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