00001 #if defined(__cplusplus)
00002 extern "C"
00003 {
00004 #endif
00005
00006 #include "log.h"
00007 #include "format.h"
00008 #include "bp.h"
00009 #include <stdio.h>
00010 #include <stdlib.h>
00011 #include <time.h>
00012
00013 #if defined(__cplusplus)
00014 }
00015 #endif
00016
00017
00018
00019 #ifndef __CINT__
00020
00021
00022
00023 #include <TROOT.h>
00024 #include <TMath.h>
00025 #include <TRandom.h>
00026
00027 TROOT root("Rint","The ROOT Interactive Interface");
00028
00029 #endif
00030
00031
00032
00033 int main (int argc, char *argv[])
00034 {
00035 Array breakPoints;
00036 BreakPoint *currBP;
00037 BreakPointRead *currBPR;
00038 int i,j;
00039 int minNumReads,minNumUniqueOffsets,minNumReadsForKS,numPossibleOffsets;
00040 double pValueCutoffForKS;
00041 Array offsets;
00042 Array randomNumbers;
00043 double *observedOffsets;
00044 double *randomOffsets;
00045
00046 if (argc != 6) {
00047 usage ("%s <minNumReads> <minNumUniqueOffsets> <minNumReadsForKS> <pValueCutoffForKS> <numPossibleOffsets>",argv[0]);
00048 }
00049 minNumReads = atoi (argv[1]);
00050 minNumUniqueOffsets = atoi (argv[2]);
00051 minNumReadsForKS = atoi (argv[3]);
00052 pValueCutoffForKS = atof (argv[4]);
00053 numPossibleOffsets = atoi (argv[5]);
00054 bp_init ("-");
00055 offsets = arrayCreate (100,int);
00056 randomNumbers = arrayCreate (100,int);
00057 breakPoints = bp_getBreakPoints ();
00058 for (i = 0; i < arrayMax (breakPoints); i++) {
00059 currBP = arrp (breakPoints,i,BreakPoint);
00060 arrayClear (offsets);
00061 for (j = 0; j < arrayMax (currBP->breakPointReads); j++) {
00062 currBPR = arrp (currBP->breakPointReads,j,BreakPointRead);
00063 array (offsets,arrayMax (offsets),int) = currBPR->offset;
00064 }
00065 arraySort (offsets,(ARRAYORDERF)arrayIntcmp);
00066 arrayUniq (offsets,NULL,(ARRAYORDERF)arrayIntcmp);
00067 if (arrayMax (currBP->breakPointReads) >= minNumReads && arrayMax (currBP->breakPointReads) < minNumReadsForKS) {
00068 if (arrayMax (offsets) >= minNumUniqueOffsets) {
00069 puts (bp_writeBreakPoint (currBP));
00070 }
00071 }
00072 else if (arrayMax (currBP->breakPointReads) >= minNumReads && arrayMax (currBP->breakPointReads) >= minNumReadsForKS) {
00073 arrayClear (randomNumbers);
00074 for (j = 0; j < arrayMax (offsets); j++) {
00075 array (randomNumbers,arrayMax (randomNumbers),int) = rand () % numPossibleOffsets;
00076 }
00077 arraySort (randomNumbers,(ARRAYORDERF)arrayIntcmp);
00078 observedOffsets = (double*)hlr_malloc (arrayMax (offsets) * sizeof (double));
00079 randomOffsets = (double*)hlr_malloc (arrayMax (offsets) * sizeof (double));
00080 for (j = 0; j < arrayMax (offsets); j++) {
00081 observedOffsets[j] = arru (offsets,j,int);
00082 randomOffsets[j] = arru (randomNumbers,j,int);
00083 }
00084 if (pValueCutoffForKS < TMath::KolmogorovTest(arrayMax (offsets),observedOffsets,arrayMax (offsets),randomOffsets,"")) {
00085 puts (bp_writeBreakPoint (currBP));
00086 }
00087 hlr_free (observedOffsets);
00088 hlr_free (randomOffsets);
00089 }
00090 }
00091 bp_deInit ();
00092 return 0;
00093 }
00094
