00001 #include "log.h"
00002 #include "format.h"
00003 #include "bowtieParser.h"
00004
00005
00006
00007 typedef struct {
00008 char *target;
00009 char *read;
00010 int position;
00011 } BreakPoint;
00012
00013
00014
00015 typedef struct {
00016 char *tileCoordinate1;
00017 char *tileCoordinate2;
00018 Array breakPoints;
00019 } SuperBreakPoint;
00020
00021
00022
00023 static int sortBreakPointsByTargetAndOffset (BreakPoint *a, BreakPoint *b)
00024 {
00025 int diff;
00026
00027 diff = strcmp (a->target,b->target);
00028 if (diff != 0) {
00029 return diff;
00030 }
00031 return b->position - a->position;
00032 }
00033
00034
00035
00036 static int sortSuperBreakPointsBySupport (SuperBreakPoint *a, SuperBreakPoint *b)
00037 {
00038 return arrayMax (b->breakPoints) - arrayMax (a->breakPoints);
00039 }
00040
00041
00042
00043 int main (int argc, char *argv[])
00044 {
00045 BowtieQuery *currBQ;
00046 BowtieEntry *currBE;
00047 Array bowtieQueries;
00048 Array breakPoints;
00049 BreakPoint *currBP,*nextBP;
00050 Array superBreakPoints;
00051 SuperBreakPoint *currSBP;
00052 char *targetCopy;
00053 char *pos;
00054 int i,j;
00055
00056 bowtieParser_initFromFile ("-");
00057 bowtieQueries = bowtieParser_getAllQueries ();
00058 bowtieParser_deInit ();
00059 breakPoints = arrayCreate (10000,BreakPoint);
00060 for (i = 0; i < arrayMax (bowtieQueries); i++) {
00061 currBQ = arrp (bowtieQueries,i,BowtieQuery);
00062 currBE = arrp (currBQ->entries,0,BowtieEntry);
00063 currBP = arrayp (breakPoints,arrayMax (breakPoints),BreakPoint);
00064 currBP->target = hlr_strdup (currBE->chromosome);
00065 currBP->read = hlr_strdup (currBE->sequence);
00066 currBP->position = currBE->position;
00067 }
00068 arraySort (breakPoints,(ARRAYORDERF)sortBreakPointsByTargetAndOffset);
00069 targetCopy = NULL;
00070 superBreakPoints = arrayCreate (1000,SuperBreakPoint);
00071 i = 0;
00072 while (i < arrayMax (breakPoints)) {
00073 currBP = arrp (breakPoints,i,BreakPoint);
00074 currSBP = arrayp (superBreakPoints,arrayMax (superBreakPoints),SuperBreakPoint);
00075 currSBP->breakPoints = arrayCreate (100,BreakPoint*);
00076 strReplace (&targetCopy,currBP->target);
00077 pos = strchr (targetCopy,'|');
00078 if (pos == NULL) {
00079 die ("Unexpected target: %s",targetCopy);
00080 }
00081 *pos = '\0';
00082 currSBP->tileCoordinate1 = hlr_strdup (targetCopy);
00083 currSBP->tileCoordinate2 = hlr_strdup (pos + 1);
00084 array (currSBP->breakPoints,arrayMax (currSBP->breakPoints),BreakPoint*) = currBP;
00085 j = i + 1;
00086 while (j < arrayMax (breakPoints)) {
00087 nextBP = arrp (breakPoints,j,BreakPoint);
00088 if (strEqual (currBP->target,nextBP->target)) {
00089 array (currSBP->breakPoints,arrayMax (currSBP->breakPoints),BreakPoint*) = nextBP;
00090 }
00091 else {
00092 break;
00093 }
00094 j++;
00095 }
00096 i = j;
00097 }
00098 arraySort (superBreakPoints,(ARRAYORDERF)sortSuperBreakPointsBySupport);
00099 for (i = 0; i < arrayMax (superBreakPoints); i++) {
00100 currSBP = arrp (superBreakPoints,i,SuperBreakPoint);
00101 printf ("%s,%s,",currSBP->tileCoordinate1,currSBP->tileCoordinate2);
00102 for (j = 0; j < arrayMax (currSBP->breakPoints); j++) {
00103 currBP = arru (currSBP->breakPoints,j,BreakPoint*);
00104 printf ("%d:%s%s",currBP->position,currBP->read, j < arrayMax (currSBP->breakPoints) - 1 ? "|" : "\n");
00105 }
00106 }
00107
00108 return 0;
00109 }
