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