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