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
  14 //**********************************************************************************************************************
  15 vector<string> SeqErrorCommand::getValidParameters(){
  16         try {
  17                 string Array[] =  {"query", "reference", "name", "qfile", "report", "threshold", "inputdir", "outputdir"};
  18                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  19                 return myArray;
  20         }
  21         catch(exception& e) {
  22                 m->errorOut(e, "SeqErrorCommand", "getValidParameters");
  23                 exit(1);
  24         }
  25 }
  26 //**********************************************************************************************************************
  27 SeqErrorCommand::SeqErrorCommand(){
  28         try {
  29                 abort = true;
  30                 //initialize outputTypes
  31                 vector<string> tempOutNames;
  32                 outputTypes["error.summary"] = tempOutNames;
  33                 outputTypes["error.seq"] = tempOutNames;
  34                 outputTypes["error.quality"] = tempOutNames;
  35                 outputTypes["error.qual.forward"] = tempOutNames;
  36                 outputTypes["error.qual.reverse"] = tempOutNames;
  37                 outputTypes["error.forward"] = tempOutNames;
  38                 outputTypes["error.reverse"] = tempOutNames;
  39                 outputTypes["error.count"] = tempOutNames;
  40                 outputTypes["error.matrix"] = tempOutNames;
  41         }
  42         catch(exception& e) {
  43                 m->errorOut(e, "SeqErrorCommand", "SeqErrorCommand");
  44                 exit(1);
  45         }
  46 }
  47 //**********************************************************************************************************************
  48 vector<string> SeqErrorCommand::getRequiredParameters(){
  49         try {
  50                 string Array[] =  {"query","reference"};
  51                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  52                 return myArray;
  53         }
  54         catch(exception& e) {
  55                 m->errorOut(e, "SeqErrorCommand", "getRequiredParameters");
  56                 exit(1);
  57         }
  58 }
  59 //**********************************************************************************************************************
  60 vector<string> SeqErrorCommand::getRequiredFiles(){
  61         try {
  62                 vector<string> myArray;
  63                 return myArray;
  64         }
  65         catch(exception& e) {
  66                 m->errorOut(e, "SeqErrorCommand", "getRequiredFiles");
  67                 exit(1);
  68         }
  69 }
  70 //***************************************************************************************************************
  71
  72 SeqErrorCommand::SeqErrorCommand(string option)  {
  73         try {
  74
  75                 abort = false;
  76
  77                 //allow user to run help
  78                 if(option == "help") { help(); abort = true; }
  79
  80                 else {
  81                         string temp;
  82
  83                         //valid paramters for this command
  84                         string AlignArray[] =  {"query", "reference", "name", "qfile", "report", "threshold", "inputdir", "outputdir"};
  85
  86                         vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
  87
  88                         OptionParser parser(option);
  89                         map<string,string> parameters = parser.getParameters();
  90
  91                         ValidParameters validParameter;
  92                         map<string,string>::iterator it;
  93
  94                         //check to make sure all parameters are valid for command
  95                         for (it = parameters.begin(); it != parameters.end(); it++) {
  96                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  97                         }
  98
  99                         //initialize outputTypes
 100                         vector<string> tempOutNames;
 101                         outputTypes["error.summary"] = tempOutNames;
 102                         outputTypes["error.seq"] = tempOutNames;
 103                         outputTypes["error.quality"] = tempOutNames;
 104                         outputTypes["error.qual.forward"] = tempOutNames;
 105                         outputTypes["error.qual.reverse"] = tempOutNames;
 106                         outputTypes["error.forward"] = tempOutNames;
 107                         outputTypes["error.reverse"] = tempOutNames;
 108                         outputTypes["error.count"] = tempOutNames;
 109                         outputTypes["error.matrix"] = tempOutNames;
 110
 111
 112                         //if the user changes the input directory command factory will send this info to us in the output parameter
 113                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 114                         if (inputDir == "not found"){   inputDir = "";          }
 115                         else {
 116                                 string path;
 117                                 it = parameters.find("query");
 118                                 //user has given a template file
 119                                 if(it != parameters.end()){
 120                                         path = m->hasPath(it->second);
 121                                         //if the user has not given a path then, add inputdir. else leave path alone.
 122                                         if (path == "") {       parameters["query"] = inputDir + it->second;            }
 123                                 }
 124
 125                                 it = parameters.find("reference");
 126                                 //user has given a template file
 127                                 if(it != parameters.end()){
 128                                         path = m->hasPath(it->second);
 129                                         //if the user has not given a path then, add inputdir. else leave path alone.
 130                                         if (path == "") {       parameters["reference"] = inputDir + it->second;                }
 131                                 }
 132
 133                                 it = parameters.find("name");
 134                                 //user has given a names file
 135                                 if(it != parameters.end()){
 136                                         path = m->hasPath(it->second);
 137                                         //if the user has not given a path then, add inputdir. else leave path alone.
 138                                         if (path == "") {       parameters["name"] = inputDir + it->second;             }
 139                                 }
 140
 141                                 it = parameters.find("qfile");
 142                                 //user has given a quality score file
 143                                 if(it != parameters.end()){
 144                                         path = m->hasPath(it->second);
 145                                         //if the user has not given a path then, add inputdir. else leave path alone.
 146                                         if (path == "") {       parameters["qfile"] = inputDir + it->second;            }
 147                                 }
 148
 149                                 it = parameters.find("report");
 150                                 //user has given a alignment report file
 151                                 if(it != parameters.end()){
 152                                         path = m->hasPath(it->second);
 153                                         //if the user has not given a path then, add inputdir. else leave path alone.
 154                                         if (path == "") {       parameters["report"] = inputDir + it->second;           }
 155                                 }
 156
 157                         }
 158                         //check for required parameters
 159                         queryFileName = validParameter.validFile(parameters, "query", true);
 160                         if (queryFileName == "not found") { m->mothurOut("query is a required parameter for the seq.error command."); m->mothurOutEndLine(); abort = true; }
 161                         else if (queryFileName == "not open") { abort = true; }
 162
 163                         referenceFileName = validParameter.validFile(parameters, "reference", true);
 164                         if (referenceFileName == "not found") { m->mothurOut("reference is a required parameter for the seq.error command."); m->mothurOutEndLine(); abort = true; }
 165                         else if (referenceFileName == "not open") { abort = true; }
 166
 167
 168                         //check for optional parameters
 169                         namesFileName = validParameter.validFile(parameters, "name", true);
 170                         if(namesFileName == "not found"){       namesFileName = "";     }
 171
 172                         qualFileName = validParameter.validFile(parameters, "qfile", true);
 173                         if(qualFileName == "not found"){        qualFileName = "";      }
 174
 175                         reportFileName = validParameter.validFile(parameters, "report", true);
 176                         if(reportFileName == "not found"){      reportFileName = "";    }
 177
 178                         if((reportFileName != "" && qualFileName == "") || (reportFileName == "" && qualFileName != "")){
 179                                 m->mothurOut("if you use either a qual file or a report file, you have to have both.");
 180                                 m->mothurOutEndLine();
 181                                 abort = true;
 182                         }
 183
 184                         outputDir = validParameter.validFile(parameters, "outputdir", false);
 185                         if (outputDir == "not found"){
 186                                 outputDir = "";
 187                                 outputDir += m->hasPath(queryFileName); //if user entered a file with a path then preserve it
 188                         }
 189
 190                         //check for optional parameter and set defaults
 191                         // ...at some point should added some additional type checking...
 192                         temp = validParameter.validFile(parameters, "threshold", false);        if (temp == "not found") { temp = "1.00"; }
 193                         convert(temp, threshold);
 194
 195                         substitutionMatrix.resize(6);
 196                         for(int i=0;i<6;i++){   substitutionMatrix[i].assign(6,0);      }
 197                 }
 198         }
 199         catch(exception& e) {
 200                 m->errorOut(e, "SeqErrorCommand", "SeqErrorCommand");
 201                 exit(1);
 202         }
 203 }
 204
 205 //**********************************************************************************************************************
 206
 207 void SeqErrorCommand::help(){
 208         try {
 209                 m->mothurOut("The seq.error command reads a query alignment file and a reference alignment file and creates .....\n");
 210                 m->mothurOut("Example seq.error(...).\n");
 211                 m->mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n");
 212                 m->mothurOut("For more details please check out the wiki http://www.mothur.org/wiki/seq.error .\n\n");
 213         }
 214         catch(exception& e) {
 215                 m->errorOut(e, "SeqErrorCommand", "help");
 216                 exit(1);
 217         }
 218 }
 219
 220 //***************************************************************************************************************
 221
 222 SeqErrorCommand::~SeqErrorCommand(){
 223
 224 }
 225
 226 //***************************************************************************************************************
 227
 228 int SeqErrorCommand::execute(){
 229         try{
 230                 if (abort == true) { return 0; }
 231
 232                 errorSummaryFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.summary";
 233                 m->openOutputFile(errorSummaryFileName, errorSummaryFile);
 234                 outputNames.push_back(errorSummaryFileName); outputTypes["error.summary"].push_back(errorSummaryFileName);
 235                 printErrorHeader();
 236
 237                 errorSeqFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq";
 238                 m->openOutputFile(errorSeqFileName, errorSeqFile);
 239                 outputNames.push_back(errorSeqFileName); outputTypes["error.seq"].push_back(errorSeqFileName);
 240                 printErrorHeader();
 241
 242                 getReferences();        //read in reference sequences - make sure there's no ambiguous bases
 243
 244                 map<string, int> weights;
 245                 if(namesFileName != ""){        weights = getWeights(); }
 246
 247                 ifstream queryFile;
 248                 m->openInputFile(queryFileName, queryFile);
 249
 250                 ifstream reportFile;
 251                 ifstream qualFile;
 252
 253                 ReportFile report;
 254                 QualityScores quality;
 255                 vector<vector<int> > qualForwardMap;
 256                 vector<vector<int> > qualReverseMap;
 257
 258                 if(qualFileName != "" && reportFileName != ""){
 259                         m->openInputFile(qualFileName, qualFile);
 260                         report = ReportFile(reportFile, reportFileName);
 261
 262                         qualForwardMap.resize(1000);
 263                         qualReverseMap.resize(1000);
 264                         for(int i=0;i<1000;i++){
 265                                 qualForwardMap[i].assign(100,0);
 266                                 qualReverseMap[i].assign(100,0);
 267                         }
 268                 }
 269
 270                 int totalBases = 0;
 271                 int totalMatches = 0;
 272
 273                 vector<int> misMatchCounts(11, 0);
 274                 int maxMismatch = 0;
 275                 int numSeqs = 0;
 276
 277                 map<string, int>::iterator it;
 278                 map<char, vector<int> > qScoreErrorMap;
 279                 qScoreErrorMap['m'].assign(41, 0);
 280                 qScoreErrorMap['s'].assign(41, 0);
 281                 qScoreErrorMap['i'].assign(41, 0);
 282                 qScoreErrorMap['a'].assign(41, 0);
 283
 284
 285
 286                 map<char, vector<int> > errorForward;
 287                 errorForward['m'].assign(1000,0);
 288                 errorForward['s'].assign(1000,0);
 289                 errorForward['i'].assign(1000,0);
 290                 errorForward['d'].assign(1000,0);
 291                 errorForward['a'].assign(1000,0);
 292
 293                 map<char, vector<int> > errorReverse;
 294                 errorReverse['m'].assign(1000,0);
 295                 errorReverse['s'].assign(1000,0);
 296                 errorReverse['i'].assign(1000,0);
 297                 errorReverse['d'].assign(1000,0);
 298                 errorReverse['a'].assign(1000,0);
 299
 300
 301
 302                 while(queryFile){
 303
 304                         if (m->control_pressed) { errorSummaryFile.close();     errorSeqFile.close(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
 305
 306                         Compare minCompare;
 307                         Sequence query(queryFile);
 308
 309                         for(int i=0;i<numRefs;i++){
 310                                 Compare currCompare = getErrors(query, referenceSeqs[i]);
 311
 312                                 if(currCompare.errorRate < minCompare.errorRate){
 313                                         minCompare = currCompare;
 314                                 }
 315                         }
 316
 317                         if(namesFileName != ""){
 318                                 it = weights.find(query.getName());
 319                                 minCompare.weight = it->second;
 320                         }
 321                         else    {       minCompare.weight = 1;  }
 322
 323                         printErrorData(minCompare);
 324
 325                         for(int i=0;i<minCompare.total;i++){
 326                                 char letter = minCompare.sequence[i];
 327                                 errorForward[letter][i] += minCompare.weight;
 328                                 errorReverse[letter][minCompare.total-i-1] += minCompare.weight;
 329                         }
 330
 331                         if(qualFileName != "" && reportFileName != ""){
 332                                 report = ReportFile(reportFile);
 333
 334                                 int origLength = report.getQueryLength();
 335                                 int startBase = report.getQueryStart();
 336                                 int endBase = report.getQueryEnd();
 337
 338                                 quality = QualityScores(qualFile, origLength);
 339                                 quality.updateQScoreErrorMap(qScoreErrorMap, minCompare.sequence, startBase, endBase, minCompare.weight);
 340                                 quality.updateForwardMap(qualForwardMap, startBase, endBase, minCompare.weight);
 341                                 quality.updateReverseMap(qualReverseMap, startBase, endBase, minCompare.weight);
 342                         }
 343
 344                         if(minCompare.errorRate < threshold){
 345                                 totalBases += (minCompare.total * minCompare.weight);
 346                                 totalMatches += minCompare.matches * minCompare.weight;
 347                                 if(minCompare.mismatches > maxMismatch){
 348                                         maxMismatch = minCompare.mismatches;
 349                                         misMatchCounts.resize(maxMismatch + 1, 0);
 350                                 }
 351                                 misMatchCounts[minCompare.mismatches] += minCompare.weight;
 352                                 numSeqs++;
 353                         }
 354
 355
 356                 }
 357                 queryFile.close();
 358
 359                 int total = 0;
 360
 361                 if(qualFileName != "" && reportFileName != ""){
 362                         string errorQualityFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.quality";
 363                         ofstream errorQualityFile;
 364                         m->openOutputFile(errorQualityFileName, errorQualityFile);
 365                         outputNames.push_back(errorQualityFileName);  outputTypes["error.quality"].push_back(errorQualityFileName);
 366
 367                         errorQualityFile << "qscore\tmatches\tsubstitutions\tinsertions\tambiguous" << endl;
 368                         for(int i=0;i<41;i++){
 369                                 errorQualityFile << i << '\t' << qScoreErrorMap['m'][i] << '\t' << qScoreErrorMap['s'][i] << '\t' << qScoreErrorMap['i'][i] << '\t'<< qScoreErrorMap['a'][i] << endl;
 370                         }
 371                         errorQualityFile.close();
 372
 373
 374
 375                         int lastRow = 0;
 376                         int lastColumn = 0;
 377
 378                         for(int i=0;i<qualForwardMap.size();i++){
 379                                 for(int j=0;j<qualForwardMap[i].size();j++){
 380                                         if(qualForwardMap[i][j] != 0){
 381                                                 if(lastRow < i)         {       lastRow = i+2;  }
 382                                                 if(lastColumn < j)      {       lastColumn = j+2;       }
 383                                         }
 384
 385                                 }
 386                         }
 387
 388
 389                         string qualityForwardFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.qual.forward";
 390                         ofstream qualityForwardFile;
 391                         m->openOutputFile(qualityForwardFileName, qualityForwardFile);
 392                         outputNames.push_back(errorQualityFileName);  outputTypes["error.qual.forward"].push_back(qualityForwardFileName);
 393
 394                         for(int i=0;i<lastColumn;i++){  qualityForwardFile << '\t' << i;        }       qualityForwardFile << endl;
 395                         for(int i=0;i<lastRow;i++){
 396
 397                                 if (m->control_pressed) { qualityForwardFile.close(); errorSummaryFile.close(); errorSeqFile.close(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
 398
 399                                 qualityForwardFile << i+1;
 400                                 for(int j=0;j<lastColumn;j++){
 401                                         qualityForwardFile << '\t' << qualForwardMap[i][j];
 402                                 }
 403                                 qualityForwardFile << endl;
 404                         }
 405                         qualityForwardFile.close();
 406
 407
 408                         string qualityReverseFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.qual.reverse";
 409                         ofstream qualityReverseFile;
 410                         m->openOutputFile(qualityReverseFileName, qualityReverseFile);
 411                         outputNames.push_back(errorQualityFileName);  outputTypes["error.qual.reverse"].push_back(qualityReverseFileName);
 412
 413                         for(int i=0;i<lastColumn;i++){  qualityReverseFile << '\t' << i;        }       qualityReverseFile << endl;
 414                         for(int i=0;i<lastRow;i++){
 415
 416                                 if (m->control_pressed) { qualityReverseFile.close(); errorSummaryFile.close(); errorSeqFile.close(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
 417
 418                                 qualityReverseFile << i+1;
 419                                 for(int j=0;j<lastColumn;j++){
 420                                         qualityReverseFile << '\t' << qualReverseMap[i][j];
 421                                 }
 422                                 qualityReverseFile << endl;
 423                         }
 424                         qualityReverseFile.close();
 425                 }
 426
 427
 428                 string errorForwardFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq.forward";
 429                 ofstream errorForwardFile;
 430                 m->openOutputFile(errorForwardFileName, errorForwardFile);
 431                 outputNames.push_back(errorForwardFileName);  outputTypes["error.forward"].push_back(errorForwardFileName);
 432
 433                 errorForwardFile << "position\ttotalseqs\tmatch\tsubstitution\tinsertion\tdeletion\tambiguous" << endl;
 434                 for(int i=0;i<1000;i++){
 435                         float match = (float)errorForward['m'][i];
 436                         float subst = (float)errorForward['s'][i];
 437                         float insert = (float)errorForward['i'][i];
 438                         float del = (float)errorForward['d'][i];
 439                         float amb = (float)errorForward['a'][i];
 440                         float total = match + subst + insert + del + amb;
 441                         if(total == 0){ break;  }
 442                         errorForwardFile << i+1 << '\t' << total << '\t' << match/total  << '\t' << subst/total  << '\t' << insert/total  << '\t' << del/total  << '\t' << amb/total << endl;
 443                 }
 444                 errorForwardFile.close();
 445
 446                 if (m->control_pressed) { errorSummaryFile.close();     errorSeqFile.close(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
 447
 448                 string errorReverseFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq.reverse";
 449                 ofstream errorReverseFile;
 450                 m->openOutputFile(errorReverseFileName, errorReverseFile);
 451                 outputNames.push_back(errorReverseFileName);  outputTypes["error.reverse"].push_back(errorReverseFileName);
 452
 453                 errorReverseFile << "position\ttotalseqs\tmatch\tsubstitution\tinsertion\tdeletion\tambiguous" << endl;
 454                 for(int i=0;i<1000;i++){
 455                         float match = (float)errorReverse['m'][i];
 456                         float subst = (float)errorReverse['s'][i];
 457                         float insert = (float)errorReverse['i'][i];
 458                         float del = (float)errorReverse['d'][i];
 459                         float amb = (float)errorReverse['a'][i];
 460                         float total = match + subst + insert + del + amb;
 461                         if(total == 0){ break;  }
 462                         errorReverseFile << i+1 << '\t' << total << '\t' << match/total  << '\t' << subst/total  << '\t' << insert/total  << '\t' << del/total  << '\t' << amb/total << endl;
 463                 }
 464                 errorReverseFile.close();
 465
 466                 if (m->control_pressed) { errorSummaryFile.close();     errorSeqFile.close(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
 467
 468                 string errorCountFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.count";
 469                 ofstream errorCountFile;
 470                 m->openOutputFile(errorCountFileName, errorCountFile);
 471                 outputNames.push_back(errorCountFileName);  outputTypes["error.count"].push_back(errorCountFileName);
 472                 m->mothurOut("Overall error rate:\t" + toString((double)(totalBases - totalMatches) / (double)totalBases) + "\n\n");
 473                 m->mothurOut("Errors\tSequences\n");
 474                 errorCountFile << "Errors\tSequences\n";
 475                 for(int i=0;i<misMatchCounts.size();i++){
 476                         m->mothurOut(toString(i) + '\t' + toString(misMatchCounts[i]) + '\n');
 477                         errorCountFile << i << '\t' << misMatchCounts[i] << endl;
 478                 }
 479                 errorCountFile.close();
 480
 481                 if (m->control_pressed) { errorSummaryFile.close();     errorSeqFile.close(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
 482
 483                 string subMatrixFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.matrix";
 484                 ofstream subMatrixFile;
 485                 m->openOutputFile(subMatrixFileName, subMatrixFile);
 486                 outputNames.push_back(subMatrixFileName);  outputTypes["error.matrix"].push_back(subMatrixFileName);
 487                 vector<string> bases(6);
 488                 bases[0] = "A";
 489                 bases[1] = "T";
 490                 bases[2] = "G";
 491                 bases[3] = "C";
 492                 bases[4] = "Gap";
 493                 bases[5] = "N";
 494                 vector<int> refSums(5,1);
 495
 496                 for(int i=0;i<5;i++){
 497                         subMatrixFile << "\tr" << bases[i];
 498
 499                         for(int j=0;j<6;j++){
 500                                 refSums[i] += substitutionMatrix[i][j];
 501                         }
 502
 503                 }
 504                 subMatrixFile << endl;
 505
 506                 for(int i=0;i<6;i++){
 507                         subMatrixFile << 'q' << bases[i];
 508                         for(int j=0;j<5;j++){
 509                                 subMatrixFile << '\t' << substitutionMatrix[j][i];
 510                         }
 511                         subMatrixFile << endl;
 512                 }
 513                 subMatrixFile << "total";
 514                 for(int i=0;i<5;i++){
 515                         subMatrixFile << '\t' << refSums[i];
 516                 }
 517                 subMatrixFile << endl;
 518                 subMatrixFile.close();
 519
 520                 errorSummaryFile.close();
 521                 errorSeqFile.close();
 522
 523                 if (m->control_pressed) {  for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
 524
 525                 m->mothurOutEndLine();
 526                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 527                 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
 528                 m->mothurOutEndLine();
 529
 530                 return 0;
 531         }
 532         catch(exception& e) {
 533                 m->errorOut(e, "SeqErrorCommand", "execute");
 534                 exit(1);
 535         }
 536 }
 537
 538 //***************************************************************************************************************
 539
 540 void SeqErrorCommand::getReferences(){
 541         try {
 542
 543                 ifstream referenceFile;
 544                 m->openInputFile(referenceFileName, referenceFile);
 545
 546                 while(referenceFile){
 547                         Sequence currentSeq(referenceFile);
 548                         int numAmbigs = currentSeq.getAmbigBases();
 549
 550                         if(numAmbigs != 0){
 551                                 m->mothurOut("Warning: " + toString(currentSeq.getName()) + " has " + toString(numAmbigs) + " ambiguous bases, these bases will be removed\n");
 552                                 currentSeq.removeAmbigBases();
 553                         }
 554                         referenceSeqs.push_back(currentSeq);
 555                         m->gobble(referenceFile);
 556                 }
 557                 numRefs = referenceSeqs.size();
 558
 559                 referenceFile.close();
 560         }
 561         catch(exception& e) {
 562                 m->errorOut(e, "SeqErrorCommand", "getReferences");
 563                 exit(1);
 564         }
 565 }
 566
 567 //***************************************************************************************************************
 568
 569 Compare SeqErrorCommand::getErrors(Sequence query, Sequence reference){
 570         try {
 571                 if(query.getAlignLength() != reference.getAlignLength()){
 572                         m->mothurOut("Warning: " + toString(query.getName()) + " and " + toString(reference.getName()) + " are different lengths\n");
 573                 }
 574                 int alignLength = query.getAlignLength();
 575
 576                 string q = query.getAligned();
 577                 string r = reference.getAligned();
 578
 579                 int started = 0;
 580                 Compare errors;
 581
 582                 for(int i=0;i<alignLength;i++){
 583                         if(q[i] != '.' && r[i] != '.' && (q[i] != '-' || r[i] != '-')){                 //      no missing data and no double gaps
 584                                 started = 1;
 585
 586                                 if(q[i] == 'A'){
 587                                         if(r[i] == 'A'){        errors.AA++;    errors.matches++;       errors.sequence += 'm'; }
 588                                         if(r[i] == 'T'){        errors.AT++;    errors.sequence += 's'; }
 589                                         if(r[i] == 'G'){        errors.AG++;    errors.sequence += 's'; }
 590                                         if(r[i] == 'C'){        errors.AC++;    errors.sequence += 's'; }
 591                                         if(r[i] == '-'){        errors.Ai++;    errors.sequence += 'i'; }
 592                                 }
 593                                 else if(q[i] == 'T'){
 594                                         if(r[i] == 'A'){        errors.TA++;    errors.sequence += 's'; }
 595                                         if(r[i] == 'T'){        errors.TT++;    errors.matches++;       errors.sequence += 'm'; }
 596                                         if(r[i] == 'G'){        errors.TG++;    errors.sequence += 's'; }
 597                                         if(r[i] == 'C'){        errors.TC++;    errors.sequence += 's'; }
 598                                         if(r[i] == '-'){        errors.Ti++;    errors.sequence += 'i'; }
 599                                 }
 600                                 else if(q[i] == 'G'){
 601                                         if(r[i] == 'A'){        errors.GA++;    errors.sequence += 's'; }
 602                                         if(r[i] == 'T'){        errors.GT++;    errors.sequence += 's'; }
 603                                         if(r[i] == 'G'){        errors.GG++;    errors.matches++;       errors.sequence += 'm'; }
 604                                         if(r[i] == 'C'){        errors.GC++;    errors.sequence += 's'; }
 605                                         if(r[i] == '-'){        errors.Gi++;    errors.sequence += 'i'; }
 606                                 }
 607                                 else if(q[i] == 'C'){
 608                                         if(r[i] == 'A'){        errors.CA++;    errors.sequence += 's'; }
 609                                         if(r[i] == 'T'){        errors.CT++;    errors.sequence += 's'; }
 610                                         if(r[i] == 'G'){        errors.CG++;    errors.sequence += 's'; }
 611                                         if(r[i] == 'C'){        errors.CC++;    errors.matches++;       errors.sequence += 'm'; }
 612                                         if(r[i] == '-'){        errors.Ci++;    errors.sequence += 'i'; }
 613                                 }
 614                                 else if(q[i] == 'N'){
 615                                         if(r[i] == 'A'){        errors.NA++;    errors.sequence += 'a'; }
 616                                         if(r[i] == 'T'){        errors.NT++;    errors.sequence += 'a'; }
 617                                         if(r[i] == 'G'){        errors.NG++;    errors.sequence += 'a'; }
 618                                         if(r[i] == 'C'){        errors.NC++;    errors.sequence += 'a'; }
 619                                         if(r[i] == '-'){        errors.Ni++;    errors.sequence += 'a'; }
 620                                 }
 621                                 else if(q[i] == '-' && r[i] != '-'){
 622                                         if(r[i] == 'A'){        errors.dA++;    errors.sequence += 'd'; }
 623                                         if(r[i] == 'T'){        errors.dT++;    errors.sequence += 'd'; }
 624                                         if(r[i] == 'G'){        errors.dG++;    errors.sequence += 'd'; }
 625                                         if(r[i] == 'C'){        errors.dC++;    errors.sequence += 'd'; }
 626                                 }
 627                                 errors.total++;
 628
 629                         }
 630                         else if(q[i] == '.' && r[i] != '.'){            //      reference extends beyond query
 631                                 if(started == 1){       break;  }
 632                         }
 633                         else if(q[i] != '.' && r[i] == '.'){            //      query extends beyond reference
 634                                 m->mothurOut("Warning: " + toString(query.getName()) + " extend beyond " + toString(reference.getName()) + ".  Ignoring the extra bases in the query\n");
 635                                 if(started == 1){       break;  }
 636                         }
 637                         else if(q[i] == '.' && r[i] == '.'){            //      both are missing data
 638                                 if(started == 1){       break;  }
 639                         }
 640
 641                 }
 642                 errors.mismatches = errors.total-errors.matches;
 643                 errors.errorRate = (double)(errors.total-errors.matches) / (double)errors.total;
 644                 errors.queryName = query.getName();
 645                 errors.refName = reference.getName();
 646
 647                 return errors;
 648         }
 649         catch(exception& e) {
 650                 m->errorOut(e, "SeqErrorCommand", "getErrors");
 651                 exit(1);
 652         }
 653 }
 654
 655 //***************************************************************************************************************
 656
 657 map<string, int> SeqErrorCommand::getWeights(){
 658         ifstream nameFile;
 659         m->openInputFile(namesFileName, nameFile);
 660
 661         string seqName;
 662         string redundantSeqs;
 663         map<string, int> nameCountMap;
 664
 665         while(nameFile){
 666                 nameFile >> seqName >> redundantSeqs;
 667                 nameCountMap[seqName] = m->getNumNames(redundantSeqs);
 668                 m->gobble(nameFile);
 669         }
 670         return nameCountMap;
 671 }
 672
 673
 674 //***************************************************************************************************************
 675
 676 void SeqErrorCommand::printErrorHeader(){
 677         try {
 678                 errorSummaryFile << "query\treference\tweight\t";
 679                 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";
 680                 errorSummaryFile << "insertions\tdeletions\tsubstitutions\tambig\tmatches\tmismatches\ttotal\terror\n";
 681
 682                 errorSummaryFile << setprecision(6);
 683                 errorSummaryFile.setf(ios::fixed);
 684         }
 685         catch(exception& e) {
 686                 m->errorOut(e, "SeqErrorCommand", "printErrorHeader");
 687                 exit(1);
 688         }
 689 }
 690
 691 //***************************************************************************************************************
 692
 693 void SeqErrorCommand::printErrorData(Compare error){
 694         try {
 695                 errorSummaryFile << error.queryName << '\t' << error.refName << '\t' << error.weight << '\t';
 696                 errorSummaryFile << error.AA << '\t' << error.AT << '\t' << error.AG << '\t' << error.AC << '\t';
 697                 errorSummaryFile << error.TA << '\t' << error.TT << '\t' << error.TG << '\t' << error.TC << '\t';
 698                 errorSummaryFile << error.GA << '\t' << error.GT << '\t' << error.GG << '\t' << error.GC << '\t';
 699                 errorSummaryFile << error.CA << '\t' << error.CT << '\t' << error.CG << '\t' << error.CC << '\t';
 700                 errorSummaryFile << error.NA << '\t' << error.NT << '\t' << error.NG << '\t' << error.NC << '\t';
 701                 errorSummaryFile << error.Ai << '\t' << error.Ti << '\t' << error.Gi << '\t' << error.Ci << '\t' << error.Ni << '\t' ;
 702                 errorSummaryFile << error.dA << '\t' << error.dT << '\t' << error.dG << '\t' << error.dC << '\t';
 703
 704                 errorSummaryFile << error.Ai + error.Ti + error.Gi + error.Ci << '\t';                  //insertions
 705                 errorSummaryFile << error.dA + error.dT + error.dG + error.dC << '\t';                  //deletions
 706                 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
 707                 errorSummaryFile << error.NA + error.NT + error.NG + error.NC + error.Ni << '\t';       //ambiguities
 708                 errorSummaryFile << error.matches << '\t' << error.mismatches << '\t' << error.total << '\t' << error.errorRate << endl;
 709
 710                 errorSeqFile << '>' << error.queryName << "\tref:" << error.refName << '\n' << error.sequence << endl;
 711
 712
 713                 int a=0;                int t=1;                int g=2;                int c=3;
 714                 int gap=4;              int n=5;
 715
 716                 substitutionMatrix[a][a] += error.weight * error.AA;
 717                 substitutionMatrix[a][t] += error.weight * error.TA;
 718                 substitutionMatrix[a][g] += error.weight * error.GA;
 719                 substitutionMatrix[a][c] += error.weight * error.CA;
 720                 substitutionMatrix[a][gap] += error.weight * error.dA;
 721                 substitutionMatrix[a][n] += error.weight * error.NA;
 722
 723                 substitutionMatrix[t][a] += error.weight * error.AT;
 724                 substitutionMatrix[t][t] += error.weight * error.TT;
 725                 substitutionMatrix[t][g] += error.weight * error.GT;
 726                 substitutionMatrix[t][c] += error.weight * error.CT;
 727                 substitutionMatrix[t][gap] += error.weight * error.dT;
 728                 substitutionMatrix[t][n] += error.weight * error.NT;
 729
 730                 substitutionMatrix[g][a] += error.weight * error.AG;
 731                 substitutionMatrix[g][t] += error.weight * error.TG;
 732                 substitutionMatrix[g][g] += error.weight * error.GG;
 733                 substitutionMatrix[g][c] += error.weight * error.CG;
 734                 substitutionMatrix[g][gap] += error.weight * error.dG;
 735                 substitutionMatrix[g][n] += error.weight * error.NG;
 736
 737                 substitutionMatrix[c][a] += error.weight * error.AC;
 738                 substitutionMatrix[c][t] += error.weight * error.TC;
 739                 substitutionMatrix[c][g] += error.weight * error.GC;
 740                 substitutionMatrix[c][c] += error.weight * error.CC;
 741                 substitutionMatrix[c][gap] += error.weight * error.dC;
 742                 substitutionMatrix[c][n] += error.weight * error.NC;
 743
 744                 substitutionMatrix[gap][a] += error.weight * error.Ai;
 745                 substitutionMatrix[gap][t] += error.weight * error.Ti;
 746                 substitutionMatrix[gap][g] += error.weight * error.Gi;
 747                 substitutionMatrix[gap][c] += error.weight * error.Ci;
 748                 substitutionMatrix[gap][n] += error.weight * error.Ni;
 749
 750         }
 751         catch(exception& e) {
 752                 m->errorOut(e, "SeqErrorCommand", "printErrorData");
 753                 exit(1);
 754         }
 755 }
 756
 757 //***************************************************************************************************************
 758
 759
 760
 761
 762
 763
 764