]> git.donarmstrong.com Git - mothur.git/blobdiff - decalc.cpp
added diffs and percent parameters to cluster.fragments command
[mothur.git] / decalc.cpp
index 4e80c65d29aa7c8d6a2354db95027c3365c7962d..8d90a78a492e63ee97ede0b577e37fde39336fb4 100644 (file)
@@ -290,10 +290,12 @@ vector<float> DeCalculator::calcFreq(vector<Sequence*> seqs, string filename) {
        try {
 
                vector<float> prob;
-               string freqfile = getRootName(filename) + "freq";
+               string freqfile = m->getRootName(filename) + "freq";
                ofstream outFreq;
                
-               openOutputFile(freqfile, outFreq);
+               m->openOutputFile(freqfile, outFreq);
+               
+               outFreq << "#" << m->getVersion() << endl;
                
                string length = toString(seqs.size());  //if there are 5000 seqs in the template then set precision to 3
                int precision = length.length() - 1;
@@ -381,9 +383,9 @@ vector<float>  DeCalculator::findQav(vector<int> window, int size, vector<float>
 }
 //***************************************************************************************************************
 //seqs have already been masked
-vector< vector<quanMember> > DeCalculator::getQuantiles(vector<Sequence*> seqs, vector<int> windowSizesTemplate, int window, vector<float> probProfile, int increment, int start, int end) {
+vector< vector<float> > DeCalculator::getQuantiles(vector<Sequence*> seqs, vector<int> windowSizesTemplate, int window, vector<float> probProfile, int increment, int start, int end) {
        try {
-               vector< vector<quanMember> > quan; 
+               vector< vector<float> > quan; 
                
                //percentage of mismatched pairs 1 to 100
                quan.resize(100);
@@ -428,9 +430,9 @@ vector< vector<quanMember> > DeCalculator::getQuantiles(vector<Sequence*> seqs,
        //cout << i << '\t' <<  j << '\t' << dist << '\t' << de << endl;                        
                                dist = ceil(dist);
                                
-                               quanMember newScore(de, i, j);
+                               //quanMember newScore(de, i, j);
                                
-                               quan[dist].push_back(newScore);
+                               quan[dist].push_back(de);
 
                                delete subject;
                        }
@@ -454,23 +456,23 @@ inline bool compareQuanMembers(quanMember left, quanMember right){
 } 
 //***************************************************************************************************************
 //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...
-void DeCalculator::removeObviousOutliers(vector< vector<quanMember> >& quantiles, int num) {
+void DeCalculator::removeObviousOutliers(vector< vector<float> >& quantiles, int num) {
        try {
                                                
                for (int i = 0; i < quantiles.size(); i++) {
                        
                        //find mean of this quantile score
-                       sort(quantiles[i].begin(), quantiles[i].end(), compareQuanMembers);
+                       sort(quantiles[i].begin(), quantiles[i].end());
                        
-                       vector<quanMember> temp;
+                       vector<float> temp;
                        if (quantiles[i].size() != 0) {
-                               float high = quantiles[i][int(quantiles[i].size() * 0.99)].score;
-                               float low =  quantiles[i][int(quantiles[i].size() * 0.01)].score;
+                               float high = quantiles[i][int(quantiles[i].size() * 0.99)];
+                               float low =  quantiles[i][int(quantiles[i].size() * 0.01)];
                        
                                //look at each value in quantiles to see if it is an outlier
                                for (int j = 0; j < quantiles[i].size(); j++) {
                                        //is this score between 1 and 99%
-                                       if ((quantiles[i][j].score > low) && (quantiles[i][j].score < high)) {
+                                       if ((quantiles[i][j] > low) && (quantiles[i][j] < high)) {
                                                temp.push_back(quantiles[i][j]);
                                        }
                                }
@@ -686,6 +688,7 @@ vector<Sequence*> DeCalculator::findClosest(Sequence* querySeq, vector<Sequence*
                indexes.clear();
                
                vector<Sequence*> seqsMatches;  
+               
                vector<SeqDist> distsLeft;
                vector<SeqDist> distsRight;
                
@@ -757,14 +760,14 @@ vector<Sequence*> DeCalculator::findClosest(Sequence* querySeq, vector<Sequence*
                        float distRight = distcalculator->getDist();
                        
                        SeqDist subjectLeft;
-                       subjectLeft.seq = db[j];
+                       subjectLeft.seq = NULL;
                        subjectLeft.dist = distLeft;
                        subjectLeft.index = j;
                        
                        distsLeft.push_back(subjectLeft);
                        
                        SeqDist subjectRight;
-                       subjectRight.seq = db[j];
+                       subjectRight.seq = NULL;
                        subjectRight.dist = distRight;
                        subjectRight.index = j;
                        
@@ -788,18 +791,18 @@ vector<Sequence*> DeCalculator::findClosest(Sequence* querySeq, vector<Sequence*
                int lasti = 0;
                for (int i = 0; i < distsLeft.size(); i++) {
                        //add left if you havent already
-                       it = seen.find(distsLeft[i].seq->getName());
+                       it = seen.find(db[distsLeft[i].index]->getName());
                        if (it == seen.end()) {  
                                dists.push_back(distsLeft[i]);
-                               seen[distsLeft[i].seq->getName()] = distsLeft[i].seq->getName();
+                               seen[db[distsLeft[i].index]->getName()] = db[distsLeft[i].index]->getName();
                                lastLeft =  distsLeft[i].dist;
                        }
 
                        //add right if you havent already
-                       it = seen.find(distsRight[i].seq->getName());
+                       it = seen.find(db[distsRight[i].index]->getName());
                        if (it == seen.end()) {  
                                dists.push_back(distsRight[i]);
-                               seen[distsRight[i].seq->getName()] = distsRight[i].seq->getName();
+                               seen[db[distsRight[i].index]->getName()] = db[distsRight[i].index]->getName();
                                lastRight =  distsRight[i].dist;
                        }
                        
@@ -827,7 +830,7 @@ vector<Sequence*> DeCalculator::findClosest(Sequence* querySeq, vector<Sequence*
 //cout << numWanted << endl;
                for (int i = 0; i < numWanted; i++) {
 //cout << dists[i].seq->getName() << '\t' << dists[i].dist << endl;
-                       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.
+                       Sequence* temp = new Sequence(db[dists[i].index]->getName(), db[dists[i].index]->getAligned()); //have to make a copy so you can trim and filter without stepping on eachother.
                        seqsMatches.push_back(temp);
                        indexes.push_back(dists[i].index);
                }
@@ -945,28 +948,35 @@ map<int, int> DeCalculator::trimSeqs(Sequence* query, vector<Sequence*> topMatch
                
                //save this spot if it is the farthest
                if (pos < rearPos) { rearPos = pos; }
-
-               //check to make sure that is not whole seq
-               if ((rearPos - frontPos - 1) <= 0) {  m->mothurOut("Error, when I trim your sequences, the entire sequence is trimmed."); m->mothurOutEndLine(); exit(1);  }
-//cout << query->getName() << " front = " << frontPos << " rear = " << rearPos << endl;                
-               //trim query
-               string newAligned = query->getAligned();
-               newAligned = newAligned.substr(frontPos, (rearPos-frontPos+1));
-               query->setAligned(newAligned);
-               
-               //trim topMatches
-               for (int i = 0; i < topMatches.size(); i++) {
-                       newAligned = topMatches[i]->getAligned();
-                       newAligned = newAligned.substr(frontPos, (rearPos-frontPos+1));
-                       topMatches[i]->setAligned(newAligned);
-               }
                
                map<int, int> trimmedPos;
-               
-               for (int i = 0; i < newAligned.length(); i++) {
-                       trimmedPos[i] = i+frontPos;
+               //check to make sure that is not whole seq
+               if ((rearPos - frontPos - 1) <= 0) {  
+                       m->mothurOut("[ERROR]: when I trim " + query->getName() + ", the entire sequence is trimmed. Skipping."); m->mothurOutEndLine();  
+                       query->setAligned("");
+                       //trim topMatches
+                       for (int i = 0; i < topMatches.size(); i++) {
+                               topMatches[i]->setAligned("");
+                       }
+                       
+               }else {
+
+                       //trim query
+                       string newAligned = query->getAligned();
+                       newAligned = newAligned.substr(frontPos, (rearPos-frontPos+1));
+                       query->setAligned(newAligned);
+                       
+                       //trim topMatches
+                       for (int i = 0; i < topMatches.size(); i++) {
+                               newAligned = topMatches[i]->getAligned();
+                               newAligned = newAligned.substr(frontPos, (rearPos-frontPos+1));
+                               topMatches[i]->setAligned(newAligned);
+                       }
+                       
+                       for (int i = 0; i < newAligned.length(); i++) {
+                               trimmedPos[i] = i+frontPos;
+                       }
                }
-               
                return trimmedPos;
        }
        catch(exception& e) {