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",false,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); parameters.push_back(pprocessors);
22 CommandParameter poutput("output", "Multiple", "column-lt-square-phylip", "column", "", "", "",false,false); 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); 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::getOutputFileNameTag(string type, string inputName=""){
70 string outputFileName = "";
71 map<string, vector<string> >::iterator it;
73 //is this a type this command creates
74 it = outputTypes.find(type);
75 if (it == outputTypes.end()) { m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); }
77 if (type == "phylip") { outputFileName = "dist"; }
78 else if (type == "column") { outputFileName = "dist"; }
79 else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true; }
81 return outputFileName;
84 m->errorOut(e, "PairwiseSeqsCommand", "getOutputFileNameTag");
88 //**********************************************************************************************************************
89 PairwiseSeqsCommand::PairwiseSeqsCommand(){
91 abort = true; calledHelp = true;
93 vector<string> tempOutNames;
94 outputTypes["phylip"] = tempOutNames;
95 outputTypes["column"] = tempOutNames;
98 m->errorOut(e, "PairwiseSeqsCommand", "PairwiseSeqsCommand");
102 //**********************************************************************************************************************
103 PairwiseSeqsCommand::PairwiseSeqsCommand(string option) {
105 abort = false; calledHelp = false;
107 //allow user to run help
108 if(option == "help") { help(); abort = true; calledHelp = true; }
109 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
112 vector<string> myArray = setParameters();
114 OptionParser parser(option);
115 map<string, string> parameters = parser.getParameters();
117 ValidParameters validParameter("pairwise.seqs");
118 map<string, string>::iterator it;
120 //check to make sure all parameters are valid for command
121 for (it = parameters.begin(); it != parameters.end(); it++) {
122 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
125 //initialize outputTypes
126 vector<string> tempOutNames;
127 outputTypes["phylip"] = tempOutNames;
128 outputTypes["column"] = tempOutNames;
130 //if the user changes the output directory command factory will send this info to us in the output parameter
131 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; }
134 //if the user changes the input directory command factory will send this info to us in the output parameter
135 string inputDir = validParameter.validFile(parameters, "inputdir", false);
137 if (inputDir == "not found"){ inputDir = ""; }
139 fastaFileName = validParameter.validFile(parameters, "fasta", false);
140 if (fastaFileName == "not found") {
141 //if there is a current fasta file, use it
142 string filename = m->getFastaFile();
143 if (filename != "") { fastaFileNames.push_back(filename); m->mothurOut("Using " + filename + " as input file for the fasta parameter."); m->mothurOutEndLine(); }
144 else { m->mothurOut("You have no current fastafile and the fasta parameter is required."); m->mothurOutEndLine(); abort = true; }
147 m->splitAtDash(fastaFileName, fastaFileNames);
149 //go through files and make sure they are good, if not, then disregard them
150 for (int i = 0; i < fastaFileNames.size(); i++) {
153 if (fastaFileNames[i] == "current") {
154 fastaFileNames[i] = m->getFastaFile();
155 if (fastaFileNames[i] != "") { m->mothurOut("Using " + fastaFileNames[i] + " as input file for the fasta parameter where you had given current."); m->mothurOutEndLine(); }
157 m->mothurOut("You have no current fastafile, ignoring current."); m->mothurOutEndLine(); ignore=true;
158 //erase from file list
159 fastaFileNames.erase(fastaFileNames.begin()+i);
166 if (inputDir != "") {
167 string path = m->hasPath(fastaFileNames[i]);
168 //if the user has not given a path then, add inputdir. else leave path alone.
169 if (path == "") { fastaFileNames[i] = inputDir + fastaFileNames[i]; }
175 ableToOpen = m->openInputFile(fastaFileNames[i], in, "noerror");
177 //if you can't open it, try default location
178 if (ableToOpen == 1) {
179 if (m->getDefaultPath() != "") { //default path is set
180 string tryPath = m->getDefaultPath() + m->getSimpleName(fastaFileNames[i]);
181 m->mothurOut("Unable to open " + fastaFileNames[i] + ". Trying default " + tryPath); m->mothurOutEndLine();
183 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
185 fastaFileNames[i] = tryPath;
189 //if you can't open it, try output location
190 if (ableToOpen == 1) {
191 if (m->getOutputDir() != "") { //default path is set
192 string tryPath = m->getOutputDir() + m->getSimpleName(fastaFileNames[i]);
193 m->mothurOut("Unable to open " + fastaFileNames[i] + ". Trying output directory " + tryPath); m->mothurOutEndLine();
195 ableToOpen = m->openInputFile(tryPath, in2, "noerror");
197 fastaFileNames[i] = tryPath;
203 if (ableToOpen == 1) {
204 m->mothurOut("Unable to open " + fastaFileNames[i] + ". It will be disregarded."); m->mothurOutEndLine();
205 //erase from file list
206 fastaFileNames.erase(fastaFileNames.begin()+i);
209 m->setFastaFile(fastaFileNames[i]);
214 //make sure there is at least one valid file left
215 if (fastaFileNames.size() == 0) { m->mothurOut("no valid files."); m->mothurOutEndLine(); abort = true; }
218 //check for optional parameter and set defaults
219 // ...at some point should added some additional type checking...
221 temp = validParameter.validFile(parameters, "match", false); if (temp == "not found"){ temp = "1.0"; }
222 m->mothurConvert(temp, match);
224 temp = validParameter.validFile(parameters, "mismatch", false); if (temp == "not found"){ temp = "-1.0"; }
225 m->mothurConvert(temp, misMatch);
226 if (misMatch > 0) { m->mothurOut("[ERROR]: mismatch must be negative.\n"); abort=true; }
228 temp = validParameter.validFile(parameters, "gapopen", false); if (temp == "not found"){ temp = "-2.0"; }
229 m->mothurConvert(temp, gapOpen);
230 if (gapOpen > 0) { m->mothurOut("[ERROR]: gapopen must be negative.\n"); abort=true; }
232 temp = validParameter.validFile(parameters, "gapextend", false); if (temp == "not found"){ temp = "-1.0"; }
233 m->mothurConvert(temp, gapExtend);
234 if (gapExtend > 0) { m->mothurOut("[ERROR]: gapextend must be negative.\n"); abort=true; }
236 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
237 m->setProcessors(temp);
238 m->mothurConvert(temp, processors);
240 temp = validParameter.validFile(parameters, "cutoff", false); if(temp == "not found"){ temp = "1.0"; }
241 m->mothurConvert(temp, cutoff);
243 temp = validParameter.validFile(parameters, "countends", false); if(temp == "not found"){ temp = "T"; }
244 countends = m->isTrue(temp);
246 temp = validParameter.validFile(parameters, "compress", false); if(temp == "not found"){ temp = "F"; }
247 compress = m->isTrue(temp);
249 align = validParameter.validFile(parameters, "align", false); if (align == "not found"){ align = "needleman"; }
251 output = validParameter.validFile(parameters, "output", false); if(output == "not found"){ output = "column"; }
252 if (output=="phylip") { output = "lt"; }
253 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"; }
255 calc = validParameter.validFile(parameters, "calc", false);
256 if (calc == "not found") { calc = "onegap"; }
258 if (calc == "default") { calc = "onegap"; }
260 m->splitAtDash(calc, Estimators);
264 catch(exception& e) {
265 m->errorOut(e, "PairwiseSeqsCommand", "PairwiseSeqsCommand");
269 //**********************************************************************************************************************
271 int PairwiseSeqsCommand::execute(){
273 if (abort == true) { if (calledHelp) { return 0; } return 2; }
275 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.
279 for (int s = 0; s < fastaFileNames.size(); s++) {
280 if (m->control_pressed) { outputTypes.clear(); return 0; }
282 m->mothurOut("Processing sequences from " + fastaFileNames[s] + " ..." ); m->mothurOutEndLine();
284 if (outputDir == "") { outputDir += m->hasPath(fastaFileNames[s]); }
287 m->openInputFile(fastaFileNames[s], inFASTA);
288 alignDB = SequenceDB(inFASTA);
291 int numSeqs = alignDB.getNumSeqs();
292 int startTime = time(NULL);
293 string outputFile = "";
295 if (output == "lt") { //does the user want lower triangle phylip formatted file
296 outputFile = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[s])) + "phylip." + getOutputFileNameTag("phylip");
297 m->mothurRemove(outputFile); outputTypes["phylip"].push_back(outputFile);
298 }else if (output == "column") { //user wants column format
299 outputFile = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[s])) + getOutputFileNameTag("column");
300 outputTypes["column"].push_back(outputFile);
301 m->mothurRemove(outputFile);
302 }else { //assume square
303 outputFile = outputDir + m->getRootName(m->getSimpleName(fastaFileNames[s])) + "square." + getOutputFileNameTag("phylip");
304 m->mothurRemove(outputFile);
305 outputTypes["phylip"].push_back(outputFile);
314 MPI_Comm_size(MPI_COMM_WORLD, &processors); //set processors to the number of mpi processes running
315 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
317 //each process gets where it should start and stop in the file
318 if (output != "square") {
319 start = int (sqrt(float(pid)/float(processors)) * numSeqs);
320 end = int (sqrt(float(pid+1)/float(processors)) * numSeqs);
322 start = int ((float(pid)/float(processors)) * numSeqs);
323 end = int ((float(pid+1)/float(processors)) * numSeqs);
326 if (output == "column") {
328 int amode=MPI_MODE_CREATE|MPI_MODE_WRONLY;
331 strcpy(filename, outputFile.c_str());
333 MPI_File_open(MPI_COMM_WORLD, filename, amode, MPI_INFO_NULL, &outMPI);
335 if (pid == 0) { //you are the root process
340 driverMPI(start, end, outMPI, cutoff);
342 if (m->control_pressed) { outputTypes.clear(); MPI_File_close(&outMPI); m->mothurRemove(outputFile); return 0; }
345 for(int i = 1; i < processors; i++) {
346 if (m->control_pressed) { outputTypes.clear(); MPI_File_close(&outMPI); m->mothurRemove(outputFile); return 0; }
349 MPI_Recv(buf, 5, MPI_CHAR, i, tag, MPI_COMM_WORLD, &status);
351 }else { //you are a child process
353 driverMPI(start, end, outMPI, cutoff);
355 if (m->control_pressed) { outputTypes.clear(); MPI_File_close(&outMPI); m->mothurRemove(outputFile); return 0; }
359 //tell parent you are done.
360 MPI_Send(buf, 5, MPI_CHAR, 0, tag, MPI_COMM_WORLD);
363 MPI_File_close(&outMPI);
365 }else { //lower triangle format
366 if (pid == 0) { //you are the root process
370 unsigned long long mySize;
372 if (output != "square"){ driverMPI(start, end, outputFile, mySize); }
373 else { driverMPI(start, end, outputFile, mySize, output); }
375 if (m->control_pressed) { outputTypes.clear(); m->mothurRemove(outputFile); return 0; }
377 int amode=MPI_MODE_APPEND|MPI_MODE_WRONLY|MPI_MODE_CREATE; //
382 strcpy(filename, outputFile.c_str());
384 MPI_File_open(MPI_COMM_SELF, filename, amode, MPI_INFO_NULL, &outMPI);
387 for(int b = 1; b < processors; b++) {
388 unsigned long long fileSize;
390 if (m->control_pressed) { outputTypes.clear(); MPI_File_close(&outMPI); m->mothurRemove(outputFile); return 0; }
392 MPI_Recv(&fileSize, 1, MPI_LONG, b, tag, MPI_COMM_WORLD, &status);
394 string outTemp = outputFile + toString(b) + ".temp";
396 char* buf = new char[outTemp.length()];
397 memcpy(buf, outTemp.c_str(), outTemp.length());
399 MPI_File_open(MPI_COMM_SELF, buf, MPI_MODE_DELETE_ON_CLOSE|MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);
403 while (count < fileSize) {
405 MPI_File_read(inMPI, buf2, 1, MPI_CHAR, &status);
406 MPI_File_write(outMPI, buf2, 1, MPI_CHAR, &status);
410 MPI_File_close(&inMPI); //deleted on close
413 MPI_File_close(&outMPI);
414 }else { //you are a child process
416 unsigned long long size;
417 if (output != "square"){ driverMPI(start, end, (outputFile + toString(pid) + ".temp"), size); }
418 else { driverMPI(start, end, (outputFile + toString(pid) + ".temp"), size, output); }
420 if (m->control_pressed) { return 0; }
422 //tell parent you are done.
423 MPI_Send(&size, 1, MPI_LONG, 0, tag, MPI_COMM_WORLD);
426 MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
429 //#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
430 //if you don't need to fork anything
432 if (output != "square") { driver(0, numSeqs, outputFile, cutoff); }
433 else { driver(0, numSeqs, outputFile, "square"); }
434 }else{ //you have multiple processors
436 for (int i = 0; i < processors; i++) {
437 distlinePair tempLine;
438 lines.push_back(tempLine);
439 if (output != "square") {
440 lines[i].start = int (sqrt(float(i)/float(processors)) * numSeqs);
441 lines[i].end = int (sqrt(float(i+1)/float(processors)) * numSeqs);
443 lines[i].start = int ((float(i)/float(processors)) * numSeqs);
444 lines[i].end = int ((float(i+1)/float(processors)) * numSeqs);
448 createProcesses(outputFile);
452 //if (output != "square") { driver(0, numSeqs, outputFile, cutoff); }
453 //else { driver(0, numSeqs, outputFile, "square"); }
457 if (m->control_pressed) { outputTypes.clear(); m->mothurRemove(outputFile); return 0; }
460 MPI_Comm_rank(MPI_COMM_WORLD, &pid);
462 if (pid == 0) { //only one process should output to screen
466 fileHandle.open(outputFile.c_str());
468 m->gobble(fileHandle);
469 if (fileHandle.eof()) { m->mothurOut(outputFile + " is blank. This can result if there are no distances below your cutoff."); m->mothurOutEndLine(); }
473 m->mothurOut("Compressing..."); m->mothurOutEndLine();
474 m->mothurOut("(Replacing " + outputFile + " with " + outputFile + ".gz)"); m->mothurOutEndLine();
475 system(("gzip -v " + outputFile).c_str());
476 outputNames.push_back(outputFile + ".gz");
477 }else { outputNames.push_back(outputFile); }
483 m->mothurOut("It took " + toString(time(NULL) - startTime) + " to calculate the distances for " + toString(numSeqs) + " sequences."); m->mothurOutEndLine();
485 if (m->control_pressed) { outputTypes.clear(); m->mothurRemove(outputFile); return 0; }
488 //set phylip file as new current phylipfile
490 itTypes = outputTypes.find("phylip");
491 if (itTypes != outputTypes.end()) {
492 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setPhylipFile(current); }
495 //set column file as new current columnfile
496 itTypes = outputTypes.find("column");
497 if (itTypes != outputTypes.end()) {
498 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setColumnFile(current); }
501 m->mothurOutEndLine();
502 m->mothurOut("Output File Name: "); m->mothurOutEndLine();
503 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
504 m->mothurOutEndLine();
509 catch(exception& e) {
510 m->errorOut(e, "PairwiseSeqsCommand", "execute");
515 /**************************************************************************************************/
516 void PairwiseSeqsCommand::createProcesses(string filename) {
521 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
524 //loop through and create all the processes you want
525 while (process != processors) {
529 processIDS.push_back(pid);
532 if (output != "square") { driver(lines[process].start, lines[process].end, filename + toString(getpid()) + ".temp", cutoff); }
533 else { driver(lines[process].start, lines[process].end, filename + toString(getpid()) + ".temp", "square"); }
536 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
537 for (int i=0;i<processIDS.size();i++) { int temp = processIDS[i]; kill (temp, SIGINT); }
543 if (output != "square") { driver(lines[0].start, lines[0].end, filename, cutoff); }
544 else { driver(lines[0].start, lines[0].end, filename, "square"); }
547 //force parent to wait until all the processes are done
548 for (int i=0;i<processIDS.size();i++) {
549 int temp = processIDS[i];
553 //////////////////////////////////////////////////////////////////////////////////////////////////////
554 //Windows version shared memory, so be careful when passing variables through the distanceData struct.
555 //Above fork() will clone, so memory is separate, but that's not the case with windows,
556 //that's why the distance calculator was moved inside of the driver to make separate copies.
557 //////////////////////////////////////////////////////////////////////////////////////////////////////
559 vector<pairwiseData*> pDataArray; //[processors-1];
560 DWORD dwThreadIdArray[processors-1];
561 HANDLE hThreadArray[processors-1];
563 //Create processor-1 worker threads.
564 for( int i=0; i<processors-1; i++ ){
565 string extension = toString(i) + ".temp";
567 // Allocate memory for thread data.
568 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);
569 pDataArray.push_back(tempDist);
570 processIDS.push_back(i);
572 if (output != "square") { hThreadArray[i] = CreateThread(NULL, 0, MyPairwiseThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]); }
573 else { hThreadArray[i] = CreateThread(NULL, 0, MyPairwiseSquareThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]); }
577 if (output != "square") { driver(lines[0].start, lines[0].end, filename, cutoff); }
578 else { driver(lines[0].start, lines[0].end, filename, "square"); }
580 //Wait until all threads have terminated.
581 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
583 //Close all thread handles and free memory allocations.
584 for(int i=0; i < pDataArray.size(); i++){
585 CloseHandle(hThreadArray[i]);
586 delete pDataArray[i];
591 //append and remove temp files
592 for (int i=0;i<processIDS.size();i++) {
593 m->appendFiles((filename + toString(processIDS[i]) + ".temp"), filename);
594 m->mothurRemove((filename + toString(processIDS[i]) + ".temp"));
598 catch(exception& e) {
599 m->errorOut(e, "PairwiseSeqsCommand", "createProcesses");
604 /**************************************************************************************************/
605 /////// need to fix to work with calcs and sequencedb
606 int PairwiseSeqsCommand::driver(int startLine, int endLine, string dFileName, float cutoff){
609 int startTime = time(NULL);
611 Alignment* alignment;
612 if(align == "gotoh") { alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase); }
613 else if(align == "needleman") { alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase); }
614 else if(align == "blast") { alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch); }
615 else if(align == "noalign") { alignment = new NoAlign(); }
617 m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
618 m->mothurOutEndLine();
619 alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
622 ValidCalculators validCalculator;
623 Dist* distCalculator;
625 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
626 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
627 else if (Estimators[0] == "eachgap") { distCalculator = new eachGapDist(); }
628 else if (Estimators[0] == "onegap") { distCalculator = new oneGapDist(); }
631 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
632 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
633 else if (Estimators[0] == "eachgap"){ distCalculator = new eachGapIgnoreTermGapDist(); }
634 else if (Estimators[0] == "onegap") { distCalculator = new oneGapIgnoreTermGapDist(); }
639 ofstream outFile(dFileName.c_str(), ios::trunc);
640 outFile.setf(ios::fixed, ios::showpoint);
641 outFile << setprecision(4);
643 if((output == "lt") && startLine == 0){ outFile << alignDB.getNumSeqs() << endl; }
645 for(int i=startLine;i<endLine;i++){
647 string name = alignDB.get(i).getName();
648 if (name.length() < 10) { //pad with spaces to make compatible
649 while (name.length() < 10) { name += " "; }
651 outFile << name << '\t';
654 for(int j=0;j<i;j++){
656 if (m->control_pressed) { outFile.close(); delete alignment; delete distCalculator; return 0; }
658 if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
659 alignment->resize(alignDB.get(i).getUnaligned().length()+1);
662 if (alignDB.get(j).getUnaligned().length() > alignment->getnRows()) {
663 alignment->resize(alignDB.get(j).getUnaligned().length()+1);
666 Sequence seqI(alignDB.get(i).getName(), alignDB.get(i).getAligned());
667 Sequence seqJ(alignDB.get(j).getName(), alignDB.get(j).getAligned());
669 alignment->align(seqI.getUnaligned(), seqJ.getUnaligned());
670 seqI.setAligned(alignment->getSeqAAln());
671 seqJ.setAligned(alignment->getSeqBAln());
674 distCalculator->calcDist(seqI, seqJ);
675 double dist = distCalculator->getDist();
678 if (output == "column") { outFile << alignDB.get(i).getName() << ' ' << alignDB.get(j).getName() << ' ' << dist << endl; }
680 if (output == "lt") { outFile << dist << '\t'; }
683 if (output == "lt") { outFile << endl; }
686 m->mothurOut(toString(i) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
690 m->mothurOut(toString(endLine-1) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
694 delete distCalculator;
698 catch(exception& e) {
699 m->errorOut(e, "PairwiseSeqsCommand", "driver");
703 /**************************************************************************************************/
704 /////// need to fix to work with calcs and sequencedb
705 int PairwiseSeqsCommand::driver(int startLine, int endLine, string dFileName, string square){
708 int startTime = time(NULL);
710 Alignment* alignment;
711 if(align == "gotoh") { alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase); }
712 else if(align == "needleman") { alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase); }
713 else if(align == "blast") { alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch); }
714 else if(align == "noalign") { alignment = new NoAlign(); }
716 m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
717 m->mothurOutEndLine();
718 alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
721 ValidCalculators validCalculator;
722 Dist* distCalculator;
724 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
725 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
726 else if (Estimators[0] == "eachgap") { distCalculator = new eachGapDist(); }
727 else if (Estimators[0] == "onegap") { distCalculator = new oneGapDist(); }
730 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
731 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
732 else if (Estimators[0] == "eachgap"){ distCalculator = new eachGapIgnoreTermGapDist(); }
733 else if (Estimators[0] == "onegap") { distCalculator = new oneGapIgnoreTermGapDist(); }
738 ofstream outFile(dFileName.c_str(), ios::trunc);
739 outFile.setf(ios::fixed, ios::showpoint);
740 outFile << setprecision(4);
742 if(startLine == 0){ outFile << alignDB.getNumSeqs() << endl; }
744 for(int i=startLine;i<endLine;i++){
746 string name = alignDB.get(i).getName();
747 //pad with spaces to make compatible
748 if (name.length() < 10) { while (name.length() < 10) { name += " "; } }
750 outFile << name << '\t';
752 for(int j=0;j<alignDB.getNumSeqs();j++){
754 if (m->control_pressed) { outFile.close(); delete alignment; delete distCalculator; return 0; }
756 if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
757 alignment->resize(alignDB.get(i).getUnaligned().length()+1);
760 if (alignDB.get(j).getUnaligned().length() > alignment->getnRows()) {
761 alignment->resize(alignDB.get(j).getUnaligned().length()+1);
764 Sequence seqI(alignDB.get(i).getName(), alignDB.get(i).getAligned());
765 Sequence seqJ(alignDB.get(j).getName(), alignDB.get(j).getAligned());
767 alignment->align(seqI.getUnaligned(), seqJ.getUnaligned());
768 seqI.setAligned(alignment->getSeqAAln());
769 seqJ.setAligned(alignment->getSeqBAln());
771 distCalculator->calcDist(seqI, seqJ);
772 double dist = distCalculator->getDist();
774 outFile << dist << '\t';
780 m->mothurOut(toString(i) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
784 m->mothurOut(toString(endLine-1) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
788 delete distCalculator;
792 catch(exception& e) {
793 m->errorOut(e, "PairwiseSeqsCommand", "driver");
798 /**************************************************************************************************/
799 /////// need to fix to work with calcs and sequencedb
800 int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, MPI_File& outMPI, float cutoff){
803 int startTime = time(NULL);
805 Alignment* alignment;
806 if(align == "gotoh") { alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase); }
807 else if(align == "needleman") { alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase); }
808 else if(align == "blast") { alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch); }
809 else if(align == "noalign") { alignment = new NoAlign(); }
811 m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
812 m->mothurOutEndLine();
813 alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
816 ValidCalculators validCalculator;
817 Dist* distCalculator;
819 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
820 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
821 else if (Estimators[0] == "eachgap") { distCalculator = new eachGapDist(); }
822 else if (Estimators[0] == "onegap") { distCalculator = new oneGapDist(); }
825 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
826 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
827 else if (Estimators[0] == "eachgap"){ distCalculator = new eachGapIgnoreTermGapDist(); }
828 else if (Estimators[0] == "onegap") { distCalculator = new oneGapIgnoreTermGapDist(); }
832 string outputString = "";
834 for(int i=startLine;i<endLine;i++){
836 for(int j=0;j<i;j++){
838 if (m->control_pressed) { delete alignment; delete distCalculator; return 0; }
840 if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
841 alignment->resize(alignDB.get(i).getUnaligned().length()+1);
844 if (alignDB.get(j).getUnaligned().length() > alignment->getnRows()) {
845 alignment->resize(alignDB.get(j).getUnaligned().length()+1);
848 Sequence seqI(alignDB.get(i).getName(), alignDB.get(i).getAligned());
849 Sequence seqJ(alignDB.get(j).getName(), alignDB.get(j).getAligned());
851 alignment->align(seqI.getUnaligned(), seqJ.getUnaligned());
852 seqI.setAligned(alignment->getSeqAAln());
853 seqJ.setAligned(alignment->getSeqBAln());
855 distCalculator->calcDist(seqI, seqJ);
856 double dist = distCalculator->getDist();
859 outputString += (alignDB.get(i).getName() + ' ' + alignDB.get(j).getName() + ' ' + toString(dist) + '\n');
864 //m->mothurOut(toString(i) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
865 cout << i << '\t' << (time(NULL) - startTime) << endl;
869 //send results to parent
870 int length = outputString.length();
872 char* buf = new char[length];
873 memcpy(buf, outputString.c_str(), length);
875 MPI_File_write_shared(outMPI, buf, length, MPI_CHAR, &status);
881 delete distCalculator;
884 catch(exception& e) {
885 m->errorOut(e, "PairwiseSeqsCommand", "driverMPI");
889 /**************************************************************************************************/
890 /////// need to fix to work with calcs and sequencedb
891 int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, string file, unsigned long long& size){
896 int amode=MPI_MODE_CREATE|MPI_MODE_WRONLY;
899 strcpy(filename, file.c_str());
901 MPI_File_open(MPI_COMM_SELF, filename, amode, MPI_INFO_NULL, &outMPI);
903 Alignment* alignment;
904 if(align == "gotoh") { alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase); }
905 else if(align == "needleman") { alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase); }
906 else if(align == "blast") { alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch); }
907 else if(align == "noalign") { alignment = new NoAlign(); }
909 m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
910 m->mothurOutEndLine();
911 alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
914 ValidCalculators validCalculator;
915 Dist* distCalculator;
917 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
918 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
919 else if (Estimators[0] == "eachgap") { distCalculator = new eachGapDist(); }
920 else if (Estimators[0] == "onegap") { distCalculator = new oneGapDist(); }
923 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
924 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
925 else if (Estimators[0] == "eachgap"){ distCalculator = new eachGapIgnoreTermGapDist(); }
926 else if (Estimators[0] == "onegap") { distCalculator = new oneGapIgnoreTermGapDist(); }
931 string outputString = "";
934 if(startLine == 0){ outputString += toString(alignDB.getNumSeqs()) + "\n"; }
936 for(int i=startLine;i<endLine;i++){
938 string name = alignDB.get(i).getName();
939 if (name.length() < 10) { //pad with spaces to make compatible
940 while (name.length() < 10) { name += " "; }
942 outputString += name + "\t";
944 for(int j=0;j<i;j++){
946 if (m->control_pressed) { delete alignment; delete distCalculator; return 0; }
948 if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
949 alignment->resize(alignDB.get(i).getUnaligned().length()+1);
952 if (alignDB.get(j).getUnaligned().length() > alignment->getnRows()) {
953 alignment->resize(alignDB.get(j).getUnaligned().length()+1);
956 Sequence seqI(alignDB.get(i).getName(), alignDB.get(i).getAligned());
957 Sequence seqJ(alignDB.get(j).getName(), alignDB.get(j).getAligned());
959 alignment->align(seqI.getUnaligned(), seqJ.getUnaligned());
960 seqI.setAligned(alignment->getSeqAAln());
961 seqJ.setAligned(alignment->getSeqBAln());
963 distCalculator->calcDist(seqI, seqJ);
964 double dist = distCalculator->getDist();
966 outputString += toString(dist) + "\t";
969 outputString += "\n";
971 //send results to parent
972 int length = outputString.length();
973 char* buf = new char[length];
974 memcpy(buf, outputString.c_str(), length);
976 MPI_File_write(outMPI, buf, length, MPI_CHAR, &status);
977 size += outputString.length();
982 MPI_File_close(&outMPI);
984 delete distCalculator;
988 catch(exception& e) {
989 m->errorOut(e, "PairwiseSeqsCommand", "driverMPI");
993 /**************************************************************************************************/
994 /////// need to fix to work with calcs and sequencedb
995 int PairwiseSeqsCommand::driverMPI(int startLine, int endLine, string file, unsigned long long& size, string square){
1000 int amode=MPI_MODE_CREATE|MPI_MODE_WRONLY;
1002 char filename[1024];
1003 strcpy(filename, file.c_str());
1005 MPI_File_open(MPI_COMM_SELF, filename, amode, MPI_INFO_NULL, &outMPI);
1007 Alignment* alignment;
1008 if(align == "gotoh") { alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase); }
1009 else if(align == "needleman") { alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase); }
1010 else if(align == "blast") { alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch); }
1011 else if(align == "noalign") { alignment = new NoAlign(); }
1013 m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
1014 m->mothurOutEndLine();
1015 alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
1018 ValidCalculators validCalculator;
1019 Dist* distCalculator;
1021 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
1022 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
1023 else if (Estimators[0] == "eachgap") { distCalculator = new eachGapDist(); }
1024 else if (Estimators[0] == "onegap") { distCalculator = new oneGapDist(); }
1027 if (validCalculator.isValidCalculator("distance", Estimators[0]) == true) {
1028 if (Estimators[0] == "nogaps") { distCalculator = new ignoreGaps(); }
1029 else if (Estimators[0] == "eachgap"){ distCalculator = new eachGapIgnoreTermGapDist(); }
1030 else if (Estimators[0] == "onegap") { distCalculator = new oneGapIgnoreTermGapDist(); }
1034 string outputString = "";
1037 if(startLine == 0){ outputString += toString(alignDB.getNumSeqs()) + "\n"; }
1039 for(int i=startLine;i<endLine;i++){
1041 string name = alignDB.get(i).getName();
1042 if (name.length() < 10) { //pad with spaces to make compatible
1043 while (name.length() < 10) { name += " "; }
1045 outputString += name + "\t";
1047 for(int j=0;j<alignDB.getNumSeqs();j++){
1049 if (m->control_pressed) { delete alignment; return 0; }
1051 if (alignDB.get(i).getUnaligned().length() > alignment->getnRows()) {
1052 alignment->resize(alignDB.get(i).getUnaligned().length()+1);
1055 if (alignDB.get(j).getUnaligned().length() > alignment->getnRows()) {
1056 alignment->resize(alignDB.get(j).getUnaligned().length()+1);
1059 Sequence seqI(alignDB.get(i).getName(), alignDB.get(i).getAligned());
1060 Sequence seqJ(alignDB.get(j).getName(), alignDB.get(j).getAligned());
1062 alignment->align(seqI.getUnaligned(), seqJ.getUnaligned());
1063 seqI.setAligned(alignment->getSeqAAln());
1064 seqJ.setAligned(alignment->getSeqBAln());
1066 distCalculator->calcDist(seqI, seqJ);
1067 double dist = distCalculator->getDist();
1069 outputString += toString(dist) + "\t";
1072 outputString += "\n";
1074 //send results to parent
1075 int length = outputString.length();
1076 char* buf = new char[length];
1077 memcpy(buf, outputString.c_str(), length);
1079 MPI_File_write(outMPI, buf, length, MPI_CHAR, &status);
1080 size += outputString.length();
1085 MPI_File_close(&outMPI);
1090 catch(exception& e) {
1091 m->errorOut(e, "PairwiseSeqsCommand", "driverMPI");
1096 /**************************************************************************************************/