+ added reverse complementation to filter reads function

git-svn-id: http://svn.tuebingen.mpg.de/ag-raetsch/projects/QPalma@9860 e1793c9e-67f9-0310-80fc-b846ff1f7b36

diff --git a/tools/data_tools/filterReads.c b/tools/data_tools/filterReads.c
index 132f3dc..0e8f563 100644 (file)
--- a/tools/data_tools/filterReads.c
+++ b/tools/data_tools/filterReads.c
@@ -46,6 +46,10 @@ Tuple join_seq(char* new_seq, char* up_seq,int u_off, int u_size, char* down_seq
 
 static char *info = "Usage is:\n./filterReads gff reads output";
 
+void reverse_complement(char** seq, int len);
+
+char current_strand;
+
 /*
  * Some constants specifying the exact behavior of the filter
  *
@@ -121,7 +125,7 @@ void process_reads(FILE* reads_fs,struct gene*** allGenes,int numGenes, FILE* ou
    }
    off_t reads_filesize = reads_stat.st_size;
    printf("Reads file is of size %lu bytes\n",(unsigned long) reads_filesize);
-   int numReads = reads_filesize / 178.0;
+   //int numReads = reads_filesize / 178.0;
 
    void *reads_area = mmap (NULL,reads_filesize,PROT_READ,MAP_PRIVATE,reads_fid,0);
    if (reads_area == MAP_FAILED) {
@@ -220,6 +224,11 @@ void process_reads(FILE* reads_fs,struct gene*** allGenes,int numGenes, FILE* ou
             goto next_read;
          }
 
+         //printf("before rc %s\n",seq);
+         if (current_strand == 'P')
+            reverse_complement(&seq,strlen(seq));
+         //printf("after rc %s\n",seq);
+
          // define read start and stop positions
          read_start = pos;
          read_stop  = pos + read_size-1;
@@ -470,6 +479,7 @@ 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;
 
@@ -538,6 +548,8 @@ int join_reads(int exon_stop, int exon_start, Read* up_read, Read* down_read, FI
    strncpy(new_chastity+u_size, down_read->chastity+d_off, d_size+additional_pos);
    new_chastity[read_size+additional_pos] = '\0';
 
+   printf("old reads: %s %s\n",up_read->seq,down_read->seq);
+   printf("new read: %s %d\n",new_seq,cut_pos);
 
    int status = 1; //fitting(up_read->prb,down_read->prb,u_off,d_off,u_size,d_size);
 
@@ -569,6 +581,35 @@ int join_reads(int exon_stop, int exon_start, Read* up_read, Read* down_read, FI
    return retval;
 }
 
+static char* translation = "T G   C            A      [ ]";
+
+void reverse_complement(char** seq, int len) {
+   int idx;
+   char *temp = malloc(sizeof(char)*len);
+   for(idx=0;idx<len;idx++)
+      temp[idx] = translation[(*seq)[idx]-65];
+
+   idx=0;
+   int temp_idx=len-1;
+   while(1) {
+      if (temp[temp_idx] == ']') {
+         (*seq)[idx] = '[';
+         (*seq)[idx+1] = temp[temp_idx-2];
+         (*seq)[idx+2] = temp[temp_idx-1];
+         (*seq)[idx+3] = ']';
+         idx      += 4;
+         temp_idx -= 4;
+      } else {
+         (*seq)[idx] = temp[temp_idx];
+         idx++;
+         temp_idx--;
+      }
+
+      if (idx == len || temp_idx == -1)
+         break;
+   }
+}
+
 
 void print_read(Read* cRead) {
    printf(line_format,
@@ -600,18 +641,22 @@ void open_file(const char* filename, const char* mode, FILE** fs) {
  *
  */
 
+
 int main(int argc, char* argv[]) {
 
-   if(argc != 7) {
+   if(argc != 8) {
       printf("%s\n",info);
       exit(EXIT_FAILURE);
    }
 
+   current_strand = argv[1][0];
+   printf("Current strand is %c\n",current_strand);
+
    int status;
    int filenameSize = 256;
    char* gff_filename = malloc(sizeof(char)*filenameSize);
 
-   strncpy(gff_filename,argv[1],filenameSize);
+   strncpy(gff_filename,argv[2],filenameSize);
 
    FILE *gff_fs;
    FILE *reads_fs;
@@ -619,15 +664,15 @@ int main(int argc, char* argv[]) {
    FILE *unfiltered_out_fs;
 
    // open file streams for all needed input/output files
-   open_file(argv[1],"r",&gff_fs);
-   open_file(argv[2],"r",&reads_fs);
-   open_file(argv[3],"w",&out_fs);
-   open_file(argv[4],"w",&unfiltered_out_fs);
+   open_file(argv[2],"r",&gff_fs);
+   open_file(argv[3],"r",&reads_fs);
+   open_file(argv[4],"w",&out_fs);
+   open_file(argv[5],"w",&unfiltered_out_fs);
 
-   read_nr = strtoul(argv[5],NULL,10);
+   read_nr = strtoul(argv[6],NULL,10);
    read_nr++;
 
-   unfiltered_read_nr = strtoul(argv[6],NULL,10);
+   unfiltered_read_nr = strtoul(argv[7],NULL,10);
    unfiltered_read_nr++;
 
    // allocate and load all genes and then close the gff file stream