From aca78ed4a47dff8672ea8fd93cef0dfbaf0f7495 Mon Sep 17 00:00:00 2001 From: Sarah Westcott Date: Wed, 2 Apr 2014 13:55:04 -0400 Subject: [PATCH] =?utf8?q?working=20on=20sra=20and=20get.mimarkspackage=20?= =?utf8?q?commands.=20=20added=20file=20parameter=20to=20cluster.split=20c?= =?utf8?q?ommand.=20seq.=20error=20not=20saving=20current=20fasta=20anymor?= =?utf8?q?e.=20sffinfo=20and=20fastq.info=20using=20=E2=80=9Cignore?= =?utf8?q?=E2=80=9D=20feature=20to=20play=20nicely=20with=20sra=20command.?= MIME-Version: 1.0 Content-Type: text/plain; charset=utf8 Content-Transfer-Encoding: 8bit --- clustersplitcommand.cpp | 278 +++++++++++++++++++++++------------ clustersplitcommand.h | 6 +- getmimarkspackagecommand.cpp | 12 +- getmimarkspackagecommand.h | 1 + parsefastaqcommand.cpp | 20 +++ seqerrorcommand.cpp | 7 - sffinfocommand.cpp | 94 +++++++----- sracommand.cpp | 56 +++++-- sracommand.h | 5 +- 9 files changed, 321 insertions(+), 158 deletions(-) diff --git a/clustersplitcommand.cpp b/clustersplitcommand.cpp index e61e107..546d633 100644 --- a/clustersplitcommand.cpp +++ b/clustersplitcommand.cpp @@ -13,6 +13,7 @@ //********************************************************************************************************************** vector ClusterSplitCommand::setParameters(){ try { + CommandParameter pfile("file", "InputTypes", "", "", "PhylipColumnFasta", "PhylipColumnFasta", "none","",false,false,true); parameters.push_back(pfile); CommandParameter ptaxonomy("taxonomy", "InputTypes", "", "", "none", "none", "FastaTaxName","",false,false,true); parameters.push_back(ptaxonomy); CommandParameter pphylip("phylip", "InputTypes", "", "", "PhylipColumnFasta", "PhylipColumnFasta", "none","list",false,false,true); parameters.push_back(pphylip); CommandParameter pfasta("fasta", "InputTypes", "", "", "PhylipColumnFasta", "PhylipColumnFasta", "FastaTaxName","list",false,false,true); parameters.push_back(pfasta); @@ -47,13 +48,14 @@ vector ClusterSplitCommand::setParameters(){ string ClusterSplitCommand::getHelpString(){ try { string helpString = ""; - helpString += "The cluster.split command parameter options are fasta, phylip, column, name, count, cutoff, precision, method, splitmethod, taxonomy, taxlevel, showabund, timing, hard, large, cluster, processors. Fasta or Phylip or column and name are required.\n"; + helpString += "The cluster.split command parameter options are file, fasta, phylip, column, name, count, cutoff, precision, method, splitmethod, taxonomy, taxlevel, showabund, timing, hard, large, cluster, processors. Fasta or Phylip or column and name are required.\n"; helpString += "The cluster.split command can split your files in 3 ways. Splitting by distance file, by classification, or by classification also using a fasta file. \n"; helpString += "For the distance file method, you need only provide your distance file and mothur will split the file into distinct groups. \n"; helpString += "For the classification method, you need to provide your distance file and taxonomy file, and set the splitmethod to classify. \n"; helpString += "You will also need to set the taxlevel you want to split by. mothur will split the sequences into distinct taxonomy groups, and split the distance file based on those groups. \n"; helpString += "For the classification method using a fasta file, you need to provide your fasta file, names file and taxonomy file. \n"; helpString += "You will also need to set the taxlevel you want to split by. mothur will split the sequence into distinct taxonomy groups, and create distance files for each grouping. \n"; + helpString += "The file option allows you to enter your file containing your list of column and names/count files as well as the singleton file. This file is mothur generated, when you run cluster.split() with the cluster=f parameter. This can be helpful when you have a large dataset that you may be able to use all your processors for the splitting step, but have to reduce them for the cluster step due to RAM constraints. For example: cluster.split(fasta=yourFasta, taxonomy=yourTax, count=yourCount, taxlevel=3, cluster=f, processors=8) then cluster.split(file=yourFile, processors=4). This allows your to maximize your processors during the splitting step. Also, if you are unsure if the cluster step will have RAM issue with multiple processors, you can avoid running the first part of the command multiple times.\n"; helpString += "The phylip and column parameter allow you to enter your distance file. \n"; helpString += "The fasta parameter allows you to enter your aligned fasta file. \n"; helpString += "The name parameter allows you to enter your name file. \n"; @@ -89,6 +91,7 @@ string ClusterSplitCommand::getOutputPattern(string type) { else if (type == "rabund") { pattern = "[filename],[clustertag],rabund"; } else if (type == "sabund") { pattern = "[filename],[clustertag],sabund"; } else if (type == "column") { pattern = "[filename],dist"; } + else if (type == "file") { pattern = "[filename],file"; } else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; } return pattern; @@ -108,6 +111,7 @@ ClusterSplitCommand::ClusterSplitCommand(){ outputTypes["rabund"] = tempOutNames; outputTypes["sabund"] = tempOutNames; outputTypes["column"] = tempOutNames; + outputTypes["file"] = tempOutNames; } catch(exception& e) { m->errorOut(e, "ClusterSplitCommand", "ClusterSplitCommand"); @@ -147,6 +151,7 @@ ClusterSplitCommand::ClusterSplitCommand(string option) { outputTypes["rabund"] = tempOutNames; outputTypes["sabund"] = tempOutNames; outputTypes["column"] = tempOutNames; + outputTypes["file"] = tempOutNames; //if the user changes the output directory command factory will send this info to us in the output parameter outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; } @@ -203,12 +208,25 @@ ClusterSplitCommand::ClusterSplitCommand(string option) { //if the user has not given a path then, add inputdir. else leave path alone. if (path == "") { parameters["count"] = inputDir + it->second; } } + + it = parameters.find("file"); + //user has given a template file + if(it != parameters.end()){ + path = m->hasPath(it->second); + //if the user has not given a path then, add inputdir. else leave path alone. + if (path == "") { parameters["file"] = inputDir + it->second; } + } } //check for required parameters - phylipfile = validParameter.validFile(parameters, "phylip", true); + file = validParameter.validFile(parameters, "file", true); + if (file == "not open") { file = ""; abort = true; } + else if (file == "not found") { file = ""; } + else { distfile = file; } + + phylipfile = validParameter.validFile(parameters, "phylip", true); if (phylipfile == "not open") { abort = true; } - else if (phylipfile == "not found") { phylipfile = ""; } + else if (phylipfile == "not found") { phylipfile = ""; } else { distfile = phylipfile; format = "phylip"; m->setPhylipFile(phylipfile); } columnfile = validParameter.validFile(parameters, "column", true); @@ -236,7 +254,7 @@ ClusterSplitCommand::ClusterSplitCommand(string option) { else if (taxFile == "not found") { taxFile = ""; } else { m->setTaxonomyFile(taxFile); if (splitmethod != "fasta") { splitmethod = "classify"; } } - if ((phylipfile == "") && (columnfile == "") && (fastafile == "")) { + if ((phylipfile == "") && (columnfile == "") && (fastafile == "") && (file == "")) { //is there are current file available for either of these? //give priority to column, then phylip, then fasta columnfile = m->getColumnFile(); @@ -248,13 +266,13 @@ ClusterSplitCommand::ClusterSplitCommand(string option) { fastafile = m->getFastaFile(); if (fastafile != "") { m->mothurOut("Using " + fastafile + " as input file for the fasta parameter."); m->mothurOutEndLine(); } else { - m->mothurOut("No valid current files. When executing a cluster.split command you must enter a phylip or a column or fastafile."); m->mothurOutEndLine(); + m->mothurOut("No valid current files. When executing a cluster.split command you must enter a file, phylip or a column or fastafile."); m->mothurOutEndLine(); abort = true; } } } } - else if ((phylipfile != "") && (columnfile != "") && (fastafile != "")) { m->mothurOut("When executing a cluster.split command you must enter ONLY ONE of the following: fasta, phylip or column."); m->mothurOutEndLine(); abort = true; } + else if ((phylipfile != "") && (columnfile != "") && (fastafile != "") && (file != "")) { m->mothurOut("When executing a cluster.split command you must enter ONLY ONE of the following: file, fasta, phylip or column."); m->mothurOutEndLine(); abort = true; } if ((countfile != "") && (namefile != "")) { m->mothurOut("When executing a cluster.split command you must enter ONLY ONE of the following: count or name."); m->mothurOutEndLine(); abort = true; } @@ -368,97 +386,104 @@ int ClusterSplitCommand::execute(){ try { if (abort == true) { if (calledHelp) { return 0; } return 2; } - - time_t estart; - vector listFileNames; - set labels; - string singletonName = ""; - double saveCutoff = cutoff; + + time_t estart; + vector listFileNames; + vector< map > distName; + set labels; + string singletonName = ""; + + double saveCutoff = cutoff; - //****************** file prep work ******************************// - #ifdef USE_MPI + if (file != "") { deleteFiles = false; estart = time(NULL); singletonName = readFile(distName); } + else { + + //****************** file prep work ******************************// +#ifdef USE_MPI int pid; int tag = 2001; - MPI_Status status; + MPI_Status status; MPI_Comm_size(MPI_COMM_WORLD, &processors); //set processors to the number of mpi processes running MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are if (pid == 0) { //only process 0 converts and splits - - #endif - - //if user gave a phylip file convert to column file - if (format == "phylip") { - estart = time(NULL); - m->mothurOut("Converting to column format..."); m->mothurOutEndLine(); - - ReadCluster* convert = new ReadCluster(distfile, cutoff, outputDir, false); - - NameAssignment* nameMap = NULL; - convert->setFormat("phylip"); - convert->read(nameMap); - - if (m->control_pressed) { delete convert; return 0; } - - distfile = convert->getOutputFile(); - - //if no names file given with phylip file, create it - ListVector* listToMakeNameFile = convert->getListVector(); - if ((namefile == "") && (countfile == "")) { //you need to make a namefile for split matrix - ofstream out; - namefile = phylipfile + ".names"; - m->openOutputFile(namefile, out); - for (int i = 0; i < listToMakeNameFile->getNumBins(); i++) { - string bin = listToMakeNameFile->get(i); - out << bin << '\t' << bin << endl; - } - out.close(); - } - delete listToMakeNameFile; - delete convert; - - m->mothurOut("It took " + toString(time(NULL) - estart) + " seconds to convert the distance file."); m->mothurOutEndLine(); - } - if (m->control_pressed) { return 0; } - - estart = time(NULL); - m->mothurOut("Splitting the file..."); m->mothurOutEndLine(); - - //split matrix into non-overlapping groups - SplitMatrix* split; - if (splitmethod == "distance") { split = new SplitMatrix(distfile, namefile, countfile, taxFile, cutoff, splitmethod, large); } - else if (splitmethod == "classify") { split = new SplitMatrix(distfile, namefile, countfile, taxFile, taxLevelCutoff, splitmethod, large); } - else if (splitmethod == "fasta") { split = new SplitMatrix(fastafile, namefile, countfile, taxFile, taxLevelCutoff, cutoff, splitmethod, processors, classic, outputDir); } - else { m->mothurOut("Not a valid splitting method. Valid splitting algorithms are distance, classify or fasta."); m->mothurOutEndLine(); return 0; } - - split->split(); - - if (m->control_pressed) { delete split; return 0; } - - singletonName = split->getSingletonNames(); - vector< map > distName = split->getDistanceFiles(); //returns map of distance files -> namefile sorted by distance file size - delete split; - - if (m->debug) { m->mothurOut("[DEBUG]: distName.size() = " + toString(distName.size()) + ".\n"); } - //output a merged distance file - //if (splitmethod == "fasta") { createMergedDistanceFile(distName); } +#endif + + //if user gave a phylip file convert to column file + if (format == "phylip") { + estart = time(NULL); + m->mothurOut("Converting to column format..."); m->mothurOutEndLine(); + + ReadCluster* convert = new ReadCluster(distfile, cutoff, outputDir, false); + + NameAssignment* nameMap = NULL; + convert->setFormat("phylip"); + convert->read(nameMap); + + if (m->control_pressed) { delete convert; return 0; } + + distfile = convert->getOutputFile(); + + //if no names file given with phylip file, create it + ListVector* listToMakeNameFile = convert->getListVector(); + if ((namefile == "") && (countfile == "")) { //you need to make a namefile for split matrix + ofstream out; + namefile = phylipfile + ".names"; + m->openOutputFile(namefile, out); + for (int i = 0; i < listToMakeNameFile->getNumBins(); i++) { + string bin = listToMakeNameFile->get(i); + out << bin << '\t' << bin << endl; + } + out.close(); + } + delete listToMakeNameFile; + delete convert; + + m->mothurOut("It took " + toString(time(NULL) - estart) + " seconds to convert the distance file."); m->mothurOutEndLine(); + } + if (m->control_pressed) { return 0; } + + estart = time(NULL); + m->mothurOut("Splitting the file..."); m->mothurOutEndLine(); + + //split matrix into non-overlapping groups + SplitMatrix* split; + if (splitmethod == "distance") { split = new SplitMatrix(distfile, namefile, countfile, taxFile, cutoff, splitmethod, large); } + else if (splitmethod == "classify") { split = new SplitMatrix(distfile, namefile, countfile, taxFile, taxLevelCutoff, splitmethod, large); } + else if (splitmethod == "fasta") { split = new SplitMatrix(fastafile, namefile, countfile, taxFile, taxLevelCutoff, cutoff, splitmethod, processors, classic, outputDir); } + else { m->mothurOut("Not a valid splitting method. Valid splitting algorithms are distance, classify or fasta."); m->mothurOutEndLine(); return 0; } + + split->split(); + + if (m->control_pressed) { delete split; return 0; } + + singletonName = split->getSingletonNames(); + distName = split->getDistanceFiles(); //returns map of distance files -> namefile sorted by distance file size + delete split; + + if (m->debug) { m->mothurOut("[DEBUG]: distName.size() = " + toString(distName.size()) + ".\n"); } + + //output a merged distance file + //if (splitmethod == "fasta") { createMergedDistanceFile(distName); } - if (m->control_pressed) { return 0; } - - m->mothurOut("It took " + toString(time(NULL) - estart) + " seconds to split the distance file."); m->mothurOutEndLine(); - estart = time(NULL); - - if (!runCluster) { - - m->mothurOutEndLine(); - m->mothurOut("Output File Names: "); m->mothurOutEndLine(); - for (int i = 0; i < distName.size(); i++) { m->mothurOut(distName[i].begin()->first); m->mothurOutEndLine(); m->mothurOut(distName[i].begin()->second); m->mothurOutEndLine(); } - m->mothurOutEndLine(); - return 0; + if (m->control_pressed) { return 0; } - } - + m->mothurOut("It took " + toString(time(NULL) - estart) + " seconds to split the distance file."); m->mothurOutEndLine(); + estart = time(NULL); + + if (!runCluster) { + printFile(singletonName, distName); + + m->mothurOutEndLine(); + m->mothurOut("Output File Names: "); m->mothurOutEndLine(); + for (int i = 0; i < distName.size(); i++) { m->mothurOut(distName[i].begin()->first); m->mothurOutEndLine(); m->mothurOut(distName[i].begin()->second); m->mothurOutEndLine(); } + m->mothurOutEndLine(); + return 0; + + } + deleteFiles = true; + } //****************** break up files between processes and cluster each file set ******************************// #ifdef USE_MPI ////you are process 0 from above//// @@ -651,6 +676,9 @@ int ClusterSplitCommand::execute(){ mergeLists(listFileNames, labelBins, listSingle); if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; } + + //delete after all are complete incase a crash happens + if (!deleteFiles) { for (int i = 0; i < distName.size(); i++) { m->mothurRemove(distName[i].begin()->first); m->mothurRemove(distName[i].begin()->second); } } m->mothurOut("It took " + toString(time(NULL) - estart) + " seconds to merge."); m->mothurOutEndLine(); @@ -1258,9 +1286,10 @@ string ClusterSplitCommand::clusterClassicFile(string thisDistFile, string thisN if(namefile != ""){ delete nameMap; } else { delete ct; } - m->mothurRemove(thisDistFile); - m->mothurRemove(thisNamefile); - + if (deleteFiles) { + m->mothurRemove(thisDistFile); + m->mothurRemove(thisNamefile); + } return listFileName; } @@ -1413,8 +1442,10 @@ string ClusterSplitCommand::clusterFile(string thisDistFile, string thisNamefile return listFileName; } - m->mothurRemove(thisDistFile); - m->mothurRemove(thisNamefile); + if (deleteFiles) { + m->mothurRemove(thisDistFile); + m->mothurRemove(thisNamefile); + } if (saveCutoff != cutoff) { if (hard) { saveCutoff = m->ceilDist(saveCutoff, precision); } @@ -1498,3 +1529,66 @@ int ClusterSplitCommand::createRabund(CountTable*& ct, ListVector*& list, RAbund } //********************************************************************************************************************** +int ClusterSplitCommand::printFile(string singleton, vector< map >& distName){ + try { + ofstream out; + map variables; + string thisOutputDir = outputDir; + if (outputDir == "") { thisOutputDir = m->hasPath(distfile); } + variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(distfile)); + string outputFileName = getOutputFileName("file", variables); + m->openOutputFile(outputFileName, out); + + out << singleton << endl; + if (namefile != "") { out << "name" << endl; } + else if (countfile != "") { out << "count" << endl; } + else { out << "unknown" << endl; } + + for (int i = 0; i < distName.size(); i++) { out << distName[i].begin()->first << '\t' << distName[i].begin()->second << endl; } + out.close(); + + return 0; + } + catch(exception& e) { + m->errorOut(e, "ClusterCommand", "printFile"); + exit(1); + } + +} +//********************************************************************************************************************** +string ClusterSplitCommand::readFile(vector< map >& distName){ + try { + string singleton, thiscolumn, thisname, type; + + ifstream in; + m->openInputFile(file, in); + + in >> singleton; m->gobble(in); + in >> type; m->gobble(in); + + if (type == "name") { namefile = "name"; } + else if (type == "count") { countfile = "count"; } + else { m->mothurOut("[ERROR]: unknown file type. Are the files in column 2 of the file name files or count files? Please change unknown to name or count.\n"); m->control_pressed = true; } + + + while(!in.eof()) { + if (m->control_pressed) { break; } + + in >> thiscolumn; m->gobble(in); + in >> thisname; m->gobble(in); + + map temp; + temp[thiscolumn] = thisname; + distName.push_back(temp); + } + in.close(); + + return singleton; + } + catch(exception& e) { + m->errorOut(e, "ClusterCommand", "readFile"); + exit(1); + } + +} +//********************************************************************************************************************** diff --git a/clustersplitcommand.h b/clustersplitcommand.h index e1f17b1..2e8ac63 100644 --- a/clustersplitcommand.h +++ b/clustersplitcommand.h @@ -48,10 +48,10 @@ private: vector processIDS; //processid vector outputNames; - string method, fileroot, tag, outputDir, phylipfile, columnfile, namefile, countfile, distfile, format, showabund, timing, splitmethod, taxFile, fastafile; + string file, method, fileroot, tag, outputDir, phylipfile, columnfile, namefile, countfile, distfile, format, showabund, timing, splitmethod, taxFile, fastafile; double cutoff, splitcutoff; int precision, length, processors, taxLevelCutoff; - bool print_start, abort, hard, large, classic, runCluster; + bool print_start, abort, hard, large, classic, runCluster, deleteFiles; time_t start; ofstream outList, outRabund, outSabund; @@ -64,6 +64,8 @@ private: map completeListFile(vector, string, set&, ListVector*&); int createMergedDistanceFile(vector< map >); int createRabund(CountTable*& ct, ListVector*& list, RAbundVector*& rabund); + string readFile(vector< map >&); + int printFile(string, vector< map >&); }; /////////////////not working for Windows//////////////////////////////////////////////////////////// diff --git a/getmimarkspackagecommand.cpp b/getmimarkspackagecommand.cpp index 2f37405..955cb79 100644 --- a/getmimarkspackagecommand.cpp +++ b/getmimarkspackagecommand.cpp @@ -34,7 +34,7 @@ vector GetMIMarksPackageCommand::setParameters(){ string GetMIMarksPackageCommand::getHelpString(){ try { string helpString = ""; - helpString += "The get.mimarkspackage command creates a mimarks package form with your groups. The required fields are flagged with * characters. Fields marked with '**' indicated they are in a group where at least one of the fields is required.\n"; + helpString += "The get.mimarkspackage command creates a mimarks package form with your groups. The required fields are flagged with * characters. \n"; helpString += "Further documentation on the different packages and required formats can be found here, http://www.mothur.org/wiki/MIMarks_Data_Packages.\n"; helpString += "The get.mimarkspackage command parameters are: oligos, group, package and requiredonly. oligos or group is required.\n"; helpString += "The oligos parameter is used to provide your oligos file so mothur can extract your group names.\n"; @@ -196,6 +196,8 @@ int GetMIMarksPackageCommand::execute(){ else if (file != "") { readFile(); } else { GroupMap groupmap(groupfile); groupmap.readMap(); Groups = groupmap.getNamesOfGroups(); } + for (set::iterator it = uniqueNames.begin(); it != uniqueNames.end(); it++) { Groups.push_back(*it); } + if (outputDir == "") { outputDir += m->hasPath(inputfile); } map variables; variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(inputfile)); @@ -207,7 +209,6 @@ int GetMIMarksPackageCommand::execute(){ out << "#This is a tab-delimited file. Additional Documentation can be found at http://www.mothur.org/wiki/MIMarks_Data_Packages." << endl; out << "#Please fill all the required fields indicated with '*'" << endl; - out << "#Fields marked with '**' indicated they are in a group where at least one of the fields is required." << endl; out << "#Unknown or inapplicable fields can be assigned NA value." << endl; out << "#You may add extra custom fields to this template. Make sure all the fields are separated by tabs." << endl; out << "#You may remove any fields not required (marked with '*'). Make sure all the fields are separated by tabs." << endl; @@ -223,9 +224,9 @@ int GetMIMarksPackageCommand::execute(){ }else if (package == "host_associated") { out << "#Environmental:MIMARKS.specimen.host-associated.3.0" << endl; if (requiredonly) { - out << "*sample_name *organism *collection_date *biome *feature *material *geo_loc_name *lat_lon *title *seq_methods *host **clone **isolate **strain" << endl; + out << "*sample_name *organism *collection_date *biome *feature *material *geo_loc_name *lat_lon *title *seq_methods *host " << endl; }else { - out << "*sample_name description bioproject_id sample_title *organism *collection_date *biome *feature *material *geo_loc_name *lat_lon *title *seq_methods **clone **isolate **strain rel_to_oxygen samp_collect_device samp_mat_process *host age altitude blood_press_diast blood_press_syst body_habitat body_product tissue chem_administration depth diet disease_stat dry_mass elev family_relationship genotype gravidity height_or_length host_body_temp host_color host_growth_cond host_shape host_subject_id host_taxid infra_specific_name infra_specific_rank last_meal life_stage misc_param organism_count oxy_stat_samp perturbation phenotype samp_size samp_salinity samp_store_dur samp_store_loc samp_store_temp sex substrate temp tot_mass" << endl; + out << "*sample_name description bioproject_id sample_title *organism *collection_date *biome *feature *material *geo_loc_name *lat_lon *title *seq_methods rel_to_oxygen samp_collect_device samp_mat_process *host age altitude blood_press_diast blood_press_syst body_habitat body_product tissue chem_administration depth diet disease_stat dry_mass elev family_relationship genotype gravidity height_or_length host_body_temp host_color host_growth_cond host_shape host_subject_id host_taxid infra_specific_name infra_specific_rank last_meal life_stage misc_param organism_count oxy_stat_samp perturbation phenotype samp_size samp_salinity samp_store_dur samp_store_loc samp_store_temp sex substrate temp tot_mass" << endl; } }else if (package == "human_associated") { out << "#Environmental:MIMARKS.specimen.human-associated.3.0" << endl; @@ -438,7 +439,7 @@ int GetMIMarksPackageCommand::readOligos(){ primerNameVector.push_back(""); } - set uniqueNames; + for(int i = 0; i < barcodeNameVector.size(); i++){ for(int j = 0; j < primerNameVector.size(); j++){ @@ -470,7 +471,6 @@ int GetMIMarksPackageCommand::readOligos(){ if (m->debug) { int count = 0; for (set::iterator it = uniqueNames.begin(); it != uniqueNames.end(); it++) { m->mothurOut("[DEBUG]: " + toString(count) + " groupName = " + *it + "\n"); count++; } } - for (set::iterator it = uniqueNames.begin(); it != uniqueNames.end(); it++) { Groups.push_back(*it); } return true; diff --git a/getmimarkspackagecommand.h b/getmimarkspackagecommand.h index ccef832..8bd2e49 100644 --- a/getmimarkspackagecommand.h +++ b/getmimarkspackagecommand.h @@ -36,6 +36,7 @@ private: string oligosfile, groupfile, package, inputfile, file; string outputDir; vector outputNames, Groups; + set uniqueNames; int readOligos(); int readFile(); diff --git a/parsefastaqcommand.cpp b/parsefastaqcommand.cpp index a509ae9..2f6f5bb 100644 --- a/parsefastaqcommand.cpp +++ b/parsefastaqcommand.cpp @@ -502,6 +502,26 @@ int ParseFastaQCommand::findGroup(fastqRead2 thisRead, int& barcode, int& primer if(!success) { trashCode += 'r'; } } + if (trashCode.length() == 0) { //is this sequence in the ignore group + string thisGroup = ""; + + if(barcodes.size() != 0){ + thisGroup = barcodeNameVector[barcode]; + if (numPrimers != 0) { + if (primerNameVector[primer] != "") { + if(thisGroup != "") { + thisGroup += "." + primerNameVector[primer]; + }else { + thisGroup = primerNameVector[primer]; + } + } + } + } + + int pos = thisGroup.find("ignore"); + if (pos != string::npos) { trashCode += "i"; } + } + return trashCode.length(); } diff --git a/seqerrorcommand.cpp b/seqerrorcommand.cpp index f7edb21..4e5abc9 100644 --- a/seqerrorcommand.cpp +++ b/seqerrorcommand.cpp @@ -430,13 +430,6 @@ int SeqErrorCommand::execute(){ m->mothurOut("It took " + toString(time(NULL) - start) + " secs to check " + toString(numSeqs) + " sequences."); m->mothurOutEndLine(); - //set fasta file as new current fastafile - string current = ""; - itTypes = outputTypes.find("errorseq"); - if (itTypes != outputTypes.end()) { - if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setFastaFile(current); } - } - m->mothurOutEndLine(); m->mothurOut("Output File Names: "); m->mothurOutEndLine(); for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); } diff --git a/sffinfocommand.cpp b/sffinfocommand.cpp index 03bcebb..26d4832 100755 --- a/sffinfocommand.cpp +++ b/sffinfocommand.cpp @@ -677,11 +677,13 @@ int SffInfoCommand::extractSffInfo(string input, string accnos, string oligos){ //append new header to reads for (int i = 0; i < filehandles.size(); i++) { for (int j = 0; j < filehandles[i].size(); j++) { - m->appendBinaryFiles(filehandles[i][j], filehandlesHeaders[i][j]); - m->renameFile(filehandlesHeaders[i][j], filehandles[i][j]); - m->mothurRemove(filehandlesHeaders[i][j]); - //cout << i << '\t' << '\t' << j << '\t' << filehandles[i][j] << " done appending headers and removing " << filehandlesHeaders[i][j] << endl; - if (numSplitReads[i][j] == 0) { m->mothurRemove(filehandles[i][j]); } + if (filehandles[i][j] != "") { + m->appendBinaryFiles(filehandles[i][j], filehandlesHeaders[i][j]); + m->renameFile(filehandlesHeaders[i][j], filehandles[i][j]); + m->mothurRemove(filehandlesHeaders[i][j]); + //cout << i << '\t' << '\t' << j << '\t' << filehandles[i][j] << " done appending headers and removing " << filehandlesHeaders[i][j] << endl; + if (numSplitReads[i][j] == 0) { m->mothurRemove(filehandles[i][j]); } + } } } //cout << "here3" << endl; @@ -1260,6 +1262,25 @@ int SffInfoCommand::findGroup(Header header, seqRead read, int& barcode, int& pr if(!success) { trashCode += 'r'; } } + if (trashCode.length() == 0) { //is this sequence in the ignore group + string thisGroup = ""; + + if(barcodes.size() != 0){ + thisGroup = barcodeNameVector[barcode]; + if (numFPrimers != 0) { + if (primerNameVector[primer] != "") { + if(thisGroup != "") { + thisGroup += "." + primerNameVector[primer]; + }else { + thisGroup = primerNameVector[primer]; + } + } + } + } + + int pos = thisGroup.find("ignore"); + if (pos != string::npos) { trashCode += "i"; } + } return trashCode.length(); } @@ -1911,36 +1932,39 @@ bool SffInfoCommand::readOligos(string oligoFile){ string primerName = primerNameVector[itPrimer->second]; string barcodeName = barcodeNameVector[itBar->second]; - string comboGroupName = ""; - string fastaFileName = ""; - string qualFileName = ""; - string nameFileName = ""; - - if(primerName == ""){ - comboGroupName = barcodeNameVector[itBar->second]; - } - else{ - if(barcodeName == ""){ - comboGroupName = primerNameVector[itPrimer->second]; - } - else{ - comboGroupName = barcodeNameVector[itBar->second] + "." + primerNameVector[itPrimer->second]; - } - } - - ofstream temp; - map variables; - variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(currentFileName)); - variables["[group]"] = comboGroupName; - string thisFilename = getOutputFileName("sff",variables); - if (uniqueNames.count(thisFilename) == 0) { - outputNames.push_back(thisFilename); - outputTypes["sff"].push_back(thisFilename); - uniqueNames.insert(thisFilename); - } - - filehandles[itBar->second][itPrimer->second] = thisFilename; - temp.open(thisFilename.c_str(), ios::binary); temp.close(); + if ((primerName == "ignore") || (barcodeName == "ignore")) { } //do nothing + else { + string comboGroupName = ""; + string fastaFileName = ""; + string qualFileName = ""; + string nameFileName = ""; + + if(primerName == ""){ + comboGroupName = barcodeNameVector[itBar->second]; + } + else{ + if(barcodeName == ""){ + comboGroupName = primerNameVector[itPrimer->second]; + } + else{ + comboGroupName = barcodeNameVector[itBar->second] + "." + primerNameVector[itPrimer->second]; + } + } + + ofstream temp; + map variables; + variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(currentFileName)); + variables["[group]"] = comboGroupName; + string thisFilename = getOutputFileName("sff",variables); + if (uniqueNames.count(thisFilename) == 0) { + outputNames.push_back(thisFilename); + outputTypes["sff"].push_back(thisFilename); + uniqueNames.insert(thisFilename); + } + + filehandles[itBar->second][itPrimer->second] = thisFilename; + temp.open(thisFilename.c_str(), ios::binary); temp.close(); + } } } } diff --git a/sracommand.cpp b/sracommand.cpp index 02e9899..cde7955 100644 --- a/sracommand.cpp +++ b/sracommand.cpp @@ -313,6 +313,8 @@ int SRACommand::execute(){ else if (sfffile != "") { parseSffFile(filesBySample); } else if (fastqfile != "") { parseFastqFile(filesBySample); } + for (set::iterator it = uniqueNames.begin(); it != uniqueNames.end(); it++) { Groups.push_back(*it); } + sanityCheckMiMarksGroups(); //checks groups and files returned from parse - removes any groups that did not get reads assigned to them, orders files. @@ -386,7 +388,7 @@ int SRACommand::execute(){ //////////////////////////////////////////////////////// for (int i = 0; i < Groups.size(); i++) { - string barcodeForThisSample = Group2Barcode[Groups[i]]; + string barcodeForThisSample = Group2Barcode[Groups[i]][0]; if (m->control_pressed) { break; } out << "\t\n"; @@ -430,7 +432,7 @@ int SRACommand::execute(){ for (int i = 0; i < Groups.size(); i++) { vector thisGroupsFiles = filesBySample[Groups[i]]; - string barcodeForThisSample = Group2Barcode[Groups[i]]; + string barcodeForThisSample = Group2Barcode[Groups[i]][0]; for (int j = 0; j < thisGroupsFiles.size(); j++) { string libId = thisGroupsFiles[j] + "." + barcodeForThisSample; @@ -444,10 +446,10 @@ int SRACommand::execute(){ out << "\t\t\t\n"; out << "\t\t\t\tgeneric-data \n"; out << "\t\t\t\n"; - vector thisBarcodes; m->splitAtChar(Group2Barcode[Groups[i]], thisBarcodes, '.'); + vector thisBarcodes; m->splitAtChar(Group2Barcode[Groups[i]][0], thisBarcodes, '.'); string forwardBarcode = thisBarcodes[0]; string reverseBarcode = thisBarcodes[1]; - vector thisPrimers; m->splitAtChar(Group2Primer[Groups[i]], thisPrimers, '.'); + vector thisPrimers; m->splitAtChar(Group2Primer[Groups[i]][0], thisPrimers, '.'); string forwardPrimer = thisPrimers[0]; string reversePrimer = thisPrimers[1]; //attributes @@ -484,8 +486,8 @@ int SRACommand::execute(){ out << "\t\t\t\n"; //attributes out << "\t\t\t" + mimarks[Groups[i]]["title"] + "\n"; - out << "\t\t\t" + Group2Barcode[Groups[i]] + "\n"; - out << "\t\t\t" + Group2Primer[Groups[i]] + "\n"; + out << "\t\t\t" + Group2Barcode[Groups[i]][0] + "\n"; + out << "\t\t\t" + Group2Primer[Groups[i]][0] + "\n"; out << "\t\t\t" + orientation + "\n"; out << "\t\t\t" + libId + "\n"; out << "\t\t\t" + libStrategy + "\n"; @@ -706,7 +708,7 @@ int SRACommand::readMIMarksFile(){ if (headers[i] == "organism") { if (!m->inUsersGroups(linePieces[i], acceptableOrganisms)) { //not an acceptable organism organismError = true; - m->mothurOut("[WARNING]: " + linePieces[i]+ " is not an acceptable organism, changing to metagenome. You can correct the issue and rerun the command, or NCBI will allow you to modify the organism after submission.\n"); linePieces[i] = "metagenome"; categories[headers[i]] = linePieces[i]; + m->mothurOut("[WARNING]: " + linePieces[i]+ " is not an acceptable organism, changing to acceptable 'metagenome'. NCBI will allow you to modify the organism after submission.\n"); linePieces[i] = "metagenome"; categories[headers[i]] = linePieces[i]; } Group2Organism[linePieces[0]] = linePieces[i]; } @@ -747,7 +749,7 @@ int SRACommand::readMIMarksFile(){ string organismTypes = ""; for (int i = 0; i < acceptableOrganisms.size()-1; i++) { organismTypes += acceptableOrganisms[i] + ", "; } organismTypes += acceptableOrganisms[acceptableOrganisms.size()-1]; - m->mothurOut("[WARNING]: The acceptable organism choices are: " + organismTypes + ".\n"); + m->mothurOut("\n[WARNING]: The acceptable organism choices are: " + organismTypes + ".\n\n\n"); } return 0; @@ -1244,7 +1246,6 @@ int SRACommand::readOligos(){ primerNameVector.push_back(""); } - set uniqueNames; //used to cleanup outputFileNames if (pairedOligos) { for(map::iterator itBar = pairedBarcodes.begin();itBar != pairedBarcodes.end();itBar++){ for(map::iterator itPrimer = pairedPrimers.begin();itPrimer != pairedPrimers.end(); itPrimer++){ @@ -1269,8 +1270,22 @@ int SRACommand::readOligos(){ } } uniqueNames.insert(comboGroupName); - Group2Barcode[comboGroupName] = (itBar->second).forward+"."+(itBar->second).reverse; - Group2Primer[comboGroupName] = (itPrimer->second).forward+"."+(itPrimer->second).reverse; + + map >::iterator itGroup2Barcode = Group2Barcode.find(comboGroupName); + if (itGroup2Barcode == Group2Barcode.end()) { + vector tempBarcodes; tempBarcodes.push_back((itBar->second).forward+"."+(itBar->second).reverse); + Group2Barcode[comboGroupName] = tempBarcodes; + }else { + Group2Barcode[comboGroupName].push_back((itBar->second).forward+"."+(itBar->second).reverse); + } + + itGroup2Barcode = Group2Primer.find(comboGroupName); + if (itGroup2Barcode == Group2Primer.end()) { + vector tempPrimers; tempPrimers.push_back((itPrimer->second).forward+"."+(itPrimer->second).reverse); + Group2Primer[comboGroupName] = tempPrimers; + }else { + Group2Primer[comboGroupName].push_back((itPrimer->second).forward+"."+(itPrimer->second).reverse); + } } } } @@ -1298,8 +1313,22 @@ int SRACommand::readOligos(){ } } uniqueNames.insert(comboGroupName); - Group2Barcode[comboGroupName] = itBar->first; - Group2Primer[comboGroupName] = itPrimer->first; + + map >::iterator itGroup2Barcode = Group2Barcode.find(comboGroupName); + if (itGroup2Barcode == Group2Barcode.end()) { + vector tempBarcodes; tempBarcodes.push_back(itBar->first); + Group2Barcode[comboGroupName] = tempBarcodes; + }else { + Group2Barcode[comboGroupName].push_back(itBar->first); + } + + itGroup2Barcode = Group2Primer.find(comboGroupName); + if (itGroup2Barcode == Group2Primer.end()) { + vector tempPrimers; tempPrimers.push_back(itPrimer->first); + Group2Primer[comboGroupName] = tempPrimers; + }else { + Group2Primer[comboGroupName].push_back(itPrimer->first); + } } } } @@ -1308,7 +1337,6 @@ int SRACommand::readOligos(){ if (m->debug) { int count = 0; for (set::iterator it = uniqueNames.begin(); it != uniqueNames.end(); it++) { m->mothurOut("[DEBUG]: " + toString(count) + " groupName = " + *it + "\n"); count++; } } - for (set::iterator it = uniqueNames.begin(); it != uniqueNames.end(); it++) { Groups.push_back(*it); } return true; diff --git a/sracommand.h b/sracommand.h index c53a17d..6c30bd9 100644 --- a/sracommand.h +++ b/sracommand.h @@ -43,10 +43,11 @@ private: vector outputNames, Groups; vector primerNameVector; vector barcodeNameVector; - map Group2Barcode; - map Group2Primer; + map > Group2Barcode; + map > Group2Primer; map Group2Organism; map > mimarks; //group -> valueForGroup> ex. F003D001 -> 42.282026 -83.733850> + set uniqueNames; bool checkCasesInstrumentModels(string&); bool checkCasesPlatforms(string&); -- 2.39.2