]> git.donarmstrong.com Git - mothur.git/blobdiff - chopseqscommand.cpp
changes to chop.seqs
[mothur.git] / chopseqscommand.cpp
index adf0910589c7d65161027bedae6cc74236e2ef5f..9020d2b4763515e1b6e1ead544247e676b043f41 100644 (file)
@@ -21,7 +21,7 @@ ChopSeqsCommand::ChopSeqsCommand(string option)  {
                
                else {
                        //valid paramters for this command
-                       string Array[] =  {"fasta","end","outputdir","inputdir"};
+                       string Array[] =  {"fasta","end","fromend","outputdir","inputdir"};
                        vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
                        
                        OptionParser parser(option);
@@ -57,13 +57,14 @@ ChopSeqsCommand::ChopSeqsCommand(string option)  {
                        if (fastafile == "not open") { abort = true; }
                        else if (fastafile == "not found") {  m->mothurOut("You must provide a fasta file."); m->mothurOutEndLine(); abort = true; }    
                        
-                       string temp = validParameter.validFile(parameters, "end", false);       
-                       if (temp == "not found") {  m->mothurOut("You must provide an end for the chops.seqs command."); m->mothurOutEndLine(); abort = true; } 
-                       else {  
-                               convert(temp, end);   
-                               if (end < 0) { m->mothurOut("End must be positive."); m->mothurOutEndLine(); abort = true;  }
-                       }
-                       
+                       string temp = validParameter.validFile(parameters, "end", false);       if (temp == "not found") { temp = "0"; } 
+                       convert(temp, end);   
+               
+                       temp = validParameter.validFile(parameters, "fromend", false);          if (temp == "not found") { temp = "0"; } 
+                       convert(temp, fromend);   
+                               
+                       if ((end == 0) && (fromend == 0))  { m->mothurOut("You must provide either end or fromend for the chops.seqs command."); m->mothurOutEndLine(); abort = true;  }
+                       if ((end != 0) && (fromend != 0))  { m->mothurOut("You must provide either end or fromend for the chops.seqs command, not both."); m->mothurOutEndLine(); abort = true;  }
                }
 
        }
@@ -77,8 +78,10 @@ ChopSeqsCommand::ChopSeqsCommand(string option)  {
 void ChopSeqsCommand::help(){
        try {
                m->mothurOut("The chop.seqs command reads a fasta file and outputs a .chop.fasta with sequences trimmed to the end position.\n");
-               m->mothurOut("The chop.seqs command parameters are fasta and end, both are required.\n");
+               m->mothurOut("The chop.seqs command parameters are fasta, end and fromend, fasta is required.\n");
                m->mothurOut("The chop.seqs command should be in the following format: chop.seqs(fasta=yourFasta, end=yourEnd).\n");
+               m->mothurOut("The end parameter allows you to specify an end base position for your sequences, default = 0.\n");
+               m->mothurOut("The fromend parameter allows you to remove the last X bases from the end of the sequence, default = 0.\n");
                m->mothurOut("Example chop.seqs(fasta=amazon.fasta, end=200).\n");
                m->mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n\n");
        }
@@ -96,38 +99,53 @@ int ChopSeqsCommand::execute(){
                if (abort == true) { return 0; }
                
                string outputFileName = outputDir + getRootName(getSimpleName(fastafile)) + "chop.fasta";
+               string outputFileNameAccnos = outputDir + getRootName(getSimpleName(fastafile)) + "chop.accnos";
                
                ofstream out;
                openOutputFile(outputFileName, out);
                
+               ofstream outAcc;
+               openOutputFile(outputFileNameAccnos, outAcc);
+               
                ifstream in;
                openInputFile(fastafile, in);
                
-               while (!in.eof()) {
+               bool wroteAccnos = false;
                
-                       Sequence seq(in, "no align");
+               while (!in.eof()) {
+                       
+                       Sequence seq(in);
                        
-                       if (m->control_pressed) { in.close(); out.close(); remove(outputFileName.c_str()); return 0;  }
+                       if (m->control_pressed) { in.close(); out.close(); outAcc.close(); remove(outputFileName.c_str()); remove(outputFileNameAccnos.c_str()); return 0;  }
                        
                        if (seq.getName() != "") {
-                               string temp = seq.getUnaligned();
-                               
-                               //output sequence name
-                               out << ">" << seq.getName() << endl;
-                               
-                               //if needed trim sequence
-                               if (temp.length() > end) {  temp = temp.substr(0, end);         }
+                               string newSeqString = "";
+                               if (seq.getIsAligned()) { //sequence is aligned
+                                       newSeqString = getChoppedAligned(seq);
+                               }else{
+                                       newSeqString = getChoppedUnaligned(seq);
+                               }
                                
-                               //output trimmed sequence       
-                               out << temp << endl;
+                               //output trimmed sequence
+                               if (newSeqString != "") {
+                                       out << ">" << seq.getName() << endl << newSeqString << endl;
+                               }else{
+                                       outAcc << seq.getName() << endl;
+                                       wroteAccnos = true;
+                               }
                        }
                }
                in.close();
                out.close();
+               outAcc.close();
                
                m->mothurOutEndLine();
                m->mothurOut("Output File Name: "); m->mothurOutEndLine();
                m->mothurOut(outputFileName); m->mothurOutEndLine();    
+               
+               if (wroteAccnos) { m->mothurOut(outputFileNameAccnos); m->mothurOutEndLine();  }
+               else {  remove(outputFileNameAccnos.c_str());  }
+               
                m->mothurOutEndLine();
                
                return 0;               
@@ -138,7 +156,74 @@ int ChopSeqsCommand::execute(){
                exit(1);
        }
 }
+//**********************************************************************************************************************
+string ChopSeqsCommand::getChoppedUnaligned(Sequence seq) {
+       try {
+               string temp = seq.getUnaligned();
+                               
+               //if needed trim sequence
+               if (end != 0) {
+                       if (temp.length() > end) {  temp = temp.substr(0, end);         }
+                       else { temp = "";  }
+               }else { //you are using fromend
+                       if (temp.length() > fromend) { temp = temp.substr(0, (temp.length()-fromend));  }
+                       else {  temp = ""; } //sequence too short
+               }
+
+               return temp;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "ChopSeqsCommand", "getChoppedUnaligned");
+               exit(1);
+       }
+}
+//**********************************************************************************************************************
+string ChopSeqsCommand::getChoppedAligned(Sequence seq) {
+       try {
+               string temp = seq.getAligned();
+               string tempUnaligned = seq.getUnaligned();
+                               
+               //if needed trim sequence
+               if (end != 0) {
+                       if (tempUnaligned.length() > end) { //you have enough bases to remove some
+                               
+                               int stopSpot = 0;
+                               int numBases = 0;
+                               
+                               for (int i = 0; i < temp.length(); i++) {
+                                       if(isalpha(temp[i])) { numBases++; }
 
+                                       if (numBases >= end) { stopSpot = i; break; }
+                               }
+                               
+                               temp = temp.substr(0, stopSpot);                
+                       }else { temp = ""; } //sequence too short
+                       
+               }else { //you are using fromend
+               
+                       if (tempUnaligned.length() > fromend) { 
+                               
+                               int stopSpot = 0;
+                               int numBases = 0;
+                               
+                               for (int i = (temp.length()-1); i >= 0; i--) {
+                                       if(isalpha(temp[i])) { numBases++; }
+
+                                       if (numBases >= fromend) { stopSpot = i; break; }
+                               }
+                               
+                               temp = temp.substr(0, stopSpot);
+                       }
+                       else {  temp = ""; } //sequence too short
+               }
+
+               return temp;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "ChopSeqsCommand", "getChoppedUnaligned");
+               exit(1);
+       }
+}
 //**********************************************************************************************************************