5 * Created by Sarah Westcott on 7/22/09.
6 * Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
13 #include "eachgapdist.h"
14 #include "ignoregaps.h"
17 //***************************************************************************************************************
18 void DeCalculator::setMask(string m) {
24 if (seqMask.length() != 0) {
25 //whereever there is a base in the mask, save that value is query and subject
26 for (int i = 0; i < seqMask.length(); i++) {
27 if (isalpha(seqMask[i])) {
35 for (int i = 0; i < alignLength; i++) {
43 errorOut(e, "DeCalculator", "setMask");
47 //***************************************************************************************************************
48 void DeCalculator::runMask(Sequence* seq) {
51 string q = seq->getAligned();
52 string tempQuery = "";
54 //whereever there is a base in the mask, save that value is query and subject
55 set<int>::iterator setit;
56 for ( setit=h.begin() ; setit != h.end(); setit++ ) {
57 tempQuery += q[*setit];
61 seq->setAligned(tempQuery);
62 seq->setUnaligned(tempQuery);
65 errorOut(e, "DeCalculator", "runMask");
69 //***************************************************************************************************************
70 //num is query's spot in querySeqs
71 void DeCalculator::trimSeqs(Sequence* query, Sequence* subject, map<int, int>& trim) {
74 string q = query->getAligned();
75 string s = subject->getAligned();
78 for (int i = 0; i < q.length(); i++) {
79 //cout << "query = " << q[i] << " subject = " << s[i] << endl;
80 if (isalpha(q[i]) && isalpha(s[i])) { front = i; break; }
82 //cout << endl << endl;
84 for (int i = q.length(); i >= 0; i--) {
85 //cout << "query = " << q[i] << " subject = " << s[i] << endl;
86 if (isalpha(q[i]) && isalpha(s[i])) { back = i; break; }
93 errorOut(e, "DeCalculator", "trimSeqs");
97 //***************************************************************************************************************
98 vector<int> DeCalculator::findWindows(Sequence* query, int front, int back, int& size, int increment) {
103 int cutoff = back - front; //back - front
105 //if window is set to default
106 if (size == 0) { if (cutoff > 1200) { size = 300; }
107 else{ size = (cutoff / 4); } }
108 else if (size > (cutoff / 4)) {
109 mothurOut("You have selected too large a window size for sequence " + query->getName() + ". I will choose an appropriate window size."); mothurOutEndLine();
113 /* string seq = query->getAligned().substr(front, cutoff);
117 for (int l = 0; l < seq.length(); l++) { if (isalpha(seq[l])) { numBases++; } }
118 //cout << "num Bases = " << numBases << endl;
120 win.push_back(front);
121 //cout << front << '\t';
122 //move ahead increment bases at a time until all bases are in a window
124 int totalBases = 0; //used to eliminate window of blanks at end of sequence
126 seq = query->getAligned();
127 for (int m = front; m < (back - size) ; m++) {
129 //count number of bases you see
130 if (isalpha(seq[m])) { countBases++; }
132 //if you have seen enough bases to make a new window
133 if (countBases >= increment) {
134 //total bases is the number of bases in a window already.
135 totalBases += countBases;
136 //cout << "total bases = " << totalBases << endl;
137 win.push_back(m); //save spot in alignment
139 countBases = 0; //reset bases you've seen in this window
142 //no need to continue if all your bases are in a window
143 if (totalBases == numBases) { break; }
147 //get last window if needed
148 if (totalBases < numBases) { win.push_back(back-size); }
149 //cout << endl << endl;
152 //this follows wigeon, but we may want to consider that it chops off the end values if the sequence cannot be evenly divided into steps
153 for (int m = front; m < (back - size) ; m+=increment) { win.push_back(m); }
160 catch(exception& e) {
161 errorOut(e, "DeCalculator", "findWindows");
166 //***************************************************************************************************************
167 vector<float> DeCalculator::calcObserved(Sequence* query, Sequence* subject, vector<int> window, int size) {
172 for (int m = 0; m < window.size(); m++) {
174 string seqFrag = query->getAligned().substr(window[m], size);
175 string seqFragsub = subject->getAligned().substr(window[m], size);
178 for (int b = 0; b < seqFrag.length(); b++) {
179 //if at least one is a base and they are not equal
180 if( (isalpha(seqFrag[b]) || isalpha(seqFragsub[b])) && (seqFrag[b] != seqFragsub[b]) ) { diff++; }
183 //if((!isalpha(seqFrag[b])) && (!isalpha(seqFragsub[b]))) { gaps++; }
186 //percentage of mismatched bases
189 //if the whole fragment is 0 distance = 0
190 //if ((seqFrag.length()-gaps) == 0) { dist = 0.0; }
192 //percentage of mismatched bases
193 //else { dist = diff / (float) (seqFrag.length()-gaps) * 100; }
195 dist = diff / (float) (seqFrag.length()) * 100;
197 temp.push_back(dist);
202 catch(exception& e) {
203 errorOut(e, "DeCalculator", "calcObserved");
207 //***************************************************************************************************************
208 float DeCalculator::calcDist(Sequence* query, Sequence* subject, int front, int back) {
211 //so you only look at the trimmed part of the sequence
212 int cutoff = back - front;
215 //from first startpoint with length back-front
216 string seqFrag = query->getAligned().substr(front, cutoff);
217 string seqFragsub = subject->getAligned().substr(front, cutoff);
220 for (int b = 0; b < seqFrag.length(); b++) {
222 if((!isalpha(seqFrag[b])) && (!isalpha(seqFragsub[b]))) { gaps++; }
223 if (seqFrag[b] != seqFragsub[b]) { diff++; }
226 //if the whole fragment is 0 distance = 0
227 if ((seqFrag.length()-gaps) == 0) { return 0.0; }
229 //percentage of mismatched bases
230 float dist = diff / (float) (seqFrag.length()-gaps) * 100;
234 catch(exception& e) {
235 errorOut(e, "DeCalculator", "calcDist");
240 //***************************************************************************************************************
241 vector<float> DeCalculator::calcExpected(vector<float> qav, float coef) {
245 vector<float> queryExpected;
247 for (int m = 0; m < qav.size(); m++) {
249 float expected = qav[m] * coef;
251 queryExpected.push_back(expected);
254 return queryExpected;
257 catch(exception& e) {
258 errorOut(e, "DeCalculator", "calcExpected");
262 //***************************************************************************************************************
263 float DeCalculator::calcDE(vector<float> obs, vector<float> exp) {
267 float sum = 0.0; //sum = sum from 1 to m of (oi-ei)^2
269 for (int m = 0; m < obs.size(); m++) {
271 //if (obs[m] != 0.0) {
272 sum += ((obs[m] - exp[m]) * (obs[m] - exp[m]));
273 //}else { numZeros++; }
277 float de = sqrt((sum / (obs.size() - 1 - numZeros)));
281 catch(exception& e) {
282 errorOut(e, "DeCalculator", "calcDE");
287 //***************************************************************************************************************
289 vector<float> DeCalculator::calcFreq(vector<Sequence*> seqs, string filename) {
293 string freqfile = getRootName(filename) + "freq";
296 openOutputFile(freqfile, outFreq);
298 string length = toString(seqs.size()); //if there are 5000 seqs in the template then set precision to 3
299 int precision = length.length() - 1;
302 outFreq.setf(ios::fixed, ios::floatfield); outFreq.setf(ios::showpoint);
304 //at each position in the sequence
305 for (int i = 0; i < seqs[0]->getAligned().length(); i++) {
307 vector<int> freq; freq.resize(4,0);
310 //find the frequency of each nucleotide
311 for (int j = 0; j < seqs.size(); j++) {
313 char value = seqs[j]->getAligned()[i];
315 if(toupper(value) == 'A') { freq[0]++; }
316 else if(toupper(value) == 'T' || toupper(value) == 'U') { freq[1]++; }
317 else if(toupper(value) == 'G') { freq[2]++; }
318 else if(toupper(value) == 'C') { freq[3]++; }
322 //find base with highest frequency
324 for (int m = 0; m < freq.size(); m++) { if (freq[m] > highest) { highest = freq[m]; } }
326 float highFreq = highest / (float) (seqs.size());
329 Pi = (highFreq - 0.25) / 0.75;
331 //cannot have probability less than 0.
332 if (Pi < 0) { Pi = 0.0; }
334 //saves this for later
335 outFreq << setprecision(precision) << i << '\t' << highFreq << endl;
337 if (h.count(i) > 0) {
347 catch(exception& e) {
348 errorOut(e, "DeCalculator", "calcFreq");
352 //***************************************************************************************************************
353 vector<float> DeCalculator::findQav(vector<int> window, int size, vector<float> probabilityProfile) {
355 vector<float> averages;
357 //for each window find average
358 for (int m = 0; m < window.size(); m++) {
362 //while you are in the window for this sequence
364 for (int j = window[m]; j < (window[m]+size); j++) {
365 average += probabilityProfile[j];
369 average = average / count;
371 //save this windows average
372 averages.push_back(average);
377 catch(exception& e) {
378 errorOut(e, "DeCalculator", "findQav");
382 //***************************************************************************************************************
383 //seqs have already been masked
384 vector< vector<quanMember> > DeCalculator::getQuantiles(vector<Sequence*> seqs, vector<int> windowSizesTemplate, int window, vector<float> probProfile, int increment, int start, int end) {
386 vector< vector<quanMember> > quan;
388 //percentage of mismatched pairs 1 to 100
391 //string out = "getQuantiles.out";
392 //openOutputFile(out, o);
395 for(int i = start; i < end; i++){
397 mothurOut("Processing sequence " + toString(i)); mothurOutEndLine();
398 Sequence* query = new Sequence(seqs[i]->getName(), seqs[i]->getAligned());
400 //compare to every other sequence in template
401 for(int j = 0; j < i; j++){
403 Sequence* subject = new Sequence(seqs[j]->getName(), seqs[j]->getAligned());
406 map<int, int>::iterator it;
408 trimSeqs(query, subject, trim);
411 int front = it->first; int back = it->second;
413 //reset window for each new comparison
414 windowSizesTemplate[i] = window;
416 vector<int> win = findWindows(query, front, back, windowSizesTemplate[i], increment);
418 vector<float> obsi = calcObserved(query, subject, win, windowSizesTemplate[i]);
420 vector<float> q = findQav(win, windowSizesTemplate[i], probProfile);
422 float alpha = getCoef(obsi, q);
424 vector<float> exp = calcExpected(q, alpha);
426 float de = calcDE(obsi, exp);
428 float dist = calcDist(query, subject, front, back);
429 //o << i << '\t' << j << '\t' << dist << '\t' << de << endl;
432 quanMember newScore(de, i, j);
434 quan[dist].push_back(newScore);
445 catch(exception& e) {
446 errorOut(e, "DeCalculator", "getQuantiles");
450 //********************************************************************************************************************
451 //sorts lowest to highest
452 inline bool compareQuanMembers(quanMember left, quanMember right){
453 return (left.score < right.score);
455 //***************************************************************************************************************
456 //this was going to be used by pintail to increase the sensitivity of the chimera detection, but it wasn't quite right. may want to revisit in the future...
457 void DeCalculator::removeObviousOutliers(vector< vector<quanMember> >& quantiles, int num) {
460 for (int i = 0; i < quantiles.size(); i++) {
462 //find mean of this quantile score
463 sort(quantiles[i].begin(), quantiles[i].end(), compareQuanMembers);
465 float high = quantiles[i][int(quantiles[i].size() * 0.99)].score;
466 float low = quantiles[i][int(quantiles[i].size() * 0.01)].score;
468 vector<quanMember> temp;
470 //look at each value in quantiles to see if it is an outlier
471 for (int j = 0; j < quantiles[i].size(); j++) {
472 //is this score between 1 and 99%
473 if ((quantiles[i][j].score > low) && (quantiles[i][j].score < high)) {
474 temp.push_back(quantiles[i][j]);
482 //find contributer with most offending score related to it
483 int largestContrib = findLargestContrib(seen);
485 //while you still have guys to eliminate
486 while (contributions.size() > 0) {
488 mothurOut("Removing scores contributed by sequence " + toString(largestContrib) + " in your template file."); mothurOutEndLine();
490 //remove from quantiles all scores that were made using this largestContrib
491 for (int i = 0; i < quantiles.size(); i++) {
492 //cout << "before remove " << quantiles[i].size() << '\t';
493 removeContrib(largestContrib, quantiles[i]);
494 //cout << "after remove " << quantiles[i].size() << endl;
496 //cout << count << " largest contrib = " << largestContrib << endl; count++;
497 //remove from contributions all scores that were made using this largestContrib
498 removeContrib(largestContrib, contributions);
500 //"erase" largestContrib
501 seen[largestContrib] = -1;
503 //get next largestContrib
504 largestContrib = findLargestContrib(seen);
507 **************************************************************************************************
509 vector<int> marked = returnObviousOutliers(quantiles, num);
510 vector<int> copyMarked = marked;
512 //find the 99% of marked
513 sort(copyMarked.begin(), copyMarked.end());
514 int high = copyMarked[int(copyMarked.size() * 0.99)];
515 cout << "high = " << high << endl;
517 for(int i = 0; i < marked.size(); i++) {
518 if (marked[i] > high) {
519 mothurOut("Removing scores contributed by sequence " + toString(marked[i]) + " in your template file."); mothurOutEndLine();
520 for (int i = 0; i < quantiles.size(); i++) {
521 removeContrib(marked[i], quantiles[i]);
529 for (int i = 0; i < quantiles.size(); i++) {
532 if (quantiles[i].size() == 0) {
533 //in case this is not a distance found in your template files
534 for (int g = 0; g < 6; g++) {
539 sort(quantiles[i].begin(), quantiles[i].end(), compareQuanMembers);
542 temp.push_back(quantiles[i][int(quantiles[i].size() * 0.10)].score);
544 temp.push_back(quantiles[i][int(quantiles[i].size() * 0.25)].score);
546 temp.push_back(quantiles[i][int(quantiles[i].size() * 0.5)].score);
548 temp.push_back(quantiles[i][int(quantiles[i].size() * 0.75)].score);
550 temp.push_back(quantiles[i][int(quantiles[i].size() * 0.95)].score);
552 temp.push_back(quantiles[i][int(quantiles[i].size() * 0.99)].score);
562 catch(exception& e) {
563 errorOut(e, "DeCalculator", "removeObviousOutliers");
567 //***************************************************************************************************************
568 //put quanMember in the vector based on how far they are from the 99% or 1%. Biggest offenders in front.
569 /*vector<quanMember> DeCalculator::sortContrib(map<quanMember*, float> quan) {
572 vector<quanMember> newQuan;
573 map<quanMember*, float>::iterator it;
575 while (quan.size() > 0) {
577 map<quanMember*, float>::iterator largest = quan.begin();
579 //find biggest member
580 for (it = quan.begin(); it != quan.end(); it++) {
581 if (it->second > largest->second) { largest = it; }
583 cout << largest->second << '\t' << largest->first->score << '\t' << largest->first->member1 << '\t' << largest->first->member2 << endl;
585 newQuan.push_back(*(largest->first));
594 catch(exception& e) {
595 errorOut(e, "DeCalculator", "sortContrib");
600 //***************************************************************************************************************
601 //used by removeObviousOutliers which was attempt to increase sensitivity of chimera detection...not currently used...
602 int DeCalculator::findLargestContrib(vector<int> seen) {
609 for (int i = 0; i < seen.size(); i++) {
611 if (seen[i] > largest) {
617 return largestContribs;
620 catch(exception& e) {
621 errorOut(e, "DeCalculator", "findLargestContribs");
625 //***************************************************************************************************************
626 void DeCalculator::removeContrib(int bad, vector<quanMember>& quan) {
629 vector<quanMember> newQuan;
630 for (int i = 0; i < quan.size(); i++) {
631 if ((quan[i].member1 != bad) && (quan[i].member2 != bad) ) {
632 newQuan.push_back(quan[i]);
639 catch(exception& e) {
640 errorOut(e, "DeCalculator", "removeContrib");
645 //***************************************************************************************************************
646 float DeCalculator::findAverage(vector<float> myVector) {
650 for (int i = 0; i < myVector.size(); i++) { total += myVector[i]; }
652 float average = total / (float) myVector.size();
657 catch(exception& e) {
658 errorOut(e, "DeCalculator", "findAverage");
663 //***************************************************************************************************************
664 float DeCalculator::getCoef(vector<float> obs, vector<float> qav) {
667 //find average prob for this seqs windows
668 float probAverage = findAverage(qav);
670 //find observed average
671 float obsAverage = findAverage(obs);
673 float coef = obsAverage / probAverage;
677 catch(exception& e) {
678 errorOut(e, "DeCalculator", "getCoef");
682 //***************************************************************************************************************
683 //gets closest matches to each end, since chimeras will most likely have different parents on each end
684 vector<Sequence*> DeCalculator::findClosest(Sequence* querySeq, vector<Sequence*> db, int& numWanted, vector<int>& indexes) {
688 vector<Sequence*> seqsMatches;
689 vector<SeqDist> distsLeft;
690 vector<SeqDist> distsRight;
692 Dist* distcalculator = new eachGapDist();
694 string queryUnAligned = querySeq->getUnaligned();
695 int numBases = int(queryUnAligned.length() * 0.33);
697 string leftQuery = ""; //first 1/3 of the sequence
698 string rightQuery = ""; //last 1/3 of the sequence
699 string queryAligned = querySeq->getAligned();
702 bool foundFirstBase = false;
705 int firstBaseSpot = 0;
706 for (int i = 0; i < queryAligned.length(); i++) {
708 if (isalpha(queryAligned[i])) {
710 if (!foundFirstBase) { foundFirstBase = true; firstBaseSpot = i; }
713 //eliminate opening .'s
714 if (foundFirstBase) { leftQuery += queryAligned[i]; }
716 if (baseCount >= numBases) { leftSpot = i; break; } //first 1/3
719 //right side - count through another 1/3, so you are at last third
722 for (int i = leftSpot; i < queryAligned.length(); i++) {
724 if (isalpha(queryAligned[i])) { baseCount++; }
726 if (baseCount >= numBases) { rightSpot = i; break; } //last 1/3
730 //find last position in query that is a non gap character
731 int lastBaseSpot = queryAligned.length()-1;
732 for (int j = queryAligned.length()-1; j >= 0; j--) {
733 if (isalpha(queryAligned[j])) {
738 rightQuery = queryAligned.substr(rightSpot, (lastBaseSpot-rightSpot)); //sequence from pos spot to end
740 Sequence queryLeft(querySeq->getName(), leftQuery);
741 Sequence queryRight(querySeq->getName(), rightQuery);
742 //cout << querySeq->getName() << '\t' << leftSpot << '\t' << rightSpot << '\t' << firstBaseSpot << '\t' << lastBaseSpot << endl;
743 //cout << queryUnAligned.length() << '\t' << queryLeft.getUnaligned().length() << '\t' << queryRight.getUnaligned().length() << endl;
744 for(int j = 0; j < db.size(); j++){
746 string dbAligned = db[j]->getAligned();
747 string leftDB = dbAligned.substr(firstBaseSpot, (leftSpot-firstBaseSpot+1)); //first 1/3 of the sequence
748 string rightDB = dbAligned.substr(rightSpot, (lastBaseSpot-rightSpot)); //last 1/3 of the sequence
750 Sequence dbLeft(db[j]->getName(), leftDB);
751 Sequence dbRight(db[j]->getName(), rightDB);
753 distcalculator->calcDist(queryLeft, dbLeft);
754 float distLeft = distcalculator->getDist();
756 distcalculator->calcDist(queryRight, dbRight);
757 float distRight = distcalculator->getDist();
760 subjectLeft.seq = db[j];
761 subjectLeft.dist = distLeft;
762 subjectLeft.index = j;
764 distsLeft.push_back(subjectLeft);
766 SeqDist subjectRight;
767 subjectRight.seq = db[j];
768 subjectRight.dist = distRight;
769 subjectRight.index = j;
771 distsRight.push_back(subjectRight);
775 delete distcalculator;
777 //sort by smallest distance
778 sort(distsRight.begin(), distsRight.end(), compareSeqDist);
779 sort(distsLeft.begin(), distsLeft.end(), compareSeqDist);
782 map<string, string> seen;
783 map<string, string>::iterator it;
785 vector<SeqDist> dists;
786 float lastRight = distsRight[0].dist;
787 float lastLeft = distsLeft[0].dist;
789 for (int i = 0; i < distsLeft.size(); i++) {
790 //add left if you havent already
791 it = seen.find(distsLeft[i].seq->getName());
792 if (it == seen.end()) {
793 dists.push_back(distsLeft[i]);
794 seen[distsLeft[i].seq->getName()] = distsLeft[i].seq->getName();
795 lastLeft = distsLeft[i].dist;
798 //add right if you havent already
799 it = seen.find(distsRight[i].seq->getName());
800 if (it == seen.end()) {
801 dists.push_back(distsRight[i]);
802 seen[distsRight[i].seq->getName()] = distsRight[i].seq->getName();
803 lastRight = distsRight[i].dist;
806 if (dists.size() > numWanted) { lasti = i; break; } //you have enough results
812 while (i < distsLeft.size()) {
813 if (distsLeft[i].dist == lastLeft) { dists.push_back(distsLeft[i]); numWanted++; }
819 while (i < distsRight.size()) {
820 if (distsRight[i].dist == lastRight) { dists.push_back(distsRight[i]); numWanted++; }
825 //cout << numWanted << endl;
826 for (int i = 0; i < numWanted; i++) {
827 //cout << dists[i].seq->getName() << '\t' << dists[i].dist << endl;
828 Sequence* temp = new Sequence(dists[i].seq->getName(), dists[i].seq->getAligned()); //have to make a copy so you can trim and filter without stepping on eachother.
829 seqsMatches.push_back(temp);
830 indexes.push_back(dists[i].index);
835 catch(exception& e) {
836 errorOut(e, "DeCalculator", "findClosest");
840 //***************************************************************************************************************
841 Sequence* DeCalculator::findClosest(Sequence* querySeq, vector<Sequence*> db) {
846 Dist* distcalculator = new eachGapDist();
848 int smallest = 1000000;
850 for(int j = 0; j < db.size(); j++){
852 distcalculator->calcDist(*querySeq, *db[j]);
853 float dist = distcalculator->getDist();
855 if (dist < smallest) {
861 delete distcalculator;
863 seqsMatch = new Sequence(db[index]->getName(), db[index]->getAligned()); //have to make a copy so you can trim and filter without stepping on eachother.
867 catch(exception& e) {
868 errorOut(e, "DeCalculator", "findClosest");
872 /***************************************************************************************************************/
873 map<int, int> DeCalculator::trimSeqs(Sequence* query, vector<Sequence*> topMatches) {
876 int frontPos = 0; //should contain first position in all seqs that is not a gap character
877 int rearPos = query->getAligned().length();
879 //********find first position in topMatches that is a non gap character***********//
880 //find first position all query seqs that is a non gap character
881 for (int i = 0; i < topMatches.size(); i++) {
883 string aligned = topMatches[i]->getAligned();
886 //find first spot in this seq
887 for (int j = 0; j < aligned.length(); j++) {
888 if (isalpha(aligned[j])) {
894 //save this spot if it is the farthest
895 if (pos > frontPos) { frontPos = pos; }
899 string aligned = query->getAligned();
902 //find first position in query that is a non gap character
903 for (int j = 0; j < aligned.length(); j++) {
904 if (isalpha(aligned[j])) {
910 //save this spot if it is the farthest
911 if (pos > frontPos) { frontPos = pos; }
914 //********find last position in topMatches that is a non gap character***********//
915 for (int i = 0; i < topMatches.size(); i++) {
917 string aligned = topMatches[i]->getAligned();
918 int pos = aligned.length();
920 //find first spot in this seq
921 for (int j = aligned.length()-1; j >= 0; j--) {
922 if (isalpha(aligned[j])) {
928 //save this spot if it is the farthest
929 if (pos < rearPos) { rearPos = pos; }
933 aligned = query->getAligned();
934 pos = aligned.length();
936 //find last position in query that is a non gap character
937 for (int j = aligned.length()-1; j >= 0; j--) {
938 if (isalpha(aligned[j])) {
944 //save this spot if it is the farthest
945 if (pos < rearPos) { rearPos = pos; }
947 //check to make sure that is not whole seq
948 if ((rearPos - frontPos - 1) <= 0) { mothurOut("Error, when I trim your sequences, the entire sequence is trimmed."); mothurOutEndLine(); exit(1); }
949 //cout << query->getName() << " front = " << frontPos << " rear = " << rearPos << endl;
951 string newAligned = query->getAligned();
952 newAligned = newAligned.substr(frontPos, (rearPos-frontPos+1));
953 query->setAligned(newAligned);
956 for (int i = 0; i < topMatches.size(); i++) {
957 newAligned = topMatches[i]->getAligned();
958 newAligned = newAligned.substr(frontPos, (rearPos-frontPos+1));
959 topMatches[i]->setAligned(newAligned);
962 map<int, int> trimmedPos;
964 for (int i = 0; i < newAligned.length(); i++) {
965 trimmedPos[i] = i+frontPos;
970 catch(exception& e) {
971 errorOut(e, "DeCalculator", "trimSequences");
976 //***************************************************************************************************************