2 // makecontigscommand.cpp
5 // Created by Sarah Westcott on 5/15/12.
6 // Copyright (c) 2012 Schloss Lab. All rights reserved.
9 #include "makecontigscommand.h"
11 //**********************************************************************************************************************
12 vector<string> MakeContigsCommand::setParameters(){
14 CommandParameter pfastq("ffastq", "InputTypes", "", "", "FastaFastqFile", "FastaFastqFile", "fastqGroup","fasta-qfile",false,false,true); parameters.push_back(pfastq);
15 CommandParameter prfastq("rfastq", "InputTypes", "", "", "none", "none", "fastqGroup","fasta-qfile",false,false,true); parameters.push_back(prfastq);
16 CommandParameter pfasta("ffasta", "InputTypes", "", "", "FastaFastqFile", "FastaFastqFile", "fastaGroup","fasta",false,false,true); parameters.push_back(pfasta);
17 CommandParameter prfasta("rfasta", "InputTypes", "", "", "none", "none", "none","fastaGroup",false,false,true); parameters.push_back(prfasta);
18 CommandParameter pfqual("fqfile", "InputTypes", "", "", "none", "none", "qfileGroup","",false,false,true); parameters.push_back(pfqual);
19 CommandParameter prqual("rqfile", "InputTypes", "", "", "none", "none", "qfileGroup","",false,false,true); parameters.push_back(prqual);
20 CommandParameter pfile("file", "InputTypes", "", "", "FastaFastqFile", "FastaFastqFile", "none","fasta-qfile",false,false,true); parameters.push_back(pfile);
21 CommandParameter poligos("oligos", "InputTypes", "", "", "none", "none", "none","group",false,false,true); parameters.push_back(poligos);
22 CommandParameter ppdiffs("pdiffs", "Number", "", "0", "", "", "","",false,false,true); parameters.push_back(ppdiffs);
23 CommandParameter pbdiffs("bdiffs", "Number", "", "0", "", "", "","",false,false,true); parameters.push_back(pbdiffs);
24 // CommandParameter pldiffs("ldiffs", "Number", "", "0", "", "", "",false,false); parameters.push_back(pldiffs);
25 // CommandParameter psdiffs("sdiffs", "Number", "", "0", "", "", "",false,false); parameters.push_back(psdiffs);
26 CommandParameter ptdiffs("tdiffs", "Number", "", "0", "", "", "","",false,false); parameters.push_back(ptdiffs);
28 CommandParameter palign("align", "Multiple", "needleman-gotoh", "needleman", "", "", "","",false,false); parameters.push_back(palign);
29 CommandParameter pallfiles("allfiles", "Boolean", "", "F", "", "", "","",false,false); parameters.push_back(pallfiles);
30 CommandParameter pmatch("match", "Number", "", "1.0", "", "", "","",false,false); parameters.push_back(pmatch);
31 CommandParameter pmismatch("mismatch", "Number", "", "-1.0", "", "", "","",false,false); parameters.push_back(pmismatch);
32 CommandParameter pgapopen("gapopen", "Number", "", "-2.0", "", "", "","",false,false); parameters.push_back(pgapopen);
33 CommandParameter pgapextend("gapextend", "Number", "", "-1.0", "", "", "","",false,false); parameters.push_back(pgapextend);
34 CommandParameter pthreshold("insert", "Number", "", "25", "", "", "","",false,false); parameters.push_back(pthreshold);
35 CommandParameter pdeltaq("deltaq", "Number", "", "6", "", "", "","",false,false); parameters.push_back(pdeltaq);
36 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors);
37 CommandParameter pformat("format", "Multiple", "sanger-illumina-solexa-illumina1.8+", "illumina1.8+", "", "", "","",false,false,true); parameters.push_back(pformat);
38 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
39 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
41 vector<string> myArray;
42 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
46 m->errorOut(e, "MakeContigsCommand", "setParameters");
50 //**********************************************************************************************************************
51 string MakeContigsCommand::getHelpString(){
53 string helpString = "";
54 helpString += "The make.contigs command reads a file, forward fastq file and a reverse fastq file or forward fasta and reverse fasta files and outputs new fasta. It will also provide new quality files if the fastq or file parameter is used.\n";
55 helpString += "If an oligos file is provided barcodes and primers will be trimmed, and a group file will be created.\n";
56 helpString += "The make.contigs command parameters are file, ffastq, rfastq, ffasta, rfasta, fqfile, rqfile, oligos, format, tdiffs, bdiffs, pdiffs, align, match, mismatch, gapopen, gapextend, insert, deltaq, allfiles and processors.\n";
57 helpString += "The ffastq and rfastq, file, or ffasta and rfasta parameters are required.\n";
58 helpString += "The file parameter is 2 column file containing the forward fastq files in the first column and their matching reverse fastq files in the second column. Mothur will process each pair and create a combined fasta and report file with all the sequences.\n";
59 helpString += "The ffastq and rfastq parameters are used to provide a forward fastq and reverse fastq file to process. If you provide one, you must provide the other.\n";
60 helpString += "The ffasta and rfasta parameters are used to provide a forward fasta and reverse fasta file to process. If you provide one, you must provide the other.\n";
61 helpString += "The fqfile and rqfile parameters are used to provide a forward quality and reverse quality files to process with the ffasta and rfasta parameters. If you provide one, you must provide the other.\n";
62 helpString += "The format parameter is used to indicate whether your sequences are sanger, solexa, illumina1.8+ or illumina, default=illumina1.8+.\n";
63 helpString += "The align parameter allows you to specify the alignment method to use. Your options are: gotoh and needleman. The default is needleman.\n";
64 helpString += "The tdiffs parameter is used to specify the total number of differences allowed in the sequence. The default is pdiffs + bdiffs + sdiffs + ldiffs.\n";
65 helpString += "The bdiffs parameter is used to specify the number of differences allowed in the barcode. The default is 0.\n";
66 helpString += "The pdiffs parameter is used to specify the number of differences allowed in the primer. The default is 0.\n";
67 //helpString += "The ldiffs parameter is used to specify the number of differences allowed in the linker. The default is 0.\n";
68 //helpString += "The sdiffs parameter is used to specify the number of differences allowed in the spacer. The default is 0.\n";
69 helpString += "The match parameter allows you to specify the bonus for having the same base. The default is 1.0.\n";
70 helpString += "The mistmatch parameter allows you to specify the penalty for having different bases. The default is -1.0.\n";
71 helpString += "The deltaq parameter allows you to specify the delta allowed between quality scores of a mismatched base. For example in the overlap, if deltaq=5 and in the alignment seqA, pos 200 has a quality score of 30 and the same position in seqB has a quality score of 20, you take the base from seqA (30-20 >= 5). If the quality score in seqB is 28 then the base in the consensus will be an N (30-28<5) The default is 6.\n";
72 helpString += "The gapopen parameter allows you to specify the penalty for opening a gap in an alignment. The default is -2.0.\n";
73 helpString += "The gapextend parameter allows you to specify the penalty for extending a gap in an alignment. The default is -1.0.\n";
74 helpString += "The insert parameter allows you to set a quality scores threshold. In the case where we are trying to decide whether to keep a base or remove it because the base is compared to a gap in the other fragment, if the base has a quality score below the threshold we eliminate it. Default=25.\n";
75 helpString += "The processors parameter allows you to specify how many processors you would like to use. The default is 1. \n";
76 helpString += "The allfiles parameter will create separate group and fasta file for each grouping. The default is F.\n";
77 helpString += "The make.contigs command should be in the following format: \n";
78 helpString += "make.contigs(ffastq=yourForwardFastqFile, rfastq=yourReverseFastqFile, align=yourAlignmentMethod) \n";
79 helpString += "Note: No spaces between parameter labels (i.e. ffastq), '=' and parameters (i.e.yourForwardFastqFile).\n";
83 m->errorOut(e, "MakeContigsCommand", "getHelpString");
87 //**********************************************************************************************************************
88 string MakeContigsCommand::getOutputPattern(string type) {
92 if (type == "fasta") { pattern = "[filename],[tag],contigs.fasta"; }
93 else if (type == "group") { pattern = "[filename],[tag],contigs.groups"; }
94 else if (type == "report") { pattern = "[filename],[tag],contigs.report"; }
95 else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; }
100 m->errorOut(e, "MakeContigsCommand", "getOutputPattern");
104 //**********************************************************************************************************************
105 MakeContigsCommand::MakeContigsCommand(){
107 abort = true; calledHelp = true;
109 vector<string> tempOutNames;
110 outputTypes["fasta"] = tempOutNames;
111 outputTypes["group"] = tempOutNames;
112 outputTypes["report"] = tempOutNames;
114 catch(exception& e) {
115 m->errorOut(e, "MakeContigsCommand", "MakeContigsCommand");
119 //**********************************************************************************************************************
120 MakeContigsCommand::MakeContigsCommand(string option) {
122 abort = false; calledHelp = false;
124 //allow user to run help
125 if(option == "help") { help(); abort = true; calledHelp = true; }
126 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
129 vector<string> myArray = setParameters();
131 OptionParser parser(option);
132 map<string, string> parameters = parser.getParameters();
134 ValidParameters validParameter("pairwise.seqs");
135 map<string, string>::iterator it;
137 //check to make sure all parameters are valid for command
138 for (it = parameters.begin(); it != parameters.end(); it++) {
139 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
142 //initialize outputTypes
143 vector<string> tempOutNames;
144 outputTypes["fasta"] = tempOutNames;
145 outputTypes["report"] = tempOutNames;
146 outputTypes["group"] = tempOutNames;
149 //if the user changes the input directory command factory will send this info to us in the output parameter
150 string inputDir = validParameter.validFile(parameters, "inputdir", false);
151 if (inputDir == "not found"){ inputDir = ""; }
154 it = parameters.find("ffastq");
155 //user has given a template file
156 if(it != parameters.end()){
157 path = m->hasPath(it->second);
158 //if the user has not given a path then, add inputdir. else leave path alone.
159 if (path == "") { parameters["ffastq"] = inputDir + it->second; }
162 it = parameters.find("rfastq");
163 //user has given a template file
164 if(it != parameters.end()){
165 path = m->hasPath(it->second);
166 //if the user has not given a path then, add inputdir. else leave path alone.
167 if (path == "") { parameters["rfastq"] = inputDir + it->second; }
170 it = parameters.find("ffasta");
171 //user has given a template file
172 if(it != parameters.end()){
173 path = m->hasPath(it->second);
174 //if the user has not given a path then, add inputdir. else leave path alone.
175 if (path == "") { parameters["ffasta"] = inputDir + it->second; }
178 it = parameters.find("rfasta");
179 //user has given a template file
180 if(it != parameters.end()){
181 path = m->hasPath(it->second);
182 //if the user has not given a path then, add inputdir. else leave path alone.
183 if (path == "") { parameters["rfasta"] = inputDir + it->second; }
186 it = parameters.find("fqfile");
187 //user has given a template file
188 if(it != parameters.end()){
189 path = m->hasPath(it->second);
190 //if the user has not given a path then, add inputdir. else leave path alone.
191 if (path == "") { parameters["fqfile"] = inputDir + it->second; }
194 it = parameters.find("rqfile");
195 //user has given a template file
196 if(it != parameters.end()){
197 path = m->hasPath(it->second);
198 //if the user has not given a path then, add inputdir. else leave path alone.
199 if (path == "") { parameters["rqfile"] = inputDir + it->second; }
202 it = parameters.find("file");
203 //user has given a template file
204 if(it != parameters.end()){
205 path = m->hasPath(it->second);
206 //if the user has not given a path then, add inputdir. else leave path alone.
207 if (path == "") { parameters["file"] = inputDir + it->second; }
210 it = parameters.find("oligos");
211 //user has given a template file
212 if(it != parameters.end()){
213 path = m->hasPath(it->second);
214 //if the user has not given a path then, add inputdir. else leave path alone.
215 if (path == "") { parameters["oligos"] = inputDir + it->second; }
219 ffastqfile = validParameter.validFile(parameters, "ffastq", true);
220 if (ffastqfile == "not open") { abort = true; }
221 else if (ffastqfile == "not found") { ffastqfile = ""; }
223 rfastqfile = validParameter.validFile(parameters, "rfastq", true);
224 if (rfastqfile == "not open") { abort = true; }
225 else if (rfastqfile == "not found") { rfastqfile = ""; }
227 ffastafile = validParameter.validFile(parameters, "ffasta", true);
228 if (ffastafile == "not open") { abort = true; }
229 else if (ffastafile == "not found") { ffastafile = ""; }
231 rfastafile = validParameter.validFile(parameters, "rfasta", true);
232 if (rfastafile == "not open") { abort = true; }
233 else if (rfastafile == "not found") { rfastafile = ""; }
235 fqualfile = validParameter.validFile(parameters, "fqfile", true);
236 if (fqualfile == "not open") { abort = true; }
237 else if (fqualfile == "not found") { fqualfile = ""; }
239 rqualfile = validParameter.validFile(parameters, "rqfile", true);
240 if (rqualfile == "not open") { abort = true; }
241 else if (rqualfile == "not found") { rqualfile = ""; }
243 file = validParameter.validFile(parameters, "file", true);
244 if (file == "not open") { abort = true; }
245 else if (file == "not found") { file = ""; }
248 if ((file == "") && (ffastafile == "") && (ffastqfile == "")) { abort = true; m->mothurOut("[ERROR]: The file, ffastq and rfastq or ffasta and rfasta parameters are required.\n"); }
249 if ((file != "") && ((ffastafile != "") || (ffastqfile != ""))) { abort = true; m->mothurOut("[ERROR]: The file, ffastq and rfastq or ffasta and rfasta parameters are required.\n"); }
250 if ((ffastqfile != "") && (rfastqfile == "")) { abort = true; m->mothurOut("[ERROR]: If you provide use the ffastq, you must provide a rfastq file.\n"); }
251 if ((ffastqfile == "") && (rfastqfile != "")) { abort = true; m->mothurOut("[ERROR]: If you provide use the rfastq, you must provide a ffastq file.\n"); }
252 if ((ffastafile != "") && (rfastafile == "")) { abort = true; m->mothurOut("[ERROR]: If you provide use the ffasta, you must provide a rfasta file.\n"); }
253 if ((ffastafile == "") && (rfastafile != "")) { abort = true; m->mothurOut("[ERROR]: If you provide use the rfasta, you must provide a ffasta file.\n"); }
254 if ((fqualfile != "") && (rqualfile == "")) { abort = true; m->mothurOut("[ERROR]: If you provide use the fqfile, you must provide a rqfile file.\n"); }
255 if ((fqualfile == "") && (rqualfile != "")) { abort = true; m->mothurOut("[ERROR]: If you provide use the rqfile, you must provide a fqfile file.\n"); }
256 if (((fqualfile != "") || (rqualfile != "")) && ((ffastafile == "") || (rfastafile == ""))) {
257 abort = true; m->mothurOut("[ERROR]: If you provide use the rqfile or fqfile file, you must provide the ffasta and rfasta parameters.\n");
260 oligosfile = validParameter.validFile(parameters, "oligos", true);
261 if (oligosfile == "not found") { oligosfile = ""; }
262 else if(oligosfile == "not open") { abort = true; }
263 else { m->setOligosFile(oligosfile); }
265 //if the user changes the output directory command factory will send this info to us in the output parameter
266 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){
271 //check for optional parameter and set defaults
272 // ...at some point should added some additional type checking...
274 temp = validParameter.validFile(parameters, "match", false); if (temp == "not found"){ temp = "1.0"; }
275 m->mothurConvert(temp, match);
277 temp = validParameter.validFile(parameters, "mismatch", false); if (temp == "not found"){ temp = "-1.0"; }
278 m->mothurConvert(temp, misMatch);
279 if (misMatch > 0) { m->mothurOut("[ERROR]: mismatch must be negative.\n"); abort=true; }
281 temp = validParameter.validFile(parameters, "gapopen", false); if (temp == "not found"){ temp = "-2.0"; }
282 m->mothurConvert(temp, gapOpen);
283 if (gapOpen > 0) { m->mothurOut("[ERROR]: gapopen must be negative.\n"); abort=true; }
285 temp = validParameter.validFile(parameters, "gapextend", false); if (temp == "not found"){ temp = "-1.0"; }
286 m->mothurConvert(temp, gapExtend);
287 if (gapExtend > 0) { m->mothurOut("[ERROR]: gapextend must be negative.\n"); abort=true; }
289 temp = validParameter.validFile(parameters, "insert", false); if (temp == "not found"){ temp = "25"; }
290 m->mothurConvert(temp, insert);
291 if ((insert < 0) || (insert > 40)) { m->mothurOut("[ERROR]: insert must be between 0 and 40.\n"); abort=true; }
293 temp = validParameter.validFile(parameters, "deltaq", false); if (temp == "not found"){ temp = "6"; }
294 m->mothurConvert(temp, deltaq);
296 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
297 m->setProcessors(temp);
298 m->mothurConvert(temp, processors);
300 temp = validParameter.validFile(parameters, "bdiffs", false); if (temp == "not found") { temp = "0"; }
301 m->mothurConvert(temp, bdiffs);
303 temp = validParameter.validFile(parameters, "pdiffs", false); if (temp == "not found") { temp = "0"; }
304 m->mothurConvert(temp, pdiffs);
306 // temp = validParameter.validFile(parameters, "ldiffs", false); if (temp == "not found") { temp = "0"; }
307 // m->mothurConvert(temp, ldiffs);
310 // temp = validParameter.validFile(parameters, "sdiffs", false); if (temp == "not found") { temp = "0"; }
311 // m->mothurConvert(temp, sdiffs);
314 temp = validParameter.validFile(parameters, "tdiffs", false); if (temp == "not found") { int tempTotal = pdiffs + bdiffs; temp = toString(tempTotal); }
315 m->mothurConvert(temp, tdiffs);
317 if(tdiffs == 0){ tdiffs = bdiffs + pdiffs; } //+ ldiffs + sdiffs;
319 temp = validParameter.validFile(parameters, "allfiles", false); if (temp == "not found") { temp = "F"; }
320 allFiles = m->isTrue(temp);
322 align = validParameter.validFile(parameters, "align", false); if (align == "not found"){ align = "needleman"; }
323 if ((align != "needleman") && (align != "gotoh")) { m->mothurOut(align + " is not a valid alignment method. Options are needleman or gotoh. I will use needleman."); m->mothurOutEndLine(); align = "needleman"; }
325 format = validParameter.validFile(parameters, "format", false); if (format == "not found"){ format = "illumina1.8+"; }
327 if ((format != "sanger") && (format != "illumina") && (format != "illumina1.8+") && (format != "solexa")) {
328 m->mothurOut(format + " is not a valid format. Your format choices are sanger, solexa, illumina1.8+ and illumina, aborting." ); m->mothurOutEndLine();
332 //fill convert table - goes from solexa to sanger. Used fq_all2std.pl as a reference.
333 for (int i = -64; i < 65; i++) {
334 char temp = (char) ((int)(33 + 10*log(1+pow(10,(i/10.0)))/log(10)+0.499));
335 convertTable.push_back(temp);
340 catch(exception& e) {
341 m->errorOut(e, "MakeContigsCommand", "MakeContigsCommand");
345 //**********************************************************************************************************************
346 int MakeContigsCommand::execute(){
348 if (abort == true) { if (calledHelp) { return 0; } return 2; }
350 //read ffastq and rfastq files creating fasta and qual files.
351 //this function will create a forward and reverse, fasta and qual files for each processor.
352 //files has an entry for each processor. files[i][0] = forwardFasta, files[i][1] = forwardQual, files[i][2] = reverseFasta, files[i][3] = reverseQual. filesToProcess is for each filepair in the file parameter file. for ffastq and rfastq this will be size 1.
353 unsigned long int numReads = 0;
354 int start = time(NULL);
356 m->mothurOut("Reading fastq data...\n");
357 vector < vector< vector<string> > > filesToProcess = preProcessData(numReads);
358 m->mothurOut("Done.\n");
360 if (m->control_pressed) { return 0; }
362 map<string, string> cvars;
363 string compOutputDir = outputDir;
364 if (outputDir == "") { compOutputDir = m->hasPath(file); }
365 cvars["[filename]"] = compOutputDir + m->getRootName(m->getSimpleName(file));
367 string compositeGroupFile = getOutputFileName("group",cvars);
368 cvars["[tag]"] = "trim";
369 string compositeFastaFile = getOutputFileName("fasta",cvars);
370 cvars["[tag]"] = "scrap";
371 string compositeScrapFastaFile = getOutputFileName("fasta",cvars);
373 string compositeMisMatchFile = getOutputFileName("report",cvars);
375 if (filesToProcess.size() > 1) { //clear files for append below
376 ofstream outCTFasta, outCTQual, outCSFasta, outCSQual, outCMisMatch;
377 m->openOutputFile(compositeFastaFile, outCTFasta); outCTFasta.close();
378 m->openOutputFile(compositeScrapFastaFile, outCSFasta); outCSFasta.close();
379 m->openOutputFile(compositeMisMatchFile, outCMisMatch); outCMisMatch.close();
380 outputNames.push_back(compositeFastaFile); outputTypes["fasta"].push_back(compositeFastaFile);
381 outputNames.push_back(compositeMisMatchFile); outputTypes["report"].push_back(compositeMisMatchFile);
382 outputNames.push_back(compositeScrapFastaFile); outputTypes["fasta"].push_back(compositeScrapFastaFile);
385 for (int l = 0; l < filesToProcess.size(); l++) {
387 m->mothurOut("\n>>>>>\tProcessing " + filesToProcess[l][0][0] + " (file " + toString(l+1) + " of " + toString(filesToProcess.size()) + ")\t<<<<<\n");
389 vector<vector<string> > fastaFileNames;
391 string outputGroupFileName;
392 map<string, string> variables;
393 string thisOutputDir = outputDir;
394 if (outputDir == "") { thisOutputDir = m->hasPath(filesToProcess[l][0][0]); }
395 variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(filesToProcess[l][0][0]));
396 variables["[tag]"] = "";
397 if(oligosfile != ""){
398 createGroup = getOligos(fastaFileNames, variables["[filename]"]);
400 outputGroupFileName = getOutputFileName("group",variables);
401 outputNames.push_back(outputGroupFileName); outputTypes["group"].push_back(outputGroupFileName);
405 variables["[tag]"] = "trim";
406 string outFastaFile = getOutputFileName("fasta",variables);
407 variables["[tag]"] = "scrap";
408 string outScrapFastaFile = getOutputFileName("fasta",variables);
409 variables["[tag]"] = "";
410 string outMisMatchFile = getOutputFileName("report",variables);
411 outputNames.push_back(outFastaFile); outputTypes["fasta"].push_back(outFastaFile);
412 outputNames.push_back(outScrapFastaFile); outputTypes["fasta"].push_back(outScrapFastaFile);
413 outputNames.push_back(outMisMatchFile); outputTypes["report"].push_back(outMisMatchFile);
415 m->mothurOut("Making contigs...\n");
416 createProcesses(filesToProcess[l], outFastaFile, outScrapFastaFile, outMisMatchFile, fastaFileNames);
417 m->mothurOut("Done.\n");
419 //remove temp fasta and qual files
420 for (int i = 0; i < processors; i++) { for(int j = 0; j < filesToProcess[l][i].size(); j++) { m->mothurRemove(filesToProcess[l][i][j]); } }
422 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
425 map<string, string> uniqueFastaNames;// so we don't add the same groupfile multiple times
426 map<string, string>::iterator it;
427 set<string> namesToRemove;
428 for(int i=0;i<fastaFileNames.size();i++){
429 for(int j=0;j<fastaFileNames[0].size();j++){
430 if (fastaFileNames[i][j] != "") {
431 if (namesToRemove.count(fastaFileNames[i][j]) == 0) {
432 if(m->isBlank(fastaFileNames[i][j])){
433 m->mothurRemove(fastaFileNames[i][j]);
434 namesToRemove.insert(fastaFileNames[i][j]);
436 it = uniqueFastaNames.find(fastaFileNames[i][j]);
437 if (it == uniqueFastaNames.end()) {
438 uniqueFastaNames[fastaFileNames[i][j]] = barcodeNameVector[i];
446 //remove names for outputFileNames, just cleans up the output
447 vector<string> outputNames2;
448 for(int i = 0; i < outputNames.size(); i++) { if (namesToRemove.count(outputNames[i]) == 0) { outputNames2.push_back(outputNames[i]); } }
449 outputNames = outputNames2;
451 for (it = uniqueFastaNames.begin(); it != uniqueFastaNames.end(); it++) {
453 m->openInputFile(it->first, in);
456 string thisGroupName = thisOutputDir + m->getRootName(m->getSimpleName(it->first));
457 thisGroupName += getOutputFileName("group",variables); outputNames.push_back(thisGroupName); outputTypes["group"].push_back(thisGroupName);
458 m->openOutputFile(thisGroupName, out);
461 if (m->control_pressed) { break; }
463 Sequence currSeq(in); m->gobble(in);
464 out << currSeq.getName() << '\t' << it->second << endl;
473 m->openOutputFile(outputGroupFileName, outGroup);
474 for (map<string, string>::iterator itGroup = groupMap.begin(); itGroup != groupMap.end(); itGroup++) {
475 outGroup << itGroup->first << '\t' << itGroup->second << endl;
480 if (filesToProcess.size() > 1) { //merge into large combo files
484 m->openOutputFile(compositeGroupFile, outCGroup); outCGroup.close();
485 outputNames.push_back(compositeGroupFile); outputTypes["group"].push_back(compositeGroupFile);
487 m->appendFiles(outputGroupFileName, compositeGroupFile);
489 m->appendFiles(outMisMatchFile, compositeMisMatchFile);
490 m->appendFiles(outFastaFile, compositeFastaFile);
491 m->appendFiles(outScrapFastaFile, compositeScrapFastaFile);
494 m->mothurOut("It took " + toString(time(NULL) - start) + " secs to process " + toString(numReads) + " sequences.\n");
496 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
498 //output group counts
499 m->mothurOutEndLine();
501 if (groupCounts.size() != 0) { m->mothurOut("Group count: \n"); }
502 for (map<string, int>::iterator it = groupCounts.begin(); it != groupCounts.end(); it++) {
503 total += it->second; m->mothurOut(it->first + "\t" + toString(it->second)); m->mothurOutEndLine();
505 if (total != 0) { m->mothurOut("Total of all groups is " + toString(total)); m->mothurOutEndLine(); }
507 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
509 string currentFasta = "";
510 itTypes = outputTypes.find("fasta");
511 if (itTypes != outputTypes.end()) {
512 if ((itTypes->second).size() != 0) { currentFasta = (itTypes->second)[0]; m->setFastaFile(currentFasta); }
515 string currentGroup = "";
516 itTypes = outputTypes.find("group");
517 if (itTypes != outputTypes.end()) {
518 if ((itTypes->second).size() != 0) { currentGroup = (itTypes->second)[0]; m->setGroupFile(currentGroup); }
521 //output files created by command
522 m->mothurOutEndLine();
523 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
524 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
525 m->mothurOutEndLine();
529 catch(exception& e) {
530 m->errorOut(e, "MakeContigsCommand", "execute");
534 //**********************************************************************************************************************
535 vector< vector< vector<string> > > MakeContigsCommand::preProcessData(unsigned long int& numReads) {
537 vector< vector< vector<string> > > filesToProcess;
539 if (ffastqfile != "") { //reading one file
540 vector< vector<string> > files = readFastqFiles(numReads, ffastqfile, rfastqfile);
541 //adjust for really large processors or really small files
542 if (numReads == 0) { m->mothurOut("[ERROR]: no good reads.\n"); m->control_pressed = true; }
543 if (numReads < processors) {
544 for (int i = numReads; i < processors; i++) { for(int j = 0; j < files[i].size(); j++) { m->mothurRemove(files[i][j]); } files[i].clear(); }
545 files.resize(numReads);
546 processors = numReads;
548 filesToProcess.push_back(files);
549 }else if (file != "") { //reading multiple files
550 //return only valid pairs
551 vector< vector<string> > filePairsToProcess = readFileNames(file);
553 if (m->control_pressed) { return filesToProcess; }
555 if (filePairsToProcess.size() != 0) {
556 for (int i = 0; i < filePairsToProcess.size(); i++) {
558 if (m->control_pressed) { for (int l = 0; l < filesToProcess.size(); l++) { for (int k = 0; k < filesToProcess[l].size(); k++) { for(int j = 0; j < filesToProcess[l][k].size(); j++) { m->mothurRemove(filesToProcess[l][k][j]); } filesToProcess[l][k].clear(); } return filesToProcess; } }
560 unsigned long int thisFilesReads;
561 vector< vector<string> > files = readFastqFiles(thisFilesReads, filePairsToProcess[i][0], filePairsToProcess[i][1]);
563 //adjust for really large processors or really small files
564 if (thisFilesReads < processors) {
565 m->mothurOut("[ERROR]: " + filePairsToProcess[i][0] + " has less than " + toString(processors) + " good reads, skipping\n");
566 for (int k = 0; k < files.size(); k++) { for(int j = 0; j < files[k].size(); j++) { m->mothurRemove(files[k][j]); } files[k].clear(); }
568 filesToProcess.push_back(files);
569 numReads += thisFilesReads;
573 if (numReads == 0) { m->control_pressed = true; }
575 }else if (ffastafile != "") {
576 vector< vector<string> > files = readFastaFiles(numReads, ffastafile, rfastafile);
577 //adjust for really large processors or really small files
578 if (numReads == 0) { m->mothurOut("[ERROR]: no good reads.\n"); m->control_pressed = true; }
579 if (numReads < processors) {
580 for (int i = numReads; i < processors; i++) { for(int j = 0; j < files[i].size(); j++) { m->mothurRemove(files[i][j]); } files[i].clear(); }
581 files.resize(numReads);
582 processors = numReads;
584 filesToProcess.push_back(files);
585 }else { m->control_pressed = true; } //should not get here
587 return filesToProcess;
589 catch(exception& e) {
590 m->errorOut(e, "MakeContigsCommand", "preProcessData");
594 //**********************************************************************************************************************
595 int MakeContigsCommand::createProcesses(vector< vector<string> > files, string outputFasta, string outputScrapFasta, string outputMisMatches, vector<vector<string> > fastaFileNames) {
598 vector<int> processIDS;
599 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
602 //loop through and create all the processes you want
603 while (process != processors-1) {
607 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
610 vector<vector<string> > tempFASTAFileNames = fastaFileNames;
615 for(int i=0;i<tempFASTAFileNames.size();i++){
616 for(int j=0;j<tempFASTAFileNames[i].size();j++){
617 if (tempFASTAFileNames[i][j] != "") {
618 tempFASTAFileNames[i][j] += toString(getpid()) + ".temp";
619 m->openOutputFile(tempFASTAFileNames[i][j], temp); temp.close();
625 num = driver(files[process],
626 outputFasta + toString(getpid()) + ".temp",
627 outputScrapFasta + toString(getpid()) + ".temp",
628 outputMisMatches + toString(getpid()) + ".temp",
629 tempFASTAFileNames, process);
631 //pass groupCounts to parent
633 string tempFile = toString(getpid()) + ".num.temp";
634 m->openOutputFile(tempFile, out);
637 out << groupCounts.size() << endl;
639 for (map<string, int>::iterator it = groupCounts.begin(); it != groupCounts.end(); it++) {
640 out << it->first << '\t' << it->second << endl;
643 out << groupMap.size() << endl;
644 for (map<string, string>::iterator it = groupMap.begin(); it != groupMap.end(); it++) {
645 out << it->first << '\t' << it->second << endl;
652 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
653 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
659 m->openOutputFile(outputFasta, temp); temp.close();
660 m->openOutputFile(outputScrapFasta, temp); temp.close();
663 num = driver(files[processors-1], outputFasta, outputScrapFasta, outputMisMatches, fastaFileNames, processors-1);
665 //force parent to wait until all the processes are done
666 for (int i=0;i<processIDS.size();i++) {
667 int temp = processIDS[i];
671 for (int i = 0; i < processIDS.size(); i++) {
673 string tempFile = toString(processIDS[i]) + ".num.temp";
674 m->openInputFile(tempFile, in);
676 in >> tempNum; num += tempNum; m->gobble(in);
680 in >> tempNum; m->gobble(in);
683 for (int j = 0; j < tempNum; j++) {
685 in >> group >> groupNum; m->gobble(in);
687 map<string, int>::iterator it = groupCounts.find(group);
688 if (it == groupCounts.end()) { groupCounts[group] = groupNum; }
689 else { groupCounts[it->first] += groupNum; }
692 in >> tempNum; m->gobble(in);
694 for (int j = 0; j < tempNum; j++) {
695 string group, seqName;
696 in >> seqName >> group; m->gobble(in);
698 map<string, string>::iterator it = groupMap.find(seqName);
699 if (it == groupMap.end()) { groupMap[seqName] = group; }
700 else { m->mothurOut("[ERROR]: " + seqName + " is in your fasta file more than once. Sequence names must be unique. please correct.\n"); }
704 in.close(); m->mothurRemove(tempFile);
708 //////////////////////////////////////////////////////////////////////////////////////////////////////
709 //Windows version shared memory, so be careful when passing variables through the contigsData struct.
710 //Above fork() will clone, so memory is separate, but that's not the case with windows,
711 //////////////////////////////////////////////////////////////////////////////////////////////////////
713 vector<contigsData*> pDataArray;
714 DWORD dwThreadIdArray[processors-1];
715 HANDLE hThreadArray[processors-1];
717 //Create processor worker threads.
718 for( int h=0; h<processors-1; h++ ){
719 string extension = "";
720 if (h != 0) { extension = toString(h) + ".temp"; processIDS.push_back(h); }
721 vector<vector<string> > tempFASTAFileNames = fastaFileNames;
726 for(int i=0;i<tempFASTAFileNames.size();i++){
727 for(int j=0;j<tempFASTAFileNames[i].size();j++){
728 if (tempFASTAFileNames[i][j] != "") {
729 tempFASTAFileNames[i][j] += extension;
730 m->openOutputFile(tempFASTAFileNames[i][j], temp); temp.close();
737 contigsData* tempcontig = new contigsData(files[h], (outputFasta + extension), (outputScrapFasta + extension), (outputMisMatches + extension), align, m, match, misMatch, gapOpen, gapExtend, insert, deltaq, barcodes, primers, tempFASTAFileNames, barcodeNameVector, primerNameVector, pdiffs, bdiffs, tdiffs, createGroup, allFiles, h);
738 pDataArray.push_back(tempcontig);
740 hThreadArray[h] = CreateThread(NULL, 0, MyContigsThreadFunction, pDataArray[h], 0, &dwThreadIdArray[h]);
743 vector<vector<string> > tempFASTAFileNames = fastaFileNames;
747 string extension = toString(processors-1) + ".temp";
749 for(int i=0;i<tempFASTAFileNames.size();i++){
750 for(int j=0;j<tempFASTAFileNames[i].size();j++){
751 if (tempFASTAFileNames[i][j] != "") {
752 tempFASTAFileNames[i][j] += extension;
753 m->openOutputFile(tempFASTAFileNames[i][j], temp); temp.close();
761 m->openOutputFile(outputFasta, temp); temp.close();
762 m->openOutputFile(outputScrapFasta, temp); temp.close();
765 processIDS.push_back(processors-1);
766 num = driver(files[processors-1], (outputFasta+ toString(processors-1) + ".temp"), (outputScrapFasta+ toString(processors-1) + ".temp"), (outputMisMatches+ toString(processors-1) + ".temp"), tempFASTAFileNames, processors-1);
768 //Wait until all threads have terminated.
769 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
771 //Close all thread handles and free memory allocations.
772 for(int i=0; i < pDataArray.size(); i++){
773 num += pDataArray[i]->count;
774 if (!pDataArray[i]->done) {
775 m->mothurOut("[ERROR]: process " + toString(i) + " only processed " + toString(pDataArray[i]->count) + " of sequences assigned to it, quitting. \n"); m->control_pressed = true;
777 for (map<string, int>::iterator it = pDataArray[i]->groupCounts.begin(); it != pDataArray[i]->groupCounts.end(); it++) {
778 map<string, int>::iterator it2 = groupCounts.find(it->first);
779 if (it2 == groupCounts.end()) { groupCounts[it->first] = it->second; }
780 else { groupCounts[it->first] += it->second; }
782 for (map<string, string>::iterator it = pDataArray[i]->groupMap.begin(); it != pDataArray[i]->groupMap.end(); it++) {
783 map<string, string>::iterator it2 = groupMap.find(it->first);
784 if (it2 == groupMap.end()) { groupMap[it->first] = it->second; }
785 else { m->mothurOut("[ERROR]: " + it->first + " is in your fasta file more than once. Sequence names must be unique. please correct.\n"); }
787 CloseHandle(hThreadArray[i]);
788 delete pDataArray[i];
793 for (int i = 0; i < processIDS.size(); i++) {
794 m->appendFiles((outputFasta + toString(processIDS[i]) + ".temp"), outputFasta);
795 m->mothurRemove((outputFasta + toString(processIDS[i]) + ".temp"));
797 m->appendFiles((outputScrapFasta + toString(processIDS[i]) + ".temp"), outputScrapFasta);
798 m->mothurRemove((outputScrapFasta + toString(processIDS[i]) + ".temp"));
800 m->appendFiles((outputMisMatches + toString(processIDS[i]) + ".temp"), outputMisMatches);
801 m->mothurRemove((outputMisMatches + toString(processIDS[i]) + ".temp"));
804 for(int j=0;j<fastaFileNames.size();j++){
805 for(int k=0;k<fastaFileNames[j].size();k++){
806 if (fastaFileNames[j][k] != "") {
807 m->appendFiles((fastaFileNames[j][k] + toString(processIDS[i]) + ".temp"), fastaFileNames[j][k]);
808 m->mothurRemove((fastaFileNames[j][k] + toString(processIDS[i]) + ".temp"));
817 catch(exception& e) {
818 m->errorOut(e, "MakeContigsCommand", "createProcesses");
822 //**********************************************************************************************************************
823 int MakeContigsCommand::driver(vector<string> files, string outputFasta, string outputScrapFasta, string outputMisMatches, vector<vector<string> > fastaFileNames, int process){
826 Alignment* alignment;
827 if(align == "gotoh") { alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase); }
828 else if(align == "needleman") { alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase); }
831 string thisffastafile = files[0];
832 string thisfqualfile = files[1];
833 string thisrfastafile = files[2];
834 string thisrqualfile = files[3];
836 if (m->debug) { m->mothurOut("[DEBUG]: ffasta = " + thisffastafile + ".\n[DEBUG]: fqual = " + thisfqualfile + ".\n[DEBUG]: rfasta = " + thisrfastafile + ".\n[DEBUG]: rqual = " + thisrqualfile + ".\n"); }
838 ifstream inFFasta, inRFasta, inFQual, inRQual;
839 ofstream outFasta, outMisMatch, outScrapFasta;
840 m->openInputFile(thisffastafile, inFFasta);
841 m->openInputFile(thisrfastafile, inRFasta);
842 if (thisfqualfile != "") {
843 m->openInputFile(thisfqualfile, inFQual);
844 m->openInputFile(thisrqualfile, inRQual);
846 m->openOutputFile(outputFasta, outFasta);
847 m->openOutputFile(outputScrapFasta, outScrapFasta);
848 m->openOutputFile(outputMisMatches, outMisMatch);
849 if (process == 0) { outMisMatch << "Name\tLength\tOverlap_Length\tOverlap_Start\tOverlap_End\tMisMatches\tNum_Ns\n"; }
851 TrimOligos trimOligos(pdiffs, bdiffs, 0, 0, primers, barcodes);
853 while ((!inFFasta.eof()) && (!inRFasta.eof())) {
855 if (m->control_pressed) { break; }
858 string trashCode = "";
859 int currentSeqsDiffs = 0;
861 //read seqs and quality info
862 Sequence fSeq(inFFasta); m->gobble(inFFasta);
863 Sequence rSeq(inRFasta); m->gobble(inRFasta);
864 QualityScores* fQual = NULL; QualityScores* rQual = NULL;
865 if (thisfqualfile != "") {
866 fQual = new QualityScores(inFQual); m->gobble(inFQual);
867 rQual = new QualityScores(inRQual); m->gobble(inRQual);
870 int barcodeIndex = 0;
873 if(barcodes.size() != 0){
874 if (thisfqualfile != "") {
875 success = trimOligos.stripBarcode(fSeq, rSeq, *fQual, *rQual, barcodeIndex);
877 success = trimOligos.stripBarcode(fSeq, rSeq, barcodeIndex);
879 if(success > bdiffs) { trashCode += 'b'; }
880 else{ currentSeqsDiffs += success; }
883 if(primers.size() != 0){
884 if (thisfqualfile != "") {
885 success = trimOligos.stripForward(fSeq, rSeq, *fQual, *rQual, primerIndex);
887 success = trimOligos.stripForward(fSeq, rSeq, primerIndex);
889 if(success > pdiffs) { trashCode += 'f'; }
890 else{ currentSeqsDiffs += success; }
893 if (currentSeqsDiffs > tdiffs) { trashCode += 't'; }
895 //flip the reverse reads
896 rSeq.reverseComplement();
897 if (thisfqualfile != "") { rQual->flipQScores(); }
900 alignment->align(fSeq.getUnaligned(), rSeq.getUnaligned());
901 map<int, int> ABaseMap = alignment->getSeqAAlnBaseMap();
902 map<int, int> BBaseMap = alignment->getSeqBAlnBaseMap();
903 fSeq.setAligned(alignment->getSeqAAln());
904 rSeq.setAligned(alignment->getSeqBAln());
905 int length = fSeq.getAligned().length();
907 //traverse alignments merging into one contiguous seq
909 int numMismatches = 0;
910 string seq1 = fSeq.getAligned();
911 string seq2 = rSeq.getAligned();
912 vector<int> scores1, scores2;
913 if (thisfqualfile != "") {
914 scores1 = fQual->getQualityScores();
915 scores2 = rQual->getQualityScores();
916 delete fQual; delete rQual;
919 // if (num < 5) { cout << fSeq.getStartPos() << '\t' << fSeq.getEndPos() << '\t' << rSeq.getStartPos() << '\t' << rSeq.getEndPos() << endl; }
920 int overlapStart = fSeq.getStartPos();
921 int seq2Start = rSeq.getStartPos();
922 //bigger of the 2 starting positions is the location of the overlapping start
923 if (overlapStart < seq2Start) { //seq2 starts later so take from 0 to seq2Start from seq1
924 overlapStart = seq2Start;
925 for (int i = 0; i < overlapStart; i++) { contig += seq1[i]; }
926 }else { //seq1 starts later so take from 0 to overlapStart from seq2
927 for (int i = 0; i < overlapStart; i++) { contig += seq2[i]; }
930 int seq1End = fSeq.getEndPos();
931 int seq2End = rSeq.getEndPos();
932 int overlapEnd = seq1End;
933 if (seq2End < overlapEnd) { overlapEnd = seq2End; } //smallest end position is where overlapping ends
935 int oStart = contig.length();
936 for (int i = overlapStart; i < overlapEnd; i++) {
937 if (seq1[i] == seq2[i]) { //match, add base and choose highest score
939 }else if (((seq1[i] == '.') || (seq1[i] == '-')) && ((seq2[i] != '-') && (seq2[i] != '.'))) { //seq1 is a gap and seq2 is a base, choose seq2, unless quality score for base is below insert. In that case eliminate base
940 if (thisfqualfile != "") {
941 if (scores2[BBaseMap[i]] < insert) { } //
942 else { contig += seq2[i]; }
943 }else { contig += seq2[i]; } //with no quality info, then we keep it?
944 }else if (((seq2[i] == '.') || (seq2[i] == '-')) && ((seq1[i] != '-') && (seq1[i] != '.'))) { //seq2 is a gap and seq1 is a base, choose seq1, unless quality score for base is below insert. In that case eliminate base
945 if (thisfqualfile != "") {
946 if (scores1[ABaseMap[i]] < insert) { } //
947 else { contig += seq1[i]; }
948 }else { contig += seq1[i]; } //with no quality info, then we keep it?
949 }else if (((seq1[i] != '-') && (seq1[i] != '.')) && ((seq2[i] != '-') && (seq2[i] != '.'))) { //both bases choose one with better quality
950 if (thisfqualfile != "") {
951 if (abs(scores1[ABaseMap[i]] - scores2[BBaseMap[i]]) >= deltaq) { //is the difference in qual scores >= deltaq, if yes choose base with higher score
953 if (scores1[ABaseMap[i]] < scores2[BBaseMap[i]]) { c = seq2[i]; }
955 }else { //if no, base becomes n
959 }else { numMismatches++; } //cant decide, so eliminate and mark as mismatch
960 }else { //should never get here
961 m->mothurOut("[ERROR]: case I didn't think of seq1 = " + toString(seq1[i]) + " and seq2 = " + toString(seq2[i]) + "\n");
964 int oend = contig.length();
965 if (seq1End < seq2End) { //seq1 ends before seq2 so take from overlap to length from seq2
966 for (int i = overlapEnd; i < length; i++) { contig += seq2[i]; }
967 }else { //seq2 ends before seq1 so take from overlap to length from seq1
968 for (int i = overlapEnd; i < length; i++) { contig += seq1[i]; }
971 if(trashCode.length() == 0){
974 if (m->debug) { m->mothurOut(fSeq.getName()); }
977 if(barcodes.size() != 0){
978 string thisGroup = barcodeNameVector[barcodeIndex];
979 if (primers.size() != 0) {
980 if (primerNameVector[primerIndex] != "") {
981 if(thisGroup != "") {
982 thisGroup += "." + primerNameVector[primerIndex];
984 thisGroup = primerNameVector[primerIndex];
989 if (m->debug) { m->mothurOut(", group= " + thisGroup + "\n"); }
991 int pos = thisGroup.find("ignore");
992 if (pos == string::npos) {
993 groupMap[fSeq.getName()] = thisGroup;
995 map<string, int>::iterator it = groupCounts.find(thisGroup);
996 if (it == groupCounts.end()) { groupCounts[thisGroup] = 1; }
997 else { groupCounts[it->first] ++; }
998 }else { ignore = true; }
1002 if (m->debug) { m->mothurOut("\n"); }
1004 if(allFiles && !ignore){
1006 m->openOutputFileAppend(fastaFileNames[barcodeIndex][primerIndex], output);
1007 output << ">" << fSeq.getName() << endl << contig << endl;
1012 outFasta << ">" << fSeq.getName() << endl << contig << endl;
1014 for (int i = 0; i < contig.length(); i++) { if (contig[i] == 'N') { numNs++; } }
1015 outMisMatch << fSeq.getName() << '\t' << contig.length() << '\t' << (oend-oStart) << '\t' << oStart << '\t' << oend << '\t' << numMismatches << '\t' << numNs << endl;
1018 outScrapFasta << ">" << fSeq.getName() << " | " << trashCode << endl << contig << endl;
1023 if((num) % 1000 == 0){ m->mothurOut(toString(num)); m->mothurOutEndLine(); }
1027 if((num) % 1000 != 0){ m->mothurOut(toString(num)); m->mothurOutEndLine(); }
1032 outScrapFasta.close();
1033 outMisMatch.close();
1034 if (thisfqualfile != "") {
1040 if (m->control_pressed) { m->mothurRemove(outputFasta); m->mothurRemove(outputScrapFasta);m->mothurRemove(outputMisMatches); }
1044 catch(exception& e) {
1045 m->errorOut(e, "MakeContigsCommand", "driver");
1049 //**********************************************************************************************************************
1050 vector< vector<string> > MakeContigsCommand::readFastqFiles(unsigned long int& count, string ffastq, string rfastq){
1052 vector< vector<string> > files;
1053 //maps processors number to file pointer
1054 map<int, vector<ofstream*> > tempfiles; //tempfiles[0] = forwardFasta, [1] = forwardQual, [2] = reverseFasta, [3] = reverseQual
1055 map<int, vector<ofstream*> >::iterator it;
1057 //create files to write to
1058 for (int i = 0; i < processors; i++) {
1059 vector<ofstream*> temp;
1060 ofstream* outFF = new ofstream; temp.push_back(outFF);
1061 ofstream* outFQ = new ofstream; temp.push_back(outFQ);
1062 ofstream* outRF = new ofstream; temp.push_back(outRF);
1063 ofstream* outRQ = new ofstream; temp.push_back(outRQ);
1064 tempfiles[i] = temp;
1066 vector<string> names;
1067 string thisOutputDir = outputDir;
1068 if (outputDir == "") { thisOutputDir = m->hasPath(ffastq); }
1069 string ffastafilename = thisOutputDir + m->getRootName(m->getSimpleName(ffastq)) + toString(i) + "ffastatemp";
1070 string rfastafilename = thisOutputDir + m->getRootName(m->getSimpleName(rfastq)) + toString(i) + "rfastatemp";
1071 string fqualfilename = thisOutputDir + m->getRootName(m->getSimpleName(ffastq)) + toString(i) + "fqualtemp";
1072 string rqualfilename = thisOutputDir + m->getRootName(m->getSimpleName(rfastq)) + toString(i) + "rqualtemp";
1073 names.push_back(ffastafilename); names.push_back(fqualfilename);
1074 names.push_back(rfastafilename); names.push_back(rqualfilename);
1075 files.push_back(names);
1077 m->openOutputFile(ffastafilename, *outFF);
1078 m->openOutputFile(rfastafilename, *outRF);
1079 m->openOutputFile(fqualfilename, *outFQ);
1080 m->openOutputFile(rqualfilename, *outRQ);
1083 if (m->control_pressed) {
1084 //close files, delete ofstreams
1085 for (it = tempfiles.begin(); it!=tempfiles.end(); it++) { for (int i = 0; i < (it->second).size(); i++) { (*(it->second)[i]).close(); delete (it->second)[i]; } }
1087 for (int i = 0; i < files.size(); i++) {
1088 for(int j = 0; j < files[i].size(); j++) { m->mothurRemove(files[i][j]); }
1093 m->openInputFile(ffastq, inForward);
1096 m->openInputFile(rfastq, inReverse);
1099 map<string, fastqRead> uniques;
1100 map<string, fastqRead>::iterator itUniques;
1101 while ((!inForward.eof()) || (!inReverse.eof())) {
1103 if (m->control_pressed) { for (it = tempfiles.begin(); it!=tempfiles.end(); it++) { for (int i = 0; i < (it->second).size(); i++) { (*(it->second)[i]).close(); delete (it->second)[i]; } } for (int i = 0; i < files.size(); i++) { for(int j = 0; j < files[i].size(); j++) { m->mothurRemove(files[i][j]); } } inForward.close(); inReverse.close(); return files; }
1105 //get a read from forward and reverse fastq files
1106 bool ignoref, ignorer;
1107 fastqRead thisFread, thisRread;
1108 if (!inForward.eof()) { thisFread = readFastq(inForward, ignoref); }
1109 else { ignoref = true; }
1110 if (!inReverse.eof()) { thisRread = readFastq(inReverse, ignorer); }
1111 else { ignorer = true; }
1113 vector<pairFastqRead> reads = getReads(ignoref, ignorer, thisFread, thisRread, uniques);
1115 for (int i = 0; i < reads.size(); i++) {
1116 fastqRead fread = reads[i].forward;
1117 fastqRead rread = reads[i].reverse;
1119 if (m->debug) { m->mothurOut(toString(count) + '\t' + fread.name + '\t' + rread.name + '\n'); }
1121 //if (checkReads(fread, rread, ffastq, rfastq)) {
1122 if (m->control_pressed) { for (it = tempfiles.begin(); it!=tempfiles.end(); it++) { for (int i = 0; i < (it->second).size(); i++) { (*(it->second)[i]).close(); delete (it->second)[i]; } } for (int i = 0; i < files.size(); i++) { for(int j = 0; j < files[i].size(); j++) { m->mothurRemove(files[i][j]); } } inForward.close(); inReverse.close(); return files; }
1124 //if the reads are okay write to output files
1125 int process = count % processors;
1127 *(tempfiles[process][0]) << ">" << fread.name << endl << fread.sequence << endl;
1128 *(tempfiles[process][1]) << ">" << fread.name << endl;
1129 for (int i = 0; i < fread.scores.size(); i++) { *(tempfiles[process][1]) << fread.scores[i] << " "; }
1130 *(tempfiles[process][1]) << endl;
1131 *(tempfiles[process][2]) << ">" << rread.name << endl << rread.sequence << endl;
1132 *(tempfiles[process][3]) << ">" << rread.name << endl;
1133 for (int i = 0; i < rread.scores.size(); i++) { *(tempfiles[process][3]) << rread.scores[i] << " "; }
1134 *(tempfiles[process][3]) << endl;
1139 if((count) % 10000 == 0){ m->mothurOut(toString(count)); m->mothurOutEndLine(); }
1144 if((count) % 10000 != 0){ m->mothurOut(toString(count)); m->mothurOutEndLine(); }
1146 if (uniques.size() != 0) {
1147 for (itUniques = uniques.begin(); itUniques != uniques.end(); itUniques++) {
1148 m->mothurOut("[WARNING]: did not find paired read for " + itUniques->first + ", ignoring.\n");
1150 m->mothurOutEndLine();
1153 //close files, delete ofstreams
1154 for (it = tempfiles.begin(); it!=tempfiles.end(); it++) { for (int i = 0; i < (it->second).size(); i++) { (*(it->second)[i]).close(); delete (it->second)[i]; } }
1160 catch(exception& e) {
1161 m->errorOut(e, "MakeContigsCommand", "readFastqFiles");
1165 //**********************************************************************************************************************
1166 vector< vector<string> > MakeContigsCommand::readFastaFiles(unsigned long int& count, string ffasta, string rfasta){
1168 vector< vector<string> > files;
1169 //maps processors number to file pointer
1170 map<int, vector<ofstream*> > tempfiles; //tempfiles[0] = forwardFasta, [1] = forwardQual, [2] = reverseFasta, [3] = reverseQual
1171 map<int, vector<ofstream*> >::iterator it;
1173 //create files to write to
1174 for (int i = 0; i < processors; i++) {
1175 vector<ofstream*> temp;
1176 ofstream* outFF = new ofstream; temp.push_back(outFF);
1177 ofstream* outFQ = new ofstream; temp.push_back(outFQ);
1178 ofstream* outRF = new ofstream; temp.push_back(outRF);
1179 ofstream* outRQ = new ofstream; temp.push_back(outRQ);
1180 tempfiles[i] = temp;
1182 vector<string> names;
1183 string thisOutputDir = outputDir;
1184 if (outputDir == "") { thisOutputDir = m->hasPath(ffasta); }
1185 string ffastafilename = thisOutputDir + m->getRootName(m->getSimpleName(ffasta)) + toString(i) + "ffastatemp";
1186 string rfastafilename = thisOutputDir + m->getRootName(m->getSimpleName(rfasta)) + toString(i) + "rfastatemp";
1187 string fqualfilename = "";
1188 if (fqualfile != "") { fqualfilename = thisOutputDir + m->getRootName(m->getSimpleName(fqualfile)) + toString(i) + "fqual.temp"; m->openOutputFile(fqualfilename, *outFQ); }
1189 string rqualfilename = "";
1190 if (rqualfile != "") { rqualfilename = thisOutputDir + m->getRootName(m->getSimpleName(rqualfile)) + toString(i) + "rqual.temp"; m->openOutputFile(rqualfilename, *outRQ); }
1191 names.push_back(ffastafilename); names.push_back(fqualfilename);
1192 names.push_back(rfastafilename); names.push_back(rqualfilename);
1193 files.push_back(names);
1195 m->openOutputFile(ffastafilename, *outFF);
1196 m->openOutputFile(rfastafilename, *outRF);
1199 if (m->control_pressed) {
1200 //close files, delete ofstreams
1201 for (it = tempfiles.begin(); it!=tempfiles.end(); it++) { for (int i = 0; i < (it->second).size(); i++) { (*(it->second)[i]).close(); delete (it->second)[i]; } }
1203 for (int i = 0; i < files.size(); i++) {
1204 for(int j = 0; j < files[i].size(); j++) { m->mothurRemove(files[i][j]); }
1208 ifstream inForwardFasta;
1209 m->openInputFile(ffasta, inForwardFasta);
1211 ifstream inReverseFasta;
1212 m->openInputFile(rfasta, inReverseFasta);
1214 ifstream inForwardQual, inReverseQual;
1215 if (fqualfile != "") { m->openInputFile(fqualfile, inForwardQual); m->openInputFile(rqualfile, inReverseQual); }
1218 map<string, fastqRead> uniques;
1219 map<string, fastqRead>::iterator itUniques;
1220 while ((!inForwardFasta.eof()) || (!inReverseFasta.eof())) {
1222 if (m->control_pressed) { for (it = tempfiles.begin(); it!=tempfiles.end(); it++) { for (int i = 0; i < (it->second).size(); i++) { (*(it->second)[i]).close(); delete (it->second)[i]; } } for (int i = 0; i < files.size(); i++) { for(int j = 0; j < files[i].size(); j++) { m->mothurRemove(files[i][j]); } } inReverseFasta.close(); inForwardFasta.close(); if (fqualfile != "") { inReverseQual.close(); inReverseQual.close(); } return files; }
1224 //get a reads from forward and reverse fasta files
1225 bool ignoref, ignorer;
1226 fastqRead thisFread, thisRread;
1227 if (!inForwardFasta.eof()) {
1229 Sequence temp(inForwardFasta); m->gobble(inForwardFasta);
1230 thisFread.name = temp.getName();
1231 thisFread.sequence = temp.getUnaligned();
1232 }else { ignoref = true; }
1233 if (!inReverseFasta.eof()) {
1235 Sequence temp(inReverseFasta); m->gobble(inReverseFasta);
1236 thisRread.name = temp.getName();
1237 thisRread.sequence = temp.getUnaligned();
1238 }else { ignorer = true; }
1240 //get qual reads if given
1241 if (fqualfile != "") {
1242 if (!inForwardQual.eof() && !ignoref) {
1243 QualityScores temp(inForwardQual); m->gobble(inForwardQual);
1244 //if forward files dont match ignore read
1245 if (thisFread.name != temp.getName()) { ignoref = true; }
1246 else { thisFread.scores = temp.getQualityScores(); }
1247 }else { ignoref = true; }
1248 if (!inReverseQual.eof() && !ignorer) {
1249 QualityScores temp(inReverseQual); m->gobble(inReverseQual);
1250 //if reverse files dont match ignore read
1251 if (thisRread.name != temp.getName()) { ignorer = true; }
1252 else { thisRread.scores = temp.getQualityScores(); }
1253 }else { ignorer = true; }
1256 vector<pairFastqRead> reads = getReads(ignoref, ignorer, thisFread, thisRread, uniques);
1258 for (int i = 0; i < reads.size(); i++) {
1259 fastqRead fread = reads[i].forward;
1260 fastqRead rread = reads[i].reverse;
1262 if (m->debug) { m->mothurOut(toString(count) + '\t' + fread.name + '\t' + rread.name + '\n'); }
1264 // if (checkReads(fread, rread, ffasta, rfasta)) {
1265 if (m->control_pressed) { for (it = tempfiles.begin(); it!=tempfiles.end(); it++) { for (int i = 0; i < (it->second).size(); i++) { (*(it->second)[i]).close(); delete (it->second)[i]; } } for (int i = 0; i < files.size(); i++) { for(int j = 0; j < files[i].size(); j++) { m->mothurRemove(files[i][j]); } } inReverseFasta.close(); inForwardFasta.close(); if (fqualfile != "") { inReverseQual.close(); inReverseQual.close(); } return files; }
1267 //if the reads are okay write to output files
1268 int process = count % processors;
1270 *(tempfiles[process][0]) << ">" << fread.name << endl << fread.sequence << endl;
1271 *(tempfiles[process][2]) << ">" << rread.name << endl << rread.sequence << endl;
1272 if (fqualfile != "") { //if you have quality info, print it
1273 *(tempfiles[process][1]) << ">" << fread.name << endl;
1274 for (int i = 0; i < fread.scores.size(); i++) { *(tempfiles[process][1]) << fread.scores[i] << " "; }
1275 *(tempfiles[process][1]) << endl;
1276 *(tempfiles[process][3]) << ">" << rread.name << endl;
1277 for (int i = 0; i < rread.scores.size(); i++) { *(tempfiles[process][3]) << rread.scores[i] << " "; }
1278 *(tempfiles[process][3]) << endl;
1283 if((count) % 10000 == 0){ m->mothurOut(toString(count)); m->mothurOutEndLine(); }
1288 if((count) % 10000 != 0){ m->mothurOut(toString(count)); m->mothurOutEndLine(); }
1290 if (uniques.size() != 0) {
1291 for (itUniques = uniques.begin(); itUniques != uniques.end(); itUniques++) {
1292 m->mothurOut("[WARNING]: did not find paired read for " + itUniques->first + ", ignoring.\n");
1294 m->mothurOutEndLine();
1297 //close files, delete ofstreams
1298 for (it = tempfiles.begin(); it!=tempfiles.end(); it++) { for (int i = 0; i < (it->second).size(); i++) { (*(it->second)[i]).close(); delete (it->second)[i]; } }
1299 inReverseFasta.close();
1300 inForwardFasta.close();
1301 if (fqualfile != "") { inReverseQual.close(); inReverseQual.close(); }
1305 catch(exception& e) {
1306 m->errorOut(e, "MakeContigsCommand", "readFastaFiles");
1310 //**********************************************************************************************************************
1311 vector<pairFastqRead> MakeContigsCommand::getReads(bool ignoref, bool ignorer, fastqRead forward, fastqRead reverse, map<string, fastqRead>& uniques){
1313 vector<pairFastqRead> reads;
1314 map<string, fastqRead>::iterator itUniques;
1316 if (!ignoref && !ignorer) {
1317 if (forward.name == reverse.name) {
1318 pairFastqRead temp(forward, reverse);
1319 reads.push_back(temp);
1322 //if no match are the names only different by 1 and 2?
1323 string tempFRead = forward.name.substr(0, forward.name.length()-1);
1324 string tempRRead = reverse.name.substr(0, reverse.name.length()-1);
1325 if (tempFRead == tempRRead) {
1326 if ((forward.name[forward.name.length()-1] == '1') && (reverse.name[reverse.name.length()-1] == '2')) {
1327 forward.name = tempFRead;
1328 reverse.name = tempRRead;
1329 pairFastqRead temp(forward, reverse);
1330 reads.push_back(temp);
1336 //look for forward pair
1337 itUniques = uniques.find(forward.name);
1338 if (itUniques != uniques.end()) { //we have the pair for this read
1339 pairFastqRead temp(forward, itUniques->second);
1340 reads.push_back(temp);
1341 uniques.erase(itUniques);
1342 }else { //save this read for later
1343 uniques[forward.name] = forward;
1346 //look for reverse pair
1347 itUniques = uniques.find(reverse.name);
1348 if (itUniques != uniques.end()) { //we have the pair for this read
1349 pairFastqRead temp(itUniques->second, reverse);
1350 reads.push_back(temp);
1351 uniques.erase(itUniques);
1352 }else { //save this read for later
1353 uniques[reverse.name] = reverse;
1358 }else if (!ignoref && ignorer) { //ignore reverse keep forward
1359 //look for forward pair
1360 itUniques = uniques.find(forward.name);
1361 if (itUniques != uniques.end()) { //we have the pair for this read
1362 pairFastqRead temp(forward, itUniques->second);
1363 reads.push_back(temp);
1364 uniques.erase(itUniques);
1365 }else { //save this read for later
1366 uniques[forward.name] = forward;
1369 }else if (ignoref && !ignorer) { //ignore forward keep reverse
1370 //look for reverse pair
1371 itUniques = uniques.find(reverse.name);
1372 if (itUniques != uniques.end()) { //we have the pair for this read
1373 pairFastqRead temp(itUniques->second, reverse);
1374 reads.push_back(temp);
1375 uniques.erase(itUniques);
1376 }else { //save this read for later
1377 uniques[reverse.name] = reverse;
1379 }//else ignore both and do nothing
1383 catch(exception& e) {
1384 m->errorOut(e, "MakeContigsCommand", "readFastqFiles");
1388 //**********************************************************************************************************************
1389 fastqRead MakeContigsCommand::readFastq(ifstream& in, bool& ignore){
1395 //read sequence name
1396 string line = m->getline(in); m->gobble(in);
1397 vector<string> pieces = m->splitWhiteSpace(line);
1398 string name = ""; if (pieces.size() != 0) { name = pieces[0]; }
1399 if (name == "") { m->mothurOut("[WARNING]: Blank fasta name, ignoring read."); m->mothurOutEndLine(); ignore=true; }
1400 else if (name[0] != '@') { m->mothurOut("[WARNING]: reading " + name + " expected a name with @ as a leading character, ignoring read."); m->mothurOutEndLine(); ignore=true; }
1401 else { name = name.substr(1); }
1404 string sequence = m->getline(in); m->gobble(in);
1405 if (sequence == "") { m->mothurOut("[WARNING]: missing sequence for " + name + ", ignoring."); ignore=true; }
1407 //read sequence name
1408 line = m->getline(in); m->gobble(in);
1409 pieces = m->splitWhiteSpace(line);
1410 string name2 = ""; if (pieces.size() != 0) { name2 = pieces[0]; }
1411 if (name2 == "") { m->mothurOut("[WARNING]: expected a name with + as a leading character, ignoring."); ignore=true; }
1412 else if (name2[0] != '+') { m->mothurOut("[WARNING]: reading " + name2 + " expected a name with + as a leading character, ignoring."); ignore=true; }
1413 else { name2 = name2.substr(1); if (name2 == "") { name2 = name; } }
1415 //read quality scores
1416 string quality = m->getline(in); m->gobble(in);
1417 if (quality == "") { m->mothurOut("[WARNING]: missing quality for " + name2 + ", ignoring."); ignore=true; }
1419 //sanity check sequence length and number of quality scores match
1420 if (name2 != "") { if (name != name2) { m->mothurOut("[WARNING]: names do not match. read " + name + " for fasta and " + name2 + " for quality, ignoring."); ignore=true; } }
1421 if (quality.length() != sequence.length()) { m->mothurOut("[WARNING]: Lengths do not match for sequence " + name + ". Read " + toString(sequence.length()) + " characters for fasta and " + toString(quality.length()) + " characters for quality scores, ignoring read."); ignore=true; }
1423 vector<int> qualScores = convertQual(quality);
1426 read.sequence = sequence;
1427 read.scores = qualScores;
1431 catch(exception& e) {
1432 m->errorOut(e, "MakeContigsCommand", "readFastq");
1436 /**********************************************************************************************************************
1437 bool MakeContigsCommand::checkReads(fastqRead& forward, fastqRead& reverse, string ffile, string rfile){
1441 //do sequence lengths match
1442 if (forward.sequence.length() != reverse.sequence.length()) {
1443 m->mothurOut("[WARNING]: For sequence " + forward.name + " I read a sequence of length " + toString(forward.sequence.length()) + " from " + ffile + ", but read a sequence of length " + toString(reverse.sequence.length()) + " from " + rfile + ", ignoring.\n");
1447 //do number of qual scores match
1448 if (forward.scores.size() != reverse.scores.size()) {
1449 m->mothurOut("[WARNING]: For sequence " + forward.name + " I read " + toString(forward.scores.size()) + " quality scores from " + ffile + ", but read " + toString(reverse.scores.size()) + " quality scores from " + rfile + ", ignoring.\n");
1455 catch(exception& e) {
1456 m->errorOut(e, "MakeContigsCommand", "checkReads");
1460 //***************************************************************************************************************
1461 vector< vector<string> > MakeContigsCommand::readFileNames(string filename){
1463 vector< vector<string> > files;
1464 string forward, reverse;
1467 m->openInputFile(filename, in);
1471 if (m->control_pressed) { return files; }
1473 in >> forward; m->gobble(in);
1474 in >> reverse; m->gobble(in);
1476 //check to make sure both are able to be opened
1478 int openForward = m->openInputFile(forward, in2, "noerror");
1480 //if you can't open it, try default location
1481 if (openForward == 1) {
1482 if (m->getDefaultPath() != "") { //default path is set
1483 string tryPath = m->getDefaultPath() + m->getSimpleName(forward);
1484 m->mothurOut("Unable to open " + forward + ". Trying default " + tryPath); m->mothurOutEndLine();
1486 openForward = m->openInputFile(tryPath, in3, "noerror");
1492 //if you can't open it, try output location
1493 if (openForward == 1) {
1494 if (m->getOutputDir() != "") { //default path is set
1495 string tryPath = m->getOutputDir() + m->getSimpleName(forward);
1496 m->mothurOut("Unable to open " + forward + ". Trying output directory " + tryPath); m->mothurOutEndLine();
1498 openForward = m->openInputFile(tryPath, in4, "noerror");
1504 if (openForward == 1) { //can't find it
1505 m->mothurOut("[WARNING]: can't find " + forward + ", ignoring pair.\n");
1506 }else{ in2.close(); }
1509 int openReverse = m->openInputFile(reverse, in3, "noerror");
1511 //if you can't open it, try default location
1512 if (openReverse == 1) {
1513 if (m->getDefaultPath() != "") { //default path is set
1514 string tryPath = m->getDefaultPath() + m->getSimpleName(reverse);
1515 m->mothurOut("Unable to open " + reverse + ". Trying default " + tryPath); m->mothurOutEndLine();
1517 openReverse = m->openInputFile(tryPath, in3, "noerror");
1523 //if you can't open it, try output location
1524 if (openReverse == 1) {
1525 if (m->getOutputDir() != "") { //default path is set
1526 string tryPath = m->getOutputDir() + m->getSimpleName(reverse);
1527 m->mothurOut("Unable to open " + reverse + ". Trying output directory " + tryPath); m->mothurOutEndLine();
1529 openReverse = m->openInputFile(tryPath, in4, "noerror");
1535 if (openReverse == 1) { //can't find it
1536 m->mothurOut("[WARNING]: can't find " + reverse + ", ignoring pair.\n");
1537 }else{ in3.close(); }
1539 if ((openForward != 1) && (openReverse != 1)) { //good pair
1540 vector<string> pair;
1541 pair.push_back(forward);
1542 pair.push_back(reverse);
1543 files.push_back(pair);
1551 catch(exception& e) {
1552 m->errorOut(e, "MakeContigsCommand", "checkReads");
1556 //***************************************************************************************************************
1557 //illumina data requires paired forward and reverse data
1558 //BARCODE atgcatgc atgcatgc groupName
1559 //PRIMER atgcatgc atgcatgc groupName
1560 //PRIMER atgcatgc atgcatgc
1561 bool MakeContigsCommand::getOligos(vector<vector<string> >& fastaFileNames, string rootname){
1564 m->openInputFile(oligosfile, in);
1568 string type, foligo, roligo, group;
1570 int indexPrimer = 0;
1571 int indexBarcode = 0;
1572 set<string> uniquePrimers;
1573 set<string> uniqueBarcodes;
1579 if (m->debug) { m->mothurOut("[DEBUG]: reading type - " + type + ".\n"); }
1582 while (!in.eof()) { char c = in.get(); if (c == 10 || c == 13){ break; } } // get rest of line if there's any crap there
1587 //make type case insensitive
1588 for(int i=0;i<type.length();i++){ type[i] = toupper(type[i]); }
1592 if (m->debug) { m->mothurOut("[DEBUG]: reading - " + foligo + ".\n"); }
1594 for(int i=0;i<foligo.length();i++){
1595 foligo[i] = toupper(foligo[i]);
1596 if(foligo[i] == 'U') { foligo[i] = 'T'; }
1599 if(type == "FORWARD"){
1604 for(int i=0;i<roligo.length();i++){
1605 roligo[i] = toupper(roligo[i]);
1606 if(roligo[i] == 'U') { roligo[i] = 'T'; }
1608 //roligo = reverseOligo(roligo);
1612 // get rest of line in case there is a primer name
1615 if (c == 10 || c == 13){ break; }
1616 else if (c == 32 || c == 9){;} //space or tab
1617 else { group += c; }
1620 oligosPair newPrimer(foligo, roligo);
1622 //check for repeat barcodes
1623 string tempPair = foligo+roligo;
1624 if (uniquePrimers.count(tempPair) != 0) { m->mothurOut("primer pair " + newPrimer.forward + " " + newPrimer.reverse + " is in your oligos file already."); m->mothurOutEndLine(); }
1625 else { uniquePrimers.insert(tempPair); }
1627 if (m->debug) { if (group != "") { m->mothurOut("[DEBUG]: reading group " + group + ".\n"); }else{ m->mothurOut("[DEBUG]: no group for primer pair " + newPrimer.forward + " " + newPrimer.reverse + ".\n"); } }
1629 primers[indexPrimer]=newPrimer; indexPrimer++;
1630 primerNameVector.push_back(group);
1631 }else if(type == "BARCODE"){
1636 for(int i=0;i<roligo.length();i++){
1637 roligo[i] = toupper(roligo[i]);
1638 if(roligo[i] == 'U') { roligo[i] = 'T'; }
1640 //roligo = reverseOligo(roligo);
1642 oligosPair newPair(foligo, roligo);
1647 if (c == 10 || c == 13){ break; }
1648 else if (c == 32 || c == 9){;} //space or tab
1649 else { group += c; }
1652 if (m->debug) { m->mothurOut("[DEBUG]: barcode pair " + newPair.forward + " " + newPair.reverse + ", and group = " + group + ".\n"); }
1654 //check for repeat barcodes
1655 string tempPair = foligo+roligo;
1656 if (uniqueBarcodes.count(tempPair) != 0) { m->mothurOut("barcode pair " + newPair.forward + " " + newPair.reverse + " is in your oligos file already, disregarding."); m->mothurOutEndLine(); }
1657 else { uniqueBarcodes.insert(tempPair); }
1659 barcodes[indexBarcode]=newPair; indexBarcode++;
1660 barcodeNameVector.push_back(group);
1661 }else if(type == "LINKER"){
1662 linker.push_back(foligo);
1663 m->mothurOut("[WARNING]: make.contigs is not setup to remove linkers, ignoring.\n");
1664 }else if(type == "SPACER"){
1665 spacer.push_back(foligo);
1666 m->mothurOut("[WARNING]: make.contigs is not setup to remove spacers, ignoring.\n");
1668 else{ m->mothurOut("[WARNING]: " + type + " is not recognized as a valid type. Choices are primer, barcode, linker and spacer. Ignoring " + foligo + "."); m->mothurOutEndLine(); }
1674 if(barcodeNameVector.size() == 0 && primerNameVector[0] == ""){ allFiles = 0; }
1676 //add in potential combos
1677 if(barcodeNameVector.size() == 0){
1678 oligosPair temp("", "");
1680 barcodeNameVector.push_back("");
1683 if(primerNameVector.size() == 0){
1684 oligosPair temp("", "");
1686 primerNameVector.push_back("");
1689 fastaFileNames.resize(barcodeNameVector.size());
1690 for(int i=0;i<fastaFileNames.size();i++){
1691 fastaFileNames[i].assign(primerNameVector.size(), "");
1695 set<string> uniqueNames; //used to cleanup outputFileNames
1696 for(map<int, oligosPair>::iterator itBar = barcodes.begin();itBar != barcodes.end();itBar++){
1697 for(map<int, oligosPair>::iterator itPrimer = primers.begin();itPrimer != primers.end(); itPrimer++){
1699 string primerName = primerNameVector[itPrimer->first];
1700 string barcodeName = barcodeNameVector[itBar->first];
1702 if ((primerName == "ignore") || (barcodeName == "ignore")) { } //do nothing
1704 string comboGroupName = "";
1705 string fastaFileName = "";
1706 string qualFileName = "";
1707 string nameFileName = "";
1708 string countFileName = "";
1710 if(primerName == ""){
1711 comboGroupName = barcodeNameVector[itBar->first];
1714 if(barcodeName == ""){
1715 comboGroupName = primerNameVector[itPrimer->first];
1718 comboGroupName = barcodeNameVector[itBar->first] + "." + primerNameVector[itPrimer->first];
1724 fastaFileName = rootname + comboGroupName + ".fasta";
1725 if (uniqueNames.count(fastaFileName) == 0) {
1726 outputNames.push_back(fastaFileName);
1727 outputTypes["fasta"].push_back(fastaFileName);
1728 uniqueNames.insert(fastaFileName);
1731 fastaFileNames[itBar->first][itPrimer->first] = fastaFileName;
1732 m->openOutputFile(fastaFileName, temp); temp.close();
1738 bool allBlank = true;
1739 for (int i = 0; i < barcodeNameVector.size(); i++) {
1740 if (barcodeNameVector[i] != "") {
1745 for (int i = 0; i < primerNameVector.size(); i++) {
1746 if (primerNameVector[i] != "") {
1753 m->mothurOut("[WARNING]: your oligos file does not contain any group names. mothur will not create a groupfile."); m->mothurOutEndLine();
1761 catch(exception& e) {
1762 m->errorOut(e, "MakeContigsCommand", "getOligos");
1766 //********************************************************************/
1767 string MakeContigsCommand::reverseOligo(string oligo){
1769 string reverse = "";
1771 for(int i=oligo.length()-1;i>=0;i--){
1773 if(oligo[i] == 'A') { reverse += 'T'; }
1774 else if(oligo[i] == 'T'){ reverse += 'A'; }
1775 else if(oligo[i] == 'U'){ reverse += 'A'; }
1777 else if(oligo[i] == 'G'){ reverse += 'C'; }
1778 else if(oligo[i] == 'C'){ reverse += 'G'; }
1780 else if(oligo[i] == 'R'){ reverse += 'Y'; }
1781 else if(oligo[i] == 'Y'){ reverse += 'R'; }
1783 else if(oligo[i] == 'M'){ reverse += 'K'; }
1784 else if(oligo[i] == 'K'){ reverse += 'M'; }
1786 else if(oligo[i] == 'W'){ reverse += 'W'; }
1787 else if(oligo[i] == 'S'){ reverse += 'S'; }
1789 else if(oligo[i] == 'B'){ reverse += 'V'; }
1790 else if(oligo[i] == 'V'){ reverse += 'B'; }
1792 else if(oligo[i] == 'D'){ reverse += 'H'; }
1793 else if(oligo[i] == 'H'){ reverse += 'D'; }
1795 else { reverse += 'N'; }
1801 catch(exception& e) {
1802 m->errorOut(e, "MakeContigsCommand", "reverseOligo");
1806 //**********************************************************************************************************************
1807 vector<int> MakeContigsCommand::convertQual(string qual) {
1809 vector<int> qualScores;
1810 bool negativeScores = false;
1812 for (int i = 0; i < qual.length(); i++) {
1815 temp = int(qual[i]);
1816 if (format == "illumina") {
1817 temp -= 64; //char '@'
1818 }else if (format == "illumina1.8+") {
1819 temp -= int('!'); //char '!'
1820 }else if (format == "solexa") {
1821 temp = int(convertTable[temp]); //convert to sanger
1822 temp -= int('!'); //char '!'
1824 temp -= int('!'); //char '!'
1827 if (temp < -5) { negativeScores = true; }
1828 qualScores.push_back(temp);
1831 if (negativeScores) { m->mothurOut("[ERROR]: finding negative quality scores, do you have the right format selected? http://en.wikipedia.org/wiki/FASTQ_format#Encoding \n"); m->control_pressed = true; }
1835 catch(exception& e) {
1836 m->errorOut(e, "MakeContigsCommand", "convertQual");
1841 //**********************************************************************************************************************