2 * filterseqscommand.cpp
5 * Created by Thomas Ryabin on 5/4/09.
6 * Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
10 #include "filterseqscommand.h"
11 #include "sequence.hpp"
13 /**************************************************************************************/
15 FilterSeqsCommand::FilterSeqsCommand(string option){
19 //allow user to run help
20 if(option == "help") { help(); abort = true; }
23 //valid paramters for this command
24 string Array[] = {"fasta", "trump", "soft", "hard", "vertical"};
25 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
27 OptionParser parser(option);
28 map<string,string> parameters = parser.getParameters();
30 ValidParameters validParameter;
32 //check to make sure all parameters are valid for command
33 for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
34 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
37 //check for required parameters
38 fastafile = validParameter.validFile(parameters, "fasta", true);
39 if (fastafile == "not found") { mothurOut("fasta is a required parameter for the filter.seqs command."); mothurOutEndLine(); abort = true; }
40 else if (fastafile == "not open") { abort = true; }
42 //check for optional parameter and set defaults
43 // ...at some point should added some additional type checking...
46 temp = validParameter.validFile(parameters, "trump", false); if (temp == "not found") { temp = "*"; }
49 temp = validParameter.validFile(parameters, "soft", false); if (temp == "not found") { soft = 0; }
50 else { soft = (float)atoi(temp.c_str()) / 100.0; }
52 hard = validParameter.validFile(parameters, "hard", true); if (hard == "not found") { hard = ""; }
53 else if (hard == "not open") { abort = true; }
55 vertical = validParameter.validFile(parameters, "vertical", false); if (vertical == "not found") { vertical = "T"; }
61 if (soft != 0) { F.setSoft(soft); }
62 if (trump != '*') { F.setTrump(trump); }
70 errorOut(e, "FilterSeqsCommand", "FilterSeqsCommand");
75 //**********************************************************************************************************************
77 void FilterSeqsCommand::help(){
79 mothurOut("The filter.seqs command reads a file containing sequences and creates a .filter and .filter.fasta file.\n");
80 mothurOut("The filter.seqs command parameters are fasta, trump, soft, hard and vertical. \n");
81 mothurOut("The fasta parameter is required.\n");
82 mothurOut("The trump parameter .... The default is ...\n");
83 mothurOut("The soft parameter .... The default is ....\n");
84 mothurOut("The hard parameter .... The default is ....\n");
85 mothurOut("The vertical parameter .... The default is T.\n");
86 mothurOut("The filter.seqs command should be in the following format: \n");
87 mothurOut("filter.seqs(fasta=yourFastaFile, trump=yourTrump, soft=yourSoft, hard=yourHard, vertical=yourVertical) \n");
88 mothurOut("Example filter.seqs(fasta=abrecovery.fasta, trump=..., soft=..., hard=..., vertical=T).\n");
89 mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n\n");
93 errorOut(e, "FilterSeqsCommand", "help");
98 /**************************************************************************************/
100 int FilterSeqsCommand::execute() {
103 if (abort == true) { return 0; }
106 openInputFile(fastafile, inFASTA);
108 Sequence testSeq(inFASTA);
109 alignmentLength = testSeq.getAlignLength();
112 F.setLength(alignmentLength);
114 if(soft != 0 || isTrue(vertical)){
118 if(hard.compare("") != 0) { F.doHard(hard); }
119 else { F.setFilter(string(alignmentLength, '1')); }
121 if(trump != '*' || isTrue(vertical) || soft != 0){
122 while(!inFASTA.eof()){ //read through and create the filter...
123 Sequence seq(inFASTA);
124 if (seq.getName() != "") {
125 if(trump != '*'){ F.doTrump(seq); }
126 if(isTrue(vertical) || soft != 0){ F.getFreqs(seq); }
134 F.setNumSeqs(numSeqs);
137 if(isTrue(vertical) == 1) { F.doVertical(); }
138 if(soft != 0) { F.doSoft(); }
140 filter = F.getFilter();
143 string filterFile = getRootName(fastafile) + "filter";
144 openOutputFile(filterFile, outFilter);
145 outFilter << filter << endl;
149 openInputFile(fastafile, inFasta2);
150 string filteredFasta = getRootName(fastafile) + "filter.fasta";
152 openOutputFile(filteredFasta, outFASTA);
155 while(!inFasta2.eof()){
156 Sequence seq(inFasta2);
157 if (seq.getName() != "") {
158 string align = seq.getAligned();
159 string filterSeq = "";
161 for(int j=0;j<alignmentLength;j++){
162 if(filter[j] == '1'){
163 filterSeq += align[j];
167 outFASTA << '>' << seq.getName() << endl << filterSeq << endl;
176 int filteredLength = 0;
177 for(int i=0;i<alignmentLength;i++){
178 if(filter[i] == '1'){ filteredLength++; }
182 mothurOut("Length of filtered alignment: " + toString(filteredLength)); mothurOutEndLine();
183 mothurOut("Number of columns removed: " + toString((alignmentLength-filteredLength))); mothurOutEndLine();
184 mothurOut("Length of the original alignment: " + toString(alignmentLength)); mothurOutEndLine();
185 mothurOut("Number of sequences used to construct filter: " + toString(numSeqs)); mothurOutEndLine();
190 catch(exception& e) {
191 errorOut(e, "FilterSeqsCommand", "execute");
196 /**************************************************************************************/