tools/data_tools/filterReads.c

   1 ////////////////////////////////////////////////////////////
   2 // The purpose of this program is to read a gff and a
   3 // solexa reads file and create a data set used by QPalma.
   4 //
   5 //
   6 ////////////////////////////////////////////////////////////
   7
   8 #include <sys/mman.h>
   9 #include <sys/stat.h>
  10 #include <stdio.h>
  11 #include <stdlib.h>
  12 #include <string.h>
  13 #include <unistd.h>
  14 #include <math.h>
  15
  16 #include "datastructures.h"
  17
  18 #define _FILE_OFFSET_BITS == 64
  19
  20 int compare_gene_struct(struct gene* a, struct gene* b) {
  21    return a->stop - b->start;
  22 }
  23
  24 int parse_gff(char* filename,FILE* fid,struct gene*** allGenes);
  25
  26 void sort_genes(struct gene*** allGenes, int numGenes);
  27
  28 void process_reads(FILE* reads_fs,struct gene*** allGenes,int numGenes,FILE* out_fs);
  29
  30 void combine_info(int exon_stop, int exon_start, void** upstream, int up_size, void** downstream, int down_size, FILE* out_fs,const char* gene_id, int exon_idx);
  31
  32 //int fitting(char* up_prb, char* up_prb_end, char* down_prb, char* down_prb_end);
  33 int fitting(char* up_prb, char* down_prb);
  34
  35 void remove_ambiguities(char * old_seq, int old_seq_size, char* new_seq);
  36
  37 static char *info = "Usage is:\n./filterReads gff reads output";
  38
  39 /*
  40  * Some constants specifying the exact behavior of the filter
  41  *
  42  */
  43
  44 const int read_size = 36;
  45
  46 const int min_overlap = 6;
  47 const int max_overlap = 30;
  48
  49 int read_nr = 1;
  50
  51 int combined_reads = 0;
  52 int main(int argc, char* argv[]) {
  53
  54    if(argc != 4) {
  55       printf("%s\n",info);
  56       exit(EXIT_FAILURE);
  57    }
  58
  59    int status;
  60    int filenameSize = 256;
  61    char* gff_filename = malloc(sizeof(char)*filenameSize);
  62    char* reads_filename = malloc(sizeof(char)*filenameSize);
  63    char* output_filename = malloc(sizeof(char)*filenameSize);
  64
  65    strncpy(gff_filename,argv[1],filenameSize);
  66    strncpy(reads_filename,argv[2],filenameSize);
  67    strncpy(output_filename,argv[3],filenameSize);
  68
  69    FILE *gff_fs = fopen(gff_filename,"r");
  70    FILE *reads_fs = fopen(reads_filename,"r");
  71    FILE *out_fs = fopen(output_filename,"w");
  72
  73    if(gff_fs == NULL) {
  74       printf("Error: Could not open file: %s",gff_filename);
  75       exit(EXIT_FAILURE);
  76    }
  77
  78    if(reads_fs == NULL) {
  79       printf("Error: Could not open file: %s",reads_filename);
  80       exit(EXIT_FAILURE);
  81    }
  82
  83    if(out_fs == NULL) {
  84       printf("Error: Could not open file: %s",output_filename);
  85       exit(EXIT_FAILURE);
  86    }
  87
  88    struct gene** allGenes;
  89    int numGenes = parse_gff(gff_filename,gff_fs,&allGenes);
  90    status = fclose(gff_fs);
  91    free(gff_filename);
  92    if(status != 0)
  93       printf("closing of gff filestream failed!\n");
  94
  95    printf("Successfully parsed gff file! Found %d genes.\n",numGenes);
  96
  97    // some entries in the gff files are not in a sorted order
  98    //printf("Sorting genes...\n");
  99    //sort_genes(&allGenes,numGenes);
 100    //qsort(allGenes,numGenes,sizeof(struct gene*),compare_gene_struct);
 101    ///printf("Genes were sorted!\n");
 102
 103    int gidx, eidx;
 104    int exon_cov = 0;
 105    for(gidx=0;gidx<numGenes;gidx++) {
 106       for(eidx=0;eidx<allGenes[gidx]->num_exons;eidx++) {
 107          exon_cov += allGenes[gidx]->exon_stops[eidx] - allGenes[gidx]->exon_starts[eidx];
 108    }}
 109    printf("Exon coverage is %f\n",(double)exon_cov/30432563);
 110
 111    process_reads(reads_fs,&allGenes,numGenes,out_fs);
 112
 113    status = fclose(reads_fs);
 114    status = fclose(out_fs);
 115    if(status != 0)
 116       perror("fclose");
 117
 118    free(reads_filename);
 119    free(output_filename);
 120    return 0;
 121 }
 122
 123 void process_reads(FILE* reads_fs,struct gene*** allGenes,int numGenes, FILE* out_fs) {
 124    int status;
 125
 126    int buffer_size= 64;
 127    int chr        = 0;
 128    int pos        = 0;
 129    char* seq      = malloc(sizeof(char)*buffer_size);
 130    unsigned long id         = 0;
 131    char strand    = ' ';
 132    int mismatch   = 0;
 133    int occurrence = 0;
 134    int size       = 0;
 135    int cut        = 0;
 136    char* prb      = malloc(sizeof(char)*buffer_size);
 137    char* cal_prb  = malloc(sizeof(char)*buffer_size);
 138    char* chastity = malloc(sizeof(char)*buffer_size);
 139
 140    int reads_fid = fileno(reads_fs);
 141    struct stat reads_stat;
 142    if ( fstat(reads_fid,&reads_stat) == -1) {
 143       perror("fstat");
 144       exit(EXIT_FAILURE);
 145    }
 146    off_t reads_filesize = reads_stat.st_size;
 147    printf("Reads file is of size %lu bytes\n",(unsigned long) reads_filesize);
 148    int numReads = reads_filesize / 178.0;
 149
 150    void *reads_area = mmap (NULL,reads_filesize,PROT_READ|PROT_WRITE,MAP_PRIVATE,reads_fid,0);
 151    if (reads_area == MAP_FAILED) {
 152       perror("mmap");
 153       exit(EXIT_FAILURE);
 154    }
 155    close(reads_fid);
 156    printf("Successfully mapped %lu bytes of reads file into memory\n",(unsigned long)reads_filesize);
 157
 158    void* linePtr = reads_area;
 159    char* current_line = malloc(sizeof(char)*256);
 160
 161    int SIZE = 500;
 162    // initialize boundary arrays
 163    void** upstream_end      = malloc(sizeof(void*)*SIZE);
 164    void** upstream_overlap  = malloc(sizeof(void*)*SIZE);
 165    void** downstream_start  = malloc(sizeof(void*)*SIZE);
 166    void** downstream_overlap= malloc(sizeof(void*)*SIZE);
 167    int ue,uo,ds,dov;
 168    ue = uo = ds = dov = 0;
 169
 170    int skippedLinesCounter = 0;
 171
 172    int prev_exon_start = -1;
 173    int prev_exon_stop = -1;
 174    int cur_exon_start = -1;
 175
 176    current_line = strncpy(current_line,linePtr,256);
 177    int gene_idx = 0;
 178    int exon_idx = 1;
 179    struct gene* currentGene = (*allGenes)[gene_idx];
 180    char* gene_id = currentGene->id;
 181
 182    //char* disamb_seq = malloc(sizeof(char)*read_size);
 183
 184    int skippedReadCtr = 0;
 185    int uselessReadCtr = 0;
 186    int exonicReadCtr = 0;
 187    int endPrevCtr = 0;
 188    int prevOverlapCtr = 0;
 189    int currentStartCtr = 0;
 190    int currentOverlapCtr = 0;
 191    int multioccurReadCtr = 0;
 192
 193    int readCtr = 0;
 194    // start of the parsing loop
 195    while(1) {
 196       if (gene_idx == numGenes || strcmp(current_line,"") == 0)
 197          break;
 198
 199       gene_id = currentGene->id;
 200
 201       if (readCtr != 0 && readCtr % 1000000 == 0)
 202          printf("Processed %d/%d genes and %d/%d reads.\n",gene_idx,numGenes,readCtr,numReads);
 203
 204       //if (gene_idx >= 1833)
 205       //   printf("currentGene start/stop: %d/%d. Positions is %d\n",currentGene->start,currentGene->stop,pos);
 206
 207       status = sscanf(current_line,"%d\t%d\t%s\t%lu\t%c\t%d\t%d\t%d\t%d\t%s\t%s\t%s\n",
 208       &chr,&pos,seq,&id,&strand,&mismatch,&occurrence,&size,&cut,prb,cal_prb,chastity);
 209       if (status < 12) {
 210          skippedLinesCounter++;
 211          goto next_read;
 212       }
 213
 214       // if the read is occurring several times elsewhere then get rid of it
 215       //if(!(occurrence >= 1 && occurrence <= 25)) {
 216       if ( occurrence != 1 ) {
 217          multioccurReadCtr++;
 218          goto next_read;
 219       }
 220
 221       if (!(currentGene->start <= pos && (pos + read_size-1) <= currentGene->stop)) { // read is not within gene borders
 222
 223          if ( currentGene->stop < (pos + read_size-1) ) { // go to next gene
 224             gene_idx++;
 225             exon_idx = 1;
 226             currentGene = (*allGenes)[gene_idx];
 227             //printf("currentGene->start / currentGene->stop %d/%d pos is %d\n",currentGene->start,currentGene->stop,pos);
 228             ue = uo = ds = dov = 0;
 229             continue;
 230          }
 231
 232          if ( pos < currentGene->start ) { // go to next read
 233             skippedReadCtr++;
 234
 235             next_read:
 236
 237             while (*(char*)linePtr != '\n') linePtr++;
 238             linePtr++;
 239             readCtr += 1;
 240             current_line = strncpy(current_line,linePtr,256);
 241             continue;
 242          }
 243
 244       } else { // read IS within gene borders
 245
 246          exon_label:
 247
 248          if (exon_idx == currentGene->num_exons) {
 249             gene_idx++;
 250             exon_idx = 1;
 251             currentGene = (*allGenes)[gene_idx];
 252             continue;
 253          }
 254
 255          prev_exon_start = currentGene->exon_starts[exon_idx-1];
 256          prev_exon_stop = currentGene->exon_stops[exon_idx-1];
 257          cur_exon_start = currentGene->exon_starts[exon_idx];
 258
 259          //printf("exon %d %d inton til %d pos %d\n",prev_exon_start,prev_exon_stop,cur_exon_start,pos);
 260
 261          if (cur_exon_start - prev_exon_stop < min_overlap || cur_exon_start < pos ) { // go to next exon
 262             if (ue != 0 && dov != 0)
 263                combine_info(prev_exon_stop,cur_exon_start,upstream_end,ue,downstream_overlap,dov,out_fs,gene_id,exon_idx);
 264
 265             if (uo != 0 && ds != 0)
 266                combine_info(prev_exon_stop,cur_exon_start,upstream_overlap,uo,downstream_start,ds,out_fs,gene_id,exon_idx);
 267
 268             exon_idx++;
 269             ue = uo = ds = dov = 0;
 270             goto exon_label;
 271          }
 272
 273          if ( prev_exon_start < pos && (pos+read_size) < prev_exon_stop ) { // read is inside previous exon
 274             exonicReadCtr++;
 275
 276    // Removed exonic reads
 277
 278             //remove_ambiguities(seq,strlen(seq),disamb_seq);
 279             //fprintf(out_fs,"%d\t%c\t%s\t%d\t%s\t%s\t%s\n",chr,strand,disamb_seq,read_size,prb,cal_prb,chastity);
 280
 281             goto next_read;
 282          }
 283
 284          if ( pos + (read_size-1) < prev_exon_stop ) // go to next read
 285             goto next_read;
 286
 287          // if this position is reached the read is somehow overlapping or
 288          // exactly on the exon boundary. now determine where exactly:
 289          if (pos + (read_size-1) == prev_exon_stop) {  // read ends at previous exon end
 290             endPrevCtr++;
 291             upstream_end[ue] = linePtr;
 292             ue++;
 293             goto next_read;
 294          }
 295
 296          if ( (prev_exon_stop - pos) >= min_overlap && (prev_exon_stop - pos) <= max_overlap) { // read overlaps with previous exon boundary
 297             prevOverlapCtr++;
 298             upstream_overlap[uo] = linePtr;
 299             uo++;
 300             goto next_read;
 301          }
 302
 303          if ( pos == cur_exon_start ) { // read starts at current exon start
 304             currentStartCtr++;
 305             downstream_start[ds] = linePtr;
 306             ds++;
 307             goto next_read;
 308          }
 309
 310          if ( (cur_exon_start - pos) >= min_overlap && (cur_exon_start - pos) <= max_overlap) { // read overlaps with current exon boundary
 311             currentOverlapCtr++;
 312             downstream_overlap[dov] = linePtr;
 313             dov++;
 314             goto next_read;
 315          }
 316
 317          uselessReadCtr++;
 318          goto next_read; // read was not useful i.e. not overlapping/starting at exon boundaries
 319       }
 320    }
 321
 322    combine_info(prev_exon_stop,cur_exon_start,upstream_end,ue,downstream_overlap,dov,out_fs,gene_id,exon_idx);
 323    combine_info(prev_exon_stop,cur_exon_start,upstream_overlap,uo,downstream_start,ds,out_fs,gene_id,exon_idx);
 324
 325    free(upstream_end);
 326    free(upstream_overlap);
 327    free(downstream_start);
 328    free(downstream_overlap);
 329
 330    printf("Processed %d reads in total while %d reads where skipped.\n",readCtr,skippedLinesCounter);
 331    printf("%d were totally intronic. %d were outside of genes and %d were occurring on several positions.\n",uselessReadCtr,skippedReadCtr,multioccurReadCtr);
 332    printf("%d reads were useless in total\n",uselessReadCtr+skippedReadCtr+multioccurReadCtr);
 333    printf("%d reads were totally exonic\n",exonicReadCtr);
 334    printf("%d reads overlap with prev. exon. %d reads overlap with current exon\n",prevOverlapCtr,currentOverlapCtr);
 335    printf("%d reads where newly combined from two original reads\n",combined_reads);
 336    printf("Total used reads: %d\n",exonicReadCtr+endPrevCtr+prevOverlapCtr+currentStartCtr+currentOverlapCtr);
 337
 338    status = munmap(reads_area,reads_filesize);
 339    if(status != 0)
 340       perror("munmap");
 341
 342    //free(current_line);
 343    free(seq);
 344    free(prb);
 345    free(cal_prb);
 346    free(chastity);
 347 }
 348
 349 void combine_info(int exon_stop, int exon_start, void** upstream, int up_size, void** downstream, int down_size,FILE* out_fs,const char* gene_id, int exon_idx) {
 350    //printf("up/down size is %d/%d\n",up_size,down_size);
 351
 352    int up_idx, down_idx, status;
 353    char* upstream_line = malloc(sizeof(char)*256);
 354    char* downstream_line = malloc(sizeof(char)*256);
 355    int p_start, p_stop;
 356
 357    int buffer_size= 64;
 358
 359    int up_chr        = 0;
 360    int up_pos        = 0;
 361    char* up_seq      = malloc(sizeof(char)*buffer_size);
 362    int up_id         = 0;
 363    char up_strand    = ' ';
 364    int up_mismatch   = 0;
 365    int up_occurrence = 0;
 366    int up_sz       = 0;
 367    int up_cut        = 0;
 368    char* up_prb      = malloc(sizeof(char)*buffer_size);
 369    char* up_cal_prb  = malloc(sizeof(char)*buffer_size);
 370    char* up_chastity = malloc(sizeof(char)*buffer_size);
 371
 372    int down_chr        = 0;
 373    int down_pos        = 0;
 374    char* down_seq      = malloc(sizeof(char)*buffer_size);
 375    int down_id         = 0;
 376    char down_strand    = ' ';
 377    int down_mismatch   = 0;
 378    int down_occurrence = 0;
 379    int down_sz         = 0;
 380    int down_cut        = 0;
 381    char* down_prb      = malloc(sizeof(char)*buffer_size);
 382    char* down_cal_prb  = malloc(sizeof(char)*buffer_size);
 383    char* down_chastity = malloc(sizeof(char)*buffer_size);
 384
 385    int new_chr        = 0;
 386    char* new_seq      = malloc(sizeof(char)*buffer_size);
 387    char new_strand    = ' ';
 388    char* new_prb        = malloc(sizeof(char)*buffer_size);
 389    char* new_cal_prb    = malloc(sizeof(char)*buffer_size);
 390    char* new_chastity   = malloc(sizeof(char)*buffer_size);
 391    char* new_up_seq = malloc(sizeof(char)*read_size);
 392    char* new_down_seq = malloc(sizeof(char)*read_size);
 393
 394    int overlap;
 395    int fit;
 396
 397    char* used_flag = calloc(down_size,sizeof(char));
 398
 399    for(up_idx=0;up_idx<up_size;up_idx++) {
 400       strncpy(upstream_line,upstream[up_idx],256);
 401       status = sscanf(upstream_line,"%d\t%d\t%s\t%d\t%c\t%d\t%d\t%d\t%d\t%s\t%s\t%s\n",
 402       &up_chr,&up_pos,up_seq,&up_id,&up_strand,&up_mismatch,&up_occurrence,&up_sz,
 403       &up_cut,up_prb,up_cal_prb,up_chastity);
 404
 405       remove_ambiguities(up_seq,strlen(up_seq),new_up_seq);
 406
 407       overlap = exon_stop - up_pos;
 408
 409       for(down_idx=0;down_idx<down_size;down_idx++) {
 410          if( used_flag[down_idx] == 1)
 411             continue;
 412
 413          strncpy(downstream_line,downstream[down_idx],256);
 414          status = sscanf(downstream_line,"%d\t%d\t%s\t%d\t%c\t%d\t%d\t%d\t%d\t%s\t%s\t%s\n",
 415          &down_chr,&down_pos,down_seq,&down_id,&down_strand,&down_mismatch,&down_occurrence,&down_sz,
 416          &down_cut,down_prb,down_cal_prb,down_chastity);
 417
 418          remove_ambiguities(down_seq,strlen(down_seq),new_down_seq);
 419
 420          if(up_strand != down_strand)
 421             continue;
 422
 423          new_seq[0] = '\0';
 424          new_prb[0] = '\0';
 425          new_cal_prb[0] = '\0';
 426          new_chastity[0] = '\0';
 427          int splitpos = 0;
 428
 429          fit = fitting(up_prb+(36-overlap)-6,down_prb+(exon_start-down_pos)-6);
 430          if (fit == -1)
 431             continue;
 432
 433          if (! (fit  < overlap ))
 434             continue;
 435
 436          new_chr     = up_chr;
 437          new_strand  = up_strand;
 438          splitpos = overlap;
 439
 440          p_start = up_pos;
 441          p_stop =  exon_start+(read_size-overlap);
 442
 443          strncat(new_seq,new_up_seq,overlap);
 444          strncat(new_prb,up_prb,overlap);
 445          strncat(new_cal_prb,up_cal_prb,overlap);
 446          strncat(new_chastity,up_chastity,overlap);
 447
 448          //strncat(new_seq,new_down_seq+fit,read_size-overlap);
 449          //strncat(new_prb,down_prb+fit,read_size-overlap);
 450          //strncat(new_cal_prb,down_cal_prb+fit,read_size-overlap);
 451          //strncat(new_chastity,down_chastity+fit,read_size-overlap);
 452
 453          int offset = exon_start-down_pos;
 454          strncat(new_seq,new_down_seq+offset,read_size-overlap);
 455          strncat(new_prb,down_prb+offset,read_size-overlap);
 456          strncat(new_cal_prb,down_cal_prb+offset,read_size-overlap);
 457          strncat(new_chastity,down_chastity+offset,read_size-overlap);
 458
 459          // The four last integers specify the positions of
 460          // p_start : the position in the dna where the first read starts
 461          // exons_stop  : the position in the dna where the first exons ends
 462          // exons_start : the position in the dna where the second exons starts
 463          // p_stop  : the position in the dna where the (truncated) second read ends
 464          fprintf(out_fs,"%08d\t%d\t%c\t%s\t%d\t%d\t%s\t%s\t%s\t%s\t%d\t%d\t%d\t%d\n",
 465          read_nr,new_chr,new_strand,new_seq,splitpos,read_size,new_prb,new_cal_prb,new_chastity,gene_id,p_start,exon_stop,exon_start,p_stop);
 466          read_nr++;
 467
 468          combined_reads++;
 469          used_flag[down_idx] = 1;
 470       }
 471    }
 472
 473    free(upstream_line);
 474    free(downstream_line);
 475    free(used_flag);
 476
 477    free(new_up_seq);
 478    free(new_down_seq);
 479
 480    free(up_seq);
 481    free(up_prb);
 482    free(up_cal_prb);
 483    free(up_chastity);
 484
 485    free(down_seq);
 486    free(down_prb);
 487    free(down_cal_prb);
 488    free(down_chastity);
 489
 490    free(new_seq);
 491    free(new_prb);
 492    free(new_cal_prb);
 493    free(new_chastity);
 494 }
 495
 496 /*
 497 int fitting(char* up_prb, char* down_prb) {
 498    int w_size = 3;
 499
 500    double current_mean_up = 0;
 501    double current_mean_down = 0;
 502
 503    int uidx;
 504    for(uidx=0;uidx<w_size;uidx++) {
 505       current_mean_up += up_prb[uidx]-50;
 506       current_mean_down += down_prb[uidx]-50;
 507    }
 508
 509    current_mean_up /= w_size;
 510    current_mean_down /= w_size;
 511
 512    double min_val = fmin(current_mean_up,current_mean_down);
 513    double max_val = fmax(current_mean_up,current_mean_down);
 514    if ( (min_val / max_val) <= 0.33 )
 515       return 1;
 516
 517    return -1;
 518 }
 519 */
 520
 521 int fitting(char* up_prb, char* down_prb) {
 522    double epsilon_mean = 30.0;
 523    double epsilon_var = 30.0;
 524    int w_size = 6;
 525
 526    double current_mean_up = 0;
 527    double current_variance_up = 0;
 528    double current_mean_down = 0;
 529    double current_variance_down = 0;
 530
 531    double* mean_up       = malloc(sizeof(double)*read_size-2*w_size);
 532    double* variance_up   = malloc(sizeof(double)*read_size-2*w_size);
 533    double* mean_down     = malloc(sizeof(double)*read_size-2*w_size);
 534    double* variance_down = malloc(sizeof(double)*read_size-2*w_size);
 535
 536    int iidx = -1;
 537    int uidx;
 538    int didx;
 539
 540    for(uidx=iidx-w_size;uidx<iidx;uidx++) {
 541       didx = uidx+w_size;
 542       current_mean_up += up_prb[uidx]-50;
 543       current_mean_down += up_prb[didx]-50;
 544    }
 545
 546    current_mean_up /= w_size;
 547    current_mean_down /= w_size;
 548
 549    for(uidx=iidx-w_size;uidx<iidx;uidx++) {
 550       didx = uidx+w_size;
 551       current_variance_up += pow(up_prb[uidx]-50 - current_mean_up,2);
 552       current_variance_down += pow(up_prb[didx]-50 - current_mean_up,2);
 553    }
 554
 555    current_variance_up /= (w_size-1);
 556    current_variance_down /= (w_size-1);
 557
 558    for(iidx=w_size;iidx<30;iidx++) {
 559       for(uidx=iidx-w_size;uidx<iidx;uidx++) {
 560          didx = uidx+w_size;
 561          mean_up[iidx-w_size] += down_prb[uidx]-50;
 562          mean_down[iidx-w_size] += down_prb[didx]-50;
 563
 564       }
 565       mean_up[iidx-w_size] /= w_size;
 566       mean_down[iidx-w_size] /= w_size;
 567
 568       for(uidx=iidx-w_size;uidx<iidx;uidx++) {
 569          didx = uidx+w_size;
 570          variance_up[iidx-w_size] += pow(down_prb[uidx]-50 - mean_up[iidx-w_size],2);
 571          variance_down[iidx-w_size] += pow(down_prb[didx]-50 - mean_down[iidx-w_size],2);
 572       }
 573       variance_up[iidx-w_size] /= (w_size-1);
 574       variance_down[iidx-w_size] /= (w_size-1);
 575
 576       //printf("means: %f %f %f %f, variances: %f %f %f %f\n",current_mean_up,current_mean_down,mean_up,mean_down,current_variance_up,current_variance_down,variance_up,variance_down);
 577       //if ( abs(current_mean_up - mean_up) < epsilon_mean && abs(current_variance_up - variance_up) < epsilon_var
 578       //&& abs(current_mean_down - mean_down) < epsilon_mean && abs(current_variance_down - variance_down) < epsilon_var )
 579       //   return iidx;
 580    }
 581
 582    int bidx;
 583    int bestIdx = -1;
 584    double min = 1000.0;
 585    for(bidx=0;bidx<read_size-2*w_size;bidx++) {
 586       if ( abs(current_mean_up - mean_up[bidx]) < epsilon_mean && abs(current_variance_up - variance_up[bidx]) < epsilon_var
 587       && abs(current_mean_down - mean_down[bidx]) < epsilon_mean && abs(current_variance_down - variance_down[bidx]) < epsilon_var ) {
 588          if ( abs(current_mean_up - mean_up[bidx]) + abs(current_variance_up - variance_up[bidx]) + abs(current_mean_down - mean_down[bidx]) + abs(current_variance_down - variance_down[bidx])) {
 589             min = abs(current_mean_up - mean_up[bidx]) + abs(current_variance_up - variance_up[bidx]) + abs(current_mean_down - mean_down[bidx]) + abs(current_variance_down - variance_down[bidx]);
 590             bestIdx = bidx;
 591          }
 592       }
 593    }
 594
 595    free(mean_up);
 596    free(variance_up);
 597    free(mean_down);
 598    free(variance_down);
 599
 600    return bestIdx;
 601 }
 602
 603 void remove_ambiguities(char * old_seq, int old_seq_size, char* new_seq) {
 604    int idx=0;
 605    int new_idx = 0;
 606
 607    while(idx<old_seq_size) {
 608       //printf("Current elem %c pos %d %d\n",old_seq[idx],idx,new_idx);
 609       if (old_seq[idx] == ']' || old_seq[idx] == '-' ) {
 610          idx++;
 611          continue;
 612       }
 613
 614       if (old_seq[idx] == '[') {
 615          idx += 2;
 616          //printf("%c %c\n",old_seq[idx-2],old_seq[idx]);
 617          continue;
 618       }
 619
 620       new_seq[new_idx] = old_seq[idx];
 621       idx++;
 622       new_idx++;
 623    }
 624
 625    if (new_idx != 36) {
 626       printf("Error: Sequence is not of length 36!\n");
 627       printf("old seq: %s\n",old_seq);
 628       printf("new seq: %s\n",new_seq);
 629       exit(EXIT_FAILURE);
 630    }
 631 }
 632
 633 /*
 634  * TODO:
 635  * - Check strand -> done simple (only if equal)
 636  * - check for [AC] and similar entries -> done simple (see function
 637  * - remove_ambiguities (exchanges [XY] by the second entry)
 638  */