00001 #include "log.h"
00002 #include "format.h"
00003 #include "intervalFind.h"
00004 #include <stdlib.h>
00005 #include <time.h>
00006 #include "geneFusionsConfig.h"
00007 #include "gfr.h"
00008 #include "mrf.h"
00009 #include "numUtil.h"
00010 #include <math.h>
00011 #include "linestream.h"
00012
00013 #define MILLION 1000000.0
00014
00015
00016 static int sortArrayByDASPER( GfrEntry* a, GfrEntry* b ) {
00017 return (int) (b->DASPER*MILLION - a->DASPER*MILLION);
00018 }
00019
00020 int main (int argc, char *argv[])
00021 {
00022 int i,j ;
00023 double inter, intra1, intra2, intras;
00024 double Npe = -1.0;
00025 char *line;
00026
00027 Array gfrA = arrayCreate(20, GfrEntry);
00028
00029 if (argc != 2) {
00030 usage ("Usage: %s <prefix>\nNote:prefix.meta is supposed to exist.",argv[0]);
00031 }
00032
00033
00034 gfr_init ( "-" );
00035 gfr_addNewColumnType("SPER");
00036 gfr_addNewColumnType("DASPER");
00037 gfr_addNewColumnType("RESPER");
00038 gfrA = gfr_parse();
00039
00040 Stringa buffer=stringCreate( 50 );
00041 stringPrintf( buffer, "%s.meta", argv[1]);
00042 LineStream ls = ls_createFromFile( string(buffer) );
00043 if( ls == NULL ) {
00044 die( "File not found: (%s.meta)", argv[1]);
00045 } else {
00046 while( line = ls_nextLine( ls ) ) {
00047 WordIter w = wordIterCreate( line, "\t", 1);
00048 if( strEqual( wordNext ( w ), "Mapped_reads") ) {
00049 Npe = atof( wordNext( w ) );
00050 }
00051 wordIterDestroy( w );
00052 }
00053 }
00054 if( Npe <= 0 ) die("Number of mapped reads cannot be less or equal than zero: %1.3f", Npe);
00055 ls_destroy(ls);
00056 double expSPER;
00057 double resperAvg = 0.0;
00058 for( j = 0; j < arrayMax(gfrA); j++) resperAvg += (double) arrp( gfrA, j, GfrEntry)->numInter;
00059 resperAvg /= ( arrayMax(gfrA) * Npe );
00060 GfrEntry* gfrE;
00061
00062 for( i=0; i<arrayMax(gfrA); i++) {
00063 gfrE = arrp( gfrA, i, GfrEntry );
00064 gfrE->SPER = ( ((double) gfrE->numInter ) / Npe ) * MILLION;
00065
00066 inter = (double) gfrE->numInter;
00067 intra1 = (double) gfrE->numIntra1;
00068 intra2 = (double) gfrE->numIntra2;
00069 intras = intra1*intra2;
00070
00071 expSPER = ( 2.0 * intra1 + inter) * (2.0 * intra2 + inter) / (4.0 * Npe * Npe ) * MILLION;
00072 gfrE->DASPER = (gfrE->SPER - expSPER);
00073
00074 gfrE->RESPER = gfrE->SPER / resperAvg / MILLION;
00075
00076 }
00077 arraySort( gfrA, (ARRAYORDERF)sortArrayByDASPER );
00078 printf("%s\n", gfr_writeHeader() );
00079 for( i=0; i<arrayMax(gfrA); i++) {
00080 GfrEntry* gfrE = arrp( gfrA, i, GfrEntry );
00081 printf("%s\n", gfr_writeGfrEntry( gfrE )) ;
00082 }
00083 gfr_deInit();
00084
00085 warn("%s_numberMappedReads: %1.0f", argv[0], Npe);
00086 warn("%s_numGfrEntries: %d", argv[0], arrayMax( gfrA ));
00087 arrayDestroy(gfrA);
00088 return 0;
00089 }
00090
