From: westcott Date: Wed, 22 Jul 2009 18:58:12 +0000 (+0000) Subject: worked on chimeras X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=commitdiff_plain;h=dbd5da8043df1cb9f5ff7c6ddb5f550ea49b52c2 worked on chimeras --- diff --git a/Mothur.xcodeproj/project.pbxproj b/Mothur.xcodeproj/project.pbxproj index 9427eec..f1b5841 100644 --- a/Mothur.xcodeproj/project.pbxproj +++ b/Mothur.xcodeproj/project.pbxproj @@ -13,6 +13,7 @@ 21E859D80FC4632E005E1A48 /* matrixoutputcommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 21E859D70FC4632E005E1A48 /* matrixoutputcommand.cpp */; }; 370B88070F8A4EE4005AB382 /* getoturepcommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 370B88060F8A4EE4005AB382 /* getoturepcommand.cpp */; }; 371B30B40FD7EE67000414CA /* screenseqscommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 371B30B20FD7EE67000414CA /* screenseqscommand.cpp */; }; + 372A55551017922B00C5194B /* decalc.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 372A55541017922B00C5194B /* decalc.cpp */; }; 372E12700F26365B0095CF7E /* readotucommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 372E126F0F26365B0095CF7E /* readotucommand.cpp */; }; 372E12960F263D5A0095CF7E /* readdistcommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 372E12950F263D5A0095CF7E /* readdistcommand.cpp */; }; 372E12ED0F264D320095CF7E /* commandfactory.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 372E12EC0F264D320095CF7E /* commandfactory.cpp */; }; @@ -195,6 +196,8 @@ 370B88060F8A4EE4005AB382 /* getoturepcommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = getoturepcommand.cpp; sourceTree = SOURCE_ROOT; }; 371B30B20FD7EE67000414CA /* screenseqscommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = screenseqscommand.cpp; sourceTree = SOURCE_ROOT; }; 371B30B30FD7EE67000414CA /* screenseqscommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = screenseqscommand.h; sourceTree = SOURCE_ROOT; }; + 372A55531017922B00C5194B /* decalc.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = decalc.h; sourceTree = ""; }; + 372A55541017922B00C5194B /* decalc.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = decalc.cpp; sourceTree = ""; }; 372E126E0F26365B0095CF7E /* readotucommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = readotucommand.h; sourceTree = SOURCE_ROOT; }; 372E126F0F26365B0095CF7E /* readotucommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = readotucommand.cpp; sourceTree = SOURCE_ROOT; }; 372E12940F263D5A0095CF7E /* readdistcommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = readdistcommand.h; sourceTree = SOURCE_ROOT; }; @@ -639,6 +642,8 @@ 374F2FD51006320200E97C66 /* chimera.h */, 374F2FE9100634B600E97C66 /* bellerophon.h */, 374F2FEA100634B600E97C66 /* bellerophon.cpp */, + 372A55531017922B00C5194B /* decalc.h */, + 372A55541017922B00C5194B /* decalc.cpp */, 374F301110063B6F00E97C66 /* pintail.h */, 374F301210063B6F00E97C66 /* pintail.cpp */, ); @@ -1157,6 +1162,7 @@ 37B73CCF1004F5E0008C4B41 /* systemcommand.cpp in Sources */, 374F2FEB100634B600E97C66 /* bellerophon.cpp in Sources */, 374F301310063B6F00E97C66 /* pintail.cpp in Sources */, + 372A55551017922B00C5194B /* decalc.cpp in Sources */, ); runOnlyForDeploymentPostprocessing = 0; }; diff --git a/chimera.h b/chimera.h index 2febd00..2dc2989 100644 --- a/chimera.h +++ b/chimera.h @@ -13,6 +13,7 @@ #include "mothur.h" #include "sparsematrix.hpp" +#include "sequence.hpp" typedef list::iterator MatData; typedef map SeqMap; //maps sequence to all distance for that seqeunce @@ -37,6 +38,29 @@ class Chimera { virtual void setCons(string) {}; virtual void setQuantiles(string) {}; + virtual vector readSeqs(string file) { + try { + ifstream in; + openInputFile(file, in); + vector container; + + //read in seqs and store in vector + while(!in.eof()){ + + Sequence* current = new Sequence(in); + container.push_back(current); + gobble(in); + } + + in.close(); + return container; + } + catch(exception& e) { + errorOut(e, "Chimera", "readSeqs"); + exit(1); + } + } + //pure functions virtual void getChimeras() = 0; diff --git a/chimeraseqscommand.cpp b/chimeraseqscommand.cpp index 08c02a8..aaed3ba 100644 --- a/chimeraseqscommand.cpp +++ b/chimeraseqscommand.cpp @@ -22,7 +22,7 @@ ChimeraSeqsCommand::ChimeraSeqsCommand(string option){ else { //valid paramters for this command - string Array[] = {"fasta", "filter", "correction", "processors", "method", "window", "increment", "template", "conservation", "quantiles" }; + string Array[] = {"fasta", "filter", "correction", "processors", "method", "window", "increment", "template", "conservation", "quantile" }; vector myArray (Array, Array+(sizeof(Array)/sizeof(string))); OptionParser parser(option); @@ -48,9 +48,9 @@ ChimeraSeqsCommand::ChimeraSeqsCommand(string option){ if (consfile == "not open") { abort = true; } else if (consfile == "not found") { consfile = ""; } - quanfile = validParameter.validFile(parameters, "quantiles", true); + quanfile = validParameter.validFile(parameters, "quantile", true); if (quanfile == "not open") { abort = true; } - else if (quanfile == "not found") { consfile = ""; } + else if (quanfile == "not found") { quanfile = ""; } string temp; @@ -71,7 +71,7 @@ ChimeraSeqsCommand::ChimeraSeqsCommand(string option){ method = validParameter.validFile(parameters, "method", false); if (method == "not found") { method = "pintail"; } - if ((method == "pintail") && (templatefile == "") && (consfile == "")) { mothurOut("You must provide a template or conservation file with the pintail method."); mothurOutEndLine(); abort = true; } + if ((method == "pintail") && (templatefile == "")) { mothurOut("You must provide a template file with the pintail method."); mothurOutEndLine(); abort = true; } } diff --git a/decalc.cpp b/decalc.cpp new file mode 100644 index 0000000..7e4a494 --- /dev/null +++ b/decalc.cpp @@ -0,0 +1,11 @@ +/* + * decalc.cpp + * Mothur + * + * Created by Sarah Westcott on 7/22/09. + * Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved. + * + */ + +#include "decalc.h" + diff --git a/decalc.h b/decalc.h new file mode 100644 index 0000000..5818539 --- /dev/null +++ b/decalc.h @@ -0,0 +1,47 @@ +#ifndef DECALC_H +#define DECALC_H +/* + * decalc.h + * Mothur + * + * Created by Sarah Westcott on 7/22/09. + * Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved. + * + */ + +#include "mothur.h" +#include "sequence.hpp" + +/***********************************************************************/ + +//This class was created using the algorythms described in the +// "At Least 1 in 20 16S rRNA Sequence Records Currently Held in the Public Repositories is Estimated To Contain Substantial Anomalies" paper +//by Kevin E. Ashelford 1, Nadia A. Chuzhanova 3, John C. Fry 1, Antonia J. Jones 2 and Andrew J. Weightman 1. + +/***********************************************************************/ + +class DeCalculator { + + public: + + void trimSeqs(Sequence*, Sequence, map&); + vector readFreq(); + vector calcFreq(vector); + + vector findPairs(int, int); + vector findWindows(Sequence*, int, int, int&); + vector calcObserved(Sequence*, Sequence, vector, int); + vector calcExpected(vector, float); + vector findQav(vector, int); + float calcDE(vector, vector); + float calcDist(Sequence*, Sequence, int, int); + float getCoef(vector, vector); + + private: + + +}; + +/***********************************************************************/ + +#endif diff --git a/pintail.cpp b/pintail.cpp index 72be679..a204294 100644 --- a/pintail.cpp +++ b/pintail.cpp @@ -31,7 +31,14 @@ void Pintail::print(ostream& out) { for (int i = 0; i < querySeqs.size(); i++) { - out << querySeqs[i]->getName() << '\t' << "div: " << deviation[i] << "\tstDev: " << DE[i] << endl; + 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"; } + + out << querySeqs[i]->getName() << '\t' << "div: " << deviation[i] << "\tstDev: " << DE[i] << "\tchimera flag: " << chimera << endl; out << "Observed\t"; for (int j = 0; j < obsDistance[i].size(); j++) { out << obsDistance[i][j] << '\t'; } @@ -71,7 +78,9 @@ void Pintail::getChimeras() { deviation.resize(numSeqs); trimmed.resize(numSeqs); windowSizes.resize(numSeqs, window); - windows.resize(numSeqs); + windowSizesTemplate.resize(templateSeqs.size(), window); + windowsForeachQuery.resize(numSeqs); + quantiles.resize(100); //one for every percent mismatch //break up file if needed int linesPerProcess = processors / numSeqs; @@ -88,8 +97,19 @@ void Pintail::getChimeras() { int i = processors - 1; lines.push_back(new linePair((i*linesPerProcess), numSeqs)); } + + //find breakup of templatefile for quantiles + if (processors == 1) { templateLines.push_back(new linePair(0, templateSeqs.size())); } + else { + for (int i = 0; i < processors; i++) { + templateLines.push_back(new linePair()); + templateLines[i]->start = int (sqrt(float(i)/float(processors)) * templateSeqs.size()); + templateLines[i]->end = int (sqrt(float(i+1)/float(processors)) * templateSeqs.size()); + } + } #else lines.push_back(new linePair(0, numSeqs)); + templateLines.push_back(new linePair(0, templateSeqs.size())); #endif distcalculator = new ignoreGaps(); @@ -98,6 +118,8 @@ void Pintail::getChimeras() { if (processors == 1) { mothurOut("Finding closest sequence in template to each sequence... "); cout.flush(); bestfit = findPairs(lines[0]->start, lines[0]->end); + + //ex.align matches from wigeon /*for (int m = 0; m < templateSeqs.size(); m++) { if (templateSeqs[m]->getName() == "159481") { bestfit[17] = *(templateSeqs[m]); } if (templateSeqs[m]->getName() == "100137") { bestfit[16] = *(templateSeqs[m]); } @@ -123,46 +145,52 @@ void Pintail::getChimeras() { for (int j = 0; j < bestfit.size(); j++) { //chops off beginning and end of sequences so they both start and end with a base - trimSeqs(querySeqs[j], bestfit[j], j); + trimSeqs(querySeqs[j], bestfit[j], trimmed[j]); } mothurOut("Done."); mothurOutEndLine(); + mothurOut("Finding window breaks... "); cout.flush(); for (int i = lines[0]->start; i < lines[0]->end; i++) { it = trimmed[i].begin(); +cout << "trimmed = " << it->first << '\t' << it->second << endl; vector win = findWindows(querySeqs[i], it->first, it->second, windowSizes[i]); - windows[i] = win; + windowsForeachQuery[i] = win; } - - - } else { createProcessesSpots(); } + mothurOut("Done."); mothurOutEndLine(); + + }else { createProcessesSpots(); } //find P - if (consfile == "") { probabilityProfile = calcFreq(templateSeqs); } - else { probabilityProfile = readFreq(); } + 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 quantiles to a .prob file so that you can input them using the conservation parameter next time you run this command. Providing the .prob file will dramatically improve speed. "); cout.flush(); + probabilityProfile = calcFreq(templateSeqs); + mothurOut("Done."); mothurOutEndLine(); + }else { probabilityProfile = readFreq(); } //make P into Q for (int i = 0; i < probabilityProfile.size(); i++) { probabilityProfile[i] = 1 - probabilityProfile[i]; } + mothurOut("Done."); mothurOutEndLine(); if (processors == 1) { mothurOut("Calculating observed distance... "); cout.flush(); for (int i = lines[0]->start; i < lines[0]->end; i++) { - vector obsi = calcObserved(querySeqs[i], bestfit[i], windows[i], windowSizes[i]); + cout << querySeqs[i]->getName() << '\t' << bestfit[i].getName() << " windows = " << windowsForeachQuery[i].size() << " size = " << windowSizes[i] << endl; + vector obsi = calcObserved(querySeqs[i], bestfit[i], windowsForeachQuery[i], windowSizes[i]); obsDistance[i] = obsi; } mothurOut("Done."); mothurOutEndLine(); - mothurOut("Finding variability... "); cout.flush(); for (int i = lines[0]->start; i < lines[0]->end; i++) { - vector q = findQav(windows[i], windowSizes[i]); + vector q = findQav(windowsForeachQuery[i], windowSizes[i]); Qav[i] = q; } mothurOut("Done."); mothurOutEndLine(); - mothurOut("Calculating alpha... "); cout.flush(); for (int i = lines[0]->start; i < lines[0]->end; i++) { float alpha = getCoef(obsDistance[i], Qav[i]); @@ -171,7 +199,6 @@ void Pintail::getChimeras() { mothurOut("Done."); mothurOutEndLine(); - mothurOut("Calculating expected distance... "); cout.flush(); for (int i = lines[0]->start; i < lines[0]->end; i++) { vector exp = calcExpected(Qav[i], seqCoef[i]); @@ -180,7 +207,6 @@ void Pintail::getChimeras() { mothurOut("Done."); mothurOutEndLine(); - mothurOut("Finding deviation... "); cout.flush(); for (int i = lines[0]->start; i < lines[0]->end; i++) { float de = calcDE(obsDistance[i], expectedDistance[i]); @@ -195,65 +221,75 @@ void Pintail::getChimeras() { } else { createProcesses(); } - delete distcalculator; //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 != "") { 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) { - quantiles = getQuantiles(lines[0]->start, lines[0]->end); + quantiles = getQuantiles(0, templateSeqs.size()); }else { createProcessesQuan(); } - } - - - //free memory - for (int i = 0; i < lines.size(); i++) { delete lines[i]; } - - } - catch(exception& e) { - errorOut(e, "Pintail", "getChimeras"); - exit(1); - } -} - -//*************************************************************************************************************** - -vector Pintail::readSeqs(string file) { - try { - - ifstream in; - openInputFile(file, in); - vector container; - - //read in seqs and store in vector - while(!in.eof()){ - Sequence* current = new Sequence(in); + ofstream out4; + string o = getRootName(templateFile) + "quan"; - //if (current->getAligned() == "") { current->setAligned(current->getUnaligned()); } - //takes out stuff is needed - //current->setUnaligned(current->getUnaligned()); + openOutputFile(o, out4); - container.push_back(current); + //adjust quantiles + for (int i = 0; i < quantiles.size(); i++) { + if (quantiles[i].size() == 0) { + //in case this is not a distance found in your template files + for (int g = 0; g < 6; g++) { + quantiles[i].push_back(0.0); + } + }else{ + + sort(quantiles[i].begin(), quantiles[i].end()); + + vector temp; + //save 10% + temp.push_back(quantiles[i][int(quantiles[i].size() * 0.10)]); + //save 25% + temp.push_back(quantiles[i][int(quantiles[i].size() * 0.25)]); + //save 50% + temp.push_back(quantiles[i][int(quantiles[i].size() * 0.5)]); + //save 75% + temp.push_back(quantiles[i][int(quantiles[i].size() * 0.75)]); + //save 95% + temp.push_back(quantiles[i][int(quantiles[i].size() * 0.95)]); + //save 99% + temp.push_back(quantiles[i][int(quantiles[i].size() * 0.99)]); + + quantiles[i] = temp; + } + + //output quan value + out4 << i+1 << '\t'; + for (int u = 0; u < quantiles[i].size(); u++) { out4 << quantiles[i][u] << '\t'; } + out4 << endl; + + } - gobble(in); + mothurOut("Done."); mothurOutEndLine(); } - - in.close(); - return container; - + + //free memory + for (int i = 0; i < lines.size(); i++) { delete lines[i]; } + for (int i = 0; i < templateLines.size(); i++) { delete templateLines[i]; } + + delete distcalculator; } catch(exception& e) { - errorOut(e, "Pintail", "readSeqs"); + errorOut(e, "Pintail", "getChimeras"); exit(1); } } - //*************************************************************************************************************** //num is query's spot in querySeqs -map Pintail::trimSeqs(Sequence* query, Sequence subject, int num) { +void Pintail::trimSeqs(Sequence* query, Sequence subject, map& trim) { try { string q = query->getAligned(); @@ -269,15 +305,7 @@ map Pintail::trimSeqs(Sequence* query, Sequence subject, int num) { if (isalpha(q[i]) && isalpha(s[i])) { back = i; break; } } - //if num = -1 then you are calling this from quantiles - if (num != -1) { - trimmed[num][front] = back; - return trimmed[num]; - } - - map temp; - temp[front] = back; - return temp; + trim[front] = back; } catch(exception& e) { @@ -328,8 +356,7 @@ vector< vector > Pintail::readQuantiles() { openInputFile(quanfile, in); vector< vector > quan; - - //read in probabilities and store in vector + int num; float ten, twentyfive, fifty, seventyfive, ninetyfive, ninetynine; while(!in.eof()){ @@ -345,7 +372,6 @@ vector< vector > Pintail::readQuantiles() { temp.push_back(ninetyfive); temp.push_back(ninetynine); - //do you want this spot quan.push_back(temp); gobble(in); @@ -360,8 +386,6 @@ vector< vector > Pintail::readQuantiles() { exit(1); } } - - //*************************************************************************************************************** //calculate the distances from each query sequence to all sequences in the template to find the closest sequence vector Pintail::findPairs(int start, int end) { @@ -457,13 +481,12 @@ vector Pintail::calcObserved(Sequence* query, Sequence subject, vector temp; - - int startpoint = 0; - for (int m = 0; m < windows.size(); m++) { +//cout << "query length = " << query->getAligned().length() << '\t' << " subject length = " << subject.getAligned().length() << endl; + for (int m = 0; m < window.size(); m++) { - string seqFrag = query->getAligned().substr(window[startpoint], size); - string seqFragsub = subject.getAligned().substr(window[startpoint], size); - + string seqFrag = query->getAligned().substr(window[m], size); + string seqFragsub = subject.getAligned().substr(window[m], size); + //cout << "start point = " << window[m] << " end point = " << window[m]+size << endl; int diff = 0; for (int b = 0; b < seqFrag.length(); b++) { if (seqFrag[b] != seqFragsub[b]) { diff++; } @@ -474,8 +497,6 @@ vector Pintail::calcObserved(Sequence* query, Sequence subject, vector obs, vector exp) { //for each window float sum = 0.0; //sum = sum from 1 to m of (oi-ei)^2 - for (int m = 0; m < obsDistance.size(); m++) { sum += ((obs[m] - exp[m]) * (obs[m] - exp[m])); } + for (int m = 0; m < obs.size(); m++) { sum += ((obs[m] - exp[m]) * (obs[m] - exp[m])); } - float de = sqrt((sum / (obsDistance.size() - 1))); + float de = sqrt((sum / (obs.size() - 1))); return de; } @@ -558,7 +579,7 @@ vector Pintail::calcFreq(vector seqs) { try { vector prob; - string freqfile = getRootName(templateFile) + "probability"; + string freqfile = getRootName(templateFile) + "prob"; ofstream outFreq; openOutputFile(freqfile, outFreq); @@ -586,11 +607,10 @@ vector Pintail::calcFreq(vector seqs) { for (int m = 0; m < freq.size(); m++) { if (freq[m] > highest) { highest = freq[m]; } } float highFreq; + //subtract gaps to "ignore them" if ( (seqs.size() - gaps) == 0 ) { highFreq = 1.0; } else { highFreq = highest / (float) (seqs.size() - gaps); } - //highFreq = highest / (float) seqs.size(); - //cout << i << '\t' << highFreq << endl; - + float Pi; Pi = (highFreq - 0.25) / 0.75; @@ -619,7 +639,7 @@ vector Pintail::findQav(vector window, int size) { vector averages; //for each window find average - for (int m = 0; m < windows.size(); m++) { + for (int m = 0; m < window.size(); m++) { float average = 0.0; @@ -645,9 +665,14 @@ vector< vector > Pintail::getQuantiles(int start, int end) { try { vector< vector > quan; + //percentage of mismatched pairs 1 to 100 + quan.resize(100); + + //for each sequence for(int i = start; i < end; i++){ + mothurOut("Processing template sequence " + toString(i)); mothurOutEndLine(); Sequence* query = templateSeqs[i]; //compare to every other sequence in template @@ -655,17 +680,37 @@ vector< vector > Pintail::getQuantiles(int start, int end) { Sequence subject = *(templateSeqs[j]); - map trim = trimSeqs(query, subject, -1); + map trim; + trimSeqs(query, subject, trim); + it = trim.begin(); + int front = it->first; int back = it->second; + //reset window for each new comparison + windowSizesTemplate[i] = window; + + vector win = findWindows(query, front, back, windowSizesTemplate[i]); + + vector obsi = calcObserved(query, subject, win, windowSizesTemplate[i]); + + vector q = findQav(win, windowSizesTemplate[i]); + + float alpha = getCoef(obsi, q); + + vector exp = calcExpected(q, alpha); + + float de = calcDE(obsi, exp); + + float dist = calcDist(query, subject, front, back); + + dist = ceil(dist); + + //dist-1 because vector indexes start at 0. + quan[dist-1].push_back(de); - - } - - - } + return quan; } @@ -682,7 +727,7 @@ float Pintail::getCoef(vector obs, vector qav) { //find average prob for this seqs windows float probAverage = 0.0; for (int j = 0; j < qav.size(); j++) { probAverage += qav[j]; } - probAverage = probAverage / (float) Qav.size(); + probAverage = probAverage / (float) qav.size(); //find observed average float obsAverage = 0.0; @@ -717,6 +762,7 @@ void Pintail::createProcessesSpots() { process++; }else if (pid == 0){ + mothurOut("Finding pairs for sequences " + toString(lines[process]->start) + " to " + toString(lines[process]->end)); mothurOutEndLine(); vector tempbest; tempbest = findPairs(lines[process]->start, lines[process]->end); int count = 0; @@ -724,15 +770,17 @@ void Pintail::createProcessesSpots() { bestfit[i] = tempbest[count]; //chops off beginning and end of sequences so they both start and end with a base - trimSeqs(querySeqs[i], bestfit[i], i); + trimSeqs(querySeqs[i], bestfit[i], trimmed[i]); count++; } + mothurOut("Done finding pairs for sequences " + toString(lines[process]->start) + " to " + toString(lines[process]->end)); mothurOutEndLine(); - + mothurOut("Finding window breaks for sequences " + toString(lines[process]->start) + " to " + toString(lines[process]->end)); mothurOutEndLine(); for (int i = lines[process]->start; i < lines[process]->end; i++) { vector temp = findWindows(querySeqs[i], it->first, it->second, windowSizes[i]); win[i] = temp; } + mothurOut("Done finding window breaks for sequences " + toString(lines[process]->start) + " to " + toString(lines[process]->end)); mothurOutEndLine(); exit(0); }else { mothurOut("unable to spawn the necessary processes."); mothurOutEndLine(); exit(0); } @@ -744,7 +792,7 @@ void Pintail::createProcessesSpots() { wait(&temp); } - windows = win; + windowsForeachQuery = win; #else bestfit = findPairs(lines[0]->start, lines[0]->end); for (int j = 0; j < bestfit.size(); j++) { @@ -790,13 +838,14 @@ void Pintail::createProcesses() { process++; }else if (pid == 0){ - //calc obs + mothurOut("Calculating observed, expected and de values for sequences " + toString(lines[process]->start) + " to " + toString(lines[process]->end)); mothurOutEndLine(); for (int i = lines[process]->start; i < lines[process]->end; i++) { - vector obsi = calcObserved(querySeqs[i], bestfit[i], windows[i], windowSizes[i]); + + vector obsi = calcObserved(querySeqs[i], bestfit[i], windowsForeachQuery[i], windowSizes[i]); obs[i] = obsi; //calc Qav - vector q = findQav(windows[i], windowSizes[i]); + vector q = findQav(windowsForeachQuery[i], windowSizes[i]); //get alpha float alpha = getCoef(obsDistance[i], q); @@ -813,6 +862,7 @@ void Pintail::createProcesses() { float dist = calcDist(querySeqs[i], bestfit[i], it->first, it->second); dev[i] = dist; } + mothurOut("Done calculating observed, expected and de values for sequences " + toString(lines[process]->start) + " to " + toString(lines[process]->end)); mothurOutEndLine(); exit(0); }else { mothurOut("unable to spawn the necessary processes."); mothurOutEndLine(); exit(0); } @@ -893,6 +943,7 @@ void Pintail::createProcessesQuan() { #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) int process = 0; vector processIDS; + vector< vector > quan; quan.resize(100); //loop through and create all the processes you want while (process != processors) { @@ -903,7 +954,19 @@ void Pintail::createProcessesQuan() { process++; }else if (pid == 0){ + vector< vector > q = getQuantiles(templateLines[process]->start, templateLines[process]->end); + + for (int i = 0; i < q.size(); i++) { + //put all values of q[i] into quan[i] + quan[i].insert(quan[i].begin(), q[i].begin(), q[i].end()); + } + for (int i = 0; i < quan.size(); i++) { + cout << i+1 << '\t'; + for (int j = 0; j < quan[i].size(); j++) { cout << quan[i][j] << '\t'; } + cout << endl; + } + exit(0); }else { mothurOut("unable to spawn the necessary processes."); mothurOutEndLine(); exit(0); } } @@ -914,8 +977,9 @@ void Pintail::createProcessesQuan() { wait(&temp); } + quantiles = quan; #else - + quantiles = getQuantiles(0, templateSeqs.size()); #endif } catch(exception& e) { diff --git a/pintail.h b/pintail.h index 289af65..8c0dce4 100644 --- a/pintail.h +++ b/pintail.h @@ -38,6 +38,7 @@ class Pintail : public Chimera { int start; int end; linePair(int i, int j) : start(i), end(j) {} + linePair(){} }; Dist* distcalculator; @@ -45,6 +46,7 @@ class Pintail : public Chimera { string fastafile, templateFile, consfile, quanfile; vector lines; + vector templateLines; vector querySeqs; vector templateSeqs; @@ -53,10 +55,11 @@ class Pintail : public Chimera { vector< vector > obsDistance; //obsDistance[0] is the vector of observed distances for queryseqs[0]... vector< vector > expectedDistance; //expectedDistance[0] is the vector of expected distances for queryseqs[0]... vector deviation; //deviation[0] is the percentage of mismatched pairs over the whole seq between querySeqs[0] and its best match. - vector< vector > windows; // windows[0] is a vector containing the starting spot in queryseqs[0] aligned sequence for each window. + vector< vector > windowsForeachQuery; // windowsForeachQuery[0] is a vector containing the starting spot in queryseqs[0] aligned sequence for each window. //this is needed so you can move by bases and not just spots in the alignment - vector windowSizes; //windowSizes[0] = window size of querySeqs[0] + vector windowSizes; //windowSizes[0] = window size of querySeqs[0] + vector windowSizesTemplate; //windowSizesTemplate[0] = window size of templateSeqs[0] vector< map > trimmed; //trimmed[0] = start and stop of trimmed sequences for querySeqs[0] map::iterator it; @@ -67,14 +70,12 @@ class Pintail : public Chimera { vector probabilityProfile; vector< vector > quantiles; //quantiles[0] is the vector of deviations with ceiling score of 1, quantiles[1] is the vector of deviations with ceiling score of 2... - vector readSeqs(string); - map trimSeqs(Sequence*, Sequence, int); + void trimSeqs(Sequence*, Sequence, map&); vector readFreq(); vector< vector > readQuantiles(); vector< vector > getQuantiles(int, int); vector calcFreq(vector); - vector findPairs(int, int); vector findWindows(Sequence*, int, int, int&); vector calcObserved(Sequence*, Sequence, vector, int);