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