getoturepcommand.cpp

   1 /*
   2  *  getoturepcommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 4/6/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "getoturepcommand.h"
  11
  12 //**********************************************************************************************************************
  13 GetOTURepCommand::GetOTURepCommand(){
  14         try{
  15                 globaldata = GlobalData::getInstance();
  16
  17                 if(globaldata->gSparseMatrix != NULL)   {       matrix = new SparseMatrix(*globaldata->gSparseMatrix);          }
  18
  19                 //listOfNames bin 0 = first name read in distance matrix, listOfNames bin 1 = second name read in distance matrix
  20                 if(globaldata->gListVector != NULL)             {
  21                         listOfNames = new ListVector(*globaldata->gListVector);
  22
  23                         //map names to rows in sparsematrix
  24                         for (int i = 0; i < listOfNames->size(); i++) {
  25                                 nameToIndex[listOfNames->get(i)] = i;
  26                         }
  27                 }else { cout << "error" << endl; }
  28
  29
  30                 fastafile = globaldata->getFastaFile();
  31                 namesfile = globaldata->getNameFile();
  32                 openInputFile(fastafile, in);
  33
  34                 fasta = new FastaMap();
  35
  36         }
  37         catch(exception& e) {
  38                 cout << "Standard Error: " << e.what() << " has occurred in the GetOTURepCommand class Function GetOTURepCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  39                 exit(1);
  40         }
  41         catch(...) {
  42                 cout << "An unknown error has occurred in the GetOTURepCommand class function GetOTURepCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  43                 exit(1);
  44         }
  45 }
  46 //**********************************************************************************************************************
  47
  48 GetOTURepCommand::~GetOTURepCommand(){
  49         delete matrix;
  50         delete list;
  51         delete input;
  52         delete read;
  53         delete fasta;
  54 }
  55
  56 //**********************************************************************************************************************
  57
  58 int GetOTURepCommand::execute(){
  59         try {
  60                 int count = 1;
  61                 string nameRep, name, sequence;
  62
  63                 //read fastafile
  64                 fasta->readFastaFile(in);
  65
  66                 //set format to list so input can get listvector
  67                 globaldata->setFormat("list");
  68
  69                 //if user gave a namesfile then use it
  70                 if (namesfile != "") {
  71                         readNamesFile();
  72                 }
  73
  74                 //read list file
  75                 read = new ReadPhilFile(globaldata->getListFile());
  76                 read->read(&*globaldata);
  77
  78                 input = globaldata->ginput;
  79                 list = globaldata->gListVector;
  80
  81                 while(list != NULL){
  82
  83                         if(globaldata->allLines == 1 || globaldata->lines.count(count) == 1 || globaldata->labels.count(list->getLabel()) == 1){
  84
  85                                 //create output file
  86                                 string outputFileName = getRootName(globaldata->getListFile()) + list->getLabel() + ".fastarep";
  87                                 openOutputFile(outputFileName, out);
  88
  89                                 cout << list->getLabel() << '\t' << count << endl;
  90
  91                                 //for each bin in the list vector
  92                                 for (int i = 0; i < list->size(); i++) {
  93                                         nameRep = FindRep(i);
  94
  95                                         //print out name and sequence for that bin
  96                                         sequence = fasta->getSequence(nameRep);
  97
  98                                         if (sequence != "not found") {
  99                                                 nameRep = nameRep + "bin" + toString(i+1);
 100                                                 out << ">" << nameRep << endl;
 101                                                 out << sequence << endl;
 102                                         }else {
 103                                                 cout << nameRep << " is missing from your fasta or name file. Please correct. " << endl;
 104                                                 remove(outputFileName.c_str());
 105                                                 return 0;
 106                                         }
 107                                 }
 108                         }
 109
 110                         list = input->getListVector();
 111                         count++;
 112                 }
 113
 114
 115                 return 0;
 116         }
 117         catch(exception& e) {
 118                 cout << "Standard Error: " << e.what() << " has occurred in the GetOTURepCommand class Function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 119                 exit(1);
 120         }
 121         catch(...) {
 122                 cout << "An unknown error has occurred in the GetOTURepCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 123                 exit(1);
 124         }
 125
 126 }
 127
 128 //**********************************************************************************************************************
 129 void GetOTURepCommand::readNamesFile() {
 130         try {
 131                 vector<string> dupNames;
 132                 openInputFile(namesfile, inNames);
 133
 134                 string name, names, sequence;
 135
 136                 while(inNames){
 137                         inNames >> name;                        //read from first column  A
 138                         inNames >> names;               //read from second column  A,B,C,D
 139
 140                         dupNames.clear();
 141
 142                         //parse names into vector
 143                         splitAtComma(names, dupNames);
 144
 145                         //store names in fasta map
 146                         sequence = fasta->getSequence(name);
 147                         for (int i = 0; i < dupNames.size(); i++) {
 148                                 fasta->push_back(dupNames[i], sequence);
 149                         }
 150
 151                         gobble(inNames);
 152                 }
 153                 inNames.close();
 154
 155         }
 156         catch(exception& e) {
 157                 cout << "Standard Error: " << e.what() << " has occurred in the GetOTURepCommand class Function readNamesFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 158                 exit(1);
 159         }
 160         catch(...) {
 161                 cout << "An unknown error has occurred in the GetOTURepCommand class function readNamesFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 162                 exit(1);
 163         }
 164 }
 165 //**********************************************************************************************************************
 166 string GetOTURepCommand::FindRep(int bin) {
 167         try{
 168                 vector<string> names;
 169                 map<string, float> sums;
 170                 map<string, float>::iterator it4;
 171                 map<int, string> binMap; //subset of namesToIndex - just member of this bin
 172                 string binnames;
 173                 float min = 10000;
 174                 string minName;
 175
 176                 binnames = list->get(bin);
 177
 178                 //parse names into vector
 179                 splitAtComma(binnames, names);
 180
 181                 //if only 1 sequence in bin then that's the rep
 182                 if (names.size() == 1) { return names[0]; }
 183                 else {
 184                         //fill binMap
 185                         for (int i = 0; i < names.size(); i++) {
 186                                 for (it3 = nameToIndex.begin(); it3 != nameToIndex.end(); it3++) {
 187                                         if (it3->first == names[i]) {
 188                                                 binMap[it3->second] = it3->first;
 189
 190                                                 //initialize sums map
 191                                                 sums[it3->first] = 0.0;
 192                                                 break;
 193                                         }
 194                                 }
 195                         }
 196
 197                         //go through each cell in the sparsematrix
 198                         for(MatData currentCell = matrix->begin(); currentCell != matrix->end(); currentCell++){
 199                                 //is this a distance between 2 members of this bin
 200                                 it = binMap.find(currentCell->row);
 201                                 it2 = binMap.find(currentCell->column);
 202
 203                                 //sum the distance of the sequences in the bin to eachother
 204                                 if ((it != binMap.end()) && (it2 != binMap.end())) {
 205                                         //this is a cell that repesents the distance between to of this bins members
 206                                         sums[it->second] += currentCell->dist;
 207                                         sums[it2->second] += currentCell->dist;
 208                                 }
 209                         }
 210
 211                         //smallest sum is the representative
 212                         for (it4 = sums.begin(); it4 != sums.end(); it4++) {
 213                                 if (it4->second < min) {
 214                                         min = it4->second;
 215                                         minName = it4->first;
 216                                 }
 217
 218                         }
 219
 220                         return minName;
 221                 }
 222
 223         }
 224         catch(exception& e) {
 225                 cout << "Standard Error: " << e.what() << " has occurred in the GetOTURepCommand class Function FindRep. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 226                 exit(1);
 227         }
 228         catch(...) {
 229                 cout << "An unknown error has occurred in the GetOTURepCommand class function FindRep. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 230                 exit(1);
 231         }
 232 }
 233
 234
 235
 236
 237