flp_seqcomp.cc

   1 //                        ----------------------------------------
   2 //                         ---------- flp_seqcomp.cc  -----------
   3 //                        ----------------------------------------
   4
   5 #include "flp.hh"
   6
   7 // Note one recording RC matches.  This version of 'seqcomp' records RC matches
   8 // as index of the window start in its nonRC form.  However, past versions of
   9 // the analysis recorded it as the index of the rear of the window, as the
  10 // following legacy comment indicates.  FR still uses this method
  11
  12 // Titus thinks the start of an RC window should be its "rear" in the
  13 // normal version of the sequence.  He's wrong, of course, and in the
  14 // future my followers will annihilate his in the Final Battle, but for
  15 // now I'll let him have his way.
  16
  17 // note seq2_i is actually the index of the leaving bp, so need to +1
  18
  19 inline void
  20 FLPs::add(int seq1_i, int seq2_i, int a_score, int i2_offset)
  21 {
  22   match a_match;
  23
  24
  25   if (a_score >= hard_threshold)
  26   {
  27     a_match.index = seq2_i + i2_offset;
  28     a_match.score = a_score;
  29
  30     all_matches[seq1_i].push_back(a_match);
  31   }
  32 }
  33
  34
  35 void
  36 FLPs::seqcomp(string sseq1, string sseq2, bool is_RC)
  37 {
  38   int start_i, seq1_i, seq2_i, win_i;      // loop variables
  39   int matches;                    // number of matches in to a window
  40   int i2_offset;
  41   char * seq1, * seq2;
  42
  43
  44   seq1 = (char *) sseq1.c_str();
  45   seq2 = (char *) sseq2.c_str();
  46
  47
  48   if (is_RC)
  49     i2_offset = window_size - seq2_length;
  50   else
  51     i2_offset = 0;
  52
  53
  54   // this does the "upper diagonals" of the search
  55
  56   // loop thru the start positions for sequence 1
  57   for(start_i = 0; start_i < seq1_win_num; start_i++)
  58   {
  59     matches = 0;
  60     // compare initial window
  61     for(win_i = 0; win_i < window_size; win_i++)
  62       if ((seq1[start_i+win_i] == seq2[win_i]) &&
  63           (seq1[start_i+win_i] != 'N'))       // N's match nothing **
  64         matches++;
  65
  66     //seq2_i is always 0 for initial win
  67     seq2_i = 0;
  68
  69     // call inlined function that adds match if it is above threshold
  70     add(start_i, seq2_i, matches, i2_offset);
  71
  72     // run through rest of seq1 and seq2 with current seq1 offset (ie start_i)
  73     // compare the bps leaving and entering the window and modify matches count
  74     seq1_i = start_i;
  75     while( (seq1_i < seq1_win_num-1) && (seq2_i < seq2_win_num-1) )
  76     {
  77       // copmpare the bp leaving the window
  78       if ((seq1[seq1_i] == seq2[seq2_i]) &&
  79           (seq1[seq1_i] != 'N'))       // N's match nothing
  80         matches = matches -1;
  81
  82       // compare the bp entering the window
  83       if ((seq1[seq1_i+window_size] == seq2[seq2_i+window_size]) &&
  84           (seq1[seq1_i+window_size] != 'N'))       // N's match nothing
  85         matches = matches + 1;
  86
  87       add(seq1_i + 1, seq2_i + 1, matches, i2_offset);
  88
  89       seq1_i = seq1_i + 1;     // increment seq1_i to next window
  90       seq2_i = seq2_i + 1;     // increment seq2_i to next window
  91     } // end of loop thru the current offset sequence
  92   } // end of loop thru the start positions of seq1 sequence
  93
  94
  95   // this does the "lower diagonals" of the search
  96
  97   // loop thru the start positions for sequence 1
  98   for(start_i = 1; start_i < seq2_win_num; start_i++)
  99   {
 100     matches = 0;
 101
 102     // compare initial window
 103     for(win_i = 0; win_i < window_size; win_i++)
 104       if ((seq2[start_i+win_i] == seq1[win_i]) &&
 105           (seq2[start_i+win_i] != 'N'))       // N's match nothing
 106         matches++;
 107
 108     //seq2_i is always start_i for initial window
 109     seq2_i = start_i;
 110
 111     add(0, seq2_i, matches, i2_offset);
 112
 113     // run through rest of seq1 and seq2 with current seq1 offset (ie start_i)
 114     // compare the bps leaving and entering the window and modify matches count
 115     seq1_i = 0;
 116     while( (seq1_i < seq1_win_num-1) && (seq2_i < seq2_win_num-1) )
 117     {
 118       // copmpare the bp leaving the window
 119       if ((seq1[seq1_i] == seq2[seq2_i]) &&
 120           (seq1[seq1_i] != 'N'))       // N's match nothing
 121         matches = matches - 1;
 122
 123       // compare the bp entering the window
 124       if ((seq1[seq1_i+window_size] == seq2[seq2_i+window_size]) &&
 125           (seq1[seq1_i+window_size] != 'N'))       // N's match nothing
 126         matches = matches + 1;
 127
 128       add(seq1_i + 1, seq2_i + 1, matches, i2_offset);
 129
 130       seq1_i = seq1_i + 1;     // increment seq1_i to next window
 131       seq2_i = seq2_i + 1;     // increment seq2_i to next window
 132     } // end of loop thru the current offset sequence
 133   } // end of loop thru the start positions of seq2 sequence
 134 }
 135
 136 /*
 137       cout << "fee\n";
 138       cout << "fie\n";
 139       cout << "foe\n";
 140       cout << "fum\n";
 141 */