#include "pintail.h"
#include "ignoregaps.h"
+#include "eachgapdist.h"
+//********************************************************************************************************************
+//sorts lowest to highest
+inline bool compareQuanMembers(quanMember left, quanMember right){
+ return (left.score < right.score);
+}
//***************************************************************************************************************
-Pintail::Pintail(string filename, string temp) { fastafile = filename; templateFile = temp; }
+Pintail::Pintail(string filename, string temp, string o) { fastafile = filename; templateFile = temp; outputDir = o; }
//***************************************************************************************************************
Pintail::~Pintail() {
for (int i = 0; i < querySeqs.size(); i++) {
int index = ceil(deviation[i]);
-
+
//is your DE value higher than the 95%
string chimera;
- if (DE[i] > quantiles[index][4]) { chimera = "Yes"; }
- else { chimera = "No"; }
+ if (quantiles[index][4] == 0.0) {
+ chimera = "Your template does not include sequences that provide quantile values at distance " + toString(index);
+ }else {
+ if (DE[i] > quantiles[index][4]) { chimera = "Yes"; }
+ else { chimera = "No"; }
+ }
out << querySeqs[i]->getName() << '\t' << "div: " << deviation[i] << "\tstDev: " << DE[i] << "\tchimera flag: " << chimera << endl;
if (chimera == "Yes") {
}
//***************************************************************************************************************
-void Pintail::getChimeras() {
+int Pintail::getChimeras() {
try {
//read in query sequences and subject sequences
int numSeqs = querySeqs.size();
+ if (unaligned) { mothurOut("Your sequences need to be aligned when you use the pintail method."); mothurOutEndLine(); return 1; }
+
obsDistance.resize(numSeqs);
expectedDistance.resize(numSeqs);
seqCoef.resize(numSeqs);
templateLines.push_back(new linePair(0, templateSeqs.size()));
#endif
- distcalculator = new ignoreGaps();
+ distcalculator = new eachGapDist();
decalc = new DeCalculator();
//if the user does enter a mask then you want to keep all the spots in the alignment
mothurOut("Done."); mothurOutEndLine();
}else { createProcessesPairs(); }
-
- for (int j = 0; j < bestfit.size(); j++) {
-
- //chops off beginning and end of sequences so they both start and end with a base
- ofstream out;
- string s = querySeqs[j]->getName();
-
- openOutputFile(s, out);
- out << ">" << querySeqs[j]->getName() << endl;
- out << querySeqs[j]->getAligned() << endl;
- out.close();
-
- string t =querySeqs[j]->getName() + ".ref";
- openOutputFile(t, out);
- out << ">" << bestfit[j]->getName() << endl;
- out << bestfit[j]->getAligned() << endl;
- out.close();
- }
+//string o = "closestmatch.eachgap.fasta";
+//ofstream out7;
+//openOutputFile(o, out7);
-
+//for (int i = 0; i < bestfit.size(); i++) {
+ //out7 << ">" << querySeqs[i]->getName() << "-"<< bestfit[i]->getName() << endl;
+ //out7 << bestfit[i]->getAligned() << endl;
+//}
+//out7.close();
//find P
mothurOut("Getting conservation... "); cout.flush();
if (consfile == "") {
mothurOut("Calculating probability of conservation for your template sequences. This can take a while... I will output the frequency of the highest base in each position to a .freq file so that you can input them using the conservation parameter next time you run this command. Providing the .freq file will improve speed. "); cout.flush();
- probabilityProfile = decalc->calcFreq(templateSeqs, templateFile);
+ probabilityProfile = decalc->calcFreq(templateSeqs, outputDir + getSimpleName(templateFile));
mothurOut("Done."); mothurOutEndLine();
}else { probabilityProfile = readFreq(); }
}
}
-
+
+ if (filter) {
+ vector<Sequence*> temp = templateSeqs;
+ for (int i = 0; i < querySeqs.size(); i++) { temp.push_back(querySeqs[i]); }
+
+ createFilter(temp);
+
+ runFilter(querySeqs);
+ runFilter(templateSeqs);
+ runFilter(bestfit);
+ }
+
+
if (processors == 1) {
for (int j = 0; j < bestfit.size(); j++) {
//quantiles are used to determine whether the de values found indicate a chimera
//if you have to calculate them, its time intensive because you are finding the de and deviation values for each
//combination of sequences in the template
- if (quanfile != "") { quantiles = readQuantiles(); }
+ if (quanfile != "") { quantiles = readQuantiles(); }
else {
mothurOut("Calculating quantiles for your template. This can take a while... I will output the quantiles to a .quan file that you can input them using the quantiles parameter next time you run this command. Providing the .quan file will dramatically improve speed. "); cout.flush();
if (processors == 1) {
- quantilesMembers = decalc->getQuantiles(templateSeqs, windowSizesTemplate, window, probabilityProfile, increment, 0, templateSeqs.size(), makeCompliant);
+ quantilesMembers = decalc->getQuantiles(templateSeqs, windowSizesTemplate, window, probabilityProfile, increment, 0, templateSeqs.size());
}else { createProcessesQuan(); }
+
-
- //decided against this because we were having trouble setting the sensitivity... may want to revisit this...
- //quantiles = decalc->removeObviousOutliers(quantilesMembers, templateSeqs.size());
-
- ofstream out4;
- string o;
+ ofstream out4, out5;
+ string noOutliers, outliers;
- o = getRootName(templateFile) + "quan";
+ if ((!filter) && (seqMask == "")) {
+ noOutliers = outputDir + getRootName(getSimpleName(templateFile)) + "pintail.quan";
+ }else if ((filter) && (seqMask == "")) {
+ noOutliers = outputDir + getRootName(getSimpleName(templateFile)) + "pintail.filtered.quan";
+ }else if ((!filter) && (seqMask != "")) {
+ noOutliers = outputDir + getRootName(getSimpleName(templateFile)) + "pintail.masked.quan";
+ }else if ((filter) && (seqMask != "")) {
+ noOutliers = outputDir + getRootName(getSimpleName(templateFile)) + "pintail.filtered.masked.quan";
+ }
+
+ //outliers = getRootName(templateFile) + "pintail.quanYESOUTLIERS";
- openOutputFile(o, out4);
+ /*openOutputFile(outliers, out4);
//adjust quantiles
- for (int i = 0; i < quantiles.size(); i++) {
+ for (int i = 0; i < quantilesMembers.size(); i++) {
vector<float> temp;
- if (quantiles[i].size() == 0) {
+ if (quantilesMembers[i].size() == 0) {
//in case this is not a distance found in your template files
for (int g = 0; g < 6; g++) {
temp.push_back(0.0);
}
}else{
- sort(quantiles[i].begin(), quantiles[i].end());
+ sort(quantilesMembers[i].begin(), quantilesMembers[i].end(), compareQuanMembers);
//save 10%
- temp.push_back(quantiles[i][int(quantiles[i].size() * 0.10)]);
+ temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.10)].score);
//save 25%
- temp.push_back(quantiles[i][int(quantiles[i].size() * 0.25)]);
+ temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.25)].score);
//save 50%
- temp.push_back(quantiles[i][int(quantiles[i].size() * 0.5)]);
+ temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.5)].score);
//save 75%
- temp.push_back(quantiles[i][int(quantiles[i].size() * 0.75)]);
+ temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.75)].score);
//save 95%
- temp.push_back(quantiles[i][int(quantiles[i].size() * 0.95)]);
+ temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.95)].score);
//save 99%
- temp.push_back(quantiles[i][int(quantiles[i].size() * 0.99)]);
+ temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.99)].score);
}
}
- out4.close();
+ out4.close();*/
+
+ decalc->removeObviousOutliers(quantilesMembers, templateSeqs.size());
+ openOutputFile(noOutliers, out5);
+
+ //adjust quantiles
+ for (int i = 0; i < quantilesMembers.size(); i++) {
+ vector<float> temp;
+
+ if (quantilesMembers[i].size() == 0) {
+ //in case this is not a distance found in your template files
+ for (int g = 0; g < 6; g++) {
+ temp.push_back(0.0);
+ }
+ }else{
+
+ sort(quantilesMembers[i].begin(), quantilesMembers[i].end(), compareQuanMembers);
+
+ //save 10%
+ temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.10)].score);
+ //save 25%
+ temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.25)].score);
+ //save 50%
+ temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.5)].score);
+ //save 75%
+ temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.75)].score);
+ //save 95%
+ temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.95)].score);
+ //save 99%
+ temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.99)].score);
+
+ }
+
+ //output quan value
+ out5 << i+1 << '\t';
+ for (int u = 0; u < temp.size(); u++) { out5 << temp[u] << '\t'; }
+ out5 << endl;
+
+ quantiles[i] = temp;
+
+ }
+
mothurOut("Done."); mothurOutEndLine();
}
delete distcalculator;
delete decalc;
+
+ return 0;
}
catch(exception& e) {
errorOut(e, "Pintail", "getChimeras");
}
}
-//***************************************************************************************************************
-
-vector< vector<float> > Pintail::readQuantiles() {
- try {
-
- ifstream in;
- openInputFile(quanfile, in);
-
- vector< vector<float> > quan;
-
- int num; float ten, twentyfive, fifty, seventyfive, ninetyfive, ninetynine;
-
- while(!in.eof()){
-
- in >> num >> ten >> twentyfive >> fifty >> seventyfive >> ninetyfive >> ninetynine;
-
- vector <float> temp;
-
- temp.push_back(ten);
- temp.push_back(twentyfive);
- temp.push_back(fifty);
- temp.push_back(seventyfive);
- temp.push_back(ninetyfive);
- temp.push_back(ninetynine);
-
- quan.push_back(temp);
-
- gobble(in);
- }
-
- in.close();
- return quan;
-
- }
- catch(exception& e) {
- errorOut(e, "Pintail", "readQuantiles");
- exit(1);
- }
-}
//***************************************************************************************************************
//calculate the distances from each query sequence to all sequences in the template to find the closest sequence
vector<Sequence*> Pintail::findPairs(int start, int end) {
vector<Sequence*> seqsMatches;
for(int i = start; i < end; i++){
-
- float smallest = 10000.0;
- Sequence query = *(querySeqs[i]);
- Sequence* match;
-
- for(int j = 0; j < templateSeqs.size(); j++){
-
- Sequence temp = *(templateSeqs[j]);
-
- distcalculator->calcDist(query, temp);
- float dist = distcalculator->getDist();
-
- if (dist < smallest) {
- match = templateSeqs[j];
- smallest = dist;
- }
- }
-
- //make copy so trimSeqs doesn't overtrim
- Sequence* copy = new Sequence(match->getName(), match->getAligned());
- seqsMatches.push_back(copy);
+ vector<Sequence*> copy = decalc->findClosest(querySeqs[i], templateSeqs, 1);
+ seqsMatches.push_back(copy[0]);
}
return seqsMatches;
process++;
}else if (pid == 0){
- quantilesMembers = decalc->getQuantiles(templateSeqs, windowSizesTemplate, window, probabilityProfile, increment, templateLines[process]->start, templateLines[process]->end, makeCompliant);
+ quantilesMembers = decalc->getQuantiles(templateSeqs, windowSizesTemplate, window, probabilityProfile, increment, templateLines[process]->start, templateLines[process]->end);
//write out data to file so parent can read it
ofstream out;
}
#else
- quantilesMembers = decalc->getQuantiles(templateSeqs, windowSizesTemplate, window, probabilityProfile, increment, 0, templateSeqs.size(), makeCompliant);
+ quantilesMembers = decalc->getQuantiles(templateSeqs, windowSizesTemplate, window, probabilityProfile, increment, 0, templateSeqs.size());
#endif
}
catch(exception& e) {