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
  14 SensSpecCommand::SensSpecCommand(string option)  {
  15         try {
  16
  17                 abort = false;
  18
  19                 //allow user to run help
  20                 if(option == "help") { help(); abort = true; }
  21
  22                 else {
  23                         string temp;
  24
  25                         //valid paramters for this command
  26                         string AlignArray[] =  {"list", "phylip", "column", "name", "hard", "label", "cutoff", "outputdir", "inputdir"};
  27
  28                         vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
  29
  30                         OptionParser parser(option);
  31                         map<string,string> parameters = parser.getParameters();
  32
  33                         ValidParameters validParameter;
  34                         map<string,string>::iterator it;
  35
  36                         //check to make sure all parameters are valid for command
  37                         for (it = parameters.begin(); it != parameters.end(); it++) {
  38                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  39                         }
  40
  41                         //if the user changes the input directory command factory will send this info to us in the output parameter
  42                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
  43                         if (inputDir == "not found"){   inputDir = "";          }
  44                         else {
  45                                 string path;
  46                                 it = parameters.find("list");
  47                                 //user has given a template file
  48                                 if(it != parameters.end()){
  49                                         path = hasPath(it->second);
  50                                         //if the user has not given a path then, add inputdir. else leave path alone.
  51                                         if (path == "") {       parameters["list"] = inputDir + it->second;             }
  52                                 }
  53
  54                                 it = parameters.find("phylip");
  55                                 //user has given a template file
  56                                 if(it != parameters.end()){
  57                                         path = hasPath(it->second);
  58                                         //if the user has not given a path then, add inputdir. else leave path alone.
  59                                         if (path == "") {       parameters["phylip"] = inputDir + it->second;           }
  60                                 }
  61
  62                                 it = parameters.find("column");
  63                                 //user has given a template file
  64                                 if(it != parameters.end()){
  65                                         path = hasPath(it->second);
  66                                         //if the user has not given a path then, add inputdir. else leave path alone.
  67                                         if (path == "") {       parameters["column"] = inputDir + it->second;           }
  68                                 }
  69
  70                                 it = parameters.find("name");
  71                                 //user has given a template file
  72                                 if(it != parameters.end()){
  73                                         path = hasPath(it->second);
  74                                         //if the user has not given a path then, add inputdir. else leave path alone.
  75                                         if (path == "") {       parameters["name"] = inputDir + it->second;             }
  76                                 }
  77
  78                         }
  79                         //check for required parameters
  80                         listFile = validParameter.validFile(parameters, "list", true);
  81                         if (listFile == "not found") { m->mothurOut("list is a required parameter for the sens.spec command."); m->mothurOutEndLine(); abort = true; }
  82                         else if (listFile == "not open") { abort = true; }
  83
  84                         distFile = validParameter.validFile(parameters, "column", true);
  85                         format = "column";
  86                         if(distFile == "not found")     {
  87                                 distFile = validParameter.validFile(parameters, "phylip", true);
  88                                 format = "phylip";
  89                         }
  90                         if(distFile == "not found") { m->mothurOut("either column or phylip are required for the sens.spec command."); m->mothurOutEndLine(); abort = true; }
  91                         else if (distFile == "not open") { abort = true; }
  92
  93                         //if the user changes the output directory command factory will send this info to us in the output parameter
  94                         outputDir = validParameter.validFile(parameters, "outputdir", false);
  95                         if (outputDir == "not found"){
  96                                 outputDir = "";
  97                                 outputDir += hasPath(listFile); //if user entered a file with a path then preserve it
  98                         }
  99
 100                         //check for optional parameter and set defaults
 101                         // ...at some point should added some additional type checking...
 102                         temp = validParameter.validFile(parameters, "hard", false);
 103                         if (temp == "not found"){       hard = 0;       }
 104                         else if(!isTrue(temp))  {       hard = 0;       }
 105                         else if(isTrue(temp))   {       hard = 1;       }
 106
 107 //                      temp = validParameter.validFile(parameters, "name", true);
 108 //                      if (temp == "not found")        {       nameFile = "";          }
 109 //                      else if(temp == "not open")     {       abort = true;           }
 110 //                      else                                            {       nameFile = temp;        }
 111 //                      cout << "name:\t" << nameFile << endl;
 112
 113                         temp = validParameter.validFile(parameters, "cutoff", false);           if (temp == "not found") { temp = ""; }
 114                         convert(temp, cutoff);
 115                         cout << "cutoff:\t" << cutoff << endl;
 116
 117 //                      cutoff = 0.0349;
 118
 119                         lineLabel = validParameter.validFile(parameters, "label", false);       if (lineLabel == "not found") { lineLabel = ""; }
 120
 121                 }
 122
 123         }
 124         catch(exception& e) {
 125                 m->errorOut(e, "SensSpecCommand", "SensSpecCommand");
 126                 exit(1);
 127         }
 128 }
 129 //**********************************************************************************************************************
 130
 131 void SensSpecCommand::help(){
 132         try {
 133                 m->mothurOut("The sens.spec command reads a fastaFile and creates .....\n");
 134
 135
 136
 137                 m->mothurOut("Example sens.spec(...).\n");
 138                 m->mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n");
 139                 m->mothurOut("For more details please check out the wiki http://www.mothur.org/wiki/Trim.seqs .\n\n");
 140
 141         }
 142         catch(exception& e) {
 143                 m->errorOut(e, "SensSpecCommand", "help");
 144                 exit(1);
 145         }
 146 }
 147
 148
 149 //***************************************************************************************************************
 150
 151 SensSpecCommand::~SensSpecCommand(){    /*      do nothing      */      }
 152
 153 //***************************************************************************************************************
 154
 155 int SensSpecCommand::execute(){
 156         try{
 157                 if (abort == true) { return 0; }
 158
 159                 if(format == "phylip")          {       processPhylip();        }
 160 //              else if(format == "column")     {       processColumn(seqMap);  }
 161
 162
 163 //              string seqList;
 164 //              map<string, int> seqMap;
 165
 166
 167
 168                 return 0;
 169         }
 170         catch(exception& e) {
 171                 m->errorOut(e, "SensSpecCommand", "execute");
 172                 exit(1);
 173         }
 174 }
 175
 176 //***************************************************************************************************************
 177
 178 void SensSpecCommand::processPhylip(){
 179         try{
 180                 ifstream inputListFile;
 181                 openInputFile(listFile, inputListFile);
 182
 183                 string label;
 184                 int numOTUs;
 185
 186                 map<string, int> seqMap;
 187                 string seqList;
 188
 189                 while(inputListFile){
 190                         inputListFile >> label >> numOTUs;
 191                         for(int i=0;i<numOTUs;i++){
 192                                 inputListFile >> seqList;
 193                                 int seqListLength = seqList.length();
 194                                 string seqName = "";
 195                                 for(int j=0;j<seqListLength;j++){
 196
 197                                         if(seqList[j] == ','){
 198                                                 seqMap[seqName] = i;
 199                                                 seqName = "";
 200                                         }
 201                                         else{
 202                                                 seqName += seqList[j];
 203                                         }
 204
 205                                 }
 206                                 seqMap[seqName] = i;
 207                         }
 208                 }
 209                 inputListFile.close();
 210
 211                 int lNumSeqs = seqMap.size();
 212                 int pNumSeqs = 0;
 213
 214                 ifstream phylipFile;
 215                 openInputFile(distFile, phylipFile);
 216                 phylipFile >> pNumSeqs;
 217                 if(pNumSeqs != lNumSeqs){       cout << "numSeq mismatch!" << endl;     }
 218
 219                 string seqName;
 220                 double distance;
 221                 vector<int> otuIndices(lNumSeqs, -1);
 222
 223                 truePositives = 0;
 224                 falsePositives = 0;
 225                 trueNegatives = 0;
 226                 falseNegatives = 0;
 227
 228
 229                 for(int i=0;i<lNumSeqs;i++){
 230                         phylipFile >> seqName;
 231                         otuIndices[i] = seqMap[seqName];
 232
 233                         for(int j=0;j<i;j++){
 234                                 phylipFile >> distance;
 235                                 if(distance <= cutoff){
 236                                         if(otuIndices[i] == otuIndices[j]){
 237                                                 truePositives++;
 238                                         }
 239                                         else{
 240                                                 falseNegatives++;
 241                                         }
 242                                 }
 243                                 else{
 244                                         if(otuIndices[i] == otuIndices[j]){
 245                                                 falsePositives++;
 246                                         }
 247                                         else{
 248                                                 trueNegatives++;
 249                                         }
 250                                 }
 251                         }
 252                 }
 253                 phylipFile.close();
 254
 255                 cout << "truePositives:\t"      << truePositives        << endl;
 256                 cout << "trueNegatives:\t"      << trueNegatives        << endl;
 257                 cout << "falsePositives:\t"     << falsePositives       << endl;
 258                 cout << "falseNegatives:\t"     << falseNegatives       << endl;
 259         }
 260         catch(exception& e) {
 261                 m->errorOut(e, "SensSpecCommand", "processPhylip");
 262                 exit(1);
 263         }
 264 }
 265
 266 //***************************************************************************************************************
 267
 268 //void SensSpecCommand::processColumn(map<string, int> seqMap){
 269 //
 270 //      truePositives = 0;
 271 //      falsePositives = 0;
 272 //      trueNegatives = numDists;
 273 //      falseNegatives = 0;
 274 //
 275 //      ifstream columnFile;
 276 //      openInputFile(distFile, columnFile);
 277 //
 278 //      string seqNameA, seqNameB, oldSeqNameA;
 279 //      int otuA, otuB, oldOTUA;
 280 //      double distance;
 281 //
 282 //      while(columnFile){
 283 //              columnFile >> seqNameA >> seqNameB >> distance;
 284 //
 285 //              if(seqNameA == oldSeqNameA)     {       otuA = oldOTUA; }
 286 //              else                                            {       otuA = seqMap[seqNameA];        oldOTUA = otuA; }
 287 //
 288 //              otuB = seqMap[seqNameB];
 289 //
 290 //              if(distance <= cutoff){
 291 //                      if(otuA == otuB){
 292 //                              truePositives++;
 293 //                      }
 294 //                      else{
 295 //                              falseNegatives++;
 296 //                      }
 297 //                      trueNegatives--;
 298 //              }
 299 //              else{
 300 //                      if(otuA == otuB){
 301 //                              falsePositives++;
 302 //                              trueNegatives--;
 303 //                      }
 304 //              }
 305 //
 306 //              gobble(columnFile);
 307 //      }
 308 //      columnFile.close();
 309 //
 310 //      cout << "truePositives:\t"      << truePositives        << endl;
 311 //      cout << "trueNegatives:\t"      << trueNegatives        << endl;
 312 //      cout << "falsePositives:\t"     << falsePositives       << endl;
 313 //      cout << "falseNegatives:\t"     << falseNegatives       << endl;
 314 //}
 315
 316 //***************************************************************************************************************