classifyseqscommand.h

   1 #ifndef CLASSIFYSEQSCOMMAND_H
   2 #define CLASSIFYSEQSCOMMAND_H
   3
   4 /*
   5  *  classifyseqscommand.h
   6  *  Mothur
   7  *
   8  *  Created by westcott on 11/2/09.
   9  *  Copyright 2009 Schloss Lab. All rights reserved.
  10  *
  11  */
  12
  13
  14 #include "command.hpp"
  15 #include "classify.h"
  16 #include "referencedb.h"
  17 #include "sequence.hpp"
  18 #include "bayesian.h"
  19 #include "phylotree.h"
  20 #include "phylosummary.h"
  21 #include "knn.h"
  22 #include "kmertree.h"
  23 #include "aligntree.h"
  24
  25
  26 //KNN and Wang methods modeled from algorithms in
  27 //Naı¨ve Bayesian Classiﬁer for Rapid Assignment of rRNA Sequences
  28 //into the New Bacterial Taxonomy􏰎†
  29 //Qiong Wang,1 George M. Garrity,1,2 James M. Tiedje,1,2 and James R. Cole1*
  30 //Center for Microbial Ecology1 and Department of Microbiology and Molecular Genetics,2 Michigan State University,
  31 //East Lansing, Michigan 48824
  32 //Received 10 January 2007/Accepted 18 June 2007
  33
  34
  35
  36 class ClassifySeqsCommand : public Command {
  37
  38 public:
  39         ClassifySeqsCommand(string);
  40         ClassifySeqsCommand();
  41         ~ClassifySeqsCommand();
  42
  43         vector<string> setParameters();
  44         string getCommandName()                 { return "classify.seqs";               }
  45         string getCommandCategory()             { return "Phylotype Analysis";  }
  46
  47         string getHelpString();
  48     string getOutputPattern(string);
  49         string getCitation() { return "Wang Q, Garrity GM, Tiedje JM, Cole JR (2007). Naive Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Appl Environ Microbiol 73: 5261-7. [ for Bayesian classifier ] \nAltschul SF, Madden TL, Schaffer AA, Zhang J, Zhang Z, Miller W, Lipman DJ (1997). Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res 25: 3389-402. [ for BLAST ] \nDeSantis TZ, Hugenholtz P, Larsen N, Rojas M, Brodie EL, Keller K, Huber T, Dalevi D, Hu P, Andersen GL (2006). Greengenes, a chimera-checked 16S rRNA gene database and workbench compatible with ARB. Appl Environ Microbiol 72: 5069-72. [ for kmer ] \nhttp://www.mothur.org/wiki/Classify.seqs"; }
  50         string getDescription()         { return "classify sequences"; }
  51
  52         int execute();
  53         void help() { m->mothurOut(getHelpString()); }
  54
  55
  56
  57 private:
  58         struct linePair {
  59                 unsigned long long start;
  60                 unsigned long long end;
  61                 linePair(unsigned long long i, unsigned long long j) : start(i), end(j) {}
  62         };
  63
  64         vector<int> processIDS;   //processid
  65         vector<linePair*> lines;
  66         vector<string> fastaFileNames;
  67         vector<string> namefileNames;
  68     vector<string> countfileNames;
  69         vector<string> groupfileNames;
  70         vector<string> outputNames;
  71         map<string, vector<string> > nameMap;
  72         map<string,  vector<string> >::iterator itNames;
  73
  74         Classify* classify;
  75         ReferenceDB* rdb;
  76
  77         string fastaFileName, templateFileName, countfile, distanceFileName, namefile, search, method, taxonomyFileName, outputDir, groupfile;
  78         int processors, kmerSize, numWanted, cutoff, iters;
  79         float match, misMatch, gapOpen, gapExtend;
  80         bool abort, probs, save, flip, hasName, hasCount, writeShortcuts;
  81
  82         int driver(linePair*, string, string, string, string);
  83         int createProcesses(string, string, string, string);
  84         string addUnclassifieds(string, int);
  85
  86         int MPIReadNamesFile(string);
  87         #ifdef USE_MPI
  88         int driverMPI(int, int, MPI_File&, MPI_File&, MPI_File&, MPI_File&, vector<unsigned long long>&);
  89         #endif
  90 };
  91
  92 /**************************************************************************************************/
  93 //custom data structure for threads to use.
  94 // This is passed by void pointer so it can be any data type
  95 // that can be passed using a single void pointer (LPVOID).
  96 struct classifyData {
  97         string taxFName;
  98         string tempTFName;
  99         string filename;
 100         string search, taxonomyFileName, templateFileName, method, accnos;
 101         unsigned long long start;
 102         unsigned long long end;
 103         MothurOut* m;
 104         float match, misMatch, gapOpen, gapExtend;
 105         int count, kmerSize, threadID, cutoff, iters, numWanted;
 106         bool probs, flip, writeShortcuts;
 107
 108         classifyData(){}
 109         classifyData(string acc, bool p, string me, string te, string tx, string a, string r, string f, string se, int ks, int i, int numW, MothurOut* mout, unsigned long long st, unsigned long long en, float ma, float misMa, float gapO, float gapE, int cut, int tid, bool fli, bool wsh) {
 110                 accnos = acc;
 111                 taxonomyFileName = tx;
 112                 templateFileName = te;
 113                 taxFName = a;
 114                 tempTFName = r;
 115                 filename = f;
 116                 search = se;
 117                 method = me;
 118                 m = mout;
 119                 start = st;
 120                 end = en;
 121                 match = ma;
 122                 misMatch = misMa;
 123                 gapOpen = gapO;
 124                 gapExtend = gapE;
 125                 kmerSize = ks;
 126                 cutoff = cut;
 127                 iters = i;
 128                 numWanted = numW;
 129                 threadID = tid;
 130                 probs = p;
 131                 count = 0;
 132                 flip = fli;
 133         writeShortcuts = wsh;
 134         }
 135 };
 136
 137 /**************************************************************************************************/
 138 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 139 #else
 140 static DWORD WINAPI MyClassThreadFunction(LPVOID lpParam){
 141         classifyData* pDataArray;
 142         pDataArray = (classifyData*)lpParam;
 143
 144         try {
 145                 ofstream outTax;
 146                 pDataArray->m->openOutputFile(pDataArray->taxFName, outTax);
 147
 148                 ofstream outTaxSimple;
 149                 pDataArray->m->openOutputFile(pDataArray->tempTFName, outTaxSimple);
 150
 151                 ofstream outAcc;
 152                 pDataArray->m->openOutputFile(pDataArray->accnos, outAcc);
 153
 154                 ifstream inFASTA;
 155                 pDataArray->m->openInputFile(pDataArray->filename, inFASTA);
 156
 157                 string taxonomy;
 158
 159                 //print header if you are process 0
 160                 if ((pDataArray->start == 0) || (pDataArray->start == 1)) {
 161                         inFASTA.seekg(0);
 162                 }else { //this accounts for the difference in line endings.
 163                         inFASTA.seekg(pDataArray->start-1); pDataArray->m->gobble(inFASTA);
 164                 }
 165
 166                 //make classify
 167                 Classify* myclassify;
 168         string outputMethodTag = pDataArray->method + ".";
 169                 if(pDataArray->method == "bayesian"){   myclassify = new Bayesian(pDataArray->taxonomyFileName, pDataArray->templateFileName, pDataArray->search, pDataArray->kmerSize, pDataArray->cutoff, pDataArray->iters, pDataArray->threadID, pDataArray->flip, pDataArray->writeShortcuts);             }
 170                 else if(pDataArray->method == "knn"){   myclassify = new Knn(pDataArray->taxonomyFileName, pDataArray->templateFileName, pDataArray->search, pDataArray->kmerSize, pDataArray->gapOpen, pDataArray->gapExtend, pDataArray->match, pDataArray->misMatch, pDataArray->numWanted, pDataArray->threadID);                           }
 171         else if(pDataArray->method == "zap"){
 172             outputMethodTag = pDataArray->search + "_" + outputMethodTag;
 173             if (pDataArray->search == "kmer") {   myclassify = new KmerTree(pDataArray->templateFileName, pDataArray->taxonomyFileName, pDataArray->kmerSize, pDataArray->cutoff); }
 174             else {  myclassify = new AlignTree(pDataArray->templateFileName, pDataArray->taxonomyFileName, pDataArray->cutoff);  }
 175         }
 176                 else {
 177                         pDataArray->m->mothurOut(pDataArray->search + " is not a valid method option. I will run the command using bayesian.");
 178                         pDataArray->m->mothurOutEndLine();
 179                         myclassify = new Bayesian(pDataArray->taxonomyFileName, pDataArray->templateFileName, pDataArray->search, pDataArray->kmerSize, pDataArray->cutoff, pDataArray->iters, pDataArray->threadID, pDataArray->flip, pDataArray->writeShortcuts);
 180                 }
 181
 182                 if (pDataArray->m->control_pressed) { delete myclassify; return 0; }
 183
 184                 pDataArray->count = 0;
 185                 for(int i = 0; i < pDataArray->end; i++){ //end is the number of sequences to process
 186
 187                         if (pDataArray->m->control_pressed) { delete myclassify; return 0; }
 188
 189                         Sequence* candidateSeq = new Sequence(inFASTA); pDataArray->m->gobble(inFASTA);
 190
 191                         if (candidateSeq->getName() != "") {
 192
 193                                 taxonomy = myclassify->getTaxonomy(candidateSeq);
 194
 195                                 if (pDataArray->m->control_pressed) { delete candidateSeq; return 0; }
 196
 197                                 if (taxonomy == "unknown;") { pDataArray->m->mothurOut("[WARNING]: " + candidateSeq->getName() + " could not be classified. You can use the remove.lineage command with taxon=unknown; to remove such sequences."); pDataArray->m->mothurOutEndLine(); }
 198
 199                                 //output confidence scores or not
 200                                 if (pDataArray->probs) {
 201                                         outTax << candidateSeq->getName() << '\t' << taxonomy << endl;
 202                                 }else{
 203                                         outTax << candidateSeq->getName() << '\t' << myclassify->getSimpleTax() << endl;
 204                                 }
 205
 206                                 outTaxSimple << candidateSeq->getName() << '\t' << myclassify->getSimpleTax() << endl;
 207
 208                                 if (myclassify->getFlipped()) { outAcc << candidateSeq->getName() << endl; }
 209
 210                                 pDataArray->count++;
 211                         }
 212                         delete candidateSeq;
 213                         //report progress
 214                         if((pDataArray->count) % 100 == 0){     pDataArray->m->mothurOut("Processing sequence: " + toString(pDataArray->count)); pDataArray->m->mothurOutEndLine();             }
 215
 216                 }
 217                 //report progress
 218                 if((pDataArray->count) % 100 != 0){     pDataArray->m->mothurOut("Processing sequence: " + toString(pDataArray->count)); pDataArray->m->mothurOutEndLine();             }
 219
 220                 delete myclassify;
 221                 inFASTA.close();
 222                 outTax.close();
 223                 outTaxSimple.close();
 224
 225         }
 226         catch(exception& e) {
 227                 pDataArray->m->errorOut(e, "ClassifySeqsCommand", "MyClassThreadFunction");
 228                 exit(1);
 229         }
 230 }
 231 #endif
 232
 233
 234
 235
 236 #endif
 237