+ }
+ }
+
+ //close files
+ MPI_File_close(&inMPI);
+ MPI_File_close(&outMPI);
+ MPI_File_close(&outMPIAccnos);
+ if (trim) { MPI_File_close(&outMPIFasta); }
+ MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
+
+
+#endif
+ return numSeqs;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "ChimeraSlayerCommand", "MPIExecute");
+ exit(1);
+ }
+}
+//**********************************************************************************************************************
+int ChimeraSlayerCommand::deconvoluteResults(map<string, string>& uniqueNames, string outputFileName, string accnosFileName, string trimFileName){
+ try {
+ map<string, string>::iterator itUnique;
+ int total = 0;
+
+ if (trimera) { //add in more potential uniqueNames
+ map<string, string> newUniqueNames = uniqueNames;
+ for (map<string, string>::iterator it = uniqueNames.begin(); it != uniqueNames.end(); it++) {
+ newUniqueNames[(it->first)+"_LEFT"] = (it->first)+"_LEFT";
+ newUniqueNames[(it->first)+"_RIGHT"] = (it->first)+"_RIGHT";
+ }
+ uniqueNames = newUniqueNames;
+ newUniqueNames.clear();
+ }
+
+ //edit accnos file
+ ifstream in2;
+ m->openInputFile(accnosFileName, in2, "no error");
+
+ ofstream out2;
+ m->openOutputFile(accnosFileName+".temp", out2);
+
+ string name; name = "";
+ set<string> chimerasInFile;
+ set<string>::iterator itChimeras;
+
+ while (!in2.eof()) {
+ if (m->control_pressed) { in2.close(); out2.close(); m->mothurRemove(outputFileName); m->mothurRemove((accnosFileName+".temp")); return 0; }
+
+ in2 >> name; m->gobble(in2);
+
+ //find unique name
+ itUnique = uniqueNames.find(name);
+
+ if (itUnique == uniqueNames.end()) { m->mothurOut("[ERROR]: trouble parsing accnos results. Cannot find "+ name + "."); m->mothurOutEndLine(); m->control_pressed = true; }
+ else {
+ itChimeras = chimerasInFile.find((itUnique->second));
+
+ if (itChimeras == chimerasInFile.end()) {
+ out2 << itUnique->second << endl;
+ chimerasInFile.insert((itUnique->second));
+ total++;
+ }
+ }
+ }
+ in2.close();
+ out2.close();
+
+ m->mothurRemove(accnosFileName);
+ rename((accnosFileName+".temp").c_str(), accnosFileName.c_str());
+
+
+ //edit chimera file
+ ifstream in;
+ m->openInputFile(outputFileName, in);
+
+ ofstream out;
+ m->openOutputFile(outputFileName+".temp", out); out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
+
+ string rest, parent1, parent2, line;
+ set<string> namesInFile; //this is so if a sequence is found to be chimera in several samples we dont write it to the results file more than once
+ set<string>::iterator itNames;
+
+ //assumptions - in file each read will always look like...
+ /*
+ F11Fcsw_92754 no
+ F11Fcsw_63104 F11Fcsw_33372 F11Fcsw_37007 0.89441 80.4469 0.2 1.03727 93.2961 52.2 no 0-241 243-369
+ */
+
+ //get header line
+ if (!in.eof()) {
+ line = m->getline(in); m->gobble(in);
+ out << line << endl;
+ }
+
+ //for the chimera file, we want to make sure if any group finds a sequence to be chimeric then all groups do,
+ //so if this is a report that did not find it to be chimeric, but it appears in the accnos file,
+ //then ignore this report and continue until we find the report that found it to be chimeric
+
+ while (!in.eof()) {
+
+ if (m->control_pressed) { in.close(); out.close(); m->mothurRemove((outputFileName+".temp")); return 0; }
+
+ in >> name; m->gobble(in);
+ in >> parent1; m->gobble(in);
+
+ if (name == "Name") { //name = "Name" because we append the header line each time we add results from the groups
+ line = m->getline(in); m->gobble(in);
+ }else {
+ if (parent1 == "no") {
+ //find unique name
+ itUnique = uniqueNames.find(name);
+
+ if (itUnique == uniqueNames.end()) { m->mothurOut("[ERROR]: trouble parsing chimera results. Cannot find "+ name + "."); m->mothurOutEndLine(); m->control_pressed = true; }
+ else {
+ //is this sequence really not chimeric??
+ itChimeras = chimerasInFile.find(itUnique->second);
+
+ if (itChimeras == chimerasInFile.end()) {
+ //is this sequence not already in the file
+ itNames = namesInFile.find((itUnique->second));
+
+ if (itNames == namesInFile.end()) { out << itUnique->second << '\t' << "no" << endl; namesInFile.insert(itUnique->second); }
+ }
+ }
+ }else { //read the rest of the line
+ double DivQLAQRB,PerIDQLAQRB,BootStrapA,DivQLBQRA,PerIDQLBQRA,BootStrapB;
+ string flag, range1, range2;
+ bool print = false;
+ in >> parent2 >> DivQLAQRB >> PerIDQLAQRB >> BootStrapA >> DivQLBQRA >> PerIDQLBQRA >> BootStrapB >> flag >> range1 >> range2; m->gobble(in);
+
+ //find unique name
+ itUnique = uniqueNames.find(name);
+
+ if (itUnique == uniqueNames.end()) { m->mothurOut("[ERROR]: trouble parsing chimera results. Cannot find "+ name + "."); m->mothurOutEndLine(); m->control_pressed = true; }
+ else {
+ name = itUnique->second;
+ //is this name already in the file
+ itNames = namesInFile.find((name));
+
+ if (itNames == namesInFile.end()) { //no not in file
+ if (flag == "no") { //are you really a no??
+ //is this sequence really not chimeric??
+ itChimeras = chimerasInFile.find(name);
+
+ //then you really are a no so print, otherwise skip
+ if (itChimeras == chimerasInFile.end()) { print = true; }
+
+ }else{ print = true; }
+ }
+ }
+
+ if (print) {
+ out << name << '\t';
+
+ namesInFile.insert(name);
+
+ //output parent1's name
+ itUnique = uniqueNames.find(parent1);
+ if (itUnique == uniqueNames.end()) { m->mothurOut("[ERROR]: trouble parsing chimera results. Cannot find parentA "+ parent1 + "."); m->mothurOutEndLine(); m->control_pressed = true; }
+ else { out << itUnique->second << '\t'; }
+
+ //output parent2's name
+ itUnique = uniqueNames.find(parent2);
+ if (itUnique == uniqueNames.end()) { m->mothurOut("[ERROR]: trouble parsing chimera results. Cannot find parentA "+ parent2 + "."); m->mothurOutEndLine(); m->control_pressed = true; }
+ else { out << itUnique->second << '\t'; }
+
+ out << DivQLAQRB << '\t' << PerIDQLAQRB << '\t' << BootStrapA << '\t' << DivQLBQRA << '\t' << PerIDQLBQRA << '\t' << BootStrapB << '\t' << flag << '\t' << range1 << '\t' << range2 << endl;
+ }
+ }
+ }
+ }
+ in.close();
+ out.close();
+
+ m->mothurRemove(outputFileName);
+ rename((outputFileName+".temp").c_str(), outputFileName.c_str());
+
+ //edit fasta file
+ if (trim) {
+ ifstream in3;
+ m->openInputFile(trimFileName, in3);
+
+ ofstream out3;
+ m->openOutputFile(trimFileName+".temp", out3);
+
+ namesInFile.clear();
+
+ while (!in3.eof()) {
+ if (m->control_pressed) { in3.close(); out3.close(); m->mothurRemove(outputFileName); m->mothurRemove(accnosFileName); m->mothurRemove((trimFileName+".temp")); return 0; }
+
+ Sequence seq(in3); m->gobble(in3);
+
+ if (seq.getName() != "") {
+ //find unique name
+ itUnique = uniqueNames.find(seq.getName());
+
+ if (itUnique == uniqueNames.end()) { m->mothurOut("[ERROR]: trouble parsing accnos results. Cannot find "+ seq.getName() + "."); m->mothurOutEndLine(); m->control_pressed = true; }
+ else {
+ itNames = namesInFile.find((itUnique->second));
+
+ if (itNames == namesInFile.end()) {
+ seq.printSequence(out3);
+ }
+ }
+ }
+ }
+ in3.close();
+ out3.close();
+
+ m->mothurRemove(trimFileName);
+ rename((trimFileName+".temp").c_str(), trimFileName.c_str());
+ }
+
+ return total;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "ChimeraSlayerCommand", "deconvoluteResults");
+ exit(1);
+ }
+}
+//**********************************************************************************************************************
+int ChimeraSlayerCommand::setUpForSelfReference(SequenceParser*& parser, map<string, string>& fileGroup, map<string, map<string, int> >& fileToPriority, int s){
+ try {
+ fileGroup.clear();
+ fileToPriority.clear();
+
+ string nameFile = "";
+ if (nameFileNames.size() != 0) { //you provided a namefile and we don't need to create one
+ nameFile = nameFileNames[s];
+ }else { nameFile = getNamesFile(fastaFileNames[s]); }
+
+ //you provided a groupfile
+ string groupFile = "";
+ if (groupFileNames.size() != 0) { groupFile = groupFileNames[s]; }
+
+ if (groupFile == "") {
+ if (processors != 1) { m->mothurOut("When using template=self, mothur can only use 1 processor, continuing."); m->mothurOutEndLine(); processors = 1; }
+
+ //sort fastafile by abundance, returns new sorted fastafile name
+ m->mothurOut("Sorting fastafile according to abundance..."); cout.flush();
+ priority = sortFastaFile(fastaFileNames[s], nameFile);
+ m->mothurOut("Done."); m->mothurOutEndLine();
+
+ fileToPriority[fastaFileNames[s]] = priority;
+ fileGroup[fastaFileNames[s]] = "noGroup";
+ }else {
+ //Parse sequences by group
+ parser = new SequenceParser(groupFile, fastaFileNames[s], nameFile);
+ vector<string> groups = parser->getNamesOfGroups();
+
+ for (int i = 0; i < groups.size(); i++) {
+ vector<Sequence> thisGroupsSeqs = parser->getSeqs(groups[i]);
+ map<string, string> thisGroupsMap = parser->getNameMap(groups[i]);
+ group2NameMap[groups[i]] = thisGroupsMap;
+ string newFastaFile = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[s])) + groups[i] + "-sortedTemp.fasta";
+ priority = sortFastaFile(thisGroupsSeqs, thisGroupsMap, newFastaFile);
+ fileToPriority[newFastaFile] = priority;
+ fileGroup[newFastaFile] = groups[i];
+ }
+ }
+
+
+ return 0;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "ChimeraSlayerCommand", "setUpForSelfReference");
+ exit(1);
+ }
+}
+//**********************************************************************************************************************
+int ChimeraSlayerCommand::setUpForSelfReference(SequenceCountParser*& parser, map<string, string>& fileGroup, map<string, map<string, int> >& fileToPriority, int s){
+ try {
+ fileGroup.clear();
+ fileToPriority.clear();
+
+ string nameFile = "";
+ if (nameFileNames.size() != 0) { //you provided a namefile and we don't need to create one
+ nameFile = nameFileNames[s];
+ }else { m->control_pressed = true; return 0; }
+
+ CountTable ct;
+ if (!ct.testGroups(nameFile)) {
+ if (processors != 1) { m->mothurOut("When using template=self, mothur can only use 1 processor, continuing."); m->mothurOutEndLine(); processors = 1; }
+
+ //sort fastafile by abundance, returns new sorted fastafile name
+ m->mothurOut("Sorting fastafile according to abundance..."); cout.flush();
+ priority = sortFastaFile(fastaFileNames[s], nameFile);
+ m->mothurOut("Done."); m->mothurOutEndLine();
+
+ fileToPriority[fastaFileNames[s]] = priority;
+ fileGroup[fastaFileNames[s]] = "noGroup";
+ }else {
+ //Parse sequences by group
+ parser = new SequenceCountParser(nameFile, fastaFileNames[s]);
+ vector<string> groups = parser->getNamesOfGroups();
+
+ for (int i = 0; i < groups.size(); i++) {
+ vector<Sequence> thisGroupsSeqs = parser->getSeqs(groups[i]);
+ map<string, int> thisGroupsMap = parser->getCountTable(groups[i]);
+ string newFastaFile = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[s])) + groups[i] + "-sortedTemp.fasta";
+ sortFastaFile(thisGroupsSeqs, thisGroupsMap, newFastaFile);
+ fileToPriority[newFastaFile] = thisGroupsMap;
+ fileGroup[newFastaFile] = groups[i];
+ }
+ }
+
+
+ return 0;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "ChimeraSlayerCommand", "setUpForSelfReference");
+ exit(1);
+ }
+}
+//**********************************************************************************************************************
+string ChimeraSlayerCommand::getNamesFile(string& inputFile){
+ try {
+ string nameFile = "";
+
+ m->mothurOutEndLine(); m->mothurOut("No namesfile given, running unique.seqs command to generate one."); m->mothurOutEndLine(); m->mothurOutEndLine();
+
+ //use unique.seqs to create new name and fastafile
+ string inputString = "fasta=" + inputFile;
+ m->mothurOut("/******************************************/"); m->mothurOutEndLine();
+ m->mothurOut("Running command: unique.seqs(" + inputString + ")"); m->mothurOutEndLine();
+ m->mothurCalling = true;
+
+ Command* uniqueCommand = new DeconvoluteCommand(inputString);
+ uniqueCommand->execute();
+
+ map<string, vector<string> > filenames = uniqueCommand->getOutputFiles();
+
+ delete uniqueCommand;
+ m->mothurCalling = false;
+ m->mothurOut("/******************************************/"); m->mothurOutEndLine();
+
+ nameFile = filenames["name"][0];
+ inputFile = filenames["fasta"][0];
+
+ return nameFile;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "ChimeraSlayerCommand", "getNamesFile");
+ exit(1);
+ }
+}
+//**********************************************************************************************************************
+
+int ChimeraSlayerCommand::driverGroups(string outputFName, string accnos, string fasta, map<string, map<string, int> >& fileToPriority, map<string, string>& fileGroup, string countlist){
+ try {
+ int totalSeqs = 0;
+ ofstream outCountList;
+
+ if (hasCount && dups) { m->openOutputFile(countlist, outCountList); }
+
+ for (map<string, map<string, int> >::iterator itFile = fileToPriority.begin(); itFile != fileToPriority.end(); itFile++) {
+
+ if (m->control_pressed) { return 0; }
+
+ int start = time(NULL);
+ string thisFastaName = itFile->first;
+ map<string, int> thisPriority = itFile->second;
+ string thisoutputFileName = outputDir + m->getRootName(m->getSimpleName(thisFastaName)) + fileGroup[thisFastaName] + "slayer.chimera";
+ string thisaccnosFileName = outputDir + m->getRootName(m->getSimpleName(thisFastaName)) + fileGroup[thisFastaName] + "slayer.accnos";
+ string thistrimFastaFileName = outputDir + m->getRootName(m->getSimpleName(thisFastaName)) + fileGroup[thisFastaName] + "slayer.fasta";
+
+ m->mothurOutEndLine(); m->mothurOut("Checking sequences from group: " + fileGroup[thisFastaName] + "."); m->mothurOutEndLine();
+
+ lines.clear();
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+ int proc = 1;
+ vector<unsigned long long> positions = m->divideFile(thisFastaName, proc);
+ lines.push_back(linePair(positions[0], positions[1]));
+#else
+ lines.push_back(linePair(0, 1000));
+#endif
+ int numSeqs = driver(lines[0], thisoutputFileName, thisFastaName, thisaccnosFileName, thistrimFastaFileName, thisPriority);
+
+ //if we provided a count file with group info and set dereplicate=t, then we want to create a *.pick.count_table
+ //This table will zero out group counts for seqs determined to be chimeric by that group.
+ if (dups) {
+ if (!m->isBlank(thisaccnosFileName)) {
+ ifstream in;
+ m->openInputFile(thisaccnosFileName, in);
+ string name;
+ if (hasCount) {
+ while (!in.eof()) {
+ in >> name; m->gobble(in);
+ outCountList << name << '\t' << fileGroup[thisFastaName] << endl;
+ }
+ in.close();
+ }else {
+ map<string, map<string, string> >::iterator itGroupNameMap = group2NameMap.find(fileGroup[thisFastaName]);
+ if (itGroupNameMap != group2NameMap.end()) {
+ map<string, string> thisnamemap = itGroupNameMap->second;
+ map<string, string>::iterator itN;
+ ofstream out;
+ m->openOutputFile(thisaccnosFileName+".temp", out);
+ while (!in.eof()) {
+ in >> name; m->gobble(in);
+ itN = thisnamemap.find(name);
+ if (itN != thisnamemap.end()) {
+ vector<string> tempNames; m->splitAtComma(itN->second, tempNames);
+ for (int j = 0; j < tempNames.size(); j++) { out << tempNames[j] << endl; }
+
+ }else { m->mothurOut("[ERROR]: parsing cannot find " + name + ".\n"); m->control_pressed = true; }
+ }
+ out.close();
+ in.close();
+ m->renameFile(thisaccnosFileName+".temp", thisaccnosFileName);
+ }else { m->mothurOut("[ERROR]: parsing cannot find " + fileGroup[thisFastaName] + ".\n"); m->control_pressed = true; }
+ }
+
+ }
+ }