summaryqualcommand.h

   1 #ifndef SUMMARYQUALCOMMAND_H
   2 #define SUMMARYQUALCOMMAND_H
   3
   4 /*
   5  *  summaryqualcommand.h
   6  *  Mothur
   7  *
   8  *  Created by westcott on 11/28/11.
   9  *  Copyright 2011 Schloss Lab. All rights reserved.
  10  *
  11  */
  12
  13
  14 #include "command.hpp"
  15 #include "qualityscores.h"
  16
  17 /**************************************************************************************************/
  18
  19 class SummaryQualCommand : public Command {
  20 public:
  21         SummaryQualCommand(string);
  22         SummaryQualCommand();
  23         ~SummaryQualCommand(){}
  24
  25         vector<string> setParameters();
  26         string getCommandName()                 { return "summary.qual";                        }
  27         string getCommandCategory()             { return "Sequence Processing";         }
  28
  29         string getHelpString();
  30     string getOutputPattern(string);
  31         string getCitation() { return "http://www.mothur.org/wiki/Summary.qual"; }
  32         string getDescription()         { return "summarize the quality of a set of sequences"; }
  33
  34         int execute();
  35         void help() { m->mothurOut(getHelpString()); }
  36
  37 private:
  38         bool abort;
  39         string qualfile, outputDir, namefile, countfile;
  40         vector<string> outputNames;
  41         map<string, int> nameMap;
  42         int processors;
  43
  44         struct linePair {
  45                 unsigned long long start;
  46                 unsigned long long end;
  47                 linePair(unsigned long long i, unsigned long long j) : start(i), end(j) {}
  48         };
  49
  50         vector<linePair> lines;
  51         vector<int> processIDS;
  52
  53         int createProcessesCreateSummary(vector<int>&, vector<int>&, vector< vector<int> >&, string);
  54         int driverCreateSummary(vector<int>&, vector<int>&, vector< vector<int> >&, string, linePair);
  55         int printQual(string, vector<int>&, vector<int>&, vector< vector<int> >&);
  56 };
  57
  58 /**************************************************************************************************/
  59 //custom data structure for threads to use.
  60 // This is passed by void pointer so it can be any data type
  61 // that can be passed using a single void pointer (LPVOID).
  62 struct seqSumQualData {
  63         vector<int> position;
  64         vector<int> averageQ;
  65         vector< vector<int> > scores;
  66         string filename;
  67         unsigned long long start;
  68         unsigned long long end;
  69         int count, numSeqs;
  70         MothurOut* m;
  71     bool hasNameMap;
  72         map<string, int> nameMap;
  73
  74         ~seqSumQualData(){}
  75         seqSumQualData(string f, MothurOut* mout, unsigned long long st, unsigned long long en, bool n, map<string, int> nam) {
  76                 filename = f;
  77                 m = mout;
  78                 start = st;
  79                 end = en;
  80                 hasNameMap = n;
  81                 nameMap = nam;
  82                 count = 0;
  83         }
  84 };
  85
  86 /**************************************************************************************************/
  87 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
  88 #else
  89 static DWORD WINAPI MySeqSumQualThreadFunction(LPVOID lpParam){
  90         seqSumQualData* pDataArray;
  91         pDataArray = (seqSumQualData*)lpParam;
  92
  93         try {
  94                 ifstream in;
  95                 pDataArray->m->openInputFile(pDataArray->filename, in);
  96
  97                 //print header if you are process 0
  98                 if ((pDataArray->start == 0) || (pDataArray->start == 1)) {
  99                         in.seekg(0);
 100                 }else { //this accounts for the difference in line endings.
 101                         in.seekg(pDataArray->start-1); pDataArray->m->gobble(in);
 102                 }
 103
 104                 pDataArray->count = 0;
 105         pDataArray->numSeqs = 0;
 106                 for(int i = 0; i < pDataArray->end; i++){ //end is the number of sequences to process
 107
 108                         if (pDataArray->m->control_pressed) { in.close(); pDataArray->count = 1; return 1; }
 109
 110                         QualityScores current(in); pDataArray->m->gobble(in);
 111
 112                         if (current.getName() != "") {
 113
 114                                 int num = 1;
 115                                 if (pDataArray->hasNameMap) {
 116                                         //make sure this sequence is in the namefile, else error
 117                                         map<string, int>::iterator it = pDataArray->nameMap.find(current.getName());
 118
 119                                         if (it == pDataArray->nameMap.end()) { pDataArray->m->mothurOut("[ERROR]: " + current.getName() + " is not in your namefile, please correct."); pDataArray->m->mothurOutEndLine(); pDataArray->m->control_pressed = true; }
 120                                         else { num = it->second; }
 121                                 }
 122
 123                                 vector<int> thisScores = current.getQualityScores();
 124
 125                                 //resize to num of positions setting number of seqs with that size to 1
 126                                 if (pDataArray->position.size() < thisScores.size()) { pDataArray->position.resize(thisScores.size(), 0); }
 127                                 if (pDataArray->averageQ.size() < thisScores.size()) { pDataArray->averageQ.resize(thisScores.size(), 0); }
 128                                 if (pDataArray->scores.size() < thisScores.size()) {
 129                                         pDataArray->scores.resize(thisScores.size());
 130                                         for (int i = 0; i < pDataArray->scores.size(); i++) { pDataArray->scores.at(i).resize(41, 0); }
 131                                 }
 132
 133                                 //increase counts of number of seqs with this position
 134                                 //average is really the total, we will average in execute
 135                                 for (int i = 0; i < thisScores.size(); i++) {
 136                                         pDataArray->position.at(i) += num;
 137                                         pDataArray->averageQ.at(i) += (thisScores[i] * num); //weighting for namesfile
 138                                         if (thisScores[i] > 40) { pDataArray->m->mothurOut("[ERROR]: " + current.getName() + " has a quality scores of " + toString(thisScores[i]) + ", expecting values to be less than 40."); pDataArray->m->mothurOutEndLine(); pDataArray->m->control_pressed = true; }
 139                                         else { pDataArray->scores.at(i)[thisScores[i]] += num; }
 140                                 }
 141
 142                                 pDataArray->numSeqs += num;
 143                 pDataArray->count++;
 144                         }
 145                 }
 146
 147                 in.close();
 148
 149                 return 0;
 150
 151         }
 152         catch(exception& e) {
 153                 pDataArray->m->errorOut(e, "SummaryQualCommand", "MySeqSumQualThreadFunction");
 154                 exit(1);
 155         }
 156 }
 157 #endif
 158
 159
 160 /**************************************************************************************************/
 161
 162
 163 #endif
 164