seqerrorcommand.cpp

   1 /*
   2  *  seqerrorcommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by Pat Schloss on 7/15/10.
   6  *  Copyright 2010 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "seqerrorcommand.h"
  11 #include "reportfile.h"
  12 #include "qualityscores.h"
  13 #include "refchimeratest.h"
  14 #include "filterseqscommand.h"
  15
  16 //**********************************************************************************************************************
  17 vector<string> SeqErrorCommand::setParameters(){
  18         try {
  19                 CommandParameter pquery("fasta", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pquery);
  20                 CommandParameter preference("reference", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(preference);
  21                 CommandParameter pqfile("qfile", "InputTypes", "", "", "none", "none", "QualReport",false,false); parameters.push_back(pqfile);
  22                 CommandParameter preport("report", "InputTypes", "", "", "none", "none", "QualReport",false,false); parameters.push_back(preport);
  23                 CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname);
  24                 CommandParameter pignorechimeras("ignorechimeras", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(pignorechimeras);
  25                 CommandParameter pthreshold("threshold", "Number", "", "1.0", "", "", "",false,false); parameters.push_back(pthreshold);
  26                 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
  27                 CommandParameter pfilter("filter", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(pfilter);
  28                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
  29                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
  30
  31                 vector<string> myArray;
  32                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  33                 return myArray;
  34         }
  35         catch(exception& e) {
  36                 m->errorOut(e, "SeqErrorCommand", "setParameters");
  37                 exit(1);
  38         }
  39 }
  40 //**********************************************************************************************************************
  41 string SeqErrorCommand::getHelpString(){
  42         try {
  43                 string helpString = "";
  44                 helpString += "The seq.error command reads a query alignment file and a reference alignment file and creates .....\n";
  45                 helpString += "Example seq.error(...).\n";
  46                 helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n";
  47                 helpString += "For more details please check out the wiki http://www.mothur.org/wiki/seq.error .\n";
  48                 return helpString;
  49         }
  50         catch(exception& e) {
  51                 m->errorOut(e, "SeqErrorCommand", "getHelpString");
  52                 exit(1);
  53         }
  54 }
  55 //**********************************************************************************************************************
  56 SeqErrorCommand::SeqErrorCommand(){
  57         try {
  58                 abort = true; calledHelp = true;
  59                 vector<string> tempOutNames;
  60                 outputTypes["error.summary"] = tempOutNames;
  61                 outputTypes["error.seq"] = tempOutNames;
  62                 outputTypes["error.quality"] = tempOutNames;
  63                 outputTypes["error.qual.forward"] = tempOutNames;
  64                 outputTypes["error.qual.reverse"] = tempOutNames;
  65                 outputTypes["error.forward"] = tempOutNames;
  66                 outputTypes["error.reverse"] = tempOutNames;
  67                 outputTypes["error.count"] = tempOutNames;
  68                 outputTypes["error.matrix"] = tempOutNames;
  69         }
  70         catch(exception& e) {
  71                 m->errorOut(e, "SeqErrorCommand", "SeqErrorCommand");
  72                 exit(1);
  73         }
  74 }
  75 //***************************************************************************************************************
  76
  77 SeqErrorCommand::SeqErrorCommand(string option)  {
  78         try {
  79
  80                 abort = false; calledHelp = false;
  81
  82                 //allow user to run help
  83                 if(option == "help") { help(); abort = true; calledHelp = true; }
  84                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
  85
  86                 else {
  87                         string temp;
  88                         vector<string> myArray = setParameters();
  89
  90                         OptionParser parser(option);
  91                         map<string,string> parameters = parser.getParameters();
  92
  93                         ValidParameters validParameter;
  94                         map<string,string>::iterator it;
  95
  96                         //check to make sure all parameters are valid for command
  97                         for (it = parameters.begin(); it != parameters.end(); it++) {
  98                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  99                         }
 100
 101                         //initialize outputTypes
 102                         vector<string> tempOutNames;
 103                         outputTypes["error.summary"] = tempOutNames;
 104                         outputTypes["error.seq"] = tempOutNames;
 105                         outputTypes["error.quality"] = tempOutNames;
 106                         outputTypes["error.qual.forward"] = tempOutNames;
 107                         outputTypes["error.qual.reverse"] = tempOutNames;
 108                         outputTypes["error.forward"] = tempOutNames;
 109                         outputTypes["error.reverse"] = tempOutNames;
 110                         outputTypes["error.count"] = tempOutNames;
 111                         outputTypes["error.matrix"] = tempOutNames;
 112
 113
 114                         //if the user changes the input directory command factory will send this info to us in the output parameter
 115                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 116                         if (inputDir == "not found"){   inputDir = "";          }
 117                         else {
 118                                 string path;
 119                                 it = parameters.find("fasta");
 120                                 //user has given a template file
 121                                 if(it != parameters.end()){
 122                                         path = m->hasPath(it->second);
 123                                         //if the user has not given a path then, add inputdir. else leave path alone.
 124                                         if (path == "") {       parameters["fasta"] = inputDir + it->second;            }
 125                                 }
 126
 127                                 it = parameters.find("reference");
 128                                 //user has given a template file
 129                                 if(it != parameters.end()){
 130                                         path = m->hasPath(it->second);
 131                                         //if the user has not given a path then, add inputdir. else leave path alone.
 132                                         if (path == "") {       parameters["reference"] = inputDir + it->second;                }
 133                                 }
 134
 135                                 it = parameters.find("name");
 136                                 //user has given a names file
 137                                 if(it != parameters.end()){
 138                                         path = m->hasPath(it->second);
 139                                         //if the user has not given a path then, add inputdir. else leave path alone.
 140                                         if (path == "") {       parameters["name"] = inputDir + it->second;             }
 141                                 }
 142
 143                                 it = parameters.find("qfile");
 144                                 //user has given a quality score file
 145                                 if(it != parameters.end()){
 146                                         path = m->hasPath(it->second);
 147                                         //if the user has not given a path then, add inputdir. else leave path alone.
 148                                         if (path == "") {       parameters["qfile"] = inputDir + it->second;            }
 149                                 }
 150
 151                                 it = parameters.find("report");
 152                                 //user has given a alignment report file
 153                                 if(it != parameters.end()){
 154                                         path = m->hasPath(it->second);
 155                                         //if the user has not given a path then, add inputdir. else leave path alone.
 156                                         if (path == "") {       parameters["report"] = inputDir + it->second;           }
 157                                 }
 158
 159                         }
 160                         //check for required parameters
 161                         queryFileName = validParameter.validFile(parameters, "fasta", true);
 162                         if (queryFileName == "not found") {
 163                                 queryFileName = m->getFastaFile();
 164                                 if (queryFileName != "") { m->mothurOut("Using " + queryFileName + " as input file for the fasta parameter."); m->mothurOutEndLine(); }
 165                                 else {  m->mothurOut("You have no current fasta file and the fasta parameter is required."); m->mothurOutEndLine(); abort = true; }
 166                         }
 167                         else if (queryFileName == "not open") { abort = true; }
 168
 169                         referenceFileName = validParameter.validFile(parameters, "reference", true);
 170                         if (referenceFileName == "not found") { m->mothurOut("reference is a required parameter for the seq.error command."); m->mothurOutEndLine(); abort = true; }
 171                         else if (referenceFileName == "not open") { abort = true; }
 172
 173
 174                         //check for optional parameters
 175                         namesFileName = validParameter.validFile(parameters, "name", true);
 176                         if(namesFileName == "not found"){       namesFileName = "";     }
 177                         else if (namesFileName == "not open") { namesFileName = ""; abort = true; }
 178
 179                         qualFileName = validParameter.validFile(parameters, "qfile", true);
 180                         if(qualFileName == "not found"){        qualFileName = "";      }
 181                         else if (qualFileName == "not open") { qualFileName = ""; abort = true; }
 182
 183                         reportFileName = validParameter.validFile(parameters, "report", true);
 184                         if(reportFileName == "not found"){      reportFileName = "";    }
 185                         else if (reportFileName == "not open") { reportFileName = ""; abort = true; }
 186
 187                         if((reportFileName != "" && qualFileName == "") || (reportFileName == "" && qualFileName != "")){
 188                                 m->mothurOut("if you use either a qual file or a report file, you have to have both.");
 189                                 m->mothurOutEndLine();
 190                                 abort = true;
 191                         }
 192
 193                         outputDir = validParameter.validFile(parameters, "outputdir", false);
 194                         if (outputDir == "not found"){
 195                                 outputDir = "";
 196                                 outputDir += m->hasPath(queryFileName); //if user entered a file with a path then preserve it
 197                         }
 198
 199                         //check for optional parameter and set defaults
 200                         // ...at some point should added some additional type checking...
 201                         temp = validParameter.validFile(parameters, "threshold", false);        if (temp == "not found") { temp = "1.00"; }
 202                         convert(temp, threshold);
 203
 204                         temp = validParameter.validFile(parameters, "ignorechimeras", false);   if (temp == "not found") { temp = "T"; }
 205                         ignoreChimeras = m->isTrue(temp);
 206
 207                         temp = validParameter.validFile(parameters, "filter", false);   if (temp == "not found") { temp = "T"; }
 208                         filter = m->isTrue(temp);
 209
 210                         temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){       temp = m->getProcessors();      }
 211                         m->setProcessors(temp);
 212                         convert(temp, processors);
 213
 214                         substitutionMatrix.resize(6);
 215                         for(int i=0;i<6;i++){   substitutionMatrix[i].resize(6,0);      }
 216                 }
 217         }
 218         catch(exception& e) {
 219                 m->errorOut(e, "SeqErrorCommand", "SeqErrorCommand");
 220                 exit(1);
 221         }
 222 }
 223 //***************************************************************************************************************
 224
 225 int SeqErrorCommand::execute(){
 226         try{
 227                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 228
 229                 int start = time(NULL);
 230                 maxLength = 2000;
 231                 totalBases = 0;
 232                 totalMatches = 0;
 233
 234                 //run vertical filter on query and reference files.
 235                 if (filter) {
 236                         string inputString = "fasta=" + queryFileName + "-" + referenceFileName;
 237                         m->mothurOut("/******************************************/"); m->mothurOutEndLine();
 238                         m->mothurOut("Running command: filter.seqs(" + inputString + ") to improve processing time."); m->mothurOutEndLine();
 239
 240                         Command* filterCommand = new FilterSeqsCommand(inputString);
 241                         filterCommand->execute();
 242
 243                         map<string, vector<string> > filenames = filterCommand->getOutputFiles();
 244
 245                         delete filterCommand;
 246
 247                         m->mothurOut("/******************************************/"); m->mothurOutEndLine();
 248
 249                         queryFileName = filenames["fasta"][0];
 250                         referenceFileName = filenames["fasta"][1];
 251                 }
 252
 253                 string errorSummaryFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.summary";
 254                 outputNames.push_back(errorSummaryFileName); outputTypes["error.summary"].push_back(errorSummaryFileName);
 255
 256                 string errorSeqFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq";
 257                 outputNames.push_back(errorSeqFileName); outputTypes["error.seq"].push_back(errorSeqFileName);
 258
 259                 string errorChimeraFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.chimera";
 260                 outputNames.push_back(errorChimeraFileName); outputTypes["error.chimera"].push_back(errorChimeraFileName);
 261
 262                 getReferences();        //read in reference sequences - make sure there's no ambiguous bases
 263
 264                 if(namesFileName != ""){        weights = getWeights(); }
 265
 266                 vector<unsigned long int> fastaFilePos;
 267                 vector<unsigned long int> qFilePos;
 268                 vector<unsigned long int> reportFilePos;
 269
 270                 setLines(queryFileName, qualFileName, reportFileName, fastaFilePos, qFilePos, reportFilePos);
 271
 272                 if (m->control_pressed) { return 0; }
 273
 274                 for (int i = 0; i < (fastaFilePos.size()-1); i++) {
 275                         lines.push_back(linePair(fastaFilePos[i], fastaFilePos[(i+1)]));
 276                         if (qualFileName != "") {  qLines.push_back(linePair(qFilePos[i], qFilePos[(i+1)]));  }
 277                         if (reportFileName != "") {  rLines.push_back(linePair(reportFilePos[i], reportFilePos[(i+1)]));  }
 278                 }
 279                 if(qualFileName == "")  {       qLines = lines; rLines = lines; } //fills with duds
 280
 281                 int numSeqs = 0;
 282 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 283                 if(processors == 1){
 284                         numSeqs = driver(queryFileName, qualFileName, reportFileName, errorSummaryFileName, errorSeqFileName, errorChimeraFileName, lines[0], qLines[0], rLines[0]);
 285                 }else{
 286                         numSeqs = createProcesses(queryFileName, qualFileName, reportFileName, errorSummaryFileName, errorSeqFileName, errorChimeraFileName);
 287                 }
 288 #else
 289                 numSeqs = driver(queryFileName, qualFileName, reportFileName, errorSummaryFileName, errorSeqFileName, errorChimeraFileName, lines[0], qLines[0], rLines[0]);
 290 #endif
 291
 292                 if(qualFileName != "" && reportFileName != ""){
 293                         printErrorQuality(qScoreErrorMap);
 294                         printQualityFR(qualForwardMap, qualReverseMap);
 295                 }
 296
 297                 printErrorFRFile(errorForward, errorReverse);
 298
 299                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
 300
 301                 string errorCountFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.count";
 302                 ofstream errorCountFile;
 303                 m->openOutputFile(errorCountFileName, errorCountFile);
 304                 outputNames.push_back(errorCountFileName);  outputTypes["error.count"].push_back(errorCountFileName);
 305                 m->mothurOut("Overall error rate:\t" + toString((double)(totalBases - totalMatches) / (double)totalBases) + "\n");
 306                 m->mothurOut("Errors\tSequences\n");
 307                 errorCountFile << "Errors\tSequences\n";
 308                 for(int i=0;i<misMatchCounts.size();i++){
 309                         m->mothurOut(toString(i) + '\t' + toString(misMatchCounts[i]) + '\n');
 310                         errorCountFile << i << '\t' << misMatchCounts[i] << endl;
 311                 }
 312                 errorCountFile.close();
 313
 314                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
 315
 316                 printSubMatrix();
 317
 318                 string megAlignmentFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.ref-query";
 319                 ofstream megAlignmentFile;
 320                 m->openOutputFile(megAlignmentFileName, megAlignmentFile);
 321                 outputNames.push_back(megAlignmentFileName);  outputTypes["error.ref-query"].push_back(megAlignmentFileName);
 322
 323                 for(int i=0;i<numRefs;i++){
 324                         megAlignmentFile << referenceSeqs[i].getInlineSeq() << endl;
 325                         megAlignmentFile << megaAlignVector[i] << endl;
 326                 }
 327
 328                 m->mothurOut("It took " + toString(time(NULL) - start) + " secs to check " + toString(numSeqs) + " sequences.");
 329                 m->mothurOutEndLine();
 330
 331                 m->mothurOutEndLine();
 332                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 333                 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
 334                 m->mothurOutEndLine();
 335
 336                 return 0;
 337         }
 338         catch(exception& e) {
 339                 m->errorOut(e, "SeqErrorCommand", "execute");
 340                 exit(1);
 341         }
 342 }
 343 //**********************************************************************************************************************
 344 int SeqErrorCommand::createProcesses(string filename, string qFileName, string rFileName, string summaryFileName, string errorOutputFileName, string chimeraOutputFileName) {
 345         try {
 346                 int process = 1;
 347                 processIDS.clear();
 348                 map<char, vector<int> >::iterator it;
 349                 int num = 0;
 350 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 351
 352                 //loop through and create all the processes you want
 353                 while (process != processors) {
 354                         int pid = fork();
 355
 356                         if (pid > 0) {
 357                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
 358                                 process++;
 359                         }else if (pid == 0){
 360
 361                                 num = driver(filename, qFileName, rFileName, summaryFileName + toString(getpid()) + ".temp", errorOutputFileName+ toString(getpid()) + ".temp", chimeraOutputFileName + toString(getpid()) + ".temp", lines[process], qLines[process], rLines[process]);
 362
 363                                 //pass groupCounts to parent
 364                                 ofstream out;
 365                                 string tempFile = filename + toString(getpid()) + ".info.temp";
 366                                 m->openOutputFile(tempFile, out);
 367
 368                                 //output totalBases and totalMatches
 369                                 out << num << '\t' << totalBases << '\t' << totalMatches << endl << endl;
 370
 371                                 //output substitutionMatrix
 372                                 for(int i = 0; i < substitutionMatrix.size(); i++) {
 373                                         for (int j = 0; j < substitutionMatrix[i].size(); j++) {
 374                                                 out << substitutionMatrix[i][j] << '\t';
 375                                         }
 376                                         out << endl;
 377                                 }
 378                                 out << endl;
 379
 380                                 //output qScoreErrorMap
 381                                 for (it = qScoreErrorMap.begin(); it != qScoreErrorMap.end(); it++) {
 382                                         vector<int> thisScoreErrorMap = it->second;
 383                                         out << it->first << '\t';
 384                                         for (int i = 0; i < thisScoreErrorMap.size(); i++) {
 385                                                 out << thisScoreErrorMap[i] << '\t';
 386                                         }
 387                                         out << endl;
 388                                 }
 389                                 out << endl;
 390
 391                                 //output qualForwardMap
 392                                 for(int i = 0; i < qualForwardMap.size(); i++) {
 393                                         for (int j = 0; j < qualForwardMap[i].size(); j++) {
 394                                                 out << qualForwardMap[i][j] << '\t';
 395                                         }
 396                                         out << endl;
 397                                 }
 398                                 out << endl;
 399
 400                                 //output qualReverseMap
 401                                 for(int i = 0; i < qualReverseMap.size(); i++) {
 402                                         for (int j = 0; j < qualReverseMap[i].size(); j++) {
 403                                                 out << qualReverseMap[i][j] << '\t';
 404                                         }
 405                                         out << endl;
 406                                 }
 407                                 out << endl;
 408
 409
 410                                 //output errorForward
 411                                 for (it = errorForward.begin(); it != errorForward.end(); it++) {
 412                                         vector<int> thisErrorForward = it->second;
 413                                         out << it->first << '\t';
 414                                         for (int i = 0; i < thisErrorForward.size(); i++) {
 415                                                 out << thisErrorForward[i] << '\t';
 416                                         }
 417                                         out << endl;
 418                                 }
 419                                 out << endl;
 420
 421                                 //output errorReverse
 422                                 for (it = errorReverse.begin(); it != errorReverse.end(); it++) {
 423                                         vector<int> thisErrorReverse = it->second;
 424                                         out << it->first << '\t';
 425                                         for (int i = 0; i < thisErrorReverse.size(); i++) {
 426                                                 out << thisErrorReverse[i] << '\t';
 427                                         }
 428                                         out << endl;
 429                                 }
 430                                 out << endl;
 431
 432                                 //output misMatchCounts
 433                                 out << misMatchCounts.size() << endl;
 434                                 for (int j = 0; j < misMatchCounts.size(); j++) {
 435                                         out << misMatchCounts[j] << '\t';
 436                                 }
 437                                 out << endl;
 438
 439
 440                                 //output megaAlignVector
 441                                 for (int j = 0; j < megaAlignVector.size(); j++) {
 442                                         out << megaAlignVector[j] << endl;
 443                                 }
 444                                 out << endl;
 445
 446                                 out.close();
 447
 448                                 exit(0);
 449                         }else {
 450                                 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
 451                                 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
 452                                 exit(0);
 453                         }
 454                 }
 455
 456                 //do my part
 457                 num = driver(filename, qFileName, rFileName, summaryFileName, errorOutputFileName, chimeraOutputFileName, lines[0], qLines[0], rLines[0]);
 458
 459                 //force parent to wait until all the processes are done
 460                 for (int i=0;i<processIDS.size();i++) {
 461                         int temp = processIDS[i];
 462                         wait(&temp);
 463                 }
 464
 465                 //append files
 466                 for(int i=0;i<processIDS.size();i++){
 467
 468                         m->mothurOut("Appending files from process " + toString(processIDS[i])); m->mothurOutEndLine();
 469
 470                         m->appendFiles((summaryFileName + toString(processIDS[i]) + ".temp"), summaryFileName);
 471                         remove((summaryFileName + toString(processIDS[i]) + ".temp").c_str());
 472                         m->appendFiles((errorOutputFileName + toString(processIDS[i]) + ".temp"), errorOutputFileName);
 473                         remove((errorOutputFileName + toString(processIDS[i]) + ".temp").c_str());
 474                         m->appendFiles((chimeraOutputFileName + toString(processIDS[i]) + ".temp"), chimeraOutputFileName);
 475                         remove((chimeraOutputFileName + toString(processIDS[i]) + ".temp").c_str());
 476
 477                         ifstream in;
 478                         string tempFile =  filename + toString(processIDS[i]) + ".info.temp";
 479                         m->openInputFile(tempFile, in);
 480
 481                         //input totalBases and totalMatches
 482                         int tempBases, tempMatches, tempNumSeqs;
 483                         in >> tempNumSeqs >> tempBases >> tempMatches; m->gobble(in);
 484                         totalBases += tempBases; totalMatches += tempMatches; num += tempNumSeqs;
 485
 486                         //input substitutionMatrix
 487                         int tempNum;
 488                         for(int i = 0; i < substitutionMatrix.size(); i++) {
 489                                 for (int j = 0; j < substitutionMatrix[i].size(); j++) {
 490                                         in >> tempNum; substitutionMatrix[i][j] += tempNum;
 491                                 }
 492                                 m->gobble(in);
 493                         }
 494                         m->gobble(in);
 495
 496                         //input qScoreErrorMap
 497                         char first;
 498                         for (int i = 0; i < qScoreErrorMap.size(); i++) {
 499                                 in >> first;
 500                                 vector<int> thisScoreErrorMap = qScoreErrorMap[first];
 501
 502                                 for (int i = 0; i < thisScoreErrorMap.size(); i++) {
 503                                         in >> tempNum; thisScoreErrorMap[i] += tempNum;
 504                                 }
 505                                 qScoreErrorMap[first] = thisScoreErrorMap;
 506                                 m->gobble(in);
 507                         }
 508                         m->gobble(in);
 509
 510                         //input qualForwardMap
 511                         for(int i = 0; i < qualForwardMap.size(); i++) {
 512                                 for (int j = 0; j < qualForwardMap[i].size(); j++) {
 513                                         in >> tempNum; qualForwardMap[i][j] += tempNum;
 514                                 }
 515                                 m->gobble(in);
 516                         }
 517                         m->gobble(in);
 518
 519                         //input qualReverseMap
 520                         for(int i = 0; i < qualReverseMap.size(); i++) {
 521                                 for (int j = 0; j < qualReverseMap[i].size(); j++) {
 522                                         in >> tempNum; qualReverseMap[i][j] += tempNum;
 523                                 }
 524                                 m->gobble(in);
 525                         }
 526                         m->gobble(in);
 527
 528                         //input errorForward
 529                         for (int i = 0; i < errorForward.size(); i++) {
 530                                 in >> first;
 531                                 vector<int> thisErrorForward = errorForward[first];
 532
 533                                 for (int i = 0; i < thisErrorForward.size(); i++) {
 534                                         in >> tempNum; thisErrorForward[i] += tempNum;
 535                                 }
 536                                 errorForward[first] = thisErrorForward;
 537                                 m->gobble(in);
 538                         }
 539                         m->gobble(in);
 540
 541                         //input errorReverse
 542                         for (int i = 0; i < errorReverse.size(); i++) {
 543                                 in >> first;
 544                                 vector<int> thisErrorReverse = errorReverse[first];
 545
 546                                 for (int i = 0; i < thisErrorReverse.size(); i++) {
 547                                         in >> tempNum; thisErrorReverse[i] += tempNum;
 548                                 }
 549                                 errorReverse[first] = thisErrorReverse;
 550                                 m->gobble(in);
 551                         }
 552                         m->gobble(in);
 553
 554                         //input misMatchCounts
 555                         int misMatchSize;
 556                         in >> misMatchSize; m->gobble(in);
 557                         if (misMatchSize > misMatchCounts.size()) {     misMatchCounts.resize(misMatchSize, 0); }
 558                         for (int j = 0; j < misMatchCounts.size(); j++) {
 559                                 in >> tempNum; misMatchCounts[j] += tempNum;
 560                         }
 561                         m->gobble(in);
 562
 563                         //input megaAlignVector
 564                         string thisLine;
 565                         for (int j = 0; j < megaAlignVector.size(); j++) {
 566                                 thisLine = m->getline(in); m->gobble(in); megaAlignVector[j] += thisLine + '\n';
 567                         }
 568                         m->gobble(in);
 569
 570                         in.close(); remove(tempFile.c_str());
 571
 572                 }
 573 #endif
 574                 return num;
 575         }
 576         catch(exception& e) {
 577                 m->errorOut(e, "SeqErrorCommand", "createProcesses");
 578                 exit(1);
 579         }
 580 }
 581 //**********************************************************************************************************************
 582 int SeqErrorCommand::driver(string filename, string qFileName, string rFileName, string summaryFileName, string errorOutputFileName, string chimeraOutputFileName, linePair line, linePair qline, linePair rline) {
 583
 584         try {
 585                 ReportFile report;
 586                 QualityScores quality;
 587
 588                 misMatchCounts.resize(11, 0);
 589                 int maxMismatch = 0;
 590                 int numSeqs = 0;
 591
 592                 map<string, int>::iterator it;
 593                 qScoreErrorMap['m'].assign(41, 0);
 594                 qScoreErrorMap['s'].assign(41, 0);
 595                 qScoreErrorMap['i'].assign(41, 0);
 596                 qScoreErrorMap['a'].assign(41, 0);
 597
 598                 errorForward['m'].assign(maxLength,0);
 599                 errorForward['s'].assign(maxLength,0);
 600                 errorForward['i'].assign(maxLength,0);
 601                 errorForward['d'].assign(maxLength,0);
 602                 errorForward['a'].assign(maxLength,0);
 603
 604                 errorReverse['m'].assign(maxLength,0);
 605                 errorReverse['s'].assign(maxLength,0);
 606                 errorReverse['i'].assign(maxLength,0);
 607                 errorReverse['d'].assign(maxLength,0);
 608                 errorReverse['a'].assign(maxLength,0);
 609
 610                 //open inputfiles and go to beginning place for this processor
 611                 ifstream queryFile;
 612                 m->openInputFile(filename, queryFile);
 613                 queryFile.seekg(line.start);
 614
 615                 ifstream reportFile;
 616                 ifstream qualFile;
 617                 if(qFileName != "" && rFileName != ""){
 618                         m->openInputFile(qFileName, qualFile);
 619                         qualFile.seekg(qline.start);
 620
 621                         //gobble headers
 622                         if (rline.start == 0) {  report = ReportFile(reportFile, rFileName); }
 623                         else{
 624                                 m->openInputFile(rFileName, reportFile);
 625                                 reportFile.seekg(rline.start);
 626                         }
 627
 628                         qualForwardMap.resize(maxLength);
 629                         qualReverseMap.resize(maxLength);
 630                         for(int i=0;i<maxLength;i++){
 631                                 qualForwardMap[i].assign(41,0);
 632                                 qualReverseMap[i].assign(41,0);
 633                         }
 634                 }
 635
 636                 ofstream outChimeraReport;
 637                 m->openOutputFile(chimeraOutputFileName, outChimeraReport);
 638                 RefChimeraTest chimeraTest(referenceSeqs);
 639                 if (line.start == 0) { chimeraTest.printHeader(outChimeraReport); }
 640
 641                 ofstream errorSummaryFile;
 642                 m->openOutputFile(summaryFileName, errorSummaryFile);
 643                 if (line.start == 0) { printErrorHeader(errorSummaryFile); }
 644
 645                 ofstream errorSeqFile;
 646                 m->openOutputFile(errorOutputFileName, errorSeqFile);
 647
 648                 megaAlignVector.resize(numRefs, "");
 649
 650                 int index = 0;
 651                 bool ignoreSeq = 0;
 652
 653                 bool moreSeqs = 1;
 654                 while (moreSeqs) {
 655
 656                         if (m->control_pressed) { queryFile.close(); if(qFileName != "" && rFileName != ""){  reportFile.close(); qualFile.close(); } outChimeraReport.close(); errorSummaryFile.close();errorSeqFile.close(); return 0; }
 657
 658                         Sequence query(queryFile);
 659
 660                         int numParentSeqs = chimeraTest.analyzeQuery(query.getName(), query.getAligned(), outChimeraReport);
 661                         int closestRefIndex = chimeraTest.getClosestRefIndex();
 662
 663                         if(numParentSeqs > 1 && ignoreChimeras == 1)    {       ignoreSeq = 1;  }
 664                         else                                                                                    {       ignoreSeq = 0;  }
 665
 666                         Compare minCompare = getErrors(query, referenceSeqs[closestRefIndex]);
 667
 668                         if(namesFileName != ""){
 669                                 it = weights.find(query.getName());
 670                                 minCompare.weight = it->second;
 671                         }
 672                         else{   minCompare.weight = 1;  }
 673
 674                         printErrorData(minCompare, numParentSeqs, errorSummaryFile, errorSeqFile);
 675
 676                         if(!ignoreSeq){
 677
 678                                 for(int i=0;i<minCompare.sequence.length();i++){
 679                                         char letter = minCompare.sequence[i];
 680
 681                                         errorForward[letter][i] += minCompare.weight;
 682                                         errorReverse[letter][minCompare.total-i-1] += minCompare.weight;
 683                                 }
 684                         }
 685
 686                         if(qualFileName != "" && reportFileName != ""){
 687                                 report = ReportFile(reportFile);
 688
 689                                 //                              int origLength = report.getQueryLength();
 690                                 int startBase = report.getQueryStart();
 691                                 int endBase = report.getQueryEnd();
 692
 693                                 quality = QualityScores(qualFile);
 694
 695                                 if(!ignoreSeq){
 696                                         quality.updateQScoreErrorMap(qScoreErrorMap, minCompare.sequence, startBase, endBase, minCompare.weight);
 697                                         quality.updateForwardMap(qualForwardMap, startBase, endBase, minCompare.weight);
 698                                         quality.updateReverseMap(qualReverseMap, startBase, endBase, minCompare.weight);
 699                                 }
 700                         }
 701
 702                         if(minCompare.errorRate < threshold && !ignoreSeq){
 703                                 totalBases += (minCompare.total * minCompare.weight);
 704                                 totalMatches += minCompare.matches * minCompare.weight;
 705                                 if(minCompare.mismatches > maxMismatch){
 706                                         maxMismatch = minCompare.mismatches;
 707                                         misMatchCounts.resize(maxMismatch + 1, 0);
 708                                 }
 709                                 misMatchCounts[minCompare.mismatches] += minCompare.weight;
 710                                 numSeqs++;
 711
 712                                 megaAlignVector[closestRefIndex] += query.getInlineSeq() + '\n';
 713                         }
 714
 715                         index++;
 716
 717                         #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 718                                 unsigned long int pos = queryFile.tellg();
 719                                 if ((pos == -1) || (pos >= line.end)) { break; }
 720                         #else
 721                                 if (queryFile.eof()) { break; }
 722                         #endif
 723
 724                         if(index % 100 == 0){   m->mothurOut(toString(index) + '\n');   }
 725                 }
 726                 queryFile.close();
 727                 if(qFileName != "" && rFileName != ""){  reportFile.close(); qualFile.close(); }
 728                 errorSummaryFile.close();
 729                 errorSeqFile.close();
 730
 731                 //report progress
 732                 if(index % 100 != 0){   m->mothurOut(toString(index) + '\n');   }
 733
 734                 return index;
 735         }
 736         catch(exception& e) {
 737                 m->errorOut(e, "SeqErrorCommand", "driver");
 738                 exit(1);
 739         }
 740 }
 741 //***************************************************************************************************************
 742
 743 void SeqErrorCommand::getReferences(){
 744         try {
 745
 746                 ifstream referenceFile;
 747                 m->openInputFile(referenceFileName, referenceFile);
 748
 749                 int numAmbigSeqs = 0;
 750
 751                 int maxStartPos = 0;
 752                 int minEndPos = 100000;
 753
 754                 while(referenceFile){
 755                         Sequence currentSeq(referenceFile);
 756                         int numAmbigs = currentSeq.getAmbigBases();
 757                         if(numAmbigs > 0){      numAmbigSeqs++; }
 758
 759 //                      int startPos = currentSeq.getStartPos();
 760 //                      if(startPos > maxStartPos)      {       maxStartPos = startPos; }
 761 //
 762 //                      int endPos = currentSeq.getEndPos();
 763 //                      if(endPos < minEndPos)          {       minEndPos = endPos;             }
 764                         referenceSeqs.push_back(currentSeq);
 765
 766                         m->gobble(referenceFile);
 767                 }
 768                 referenceFile.close();
 769                 numRefs = referenceSeqs.size();
 770
 771
 772                 for(int i=0;i<numRefs;i++){
 773                         referenceSeqs[i].padToPos(maxStartPos);
 774                         referenceSeqs[i].padFromPos(minEndPos);
 775                 }
 776
 777                 if(numAmbigSeqs != 0){
 778                         m->mothurOut("Warning: " + toString(numAmbigSeqs) + " reference sequences have ambiguous bases, these bases will be ignored\n");
 779                 }
 780
 781         }
 782         catch(exception& e) {
 783                 m->errorOut(e, "SeqErrorCommand", "getReferences");
 784                 exit(1);
 785         }
 786 }
 787
 788 //***************************************************************************************************************
 789
 790 Compare SeqErrorCommand::getErrors(Sequence query, Sequence reference){
 791         try {
 792                 if(query.getAlignLength() != reference.getAlignLength()){
 793                         m->mothurOut("Warning: " + toString(query.getName()) + " and " + toString(reference.getName()) + " are different lengths\n");
 794                 }
 795                 int alignLength = query.getAlignLength();
 796
 797                 string q = query.getAligned();
 798                 string r = reference.getAligned();
 799
 800                 int started = 0;
 801                 Compare errors;
 802
 803                 for(int i=0;i<alignLength;i++){
 804                         if(r[i] != 'N' && q[i] != '.' && r[i] != '.' && (q[i] != '-' || r[i] != '-')){                  //      no missing data and no double gaps
 805                                 started = 1;
 806
 807                                 if(q[i] == 'A'){
 808                                         if(r[i] == 'A'){        errors.AA++;    errors.matches++;       errors.sequence += 'm'; }
 809                                         if(r[i] == 'T'){        errors.AT++;    errors.sequence += 's'; }
 810                                         if(r[i] == 'G'){        errors.AG++;    errors.sequence += 's'; }
 811                                         if(r[i] == 'C'){        errors.AC++;    errors.sequence += 's'; }
 812                                         if(r[i] == '-'){        errors.Ai++;    errors.sequence += 'i'; }
 813                                 }
 814                                 else if(q[i] == 'T'){
 815                                         if(r[i] == 'A'){        errors.TA++;    errors.sequence += 's'; }
 816                                         if(r[i] == 'T'){        errors.TT++;    errors.matches++;       errors.sequence += 'm'; }
 817                                         if(r[i] == 'G'){        errors.TG++;    errors.sequence += 's'; }
 818                                         if(r[i] == 'C'){        errors.TC++;    errors.sequence += 's'; }
 819                                         if(r[i] == '-'){        errors.Ti++;    errors.sequence += 'i'; }
 820                                 }
 821                                 else if(q[i] == 'G'){
 822                                         if(r[i] == 'A'){        errors.GA++;    errors.sequence += 's'; }
 823                                         if(r[i] == 'T'){        errors.GT++;    errors.sequence += 's'; }
 824                                         if(r[i] == 'G'){        errors.GG++;    errors.matches++;       errors.sequence += 'm'; }
 825                                         if(r[i] == 'C'){        errors.GC++;    errors.sequence += 's'; }
 826                                         if(r[i] == '-'){        errors.Gi++;    errors.sequence += 'i'; }
 827                                 }
 828                                 else if(q[i] == 'C'){
 829                                         if(r[i] == 'A'){        errors.CA++;    errors.sequence += 's'; }
 830                                         if(r[i] == 'T'){        errors.CT++;    errors.sequence += 's'; }
 831                                         if(r[i] == 'G'){        errors.CG++;    errors.sequence += 's'; }
 832                                         if(r[i] == 'C'){        errors.CC++;    errors.matches++;       errors.sequence += 'm'; }
 833                                         if(r[i] == '-'){        errors.Ci++;    errors.sequence += 'i'; }
 834                                 }
 835                                 else if(q[i] == 'N'){
 836                                         if(r[i] == 'A'){        errors.NA++;    errors.sequence += 'a'; }
 837                                         if(r[i] == 'T'){        errors.NT++;    errors.sequence += 'a'; }
 838                                         if(r[i] == 'G'){        errors.NG++;    errors.sequence += 'a'; }
 839                                         if(r[i] == 'C'){        errors.NC++;    errors.sequence += 'a'; }
 840                                         if(r[i] == '-'){        errors.Ni++;    errors.sequence += 'a'; }
 841                                 }
 842                                 else if(q[i] == '-' && r[i] != '-'){
 843                                         if(r[i] == 'A'){        errors.dA++;    errors.sequence += 'd'; }
 844                                         if(r[i] == 'T'){        errors.dT++;    errors.sequence += 'd'; }
 845                                         if(r[i] == 'G'){        errors.dG++;    errors.sequence += 'd'; }
 846                                         if(r[i] == 'C'){        errors.dC++;    errors.sequence += 'd'; }
 847                                 }
 848                                 errors.total++;
 849
 850                         }
 851                         else if(q[i] == '.' && r[i] != '.'){            //      reference extends beyond query
 852                                 if(started == 1){       break;  }
 853                         }
 854                         else if(q[i] != '.' && r[i] == '.'){            //      query extends beyond reference
 855                                 if(started == 1){       break;  }
 856                         }
 857                         else if(q[i] == '.' && r[i] == '.'){            //      both are missing data
 858                                 if(started == 1){       break;  }
 859                         }
 860
 861                 }
 862                 errors.mismatches = errors.total-errors.matches;
 863                 errors.errorRate = (double)(errors.total-errors.matches) / (double)errors.total;
 864                 errors.queryName = query.getName();
 865                 errors.refName = reference.getName();
 866
 867                 return errors;
 868         }
 869         catch(exception& e) {
 870                 m->errorOut(e, "SeqErrorCommand", "getErrors");
 871                 exit(1);
 872         }
 873 }
 874
 875 //***************************************************************************************************************
 876
 877 map<string, int> SeqErrorCommand::getWeights(){
 878         ifstream nameFile;
 879         m->openInputFile(namesFileName, nameFile);
 880
 881         string seqName;
 882         string redundantSeqs;
 883         map<string, int> nameCountMap;
 884
 885         while(nameFile){
 886                 nameFile >> seqName >> redundantSeqs;
 887                 nameCountMap[seqName] = m->getNumNames(redundantSeqs);
 888                 m->gobble(nameFile);
 889         }
 890         return nameCountMap;
 891 }
 892
 893
 894 //***************************************************************************************************************
 895
 896 void SeqErrorCommand::printErrorHeader(ofstream& errorSummaryFile){
 897         try {
 898                 errorSummaryFile << "query\treference\tweight\t";
 899                 errorSummaryFile << "AA\tAT\tAG\tAC\tTA\tTT\tTG\tTC\tGA\tGT\tGG\tGC\tCA\tCT\tCG\tCC\tNA\tNT\tNG\tNC\tAi\tTi\tGi\tCi\tNi\tdA\tdT\tdG\tdC\t";
 900                 errorSummaryFile << "insertions\tdeletions\tsubstitutions\tambig\tmatches\tmismatches\ttotal\terror\tnumparents\n";
 901
 902                 errorSummaryFile << setprecision(6);
 903                 errorSummaryFile.setf(ios::fixed);
 904         }
 905         catch(exception& e) {
 906                 m->errorOut(e, "SeqErrorCommand", "printErrorHeader");
 907                 exit(1);
 908         }
 909 }
 910
 911 //***************************************************************************************************************
 912
 913 void SeqErrorCommand::printErrorData(Compare error, int numParentSeqs, ofstream& errorSummaryFile, ofstream& errorSeqFile){
 914         try {
 915
 916                 errorSummaryFile << error.queryName << '\t' << error.refName << '\t' << error.weight << '\t';
 917                 errorSummaryFile << error.AA << '\t' << error.AT << '\t' << error.AG << '\t' << error.AC << '\t';
 918                 errorSummaryFile << error.TA << '\t' << error.TT << '\t' << error.TG << '\t' << error.TC << '\t';
 919                 errorSummaryFile << error.GA << '\t' << error.GT << '\t' << error.GG << '\t' << error.GC << '\t';
 920                 errorSummaryFile << error.CA << '\t' << error.CT << '\t' << error.CG << '\t' << error.CC << '\t';
 921                 errorSummaryFile << error.NA << '\t' << error.NT << '\t' << error.NG << '\t' << error.NC << '\t';
 922                 errorSummaryFile << error.Ai << '\t' << error.Ti << '\t' << error.Gi << '\t' << error.Ci << '\t' << error.Ni << '\t';
 923                 errorSummaryFile << error.dA << '\t' << error.dT << '\t' << error.dG << '\t' << error.dC << '\t';
 924
 925                 errorSummaryFile << error.Ai + error.Ti + error.Gi + error.Ci << '\t';                  //insertions
 926                 errorSummaryFile << error.dA + error.dT + error.dG + error.dC << '\t';                  //deletions
 927                 errorSummaryFile << error.mismatches - (error.Ai + error.Ti + error.Gi + error.Ci) - (error.dA + error.dT + error.dG + error.dC) - (error.NA + error.NT + error.NG + error.NC + error.Ni) << '\t';      //substitutions
 928                 errorSummaryFile << error.NA + error.NT + error.NG + error.NC + error.Ni << '\t';       //ambiguities
 929                 errorSummaryFile << error.matches << '\t' << error.mismatches << '\t' << error.total << '\t' << error.errorRate << '\t' << numParentSeqs << endl;
 930
 931                 errorSeqFile << '>' << error.queryName << "\tref:" << error.refName << '\n' << error.sequence << endl;
 932
 933                 int a=0;                int t=1;                int g=2;                int c=3;
 934                 int gap=4;              int n=5;
 935
 936                 if(numParentSeqs == 1 || ignoreChimeras == 0){
 937                         substitutionMatrix[a][a] += error.weight * error.AA;
 938                         substitutionMatrix[a][t] += error.weight * error.TA;
 939                         substitutionMatrix[a][g] += error.weight * error.GA;
 940                         substitutionMatrix[a][c] += error.weight * error.CA;
 941                         substitutionMatrix[a][gap] += error.weight * error.dA;
 942                         substitutionMatrix[a][n] += error.weight * error.NA;
 943
 944                         substitutionMatrix[t][a] += error.weight * error.AT;
 945                         substitutionMatrix[t][t] += error.weight * error.TT;
 946                         substitutionMatrix[t][g] += error.weight * error.GT;
 947                         substitutionMatrix[t][c] += error.weight * error.CT;
 948                         substitutionMatrix[t][gap] += error.weight * error.dT;
 949                         substitutionMatrix[t][n] += error.weight * error.NT;
 950
 951                         substitutionMatrix[g][a] += error.weight * error.AG;
 952                         substitutionMatrix[g][t] += error.weight * error.TG;
 953                         substitutionMatrix[g][g] += error.weight * error.GG;
 954                         substitutionMatrix[g][c] += error.weight * error.CG;
 955                         substitutionMatrix[g][gap] += error.weight * error.dG;
 956                         substitutionMatrix[g][n] += error.weight * error.NG;
 957
 958                         substitutionMatrix[c][a] += error.weight * error.AC;
 959                         substitutionMatrix[c][t] += error.weight * error.TC;
 960                         substitutionMatrix[c][g] += error.weight * error.GC;
 961                         substitutionMatrix[c][c] += error.weight * error.CC;
 962                         substitutionMatrix[c][gap] += error.weight * error.dC;
 963                         substitutionMatrix[c][n] += error.weight * error.NC;
 964
 965                         substitutionMatrix[gap][a] += error.weight * error.Ai;
 966                         substitutionMatrix[gap][t] += error.weight * error.Ti;
 967                         substitutionMatrix[gap][g] += error.weight * error.Gi;
 968                         substitutionMatrix[gap][c] += error.weight * error.Ci;
 969                         substitutionMatrix[gap][n] += error.weight * error.Ni;
 970                 }
 971         }
 972         catch(exception& e) {
 973                 m->errorOut(e, "SeqErrorCommand", "printErrorData");
 974                 exit(1);
 975         }
 976 }
 977
 978 //***************************************************************************************************************
 979
 980 void SeqErrorCommand::printSubMatrix(){
 981         try {
 982                 string subMatrixFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.matrix";
 983                 ofstream subMatrixFile;
 984                 m->openOutputFile(subMatrixFileName, subMatrixFile);
 985                 outputNames.push_back(subMatrixFileName);  outputTypes["error.matrix"].push_back(subMatrixFileName);
 986                 vector<string> bases(6);
 987                 bases[0] = "A";
 988                 bases[1] = "T";
 989                 bases[2] = "G";
 990                 bases[3] = "C";
 991                 bases[4] = "Gap";
 992                 bases[5] = "N";
 993                 vector<int> refSums(5,1);
 994
 995                 for(int i=0;i<5;i++){
 996                         subMatrixFile << "\tr" << bases[i];
 997
 998                         for(int j=0;j<6;j++){
 999                                 refSums[i] += substitutionMatrix[i][j];
1000                         }
1001                 }
1002                 subMatrixFile << endl;
1003
1004                 for(int i=0;i<6;i++){
1005                         subMatrixFile << 'q' << bases[i];
1006                         for(int j=0;j<5;j++){
1007                                 subMatrixFile << '\t' << substitutionMatrix[j][i];
1008                         }
1009                         subMatrixFile << endl;
1010                 }
1011
1012                 subMatrixFile << "total";
1013                 for(int i=0;i<5;i++){
1014                         subMatrixFile << '\t' << refSums[i];
1015                 }
1016                 subMatrixFile << endl;
1017                 subMatrixFile.close();
1018         }
1019         catch(exception& e) {
1020                 m->errorOut(e, "SeqErrorCommand", "printSubMatrix");
1021                 exit(1);
1022         }
1023 }
1024 //***************************************************************************************************************
1025
1026 void SeqErrorCommand::printErrorFRFile(map<char, vector<int> > errorForward, map<char, vector<int> > errorReverse){
1027         try{
1028                 string errorForwardFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq.forward";
1029                 ofstream errorForwardFile;
1030                 m->openOutputFile(errorForwardFileName, errorForwardFile);
1031                 outputNames.push_back(errorForwardFileName);  outputTypes["error.forward"].push_back(errorForwardFileName);
1032
1033                 errorForwardFile << "position\ttotalseqs\tmatch\tsubstitution\tinsertion\tdeletion\tambiguous" << endl;
1034                 for(int i=0;i<maxLength;i++){
1035                         float match = (float)errorForward['m'][i];
1036                         float subst = (float)errorForward['s'][i];
1037                         float insert = (float)errorForward['i'][i];
1038                         float del = (float)errorForward['d'][i];
1039                         float amb = (float)errorForward['a'][i];
1040                         float total = match + subst + insert + del + amb;
1041                         if(total == 0){ break;  }
1042                         errorForwardFile << i+1 << '\t' << total << '\t' << match/total  << '\t' << subst/total  << '\t' << insert/total  << '\t' << del/total  << '\t' << amb/total << endl;
1043                 }
1044                 errorForwardFile.close();
1045
1046                 string errorReverseFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq.reverse";
1047                 ofstream errorReverseFile;
1048                 m->openOutputFile(errorReverseFileName, errorReverseFile);
1049                 outputNames.push_back(errorReverseFileName);  outputTypes["error.reverse"].push_back(errorReverseFileName);
1050
1051                 errorReverseFile << "position\ttotalseqs\tmatch\tsubstitution\tinsertion\tdeletion\tambiguous" << endl;
1052                 for(int i=0;i<maxLength;i++){
1053                         float match = (float)errorReverse['m'][i];
1054                         float subst = (float)errorReverse['s'][i];
1055                         float insert = (float)errorReverse['i'][i];
1056                         float del = (float)errorReverse['d'][i];
1057                         float amb = (float)errorReverse['a'][i];
1058                         float total = match + subst + insert + del + amb;
1059                         if(total == 0){ break;  }
1060                         errorReverseFile << i+1 << '\t' << total << '\t' << match/total  << '\t' << subst/total  << '\t' << insert/total  << '\t' << del/total  << '\t' << amb/total << endl;
1061                 }
1062                 errorReverseFile.close();
1063         }
1064         catch(exception& e) {
1065                 m->errorOut(e, "SeqErrorCommand", "printErrorFRFile");
1066                 exit(1);
1067         }
1068 }
1069
1070 //***************************************************************************************************************
1071
1072 void SeqErrorCommand::printErrorQuality(map<char, vector<int> > qScoreErrorMap){
1073         try{
1074
1075                 string errorQualityFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.quality";
1076                 ofstream errorQualityFile;
1077                 m->openOutputFile(errorQualityFileName, errorQualityFile);
1078                 outputNames.push_back(errorQualityFileName);  outputTypes["error.quality"].push_back(errorQualityFileName);
1079
1080                 errorQualityFile << "qscore\tmatches\tsubstitutions\tinsertions\tambiguous" << endl;
1081                 for(int i=0;i<41;i++){
1082                         errorQualityFile << i << '\t' << qScoreErrorMap['m'][i] << '\t' << qScoreErrorMap['s'][i] << '\t' << qScoreErrorMap['i'][i] << '\t'<< qScoreErrorMap['a'][i] << endl;
1083                 }
1084                 errorQualityFile.close();
1085         }
1086         catch(exception& e) {
1087                 m->errorOut(e, "SeqErrorCommand", "printErrorFRFile");
1088                 exit(1);
1089         }
1090 }
1091
1092
1093 //***************************************************************************************************************
1094
1095 void SeqErrorCommand::printQualityFR(vector<vector<int> > qualForwardMap, vector<vector<int> > qualReverseMap){
1096
1097         try{
1098                 int numRows = 0;
1099                 int numColumns = qualForwardMap[0].size();
1100
1101                 for(int i=0;i<qualForwardMap.size();i++){
1102                         for(int j=0;j<numColumns;j++){
1103                                 if(qualForwardMap[i][j] != 0){
1104                                         if(numRows < i)         {       numRows = i+20;         }
1105                                 }
1106                         }
1107                 }
1108
1109                 string qualityForwardFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.qual.forward";
1110                 ofstream qualityForwardFile;
1111                 m->openOutputFile(qualityForwardFileName, qualityForwardFile);
1112                 outputNames.push_back(qualityForwardFileName);  outputTypes["error.qual.forward"].push_back(qualityForwardFileName);
1113
1114                 for(int i=0;i<numColumns;i++){  qualityForwardFile << '\t' << i;        }       qualityForwardFile << endl;
1115
1116                 for(int i=0;i<numRows;i++){
1117                         qualityForwardFile << i+1;
1118                         for(int j=0;j<numColumns;j++){
1119                                 qualityForwardFile << '\t' << qualForwardMap[i][j];
1120                         }
1121
1122                         qualityForwardFile << endl;
1123                 }
1124                 qualityForwardFile.close();
1125
1126
1127                 string qualityReverseFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.qual.reverse";
1128                 ofstream qualityReverseFile;
1129                 m->openOutputFile(qualityReverseFileName, qualityReverseFile);
1130                 outputNames.push_back(qualityReverseFileName);  outputTypes["error.qual.reverse"].push_back(qualityReverseFileName);
1131
1132                 for(int i=0;i<numColumns;i++){  qualityReverseFile << '\t' << i;        }       qualityReverseFile << endl;
1133                 for(int i=0;i<numRows;i++){
1134
1135                         qualityReverseFile << i+1;
1136                         for(int j=0;j<numColumns;j++){
1137                                 qualityReverseFile << '\t' << qualReverseMap[i][j];
1138                         }
1139                         qualityReverseFile << endl;
1140                 }
1141                 qualityReverseFile.close();
1142         }
1143         catch(exception& e) {
1144                 m->errorOut(e, "SeqErrorCommand", "printErrorFRFile");
1145                 exit(1);
1146         }
1147
1148 }
1149 /**************************************************************************************************/
1150
1151 int SeqErrorCommand::setLines(string filename, string qfilename, string rfilename, vector<unsigned long int>& fastaFilePos, vector<unsigned long int>& qfileFilePos, vector<unsigned long int>& rfileFilePos) {
1152         try {
1153 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1154                 //set file positions for fasta file
1155                 fastaFilePos = m->divideFile(filename, processors);
1156
1157                 if (qfilename == "") { return processors; }
1158
1159                 //get name of first sequence in each chunk
1160                 map<string, int> firstSeqNames;
1161                 for (int i = 0; i < (fastaFilePos.size()-1); i++) {
1162                         ifstream in;
1163                         m->openInputFile(filename, in);
1164                         in.seekg(fastaFilePos[i]);
1165
1166                         Sequence temp(in);
1167                         firstSeqNames[temp.getName()] = i;
1168
1169                         in.close();
1170                 }
1171
1172                 //make copy to use below
1173                 map<string, int> firstSeqNamesReport = firstSeqNames;
1174
1175                 //seach for filePos of each first name in the qfile and save in qfileFilePos
1176                 ifstream inQual;
1177                 m->openInputFile(qfilename, inQual);
1178
1179                 string input;
1180                 while(!inQual.eof()){
1181                         input = m->getline(inQual);
1182
1183                         if (input.length() != 0) {
1184                                 if(input[0] == '>'){ //this is a sequence name line
1185                                         istringstream nameStream(input);
1186
1187                                         string sname = "";  nameStream >> sname;
1188                                         sname = sname.substr(1);
1189
1190                                         map<string, int>::iterator it = firstSeqNames.find(sname);
1191
1192                                         if(it != firstSeqNames.end()) { //this is the start of a new chunk
1193                                                 unsigned long int pos = inQual.tellg();
1194                                                 qfileFilePos.push_back(pos - input.length() - 1);
1195                                                 firstSeqNames.erase(it);
1196                                         }
1197                                 }
1198                         }
1199
1200                         if (firstSeqNames.size() == 0) { break; }
1201                 }
1202                 inQual.close();
1203
1204                 if (firstSeqNames.size() != 0) {
1205                         for (map<string, int>::iterator it = firstSeqNames.begin(); it != firstSeqNames.end(); it++) {
1206                                 m->mothurOut(it->first + " is in your fasta file and not in your quality file, aborting."); m->mothurOutEndLine();
1207                         }
1208                         m->control_pressed = true;
1209                         return processors;
1210                 }
1211
1212                 //get last file position of qfile
1213                 FILE * pFile;
1214                 unsigned long int size;
1215
1216                 //get num bytes in file
1217                 pFile = fopen (qfilename.c_str(),"rb");
1218                 if (pFile==NULL) perror ("Error opening file");
1219                 else{
1220                         fseek (pFile, 0, SEEK_END);
1221                         size=ftell (pFile);
1222                         fclose (pFile);
1223                 }
1224
1225                 qfileFilePos.push_back(size);
1226
1227                 //seach for filePos of each first name in the rfile and save in rfileFilePos
1228                 string junk;
1229                 ifstream inR;
1230                 m->openInputFile(rfilename, inR);
1231
1232                 //read column headers
1233                 for (int i = 0; i < 16; i++) {
1234                         if (!inR.eof()) {       inR >> junk;    }
1235                         else                    {       break;                  }
1236                 }
1237
1238                 while(!inR.eof()){
1239
1240                         if (m->control_pressed) { inR.close();  return processors; }
1241
1242                         input = m->getline(inR);
1243
1244                         if (input.length() != 0) {
1245
1246                                 istringstream nameStream(input);
1247                                 string sname = "";  nameStream >> sname;
1248
1249                                 map<string, int>::iterator it = firstSeqNamesReport.find(sname);
1250
1251                                 if(it != firstSeqNamesReport.end()) { //this is the start of a new chunk
1252                                         unsigned long int pos = inR.tellg();
1253                                         rfileFilePos.push_back(pos - input.length() - 1);
1254                                         firstSeqNamesReport.erase(it);
1255                                 }
1256                         }
1257
1258                         if (firstSeqNamesReport.size() == 0) { break; }
1259                         m->gobble(inR);
1260                 }
1261                 inR.close();
1262
1263                 if (firstSeqNamesReport.size() != 0) {
1264                         for (map<string, int>::iterator it = firstSeqNamesReport.begin(); it != firstSeqNamesReport.end(); it++) {
1265                                 m->mothurOut(it->first + " is in your fasta file and not in your report file, aborting."); m->mothurOutEndLine();
1266                         }
1267                         m->control_pressed = true;
1268                         return processors;
1269                 }
1270
1271                 //get last file position of qfile
1272                 FILE * rFile;
1273                 unsigned long int sizeR;
1274
1275                 //get num bytes in file
1276                 rFile = fopen (rfilename.c_str(),"rb");
1277                 if (rFile==NULL) perror ("Error opening file");
1278                 else{
1279                         fseek (rFile, 0, SEEK_END);
1280                         sizeR=ftell (rFile);
1281                         fclose (rFile);
1282                 }
1283
1284                 rfileFilePos.push_back(sizeR);
1285
1286                 return processors;
1287
1288 #else
1289
1290                 fastaFilePos.push_back(0); qfileFilePos.push_back(0);
1291                 //get last file position of fastafile
1292                 FILE * pFile;
1293                 unsigned long int size;
1294
1295                 //get num bytes in file
1296                 pFile = fopen (filename.c_str(),"rb");
1297                 if (pFile==NULL) perror ("Error opening file");
1298                 else{
1299                         fseek (pFile, 0, SEEK_END);
1300                         size=ftell (pFile);
1301                         fclose (pFile);
1302                 }
1303                 fastaFilePos.push_back(size);
1304
1305                 //get last file position of fastafile
1306                 FILE * qFile;
1307
1308                 //get num bytes in file
1309                 qFile = fopen (qfilename.c_str(),"rb");
1310                 if (qFile==NULL) perror ("Error opening file");
1311                 else{
1312                         fseek (qFile, 0, SEEK_END);
1313                         size=ftell (qFile);
1314                         fclose (qFile);
1315                 }
1316                 qfileFilePos.push_back(size);
1317
1318                 return 1;
1319
1320 #endif
1321         }
1322         catch(exception& e) {
1323                 m->errorOut(e, "SeqErrorCommand", "setLines");
1324                 exit(1);
1325         }
1326 }
1327 //***************************************************************************************************************