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 ms) {
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 m->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 m->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 m->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 m->mothurOut("You have selected too large a window size for sequence " + query->getName() + ". I will choose an appropriate window size."); m->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 m->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 m->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 m->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 m->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 m->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 m->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 m->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 m->mothurOut("Processing sequence " + toString(i)); m->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());
405 if (m->control_pressed) { delete query; delete subject; return quan; }
408 map<int, int>::iterator it;
410 trimSeqs(query, subject, trim);
413 int front = it->first; int back = it->second;
415 //reset window for each new comparison
416 windowSizesTemplate[i] = window;
418 vector<int> win = findWindows(query, front, back, windowSizesTemplate[i], increment);
420 vector<float> obsi = calcObserved(query, subject, win, windowSizesTemplate[i]);
422 vector<float> q = findQav(win, windowSizesTemplate[i], probProfile);
424 float alpha = getCoef(obsi, q);
426 vector<float> exp = calcExpected(q, alpha);
428 float de = calcDE(obsi, exp);
430 float dist = calcDist(query, subject, front, back);
431 //o << i << '\t' << j << '\t' << dist << '\t' << de << endl;
434 quanMember newScore(de, i, j);
436 quan[dist].push_back(newScore);
447 catch(exception& e) {
448 m->errorOut(e, "DeCalculator", "getQuantiles");
452 //********************************************************************************************************************
453 //sorts lowest to highest
454 inline bool compareQuanMembers(quanMember left, quanMember right){
455 return (left.score < right.score);
457 //***************************************************************************************************************
458 //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...
459 void DeCalculator::removeObviousOutliers(vector< vector<quanMember> >& quantiles, int num) {
462 for (int i = 0; i < quantiles.size(); i++) {
464 //find mean of this quantile score
465 sort(quantiles[i].begin(), quantiles[i].end(), compareQuanMembers);
467 float high = quantiles[i][int(quantiles[i].size() * 0.99)].score;
468 float low = quantiles[i][int(quantiles[i].size() * 0.01)].score;
470 vector<quanMember> temp;
472 //look at each value in quantiles to see if it is an outlier
473 for (int j = 0; j < quantiles[i].size(); j++) {
474 //is this score between 1 and 99%
475 if ((quantiles[i][j].score > low) && (quantiles[i][j].score < high)) {
476 temp.push_back(quantiles[i][j]);
484 //find contributer with most offending score related to it
485 int largestContrib = findLargestContrib(seen);
487 //while you still have guys to eliminate
488 while (contributions.size() > 0) {
490 m->mothurOut("Removing scores contributed by sequence " + toString(largestContrib) + " in your template file."); m->mothurOutEndLine();
492 //remove from quantiles all scores that were made using this largestContrib
493 for (int i = 0; i < quantiles.size(); i++) {
494 //cout << "before remove " << quantiles[i].size() << '\t';
495 removeContrib(largestContrib, quantiles[i]);
496 //cout << "after remove " << quantiles[i].size() << endl;
498 //cout << count << " largest contrib = " << largestContrib << endl; count++;
499 //remove from contributions all scores that were made using this largestContrib
500 removeContrib(largestContrib, contributions);
502 //"erase" largestContrib
503 seen[largestContrib] = -1;
505 //get next largestContrib
506 largestContrib = findLargestContrib(seen);
509 **************************************************************************************************
511 vector<int> marked = returnObviousOutliers(quantiles, num);
512 vector<int> copyMarked = marked;
514 //find the 99% of marked
515 sort(copyMarked.begin(), copyMarked.end());
516 int high = copyMarked[int(copyMarked.size() * 0.99)];
517 cout << "high = " << high << endl;
519 for(int i = 0; i < marked.size(); i++) {
520 if (marked[i] > high) {
521 m->mothurOut("Removing scores contributed by sequence " + toString(marked[i]) + " in your template file."); m->mothurOutEndLine();
522 for (int i = 0; i < quantiles.size(); i++) {
523 removeContrib(marked[i], quantiles[i]);
531 for (int i = 0; i < quantiles.size(); i++) {
534 if (quantiles[i].size() == 0) {
535 //in case this is not a distance found in your template files
536 for (int g = 0; g < 6; g++) {
541 sort(quantiles[i].begin(), quantiles[i].end(), compareQuanMembers);
544 temp.push_back(quantiles[i][int(quantiles[i].size() * 0.10)].score);
546 temp.push_back(quantiles[i][int(quantiles[i].size() * 0.25)].score);
548 temp.push_back(quantiles[i][int(quantiles[i].size() * 0.5)].score);
550 temp.push_back(quantiles[i][int(quantiles[i].size() * 0.75)].score);
552 temp.push_back(quantiles[i][int(quantiles[i].size() * 0.95)].score);
554 temp.push_back(quantiles[i][int(quantiles[i].size() * 0.99)].score);
564 catch(exception& e) {
565 m->errorOut(e, "DeCalculator", "removeObviousOutliers");
569 //***************************************************************************************************************
570 //put quanMember in the vector based on how far they are from the 99% or 1%. Biggest offenders in front.
571 /*vector<quanMember> DeCalculator::sortContrib(map<quanMember*, float> quan) {
574 vector<quanMember> newQuan;
575 map<quanMember*, float>::iterator it;
577 while (quan.size() > 0) {
579 map<quanMember*, float>::iterator largest = quan.begin();
581 //find biggest member
582 for (it = quan.begin(); it != quan.end(); it++) {
583 if (it->second > largest->second) { largest = it; }
585 cout << largest->second << '\t' << largest->first->score << '\t' << largest->first->member1 << '\t' << largest->first->member2 << endl;
587 newQuan.push_back(*(largest->first));
596 catch(exception& e) {
597 m->errorOut(e, "DeCalculator", "sortContrib");
602 //***************************************************************************************************************
603 //used by removeObviousOutliers which was attempt to increase sensitivity of chimera detection...not currently used...
604 int DeCalculator::findLargestContrib(vector<int> seen) {
611 for (int i = 0; i < seen.size(); i++) {
613 if (seen[i] > largest) {
619 return largestContribs;
622 catch(exception& e) {
623 m->errorOut(e, "DeCalculator", "findLargestContribs");
627 //***************************************************************************************************************
628 void DeCalculator::removeContrib(int bad, vector<quanMember>& quan) {
631 vector<quanMember> newQuan;
632 for (int i = 0; i < quan.size(); i++) {
633 if ((quan[i].member1 != bad) && (quan[i].member2 != bad) ) {
634 newQuan.push_back(quan[i]);
641 catch(exception& e) {
642 m->errorOut(e, "DeCalculator", "removeContrib");
647 //***************************************************************************************************************
648 float DeCalculator::findAverage(vector<float> myVector) {
652 for (int i = 0; i < myVector.size(); i++) { total += myVector[i]; }
654 float average = total / (float) myVector.size();
659 catch(exception& e) {
660 m->errorOut(e, "DeCalculator", "findAverage");
665 //***************************************************************************************************************
666 float DeCalculator::getCoef(vector<float> obs, vector<float> qav) {
669 //find average prob for this seqs windows
670 float probAverage = findAverage(qav);
672 //find observed average
673 float obsAverage = findAverage(obs);
675 float coef = obsAverage / probAverage;
679 catch(exception& e) {
680 m->errorOut(e, "DeCalculator", "getCoef");
684 //***************************************************************************************************************
685 //gets closest matches to each end, since chimeras will most likely have different parents on each end
686 vector<Sequence*> DeCalculator::findClosest(Sequence* querySeq, vector<Sequence*> db, int& numWanted, vector<int>& indexes) {
690 vector<Sequence*> seqsMatches;
691 vector<SeqDist> distsLeft;
692 vector<SeqDist> distsRight;
694 Dist* distcalculator = new eachGapDist();
696 string queryUnAligned = querySeq->getUnaligned();
697 int numBases = int(queryUnAligned.length() * 0.33);
699 string leftQuery = ""; //first 1/3 of the sequence
700 string rightQuery = ""; //last 1/3 of the sequence
701 string queryAligned = querySeq->getAligned();
704 bool foundFirstBase = false;
707 int firstBaseSpot = 0;
708 for (int i = 0; i < queryAligned.length(); i++) {
710 if (isalpha(queryAligned[i])) {
712 if (!foundFirstBase) { foundFirstBase = true; firstBaseSpot = i; }
715 //eliminate opening .'s
716 if (foundFirstBase) { leftQuery += queryAligned[i]; }
718 if (baseCount >= numBases) { leftSpot = i; break; } //first 1/3
721 //right side - count through another 1/3, so you are at last third
724 for (int i = leftSpot; i < queryAligned.length(); i++) {
726 if (isalpha(queryAligned[i])) { baseCount++; }
728 if (baseCount >= numBases) { rightSpot = i; break; } //last 1/3
732 //find last position in query that is a non gap character
733 int lastBaseSpot = queryAligned.length()-1;
734 for (int j = queryAligned.length()-1; j >= 0; j--) {
735 if (isalpha(queryAligned[j])) {
740 rightQuery = queryAligned.substr(rightSpot, (lastBaseSpot-rightSpot)); //sequence from pos spot to end
742 Sequence queryLeft(querySeq->getName(), leftQuery);
743 Sequence queryRight(querySeq->getName(), rightQuery);
744 //cout << querySeq->getName() << '\t' << leftSpot << '\t' << rightSpot << '\t' << firstBaseSpot << '\t' << lastBaseSpot << endl;
745 //cout << queryUnAligned.length() << '\t' << queryLeft.getUnaligned().length() << '\t' << queryRight.getUnaligned().length() << endl;
746 for(int j = 0; j < db.size(); j++){
748 string dbAligned = db[j]->getAligned();
749 string leftDB = dbAligned.substr(firstBaseSpot, (leftSpot-firstBaseSpot+1)); //first 1/3 of the sequence
750 string rightDB = dbAligned.substr(rightSpot, (lastBaseSpot-rightSpot)); //last 1/3 of the sequence
752 Sequence dbLeft(db[j]->getName(), leftDB);
753 Sequence dbRight(db[j]->getName(), rightDB);
755 distcalculator->calcDist(queryLeft, dbLeft);
756 float distLeft = distcalculator->getDist();
758 distcalculator->calcDist(queryRight, dbRight);
759 float distRight = distcalculator->getDist();
762 subjectLeft.seq = db[j];
763 subjectLeft.dist = distLeft;
764 subjectLeft.index = j;
766 distsLeft.push_back(subjectLeft);
768 SeqDist subjectRight;
769 subjectRight.seq = db[j];
770 subjectRight.dist = distRight;
771 subjectRight.index = j;
773 distsRight.push_back(subjectRight);
777 delete distcalculator;
779 //sort by smallest distance
780 sort(distsRight.begin(), distsRight.end(), compareSeqDist);
781 sort(distsLeft.begin(), distsLeft.end(), compareSeqDist);
784 map<string, string> seen;
785 map<string, string>::iterator it;
787 vector<SeqDist> dists;
788 float lastRight = distsRight[0].dist;
789 float lastLeft = distsLeft[0].dist;
791 for (int i = 0; i < distsLeft.size(); i++) {
792 //add left if you havent already
793 it = seen.find(distsLeft[i].seq->getName());
794 if (it == seen.end()) {
795 dists.push_back(distsLeft[i]);
796 seen[distsLeft[i].seq->getName()] = distsLeft[i].seq->getName();
797 lastLeft = distsLeft[i].dist;
800 //add right if you havent already
801 it = seen.find(distsRight[i].seq->getName());
802 if (it == seen.end()) {
803 dists.push_back(distsRight[i]);
804 seen[distsRight[i].seq->getName()] = distsRight[i].seq->getName();
805 lastRight = distsRight[i].dist;
808 if (dists.size() > numWanted) { lasti = i; break; } //you have enough results
814 while (i < distsLeft.size()) {
815 if (distsLeft[i].dist == lastLeft) { dists.push_back(distsLeft[i]); numWanted++; }
821 while (i < distsRight.size()) {
822 if (distsRight[i].dist == lastRight) { dists.push_back(distsRight[i]); numWanted++; }
827 //cout << numWanted << endl;
828 for (int i = 0; i < numWanted; i++) {
829 //cout << dists[i].seq->getName() << '\t' << dists[i].dist << endl;
830 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.
831 seqsMatches.push_back(temp);
832 indexes.push_back(dists[i].index);
837 catch(exception& e) {
838 m->errorOut(e, "DeCalculator", "findClosest");
842 //***************************************************************************************************************
843 Sequence* DeCalculator::findClosest(Sequence* querySeq, vector<Sequence*> db) {
848 Dist* distcalculator = new eachGapDist();
850 int smallest = 1000000;
852 for(int j = 0; j < db.size(); j++){
854 distcalculator->calcDist(*querySeq, *db[j]);
855 float dist = distcalculator->getDist();
857 if (dist < smallest) {
863 delete distcalculator;
865 seqsMatch = new Sequence(db[index]->getName(), db[index]->getAligned()); //have to make a copy so you can trim and filter without stepping on eachother.
869 catch(exception& e) {
870 m->errorOut(e, "DeCalculator", "findClosest");
874 /***************************************************************************************************************/
875 map<int, int> DeCalculator::trimSeqs(Sequence* query, vector<Sequence*> topMatches) {
878 int frontPos = 0; //should contain first position in all seqs that is not a gap character
879 int rearPos = query->getAligned().length();
881 //********find first position in topMatches that is a non gap character***********//
882 //find first position all query seqs that is a non gap character
883 for (int i = 0; i < topMatches.size(); i++) {
885 string aligned = topMatches[i]->getAligned();
888 //find first spot in this seq
889 for (int j = 0; j < aligned.length(); j++) {
890 if (isalpha(aligned[j])) {
896 //save this spot if it is the farthest
897 if (pos > frontPos) { frontPos = pos; }
901 string aligned = query->getAligned();
904 //find first position in query that is a non gap character
905 for (int j = 0; j < aligned.length(); j++) {
906 if (isalpha(aligned[j])) {
912 //save this spot if it is the farthest
913 if (pos > frontPos) { frontPos = pos; }
916 //********find last position in topMatches that is a non gap character***********//
917 for (int i = 0; i < topMatches.size(); i++) {
919 string aligned = topMatches[i]->getAligned();
920 int pos = aligned.length();
922 //find first spot in this seq
923 for (int j = aligned.length()-1; j >= 0; j--) {
924 if (isalpha(aligned[j])) {
930 //save this spot if it is the farthest
931 if (pos < rearPos) { rearPos = pos; }
935 aligned = query->getAligned();
936 pos = aligned.length();
938 //find last position in query that is a non gap character
939 for (int j = aligned.length()-1; j >= 0; j--) {
940 if (isalpha(aligned[j])) {
946 //save this spot if it is the farthest
947 if (pos < rearPos) { rearPos = pos; }
949 //check to make sure that is not whole seq
950 if ((rearPos - frontPos - 1) <= 0) { m->mothurOut("Error, when I trim your sequences, the entire sequence is trimmed."); m->mothurOutEndLine(); exit(1); }
951 //cout << query->getName() << " front = " << frontPos << " rear = " << rearPos << endl;
953 string newAligned = query->getAligned();
954 newAligned = newAligned.substr(frontPos, (rearPos-frontPos+1));
955 query->setAligned(newAligned);
958 for (int i = 0; i < topMatches.size(); i++) {
959 newAligned = topMatches[i]->getAligned();
960 newAligned = newAligned.substr(frontPos, (rearPos-frontPos+1));
961 topMatches[i]->setAligned(newAligned);
964 map<int, int> trimmedPos;
966 for (int i = 0; i < newAligned.length(); i++) {
967 trimmedPos[i] = i+frontPos;
972 catch(exception& e) {
973 m->errorOut(e, "DeCalculator", "trimSequences");
978 //***************************************************************************************************************