00001 #include "log.h"
00002 #include "format.h"
00003 #include "fasta.h"
00004 #include "bp.h"
00005 #include <stdio.h>
00006 #include <unistd.h>
00007 #include "geneFusionsConfig.h"
00008
00009
00010
00011 #define TILE_SEPARATOR " "
00012
00013
00014
00015 static char* getBreakPointSequence (char *tileCoordinate1, char *tileCoordinate2)
00016 {
00017 Stringa buffer;
00018 Stringa targetsFile;
00019 FILE *fp;
00020 Array targetSeqs;
00021 int i;
00022 Seq *currSeq;
00023 static Stringa sequence = NULL;
00024
00025 buffer = stringCreate (100);
00026 targetsFile = stringCreate (100);
00027 stringPrintf (targetsFile,"targets_%d.txt",getpid ());
00028 if (!(fp = fopen (string (targetsFile),"w")) ){
00029 die ("Unable to open target file: %s",string (targetsFile));
00030 }
00031 fprintf (fp,"%s\n%s",tileCoordinate1,tileCoordinate2);
00032 fclose (fp);
00033 stringPrintf (buffer,"%s %s/%s stdout -noMask -seqList=%s",
00034 BLAT_TWO_BIT_TO_FA,BLAT_DATA_DIR,BLAT_TWO_BIT_DATA_FILENAME,string (targetsFile));
00035 fasta_initFromPipe (string (buffer));
00036 targetSeqs = fasta_readAllSequences (0);
00037 fasta_deInit ();
00038 if (arrayMax (targetSeqs) != 2) {
00039 die ("Expected only two target sequences");
00040 }
00041 stringCreateClear (sequence,100);
00042 for (i = 0; i < arrayMax (targetSeqs); i++) {
00043 currSeq = arrp (targetSeqs,i,Seq);
00044 stringAppendf (sequence,"%s",currSeq->sequence);
00045 hlr_free (currSeq->name);
00046 hlr_free (currSeq->sequence);
00047 }
00048 arrayDestroy (targetSeqs);
00049 stringPrintf (buffer,"rm -rf %s",string (targetsFile));
00050 hlr_system (string (buffer),0);
00051 stringDestroy (targetsFile);
00052 stringDestroy (buffer);
00053 return string (sequence);
00054 }
00055
00056
00057
00058 static void extractCoordinates (char *string, int *start, int *end)
00059 {
00060 static char* stringCopy = NULL;
00061 char *pos1,*pos2;
00062
00063 strReplace (&stringCopy,string);
00064 pos1 = strchr (stringCopy,':');
00065 pos2 = strchr (stringCopy,'-');
00066 *pos2 = '\0';
00067 *start = atoi (pos1 + 1);
00068 *end = atoi (pos2 + 1);
00069 }
00070
00071
00072
00073 static int getTileSize (char *tileCoordinate1, char *tileCoordinate2)
00074 {
00075 int startTile1,endTile1,tileSize1;
00076 int startTile2,endTile2,tileSize2;
00077
00078 extractCoordinates (tileCoordinate1,&startTile1,&endTile1);
00079 tileSize1 = endTile1 - startTile1;
00080 extractCoordinates (tileCoordinate2,&startTile2,&endTile2);
00081 tileSize2 = endTile2 - startTile2;
00082 if (tileSize1 != tileSize2) {
00083
00084 }
00085 int junctionSize= tileSize1 + tileSize2;
00086 return (junctionSize/2 + junctionSize%2);
00087 }
00088
00089
00090
00091 int main (int argc, char *argv[])
00092 {
00093 Array breakPoints;
00094 BreakPoint *currBP;
00095 BreakPointRead *currBPR;
00096 int i,j,k;
00097 int readLength;
00098 int tileSize;
00099 char *breakPointSequence;
00100
00101 bp_init ("-");
00102 breakPoints = bp_getBreakPoints ();
00103 for (i = 0; i < arrayMax (breakPoints); i++) {
00104 currBP = arrp (breakPoints,i,BreakPoint);
00105 tileSize = getTileSize (currBP->tileCoordinate1,currBP->tileCoordinate2);
00106 breakPointSequence = getBreakPointSequence (currBP->tileCoordinate1,currBP->tileCoordinate2);
00107 printf ("Tile 1: %s\n",currBP->tileCoordinate1);
00108 printf ("Tile 2: %s\n",currBP->tileCoordinate2);
00109 printf ("Number of reads spanning breakpoint: %d\n\n\n",arrayMax (currBP->breakPointReads));
00110 for (j = 0; j < arrayMax (currBP->breakPointReads); j++) {
00111 currBPR = arrp (currBP->breakPointReads,j,BreakPointRead);
00112 readLength = strlen (currBPR->read);
00113 for (k = 0; k < currBPR->offset; k++) {
00114 printf (" ");
00115 }
00116 for (k = 0; k < readLength; k++) {
00117 if (((currBPR->offset + k) % tileSize) == 0 && (currBPR->offset + k) != 0) {
00118 printf ("%s",TILE_SEPARATOR);
00119 }
00120 printf ("%c",currBPR->read[k]);
00121 }
00122 printf ("\n");
00123 }
00124 for (k = 0; k < (2 * tileSize); k++) {
00125 if ((k % tileSize) == 0 && k != 0) {
00126 printf ("%s",TILE_SEPARATOR);
00127 }
00128 printf ("%c",breakPointSequence[k]);
00129 }
00130 printf ("\n\n\n\n\n");
00131 }
00132 bp_deInit ();
00133 return 0;
00134 }
00135
00136
