alg/flp_seqcomp.cpp

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