makecontigscommand.h

   1 #ifndef Mothur_makecontigscommand_h
   2 #define Mothur_makecontigscommand_h
   3
   4 //
   5 //  makecontigscommand.h
   6 //  Mothur
   7 //
   8 //  Created by Sarah Westcott on 5/15/12.
   9 //  Copyright (c) 2012 Schloss Lab. All rights reserved.
  10 //
  11
  12 #include "command.hpp"
  13 #include "sequence.hpp"
  14 #include "qualityscores.h"
  15 #include "alignment.hpp"
  16 #include "gotohoverlap.hpp"
  17 #include "needlemanoverlap.hpp"
  18 #include "blastalign.hpp"
  19 #include "noalign.hpp"
  20
  21
  22 struct fastqRead {
  23         vector<int> scores;
  24         string name;
  25         string sequence;
  26
  27         fastqRead() { name = ""; sequence = ""; scores.clear(); };
  28         fastqRead(string n, string s, vector<int> sc) : name(n), sequence(s), scores(sc) {};
  29         ~fastqRead() {};
  30 };
  31
  32 /**************************************************************************************************/
  33
  34 class MakeContigsCommand : public Command {
  35 public:
  36     MakeContigsCommand(string);
  37     MakeContigsCommand();
  38     ~MakeContigsCommand(){}
  39
  40     vector<string> setParameters();
  41     string getCommandName()                     { return "make.contigs";                        }
  42     string getCommandCategory()         { return "Sequence Processing";         }
  43     //commmand category choices: Sequence Processing, OTU-Based Approaches, Hypothesis Testing, Phylotype Analysis, General, Clustering and Hidden
  44     string getOutputFileNameTag(string, string);
  45         string getHelpString();
  46     string getCitation() { return "http://www.mothur.org/wiki/Make.contigs"; }
  47     string getDescription()             { return "description"; }
  48
  49     int execute();
  50     void help() { m->mothurOut(getHelpString()); }
  51
  52 private:
  53     bool abort;
  54     string outputDir, ffastqfile, rfastqfile, align;
  55         float match, misMatch, gapOpen, gapExtend;
  56         int processors, longestBase, threshold;
  57     vector<string> outputNames;
  58
  59     fastqRead readFastq(ifstream&);
  60     vector< vector<string> > readFastqFiles(int&);
  61     bool checkReads(fastqRead&, fastqRead&);
  62     int createProcesses(vector< vector<string> >, string, string, string);
  63     int driver(vector<string>, string, string, string);
  64 };
  65
  66 /**************************************************************************************************/
  67
  68 /**************************************************************************************************/
  69 //custom data structure for threads to use.
  70 // This is passed by void pointer so it can be any data type
  71 // that can be passed using a single void pointer (LPVOID).
  72 struct contigsData {
  73         string outputFasta;
  74         string outputQual;
  75         string outputMisMatches;
  76         string align;
  77     vector<string> files;
  78         MothurOut* m;
  79         float match, misMatch, gapOpen, gapExtend;
  80         int count, threshold, threadID;
  81
  82         contigsData(){}
  83         contigsData(vector<string> f, string of, string oq, string om, string al, MothurOut* mout, float ma, float misMa, float gapO, float gapE, int thr, int tid) {
  84         files = f;
  85                 outputFasta = of;
  86         outputQual = oq;
  87         outputMisMatches = om;
  88         m = mout;
  89                 match = ma;
  90                 misMatch = misMa;
  91                 gapOpen = gapO;
  92                 gapExtend = gapE;
  93         threshold = thr;
  94                 align = al;
  95                 count = 0;
  96                 threadID = tid;
  97         }
  98 };
  99
 100 /**************************************************************************************************/
 101 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 102 #else
 103 static DWORD WINAPI MyContigsThreadFunction(LPVOID lpParam){
 104         contigsData* pDataArray;
 105         pDataArray = (contigsData*)lpParam;
 106
 107         try {
 108         int longestBase = 1000;
 109         Alignment* alignment;
 110         if(pDataArray->align == "gotoh")                        {       alignment = new GotohOverlap(pDataArray->gapOpen, pDataArray->gapExtend, pDataArray->match, pDataArray->misMatch, longestBase);                 }
 111                 else if(pDataArray->align == "needleman")       {       alignment = new NeedlemanOverlap(pDataArray->gapOpen, pDataArray->match, pDataArray->misMatch, longestBase);                            }
 112
 113         int num = 0;
 114         string thisffastafile = pDataArray->files[0];
 115         string thisfqualfile = pDataArray->files[1];
 116         string thisrfastafile = pDataArray->files[2];
 117         string thisrqualfile = pDataArray->files[3];
 118
 119         if (pDataArray->m->debug) {  pDataArray->m->mothurOut("[DEBUG]: ffasta = " + thisffastafile + ".\n[DEBUG]: fqual = " + thisfqualfile + ".\n[DEBUG]: rfasta = " + thisrfastafile + ".\n[DEBUG]: rqual = " + thisrqualfile + ".\n"); }
 120
 121         ifstream inFFasta, inRFasta, inFQual, inRQual;
 122         pDataArray->m->openInputFile(thisffastafile, inFFasta);
 123         pDataArray->m->openInputFile(thisfqualfile, inFQual);
 124         pDataArray->m->openInputFile(thisrfastafile, inRFasta);
 125         pDataArray->m->openInputFile(thisrqualfile, inRQual);
 126
 127         ofstream outFasta, outQual, outMisMatch;
 128         pDataArray->m->openOutputFile(pDataArray->outputFasta, outFasta);
 129         pDataArray->m->openOutputFile(pDataArray->outputQual, outQual);
 130         pDataArray->m->openOutputFile(pDataArray->outputMisMatches, outMisMatch);
 131         outMisMatch << "Name\tLength\tMisMatches\n";
 132
 133         while ((!inFQual.eof()) && (!inFFasta.eof()) && (!inRFasta.eof()) && (!inRQual.eof())) {
 134
 135             if (pDataArray->m->control_pressed) { break; }
 136
 137             //read seqs and quality info
 138             Sequence fSeq(inFFasta); pDataArray->m->gobble(inFFasta);
 139             Sequence rSeq(inRFasta); pDataArray->m->gobble(inRFasta);
 140             QualityScores fQual(inFQual); pDataArray->m->gobble(inFQual);
 141             QualityScores rQual(inRQual); pDataArray->m->gobble(inRQual);
 142
 143             //flip the reverse reads
 144             rSeq.reverseComplement();
 145             rQual.flipQScores();
 146
 147             //pairwise align
 148             alignment->align(fSeq.getUnaligned(), rSeq.getUnaligned());
 149             map<int, int> ABaseMap = alignment->getSeqAAlnBaseMap();
 150             map<int, int> BBaseMap = alignment->getSeqBAlnBaseMap();
 151             fSeq.setAligned(alignment->getSeqAAln());
 152             rSeq.setAligned(alignment->getSeqBAln());
 153             int length = fSeq.getAligned().length();
 154
 155             //traverse alignments merging into one contiguous seq
 156             string contig = "";
 157             vector<int> contigScores;
 158             int numMismatches = 0;
 159             string seq1 = fSeq.getAligned();
 160             string seq2 = rSeq.getAligned();
 161
 162             vector<int> scores1 = fQual.getQualityScores();
 163             vector<int> scores2 = rQual.getQualityScores();
 164
 165             for (int i = 0; i < length; i++) {
 166                 if (seq1[i] == seq2[i]) { //match, add base and choose highest score
 167                     contig += seq1[i];
 168                     contigScores.push_back(scores1[ABaseMap[i]]);
 169                     if (scores1[ABaseMap[i]] < scores2[BBaseMap[i]]) { contigScores[i] = scores2[BBaseMap[i]]; }
 170                 }else if (((seq1[i] == '.') || (seq1[i] == '-')) && ((seq2[i] != '-') && (seq2[i] != '.'))) { //seq1 is a gap and seq2 is a base, choose seq2, unless quality score for base is below threshold. In that case eliminate base
 171                     if (scores2[BBaseMap[i]] >= pDataArray->threshold) {
 172                         contig += seq2[i];
 173                         contigScores.push_back(scores2[BBaseMap[i]]);
 174                     }
 175                 }else if (((seq2[i] == '.') || (seq2[i] == '-')) && ((seq1[i] != '-') && (seq1[i] != '.'))) { //seq2 is a gap and seq1 is a base, choose seq1, unless quality score for base is below threshold. In that case eliminate base
 176                     if (scores1[ABaseMap[i]] >= pDataArray->threshold) {
 177                         contig += seq1[i];
 178                         contigScores.push_back(scores1[ABaseMap[i]]);
 179                     }
 180                 }else if (((seq1[i] != '-') && (seq1[i] != '.')) && ((seq2[i] != '-') && (seq2[i] != '.'))) { //both bases choose one with better quality
 181                     char c = seq1[i];
 182                     contigScores.push_back(scores1[ABaseMap[i]]);
 183                     if (scores1[ABaseMap[i]] < scores2[BBaseMap[i]]) { contigScores[i] = scores2[BBaseMap[i]]; c = seq2[i]; }
 184                     contig += c;
 185                     numMismatches++;
 186                 }else { //should never get here
 187                     pDataArray->m->mothurOut("[ERROR]: case I didn't think of seq1 = " + toString(seq1[i]) + " and seq2 = " + toString(seq2[i]) + "\n");
 188                 }
 189             }
 190
 191             //output
 192             outFasta << ">" << fSeq.getName() << endl << contig << endl;
 193             outQual << ">" << fSeq.getName() << endl;
 194             for (int i = 0; i < contigScores.size(); i++) { outQual << contigScores[i] << ' '; }
 195             outQual << endl;
 196             outMisMatch << fSeq.getName() << '\t' << contig.length() << '\t' << numMismatches << endl;
 197
 198             num++;
 199
 200                         //report progress
 201                         if((num) % 1000 == 0){  pDataArray->m->mothurOut(toString(num)); pDataArray->m->mothurOutEndLine();             }
 202                 }
 203
 204                 //report progress
 205                 if((num) % 1000 != 0){  pDataArray->m->mothurOut(toString(num)); pDataArray->m->mothurOutEndLine();             }
 206
 207         inFFasta.close();
 208         inFQual.close();
 209         inRFasta.close();
 210         inRQual.close();
 211         outFasta.close();
 212         outQual.close();
 213         outMisMatch.close();
 214         delete alignment;
 215
 216         if (pDataArray->m->control_pressed) { pDataArray->m->mothurRemove(pDataArray->outputQual); pDataArray->m->mothurRemove(pDataArray->outputFasta);  pDataArray->m->mothurRemove(pDataArray->outputMisMatches);}
 217
 218         return 0;
 219
 220         }
 221         catch(exception& e) {
 222                 pDataArray->m->errorOut(e, "AlignCommand", "MyContigsThreadFunction");
 223                 exit(1);
 224         }
 225 }
 226 #endif
 227
 228
 229 #endif