+ fixed some issues with the splice site scores and ugly code fragments
[qpalma.git] / tools / data_tools / filterReads.c
index 84e45ae..139c681 100644 (file)
@@ -2,13 +2,10 @@
 // The purpose of this program is to read a gff and a    
 // solexa reads file and create a data set used by QPalma. 
 //
-//
-//
 // Notes:
 //
 // - Both the read indices and the gff gene indices are one-based
 //
-//
 ////////////////////////////////////////////////////////////////////////////////
 
 #include <sys/mman.h>
@@ -58,8 +55,6 @@ char current_strand;
 //#define _FDEBUG_ 0
 //#define DEBUG_READ 38603
 
-
-//
 /*
  * TODO:
  * - Check strand -> done simple (only if equal)
@@ -453,6 +448,7 @@ void combine_info(int exon_stop, int exon_start, Read** upstream, int up_size, R
       }
    }
 
+   exit(0);
    free(up_used_flag);
    free(down_used_flag);
 }
@@ -481,7 +477,6 @@ int join_reads(int exon_stop, int exon_start, Read* up_read, Read* down_read, FI
    int down_range = down_read_stop - exon_start + 1;
    int retval;
 
-
    int u_size, d_size;
    u_size = d_size = -1;
 
@@ -505,9 +500,7 @@ int join_reads(int exon_stop, int exon_start, Read* up_read, Read* down_read, FI
    int u_off = p_start - up_read_start;
    int d_off = exon_start - down_read_start;
 
-   FA(u_off >= 0);
-   FA(d_off >= 0);
-
+   FA( u_off >= 0 && d_off >= 0 );
    FA( exon_stop - p_start + p_stop - exon_start + 2 == read_size);
    FA( u_size + d_size == read_size );
 
@@ -517,51 +510,84 @@ int join_reads(int exon_stop, int exon_start, Read* up_read, Read* down_read, FI
    char* new_seq = malloc(sizeof(char)*buf_size);
    memset(new_seq,'z',sizeof(char)*buf_size);
 
+   if ( current_strand == 'P' ) {
+      printf("flipping read sequences...\n");
+      printf("%s %s\n",up_read->seq,down_read->seq);
+      char *tmp = malloc(sizeof(char)*(strlen(up_read->seq)+1));
+      strncpy(tmp,up_read->seq,strlen(up_read->seq));
+      tmp[strlen(up_read->seq)]='\0';
+
+      realloc(up_read->seq,sizeof(char)*(strlen(down_read->seq)+1));
+      strncpy(up_read->seq,down_read->seq,strlen(down_read->seq));
+      up_read->seq[strlen(down_read->seq)] = '\0';
+
+      if (up_read->seq == NULL)
+         perror("realloc\n");
+
+      realloc(down_read->seq,sizeof(char)*strlen(tmp)+1);
+      strncpy(down_read->seq,tmp,strlen(tmp));
+      down_read->seq[strlen(tmp)] = '\0';
+      
+      if (down_read->seq == NULL)
+         perror("realloc\n");
+      
+      free(tmp);
+      printf("flipping done...\n");
+      printf("%s %s\n",up_read->seq,down_read->seq);
+   }
+   
+   printf("start joining...\n");
    Tuple jinfo = join_seq(new_seq,up_read->seq,u_off,u_size,down_read->seq,d_off,d_size,down_range);
+   printf("end of joining...\n");
+
+   printf("new_seq is %s\n",new_seq);
+
 
    int cut_pos = jinfo.first;
    int additional_pos = jinfo.second;
+   printf("jinfo contains %d/%d\n",jinfo.first,jinfo.second);
 
    buf_size = read_size+1+additional_pos;
 
+   printf("allocating quality arrays (size=%d)...\n",buf_size);
    char* new_prb        = malloc(sizeof(char)*buf_size);
    char* new_cal_prb    = malloc(sizeof(char)*buf_size);
    char* new_chastity   = malloc(sizeof(char)*buf_size);
 
-   if( jinfo.first == -1 ) {
-      retval = 0;
-      goto free;
-   }
-
-   if( additional_pos > (down_range - d_size) ) {
+   if (new_prb == NULL || new_cal_prb == NULL || new_chastity == NULL)
+      perror("malloc\n");
+   
+   if( jinfo.first == -1 || (additional_pos > (down_range - d_size)) ) {
       retval = 0;
       goto free;
    }
 
+   printf("joining qualities...\n");
    strncpy(new_prb, up_read->prb+u_off, u_size);
    strncpy(new_prb+u_size, down_read->prb+d_off, d_size+additional_pos);
-   new_prb[read_size+additional_pos] = '\0';
+   new_prb[buf_size-1] = '\0';
 
    strncpy(new_cal_prb, up_read->cal_prb+u_off, u_size);
    strncpy(new_cal_prb+u_size, down_read->cal_prb+d_off, d_size+additional_pos);
-   new_cal_prb[read_size+additional_pos] = '\0';
+   new_cal_prb[buf_size-1] = '\0';
 
    strncpy(new_chastity, up_read->chastity+u_off, u_size);
    strncpy(new_chastity+u_size, down_read->chastity+d_off, d_size+additional_pos);
-   new_chastity[read_size+additional_pos] = '\0';
+   new_chastity[buf_size-1] = '\0';
+   printf("end of joining qualities...\n");
 
    //printf("old reads: %s %s (%d %d %d/%d)\n",up_read->seq,down_read->seq,up_read->pos,down_read->pos,u_off,d_off);
    //printf("new read: %s %d %d\n",new_seq,cut_pos,u_size);
-
-   int alpha   = read_size - u_off - 1;
-   int beta    = alpha - u_size ;
-   p_start     = up_read->pos + beta + 1;
-   exon_stop   = up_read->pos + alpha;
-
-   alpha   = read_size - d_off - 1;
-   beta    = alpha - (read_size - u_size);
-   exon_start  = down_read->pos + beta + 1; 
-   p_stop      = down_read->pos + alpha;
+   //if ( current_strand == 'P' ) {
+   //   int alpha   = read_size - u_off - 1;
+   //   int beta    = alpha - u_size ;
+   //   p_start     = up_read->pos + beta + 1;
+   //   exon_stop   = up_read->pos + alpha;
+   //   alpha   = read_size - d_off - 1;
+   //   beta    = alpha - (read_size - u_size);
+   //   exon_start  = down_read->pos + beta + 1; 
+   //   p_stop      = down_read->pos + alpha;
+   //}
 
    int status = 1; //fitting(up_read->prb,down_read->prb,u_off,d_off,u_size,d_size);
 
@@ -614,6 +640,7 @@ void reverse_complement(char** seq, int len) {
       if (idx == len || temp_idx == -1)
          break;
    }
+   free(temp);
 }
 
 
@@ -647,7 +674,6 @@ void open_file(const char* filename, const char* mode, FILE** fs) {
  *
  */
 
-
 int main(int argc, char* argv[]) {
 
    if(argc != 8) {
@@ -706,6 +732,8 @@ int main(int argc, char* argv[]) {
       if (! (old_gene_stop <  currentGene->start ) ) {
          printf("Gene %s is overlapping\n",currentGene->id);
          old_gene_stop = currentGene->stop;
+         free_gene(allGenes[g_idx]);
+         free(allGenes[g_idx]);
          allGenes[g_idx] = 0;
          nulled_genes++;
          continue;