sensspeccommand.cpp

   1 /*
   2  *  sensspeccommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by Pat Schloss on 7/6/10.
   6  *  Copyright 2010 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "sensspeccommand.h"
  11
  12 //**********************************************************************************************************************
  13 vector<string> SensSpecCommand::setParameters(){
  14         try {
  15                 CommandParameter plist("list", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(plist);
  16                 CommandParameter pphylip("phylip", "InputTypes", "", "", "PhylipColumn", "PhylipColumn", "none",false,false); parameters.push_back(pphylip);
  17                 //CommandParameter pname("name", "InputTypes", "", "", "none", "none", "ColumnName",false,false); parameters.push_back(pname);
  18                 CommandParameter pcolumn("column", "InputTypes", "", "", "PhylipColumn", "PhylipColumn", "none",false,false); parameters.push_back(pcolumn);
  19                 CommandParameter plabel("label", "String", "", "", "", "", "",false,false); parameters.push_back(plabel);
  20                 CommandParameter pcutoff("cutoff", "Number", "", "-1.00", "", "", "",false,false); parameters.push_back(pcutoff);
  21                 CommandParameter pprecision("precision", "Number", "", "100", "", "", "",false,false); parameters.push_back(pprecision);
  22                 CommandParameter phard("hard", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(phard);
  23                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
  24                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
  25
  26                 vector<string> myArray;
  27                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  28                 return myArray;
  29         }
  30         catch(exception& e) {
  31                 m->errorOut(e, "SensSpecCommand", "setParameters");
  32                 exit(1);
  33         }
  34 }
  35 //**********************************************************************************************************************
  36 string SensSpecCommand::getHelpString(){
  37         try {
  38                 string helpString = "";
  39                 helpString += "The sens.spec command....\n";
  40                 return helpString;
  41         }
  42         catch(exception& e) {
  43                 m->errorOut(e, "SensSpecCommand", "getHelpString");
  44                 exit(1);
  45         }
  46 }
  47 //**********************************************************************************************************************
  48 string SensSpecCommand::getOutputFileNameTag(string type, string inputName=""){
  49         try {
  50         string outputFileName = "";
  51                 map<string, vector<string> >::iterator it;
  52
  53         //is this a type this command creates
  54         it = outputTypes.find(type);
  55         if (it == outputTypes.end()) {  m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); }
  56         else {
  57             if (type == "sensspec")            {   outputFileName =  "sensspec";   }
  58             else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true;  }
  59         }
  60         return outputFileName;
  61         }
  62         catch(exception& e) {
  63                 m->errorOut(e, "SensSpecCommand", "getOutputFileNameTag");
  64                 exit(1);
  65         }
  66 }
  67 //**********************************************************************************************************************
  68 SensSpecCommand::SensSpecCommand(){
  69         try {
  70                 abort = true; calledHelp = true;
  71                 setParameters();
  72                 vector<string> tempOutNames;
  73                 outputTypes["sensspec"] = tempOutNames;
  74         }
  75         catch(exception& e) {
  76                 m->errorOut(e, "SensSpecCommand", "SensSpecCommand");
  77                 exit(1);
  78         }
  79 }
  80 //***************************************************************************************************************
  81
  82 SensSpecCommand::SensSpecCommand(string option)  {
  83         try {
  84
  85                 abort = false; calledHelp = false;
  86                 allLines = 1;
  87
  88                 //allow user to run help
  89                 if(option == "help") { help(); abort = true; calledHelp = true; }
  90                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
  91
  92                 else {
  93                         string temp;
  94
  95                         vector<string> myArray = setParameters();
  96
  97                         OptionParser parser(option);
  98                         map<string,string> parameters = parser.getParameters();
  99
 100                         ValidParameters validParameter;
 101                         map<string,string>::iterator it;
 102
 103                         //check to make sure all parameters are valid for command
 104                         for (it = parameters.begin(); it != parameters.end(); it++) {
 105                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
 106                         }
 107
 108                         //initialize outputTypes
 109                         vector<string> tempOutNames;
 110                         outputTypes["sensspec"] = tempOutNames;
 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("list");
 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["list"] = inputDir + it->second;             }
 123                                 }
 124
 125                                 it = parameters.find("phylip");
 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["phylip"] = inputDir + it->second;           }
 131                                 }
 132
 133                                 it = parameters.find("column");
 134                                 //user has given a template 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["column"] = inputDir + it->second;           }
 139                                 }
 140
 141                                 //it = parameters.find("name");
 142                                 //user has given a template 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["name"] = inputDir + it->second;             }
 147                                 //}
 148
 149                         }
 150                         //check for required parameters
 151                         listFile = validParameter.validFile(parameters, "list", true);
 152                         if (listFile == "not found") {
 153                                 listFile = m->getListFile();
 154                                 if (listFile != "") { m->mothurOut("Using " + listFile + " as input file for the list parameter."); m->mothurOutEndLine(); }
 155                                 else {  m->mothurOut("You have no current list file and the list parameter is required."); m->mothurOutEndLine(); abort = true; }
 156                         }
 157                         else if (listFile == "not open") { abort = true; }
 158                         else { m->setListFile(listFile); }
 159
 160                         phylipfile = validParameter.validFile(parameters, "phylip", true);
 161                         if (phylipfile == "not found") { phylipfile = "";  }
 162                         else if (phylipfile == "not open") { abort = true; }
 163                         else { distFile = phylipfile; format = "phylip"; m->setPhylipFile(phylipfile);  }
 164
 165                         columnfile = validParameter.validFile(parameters, "column", true);
 166                         if (columnfile == "not found") { columnfile = ""; }
 167                         else if (columnfile == "not open") { abort = true; }
 168                         else { distFile = columnfile; format = "column";   m->setColumnFile(columnfile); }
 169
 170                         if ((phylipfile == "") && (columnfile == "")) { //is there are current file available for either of these?
 171                                 //give priority to column, then phylip
 172                                 columnfile = m->getColumnFile();
 173                                 if (columnfile != "") {  distFile = columnfile; format = "column";  m->mothurOut("Using " + columnfile + " as input file for the column parameter."); m->mothurOutEndLine(); }
 174                                 else {
 175                                         phylipfile = m->getPhylipFile();
 176                                         if (phylipfile != "") {  distFile = phylipfile; format = "phylip"; m->mothurOut("Using " + phylipfile + " as input file for the phylip parameter."); m->mothurOutEndLine(); }
 177                                         else {
 178                                                 m->mothurOut("No valid current files. You must provide a phylip or column file."); m->mothurOutEndLine();
 179                                                 abort = true;
 180                                         }
 181                                 }
 182                         }else if ((phylipfile != "") && (columnfile != "")) { m->mothurOut("When executing a sens.spec command you must enter ONLY ONE of the following: phylip or column."); m->mothurOutEndLine(); abort = true; }
 183
 184
 185                         //if the user changes the output directory command factory will send this info to us in the output parameter
 186                         outputDir = validParameter.validFile(parameters, "outputdir", false);
 187                         if (outputDir == "not found"){
 188                                 outputDir = "";
 189                                 outputDir += m->hasPath(listFile); //if user entered a file with a path then preserve it
 190                         }
 191
 192                         //check for optional parameter and set defaults
 193                         // ...at some point should added some additional type checking...
 194                         temp = validParameter.validFile(parameters, "hard", false);
 195                         if (temp == "not found"){       hard = 0;       }
 196                         else if(!m->isTrue(temp))       {       hard = 0;       }
 197                         else if(m->isTrue(temp))        {       hard = 1;       }
 198
 199 //                      temp = validParameter.validFile(parameters, "name", true);
 200 //                      if (temp == "not found")        {       nameFile = "";          }
 201 //                      else if(temp == "not open")     {       abort = true;           }
 202 //                      else                                            {       nameFile = temp;        }
 203 //                      cout << "name:\t" << nameFile << endl;
 204
 205                         temp = validParameter.validFile(parameters, "cutoff", false);           if (temp == "not found") { temp = "-1.00"; }
 206                         m->mothurConvert(temp, cutoff);
 207 //                      cout << cutoff << endl;
 208
 209                         temp = validParameter.validFile(parameters, "precision", false);        if (temp == "not found") { temp = "100"; }
 210                         m->mothurConvert(temp, precision);
 211 //                      cout << precision << endl;
 212
 213                         string label = validParameter.validFile(parameters, "label", false);
 214                         if (label == "not found") { label = ""; }
 215                         else {
 216                                 if(label != "all") {  m->splitAtDash(label, labels);  allLines = 0;  }
 217                                 else { allLines = 1;  }
 218                         }
 219
 220                         sensSpecFileName = outputDir + m->getRootName(m->getSimpleName(listFile)) + getOutputFileNameTag("sensspec");
 221                 }
 222         }
 223         catch(exception& e) {
 224                 m->errorOut(e, "SensSpecCommand", "SensSpecCommand");
 225                 exit(1);
 226         }
 227 }
 228 //***************************************************************************************************************
 229
 230 int SensSpecCommand::execute(){
 231         try{
 232                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 233
 234                 setUpOutput();
 235                 outputNames.push_back(sensSpecFileName); outputTypes["sensspec"].push_back(sensSpecFileName);
 236                 if(format == "phylip")          {       processPhylip();        }
 237                 else if(format == "column")     {       processColumn();        }
 238
 239                 if (m->control_pressed) { m->mothurRemove(sensSpecFileName); return 0; }
 240
 241                 m->mothurOutEndLine();
 242                 m->mothurOut("Output File Name: "); m->mothurOutEndLine();
 243                 m->mothurOut(sensSpecFileName); m->mothurOutEndLine();
 244                 m->mothurOutEndLine();
 245
 246
 247                 return 0;
 248         }
 249         catch(exception& e) {
 250                 m->errorOut(e, "SensSpecCommand", "execute");
 251                 exit(1);
 252         }
 253 }
 254
 255 //***************************************************************************************************************
 256
 257 int SensSpecCommand::processPhylip(){
 258         try{
 259                 //probably need some checking to confirm that the names in the distance matrix are the same as those in the list file
 260                 string origCutoff = "";
 261                 bool getCutoff = 0;
 262                 if(cutoff == -1.00)     {       getCutoff = 1;                                                                                                                  }
 263                 else                            {       origCutoff = toString(cutoff);  cutoff += (0.49 / double(precision));   }
 264
 265                 map<string, int> seqMap;
 266                 string seqList;
 267
 268                 InputData input(listFile, "list");
 269                 ListVector* list = input.getListVector();
 270                 string lastLabel = list->getLabel();
 271
 272                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 273                 set<string> processedLabels;
 274                 set<string> userLabels = labels;
 275
 276
 277                 while((list != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 278
 279                         if(m->control_pressed){
 280                 for (int i = 0; i < outputNames.size(); i++){   m->mothurRemove(outputNames[i]);  }  delete list;  return 0;
 281             }
 282
 283                         if(allLines == 1 || labels.count(list->getLabel()) == 1){
 284
 285                                 processedLabels.insert(list->getLabel());
 286                                 userLabels.erase(list->getLabel());
 287
 288                                 //process
 289                                 fillSeqMap(seqMap, list);
 290                                 process(seqMap, list->getLabel(), getCutoff, origCutoff);
 291                         }
 292
 293                         if ((m->anyLabelsToProcess(list->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 294                                 string saveLabel = list->getLabel();
 295
 296                                 delete list;
 297                                 list = input.getListVector(lastLabel);
 298
 299                                 processedLabels.insert(list->getLabel());
 300                                 userLabels.erase(list->getLabel());
 301
 302                                 //process
 303                                 fillSeqMap(seqMap, list);
 304                                 process(seqMap, list->getLabel(), getCutoff, origCutoff);
 305
 306                                 //restore real lastlabel to save below
 307                                 list->setLabel(saveLabel);
 308                         }
 309
 310                         lastLabel = list->getLabel();
 311
 312                         delete list;
 313                         list = input.getListVector();
 314                 }
 315
 316
 317                 //output error messages about any remaining user labels
 318                 set<string>::iterator it;
 319                 bool needToRun = false;
 320                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 321                         m->mothurOut("Your file does not include the label " + *it);
 322                         if (processedLabels.count(lastLabel) != 1) {
 323                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 324                                 needToRun = true;
 325                         }else {
 326                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 327                         }
 328                 }
 329
 330                 //run last label if you need to
 331                 if (needToRun == true)  {
 332                         if (list != NULL) {     delete list;    }
 333                         list = input.getListVector(lastLabel);
 334
 335                         //process
 336                         fillSeqMap(seqMap, list);
 337                         process(seqMap, list->getLabel(), getCutoff, origCutoff);
 338
 339                         delete list;
 340                 }
 341
 342                 return 0;
 343         }
 344         catch(exception& e) {
 345                 m->errorOut(e, "SensSpecCommand", "processPhylip");
 346                 exit(1);
 347         }
 348 }
 349
 350 //***************************************************************************************************************
 351
 352 int SensSpecCommand::fillSeqMap(map<string, int>& seqMap, ListVector*& list){
 353         try {
 354                 //for each otu
 355                 for(int i=0;i<list->getNumBins();i++){
 356
 357                         if (m->control_pressed) { return 0; }
 358
 359                         string seqList = list->get(i);
 360                         int seqListLength = seqList.length();
 361                         string seqName = "";
 362
 363                         //parse bin by name, mapping each name to its otu number
 364                         for(int j=0;j<seqListLength;j++){
 365
 366                                 if(seqList[j] == ','){
 367                                         seqMap[seqName] = i;
 368                                         seqName = "";
 369                                 }
 370                                 else{
 371                                         seqName += seqList[j];
 372                                 }
 373
 374                         }
 375                         seqMap[seqName] = i;
 376                 }
 377
 378                 return 0;
 379         }
 380         catch(exception& e) {
 381                 m->errorOut(e, "SensSpecCommand", "fillSeqMap");
 382                 exit(1);
 383         }
 384 }
 385 //***************************************************************************************************************
 386 int SensSpecCommand::fillSeqPairSet(set<string>& seqPairSet, ListVector*& list){
 387         try {
 388                 int numSeqs = 0;
 389
 390                 //for each otu
 391                 for(int i=0;i<list->getNumBins();i++){
 392
 393                         if (m->control_pressed) { return 0; }
 394
 395                         vector<string> seqNameVector;
 396                         string bin = list->get(i);
 397                         m->splitAtComma(bin, seqNameVector);
 398
 399                         numSeqs += seqNameVector.size();
 400
 401                         for(int j=0;j<seqNameVector.size();j++){
 402                                 string seqPairString = "";
 403                                 for(int k=0;k<j;k++){
 404                                         if(seqNameVector[j] < seqNameVector[k]) {       seqPairString = seqNameVector[j] + '\t' + seqNameVector[k];     }
 405                                         else                                                                    {       seqPairString = seqNameVector[k] + '\t' + seqNameVector[j];     }
 406                                         seqPairSet.insert(seqPairString);
 407                                 }
 408                         }
 409                 }
 410
 411                 return numSeqs;
 412         }
 413         catch(exception& e) {
 414                 m->errorOut(e, "SensSpecCommand", "fillSeqPairSet");
 415                 exit(1);
 416         }
 417 }
 418 //***************************************************************************************************************
 419 int SensSpecCommand::process(map<string, int>& seqMap, string label, bool& getCutoff, string& origCutoff){
 420         try {
 421
 422                 int lNumSeqs = seqMap.size();
 423                 int pNumSeqs = 0;
 424
 425                 ifstream phylipFile;
 426                 m->openInputFile(distFile, phylipFile);
 427                 phylipFile >> pNumSeqs;
 428                 if(pNumSeqs != lNumSeqs){       m->mothurOut("numSeq mismatch!\n"); /*m->control_pressed = true;*/ }
 429
 430                 string seqName;
 431                 double distance;
 432                 vector<int> otuIndices(lNumSeqs, -1);
 433
 434                 truePositives = 0;
 435                 falsePositives = 0;
 436                 trueNegatives = 0;
 437                 falseNegatives = 0;
 438
 439                 if(getCutoff == 1){
 440                         if(label != "unique"){
 441                                 origCutoff = label;
 442                                 convert(label, cutoff);
 443                                 if(hard == 0){  cutoff += (0.49 / double(precision));   }
 444                         }
 445                         else{
 446                                 origCutoff = "unique";
 447                                 cutoff = 0.0000;
 448                         }
 449                 }
 450
 451                 m->mothurOut(label); m->mothurOutEndLine();
 452
 453                 for(int i=0;i<pNumSeqs;i++){
 454
 455                         if (m->control_pressed) { return 0; }
 456
 457                         phylipFile >> seqName;
 458                         otuIndices[i] = seqMap[seqName];
 459
 460                         for(int j=0;j<i;j++){
 461                                 phylipFile >> distance;
 462
 463                                 if(distance <= cutoff){
 464                                         if(otuIndices[i] == otuIndices[j])      {       truePositives++;        }
 465                                         else                                                            {       falseNegatives++;       }
 466                                 }
 467                                 else{
 468                                         if(otuIndices[i] == otuIndices[j])      {       falsePositives++;       }
 469                                         else                                                            {       trueNegatives++;        }
 470                                 }
 471                         }
 472                 }
 473                 phylipFile.close();
 474
 475                 outputStatistics(label, origCutoff);
 476
 477                 return 0;
 478         }
 479         catch(exception& e) {
 480                 m->errorOut(e, "SensSpecCommand", "process");
 481                 exit(1);
 482         }
 483 }
 484 //***************************************************************************************************************
 485 int SensSpecCommand::process(set<string>& seqPairSet, string label, bool& getCutoff, string& origCutoff, int numSeqs){
 486         try {
 487                 int numDists = (numSeqs * (numSeqs-1) / 2);
 488
 489                 ifstream columnFile;
 490                 m->openInputFile(distFile, columnFile);
 491                 string seqNameA, seqNameB, seqPairString;
 492                 double distance;
 493
 494                 truePositives = 0;
 495                 falsePositives = 0;
 496                 trueNegatives = numDists;
 497                 falseNegatives = 0;
 498
 499                 if(getCutoff == 1){
 500                         if(label != "unique"){
 501                                 origCutoff = label;
 502                                 convert(label, cutoff);
 503                                 if(hard == 0){  cutoff += (0.49 / double(precision));   }
 504                         }
 505                         else{
 506                                 origCutoff = "unique";
 507                                 cutoff = 0.0000;
 508                         }
 509                 }
 510
 511                 m->mothurOut(label); m->mothurOutEndLine();
 512
 513                 while(columnFile){
 514                         columnFile >> seqNameA >> seqNameB >> distance;
 515                         if(seqNameA < seqNameB) {       seqPairString = seqNameA + '\t' + seqNameB;     }
 516                         else                                    {       seqPairString = seqNameB + '\t' + seqNameA;     }
 517
 518                         set<string>::iterator it = seqPairSet.find(seqPairString);
 519
 520                         if(distance <= cutoff){
 521                                 if(it != seqPairSet.end()){
 522                                         truePositives++;
 523                                         seqPairSet.erase(it);
 524                                 }
 525                                 else{
 526                                         falseNegatives++;
 527                                 }
 528                                 trueNegatives--;
 529                         }
 530                         else if(it != seqPairSet.end()){
 531                                 falsePositives++;
 532                                 trueNegatives--;
 533                                 seqPairSet.erase(it);
 534                         }
 535
 536                         m->gobble(columnFile);
 537                 }
 538                 falsePositives += seqPairSet.size();
 539
 540                 outputStatistics(label, origCutoff);
 541
 542
 543                 return 0;
 544         }
 545         catch(exception& e) {
 546                 m->errorOut(e, "SensSpecCommand", "process");
 547                 exit(1);
 548         }
 549 }
 550 //***************************************************************************************************************
 551
 552 int SensSpecCommand::processColumn(){
 553         try{
 554                 string origCutoff = "";
 555                 bool getCutoff = 0;
 556                 if(cutoff == -1.00)     {       getCutoff = 1;                                                                                                                  }
 557                 else                            {       origCutoff = toString(cutoff);  cutoff += (0.49 / double(precision));   }
 558
 559                 set<string> seqPairSet;
 560                 int numSeqs = 0;
 561
 562                 InputData input(listFile, "list");
 563                 ListVector* list = input.getListVector();
 564                 string lastLabel = list->getLabel();
 565
 566                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 567                 set<string> processedLabels;
 568                 set<string> userLabels = labels;
 569
 570
 571                 while((list != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 572
 573                         if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        m->mothurRemove(outputNames[i]);  }  delete list;  return 0;  }
 574
 575                         if(allLines == 1 || labels.count(list->getLabel()) == 1){
 576
 577                                 processedLabels.insert(list->getLabel());
 578                                 userLabels.erase(list->getLabel());
 579
 580                                 //process
 581                                 numSeqs = fillSeqPairSet(seqPairSet, list);
 582                                 process(seqPairSet, list->getLabel(), getCutoff, origCutoff, numSeqs);
 583                         }
 584
 585                         if ((m->anyLabelsToProcess(list->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 586                                 string saveLabel = list->getLabel();
 587
 588                                 delete list;
 589                                 list = input.getListVector(lastLabel);
 590
 591                                 processedLabels.insert(list->getLabel());
 592                                 userLabels.erase(list->getLabel());
 593
 594                                 //process
 595                                 numSeqs = fillSeqPairSet(seqPairSet, list);
 596                                 process(seqPairSet, list->getLabel(), getCutoff, origCutoff, numSeqs);
 597
 598                                 //restore real lastlabel to save below
 599                                 list->setLabel(saveLabel);
 600                         }
 601
 602                         lastLabel = list->getLabel();
 603
 604                         delete list;
 605                         list = input.getListVector();
 606                 }
 607
 608
 609                 //output error messages about any remaining user labels
 610                 set<string>::iterator it;
 611                 bool needToRun = false;
 612                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 613                         m->mothurOut("Your file does not include the label " + *it);
 614                         if (processedLabels.count(lastLabel) != 1) {
 615                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 616                                 needToRun = true;
 617                         }else {
 618                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 619                         }
 620                 }
 621
 622                 //run last label if you need to
 623                 if (needToRun == true)  {
 624                         if (list != NULL) {     delete list;    }
 625                         list = input.getListVector(lastLabel);
 626
 627                         //process
 628                         numSeqs = fillSeqPairSet(seqPairSet, list);
 629                         delete list;
 630                         process(seqPairSet, list->getLabel(), getCutoff, origCutoff, numSeqs);
 631                 }
 632
 633                 return 0;
 634         }
 635         catch(exception& e) {
 636                 m->errorOut(e, "SensSpecCommand", "processColumn");
 637                 exit(1);
 638         }
 639 }
 640
 641 //***************************************************************************************************************
 642
 643 void SensSpecCommand::setUpOutput(){
 644         try{
 645                 ofstream sensSpecFile;
 646                 m->openOutputFile(sensSpecFileName, sensSpecFile);
 647
 648                 sensSpecFile << "label\tcutoff\ttp\ttn\tfp\tfn\tsensitivity\tspecificity\tppv\tnpv\tfdr\taccuracy\tmcc\tf1score\n";
 649
 650                 sensSpecFile.close();
 651         }
 652         catch(exception& e) {
 653                 m->errorOut(e, "SensSpecCommand", "setUpOutput");
 654                 exit(1);
 655         }
 656 }
 657
 658 //***************************************************************************************************************
 659
 660 void SensSpecCommand::outputStatistics(string label, string cutoff){
 661         try{
 662                 double tp = (double) truePositives;
 663                 double fp = (double) falsePositives;
 664                 double tn = (double) trueNegatives;
 665                 double fn = (double) falseNegatives;
 666
 667                 double p = tp + fn;
 668                 double n = fp + tn;
 669                 double pPrime = tp + fp;
 670                 double nPrime = tn + fn;
 671
 672                 double sensitivity = tp / p;
 673                 double specificity = tn / n;
 674                 double positivePredictiveValue = tp / pPrime;
 675                 double negativePredictiveValue = tn / nPrime;
 676                 double falseDiscoveryRate = fp / pPrime;
 677
 678                 double accuracy = (tp + tn) / (p + n);
 679                 double matthewsCorrCoef = (tp * tn - fp * fn) / sqrt(p * n * pPrime * nPrime);  if(p == 0 || n == 0){   matthewsCorrCoef = 0;   }
 680                 double f1Score = 2.0 * tp / (p + pPrime);
 681
 682
 683                 if(p == 0)                      {       sensitivity = 0;        matthewsCorrCoef = 0;   }
 684                 if(n == 0)                      {       specificity = 0;        matthewsCorrCoef = 0;   }
 685                 if(p + n == 0)          {       accuracy = 0;                                                           }
 686                 if(p + pPrime == 0)     {       f1Score = 0;                                                            }
 687                 if(pPrime == 0)         {       positivePredictiveValue = 0;    falseDiscoveryRate = 0; matthewsCorrCoef = 0;   }
 688                 if(nPrime == 0)         {       negativePredictiveValue = 0;    matthewsCorrCoef = 0;                                                   }
 689
 690                 ofstream sensSpecFile;
 691                 m->openOutputFileAppend(sensSpecFileName, sensSpecFile);
 692
 693                 sensSpecFile << label << '\t' << cutoff << '\t';
 694                 sensSpecFile << truePositives << '\t' << trueNegatives << '\t' << falsePositives << '\t' << falseNegatives << '\t';
 695                 sensSpecFile << setprecision(4);
 696                 sensSpecFile << sensitivity << '\t' << specificity << '\t' << positivePredictiveValue << '\t' << negativePredictiveValue << '\t';
 697                 sensSpecFile << falseDiscoveryRate << '\t' << accuracy << '\t' << matthewsCorrCoef << '\t' << f1Score << endl;
 698
 699                 sensSpecFile.close();
 700         }
 701         catch(exception& e) {
 702                 m->errorOut(e, "SensSpecCommand", "outputStatistics");
 703                 exit(1);
 704         }
 705 }
 706
 707 //***************************************************************************************************************
 708
 709
 710