getmetacommunitycommand.h

   1 //
   2 //  getmetacommunitycommand.h
   3 //  Mothur
   4 //
   5 //  Created by SarahsWork on 4/9/13.
   6 //  Copyright (c) 2013 Schloss Lab. All rights reserved.
   7 //
   8
   9 #ifndef Mothur_getmetacommunitycommand_h
  10 #define Mothur_getmetacommunitycommand_h
  11
  12 #include "command.hpp"
  13 #include "inputdata.h"
  14 #include "qFinderDMM.h"
  15
  16 /**************************************************************************************************/
  17
  18 class GetMetaCommunityCommand : public Command {
  19 public:
  20     GetMetaCommunityCommand(string);
  21     GetMetaCommunityCommand();
  22     ~GetMetaCommunityCommand(){}
  23
  24     vector<string> setParameters();
  25     string getCommandName()                     { return "get.communitytype";           }
  26     string getCommandCategory()         { return "OTU-Based Approaches";         }
  27
  28     string getOutputPattern(string);
  29
  30         string getHelpString();
  31     string getCitation() { return "Holmes I, Harris K, Quince C (2012) Dirichlet Multinomial Mixtures: Generative Models for Microbial Metagenomics. PLoS ONE 7(2): e30126. doi:10.1371/journal.pone.0030126 http://www.mothur.org/wiki/get.communitytype"; }
  32     string getDescription()             { return "Assigns samples to bins using a Dirichlet multinomial mixture model"; }
  33
  34     int execute();
  35     void help() { m->mothurOut(getHelpString()); }
  36
  37 private:
  38     struct linePair {
  39                 unsigned long long start;
  40                 unsigned long long end;
  41                 linePair(unsigned long long i, unsigned long long j) : start(i), end(j) {}
  42         };
  43     bool abort, allLines;
  44     string outputDir;
  45     vector<string> outputNames;
  46     string sharedfile;
  47     int minpartitions, maxpartitions, optimizegap, processors;
  48     vector<string> Groups;
  49     set<string> labels;
  50
  51     int processDriver(vector<SharedRAbundVector*>&, vector<int>&, string, vector<string>, vector<string>, vector<string>, int);
  52     int createProcesses(vector<SharedRAbundVector*>&);
  53     vector<double> generateDesignFile(int, map<string,string>);
  54     int generateSummaryFile(int, map<string,string>);
  55
  56 };
  57
  58 /**************************************************************************************************/
  59 struct summaryData {
  60
  61     string name;
  62     double refMean, difference;
  63     vector<double> partMean, partLCI, partUCI;
  64
  65 };
  66 /**************************************************************************************************/
  67
  68 struct metaCommunityData {
  69     vector<SharedRAbundVector*> thislookup;
  70         MothurOut* m;
  71         string outputFileName;
  72     vector<string> relabunds, matrix, outputNames;
  73     int minpartitions, maxpartitions, optimizegap;
  74     vector<int> parts;
  75     int minPartition;
  76
  77         metaCommunityData(){}
  78         metaCommunityData(vector<SharedRAbundVector*> lu, MothurOut* mout, vector<int> dp, string fit, vector<string> rels, vector<string> mat, int minp, int maxp, int opg) {
  79                 m = mout;
  80         thislookup = lu;
  81         parts = dp;
  82         outputFileName = fit;
  83         relabunds = rels;
  84         matrix = mat;
  85         minpartitions = minp;
  86         maxpartitions = maxp;
  87         optimizegap = opg;
  88         minPartition = 0;
  89         }
  90 };
  91 /**************************************************************************************************/
  92 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
  93 #else
  94 static DWORD WINAPI MyMetaCommunityThreadFunction(LPVOID lpParam){
  95         metaCommunityData* pDataArray;
  96         pDataArray = (metaCommunityData*)lpParam;
  97
  98         try {
  99
 100         double minLaplace = 1e10;
 101
 102                 ofstream fitData;
 103                 pDataArray->m->openOutputFile(pDataArray->outputFileName, fitData);
 104         fitData.setf(ios::fixed, ios::floatfield);
 105         fitData.setf(ios::showpoint);
 106         cout.setf(ios::fixed, ios::floatfield);
 107         cout.setf(ios::showpoint);
 108
 109         vector< vector<int> > sharedMatrix;
 110         for (int i = 0; i < pDataArray->thislookup.size(); i++) { sharedMatrix.push_back(pDataArray->thislookup[i]->getAbundances()); }
 111
 112         pDataArray->m->mothurOut("K\tNLE\t\tlogDet\tBIC\t\tAIC\t\tLaplace\n");
 113         fitData << "K\tNLE\tlogDet\tBIC\tAIC\tLaplace" << endl;
 114
 115         for(int i=0;i<pDataArray->parts.size();i++){
 116
 117             int numPartitions = pDataArray->parts[i];
 118
 119             if (pDataArray->m->debug) { pDataArray->m->mothurOut("[DEBUG]: running partition " + toString(numPartitions) + "\n"); }
 120
 121             if (pDataArray->m->control_pressed) { break; }
 122
 123             qFinderDMM* findQ = new qFinderDMM(sharedMatrix, numPartitions);
 124
 125             if (pDataArray->m->debug) { pDataArray->m->mothurOut("[DEBUG]: done finding Q " + toString(numPartitions) + "\n"); }
 126
 127             double laplace = findQ->getLaplace();
 128             pDataArray->m->mothurOut(toString(numPartitions) + '\t');
 129             cout << setprecision (2) << findQ->getNLL() << '\t' << findQ->getLogDet() << '\t';
 130             pDataArray->m->mothurOutJustToLog(toString(findQ->getNLL()) + '\t' + toString(findQ->getLogDet()) + '\t');
 131             cout << findQ->getBIC() << '\t' << findQ->getAIC() << '\t' << laplace;
 132             pDataArray->m->mothurOutJustToLog(toString(findQ->getBIC()) + '\t' + toString(findQ->getAIC()) + '\t' + toString(laplace));
 133
 134             fitData << numPartitions << '\t';
 135             fitData << setprecision (2) << findQ->getNLL() << '\t' << findQ->getLogDet() << '\t';
 136             fitData << findQ->getBIC() << '\t' << findQ->getAIC() << '\t' << laplace << endl;
 137
 138             if(laplace < minLaplace){
 139                 pDataArray->minPartition = numPartitions;
 140                 minLaplace = laplace;
 141                 pDataArray->m->mothurOut("***");
 142             }
 143             pDataArray->m->mothurOutEndLine();
 144
 145             pDataArray->outputNames.push_back(pDataArray->relabunds[i]);
 146             pDataArray->outputNames.push_back(pDataArray->matrix[i]);
 147
 148             findQ->printZMatrix(pDataArray->matrix[i], pDataArray->m->getGroups());
 149             findQ->printRelAbund(pDataArray->relabunds[i], pDataArray->m->currentBinLabels);
 150
 151             if(pDataArray->optimizegap != -1 && (numPartitions - pDataArray->minPartition) >= pDataArray->optimizegap && numPartitions >= pDataArray->minpartitions){ break; }
 152
 153             delete findQ;
 154         }
 155         fitData.close();
 156
 157         //minPartition = 4;
 158
 159         if (pDataArray->m->control_pressed) { return 0; }
 160
 161         return 0;
 162
 163         }
 164         catch(exception& e) {
 165                 pDataArray->m->errorOut(e, "GetMetaCommunityCommand", "MyMetaCommunityThreadFunction");
 166                 exit(1);
 167         }
 168 }
 169 #endif
 170
 171
 172
 173 #endif