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, cutoff, 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; }
647 for(int i=startLine;i<endLine;i++){
649 string name = alignDB.get(i).getName();
650 if (name.length() < 10) { //pad with spaces to make compatible
651 while (name.length() < 10) { name += " "; }
653 outFile << name << '\t';
656 for(int j=0;j<i;j++){
658 if (m->control_pressed) { outFile.close(); delete alignment; delete distCalculator; return 0; }
660 if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
661 alignment->resize(alignDB.get(i).getUnaligned().length()+1);
664 if (alignDB.get(j).getUnaligned().length() > alignment->getnRows()) {
665 alignment->resize(alignDB.get(j).getUnaligned().length()+1);
668 Sequence seqI(alignDB.get(i).getName(), alignDB.get(i).getAligned());
669 Sequence seqJ(alignDB.get(j).getName(), alignDB.get(j).getAligned());
671 alignment->align(seqI.getUnaligned(), seqJ.getUnaligned());
672 seqI.setAligned(alignment->getSeqAAln());
673 seqJ.setAligned(alignment->getSeqBAln());
675 //cout << seqI.getName() << '\t' << seqJ.getName() << endl;
676 //cout << alignment->getSeqAAln() << endl << alignment->getSeqBAln() << endl;
678 distCalculator->calcDist(seqI, seqJ);
679 double dist = distCalculator->getDist();
681 //cout << "dist = " << dist << endl;
684 if (output == "column") { outFile << alignDB.get(i).getName() << ' ' << alignDB.get(j).getName() << ' ' << dist << endl; }
686 if (output == "lt") { outFile << dist << '\t'; }
689 if (output == "lt") { outFile << endl; }
692 m->mothurOut(toString(i) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
696 m->mothurOut(toString(endLine-1) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
700 delete distCalculator;
704 catch(exception& e) {
705 m->errorOut(e, "PairwiseSeqsCommand", "driver");
709 /**************************************************************************************************/
710 /////// need to fix to work with calcs and sequencedb
711 int PairwiseSeqsCommand::driver(int startLine, int endLine, string dFileName, string square){
714 int startTime = time(NULL);
716 Alignment* alignment;
717 if(align == "gotoh") { alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase); }
718 else if(align == "needleman") { alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase); }
719 else if(align == "blast") { alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch); }
720 else if(align == "noalign") { alignment = new NoAlign(); }
722 m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
723 m->mothurOutEndLine();
724 alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
727 ValidCalculators validCalculator;
728 Dist* distCalculator;
730 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
731 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
732 else if (Estimators[0] == "eachgap") { distCalculator = new eachGapDist(); }
733 else if (Estimators[0] == "onegap") { distCalculator = new oneGapDist(); }
736 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
737 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
738 else if (Estimators[0] == "eachgap"){ distCalculator = new eachGapIgnoreTermGapDist(); }
739 else if (Estimators[0] == "onegap") { distCalculator = new oneGapIgnoreTermGapDist(); }
744 ofstream outFile(dFileName.c_str(), ios::trunc);
745 outFile.setf(ios::fixed, ios::showpoint);
746 outFile << setprecision(4);
748 if(startLine == 0){ outFile << alignDB.getNumSeqs() << endl; }
750 for(int i=startLine;i<endLine;i++){
752 string name = alignDB.get(i).getName();
753 //pad with spaces to make compatible
754 if (name.length() < 10) { while (name.length() < 10) { name += " "; } }
756 outFile << name << '\t';
758 for(int j=0;j<alignDB.getNumSeqs();j++){
760 if (m->control_pressed) { outFile.close(); delete alignment; delete distCalculator; return 0; }
762 if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
763 alignment->resize(alignDB.get(i).getUnaligned().length()+1);
766 if (alignDB.get(j).getUnaligned().length() > alignment->getnRows()) {
767 alignment->resize(alignDB.get(j).getUnaligned().length()+1);
770 Sequence seqI(alignDB.get(i).getName(), alignDB.get(i).getAligned());
771 Sequence seqJ(alignDB.get(j).getName(), alignDB.get(j).getAligned());
773 alignment->align(seqI.getUnaligned(), seqJ.getUnaligned());
774 seqI.setAligned(alignment->getSeqAAln());
775 seqJ.setAligned(alignment->getSeqBAln());
777 distCalculator->calcDist(seqI, seqJ);
778 double dist = distCalculator->getDist();
780 outFile << dist << '\t';
786 m->mothurOut(toString(i) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
790 m->mothurOut(toString(endLine-1) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
794 delete distCalculator;
798 catch(exception& e) {
799 m->errorOut(e, "PairwiseSeqsCommand", "driver");
804 /**************************************************************************************************/
805 /////// need to fix to work with calcs and sequencedb
806 int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, MPI_File& outMPI, float cutoff){
809 int startTime = time(NULL);
811 Alignment* alignment;
812 if(align == "gotoh") { alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase); }
813 else if(align == "needleman") { alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase); }
814 else if(align == "blast") { alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch); }
815 else if(align == "noalign") { alignment = new NoAlign(); }
817 m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
818 m->mothurOutEndLine();
819 alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
822 ValidCalculators validCalculator;
823 Dist* distCalculator;
825 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
826 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
827 else if (Estimators[0] == "eachgap") { distCalculator = new eachGapDist(); }
828 else if (Estimators[0] == "onegap") { distCalculator = new oneGapDist(); }
831 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
832 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
833 else if (Estimators[0] == "eachgap"){ distCalculator = new eachGapIgnoreTermGapDist(); }
834 else if (Estimators[0] == "onegap") { distCalculator = new oneGapIgnoreTermGapDist(); }
838 string outputString = "";
840 for(int i=startLine;i<endLine;i++){
842 for(int j=0;j<i;j++){
844 if (m->control_pressed) { delete alignment; delete distCalculator; return 0; }
846 if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
847 alignment->resize(alignDB.get(i).getUnaligned().length()+1);
850 if (alignDB.get(j).getUnaligned().length() > alignment->getnRows()) {
851 alignment->resize(alignDB.get(j).getUnaligned().length()+1);
854 Sequence seqI(alignDB.get(i).getName(), alignDB.get(i).getAligned());
855 Sequence seqJ(alignDB.get(j).getName(), alignDB.get(j).getAligned());
857 alignment->align(seqI.getUnaligned(), seqJ.getUnaligned());
858 seqI.setAligned(alignment->getSeqAAln());
859 seqJ.setAligned(alignment->getSeqBAln());
861 distCalculator->calcDist(seqI, seqJ);
862 double dist = distCalculator->getDist();
865 outputString += (alignDB.get(i).getName() + ' ' + alignDB.get(j).getName() + ' ' + toString(dist) + '\n');
870 //m->mothurOut(toString(i) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
871 cout << i << '\t' << (time(NULL) - startTime) << endl;
875 //send results to parent
876 int length = outputString.length();
878 char* buf = new char[length];
879 memcpy(buf, outputString.c_str(), length);
881 MPI_File_write_shared(outMPI, buf, length, MPI_CHAR, &status);
887 delete distCalculator;
890 catch(exception& e) {
891 m->errorOut(e, "PairwiseSeqsCommand", "driverMPI");
895 /**************************************************************************************************/
896 /////// need to fix to work with calcs and sequencedb
897 int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, string file, unsigned long long& size){
902 int amode=MPI_MODE_CREATE|MPI_MODE_WRONLY;
905 strcpy(filename, file.c_str());
907 MPI_File_open(MPI_COMM_SELF, filename, amode, MPI_INFO_NULL, &outMPI);
909 Alignment* alignment;
910 if(align == "gotoh") { alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase); }
911 else if(align == "needleman") { alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase); }
912 else if(align == "blast") { alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch); }
913 else if(align == "noalign") { alignment = new NoAlign(); }
915 m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
916 m->mothurOutEndLine();
917 alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
920 ValidCalculators validCalculator;
921 Dist* distCalculator;
923 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
924 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
925 else if (Estimators[0] == "eachgap") { distCalculator = new eachGapDist(); }
926 else if (Estimators[0] == "onegap") { distCalculator = new oneGapDist(); }
929 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
930 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
931 else if (Estimators[0] == "eachgap"){ distCalculator = new eachGapIgnoreTermGapDist(); }
932 else if (Estimators[0] == "onegap") { distCalculator = new oneGapIgnoreTermGapDist(); }
937 string outputString = "";
940 if(startLine == 0){ outputString += toString(alignDB.getNumSeqs()) + "\n"; }
942 for(int i=startLine;i<endLine;i++){
944 string name = alignDB.get(i).getName();
945 if (name.length() < 10) { //pad with spaces to make compatible
946 while (name.length() < 10) { name += " "; }
948 outputString += name + "\t";
950 for(int j=0;j<i;j++){
952 if (m->control_pressed) { delete alignment; delete distCalculator; return 0; }
954 if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
955 alignment->resize(alignDB.get(i).getUnaligned().length()+1);
958 if (alignDB.get(j).getUnaligned().length() > alignment->getnRows()) {
959 alignment->resize(alignDB.get(j).getUnaligned().length()+1);
962 Sequence seqI(alignDB.get(i).getName(), alignDB.get(i).getAligned());
963 Sequence seqJ(alignDB.get(j).getName(), alignDB.get(j).getAligned());
965 alignment->align(seqI.getUnaligned(), seqJ.getUnaligned());
966 seqI.setAligned(alignment->getSeqAAln());
967 seqJ.setAligned(alignment->getSeqBAln());
969 distCalculator->calcDist(seqI, seqJ);
970 double dist = distCalculator->getDist();
972 outputString += toString(dist) + "\t";
975 outputString += "\n";
977 //send results to parent
978 int length = outputString.length();
979 char* buf = new char[length];
980 memcpy(buf, outputString.c_str(), length);
982 MPI_File_write(outMPI, buf, length, MPI_CHAR, &status);
983 size += outputString.length();
988 MPI_File_close(&outMPI);
990 delete distCalculator;
994 catch(exception& e) {
995 m->errorOut(e, "PairwiseSeqsCommand", "driverMPI");
999 /**************************************************************************************************/
1000 /////// need to fix to work with calcs and sequencedb
1001 int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, string file, unsigned long long& size, string square){
1006 int amode=MPI_MODE_CREATE|MPI_MODE_WRONLY;
1008 char filename[1024];
1009 strcpy(filename, file.c_str());
1011 MPI_File_open(MPI_COMM_SELF, filename, amode, MPI_INFO_NULL, &outMPI);
1013 Alignment* alignment;
1014 if(align == "gotoh") { alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase); }
1015 else if(align == "needleman") { alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase); }
1016 else if(align == "blast") { alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch); }
1017 else if(align == "noalign") { alignment = new NoAlign(); }
1019 m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
1020 m->mothurOutEndLine();
1021 alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
1024 ValidCalculators validCalculator;
1025 Dist* distCalculator;
1027 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
1028 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
1029 else if (Estimators[0] == "eachgap") { distCalculator = new eachGapDist(); }
1030 else if (Estimators[0] == "onegap") { distCalculator = new oneGapDist(); }
1033 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
1034 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
1035 else if (Estimators[0] == "eachgap"){ distCalculator = new eachGapIgnoreTermGapDist(); }
1036 else if (Estimators[0] == "onegap") { distCalculator = new oneGapIgnoreTermGapDist(); }
1040 string outputString = "";
1043 if(startLine == 0){ outputString += toString(alignDB.getNumSeqs()) + "\n"; }
1045 for(int i=startLine;i<endLine;i++){
1047 string name = alignDB.get(i).getName();
1048 if (name.length() < 10) { //pad with spaces to make compatible
1049 while (name.length() < 10) { name += " "; }
1051 outputString += name + "\t";
1053 for(int j=0;j<alignDB.getNumSeqs();j++){
1055 if (m->control_pressed) { delete alignment; return 0; }
1057 if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
1058 alignment->resize(alignDB.get(i).getUnaligned().length()+1);
1061 if (alignDB.get(j).getUnaligned().length() > alignment->getnRows()) {
1062 alignment->resize(alignDB.get(j).getUnaligned().length()+1);
1065 Sequence seqI(alignDB.get(i).getName(), alignDB.get(i).getAligned());
1066 Sequence seqJ(alignDB.get(j).getName(), alignDB.get(j).getAligned());
1068 alignment->align(seqI.getUnaligned(), seqJ.getUnaligned());
1069 seqI.setAligned(alignment->getSeqAAln());
1070 seqJ.setAligned(alignment->getSeqBAln());
1072 distCalculator->calcDist(seqI, seqJ);
1073 double dist = distCalculator->getDist();
1075 outputString += toString(dist) + "\t";
1078 outputString += "\n";
1080 //send results to parent
1081 int length = outputString.length();
1082 char* buf = new char[length];
1083 memcpy(buf, outputString.c_str(), length);
1085 MPI_File_write(outMPI, buf, length, MPI_CHAR, &status);
1086 size += outputString.length();
1091 MPI_File_close(&outMPI);
1096 catch(exception& e) {
1097 m->errorOut(e, "PairwiseSeqsCommand", "driverMPI");
1102 /**************************************************************************************************/