filterseqscommand.cpp

   1 /*
   2  *  filterseqscommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by Thomas Ryabin on 5/4/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "filterseqscommand.h"
  11
  12 /**************************************************************************************/
  13 void FilterSeqsCommand::doTrump() {
  14         trump = globaldata->getTrump();
  15         for(int i = 0; i < db->size(); i++) {
  16                 Sequence cur = db->get(i);
  17                 string curAligned = cur.getAligned();
  18                 for(int j = 0; j < curAligned.length(); j++) {
  19                         string curChar = curAligned.substr(j, 1);
  20                         if(curChar.compare(trump) == 0)
  21                                 columnsToRemove[j] = true;
  22                 }
  23         }
  24 }
  25
  26 /**************************************************************************************/
  27 void FilterSeqsCommand::doSoft() {
  28         soft = atoi(globaldata->getSoft().c_str());
  29         vector<vector<int> > columnSymbolSums;
  30         vector<vector<string> > columnSymbols;
  31         for(int i = 0; i < db->get(0).getLength(); i++) {
  32                 vector<string> symbols;
  33                 vector<int> sums;
  34                 columnSymbols.push_back(symbols);
  35                 columnSymbolSums.push_back(sums);
  36         }
  37
  38         for(int i = 0; i < db->size(); i++) {
  39                 Sequence cur = db->get(i);
  40                 string curAligned = cur.getAligned();
  41
  42                 for(int j = 0; j < curAligned.length(); j++) {
  43                         string curChar = curAligned.substr(j, 1);
  44                         vector<string> curColumnSymbols = columnSymbols[j];
  45                         bool newSymbol = true;
  46
  47                         for(int k = 0; k < curColumnSymbols.size(); k++)
  48                                 if(curChar.compare(curColumnSymbols[k]) == 0) {
  49                                         newSymbol = false;
  50                                         columnSymbolSums[j][k]++;
  51                                 }
  52
  53                         if(newSymbol) {
  54                                 columnSymbols[j].push_back(curChar);
  55                                 columnSymbolSums[j].push_back(1);
  56                         }
  57                 }
  58         }
  59
  60
  61         for(int i = 0; i < columnSymbolSums.size(); i++) {
  62                 int totalSum = 0;
  63                 int max = 0;
  64                 vector<int> curColumnSymbols = columnSymbolSums[i];
  65
  66                 for(int j = 0; j < curColumnSymbols.size(); j++) {
  67                         int curSum = curColumnSymbols[j];
  68                         //cout << columnSymbols[i][j] << ": " << curSum << "\n";
  69                         if(curSum > max)
  70                                 max = curSum;
  71                         totalSum += curSum;
  72                 }
  73                 //cout << "\n";
  74
  75                 if((double)max/(double)totalSum * 100 < soft)
  76                         columnsToRemove[i] = true;
  77         }
  78 }
  79
  80 /**************************************************************************************/
  81 void FilterSeqsCommand::doFilter() {
  82         filter = globaldata->getFilter();
  83         ifstream filehandle;
  84         openInputFile(filter, filehandle);
  85
  86         char c;
  87         int count = 0;
  88         while(!filehandle.eof()) {
  89                 c = filehandle.get();
  90                 if(c == '0')
  91                         columnsToRemove[count] = true;
  92                 count++;
  93         }
  94 }
  95
  96 /**************************************************************************************/
  97 int FilterSeqsCommand::execute() {
  98         try {
  99                 globaldata = GlobalData::getInstance();
 100                 filename = globaldata->inputFileName;
 101
 102                 if(globaldata->getFastaFile() != "") {
 103                         readSeqs =  new ReadFasta(filename); }
 104                 else if(globaldata->getNexusFile() != "") {
 105                         readSeqs = new ReadNexus(filename); }
 106                 else if(globaldata->getClustalFile() != "") {
 107                         readSeqs = new ReadClustal(filename); }
 108                 else if(globaldata->getPhylipFile() != "") {
 109                         readSeqs = new ReadPhylip(filename); }
 110
 111                 readSeqs->read();
 112                 db = readSeqs->getDB();
 113
 114                 //for(int i = 0; i < db->size(); i++) {
 115 //                      cout << db->get(i).getLength() << "\n" << db->get(i).getName() << ": " << db->get(i).getAligned() << "\n\n";
 116 //              }
 117
 118                 for(int i = 0; i < db->get(0).getLength(); i++)
 119                         columnsToRemove.push_back(false);
 120
 121
 122                 if(globaldata->getTrump().compare("") != 0)
 123                         doTrump();
 124                 else if(globaldata->getSoft().compare("") != 0)
 125                         doSoft();
 126                 else if(globaldata->getFilter().compare("") != 0)
 127                         doFilter();
 128
 129                 //for(int i = 0; i < columnsToRemove.size(); i++)
 130 //              {
 131 //                      cout << "Remove Column " << i << " = ";
 132 //                      if(columnsToRemove[i])
 133 //                              cout << "true\n";
 134 //                      else
 135 //                              cout << "false\n";
 136 //              }
 137
 138
 139                 //Creating the new SequenceDB
 140                 SequenceDB newDB;
 141                 for(int i = 0; i < db->size(); i++) {
 142                         Sequence curSeq = db->get(i);
 143                         string curAligned = curSeq.getAligned();
 144                         string curName = curSeq.getName();
 145                         string newAligned = "";
 146                         for(int j = 0; j < curAligned.length(); j++)
 147                                 if(!columnsToRemove[j])
 148                                         newAligned += curAligned.substr(j, 1);
 149
 150                         Sequence newSeq(curName, newAligned);
 151                         newDB.add(newSeq);
 152                 }
 153
 154                 string newFileName = getRootName(filename) + "filter.fa";
 155                 ofstream outfile;
 156                 outfile.open(newFileName.c_str());
 157                 newDB.print(outfile);
 158                 outfile.close();
 159
 160                 globaldata->clear();
 161                 //delete db;
 162                 //delete newDB;
 163
 164                 return 0;
 165         }
 166         catch(exception& e) {
 167                 cout << "Standard Error: " << e.what() << " has occurred in the FilterSeqsCommand class Function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 168                 exit(1);
 169         }
 170         catch(...) {
 171                 cout << "An unknown error has occurred in the FilterSeqsCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 172                 exit(1);
 173         }
 174 }