flp_seqcomp.cc

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