flp_seqcomp.cc

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