2 * pairwiseseqscommand.cpp
5 * Created by westcott on 10/20/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "pairwiseseqscommand.h"
12 //**********************************************************************************************************************
13 vector<string> PairwiseSeqsCommand::setParameters(){
15 CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "none", "none","phylip-column",false,true,true); parameters.push_back(pfasta);
16 CommandParameter palign("align", "Multiple", "needleman-gotoh-blast-noalign", "needleman", "", "", "","",false,false); parameters.push_back(palign);
17 CommandParameter pmatch("match", "Number", "", "1.0", "", "", "","",false,false); parameters.push_back(pmatch);
18 CommandParameter pmismatch("mismatch", "Number", "", "-1.0", "", "", "","",false,false); parameters.push_back(pmismatch);
19 CommandParameter pgapopen("gapopen", "Number", "", "-2.0", "", "", "","",false,false); parameters.push_back(pgapopen);
20 CommandParameter pgapextend("gapextend", "Number", "", "-1.0", "", "", "","",false,false); parameters.push_back(pgapextend);
21 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors);
22 CommandParameter poutput("output", "Multiple", "column-lt-square-phylip", "column", "", "", "","phylip-column",false,false,true); parameters.push_back(poutput);
23 CommandParameter pcalc("calc", "Multiple", "nogaps-eachgap-onegap", "onegap", "", "", "","",false,false); parameters.push_back(pcalc);
24 CommandParameter pcountends("countends", "Boolean", "", "T", "", "", "","",false,false); parameters.push_back(pcountends);
25 CommandParameter pcompress("compress", "Boolean", "", "F", "", "", "","",false,false); parameters.push_back(pcompress);
26 CommandParameter pcutoff("cutoff", "Number", "", "1.0", "", "", "","",false,false,true); parameters.push_back(pcutoff);
27 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
28 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
30 vector<string> myArray;
31 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
35 m->errorOut(e, "PairwiseSeqsCommand", "setParameters");
39 //**********************************************************************************************************************
40 string PairwiseSeqsCommand::getHelpString(){
42 string helpString = "";
43 helpString += "The pairwise.seqs command reads a fasta file and creates distance matrix.\n";
44 helpString += "The pairwise.seqs command parameters are fasta, align, match, mismatch, gapopen, gapextend, calc, output, cutoff and processors.\n";
45 helpString += "The fasta parameter is required. You may enter multiple fasta files by separating their names with dashes. ie. fasta=abrecovery.fasta-amzon.fasta \n";
46 helpString += "The align parameter allows you to specify the alignment method to use. Your options are: gotoh, needleman, blast and noalign. The default is needleman.\n";
47 helpString += "The match parameter allows you to specify the bonus for having the same base. The default is 1.0.\n";
48 helpString += "The mistmatch parameter allows you to specify the penalty for having different bases. The default is -1.0.\n";
49 helpString += "The gapopen parameter allows you to specify the penalty for opening a gap in an alignment. The default is -2.0.\n";
50 helpString += "The gapextend parameter allows you to specify the penalty for extending a gap in an alignment. The default is -1.0.\n";
51 helpString += "The calc parameter allows you to specify the method of calculating the distances. Your options are: nogaps, onegap or eachgap. The default is onegap.\n";
52 helpString += "The countends parameter allows you to specify whether to include terminal gaps in distance. Your options are: T or F. The default is T.\n";
53 helpString += "The cutoff parameter allows you to specify maximum distance to keep. The default is 1.0.\n";
54 helpString += "The output parameter allows you to specify format of your distance matrix. Options are column, lt, and square. The default is column.\n";
55 helpString += "The compress parameter allows you to indicate that you want the resulting distance file compressed. The default is false.\n";
56 helpString += "The pairwise.seqs command should be in the following format: \n";
57 helpString += "pairwise.seqs(fasta=yourfastaFile, align=yourAlignmentMethod) \n";
58 helpString += "Example pairwise.seqs(fasta=candidate.fasta, align=blast)\n";
59 helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFastaFile).\n";
63 m->errorOut(e, "PairwiseSeqsCommand", "getHelpString");
67 //**********************************************************************************************************************
68 string PairwiseSeqsCommand::getOutputPattern(string type) {
72 if (type == "phylip") { pattern = "[filename],[outputtag],dist"; }
73 else if (type == "column") { pattern = "[filename],dist"; }
74 else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; }
79 m->errorOut(e, "PairwiseSeqsCommand", "getOutputPattern");
83 //**********************************************************************************************************************
84 PairwiseSeqsCommand::PairwiseSeqsCommand(){
86 abort = true; calledHelp = true;
88 vector<string> tempOutNames;
89 outputTypes["phylip"] = tempOutNames;
90 outputTypes["column"] = tempOutNames;
93 m->errorOut(e, "PairwiseSeqsCommand", "PairwiseSeqsCommand");
97 //**********************************************************************************************************************
98 PairwiseSeqsCommand::PairwiseSeqsCommand(string option) {
100 abort = false; calledHelp = false;
102 //allow user to run help
103 if(option == "help") { help(); abort = true; calledHelp = true; }
104 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
107 vector<string> myArray = setParameters();
109 OptionParser parser(option);
110 map<string, string> parameters = parser.getParameters();
112 ValidParameters validParameter("pairwise.seqs");
113 map<string, string>::iterator it;
115 //check to make sure all parameters are valid for command
116 for (it = parameters.begin(); it != parameters.end(); it++) {
117 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
120 //initialize outputTypes
121 vector<string> tempOutNames;
122 outputTypes["phylip"] = tempOutNames;
123 outputTypes["column"] = tempOutNames;
125 //if the user changes the output directory command factory will send this info to us in the output parameter
126 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; }
129 //if the user changes the input directory command factory will send this info to us in the output parameter
130 string inputDir = validParameter.validFile(parameters, "inputdir", false);
132 if (inputDir == "not found"){ inputDir = ""; }
134 fastaFileName = validParameter.validFile(parameters, "fasta", false);
135 if (fastaFileName == "not found") {
136 //if there is a current fasta file, use it
137 string filename = m->getFastaFile();
138 if (filename != "") { fastaFileNames.push_back(filename); m->mothurOut("Using " + filename + " as input file for the fasta parameter."); m->mothurOutEndLine(); }
139 else { m->mothurOut("You have no current fastafile and the fasta parameter is required."); m->mothurOutEndLine(); abort = true; }
142 m->splitAtDash(fastaFileName, fastaFileNames);
144 //go through files and make sure they are good, if not, then disregard them
145 for (int i = 0; i < fastaFileNames.size(); i++) {
148 if (fastaFileNames[i] == "current") {
149 fastaFileNames[i] = m->getFastaFile();
150 if (fastaFileNames[i] != "") { m->mothurOut("Using " + fastaFileNames[i] + " as input file for the fasta parameter where you had given current."); m->mothurOutEndLine(); }
152 m->mothurOut("You have no current fastafile, ignoring current."); m->mothurOutEndLine(); ignore=true;
153 //erase from file list
154 fastaFileNames.erase(fastaFileNames.begin()+i);
161 if (inputDir != "") {
162 string path = m->hasPath(fastaFileNames[i]);
163 //if the user has not given a path then, add inputdir. else leave path alone.
164 if (path == "") { fastaFileNames[i] = inputDir + fastaFileNames[i]; }
170 ableToOpen = m->openInputFile(fastaFileNames[i], in, "noerror");
172 //if you can't open it, try default location
173 if (ableToOpen == 1) {
174 if (m->getDefaultPath() != "") { //default path is set
175 string tryPath = m->getDefaultPath() + m->getSimpleName(fastaFileNames[i]);
176 m->mothurOut("Unable to open " + fastaFileNames[i] + ". Trying default " + tryPath); m->mothurOutEndLine();
178 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
180 fastaFileNames[i] = tryPath;
184 //if you can't open it, try output location
185 if (ableToOpen == 1) {
186 if (m->getOutputDir() != "") { //default path is set
187 string tryPath = m->getOutputDir() + m->getSimpleName(fastaFileNames[i]);
188 m->mothurOut("Unable to open " + fastaFileNames[i] + ". Trying output directory " + tryPath); m->mothurOutEndLine();
190 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
192 fastaFileNames[i] = tryPath;
198 if (ableToOpen == 1) {
199 m->mothurOut("Unable to open " + fastaFileNames[i] + ". It will be disregarded."); m->mothurOutEndLine();
200 //erase from file list
201 fastaFileNames.erase(fastaFileNames.begin()+i);
204 m->setFastaFile(fastaFileNames[i]);
209 //make sure there is at least one valid file left
210 if (fastaFileNames.size() == 0) { m->mothurOut("no valid files."); m->mothurOutEndLine(); abort = true; }
213 //check for optional parameter and set defaults
214 // ...at some point should added some additional type checking...
216 temp = validParameter.validFile(parameters, "match", false); if (temp == "not found"){ temp = "1.0"; }
217 m->mothurConvert(temp, match);
219 temp = validParameter.validFile(parameters, "mismatch", false); if (temp == "not found"){ temp = "-1.0"; }
220 m->mothurConvert(temp, misMatch);
221 if (misMatch > 0) { m->mothurOut("[ERROR]: mismatch must be negative.\n"); abort=true; }
223 temp = validParameter.validFile(parameters, "gapopen", false); if (temp == "not found"){ temp = "-2.0"; }
224 m->mothurConvert(temp, gapOpen);
225 if (gapOpen > 0) { m->mothurOut("[ERROR]: gapopen must be negative.\n"); abort=true; }
227 temp = validParameter.validFile(parameters, "gapextend", false); if (temp == "not found"){ temp = "-1.0"; }
228 m->mothurConvert(temp, gapExtend);
229 if (gapExtend > 0) { m->mothurOut("[ERROR]: gapextend must be negative.\n"); abort=true; }
231 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
232 m->setProcessors(temp);
233 m->mothurConvert(temp, processors);
235 temp = validParameter.validFile(parameters, "cutoff", false); if(temp == "not found"){ temp = "1.0"; }
236 m->mothurConvert(temp, cutoff);
238 temp = validParameter.validFile(parameters, "countends", false); if(temp == "not found"){ temp = "T"; }
239 countends = m->isTrue(temp);
241 temp = validParameter.validFile(parameters, "compress", false); if(temp == "not found"){ temp = "F"; }
242 compress = m->isTrue(temp);
244 align = validParameter.validFile(parameters, "align", false); if (align == "not found"){ align = "needleman"; }
246 output = validParameter.validFile(parameters, "output", false); if(output == "not found"){ output = "column"; }
247 if (output=="phylip") { output = "lt"; }
248 if ((output != "column") && (output != "lt") && (output != "square")) { m->mothurOut(output + " is not a valid output form. Options are column, lt and square. I will use column."); m->mothurOutEndLine(); output = "column"; }
250 calc = validParameter.validFile(parameters, "calc", false);
251 if (calc == "not found") { calc = "onegap"; }
253 if (calc == "default") { calc = "onegap"; }
255 m->splitAtDash(calc, Estimators);
259 catch(exception& e) {
260 m->errorOut(e, "PairwiseSeqsCommand", "PairwiseSeqsCommand");
264 //**********************************************************************************************************************
266 int PairwiseSeqsCommand::execute(){
268 if (abort == true) { if (calledHelp) { return 0; } return 2; }
270 longestBase = 2000; //will need to update this in driver if we find sequences with more bases. hardcoded so we don't have the pre-read user fasta file.
274 for (int s = 0; s < fastaFileNames.size(); s++) {
275 if (m->control_pressed) { outputTypes.clear(); return 0; }
277 m->mothurOut("Processing sequences from " + fastaFileNames[s] + " ..." ); m->mothurOutEndLine();
279 if (outputDir == "") { outputDir += m->hasPath(fastaFileNames[s]); }
282 m->openInputFile(fastaFileNames[s], inFASTA);
283 alignDB = SequenceDB(inFASTA);
286 int numSeqs = alignDB.getNumSeqs();
287 int startTime = time(NULL);
288 string outputFile = "";
290 map<string, string> variables;
291 variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[s]));
292 if (output == "lt") { //does the user want lower triangle phylip formatted file
293 variables["[outputtag]"] = "phylip";
294 outputFile = getOutputFileName("phylip", variables);
295 m->mothurRemove(outputFile); outputTypes["phylip"].push_back(outputFile);
296 }else if (output == "column") { //user wants column format
297 outputFile = getOutputFileName("column", variables);
298 outputTypes["column"].push_back(outputFile);
299 m->mothurRemove(outputFile);
300 }else { //assume square
301 variables["[outputtag]"] = "square";
302 outputFile = getOutputFileName("phylip", variables);
303 m->mothurRemove(outputFile);
304 outputTypes["phylip"].push_back(outputFile);
313 MPI_Comm_size(MPI_COMM_WORLD, &processors); //set processors to the number of mpi processes running
314 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
316 //each process gets where it should start and stop in the file
317 if (output != "square") {
318 start = int (sqrt(float(pid)/float(processors)) * numSeqs);
319 end = int (sqrt(float(pid+1)/float(processors)) * numSeqs);
321 start = int ((float(pid)/float(processors)) * numSeqs);
322 end = int ((float(pid+1)/float(processors)) * numSeqs);
325 if (output == "column") {
327 int amode=MPI_MODE_CREATE|MPI_MODE_WRONLY;
330 strcpy(filename, outputFile.c_str());
332 MPI_File_open(MPI_COMM_WORLD, filename, amode, MPI_INFO_NULL, &outMPI);
334 if (pid == 0) { //you are the root process
339 driverMPI(start, end, outMPI, cutoff);
341 if (m->control_pressed) { outputTypes.clear(); MPI_File_close(&outMPI); m->mothurRemove(outputFile); return 0; }
344 for(int i = 1; i < processors; i++) {
345 if (m->control_pressed) { outputTypes.clear(); MPI_File_close(&outMPI); m->mothurRemove(outputFile); return 0; }
348 MPI_Recv(buf, 5, MPI_CHAR, i, tag, MPI_COMM_WORLD, &status);
350 }else { //you are a child process
352 driverMPI(start, end, outMPI, cutoff);
354 if (m->control_pressed) { outputTypes.clear(); MPI_File_close(&outMPI); m->mothurRemove(outputFile); return 0; }
358 //tell parent you are done.
359 MPI_Send(buf, 5, MPI_CHAR, 0, tag, MPI_COMM_WORLD);
362 MPI_File_close(&outMPI);
364 }else { //lower triangle format
365 if (pid == 0) { //you are the root process
369 unsigned long long mySize;
371 if (output != "square"){ driverMPI(start, end, outputFile, mySize); }
372 else { driverMPI(start, end, outputFile, mySize, output); }
374 if (m->control_pressed) { outputTypes.clear(); m->mothurRemove(outputFile); return 0; }
376 int amode=MPI_MODE_APPEND|MPI_MODE_WRONLY|MPI_MODE_CREATE; //
381 strcpy(filename, outputFile.c_str());
383 MPI_File_open(MPI_COMM_SELF, filename, amode, MPI_INFO_NULL, &outMPI);
386 for(int b = 1; b < processors; b++) {
387 unsigned long long fileSize;
389 if (m->control_pressed) { outputTypes.clear(); MPI_File_close(&outMPI); m->mothurRemove(outputFile); return 0; }
391 MPI_Recv(&fileSize, 1, MPI_LONG, b, tag, MPI_COMM_WORLD, &status);
393 string outTemp = outputFile + toString(b) + ".temp";
395 char* buf = new char[outTemp.length()];
396 memcpy(buf, outTemp.c_str(), outTemp.length());
398 MPI_File_open(MPI_COMM_SELF, buf, MPI_MODE_DELETE_ON_CLOSE|MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);
402 while (count < fileSize) {
404 MPI_File_read(inMPI, buf2, 1, MPI_CHAR, &status);
405 MPI_File_write(outMPI, buf2, 1, MPI_CHAR, &status);
409 MPI_File_close(&inMPI); //deleted on close
412 MPI_File_close(&outMPI);
413 }else { //you are a child process
415 unsigned long long size;
416 if (output != "square"){ driverMPI(start, end, (outputFile + toString(pid) + ".temp"), size); }
417 else { driverMPI(start, end, (outputFile + toString(pid) + ".temp"), size, output); }
419 if (m->control_pressed) { return 0; }
421 //tell parent you are done.
422 MPI_Send(&size, 1, MPI_LONG, 0, tag, MPI_COMM_WORLD);
425 MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
428 //#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
429 //if you don't need to fork anything
431 if (output != "square") { driver(0, numSeqs, outputFile, cutoff); }
432 else { driver(0, numSeqs, outputFile, "square"); }
433 }else{ //you have multiple processors
435 for (int i = 0; i < processors; i++) {
436 distlinePair tempLine;
437 lines.push_back(tempLine);
438 if (output != "square") {
439 lines[i].start = int (sqrt(float(i)/float(processors)) * numSeqs);
440 lines[i].end = int (sqrt(float(i+1)/float(processors)) * numSeqs);
442 lines[i].start = int ((float(i)/float(processors)) * numSeqs);
443 lines[i].end = int ((float(i+1)/float(processors)) * numSeqs);
447 createProcesses(outputFile);
451 //if (output != "square") { driver(0, numSeqs, outputFile, cutoff); }
452 //else { driver(0, numSeqs, outputFile, "square"); }
456 if (m->control_pressed) { outputTypes.clear(); m->mothurRemove(outputFile); return 0; }
459 MPI_Comm_rank(MPI_COMM_WORLD, &pid);
461 if (pid == 0) { //only one process should output to screen
465 fileHandle.open(outputFile.c_str());
467 m->gobble(fileHandle);
468 if (fileHandle.eof()) { m->mothurOut(outputFile + " is blank. This can result if there are no distances below your cutoff."); m->mothurOutEndLine(); }
472 m->mothurOut("Compressing..."); m->mothurOutEndLine();
473 m->mothurOut("(Replacing " + outputFile + " with " + outputFile + ".gz)"); m->mothurOutEndLine();
474 system(("gzip -v " + outputFile).c_str());
475 outputNames.push_back(outputFile + ".gz");
476 }else { outputNames.push_back(outputFile); }
482 m->mothurOut("It took " + toString(time(NULL) - startTime) + " to calculate the distances for " + toString(numSeqs) + " sequences."); m->mothurOutEndLine();
484 if (m->control_pressed) { outputTypes.clear(); m->mothurRemove(outputFile); return 0; }
487 //set phylip file as new current phylipfile
489 itTypes = outputTypes.find("phylip");
490 if (itTypes != outputTypes.end()) {
491 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setPhylipFile(current); }
494 //set column file as new current columnfile
495 itTypes = outputTypes.find("column");
496 if (itTypes != outputTypes.end()) {
497 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setColumnFile(current); }
500 m->mothurOutEndLine();
501 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
502 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
503 m->mothurOutEndLine();
508 catch(exception& e) {
509 m->errorOut(e, "PairwiseSeqsCommand", "execute");
514 /**************************************************************************************************/
515 void PairwiseSeqsCommand::createProcesses(string filename) {
520 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
523 //loop through and create all the processes you want
524 while (process != processors) {
528 processIDS.push_back(pid);
531 if (output != "square") { driver(lines[process].start, lines[process].end, filename + toString(getpid()) + ".temp", cutoff); }
532 else { driver(lines[process].start, lines[process].end, filename + toString(getpid()) + ".temp", "square"); }
535 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
536 for (int i=0;i<processIDS.size();i++) { int temp = processIDS[i]; kill (temp, SIGINT); }
542 if (output != "square") { driver(lines[0].start, lines[0].end, filename, cutoff); }
543 else { driver(lines[0].start, lines[0].end, filename, "square"); }
546 //force parent to wait until all the processes are done
547 for (int i=0;i<processIDS.size();i++) {
548 int temp = processIDS[i];
552 //////////////////////////////////////////////////////////////////////////////////////////////////////
553 //Windows version shared memory, so be careful when passing variables through the distanceData struct.
554 //Above fork() will clone, so memory is separate, but that's not the case with windows,
555 //that's why the distance calculator was moved inside of the driver to make separate copies.
556 //////////////////////////////////////////////////////////////////////////////////////////////////////
558 vector<pairwiseData*> pDataArray; //[processors-1];
559 DWORD dwThreadIdArray[processors-1];
560 HANDLE hThreadArray[processors-1];
562 //Create processor-1 worker threads.
563 for( int i=0; i<processors-1; i++ ){
564 string extension = toString(i) + ".temp";
566 // Allocate memory for thread data.
567 pairwiseData* tempDist = new pairwiseData((filename+extension), align, "square", Estimators[0], countends, output, alignDB, m, lines[i+1].start, lines[i+1].end, match, misMatch, gapOpen, gapExtend, longestBase, i);
568 pDataArray.push_back(tempDist);
569 processIDS.push_back(i);
571 if (output != "square") { hThreadArray[i] = CreateThread(NULL, 0, MyPairwiseThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]); }
572 else { hThreadArray[i] = CreateThread(NULL, 0, MyPairwiseSquareThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]); }
576 if (output != "square") { driver(lines[0].start, lines[0].end, filename, cutoff); }
577 else { driver(lines[0].start, lines[0].end, filename, "square"); }
579 //Wait until all threads have terminated.
580 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
582 //Close all thread handles and free memory allocations.
583 for(int i=0; i < pDataArray.size(); i++){
584 if (pDataArray[i]->count != (pDataArray[i]->end-pDataArray[i]->start)) {
585 m->mothurOut("[ERROR]: process " + toString(i) + " only processed " + toString(pDataArray[i]->count) + " of " + toString(pDataArray[i]->end-pDataArray[i]->start) + " sequences assigned to it, quitting. \n"); m->control_pressed = true;
587 CloseHandle(hThreadArray[i]);
588 delete pDataArray[i];
593 //append and remove temp files
594 for (int i=0;i<processIDS.size();i++) {
595 m->appendFiles((filename + toString(processIDS[i]) + ".temp"), filename);
596 m->mothurRemove((filename + toString(processIDS[i]) + ".temp"));
600 catch(exception& e) {
601 m->errorOut(e, "PairwiseSeqsCommand", "createProcesses");
606 /**************************************************************************************************/
607 /////// need to fix to work with calcs and sequencedb
608 int PairwiseSeqsCommand::driver(int startLine, int endLine, string dFileName, float cutoff){
611 int startTime = time(NULL);
613 Alignment* alignment;
614 if(align == "gotoh") { alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase); }
615 else if(align == "needleman") { alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase); }
616 else if(align == "blast") { alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch); }
617 else if(align == "noalign") { alignment = new NoAlign(); }
619 m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
620 m->mothurOutEndLine();
621 alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
624 ValidCalculators validCalculator;
625 Dist* distCalculator;
627 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
628 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
629 else if (Estimators[0] == "eachgap") { distCalculator = new eachGapDist(); }
630 else if (Estimators[0] == "onegap") { distCalculator = new oneGapDist(); }
633 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
634 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
635 else if (Estimators[0] == "eachgap"){ distCalculator = new eachGapIgnoreTermGapDist(); }
636 else if (Estimators[0] == "onegap") { distCalculator = new oneGapIgnoreTermGapDist(); }
641 ofstream outFile(dFileName.c_str(), ios::trunc);
642 outFile.setf(ios::fixed, ios::showpoint);
643 outFile << setprecision(4);
645 if((output == "lt") && startLine == 0){ outFile << alignDB.getNumSeqs() << endl; }
648 for(int i=startLine;i<endLine;i++){
650 string name = alignDB.get(i).getName();
651 if (name.length() < 10) { //pad with spaces to make compatible
652 while (name.length() < 10) { name += " "; }
654 outFile << name << '\t';
657 for(int j=0;j<i;j++){
659 if (m->control_pressed) { outFile.close(); delete alignment; delete distCalculator; return 0; }
661 if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
662 alignment->resize(alignDB.get(i).getUnaligned().length()+1);
665 if (alignDB.get(j).getUnaligned().length() > alignment->getnRows()) {
666 alignment->resize(alignDB.get(j).getUnaligned().length()+1);
669 Sequence seqI(alignDB.get(i).getName(), alignDB.get(i).getAligned());
670 Sequence seqJ(alignDB.get(j).getName(), alignDB.get(j).getAligned());
672 alignment->align(seqI.getUnaligned(), seqJ.getUnaligned());
673 seqI.setAligned(alignment->getSeqAAln());
674 seqJ.setAligned(alignment->getSeqBAln());
676 //cout << seqI.getName() << '\t' << seqJ.getName() << endl;
677 //cout << alignment->getSeqAAln() << endl << alignment->getSeqBAln() << endl;
679 distCalculator->calcDist(seqI, seqJ);
680 double dist = distCalculator->getDist();
682 //cout << "dist = " << dist << endl;
685 if (dist < 0.01) { countSmall++; }
688 if (output == "column") { outFile << alignDB.get(i).getName() << ' ' << alignDB.get(j).getName() << ' ' << dist << endl; }
690 if (output == "lt") { outFile << dist << '\t'; }
693 if (output == "lt") { outFile << endl; }
696 m->mothurOut(toString(i) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
700 m->mothurOut(toString(endLine-1) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
701 cout << "num less than 0.01 = " << countSmall << " of " << countAll << endl;
704 delete distCalculator;
708 catch(exception& e) {
709 m->errorOut(e, "PairwiseSeqsCommand", "driver");
713 /**************************************************************************************************/
714 /////// need to fix to work with calcs and sequencedb
715 int PairwiseSeqsCommand::driver(int startLine, int endLine, string dFileName, string square){
718 int startTime = time(NULL);
720 Alignment* alignment;
721 if(align == "gotoh") { alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase); }
722 else if(align == "needleman") { alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase); }
723 else if(align == "blast") { alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch); }
724 else if(align == "noalign") { alignment = new NoAlign(); }
726 m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
727 m->mothurOutEndLine();
728 alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
731 ValidCalculators validCalculator;
732 Dist* distCalculator;
734 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
735 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
736 else if (Estimators[0] == "eachgap") { distCalculator = new eachGapDist(); }
737 else if (Estimators[0] == "onegap") { distCalculator = new oneGapDist(); }
740 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
741 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
742 else if (Estimators[0] == "eachgap"){ distCalculator = new eachGapIgnoreTermGapDist(); }
743 else if (Estimators[0] == "onegap") { distCalculator = new oneGapIgnoreTermGapDist(); }
748 ofstream outFile(dFileName.c_str(), ios::trunc);
749 outFile.setf(ios::fixed, ios::showpoint);
750 outFile << setprecision(4);
752 if(startLine == 0){ outFile << alignDB.getNumSeqs() << endl; }
754 for(int i=startLine;i<endLine;i++){
756 string name = alignDB.get(i).getName();
757 //pad with spaces to make compatible
758 if (name.length() < 10) { while (name.length() < 10) { name += " "; } }
760 outFile << name << '\t';
762 for(int j=0;j<alignDB.getNumSeqs();j++){
764 if (m->control_pressed) { outFile.close(); delete alignment; delete distCalculator; return 0; }
766 if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
767 alignment->resize(alignDB.get(i).getUnaligned().length()+1);
770 if (alignDB.get(j).getUnaligned().length() > alignment->getnRows()) {
771 alignment->resize(alignDB.get(j).getUnaligned().length()+1);
774 Sequence seqI(alignDB.get(i).getName(), alignDB.get(i).getAligned());
775 Sequence seqJ(alignDB.get(j).getName(), alignDB.get(j).getAligned());
777 alignment->align(seqI.getUnaligned(), seqJ.getUnaligned());
778 seqI.setAligned(alignment->getSeqAAln());
779 seqJ.setAligned(alignment->getSeqBAln());
781 distCalculator->calcDist(seqI, seqJ);
782 double dist = distCalculator->getDist();
784 outFile << dist << '\t';
790 m->mothurOut(toString(i) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
794 m->mothurOut(toString(endLine-1) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
798 delete distCalculator;
802 catch(exception& e) {
803 m->errorOut(e, "PairwiseSeqsCommand", "driver");
808 /**************************************************************************************************/
809 /////// need to fix to work with calcs and sequencedb
810 int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, MPI_File& outMPI, float cutoff){
813 int startTime = time(NULL);
815 Alignment* alignment;
816 if(align == "gotoh") { alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase); }
817 else if(align == "needleman") { alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase); }
818 else if(align == "blast") { alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch); }
819 else if(align == "noalign") { alignment = new NoAlign(); }
821 m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
822 m->mothurOutEndLine();
823 alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
826 ValidCalculators validCalculator;
827 Dist* distCalculator;
829 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
830 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
831 else if (Estimators[0] == "eachgap") { distCalculator = new eachGapDist(); }
832 else if (Estimators[0] == "onegap") { distCalculator = new oneGapDist(); }
835 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
836 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
837 else if (Estimators[0] == "eachgap"){ distCalculator = new eachGapIgnoreTermGapDist(); }
838 else if (Estimators[0] == "onegap") { distCalculator = new oneGapIgnoreTermGapDist(); }
842 string outputString = "";
844 for(int i=startLine;i<endLine;i++){
846 for(int j=0;j<i;j++){
848 if (m->control_pressed) { delete alignment; delete distCalculator; return 0; }
850 if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
851 alignment->resize(alignDB.get(i).getUnaligned().length()+1);
854 if (alignDB.get(j).getUnaligned().length() > alignment->getnRows()) {
855 alignment->resize(alignDB.get(j).getUnaligned().length()+1);
858 Sequence seqI(alignDB.get(i).getName(), alignDB.get(i).getAligned());
859 Sequence seqJ(alignDB.get(j).getName(), alignDB.get(j).getAligned());
861 alignment->align(seqI.getUnaligned(), seqJ.getUnaligned());
862 seqI.setAligned(alignment->getSeqAAln());
863 seqJ.setAligned(alignment->getSeqBAln());
865 distCalculator->calcDist(seqI, seqJ);
866 double dist = distCalculator->getDist();
869 outputString += (alignDB.get(i).getName() + ' ' + alignDB.get(j).getName() + ' ' + toString(dist) + '\n');
874 //m->mothurOut(toString(i) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
875 cout << i << '\t' << (time(NULL) - startTime) << endl;
879 //send results to parent
880 int length = outputString.length();
882 char* buf = new char[length];
883 memcpy(buf, outputString.c_str(), length);
885 MPI_File_write_shared(outMPI, buf, length, MPI_CHAR, &status);
891 delete distCalculator;
894 catch(exception& e) {
895 m->errorOut(e, "PairwiseSeqsCommand", "driverMPI");
899 /**************************************************************************************************/
900 /////// need to fix to work with calcs and sequencedb
901 int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, string file, unsigned long long& size){
906 int amode=MPI_MODE_CREATE|MPI_MODE_WRONLY;
909 strcpy(filename, file.c_str());
911 MPI_File_open(MPI_COMM_SELF, filename, amode, MPI_INFO_NULL, &outMPI);
913 Alignment* alignment;
914 if(align == "gotoh") { alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase); }
915 else if(align == "needleman") { alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase); }
916 else if(align == "blast") { alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch); }
917 else if(align == "noalign") { alignment = new NoAlign(); }
919 m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
920 m->mothurOutEndLine();
921 alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
924 ValidCalculators validCalculator;
925 Dist* distCalculator;
927 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
928 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
929 else if (Estimators[0] == "eachgap") { distCalculator = new eachGapDist(); }
930 else if (Estimators[0] == "onegap") { distCalculator = new oneGapDist(); }
933 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
934 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
935 else if (Estimators[0] == "eachgap"){ distCalculator = new eachGapIgnoreTermGapDist(); }
936 else if (Estimators[0] == "onegap") { distCalculator = new oneGapIgnoreTermGapDist(); }
941 string outputString = "";
944 if(startLine == 0){ outputString += toString(alignDB.getNumSeqs()) + "\n"; }
946 for(int i=startLine;i<endLine;i++){
948 string name = alignDB.get(i).getName();
949 if (name.length() < 10) { //pad with spaces to make compatible
950 while (name.length() < 10) { name += " "; }
952 outputString += name + "\t";
954 for(int j=0;j<i;j++){
956 if (m->control_pressed) { delete alignment; delete distCalculator; return 0; }
958 if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
959 alignment->resize(alignDB.get(i).getUnaligned().length()+1);
962 if (alignDB.get(j).getUnaligned().length() > alignment->getnRows()) {
963 alignment->resize(alignDB.get(j).getUnaligned().length()+1);
966 Sequence seqI(alignDB.get(i).getName(), alignDB.get(i).getAligned());
967 Sequence seqJ(alignDB.get(j).getName(), alignDB.get(j).getAligned());
969 alignment->align(seqI.getUnaligned(), seqJ.getUnaligned());
970 seqI.setAligned(alignment->getSeqAAln());
971 seqJ.setAligned(alignment->getSeqBAln());
973 distCalculator->calcDist(seqI, seqJ);
974 double dist = distCalculator->getDist();
976 outputString += toString(dist) + "\t";
979 outputString += "\n";
981 //send results to parent
982 int length = outputString.length();
983 char* buf = new char[length];
984 memcpy(buf, outputString.c_str(), length);
986 MPI_File_write(outMPI, buf, length, MPI_CHAR, &status);
987 size += outputString.length();
992 MPI_File_close(&outMPI);
994 delete distCalculator;
998 catch(exception& e) {
999 m->errorOut(e, "PairwiseSeqsCommand", "driverMPI");
1003 /**************************************************************************************************/
1004 /////// need to fix to work with calcs and sequencedb
1005 int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, string file, unsigned long long& size, string square){
1010 int amode=MPI_MODE_CREATE|MPI_MODE_WRONLY;
1012 char filename[1024];
1013 strcpy(filename, file.c_str());
1015 MPI_File_open(MPI_COMM_SELF, filename, amode, MPI_INFO_NULL, &outMPI);
1017 Alignment* alignment;
1018 if(align == "gotoh") { alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase); }
1019 else if(align == "needleman") { alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase); }
1020 else if(align == "blast") { alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch); }
1021 else if(align == "noalign") { alignment = new NoAlign(); }
1023 m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
1024 m->mothurOutEndLine();
1025 alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
1028 ValidCalculators validCalculator;
1029 Dist* distCalculator;
1031 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
1032 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
1033 else if (Estimators[0] == "eachgap") { distCalculator = new eachGapDist(); }
1034 else if (Estimators[0] == "onegap") { distCalculator = new oneGapDist(); }
1037 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
1038 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
1039 else if (Estimators[0] == "eachgap"){ distCalculator = new eachGapIgnoreTermGapDist(); }
1040 else if (Estimators[0] == "onegap") { distCalculator = new oneGapIgnoreTermGapDist(); }
1044 string outputString = "";
1047 if(startLine == 0){ outputString += toString(alignDB.getNumSeqs()) + "\n"; }
1049 for(int i=startLine;i<endLine;i++){
1051 string name = alignDB.get(i).getName();
1052 if (name.length() < 10) { //pad with spaces to make compatible
1053 while (name.length() < 10) { name += " "; }
1055 outputString += name + "\t";
1057 for(int j=0;j<alignDB.getNumSeqs();j++){
1059 if (m->control_pressed) { delete alignment; return 0; }
1061 if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
1062 alignment->resize(alignDB.get(i).getUnaligned().length()+1);
1065 if (alignDB.get(j).getUnaligned().length() > alignment->getnRows()) {
1066 alignment->resize(alignDB.get(j).getUnaligned().length()+1);
1069 Sequence seqI(alignDB.get(i).getName(), alignDB.get(i).getAligned());
1070 Sequence seqJ(alignDB.get(j).getName(), alignDB.get(j).getAligned());
1072 alignment->align(seqI.getUnaligned(), seqJ.getUnaligned());
1073 seqI.setAligned(alignment->getSeqAAln());
1074 seqJ.setAligned(alignment->getSeqBAln());
1076 distCalculator->calcDist(seqI, seqJ);
1077 double dist = distCalculator->getDist();
1079 outputString += toString(dist) + "\t";
1082 outputString += "\n";
1084 //send results to parent
1085 int length = outputString.length();
1086 char* buf = new char[length];
1087 memcpy(buf, outputString.c_str(), length);
1089 MPI_File_write(outMPI, buf, length, MPI_CHAR, &status);
1090 size += outputString.length();
1095 MPI_File_close(&outMPI);
1100 catch(exception& e) {
1101 m->errorOut(e, "PairwiseSeqsCommand", "driverMPI");
1106 /**************************************************************************************************/