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