]> git.donarmstrong.com Git - mothur.git/blobdiff - classifyseqscommand.cpp
reworked the classifiers summary file added groupfile option to classify.seqs, added...
[mothur.git] / classifyseqscommand.cpp
index a9f0a36a2a440e17e8f4f433954c766abba041ed..38378491c072c1ed33cc257ec954de6a81f3bb7d 100644 (file)
@@ -11,6 +11,7 @@
 #include "sequence.hpp"
 #include "bayesian.h"
 #include "phylotree.h"
+#include "phylosummary.h"
 #include "knn.h"
 
 //**********************************************************************************************************************
@@ -25,7 +26,7 @@ ClassifySeqsCommand::ClassifySeqsCommand(string option)  {
                else {
                        
                        //valid paramters for this command
-                       string AlignArray[] =  {"template","fasta","name","search","ksize","method","processors","taxonomy","match","mismatch","gapopen","gapextend","numwanted","cutoff","probs","iters", "outputdir","inputdir"};
+                       string AlignArray[] =  {"template","fasta","name","group","search","ksize","method","processors","taxonomy","match","mismatch","gapopen","gapextend","numwanted","cutoff","probs","iters", "outputdir","inputdir"};
                        vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
                        
                        OptionParser parser(option);
@@ -62,6 +63,14 @@ ClassifySeqsCommand::ClassifySeqsCommand(string option)  {
                                        //if the user has not given a path then, add inputdir. else leave path alone.
                                        if (path == "") {       parameters["taxonomy"] = inputDir + it->second;         }
                                }
+                               
+                               it = parameters.find("group");
+                               //user has given a template file
+                               if(it != parameters.end()){ 
+                                       path = hasPath(it->second);
+                                       //if the user has not given a path then, add inputdir. else leave path alone.
+                                       if (path == "") {       parameters["group"] = inputDir + it->second;            }
+                               }
                        }
 
                        //check for required parameters
@@ -73,6 +82,10 @@ ClassifySeqsCommand::ClassifySeqsCommand(string option)  {
                        }
                        else if (templateFileName == "not open") { abort = true; }      
                        
+                       groupfile = validParameter.validFile(parameters, "group", true);
+                       if (groupfile == "not open") { abort = true; }  
+                       else if (groupfile == "not found") { groupfile = ""; }
+                       
                        fastaFileName = validParameter.validFile(parameters, "fasta", false);
                        if (fastaFileName == "not found") { m->mothurOut("fasta is a required parameter for the classify.seqs command."); m->mothurOutEndLine(); abort = true;  }
                        else { 
@@ -250,6 +263,7 @@ void ClassifySeqsCommand::help(){
                m->mothurOut("The template, fasta and taxonomy parameters are required. You may enter multiple fasta files by separating their names with dashes. ie. fasta=abrecovery.fasta-amzon.fasta \n");
                m->mothurOut("The search parameter allows you to specify the method to find most similar template.  Your options are: suffix, kmer, blast and distance. The default is kmer.\n");
                m->mothurOut("The name parameter allows you add a names file with your fasta file, if you enter multiple fasta files, you must enter matching names files for them.\n");
+               m->mothurOut("The group parameter allows you add a group file so you can have the summary totals broken up by group.\n");
                m->mothurOut("The method parameter allows you to specify classification method to use.  Your options are: bayesian and knn. The default is bayesian.\n");
                m->mothurOut("The ksize parameter allows you to specify the kmer size for finding most similar template to candidate.  The default is 8.\n");
                m->mothurOut("The processors parameter allows you to specify the number of processors to use. The default is 1.\n");
@@ -327,20 +341,33 @@ int ClassifySeqsCommand::execute(){
                                                        
                                int outMode=MPI_MODE_CREATE|MPI_MODE_WRONLY; 
                                int inMode=MPI_MODE_RDONLY; 
-                                                               
-                               char outNewTax[newTaxonomyFile.length()];
+                               
+                               //char* outNewTax = new char[newTaxonomyFile.length()];
+                               //memcpy(outNewTax, newTaxonomyFile.c_str(), newTaxonomyFile.length());
+                               
+                               char outNewTax[1024];
                                strcpy(outNewTax, newTaxonomyFile.c_str());
+
+                               //char* outTempTax = new char[tempTaxonomyFile.length()];
+                               //memcpy(outTempTax, tempTaxonomyFile.c_str(), tempTaxonomyFile.length());
                                
-                               char outTempTax[tempTaxonomyFile.length()];
+                               char outTempTax[1024];
                                strcpy(outTempTax, tempTaxonomyFile.c_str());
+
+                               //char* inFileName = new char[fastaFileNames[s].length()];
+                               //memcpy(inFileName, fastaFileNames[s].c_str(), fastaFileNames[s].length());
                                
-                               char inFileName[fastaFileNames[s].length()];
+                               char inFileName[1024];
                                strcpy(inFileName, fastaFileNames[s].c_str());
 
                                MPI_File_open(MPI_COMM_WORLD, inFileName, inMode, MPI_INFO_NULL, &inMPI);  //comm, filename, mode, info, filepointer
                                MPI_File_open(MPI_COMM_WORLD, outNewTax, outMode, MPI_INFO_NULL, &outMPINewTax);
                                MPI_File_open(MPI_COMM_WORLD, outTempTax, outMode, MPI_INFO_NULL, &outMPITempTax);
                                
+                               //delete outNewTax;
+                               //delete outTempTax;
+                               //delete inFileName;
+
                                if (m->control_pressed) {  MPI_File_close(&inMPI);  MPI_File_close(&outMPINewTax);   MPI_File_close(&outMPITempTax);  delete classify; return 0;  }
 
                                if(namefile != "") {  MPIReadNamesFile(namefileNames[s]);  }
@@ -474,50 +501,44 @@ int ClassifySeqsCommand::execute(){
                        if (pid == 0) {  //this part does not need to be paralellized
                #endif
 
-                       //make taxonomy tree from new taxonomy file 
-                       PhyloTree taxaBrowser;
+                       m->mothurOutEndLine();
+                       m->mothurOut("It took " + toString(time(NULL) - start) + " secs to classify " + toString(numFastaSeqs) + " sequences."); m->mothurOutEndLine(); m->mothurOutEndLine();
+                       start = time(NULL);
+                       
+                       PhyloSummary taxaSum(taxonomyFileName, groupfile);
                        
                        if (m->control_pressed) {  for (int i = 0; i < outputNames.size(); i++) {       remove(outputNames[i].c_str()); } delete classify; return 0; }
                        
-                       ifstream in;
-                       openInputFile(tempTaxonomyFile, in);
-               
-                       //read in users taxonomy file and add sequences to tree
-                       string name, taxon;
-                       while(!in.eof()){
-                               in >> name >> taxon; gobble(in);
-                               
-                               if (m->control_pressed) {  for (int i = 0; i < outputNames.size(); i++) {       remove(outputNames[i].c_str()); } remove(tempTaxonomyFile.c_str()); delete classify; return 0; }
+                       if (namefile == "") {  taxaSum.summarize(tempTaxonomyFile);  }
+                       else {
+                               ifstream in;
+                               openInputFile(tempTaxonomyFile, in);
                                
-                               if (namefile != "") {
+                               //read in users taxonomy file and add sequences to tree
+                               string name, taxon;
+                               while(!in.eof()){
+                                       in >> name >> taxon; gobble(in);
+                                       
                                        itNames = nameMap.find(name);
                
                                        if (itNames == nameMap.end()) { 
                                                m->mothurOut(name + " is not in your name file please correct."); m->mothurOutEndLine(); exit(1);
                                        }else{
                                                for (int i = 0; i < itNames->second; i++) { 
-                                                       taxaBrowser.addSeqToTree(name+toString(i), taxon);  //add it as many times as there are identical seqs
+                                                       taxaSum.addSeqToTree(name+toString(i), taxon);  //add it as many times as there are identical seqs
                                                }
                                        }
-                               }else {  taxaBrowser.addSeqToTree(name, taxon);  } //add it once
+                               }
+                               in.close();
                        }
-                       in.close();
-       
-                       taxaBrowser.assignHeirarchyIDs(0);
-                       
-                       if (m->control_pressed) {  for (int i = 0; i < outputNames.size(); i++) {       remove(outputNames[i].c_str()); } remove(tempTaxonomyFile.c_str()); delete classify; return 0; }
-
-                       taxaBrowser.binUnclassified();
-                       
                        remove(tempTaxonomyFile.c_str());
                        
                        if (m->control_pressed) {  for (int i = 0; i < outputNames.size(); i++) {       remove(outputNames[i].c_str()); } delete classify; return 0; }
-
                        
                        //print summary file
                        ofstream outTaxTree;
                        openOutputFile(taxSummary, outTaxTree);
-                       taxaBrowser.print(outTaxTree);
+                       taxaSum.print(outTaxTree);
                        outTaxTree.close();
                        
                        //output taxonomy with the unclassified bins added
@@ -529,9 +550,10 @@ int ClassifySeqsCommand::execute(){
                        openOutputFile(unclass, outTax);
                        
                        //get maxLevel from phylotree so you know how many 'unclassified's to add
-                       int maxLevel = taxaBrowser.getMaxLevel();
+                       int maxLevel = taxaSum.getMaxLevel();
                        
                        //read taxfile - this reading and rewriting is done to preserve the confidence scores.
+                       string name, taxon;
                        while (!inTax.eof()) {
                                if (m->control_pressed) {  for (int i = 0; i < outputNames.size(); i++) {       remove(outputNames[i].c_str()); } remove(unclass.c_str()); delete classify; return 0; }
 
@@ -547,6 +569,9 @@ int ClassifySeqsCommand::execute(){
                        remove(newTaxonomyFile.c_str());
                        rename(unclass.c_str(), newTaxonomyFile.c_str());
                        
+                       m->mothurOutEndLine();
+                       m->mothurOut("It took " + toString(time(NULL) - start) + " secs to create the summary file for  " + toString(numFastaSeqs) + " sequences."); m->mothurOutEndLine(); m->mothurOutEndLine();
+                       
                        #ifdef USE_MPI  
                                }
                        #endif
@@ -555,10 +580,7 @@ int ClassifySeqsCommand::execute(){
                        m->mothurOut("Output File Names: "); m->mothurOutEndLine();
                        for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
                        m->mothurOutEndLine();
-
                        
-                       m->mothurOutEndLine();
-                       m->mothurOut("It took " + toString(time(NULL) - start) + " secs to classify " + toString(numFastaSeqs) + " sequences."); m->mothurOutEndLine(); m->mothurOutEndLine();
                }
                
                delete classify;
@@ -732,12 +754,13 @@ int ClassifySeqsCommand::driverMPI(int start, int num, MPI_File& inMPI, MPI_File
                
                        //read next sequence
                        int length = MPIPos[start+i+1] - MPIPos[start+i];
-                       char buf4[length];
+                       char* buf4 = new char[length];
                        MPI_File_read_at(inMPI, MPIPos[start+i], buf4, length, MPI_CHAR, &status);
                        
                        string tempBuf = buf4;
                        if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length);  }
                        istringstream iss (tempBuf,istringstream::in);
+                       delete buf4;
 
                        Sequence* candidateSeq = new Sequence(iss);
                        
@@ -753,17 +776,19 @@ int ClassifySeqsCommand::driverMPI(int start, int num, MPI_File& inMPI, MPI_File
                                        }
                                        
                                        int length = outputString.length();
-                                       char buf2[length];
-                                       strcpy(buf2, outputString.c_str()); 
+                                       char* buf2 = new char[length];
+                                       memcpy(buf2, outputString.c_str(), length);
                                
                                        MPI_File_write_shared(newFile, buf2, length, MPI_CHAR, &statusNew);
-                                       
+                                       delete buf2;
+
                                        outputString =  candidateSeq->getName() + "\t" + classify->getSimpleTax() + "\n";
                                        length = outputString.length();
-                                       char buf[length];
-                                       strcpy(buf, outputString.c_str()); 
+                                       char* buf = new char[length];
+                                       memcpy(buf, outputString.c_str(), length);
                                
                                        MPI_File_write_shared(tempFile, buf, length, MPI_CHAR, &statusTemp);
+                                       delete buf;
                                }
                        }                               
                        delete candidateSeq;
@@ -791,19 +816,24 @@ int ClassifySeqsCommand::MPIReadNamesFile(string nameFilename){
                MPI_File inMPI;
                MPI_Offset size;
                MPI_Status status;
+
+               //char* inFileName = new char[nameFilename.length()];
+               //memcpy(inFileName, nameFilename.c_str(), nameFilename.length());
                
-               char inFileName[nameFilename.length()];
+               char inFileName[1024];
                strcpy(inFileName, nameFilename.c_str());
 
                MPI_File_open(MPI_COMM_WORLD, inFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);  
                MPI_File_get_size(inMPI, &size);
+               //delete inFileName;
 
-               char buffer[size];
+               char* buffer = new char[size];
                MPI_File_read(inMPI, buffer, size, MPI_CHAR, &status);
 
                string tempBuf = buffer;
                if (tempBuf.length() > size) { tempBuf = tempBuf.substr(0, size);  }
                istringstream iss (tempBuf,istringstream::in);
+               delete buffer;
                
                string firstCol, secondCol;
                while(!iss.eof()) {