5 * Created by westcott on 5/10/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "chopseqscommand.h"
11 #include "sequence.hpp"
13 //**********************************************************************************************************************
14 vector<string> ChopSeqsCommand::setParameters(){
16 CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pfasta);
17 CommandParameter pnumbases("numbases", "Number", "", "0", "", "", "",false,true); parameters.push_back(pnumbases);
18 CommandParameter pcountgaps("countgaps", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pcountgaps);
19 CommandParameter pshort("short", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pshort);
20 CommandParameter pkeep("keep", "Multiple", "front-back", "front", "", "", "",false,false); parameters.push_back(pkeep);
21 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
22 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
24 vector<string> myArray;
25 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
29 m->errorOut(e, "ChopSeqsCommand", "setParameters");
33 //**********************************************************************************************************************
34 string ChopSeqsCommand::getHelpString(){
36 string helpString = "";
37 helpString += "The chop.seqs command reads a fasta file and outputs a .chop.fasta containing the trimmed sequences. Note: If a sequence is completely 'chopped', an accnos file will be created with the names of the sequences removed. \n";
38 helpString += "The chop.seqs command parameters are fasta, numbases, countgaps and keep. fasta is required unless you have a valid current fasta file. numbases is required.\n";
39 helpString += "The chop.seqs command should be in the following format: chop.seqs(fasta=yourFasta, numbases=yourNum, keep=yourKeep).\n";
40 helpString += "The numbases parameter allows you to specify the number of bases you want to keep.\n";
41 helpString += "The keep parameter allows you to specify whether you want to keep the front or the back of your sequence, default=front.\n";
42 helpString += "The countgaps parameter allows you to specify whether you want to count gaps as bases, default=false.\n";
43 helpString += "The short parameter allows you to specify you want to keep sequences that are too short to chop, default=false.\n";
44 helpString += "For example, if you ran chop.seqs with numbases=200 and short=t, if a sequence had 100 bases mothur would keep the sequence rather than eliminate it.\n";
45 helpString += "Example chop.seqs(fasta=amazon.fasta, numbases=200, keep=front).\n";
46 helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n";
50 m->errorOut(e, "ChopSeqsCommand", "getHelpString");
55 //**********************************************************************************************************************
56 ChopSeqsCommand::ChopSeqsCommand(){
58 abort = true; calledHelp = true;
60 vector<string> tempOutNames;
61 outputTypes["fasta"] = tempOutNames;
62 outputTypes["accnos"] = tempOutNames;
65 m->errorOut(e, "ChopSeqsCommand", "ChopSeqsCommand");
69 //**********************************************************************************************************************
70 ChopSeqsCommand::ChopSeqsCommand(string option) {
72 abort = false; calledHelp = false;
74 //allow user to run help
75 if(option == "help") { help(); abort = true; calledHelp = true; }
78 vector<string> myArray = setParameters();
80 OptionParser parser(option);
81 map<string,string> parameters = parser.getParameters();
83 ValidParameters validParameter;
84 map<string,string>::iterator it;
86 //check to make sure all parameters are valid for command
87 for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
88 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
91 //initialize outputTypes
92 vector<string> tempOutNames;
93 outputTypes["fasta"] = tempOutNames;
94 outputTypes["accnos"] = tempOutNames;
96 //if the user changes the input directory command factory will send this info to us in the output parameter
97 string inputDir = validParameter.validFile(parameters, "inputdir", false);
98 if (inputDir == "not found"){ inputDir = ""; }
101 it = parameters.find("fasta");
102 //user has given a template file
103 if(it != parameters.end()){
104 path = m->hasPath(it->second);
105 //if the user has not given a path then, add inputdir. else leave path alone.
106 if (path == "") { parameters["fasta"] = inputDir + it->second; }
110 //check for required parameters
111 fastafile = validParameter.validFile(parameters, "fasta", true);
112 if (fastafile == "not open") { abort = true; }
113 else if (fastafile == "not found") { //if there is a current fasta file, use it
114 fastafile = m->getFastaFile();
115 if (fastafile != "") { m->mothurOut("Using " + fastafile + " as input file for the fasta parameter."); m->mothurOutEndLine(); }
116 else { m->mothurOut("You have no current fastafile and the fasta parameter is required."); m->mothurOutEndLine(); abort = true; }
119 //if the user changes the output directory command factory will send this info to us in the output parameter
120 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = m->hasPath(fastafile); }
122 string temp = validParameter.validFile(parameters, "numbases", false); if (temp == "not found") { temp = "0"; }
123 convert(temp, numbases);
125 temp = validParameter.validFile(parameters, "countgaps", false); if (temp == "not found") { temp = "f"; }
126 countGaps = m->isTrue(temp);
128 temp = validParameter.validFile(parameters, "short", false); if (temp == "not found") { temp = "f"; }
129 Short = m->isTrue(temp);
131 keep = validParameter.validFile(parameters, "keep", false); if (keep == "not found") { keep = "front"; }
133 if (numbases == 0) { m->mothurOut("You must provide the number of bases you want to keep for the chops.seqs command."); m->mothurOutEndLine(); abort = true; }
137 catch(exception& e) {
138 m->errorOut(e, "ChopSeqsCommand", "ChopSeqsCommand");
142 //**********************************************************************************************************************
144 int ChopSeqsCommand::execute(){
147 if (abort == true) { if (calledHelp) { return 0; } return 2; }
149 string outputFileName = outputDir + m->getRootName(m->getSimpleName(fastafile)) + "chop.fasta";
150 string outputFileNameAccnos = outputDir + m->getRootName(m->getSimpleName(fastafile)) + "chop.accnos";
153 m->openOutputFile(outputFileName, out);
156 m->openOutputFile(outputFileNameAccnos, outAcc);
159 m->openInputFile(fastafile, in);
161 bool wroteAccnos = false;
167 if (m->control_pressed) { outputTypes.clear(); in.close(); out.close(); outAcc.close(); remove(outputFileName.c_str()); remove(outputFileNameAccnos.c_str()); return 0; }
169 if (seq.getName() != "") {
170 string newSeqString = getChopped(seq);
172 //output trimmed sequence
173 if (newSeqString != "") {
174 out << ">" << seq.getName() << endl << newSeqString << endl;
176 outAcc << seq.getName() << endl;
185 m->mothurOutEndLine();
186 m->mothurOut("Output File Name: "); m->mothurOutEndLine();
187 m->mothurOut(outputFileName); m->mothurOutEndLine(); outputNames.push_back(outputFileName); outputTypes["fasta"].push_back(outputFileName);
189 if (wroteAccnos) { m->mothurOut(outputFileNameAccnos); m->mothurOutEndLine(); outputNames.push_back(outputFileNameAccnos); outputTypes["accnos"].push_back(outputFileNameAccnos); }
190 else { remove(outputFileNameAccnos.c_str()); }
192 m->mothurOutEndLine();
194 //set fasta file as new current fastafile
196 itTypes = outputTypes.find("fasta");
197 if (itTypes != outputTypes.end()) {
198 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setFastaFile(current); }
201 if (wroteAccnos) { //set accnos file as new current accnosfile
202 itTypes = outputTypes.find("accnos");
203 if (itTypes != outputTypes.end()) {
204 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setAccnosFile(current); }
212 catch(exception& e) {
213 m->errorOut(e, "ChopSeqsCommand", "execute");
217 //**********************************************************************************************************************
218 string ChopSeqsCommand::getChopped(Sequence seq) {
220 string temp = seq.getAligned();
221 string tempUnaligned = seq.getUnaligned();
224 //if needed trim sequence
225 if (keep == "front") {//you want to keep the beginning
226 int tempLength = temp.length();
228 if (tempLength > numbases) { //you have enough bases to remove some
231 int numBasesCounted = 0;
233 for (int i = 0; i < temp.length(); i++) {
235 if (toupper(temp[i]) == 'N') { temp[i] == '.'; }
239 if (numBasesCounted >= numbases) { stopSpot = i; break; }
242 if (stopSpot == 0) { temp = ""; }
243 else { temp = temp.substr(0, stopSpot); }
246 if (!Short) { temp = ""; } //sequence too short
248 }else { //you are keeping the back
249 int tempLength = temp.length();
250 if (tempLength > numbases) { //you have enough bases to remove some
253 int numBasesCounted = 0;
255 for (int i = (temp.length()-1); i >= 0; i--) {
257 if (toupper(temp[i]) == 'N') { temp[i] == '.'; }
261 if (numBasesCounted >= numbases) { stopSpot = i; break; }
264 if (stopSpot == 0) { temp = ""; }
265 else { temp = temp.substr(stopSpot+1); }
267 if (!Short) { temp = ""; } //sequence too short
273 //if needed trim sequence
274 if (keep == "front") {//you want to keep the beginning
275 int tempLength = tempUnaligned.length();
277 if (tempLength > numbases) { //you have enough bases to remove some
280 int numBasesCounted = 0;
282 for (int i = 0; i < temp.length(); i++) {
284 if (toupper(temp[i]) == 'N') {
287 if (tempLength < numbases) { stopSpot = 0; break; }
290 if(isalpha(temp[i])) { numBasesCounted++; }
292 if (numBasesCounted >= numbases) { stopSpot = i; break; }
295 if (stopSpot == 0) { temp = ""; }
296 else { temp = temp.substr(0, stopSpot); }
299 if (!Short) { temp = ""; } //sequence too short
301 }else { //you are keeping the back
302 int tempLength = tempUnaligned.length();
303 if (tempLength > numbases) { //you have enough bases to remove some
306 int numBasesCounted = 0;
308 for (int i = (temp.length()-1); i >= 0; i--) {
310 if (toupper(temp[i]) == 'N') {
313 if (tempLength < numbases) { stopSpot = 0; break; }
316 if(isalpha(temp[i])) { numBasesCounted++; }
318 if (numBasesCounted >= numbases) { stopSpot = i; break; }
321 if (stopSpot == 0) { temp = ""; }
322 else { temp = temp.substr(stopSpot+1); }
324 if (!Short) { temp = ""; } //sequence too short
331 catch(exception& e) {
332 m->errorOut(e, "ChopSeqsCommand", "getChopped");
336 //**********************************************************************************************************************