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