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