chimeraperseuscommand.h

   1 #ifndef CHIMERAPERSEUSCOMMAND_H
   2 #define CHIMERAPERSEUSCOMMAND_H
   3
   4
   5 /*
   6  *  chimeraperseuscommand.h
   7  *  Mothur
   8  *
   9  *  Created by westcott on 10/26/11.
  10  *  Copyright 2011 Schloss Lab. All rights reserved.
  11  *
  12  */
  13
  14
  15
  16 #include "mothur.h"
  17 #include "command.hpp"
  18 #include "sequenceparser.h"
  19 #include "myPerseus.h"
  20
  21 /***********************************************************/
  22 class ChimeraPerseusCommand : public Command {
  23 public:
  24         ChimeraPerseusCommand(string);
  25         ChimeraPerseusCommand();
  26         ~ChimeraPerseusCommand() {}
  27
  28         vector<string> setParameters();
  29         string getCommandName()                 { return "chimera.perseus";             }
  30         string getCommandCategory()             { return "Sequence Processing"; }
  31         string getHelpString();
  32         string getCitation() { return "http://www.mothur.org/wiki/Chimera.perseus\n"; }
  33         string getDescription()         { return "detect chimeric sequences"; }
  34
  35         int execute();
  36         void help() { m->mothurOut(getHelpString()); }
  37
  38 private:
  39         struct linePair {
  40                 int start;
  41                 int end;
  42                 linePair(int i, int j) : start(i), end(j) {}
  43         };
  44
  45         bool abort;
  46         string fastafile, groupfile, outputDir, namefile;
  47         int processors;
  48         double cutoff, alpha, beta;
  49
  50         vector<string> outputNames;
  51         vector<string> fastaFileNames;
  52         vector<string> nameFileNames;
  53         vector<string> groupFileNames;
  54
  55         string getNamesFile(string&);
  56         int driver(string, vector<seqData>&, string, int&);
  57         vector<seqData> readFiles(string, string);
  58         vector<seqData> loadSequences(SequenceParser&, string);
  59         int deconvoluteResults(SequenceParser&, string, string);
  60         int driverGroups(SequenceParser&, string, string, int, int, vector<string>);
  61         int createProcessesGroups(SequenceParser&, string, string, vector<string>, string, string, string);
  62 };
  63
  64 /**************************************************************************************************/
  65 //custom data structure for threads to use.
  66 // This is passed by void pointer so it can be any data type
  67 // that can be passed using a single void pointer (LPVOID).
  68 struct perseusData {
  69         string fastafile;
  70         string namefile;
  71         string groupfile;
  72         string outputFName;
  73         string accnos;
  74         MothurOut* m;
  75         int start;
  76         int end;
  77         int threadID, count, numChimeras;
  78         double alpha, beta, cutoff;
  79         vector<string> groups;
  80
  81         perseusData(){}
  82         perseusData(double a, double b, double c, string o,  string f, string n, string g, string ac, vector<string> gr, MothurOut* mout, int st, int en, int tid) {
  83                 alpha = a;
  84                 beta = b;
  85                 cutoff = c;
  86                 fastafile = f;
  87                 namefile = n;
  88                 groupfile = g;
  89                 outputFName = o;
  90                 accnos = ac;
  91                 m = mout;
  92                 start = st;
  93                 end = en;
  94                 threadID = tid;
  95                 groups = gr;
  96                 count = 0;
  97                 numChimeras = 0;
  98         }
  99 };
 100 /**************************************************************************************************/
 101 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 102 #else
 103 static DWORD WINAPI MyPerseusThreadFunction(LPVOID lpParam){
 104         perseusData* pDataArray;
 105         pDataArray = (perseusData*)lpParam;
 106
 107         try {
 108
 109                 //clears files
 110                 ofstream out, out1, out2;
 111                 pDataArray->m->openOutputFile(pDataArray->outputFName, out); out.close();
 112                 pDataArray->m->openOutputFile(pDataArray->accnos, out1); out1.close();
 113
 114                 //parse fasta and name file by group
 115                 SequenceParser* parser;
 116                 if (pDataArray->namefile != "") { parser = new SequenceParser(pDataArray->groupfile, pDataArray->fastafile, pDataArray->namefile);      }
 117                 else                                                    { parser = new SequenceParser(pDataArray->groupfile, pDataArray->fastafile);                                            }
 118
 119                 int totalSeqs = 0;
 120                 int numChimeras = 0;
 121
 122                 for (int i = pDataArray->start; i < pDataArray->end; i++) {
 123
 124                         int start = time(NULL);  if (pDataArray->m->control_pressed) {  delete parser; pDataArray->m->mothurRemove(pDataArray->outputFName); pDataArray->m->mothurRemove(pDataArray->accnos); return 0; }
 125
 126                         pDataArray->m->mothurOutEndLine(); pDataArray->m->mothurOut("Checking sequences from group " + pDataArray->groups[i] + "...");  pDataArray->m->mothurOutEndLine();
 127
 128                         //vector<seqData> sequences = loadSequences(parser, groups[i]); - same function below
 129                         ////////////////////////////////////////////////////////////////////////////////////////
 130                         vector<Sequence> thisGroupsSeqs = parser->getSeqs(pDataArray->groups[i]);
 131                         map<string, string> nameMap = parser->getNameMap(pDataArray->groups[i]);
 132                         map<string, string>::iterator it;
 133
 134                         vector<seqData> sequences;
 135                         bool error = false;
 136
 137                         for (int j = 0; j < thisGroupsSeqs.size(); j++) {
 138
 139                                 if (pDataArray->m->control_pressed) {  delete parser; pDataArray->m->mothurRemove(pDataArray->outputFName); pDataArray->m->mothurRemove(pDataArray->accnos); return 0; }
 140
 141                                 it = nameMap.find(thisGroupsSeqs[j].getName());
 142                                 if (it == nameMap.end()) { error = true; pDataArray->m->mothurOut("[ERROR]: " + thisGroupsSeqs[j].getName() + " is in your fasta file and not in your namefile, please correct."); pDataArray->m->mothurOutEndLine(); }
 143                                 else {
 144                                         int num = pDataArray->m->getNumNames(it->second);
 145                                         sequences.push_back(seqData(thisGroupsSeqs[j].getName(), thisGroupsSeqs[j].getUnaligned(), num));
 146                                 }
 147                         }
 148
 149                         if (error) { pDataArray->m->control_pressed = true; }
 150
 151                         //sort by frequency
 152                         sort(sequences.rbegin(), sequences.rend());
 153                         ////////////////////////////////////////////////////////////////////////////////////////
 154
 155                         if (pDataArray->m->control_pressed) { delete parser; pDataArray->m->mothurRemove(pDataArray->outputFName); pDataArray->m->mothurRemove(pDataArray->accnos); return 0; }
 156
 157                         //int numSeqs = driver((outputFName + groups[i]), sequences, (accnos+groups[i]), numChimeras); - same function below
 158                         ////////////////////////////////////////////////////////////////////////////////////////
 159                         string chimeraFileName = pDataArray->outputFName+pDataArray->groups[i];
 160                         string accnosFileName = pDataArray->accnos+pDataArray->groups[i];
 161
 162                         vector<vector<double> > correctModel(4);        //could be an option in the future to input own model matrix
 163                         for(int j=0;j<4;j++){   correctModel[j].resize(4);      }
 164
 165                         correctModel[0][0] = 0.000000;  //AA
 166                         correctModel[1][0] = 11.619259; //CA
 167                         correctModel[2][0] = 11.694004; //TA
 168                         correctModel[3][0] = 7.748623;  //GA
 169
 170                         correctModel[1][1] = 0.000000;  //CC
 171                         correctModel[2][1] = 7.619657;  //TC
 172                         correctModel[3][1] = 12.852562; //GC
 173
 174                         correctModel[2][2] = 0.000000;  //TT
 175                         correctModel[3][2] = 10.964048; //TG
 176
 177                         correctModel[3][3] = 0.000000;  //GG
 178
 179                         for(int k=0;k<4;k++){
 180                                 for(int j=0;j<k;j++){
 181                                         correctModel[j][k] = correctModel[k][j];
 182                                 }
 183                         }
 184
 185                         int numSeqs = sequences.size();
 186                         int alignLength = sequences[0].sequence.size();
 187
 188                         ofstream chimeraFile;
 189                         ofstream accnosFile;
 190                         pDataArray->m->openOutputFile(chimeraFileName, chimeraFile);
 191                         pDataArray->m->openOutputFile(accnosFileName, accnosFile);
 192
 193                         Perseus myPerseus;
 194                         vector<vector<double> > binMatrix = myPerseus.binomial(alignLength);
 195
 196                         chimeraFile << "SequenceIndex\tName\tDiffsToBestMatch\tBestMatchIndex\tBestMatchName\tDiffstToChimera\tIndexofLeftParent\tIndexOfRightParent\tNameOfLeftParent\tNameOfRightParent\tDistanceToBestMatch\tcIndex\t(cIndex - singleDist)\tloonIndex\tMismatchesToChimera\tMismatchToTrimera\tChimeraBreakPoint\tLogisticProbability\tTypeOfSequence\n";
 197
 198                         vector<bool> chimeras(numSeqs, 0);
 199
 200                         for(int j=0;j<numSeqs;j++){
 201
 202                                 if (pDataArray->m->control_pressed) { delete parser; pDataArray->m->mothurRemove(pDataArray->outputFName); pDataArray->m->mothurRemove(pDataArray->accnos); chimeraFile.close(); pDataArray->m->mothurRemove(chimeraFileName); accnosFile.close(); pDataArray->m->mothurRemove(accnosFileName); return 0; }
 203
 204                                 vector<bool> restricted = chimeras;
 205
 206                                 vector<vector<int> > leftDiffs(numSeqs);
 207                                 vector<vector<int> > leftMaps(numSeqs);
 208                                 vector<vector<int> > rightDiffs(numSeqs);
 209                                 vector<vector<int> > rightMaps(numSeqs);
 210
 211                                 vector<int> singleLeft, bestLeft;
 212                                 vector<int> singleRight, bestRight;
 213
 214                                 int bestSingleIndex, bestSingleDiff;
 215                                 vector<pwAlign> alignments(numSeqs);
 216
 217                                 int comparisons = myPerseus.getAlignments(j, sequences, alignments, leftDiffs, leftMaps, rightDiffs, rightMaps, bestSingleIndex, bestSingleDiff, restricted);
 218
 219                                 if (pDataArray->m->control_pressed) { delete parser; pDataArray->m->mothurRemove(pDataArray->outputFName); pDataArray->m->mothurRemove(pDataArray->accnos); chimeraFile.close(); pDataArray->m->mothurRemove(chimeraFileName); accnosFile.close(); pDataArray->m->mothurRemove(accnosFileName); return 0; }
 220
 221                                 int minMismatchToChimera, leftParentBi, rightParentBi, breakPointBi;
 222
 223                                 string dummyA, dummyB;
 224
 225                                 if(comparisons >= 2){
 226                                         minMismatchToChimera = myPerseus.getChimera(sequences, leftDiffs, rightDiffs, leftParentBi, rightParentBi, breakPointBi, singleLeft, bestLeft, singleRight, bestRight, restricted);
 227
 228                                         if (pDataArray->m->control_pressed) { delete parser;  pDataArray->m->mothurRemove(pDataArray->outputFName); pDataArray->m->mothurRemove(pDataArray->accnos); chimeraFile.close(); pDataArray->m->mothurRemove(chimeraFileName); accnosFile.close(); pDataArray->m->mothurRemove(accnosFileName); return 0; }
 229
 230                                         int minMismatchToTrimera = numeric_limits<int>::max();
 231                                         int leftParentTri, middleParentTri, rightParentTri, breakPointTriA, breakPointTriB;
 232
 233                                         if(minMismatchToChimera >= 3 && comparisons >= 3){
 234                                                 minMismatchToTrimera = myPerseus.getTrimera(sequences, leftDiffs, leftParentTri, middleParentTri, rightParentTri, breakPointTriA, breakPointTriB, singleLeft, bestLeft, singleRight, bestRight, restricted);
 235
 236                                                 if (pDataArray->m->control_pressed) { delete parser;  pDataArray->m->mothurRemove(pDataArray->outputFName); pDataArray->m->mothurRemove(pDataArray->accnos); chimeraFile.close(); pDataArray->m->mothurRemove(chimeraFileName); accnosFile.close(); pDataArray->m->mothurRemove(accnosFileName); return 0; }
 237                                         }
 238
 239                                         double singleDist = myPerseus.modeledPairwiseAlignSeqs(sequences[j].sequence, sequences[bestSingleIndex].sequence, dummyA, dummyB, correctModel);
 240
 241                                         if (pDataArray->m->control_pressed) { delete parser;  pDataArray->m->mothurRemove(pDataArray->outputFName); pDataArray->m->mothurRemove(pDataArray->accnos); chimeraFile.close(); pDataArray->m->mothurRemove(chimeraFileName); accnosFile.close(); pDataArray->m->mothurRemove(accnosFileName); return 0; }
 242
 243                                         string type;
 244                                         string chimeraRefSeq;
 245
 246                                         if(minMismatchToChimera - minMismatchToTrimera >= 3){
 247                                                 type = "trimera";
 248                                                 chimeraRefSeq = myPerseus.stitchTrimera(alignments, leftParentTri, middleParentTri, rightParentTri, breakPointTriA, breakPointTriB, leftMaps, rightMaps);
 249                                         }
 250                                         else{
 251                                                 type = "chimera";
 252                                                 chimeraRefSeq = myPerseus.stitchBimera(alignments, leftParentBi, rightParentBi, breakPointBi, leftMaps, rightMaps);
 253                                         }
 254
 255                                         if (pDataArray->m->control_pressed) { delete parser; pDataArray->m->mothurRemove(pDataArray->outputFName); pDataArray->m->mothurRemove(pDataArray->accnos); chimeraFile.close(); pDataArray->m->mothurRemove(chimeraFileName); accnosFile.close(); pDataArray->m->mothurRemove(accnosFileName); return 0; }
 256
 257                                         double chimeraDist = myPerseus.modeledPairwiseAlignSeqs(sequences[j].sequence, chimeraRefSeq, dummyA, dummyB, correctModel);
 258
 259                                         if (pDataArray->m->control_pressed) { delete parser; pDataArray->m->mothurRemove(pDataArray->outputFName); pDataArray->m->mothurRemove(pDataArray->accnos); chimeraFile.close(); pDataArray->m->mothurRemove(chimeraFileName); accnosFile.close(); pDataArray->m->mothurRemove(accnosFileName); return 0; }
 260
 261                                         double cIndex = chimeraDist;//modeledPairwiseAlignSeqs(sequences[j].sequence, chimeraRefSeq);
 262                                         double loonIndex = myPerseus.calcLoonIndex(sequences[j].sequence, sequences[leftParentBi].sequence, sequences[rightParentBi].sequence, breakPointBi, binMatrix);
 263
 264                                         if (pDataArray->m->control_pressed) { delete parser; pDataArray->m->mothurRemove(pDataArray->outputFName); pDataArray->m->mothurRemove(pDataArray->accnos); chimeraFile.close(); pDataArray->m->mothurRemove(chimeraFileName); accnosFile.close(); pDataArray->m->mothurRemove(accnosFileName); return 0; }
 265
 266                                         chimeraFile << j << '\t' << sequences[j].seqName << '\t' << bestSingleDiff << '\t' << bestSingleIndex << '\t' << sequences[bestSingleIndex].seqName << '\t';
 267                                         chimeraFile << minMismatchToChimera << '\t' << leftParentBi << '\t' << rightParentBi << '\t' << sequences[leftParentBi].seqName << '\t' << sequences[rightParentBi].seqName << '\t';
 268                                         chimeraFile << singleDist << '\t' << cIndex << '\t' << (cIndex - singleDist) << '\t' << loonIndex << '\t';
 269                                         chimeraFile << minMismatchToChimera << '\t' << minMismatchToTrimera << '\t' << breakPointBi << '\t';
 270
 271                                         double probability = myPerseus.classifyChimera(singleDist, cIndex, loonIndex, pDataArray->alpha, pDataArray->beta);
 272
 273                                         chimeraFile << probability << '\t';
 274
 275                                         if(probability > pDataArray->cutoff){
 276                                                 chimeraFile << type << endl;
 277                                                 accnosFile << sequences[j].seqName << endl;
 278                                                 chimeras[j] = 1;
 279                                                 numChimeras++;
 280                                         }
 281                                         else{
 282                                                 chimeraFile << "good" << endl;
 283                                         }
 284
 285                                 }
 286                                 else{
 287                                         chimeraFile << j << '\t' << sequences[j].seqName << "\t0\t0\tNull\t0\t0\t0\tNull\tNull\t0.0\t0.0\t0.0\t0\t0\t0\t0.0\t0.0\tgood" << endl;
 288                                 }
 289                                 //report progress
 290                                 if((j+1) % 100 == 0){   pDataArray->m->mothurOut("Processing sequence: " + toString(j+1) + "\n");               }
 291                         }
 292
 293                         if((numSeqs) % 100 != 0){       pDataArray->m->mothurOut("Processing sequence: " + toString(numSeqs) + "\n");           }
 294
 295                         chimeraFile.close();
 296                         accnosFile.close();
 297                         ////////////////////////////////////////////////////////////////////////////////////////
 298
 299                         totalSeqs += numSeqs;
 300
 301                         //append files
 302                         pDataArray->m->appendFiles(chimeraFileName, pDataArray->outputFName); pDataArray->m->mothurRemove(chimeraFileName);
 303                         pDataArray->m->appendFiles(accnosFileName, pDataArray->accnos); pDataArray->m->mothurRemove(accnosFileName);
 304                         pDataArray->m->mothurOutEndLine(); pDataArray->m->mothurOut("It took " + toString(time(NULL) - start) + " secs to check " + toString(numSeqs) + " sequences from group " + pDataArray->groups[i] + ".");        pDataArray->m->mothurOutEndLine();
 305
 306                         if (pDataArray->m->control_pressed) { delete parser; pDataArray->m->mothurRemove(pDataArray->outputFName); pDataArray->m->mothurRemove(pDataArray->accnos); return 0; }
 307                 }
 308
 309                 pDataArray->count = totalSeqs;
 310                 delete parser;
 311                 return totalSeqs;
 312
 313         }
 314         catch(exception& e) {
 315                 pDataArray->m->errorOut(e, "ChimeraUchimeCommand", "MyPerseusThreadFunction");
 316                 exit(1);
 317         }
 318 }
 319 /**************************************************************************************************/
 320
 321 #endif
 322
 323 #endif
 324
 325