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 getHelpString();
  45     string getCitation() { return "http://www.mothur.org/wiki/Make.contigs"; }
  46     string getDescription()             { return "description"; }
  47
  48     int execute();
  49     void help() { m->mothurOut(getHelpString()); }
  50
  51 private:
  52     bool abort;
  53     string outputDir, ffastqfile, rfastqfile, align;
  54         float match, misMatch, gapOpen, gapExtend;
  55         int processors, longestBase, threshold;
  56     vector<string> outputNames;
  57
  58     fastqRead readFastq(ifstream&);
  59     vector< vector<string> > readFastqFiles(int&);
  60     bool checkReads(fastqRead&, fastqRead&);
  61     int createProcesses(vector< vector<string> >, string, string, string);
  62     int driver(vector<string>, string, string, string);
  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 contigsData {
  72         string outputFasta;
  73         string outputQual;
  74         string outputMisMatches;
  75         string align;
  76     vector<string> files;
  77         MothurOut* m;
  78         float match, misMatch, gapOpen, gapExtend;
  79         int count, threshold, threadID;
  80
  81         contigsData(){}
  82         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) {
  83         files = f;
  84                 outputFasta = of;
  85         outputQual = oq;
  86         outputMisMatches = om;
  87         m = mout;
  88                 match = ma;
  89                 misMatch = misMa;
  90                 gapOpen = gapO;
  91                 gapExtend = gapE;
  92         threshold = thr;
  93                 align = al;
  94                 count = 0;
  95                 threadID = tid;
  96         }
  97 };
  98
  99 /**************************************************************************************************/
 100 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 101 #else
 102 static DWORD WINAPI MyContigsThreadFunction(LPVOID lpParam){
 103         contigsData* pDataArray;
 104         pDataArray = (contigsData*)lpParam;
 105
 106         try {
 107         int longestBase = 1000;
 108         Alignment* alignment;
 109         if(pDataArray->align == "gotoh")                        {       alignment = new GotohOverlap(pDataArray->gapOpen, pDataArray->gapExtend, pDataArray->match, pDataArray->misMatch, longestBase);                 }
 110                 else if(pDataArray->align == "needleman")       {       alignment = new NeedlemanOverlap(pDataArray->gapOpen, pDataArray->match, pDataArray->misMatch, longestBase);                            }
 111
 112         int num = 0;
 113         string thisffastafile = pDataArray->files[0];
 114         string thisfqualfile = pDataArray->files[1];
 115         string thisrfastafile = pDataArray->files[2];
 116         string thisrqualfile = pDataArray->files[3];
 117
 118         if (pDataArray->m->debug) {  pDataArray->m->mothurOut("[DEBUG]: ffasta = " + thisffastafile + ".\n[DEBUG]: fqual = " + thisfqualfile + ".\n[DEBUG]: rfasta = " + thisrfastafile + ".\n[DEBUG]: rqual = " + thisrqualfile + ".\n"); }
 119
 120         ifstream inFFasta, inRFasta, inFQual, inRQual;
 121         pDataArray->m->openInputFile(thisffastafile, inFFasta);
 122         pDataArray->m->openInputFile(thisfqualfile, inFQual);
 123         pDataArray->m->openInputFile(thisrfastafile, inRFasta);
 124         pDataArray->m->openInputFile(thisrqualfile, inRQual);
 125
 126         ofstream outFasta, outQual, outMisMatch;
 127         pDataArray->m->openOutputFile(pDataArray->outputFasta, outFasta);
 128         pDataArray->m->openOutputFile(pDataArray->outputQual, outQual);
 129         pDataArray->m->openOutputFile(pDataArray->outputMisMatches, outMisMatch);
 130         outMisMatch << "Name\tLength\tMisMatches\n";
 131
 132         while ((!inFQual.eof()) && (!inFFasta.eof()) && (!inRFasta.eof()) && (!inRQual.eof())) {
 133
 134             if (pDataArray->m->control_pressed) { break; }
 135
 136             //read seqs and quality info
 137             Sequence fSeq(inFFasta); pDataArray->m->gobble(inFFasta);
 138             Sequence rSeq(inRFasta); pDataArray->m->gobble(inRFasta);
 139             QualityScores fQual(inFQual); pDataArray->m->gobble(inFQual);
 140             QualityScores rQual(inRQual); pDataArray->m->gobble(inRQual);
 141
 142             //flip the reverse reads
 143             rSeq.reverseComplement();
 144             rQual.flipQScores();
 145
 146             //pairwise align
 147             alignment->align(fSeq.getUnaligned(), rSeq.getUnaligned());
 148             map<int, int> ABaseMap = alignment->getSeqAAlnBaseMap();
 149             map<int, int> BBaseMap = alignment->getSeqBAlnBaseMap();
 150             fSeq.setAligned(alignment->getSeqAAln());
 151             rSeq.setAligned(alignment->getSeqBAln());
 152             int length = fSeq.getAligned().length();
 153
 154             //traverse alignments merging into one contiguous seq
 155             string contig = "";
 156             vector<int> contigScores;
 157             int numMismatches = 0;
 158             string seq1 = fSeq.getAligned();
 159             string seq2 = rSeq.getAligned();
 160
 161             vector<int> scores1 = fQual.getQualityScores();
 162             vector<int> scores2 = rQual.getQualityScores();
 163
 164             for (int i = 0; i < length; i++) {
 165                 if (seq1[i] == seq2[i]) { //match, add base and choose highest score
 166                     contig += seq1[i];
 167                     contigScores.push_back(scores1[ABaseMap[i]]);
 168                     if (scores1[ABaseMap[i]] < scores2[BBaseMap[i]]) { contigScores[i] = scores2[BBaseMap[i]]; }
 169                 }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
 170                     if (scores2[BBaseMap[i]] >= pDataArray->threshold)) {
 171                         contig += seq2[i];
 172                         contigScores.push_back(scores2[BBaseMap[i]]);
 173                     }
 174                 }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
 175                     if (scores1[ABaseMap[i]] >= pDataArray->threshold) {
 176                         contig += seq1[i];
 177                         contigScores.push_back(scores1[ABaseMap[i]]);
 178                     }
 179                 }else if (((seq1[i] != '-') && (seq1[i] != '.')) && ((seq2[i] != '-') && (seq2[i] != '.'))) { //both bases choose one with better quality
 180                     char c = seq1[i];
 181                     contigScores.push_back(scores1[ABaseMap[i]]);
 182                     if (scores1[ABaseMap[i]] < scores2[BBaseMap[i]]) { contigScores[i] = scores2[BBaseMap[i]]; c = seq2[i]; }
 183                     contig += c;
 184                     numMismatches++;
 185                 }else { //should never get here
 186                     pDataArray->m->mothurOut("[ERROR]: case I didn't think of seq1 = " + toString(seq1[i]) + " and seq2 = " + toString(seq2[i]) + "\n");
 187                 }
 188             }
 189
 190             //output
 191             outFasta << ">" << fSeq.getName() << endl << contig << endl;
 192             outQual << ">" << fSeq.getName() << endl;
 193             for (int i = 0; i < contigScores.size(); i++) { outQual << contigScores[i] << ' '; }
 194             outQual << endl;
 195             outMisMatch << fSeq.getName() << '\t' << contig.length() << '\t' << numMismatches << endl;
 196
 197             num++;
 198
 199                         //report progress
 200                         if((num) % 1000 == 0){  pDataArray->m->mothurOut(toString(num)); pDataArray->m->mothurOutEndLine();             }
 201                 }
 202
 203                 //report progress
 204                 if((num) % 1000 != 0){  pDataArray->m->mothurOut(toString(num)); pDataArray->m->mothurOutEndLine();             }
 205
 206         inFFasta.close();
 207         inFQual.close();
 208         inRFasta.close();
 209         inRQual.close();
 210         outFasta.close();
 211         outQual.close();
 212         outMisMatch.close();
 213         delete alignment;
 214
 215         if (pDataArray->m->control_pressed) { pDataArray->m->mothurRemove(pDataArray->outputQual); pDataArray->m->mothurRemove(pDataArray->outputFasta);  pDataArray->m->mothurRemove(pDataArray->outputMisMatches);}
 216
 217         return 0;
 218
 219         }
 220         catch(exception& e) {
 221                 pDataArray->m->errorOut(e, "AlignCommand", "MyContigsThreadFunction");
 222                 exit(1);
 223         }
 224 }
 225 #endif
 226
 227
 228 #endif