+ minor changes in the paths
[qpalma.git] / tools / run_specific_scripts / transcriptome_analysis / compare_predictions / compare.c
index fe51eaa..86e174e 100644 (file)
@@ -25,19 +25,90 @@ int* pred_intron_starts;
 int* pred_intron_stops;
 int pred_size;
 
-int matching_introns;
+int matching_values;
 
 static char* info = "Usage: ./<app name> <input 1 file name> <input 2 file name> <result file name>";
 
+static int chromo_sizes[5] = {30432563, 19705359, 23470805, 18585042, 26992728};
+
+short strand_info;
+
+
+void check_seeds() {
+
+   matching_values = 0;
+   int pred_idx, gt_idx, chromo, tmp;
+   int pos,gt_start,gt_stop,pred_start,pred_stop;
+
+   printf("sizes are %d/%d (gt/pred)\n",gt_size,pred_size);
+
+   for(pred_idx=0;pred_idx<pred_size;pred_idx++) {
+
+      pred_start = pos-1500;
+      pred_stop = pos+1500;
+      pos = pred_intron_stops[pred_idx];
+      chromo = pred_intron_starts[pred_idx];
+
+      for(gt_idx=0;gt_idx<gt_size;gt_idx++) {
+
+         gt_start = gt_intron_starts[gt_idx];
+         gt_stop = gt_intron_stops[gt_idx];
+
+         if(strand_info == 0) {
+            tmp = gt_start;
+            gt_start = chromo_sizes[chromo] - gt_stop;
+            gt_stop =  chromo_sizes[chromo] - tmp;
+         }
+
+         //printf("pos: %d/%d %d/%d\n",gt_start,gt_stop,pred_start,pred_stop);
+         if ( pred_start <= gt_start && gt_stop <= pred_stop ) {
+            matching_values++;
+            break;
+         }
+
+         //if (gt_start == pred_start && gt_stop == pred_stop)
+         //   matching_values++;
+   }}
+}
+
+
+void compare_introns() {
+
+   matching_values = 0;
+   int pred_idx, gt_idx;
+   int gt_start,gt_stop,pred_start,pred_stop;
+
+   printf("sizes are %d/%d (gt/pred)\n",gt_size,pred_size);
+
+   for(pred_idx=0;pred_idx<pred_size;pred_idx++) {
+
+      pred_start = pred_intron_starts[pred_idx];
+      pred_stop = pred_intron_stops[pred_idx];
+
+      for(gt_idx=0;gt_idx<gt_size;gt_idx++) {
+
+         gt_start = gt_intron_starts[gt_idx];
+         gt_stop = gt_intron_stops[gt_idx];
+
+         //printf("pos: %d/%d %d/%d\n",gt_start,gt_stop,pred_start,pred_stop);
+         if ( fabs(gt_start - pred_start) <= 2 && fabs(gt_stop - pred_stop) <= 2) {
+            matching_values++;
+            break;
+         }
+         //if (gt_start == pred_start && gt_stop == pred_stop)
+         //   matching_values++;
+   }}
+}
+
 
 /* 
  * This function compares the intron boundaries of the two files.
  *
  */
 
-void compare_introns() {
+void _compare_introns() {
 
-   matching_introns = 0;
+   matching_values = 0;
    
    int gt_idx = 0;
    int pred_idx = 0;
@@ -57,7 +128,18 @@ void compare_introns() {
 
       //printf("pos: %d/%d %d/%d\n",gt_start,gt_stop,pred_start,pred_stop);
       if (gt_start == pred_start && gt_stop == pred_stop) {
-         matching_introns++;
+         matching_values++;
+
+         // take doubles into account
+#if 1
+         int idx = pred_idx+1;
+         for(;idx<pred_size;idx++) {
+            if (pred_intron_starts[idx] == pred_start && pred_intron_stops[idx] == pred_stop)
+               matching_values++;
+            else
+               break;
+         }
+#endif
          gt_idx++;
          pred_idx++;
          continue;
@@ -110,6 +192,7 @@ void load_introns(char* fn, int** starts, int** stops, int* size) {
    while( getline(&current_line,&line_size,fs) != -1 ) (*size)++;
 
    // now we allocate memory to store all positions
+   //*id      = malloc(sizeof(char)*(*size));
    *starts  = malloc(sizeof(int)*(*size));
    *stops   = malloc(sizeof(int)*(*size));
 
@@ -122,6 +205,7 @@ void load_introns(char* fn, int** starts, int** stops, int* size) {
       exit(EXIT_FAILURE);
    }
 
+   //char* id_string  = malloc(sizeof(char)*line_size);
    char* first_num  = malloc(sizeof(char)*line_size);
    char* second_num = malloc(sizeof(char)*line_size);
 
@@ -152,7 +236,8 @@ void load_introns(char* fn, int** starts, int** stops, int* size) {
       ctr++;
    }
 
-   printf("ctr is %d\n",ctr);
+   //printf("ctr is %d\n",ctr);
+
    if (fclose(fs) != 0)
       perror("Closing of filestream failed!");
 }
@@ -160,7 +245,7 @@ void load_introns(char* fn, int** starts, int** stops, int* size) {
 
 int main(int argc, char* argv[]) {
 
-   if(argc != 4) {
+   if(argc != 6) {
       printf("%s\n",info);
       exit(EXIT_FAILURE);
    }
@@ -173,9 +258,29 @@ int main(int argc, char* argv[]) {
    if ( gt_fn == NULL || pred_fn == NULL || result_fn == NULL)
       perror("Could not allocate memory for filenames");
 
-   strncpy(gt_fn,argv[1],strlen(argv[1]));
-   strncpy(pred_fn,argv[2],strlen(argv[2]));
-   strncpy(result_fn,argv[3],strlen(argv[3]));
+   short intron_mode = -1;
+   if ( strcmp(argv[1],"-intron") == 0 )
+      intron_mode = 1;
+
+   if ( strcmp(argv[1],"-seed") == 0 )
+      intron_mode = 0;
+
+   if (intron_mode == -1)
+      perror("You have to choose between -intron or -seed mode!");
+
+   strand_info = -1;
+   if ( strcmp(argv[2],"-strand=+") == 0 )
+      strand_info = 1;
+
+   if ( strcmp(argv[2],"-strand=-") == 0 )
+      strand_info = 0;
+
+   if (strand_info == -1)
+      perror("You have to choose between -strand=+ or -strand=- !");
+
+   strncpy(gt_fn,argv[3],strlen(argv[3]));
+   strncpy(pred_fn,argv[4],strlen(argv[4]));
+   strncpy(result_fn,argv[5],strlen(argv[5]));
 
    FILE *result_fs = fopen(result_fn,"w+");
    if(result_fs == NULL) {
@@ -186,7 +291,10 @@ int main(int argc, char* argv[]) {
    load_introns(gt_fn,&gt_intron_starts,&gt_intron_stops,&gt_size);
    load_introns(pred_fn,&pred_intron_starts,&pred_intron_stops,&pred_size);
 
-   compare_introns();
+   if (intron_mode == 1)
+      compare_introns();
+   else
+      check_seeds();
 
    int f_status = fclose(result_fs);
    if(f_status != 0)
@@ -196,6 +304,6 @@ int main(int argc, char* argv[]) {
    free(pred_fn);
    free(result_fn);
 
-   printf("Found %d matching intron(s).\n",matching_introns);
+   printf("Found %d matching intron(s)/seed(s).\n",matching_values);
    exit(EXIT_SUCCESS);
 }