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