#include "command.hpp"
#include "sequence.hpp"
#include "sequenceparser.h"
+#include "sequencecountparser.h"
/************************************************************/
struct seqPNode {
Sequence seq;
string names;
bool active;
+ int diffs;
seqPNode() {}
- seqPNode(int n, Sequence s, string nm) : numIdentical(n), seq(s), names(nm), active(1) {}
+ seqPNode(int n, Sequence s, string nm) : numIdentical(n), seq(s), names(nm), active(1) { diffs = 0; }
~seqPNode() {}
};
/************************************************************/
-inline bool comparePriority(seqPNode first, seqPNode second) { return (first.numIdentical > second.numIdentical); }
+inline bool comparePriorityTopDown(seqPNode first, seqPNode second) {
+ if (first.numIdentical > second.numIdentical) { return true; }
+ else if (first.numIdentical == second.numIdentical) {
+ if (first.seq.getName() > second.seq.getName()) { return true; }
+ }
+ return false;
+}
+/************************************************************/
+inline bool comparePriorityDownTop(seqPNode first, seqPNode second) {
+ if (first.numIdentical < second.numIdentical) { return true; }
+ else if (first.numIdentical == second.numIdentical) {
+ if (first.seq.getName() > second.seq.getName()) { return true; }
+ }
+ return false;
+}
//************************************************************/
class PreClusterCommand : public Command {
vector<string> setParameters();
string getCommandName() { return "pre.cluster"; }
string getCommandCategory() { return "Sequence Processing"; }
+
string getHelpString();
- string getCitation() { return "http://www.mothur.org/wiki/Pre.cluster"; }
+ string getOutputPattern(string);
+ string getCitation() { return "Schloss PD, Gevers D, Westcott SL (2011). Reducing the effects of PCR amplification and sequencing artifacts on 16S rRNA-based studies. PLoS ONE. 6:e27310.\nhttp://www.mothur.org/wiki/Pre.cluster"; }
string getDescription() { return "implements a pseudo-single linkage algorithm with the goal of removing sequences that are likely due to pyrosequencing errors"; }
linePair(int i, int j) : start(i), end(j) {}
};
+ SequenceParser* parser;
+ SequenceCountParser* cparser;
+ CountTable ct;
+
int diffs, length, processors;
- bool abort, bygroup;
- string fastafile, namefile, outputDir, groupfile;
+ bool abort, bygroup, topdown;
+ string fastafile, namefile, outputDir, groupfile, countfile;
vector<seqPNode> alignSeqs; //maps the number of identical seqs to a sequence
map<string, string> names; //represents the names file first column maps to second column
map<string, int> sizes; //this map a seq name to the number of identical seqs in the names file
map<string, int>::iterator itSize;
// map<string, bool> active; //maps sequence name to whether it has already been merged or not.
vector<string> outputNames;
- map<string, vector<string> > outputTypes;
int readFASTA();
void readNameFile();
//int readNamesFASTA();
int calcMisMatches(string, string);
- void printData(string, string); //fasta filename, names file name
- int process();
- int loadSeqs(map<string, string>&, vector<Sequence>&);
- int driverGroups(SequenceParser*, string, string, int, int, vector<string> groups);
- int createProcessesGroups(SequenceParser*, string, string, vector<string>);
+ void printData(string, string, string); //fasta filename, names file name
+ int process(string);
+ int loadSeqs(map<string, string>&, vector<Sequence>&, string);
+ int driverGroups(string, string, string, int, int, vector<string> groups);
+ int createProcessesGroups(string, string, string, vector<string>);
+ int mergeGroupCounts(string, string, string);
};
/**************************************************************************************************/
struct preClusterData {
string fastafile;
string namefile;
- string groupfile;
- string newFName, newNName;
+ string groupfile, countfile;
+ string newFName, newNName, newMName;
MothurOut* m;
int start;
int end;
int diffs, threadID;
vector<string> groups;
+ vector<string> mapFileNames;
+ bool topdown;
preClusterData(){}
- preClusterData(string f, string n, string g, string nff, string nnf, vector<string> gr, MothurOut* mout, int st, int en, int d, int tid) {
+ preClusterData(string f, string n, string g, string c, string nff, string nnf, string nmf, vector<string> gr, MothurOut* mout, int st, int en, int d, bool td, int tid) {
fastafile = f;
namefile = n;
groupfile = g;
newFName = nff;
newNName = nnf;
+ newMName = nmf;
m = mout;
start = st;
end = en;
diffs = d;
threadID = tid;
groups = gr;
+ countfile = c;
+ topdown = td;
}
};
/**************************************************************************************************/
-#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
#else
static DWORD WINAPI MyPreclusterThreadFunction(LPVOID lpParam){
preClusterData* pDataArray;
//parse fasta and name file by group
SequenceParser* parser;
- if (pDataArray->namefile != "") { parser = new SequenceParser(pDataArray->groupfile, pDataArray->fastafile, pDataArray->namefile); }
- else { parser = new SequenceParser(pDataArray->groupfile, pDataArray->fastafile); }
-
- int numSeqs = 0;
+ SequenceCountParser* cparser;
+ if (pDataArray->countfile != "") {
+ cparser = new SequenceCountParser(pDataArray->countfile, pDataArray->fastafile);
+ }else {
+ if (pDataArray->namefile != "") { parser = new SequenceParser(pDataArray->groupfile, pDataArray->fastafile, pDataArray->namefile); }
+ else { parser = new SequenceParser(pDataArray->groupfile, pDataArray->fastafile); }
+ }
+
+ int numSeqs = 0;
vector<seqPNode> alignSeqs;
//clear out old files
ofstream outF; pDataArray->m->openOutputFile(pDataArray->newFName, outF); outF.close();
pDataArray->m->mothurOutEndLine(); pDataArray->m->mothurOut("Processing group " + pDataArray->groups[k] + ":"); pDataArray->m->mothurOutEndLine();
map<string, string> thisNameMap;
- if (pDataArray->namefile != "") { thisNameMap = parser->getNameMap(pDataArray->groups[k]); }
- vector<Sequence> thisSeqs = parser->getSeqs(pDataArray->groups[k]);
+ vector<Sequence> thisSeqs;
+ if (pDataArray->groupfile != "") {
+ thisSeqs = parser->getSeqs(pDataArray->groups[k]);
+ }else if (pDataArray->countfile != "") {
+ thisSeqs = cparser->getSeqs(pDataArray->groups[k]);
+ }
+ if (pDataArray->namefile != "") { thisNameMap = parser->getNameMap(pDataArray->groups[k]); }
//fill alignSeqs with this groups info.
////////////////////////////////////////////////////
alignSeqs.clear();
map<string, string>::iterator it;
bool error = false;
+ map<string, int> thisCount;
+ if (pDataArray->countfile != "") { thisCount = cparser->getCountTable(pDataArray->groups[k]); }
+
for (int i = 0; i < thisSeqs.size(); i++) {
if (thisSeqs[i].getAligned().length() > length) { length = thisSeqs[i].getAligned().length(); }
}
}else { //no names file, you are identical to yourself
- seqPNode tempNode(1, thisSeqs[i], thisSeqs[i].getName());
- alignSeqs.push_back(tempNode);
+ int numRep = 1;
+ if (pDataArray->countfile != "") {
+ map<string, int>::iterator it2 = thisCount.find(thisSeqs[i].getName());
+
+ //should never be true since parser checks for this
+ if (it2 == thisCount.end()) { pDataArray->m->mothurOut(thisSeqs[i].getName() + " is not in your count file, please correct."); pDataArray->m->mothurOutEndLine(); error = true; }
+ else { numRep = it2->second; }
+ }
+ seqPNode tempNode(numRep, thisSeqs[i], thisSeqs[i].getName());
+ alignSeqs.push_back(tempNode);
if (thisSeqs[i].getAligned().length() > length) { length = thisSeqs[i].getAligned().length(); }
}
}
////////////////////////////////////////////////////
//int count = process(); - same function below
- //sort seqs by number of identical seqs
- sort(alignSeqs.begin(), alignSeqs.end(), comparePriority);
+ ofstream out;
+ pDataArray->m->openOutputFile(pDataArray->newMName+pDataArray->groups[k]+".map", out);
+ pDataArray->mapFileNames.push_back(pDataArray->newMName+pDataArray->groups[k]+".map");
+ //sort seqs by number of identical seqs
+ if (pDataArray->topdown) { sort(alignSeqs.begin(), alignSeqs.end(), comparePriorityTopDown); }
+ else { sort(alignSeqs.begin(), alignSeqs.end(), comparePriorityDownTop); }
+
int count = 0;
//think about running through twice...
if (alignSeqs[i].active) { //this sequence has not been merged yet
+ string chunk = alignSeqs[i].seq.getName() + "\t" + toString(alignSeqs[i].numIdentical) + "\t" + toString(0) + "\t" + alignSeqs[i].seq.getAligned() + "\n";
+
//try to merge it with all smaller seqs
for (int j = i+1; j < numSeqs; j++) {
alignSeqs[j].active = 0;
alignSeqs[j].numIdentical = 0;
+ alignSeqs[j].diffs = mismatch;
count++;
+ chunk += alignSeqs[j].seq.getName() + "\t" + toString(alignSeqs[j].numIdentical) + "\t" + toString(mismatch) + "\t" + alignSeqs[j].seq.getAligned() + "\n";
}
}//end if j active
- }//end if i != j
+ }//end for loop j
//remove from active list
alignSeqs[i].active = 0;
+ out << "ideal_seq_" << (i+1) << '\t' << alignSeqs[i].numIdentical << endl << chunk << endl;
+
}//end if active i
if(i % 100 == 0) { pDataArray->m->mothurOut(toString(i) + "\t" + toString(numSeqs - count) + "\t" + toString(count)); pDataArray->m->mothurOutEndLine(); }
}
-
+ out.close();
if(numSeqs % 100 != 0) { pDataArray->m->mothurOut(toString(numSeqs) + "\t" + toString(numSeqs - count) + "\t" + toString(count)); pDataArray->m->mothurOutEndLine(); }
////////////////////////////////////////////////////
for (int i = 0; i < alignSeqs.size(); i++) {
if (alignSeqs[i].numIdentical != 0) {
alignSeqs[i].seq.printSequence(outFasta);
- outNames << alignSeqs[i].seq.getName() << '\t' << alignSeqs[i].names << endl;
+ if (pDataArray->countfile != "") { outNames << pDataArray->groups[k] << '\t' << alignSeqs[i].seq.getName() << '\t' << alignSeqs[i].names << endl;
+ }else { outNames << alignSeqs[i].seq.getName() << '\t' << alignSeqs[i].names << endl; }
+
}
}
projects / mothur.git / blobdiff