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