00001 #include "log.h"
00002 #include "format.h"
00003 #include "gfr.h"
00004
00005
00006 int main (int argc, char *argv[])
00007 {
00008 GfrEntry *currGE;
00009 int count;
00010 int countRemoved;
00011 int i;
00012
00013 if (argc != 3) {
00014 usage ("%s <offsetCutoff> <minNumUniqueReads>",argv[0]);
00015 }
00016 count = 0;
00017 countRemoved = 0;
00018 int offsetCutOff = atoi (argv[1]);
00019 int minNumUniqueReads = atoi (argv[2]);
00020 gfr_init ("-");
00021 puts (gfr_writeHeader ());
00022 while (currGE = gfr_nextEntry ()){
00023 Array starts = arrayCreate( 100, int);
00024 for( i = 0; i < arrayMax( currGE->interReads ); i++) {
00025 int currStart = arrp(currGE->interReads, i, GfrInterRead)->readStart1 + arrp(currGE->interReads, i, GfrInterRead)->readStart2;
00026 array(starts, arrayMax(starts), int) = currStart;
00027 }
00028 arraySort( starts, (ARRAYORDERF) arrayIntcmp );
00029 arrayUniq( starts, NULL, (ARRAYORDERF) arrayIntcmp ) ;
00030 int numUniqeOffsets = arrayMax( starts );
00031 arrayDestroy( starts );
00032
00033 if( arrayMax( currGE->readsTranscript1 ) != arrayMax( currGE->readsTranscript2 ) )
00034 die( "The two ends have a different number of reads");
00035 Texta reads = textCreate( arrayMax( currGE->readsTranscript1 ) );
00036 for( i=0; i<arrayMax( currGE->readsTranscript1 ); i++) {
00037 Stringa strA = stringCreate( strlen(textItem( currGE->readsTranscript1, i) ) * 2 + 1);
00038 stringAppendf( strA, textItem( currGE->readsTranscript1,i));
00039 stringAppendf( strA, textItem( currGE->readsTranscript2,i));
00040 textAdd( reads, string(strA));
00041 stringDestroy( strA );
00042 }
00043 textUniqKeepOrder( reads );
00044 int numRemaining = arrayMax( reads );
00045 textDestroy ( reads );
00046
00047 if( numRemaining <= minNumUniqueReads || numUniqeOffsets <= offsetCutOff )
00048 {
00049 countRemoved++;
00050 continue;
00051 }
00052 puts (gfr_writeGfrEntry (currGE));
00053 count++;
00054 }
00055 gfr_deInit ();
00056 warn ("%s_PCRFilter: offset=%d minNumUniqueReads=%d",argv[0],offsetCutOff, minNumUniqueReads);
00057 warn ("%s_numRemoved: %d",argv[0],countRemoved);
00058 warn ("%s_numGfrEntries: %d",argv[0],count);
00059 return 0;
00060 }
