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;
  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, 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 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                 align = al;
  93                 count = 0;
  94                 threadID = tid;
  95         }
  96 };
  97
  98 /**************************************************************************************************/
  99 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 100 #else
 101 static DWORD WINAPI MyContigsThreadFunction(LPVOID lpParam){
 102         contigsData* pDataArray;
 103         pDataArray = (contigsData*)lpParam;
 104
 105         try {
 106         int longestBase = 1000;
 107         Alignment* alignment;
 108         if(pDataArray->align == "gotoh")                        {       alignment = new GotohOverlap(pDataArray->gapOpen, pDataArray->gapExtend, pDataArray->match, pDataArray->misMatch, longestBase);                 }
 109                 else if(pDataArray->align == "needleman")       {       alignment = new NeedlemanOverlap(pDataArray->gapOpen, pDataArray->match, pDataArray->misMatch, longestBase);                            }
 110
 111         int num = 0;
 112         string thisffastafile = pDataArray->files[0];
 113         string thisfqualfile = pDataArray->files[1];
 114         string thisrfastafile = pDataArray->files[2];
 115         string thisrqualfile = pDataArray->files[3];
 116
 117         if (pDataArray->m->debug) {  pDataArray->m->mothurOut("[DEBUG]: ffasta = " + thisffastafile + ".\n[DEBUG]: fqual = " + thisfqualfile + ".\n[DEBUG]: rfasta = " + thisrfastafile + ".\n[DEBUG]: rqual = " + thisrqualfile + ".\n"); }
 118
 119         ifstream inFFasta, inRFasta, inFQual, inRQual;
 120         pDataArray->m->openInputFile(thisffastafile, inFFasta);
 121         pDataArray->m->openInputFile(thisfqualfile, inFQual);
 122         pDataArray->m->openInputFile(thisrfastafile, inRFasta);
 123         pDataArray->m->openInputFile(thisrqualfile, inRQual);
 124
 125         ofstream outFasta, outQual, outMisMatch;
 126         pDataArray->m->openOutputFile(pDataArray->outputFasta, outFasta);
 127         pDataArray->m->openOutputFile(pDataArray->outputQual, outQual);
 128         pDataArray->m->openOutputFile(pDataArray->outputMisMatches, outMisMatch);
 129         outMisMatch << "Name\tLength\tMisMatches\n";
 130
 131         while ((!inFQual.eof()) && (!inFFasta.eof()) && (!inRFasta.eof()) && (!inRQual.eof())) {
 132
 133             if (pDataArray->m->control_pressed) { break; }
 134
 135             //read seqs and quality info
 136             Sequence fSeq(inFFasta); pDataArray->m->gobble(inFFasta);
 137             Sequence rSeq(inRFasta); pDataArray->m->gobble(inRFasta);
 138             QualityScores fQual(inFQual); pDataArray->m->gobble(inFQual);
 139             QualityScores rQual(inRQual); pDataArray->m->gobble(inRQual);
 140
 141             //flip the reverse reads
 142             rSeq.reverseComplement();
 143             rQual.flipQScores();
 144
 145             //pairwise align
 146             alignment->align(fSeq.getUnaligned(), rSeq.getUnaligned());
 147             map<int, int> ABaseMap = alignment->getSeqAAlnBaseMap();
 148             map<int, int> BBaseMap = alignment->getSeqBAlnBaseMap();
 149             fSeq.setAligned(alignment->getSeqAAln());
 150             rSeq.setAligned(alignment->getSeqBAln());
 151             int length = fSeq.getAligned().length();
 152
 153             //traverse alignments merging into one contiguous seq
 154             string contig = "";
 155             vector<int> contigScores;
 156             int numMismatches = 0;
 157             string seq1 = fSeq.getAligned();
 158             string seq2 = rSeq.getAligned();
 159
 160             vector<int> scores1 = fQual.getQualityScores();
 161             vector<int> scores2 = rQual.getQualityScores();
 162
 163             for (int i = 0; i < length; i++) {
 164                 if (seq1[i] == seq2[i]) { //match, add base and choose highest score
 165                     contig += seq1[i];
 166                     contigScores.push_back(scores1[ABaseMap[i]]);
 167                     if (scores1[ABaseMap[i]] < scores2[BBaseMap[i]]) { contigScores[i] = scores2[BBaseMap[i]]; }
 168                 }else if (((seq1[i] == '.') || (seq1[i] == '-')) && ((seq2[i] != '-') && (seq2[i] != '.'))) { //seq1 is a gap and seq2 is a base, choose seq2
 169                     contig += seq2[i];
 170                     contigScores.push_back(scores2[BBaseMap[i]]);
 171                 }else if (((seq2[i] == '.') || (seq2[i] == '-')) && ((seq1[i] != '-') && (seq1[i] != '.'))) { //seq2 is a gap and seq1 is a base, choose seq1
 172                     contig += seq1[i];
 173                     contigScores.push_back(scores1[ABaseMap[i]]);
 174                 }else if (((seq1[i] != '-') && (seq1[i] != '.')) && ((seq2[i] != '-') && (seq2[i] != '.'))) { //both bases choose one with better quality
 175                     char c = seq1[i];
 176                     contigScores.push_back(scores1[ABaseMap[i]]);
 177                     if (scores1[ABaseMap[i]] < scores2[BBaseMap[i]]) { contigScores[i] = scores2[BBaseMap[i]]; c = seq2[i]; }
 178                     contig += c;
 179                     numMismatches++;
 180                 }else { //should never get here
 181                     pDataArray->m->mothurOut("[ERROR]: case I didn't think of seq1 = " + toString(seq1[i]) + " and seq2 = " + toString(seq2[i]) + "\n");
 182                 }
 183             }
 184
 185             //output
 186             outFasta << ">" << fSeq.getName() << endl << contig << endl;
 187             outQual << ">" << fSeq.getName() << endl;
 188             for (int i = 0; i < contigScores.size(); i++) { outQual << contigScores[i] << ' '; }
 189             outQual << endl;
 190             outMisMatch << fSeq.getName() << '\t' << contig.length() << '\t' << numMismatches << endl;
 191
 192             num++;
 193
 194                         //report progress
 195                         if((num) % 1000 == 0){  pDataArray->m->mothurOut(toString(num)); pDataArray->m->mothurOutEndLine();             }
 196                 }
 197
 198                 //report progress
 199                 if((num) % 1000 != 0){  pDataArray->m->mothurOut(toString(num)); pDataArray->m->mothurOutEndLine();             }
 200
 201         inFFasta.close();
 202         inFQual.close();
 203         inRFasta.close();
 204         inRQual.close();
 205         outFasta.close();
 206         outQual.close();
 207         outMisMatch.close();
 208         delete alignment;
 209
 210         if (pDataArray->m->control_pressed) { pDataArray->m->mothurRemove(pDataArray->outputQual); pDataArray->m->mothurRemove(pDataArray->outputFasta);  pDataArray->m->mothurRemove(pDataArray->outputMisMatches);}
 211
 212         return 0;
 213
 214         }
 215         catch(exception& e) {
 216                 pDataArray->m->errorOut(e, "AlignCommand", "MyContigsThreadFunction");
 217                 exit(1);
 218         }
 219 }
 220 #endif
 221
 222
 223 #endif