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