X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=decalc.cpp;h=2e214bffe969181c9f83061edc69a491610fe8a8;hb=2f17ce987b15dca9a309116d9df5131a1d3a35ee;hp=74cd8652909d17f0cd0e4e899b1472c6eea2ef93;hpb=cad05a21b084833b07808c1586e755be48fe7e1a;p=mothur.git diff --git a/decalc.cpp b/decalc.cpp index 74cd865..2e214bf 100644 --- a/decalc.cpp +++ b/decalc.cpp @@ -150,7 +150,7 @@ vector DeCalculator::findWindows(Sequence* query, int front, int back, int */ //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 - for (int m = front; m < (back - size) ; m+=increment) { win.push_back(m); } + for (int i = front; i < (back - size) ; i+=increment) { win.push_back(i); } @@ -169,10 +169,10 @@ vector DeCalculator::calcObserved(Sequence* query, Sequence* subject, vec vector temp; //int gaps = 0; - for (int m = 0; m < window.size(); m++) { + for (int i = 0; i < window.size(); i++) { - string seqFrag = query->getAligned().substr(window[m], size); - string seqFragsub = subject->getAligned().substr(window[m], size); + string seqFrag = query->getAligned().substr(window[i], size); + string seqFragsub = subject->getAligned().substr(window[i], size); int diff = 0; for (int b = 0; b < seqFrag.length(); b++) { @@ -244,9 +244,9 @@ vector DeCalculator::calcExpected(vector qav, float coef) { //for each window vector queryExpected; - for (int m = 0; m < qav.size(); m++) { + for (int j = 0; j < qav.size(); j++) { - float expected = qav[m] * coef; + float expected = qav[j] * coef; queryExpected.push_back(expected); } @@ -264,12 +264,12 @@ float DeCalculator::calcDE(vector obs, vector exp) { try { //for each window - float sum = 0.0; //sum = sum from 1 to m of (oi-ei)^2 + float sum = 0.0; //sum = sum from 1 to i of (oi-ei)^2 int numZeros = 0; - for (int m = 0; m < obs.size(); m++) { + for (int j = 0; j < obs.size(); j++) { - //if (obs[m] != 0.0) { - sum += ((obs[m] - exp[m]) * (obs[m] - exp[m])); + //if (obs[j] != 0.0) { + sum += ((obs[j] - exp[j]) * (obs[j] - exp[j])); //}else { numZeros++; } } @@ -321,7 +321,7 @@ vector DeCalculator::calcFreq(vector seqs, string filename) { //find base with highest frequency int highest = 0; - for (int m = 0; m < freq.size(); m++) { if (freq[m] > highest) { highest = freq[m]; } } + for (int j = 0; j < freq.size(); j++) { if (freq[j] > highest) { highest = freq[j]; } } float highFreq = highest / (float) (seqs.size()); @@ -355,13 +355,13 @@ vector DeCalculator::findQav(vector window, int size, vector vector averages; //for each window find average - for (int m = 0; m < window.size(); m++) { + for (int i = 0; i < window.size(); i++) { float average = 0.0; //while you are in the window for this sequence int count = 0; - for (int j = window[m]; j < (window[m]+size); j++) { + for (int j = window[i]; j < (window[i]+size); j++) { average += probabilityProfile[j]; count++; } @@ -381,16 +381,13 @@ vector DeCalculator::findQav(vector window, int size, vector } //*************************************************************************************************************** //seqs have already been masked -vector< vector > DeCalculator::getQuantiles(vector seqs, vector windowSizesTemplate, int window, vector probProfile, int increment, int start, int end) { +vector< vector > DeCalculator::getQuantiles(vector seqs, vector windowSizesTemplate, int window, vector probProfile, int increment, int start, int end) { try { - vector< vector > quan; + vector< vector > quan; //percentage of mismatched pairs 1 to 100 quan.resize(100); -//ofstream o; -//string out = "getQuantiles.out"; -//openOutputFile(out, o); - + //for each sequence for(int i = start; i < end; i++){ @@ -431,15 +428,16 @@ vector< vector > DeCalculator::getQuantiles(vector seqs, //cout << i << '\t' << j << '\t' << dist << '\t' << de << endl; dist = ceil(dist); - quanMember newScore(de, i, j); - - quan[dist].push_back(newScore); + //quanMember newScore(de, i, j); + quan[dist].push_back(de); + delete subject; } delete query; } + return quan; @@ -456,27 +454,27 @@ 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 >& quantiles, int num) { +void DeCalculator::removeObviousOutliers(vector< vector >& 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); - - float high = quantiles[i][int(quantiles[i].size() * 0.99)].score; - float low = quantiles[i][int(quantiles[i].size() * 0.01)].score; - vector temp; + //find mean of this quantile score + sort(quantiles[i].begin(), quantiles[i].end()); - //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)) { - temp.push_back(quantiles[i][j]); + vector temp; + if (quantiles[i].size() != 0) { + 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] > low) && (quantiles[i][j] < high)) { + temp.push_back(quantiles[i][j]); + } } } - quantiles[i] = temp; }