shhhseqscommand.h

   1 #ifndef SHHHSEQSCOMMAND_H
   2 #define SHHHSEQSCOMMAND_H
   3
   4 /*
   5  *  shhhseqscommand.h
   6  *  Mothur
   7  *
   8  *  Created by westcott on 11/8/11.
   9  *  Copyright 2011 Schloss Lab. All rights reserved.
  10  *
  11  */
  12
  13 #include "command.hpp"
  14 #include "myseqdist.h"
  15 #include "seqnoise.h"
  16 #include "sequenceparser.h"
  17 #include "deconvolutecommand.h"
  18 #include "clustercommand.h"
  19
  20 //**********************************************************************************************************************
  21
  22 class ShhhSeqsCommand : public Command {
  23
  24 public:
  25         ShhhSeqsCommand(string);
  26         ShhhSeqsCommand();
  27         ~ShhhSeqsCommand() {}
  28
  29         vector<string> setParameters();
  30         string getCommandName()                 { return "shhh.seqs";   }
  31         string getCommandCategory()             { return "Sequence Processing";         }
  32         string getOutputFileNameTag(string, string);
  33         string getHelpString();
  34         string getCitation() { return "Schloss PD, Gevers D, Westcott SL (2011).  Reducing the effects of PCR amplification and sequencing artifacts on 16S rRNA-based studies.  PLoS ONE.  6:e27310.\nQuince C, Lanzen A, Davenport RJ, Turnbaugh PJ (2011).  Removing noise from pyrosequenced amplicons.  BMC Bioinformatics  12:38.\nhttp://www.mothur.org/wiki/Shhh.seqs"; }
  35         string getDescription()         { return "shhh.seqs"; }
  36
  37
  38         int execute();
  39         void help() { m->mothurOut(getHelpString()); }
  40
  41 private:
  42
  43         struct linePair {
  44                 int start;
  45                 int end;
  46                 linePair(int i, int j) : start(i), end(j) {}
  47         };
  48
  49         bool abort;
  50         string outputDir, fastafile, namefile, groupfile;
  51         int processors;
  52         double sigma;
  53         vector<string> outputNames;
  54
  55         int readData(correctDist*, seqNoise&, vector<string>&, vector<string>&, vector<string>&, vector<int>&);
  56         int loadData(correctDist*, seqNoise&, vector<string>&, vector<string>&, vector<string>&, vector<int>&, map<string, string>&, vector<Sequence>&);
  57
  58         int driver(seqNoise&, vector<string>&, vector<string>&, vector<string>&, vector<int>&, string, string, string, string);
  59         vector<string> driverGroups(SequenceParser&, string, string, string, int, int, vector<string>);
  60         vector<string> createProcessesGroups(SequenceParser&, string, string, string, vector<string>);
  61         int deconvoluteResults(string, string);
  62
  63
  64
  65 };
  66
  67 /**************************************************************************************************/
  68 //custom data structure for threads to use.
  69 // This is passed by void pointer so it can be any data type
  70 // that can be passed using a single void pointer (LPVOID).
  71 struct shhhseqsData {
  72         string fastafile;
  73         string namefile;
  74         string groupfile;
  75         string newFName, newNName, newMName;
  76         MothurOut* m;
  77         int start;
  78         int end;
  79         int sigma, threadID;
  80         vector<string> groups;
  81         vector<string> mapfileNames;
  82
  83         shhhseqsData(){}
  84         shhhseqsData(string f, string n, string g, string nff,  string nnf, string nmf, vector<string> gr, MothurOut* mout, int st, int en, int s, int tid) {
  85                 fastafile = f;
  86                 namefile = n;
  87                 groupfile = g;
  88                 newFName = nff;
  89                 newNName = nnf;
  90                 newMName = nmf;
  91                 m = mout;
  92                 start = st;
  93                 end = en;
  94                 sigma = s;
  95                 threadID = tid;
  96                 groups = gr;
  97         }
  98 };
  99
 100 /**************************************************************************************************/
 101 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 102 #else
 103 static DWORD WINAPI MyShhhSeqsThreadFunction(LPVOID lpParam){
 104         shhhseqsData* pDataArray;
 105         pDataArray = (shhhseqsData*)lpParam;
 106
 107         try {
 108
 109                 //parse fasta and name file by group
 110                 SequenceParser parser(pDataArray->groupfile, pDataArray->fastafile, pDataArray->namefile);
 111
 112                 //precluster each group
 113                 for (int k = pDataArray->start; k < pDataArray->end; k++) {
 114
 115                         int start = time(NULL);
 116
 117                         if (pDataArray->m->control_pressed) {  return 0; }
 118
 119                         pDataArray->m->mothurOutEndLine(); pDataArray->m->mothurOut("Processing group " + pDataArray->groups[k] + ":"); pDataArray->m->mothurOutEndLine();
 120
 121                         map<string, string> thisNameMap;
 122                         thisNameMap = parser.getNameMap(pDataArray->groups[k]);
 123                         vector<Sequence> thisSeqs = parser.getSeqs(pDataArray->groups[k]);
 124
 125                         if (pDataArray->m->control_pressed) {  return 0; }
 126
 127                         vector<string> sequences;
 128                         vector<string> uniqueNames;
 129                         vector<string> redundantNames;
 130                         vector<int> seqFreq;
 131
 132                         seqNoise noise;
 133                         correctDist* correct = new correctDist(1); //we use one processor since we already split up the work load.
 134
 135                         //load this groups info in order
 136                         //loadData(correct, noise, sequences, uniqueNames, redundantNames, seqFreq, thisNameMap, thisSeqs);
 137                         //////////////////////////////////////////////////////////////////////////////////////////////////
 138                         bool error = false;
 139                         map<string, string>::iterator it;
 140
 141                         for (int i = 0; i < thisSeqs.size(); i++) {
 142
 143                                 if (pDataArray->m->control_pressed) { return 0; }
 144
 145                                 if (thisSeqs[i].getName() != "") {
 146                                         correct->addSeq(thisSeqs[i].getName(), thisSeqs[i].getAligned());
 147
 148                                         it = thisNameMap.find(thisSeqs[i].getName());
 149                                         if (it != thisNameMap.end()) {
 150                                                 noise.addSeq(thisSeqs[i].getAligned(), sequences);
 151                                                 noise.addRedundantName(it->first, it->second, uniqueNames, redundantNames, seqFreq);
 152                                         }else {
 153                                                 pDataArray->m->mothurOut("[ERROR]: " + thisSeqs[i].getName() + " is in your fasta file and not in your namefile, please correct.");
 154                                                 error = true;
 155                                         }
 156                                 }
 157                         }
 158
 159                         if (error) { return 0; }
 160                         //////////////////////////////////////////////////////////////////////////////////////////////////
 161
 162                         if (pDataArray->m->control_pressed) { return 0; }
 163
 164                         //calc distances for cluster
 165                         string distFileName = pDataArray->m->getRootName(pDataArray->m->getSimpleName(pDataArray->fastafile)) + pDataArray->groups[k] + ".shhh.dist";
 166                         correct->execute(distFileName);
 167                         delete correct;
 168
 169                         if (pDataArray->m->control_pressed) { pDataArray->m->mothurRemove(distFileName); return 0; }
 170
 171                         //driver(noise, sequences, uniqueNames, redundantNames, seqFreq, distFileName, newFFile+groups[i], newNFile+groups[i], newMFile+groups[i]+".map");
 172                         ///////////////////////////////////////////////////////////////////////////////////////////////////
 173                         double cutOff = 0.08;
 174                         int minIter = 10;
 175                         int maxIter = 1000;
 176                         double minDelta = 1e-6;
 177                         int numIters = 0;
 178                         double maxDelta = 1e6;
 179                         int numSeqs = sequences.size();
 180
 181                         //run cluster command
 182                         string inputString = "phylip=" + distFileName + ", method=furthest, cutoff=0.08";
 183                         pDataArray->m->mothurOut("/******************************************/"); pDataArray->m->mothurOutEndLine();
 184                         pDataArray->m->mothurOut("Running command: cluster(" + inputString + ")"); pDataArray->m->mothurOutEndLine();
 185
 186                         Command* clusterCommand = new ClusterCommand(inputString);
 187                         clusterCommand->execute();
 188
 189                         map<string, vector<string> > filenames = clusterCommand->getOutputFiles();
 190                         string listFileName = filenames["list"][0];
 191                         string rabundFileName = filenames["rabund"][0]; pDataArray->m->mothurRemove(rabundFileName);
 192                         string sabundFileName = filenames["sabund"][0]; pDataArray->m->mothurRemove(sabundFileName);
 193
 194                         delete clusterCommand;
 195                         pDataArray->m->mothurOut("/******************************************/"); pDataArray->m->mothurOutEndLine();
 196
 197                         if (pDataArray->m->control_pressed) { pDataArray->m->mothurRemove(distFileName); pDataArray->m->mothurRemove(listFileName); return 0; }
 198
 199                         vector<double> distances(numSeqs * numSeqs);
 200                         noise.getDistanceData(distFileName, distances);
 201                         pDataArray->m->mothurRemove(distFileName);
 202                         if (pDataArray->m->control_pressed) { pDataArray->m->mothurRemove(listFileName); return 0; }
 203
 204                         vector<int> otuData(numSeqs);
 205                         vector<int> otuFreq;
 206                         vector<vector<int> > otuBySeqLookUp;
 207                         noise.getListData(listFileName, cutOff, otuData, otuFreq, otuBySeqLookUp);
 208                         pDataArray->m->mothurRemove(listFileName);
 209                         if (pDataArray->m->control_pressed) { return 0; }
 210
 211                         int numOTUs = otuFreq.size();
 212
 213                         vector<double> weights(numOTUs, 0);
 214                         vector<int> change(numOTUs, 1);
 215                         vector<int> centroids(numOTUs, -1);
 216                         vector<int> cumCount(numOTUs, 0);
 217
 218                         vector<double> tau(numSeqs, 1);
 219                         vector<int> anP(numSeqs, 0);
 220                         vector<int> anI(numSeqs, 0);
 221                         vector<int> anN(numSeqs, 0);
 222                         vector<vector<int> > aanI = otuBySeqLookUp;
 223
 224                         while(numIters < minIter || ((maxDelta > minDelta) && (numIters < maxIter))){
 225
 226                                 if (pDataArray->m->control_pressed) { return 0; }
 227
 228                                 noise.updateOTUCountData(otuFreq, otuBySeqLookUp, aanI, anP, anI, cumCount); if (pDataArray->m->control_pressed) { return 0; }
 229                                 maxDelta = noise.calcNewWeights(weights, seqFreq, anI, cumCount, anP, otuFreq, tau);  if (pDataArray->m->control_pressed) { return 0; }
 230
 231                                 noise.calcCentroids(anI, anP, change, centroids, cumCount, distances, seqFreq, otuFreq, tau); if (pDataArray->m->control_pressed) { return 0; }
 232                                 noise.checkCentroids(weights, centroids); if (pDataArray->m->control_pressed) { return 0; }
 233
 234                                 otuFreq.assign(numOTUs, 0);
 235
 236                                 int total = 0;
 237
 238                                 for(int i=0;i<numSeqs;i++){
 239                                         if (pDataArray->m->control_pressed) { return 0; }
 240
 241                                         double offset = 1e6;
 242                                         double norm = 0.0000;
 243                                         double minWeight = 0.1;
 244                                         vector<double> currentTau(numOTUs);
 245
 246                                         for(int j=0;j<numOTUs;j++){
 247                                                 if (pDataArray->m->control_pressed) { return 0; }
 248                                                 if(weights[j] > minWeight && distances[i * numSeqs+centroids[j]] < offset){
 249                                                         offset = distances[i * numSeqs+centroids[j]];
 250                                                 }
 251                                         }
 252
 253                                         for(int j=0;j<numOTUs;j++){
 254                                                 if (pDataArray->m->control_pressed) { return 0; }
 255                                                 if(weights[j] > minWeight){
 256                                                         currentTau[j] = exp(pDataArray->sigma * (-distances[(i * numSeqs + centroids[j])] + offset)) * weights[j];
 257                                                         norm += currentTau[j];
 258                                                 }
 259                                                 else{
 260                                                         currentTau[j] = 0.0000;
 261                                                 }
 262                                         }
 263
 264                                         for(int j=0;j<numOTUs;j++){
 265                                                 if (pDataArray->m->control_pressed) { return 0; }
 266                                                 currentTau[j] /= norm;
 267                                         }
 268
 269                                         for(int j=0;j<numOTUs;j++){
 270                                                 if (pDataArray->m->control_pressed) { return 0; }
 271
 272                                                 if(currentTau[j] > 1.0e-4){
 273                                                         int oldTotal = total;
 274                                                         total++;
 275
 276                                                         tau.resize(oldTotal+1);
 277                                                         tau[oldTotal] = currentTau[j];
 278                                                         otuBySeqLookUp[j][otuFreq[j]] = oldTotal;
 279                                                         aanI[j][otuFreq[j]] = i;
 280                                                         otuFreq[j]++;
 281
 282                                                 }
 283                                         }
 284
 285                                         anP.resize(total);
 286                                         anI.resize(total);
 287                                 }
 288
 289                                 numIters++;
 290                         }
 291
 292                         noise.updateOTUCountData(otuFreq, otuBySeqLookUp, aanI, anP, anI, cumCount);  if (pDataArray->m->control_pressed) { return 0; }
 293
 294                         vector<double> percentage(numSeqs);
 295                         noise.setUpOTUData(otuData, percentage, cumCount, tau, otuFreq, anP, anI);  if (pDataArray->m->control_pressed) { return 0; }
 296                         noise.finishOTUData(otuData, otuFreq, anP, anI, cumCount, otuBySeqLookUp, aanI, tau);  if (pDataArray->m->control_pressed) { return 0; }
 297
 298                         change.assign(numOTUs, 1);
 299                         noise.calcCentroids(anI, anP, change, centroids, cumCount, distances, seqFreq, otuFreq, tau); if (pDataArray->m->control_pressed) { return 0; }
 300
 301
 302                         vector<int> finalTau(numOTUs, 0);
 303                         for(int i=0;i<numSeqs;i++){
 304                                 if (pDataArray->m->control_pressed) { return 0; }
 305                                 finalTau[otuData[i]] += int(seqFreq[i]);
 306                         }
 307
 308                         noise.writeOutput(pDataArray->newFName+pDataArray->groups[k], pDataArray->newNName+pDataArray->groups[k], pDataArray->newMName+pDataArray->groups[k]+".map", finalTau, centroids, otuData, sequences, uniqueNames, redundantNames, seqFreq, distances);
 309
 310                         ///////////////////////////////////////////////////////////////////////////////////////////////////
 311
 312                         if (pDataArray->m->control_pressed) { return 0; }
 313
 314                         pDataArray->m->appendFiles(pDataArray->newFName+pDataArray->groups[k], pDataArray->newFName); pDataArray->m->mothurRemove(pDataArray->newFName+pDataArray->groups[k]);
 315                         pDataArray->m->appendFiles(pDataArray->newNName+pDataArray->groups[k], pDataArray->newNName); pDataArray->m->mothurRemove(pDataArray->newNName+pDataArray->groups[k]);
 316                         pDataArray->mapfileNames.push_back(pDataArray->newMName+pDataArray->groups[k]+".map");
 317
 318                         pDataArray->m->mothurOut("It took " + toString(time(NULL) - start) + " secs to process group " + pDataArray->groups[k] + "."); pDataArray->m->mothurOutEndLine();
 319                 }
 320
 321                 return 0;
 322
 323
 324         }
 325         catch(exception& e) {
 326                 pDataArray->m->errorOut(e, "ShhhSeqsCommand", "MyShhhSeqsThreadFunction");
 327                 exit(1);
 328         }
 329 }
 330 #endif
 331 /**************************************************************************************************/
 332
 333 #endif