chimeraseqscommand.cpp

   1 /*
   2  *  chimeraseqscommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 6/29/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "chimeraseqscommand.h"
  11
  12 //***************************************************************************************************************
  13
  14 ChimeraSeqsCommand::ChimeraSeqsCommand(string option){
  15         try {
  16                 abort = false;
  17
  18                 //allow user to run help
  19                 if(option == "help") { help(); abort = true; }
  20
  21                 else {
  22                         //valid paramters for this command
  23                         string Array[] =  {"fasta", "filter", "correction", "processors", "method" };
  24                         vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  25
  26                         OptionParser parser(option);
  27                         map<string,string> parameters = parser.getParameters();
  28
  29                         ValidParameters validParameter;
  30
  31                         //check to make sure all parameters are valid for command
  32                         for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
  33                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  34                         }
  35
  36                         //check for required parameters
  37                         fastafile = validParameter.validFile(parameters, "fasta", true);
  38                         if (fastafile == "not open") { abort = true; }
  39                         else if (fastafile == "not found") { fastafile = ""; mothurOut("fasta is a required parameter for the chimera.seqs command."); mothurOutEndLine(); abort = true;  }
  40
  41                         string temp;
  42                         temp = validParameter.validFile(parameters, "filter", false);                   if (temp == "not found") { temp = "F"; }
  43                         filter = isTrue(temp);
  44
  45                         temp = validParameter.validFile(parameters, "correction", false);               if (temp == "not found") { temp = "T"; }
  46                         correction = isTrue(temp);
  47
  48                         temp = validParameter.validFile(parameters, "processors", true);                if (temp == "not found") { temp = "1"; }
  49                         convert(temp, processors);
  50
  51                         method = validParameter.validFile(parameters, "method", false);         if (method == "not found") { method = "bellerophon"; }
  52
  53                 }
  54         }
  55         catch(exception& e) {
  56                 errorOut(e, "ChimeraSeqsCommand", "ChimeraSeqsCommand");
  57                 exit(1);
  58         }
  59 }
  60 //**********************************************************************************************************************
  61
  62 void ChimeraSeqsCommand::help(){
  63         try {
  64                 mothurOut("The chimera.seqs command reads a fastafile and creates a sorted priority score list of potentially chimeric sequences (ideally, the sequences should already be aligned).\n");
  65                 mothurOut("The chimera.seqs command parameters are fasta, filter, correction, processors and method.  fasta is required.\n");
  66                 mothurOut("The filter parameter allows you to specify if you would like to apply a 50% soft filter.  The default is false. \n");
  67                 mothurOut("The correction parameter allows you to .....  The default is true. \n");
  68                 mothurOut("The processors parameter allows you to specify how many processors you would like to use.  The default is 1. \n");
  69                 mothurOut("The method parameter allows you to specify the method for finding chimeric sequences.  The default is bellerophon. \n");
  70                 mothurOut("The chimera.seqs command should be in the following format: \n");
  71                 mothurOut("chimera.seqs(fasta=yourFastaFile, filter=yourFilter, correction=yourCorrection, processors=yourProcessors, method=bellerophon) \n");
  72                 mothurOut("Example: chimera.seqs(fasta=AD.align, filter=True, correction=true, processors=2, method=yourMethod) \n");
  73                 mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFastaFile).\n\n");
  74         }
  75         catch(exception& e) {
  76                 errorOut(e, "ChimeraSeqsCommand", "help");
  77                 exit(1);
  78         }
  79 }
  80
  81 //***************************************************************************************************************
  82
  83 ChimeraSeqsCommand::~ChimeraSeqsCommand(){      /*      do nothing      */      }
  84
  85 //***************************************************************************************************************
  86
  87 int ChimeraSeqsCommand::execute(){
  88         try{
  89
  90                 if (abort == true) { return 0; }
  91
  92                 //do soft filter
  93                 if (filter)  {
  94                         string optionString = "fasta=" + fastafile + ", soft=50.0, vertical=F";
  95                         filterSeqs = new FilterSeqsCommand(optionString);
  96                         filterSeqs->execute();
  97                         delete filterSeqs;
  98
  99                         //reset fastafile to filtered file
 100                         fastafile = getRootName(fastafile) + "filter.fasta";
 101                 }
 102
 103                 //read in sequences
 104                 readSeqs();
 105
 106                 //int numSeqs = seqs.size();
 107
 108                 //find average midpoint of seqs
 109                 midpoint = findAverageMidPoint();
 110
 111                 //this should be parallelized
 112                 //generatePreferences();
 113
 114
 115                 //output results to screen
 116                 mothurOutEndLine();
 117                 mothurOut("\t\t"); mothurOutEndLine();
 118                 //mothurOut("Minimum:\t" + toString(startPosition[0]) + "\t" + toString(endPosition[0]) + "\t" + toString(seqLength[0]) + "\t" + toString(ambigBases[0]) + "\t" + toString(longHomoPolymer[0])); mothurOutEndLine();
 119                 //mothurOut("2.5%-tile:\t" + toString(startPosition[ptile0_25]) + "\t" + toString(endPosition[ptile0_25]) + "\t" + toString(seqLength[ptile0_25]) + "\t" + toString(ambigBases[ptile0_25]) + "\t"+ toString(longHomoPolymer[ptile0_25])); mothurOutEndLine();
 120                 //mothurOut("25%-tile:\t" + toString(startPosition[ptile25]) + "\t" + toString(endPosition[ptile25]) + "\t" + toString(seqLength[ptile25]) + "\t" + toString(ambigBases[ptile25]) + "\t" + toString(longHomoPolymer[ptile25])); mothurOutEndLine();
 121                 //mothurOut("Median: \t" + toString(startPosition[ptile50]) + "\t" + toString(endPosition[ptile50]) + "\t" + toString(seqLength[ptile50]) + "\t" + toString(ambigBases[ptile50]) + "\t" + toString(longHomoPolymer[ptile50])); mothurOutEndLine();
 122                 //mothurOut("75%-tile:\t" + toString(startPosition[ptile75]) + "\t" + toString(endPosition[ptile75]) + "\t" + toString(seqLength[ptile75]) + "\t" + toString(ambigBases[ptile75]) + "\t" + toString(longHomoPolymer[ptile75])); mothurOutEndLine();
 123                 //mothurOut("97.5%-tile:\t" + toString(startPosition[ptile97_5]) + "\t" + toString(endPosition[ptile97_5]) + "\t" + toString(seqLength[ptile97_5]) + "\t" + toString(ambigBases[ptile97_5]) + "\t" + toString(longHomoPolymer[ptile97_5])); mothurOutEndLine();
 124                 //mothurOut("Maximum:\t" + toString(startPosition[ptile100]) + "\t" + toString(endPosition[ptile100]) + "\t" + toString(seqLength[ptile100]) + "\t" + toString(ambigBases[ptile100]) + "\t" + toString(longHomoPolymer[ptile100])); mothurOutEndLine();
 125                 //mothurOut("# of Seqs:\t" + toString(numSeqs)); mothurOutEndLine();
 126
 127                 //outSummary.close();
 128                 return 0;
 129         }
 130         catch(exception& e) {
 131                 errorOut(e, "ChimeraSeqsCommand", "execute");
 132                 exit(1);
 133         }
 134 }
 135
 136 //***************************************************************************************************************
 137 void ChimeraSeqsCommand::readSeqs(){
 138         try {
 139                 ifstream inFASTA;
 140                 openInputFile(fastafile, inFASTA);
 141
 142                 //read in seqs and store in vector
 143                 while(!inFASTA.eof()){
 144                         Sequence current(inFASTA);
 145
 146                         seqs.push_back(current);
 147
 148                         gobble(inFASTA);
 149                 }
 150                 inFASTA.close();
 151
 152         }
 153         catch(exception& e) {
 154                 errorOut(e, "ChimeraSeqsCommand", "readSeqs");
 155                 exit(1);
 156         }
 157 }
 158
 159
 160 //***************************************************************************************************************
 161 int ChimeraSeqsCommand::findAverageMidPoint(){
 162         try {
 163                 int totalMids = 0;
 164                 int averageMid = 0;
 165
 166                 //loop through the seqs and find midpoint
 167                 for (int i = 0; i < seqs.size(); i++) {
 168
 169                         //get unaligned sequence
 170                         seqs[i].setUnaligned(seqs[i].getUnaligned());  //if you read an aligned file the unaligned is really aligned, so we need to make sure its unaligned
 171
 172                         string unaligned = seqs[i].getUnaligned();
 173                         string aligned = seqs[i].getAligned();
 174
 175                         //find midpoint of this seq
 176                         int count = 0;
 177                         int thismid = 0;
 178                         for (int j = 0; j < aligned.length(); j++) {
 179
 180                                 thismid++;
 181
 182                                 //if you are part of the unaligned sequence increment
 183                                 if (isalpha(aligned[j])) {  count++;  }
 184
 185                                 //if you have reached the halfway point stop
 186                                 if (count >= (unaligned.length() / 2)) { break; }
 187                         }
 188
 189                         //add this mid to total
 190                         totalMids += thismid;
 191
 192                 }
 193
 194                 averageMid = (totalMids / seqs.size());
 195
 196                 return averageMid;
 197
 198
 199         }
 200         catch(exception& e) {
 201                 errorOut(e, "ChimeraSeqsCommand", "findAverageMidPoint");
 202                 exit(1);
 203         }
 204 }
 205
 206 /***************************************************************************************************************
 207 int ChimeraSeqsCommand::createSparseMatrix(int startLine, int endLine, SparseMatrix* sparse){
 208         try {
 209
 210                 for(int i=startLine; i<endLine; i++){
 211
 212                         for(int j=0;j<i;j++){
 213
 214                                 distCalculator->calcDist(seqs.get(i), seqs.get(j));
 215                                 double dist = distCalculator->getDist();
 216
 217
 218
 219                         }
 220
 221
 222                 return 1;
 223         }
 224         catch(exception& e) {
 225                 errorOut(e, "ChimeraSeqsCommand", "createSparseMatrix");
 226                 exit(1);
 227         }
 228 }
 229 /**************************************************************************************************/
 230