#include "trimseqscommand.h"
#include "needlemanoverlap.hpp"
+//**********************************************************************************************************************
+
+vector<string> TrimSeqsCommand::getValidParameters(){
+ try {
+ string Array[] = {"fasta", "flip", "oligos", "maxambig", "maxhomop","minlength", "maxlength", "qfile",
+ "qthreshold", "qwindowaverage", "qstepsize", "qwindowsize", "qaverage", "rollaverage",
+ "keepfirst", "removelast",
+ "allfiles", "qtrim","tdiffs", "pdiffs", "bdiffs", "processors", "outputdir","inputdir"};
+ vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
+ return myArray;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "TrimSeqsCommand", "getValidParameters");
+ exit(1);
+ }
+}
+
+//**********************************************************************************************************************
+
+TrimSeqsCommand::TrimSeqsCommand(){
+ try {
+ abort = true; calledHelp = true;
+ vector<string> tempOutNames;
+ outputTypes["fasta"] = tempOutNames;
+ outputTypes["qfile"] = tempOutNames;
+ outputTypes["group"] = tempOutNames;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "TrimSeqsCommand", "TrimSeqsCommand");
+ exit(1);
+ }
+}
+
+//**********************************************************************************************************************
+
+vector<string> TrimSeqsCommand::getRequiredParameters(){
+ try {
+ string Array[] = {"fasta"};
+ vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
+ return myArray;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "TrimSeqsCommand", "getRequiredParameters");
+ exit(1);
+ }
+}
+
+//**********************************************************************************************************************
+
+vector<string> TrimSeqsCommand::getRequiredFiles(){
+ try {
+ vector<string> myArray;
+ return myArray;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "TrimSeqsCommand", "getRequiredFiles");
+ exit(1);
+ }
+}
+
//***************************************************************************************************************
TrimSeqsCommand::TrimSeqsCommand(string option) {
try {
- abort = false;
+ abort = false; calledHelp = false;
+ comboStarts = 0;
//allow user to run help
- if(option == "help") { help(); abort = true; }
+ if(option == "help") { help(); abort = true; calledHelp = true; }
else {
//valid paramters for this command
- string AlignArray[] = {"fasta", "flip", "oligos", "maxambig", "maxhomop", "minlength", "maxlength", "qfile",
- "qthreshold", "qaverage", "allfiles", "qtrim","tdiffs", "pdiffs", "bdiffs", "processors", "outputdir","inputdir"};
+ string AlignArray[] = { "fasta", "flip", "oligos", "maxambig", "maxhomop", "minlength", "maxlength", "qfile",
+ "qthreshold", "qwindowaverage", "qstepsize", "qwindowsize", "qaverage", "rollaverage",
+ "keepfirst", "removelast",
+ "allfiles", "qtrim","tdiffs", "pdiffs", "bdiffs", "processors", "outputdir","inputdir"};
vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
}
+ //initialize outputTypes
+ vector<string> tempOutNames;
+ outputTypes["fasta"] = tempOutNames;
+ outputTypes["qfile"] = tempOutNames;
+ outputTypes["group"] = tempOutNames;
+
//if the user changes the input directory command factory will send this info to us in the output parameter
string inputDir = validParameter.validFile(parameters, "inputdir", false);
if (inputDir == "not found"){ inputDir = ""; }
it = parameters.find("fasta");
//user has given a template file
if(it != parameters.end()){
- path = hasPath(it->second);
+ path = m->hasPath(it->second);
//if the user has not given a path then, add inputdir. else leave path alone.
if (path == "") { parameters["fasta"] = inputDir + it->second; }
}
it = parameters.find("oligos");
//user has given a template file
if(it != parameters.end()){
- path = hasPath(it->second);
+ path = m->hasPath(it->second);
//if the user has not given a path then, add inputdir. else leave path alone.
if (path == "") { parameters["oligos"] = inputDir + it->second; }
}
it = parameters.find("qfile");
//user has given a template file
if(it != parameters.end()){
- path = hasPath(it->second);
+ path = m->hasPath(it->second);
//if the user has not given a path then, add inputdir. else leave path alone.
if (path == "") { parameters["qfile"] = inputDir + it->second; }
}
+
}
//check for required parameters
fastaFile = validParameter.validFile(parameters, "fasta", true);
- if (fastaFile == "not found") { m->mothurOut("fasta is a required parameter for the screen.seqs command."); m->mothurOutEndLine(); abort = true; }
+ if (fastaFile == "not found") { m->mothurOut("fasta is a required parameter for the trim.seqs command."); m->mothurOutEndLine(); abort = true; }
else if (fastaFile == "not open") { abort = true; }
//if the user changes the output directory command factory will send this info to us in the output parameter
outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){
outputDir = "";
- outputDir += hasPath(fastaFile); //if user entered a file with a path then preserve it
+ outputDir += m->hasPath(fastaFile); //if user entered a file with a path then preserve it
}
+
//check for optional parameter and set defaults
// ...at some point should added some additional type checking...
string temp;
temp = validParameter.validFile(parameters, "flip", false);
if (temp == "not found"){ flip = 0; }
- else if(isTrue(temp)) { flip = 1; }
+ else if(m->isTrue(temp)) { flip = 1; }
temp = validParameter.validFile(parameters, "oligos", true);
if (temp == "not found"){ oligoFile = ""; }
else if(temp == "not open"){ abort = true; }
else { oligoFile = temp; }
+
temp = validParameter.validFile(parameters, "maxambig", false); if (temp == "not found") { temp = "-1"; }
convert(temp, maxAmbig);
temp = validParameter.validFile(parameters, "maxlength", false); if (temp == "not found") { temp = "0"; }
convert(temp, maxLength);
- temp = validParameter.validFile(parameters, "tdiffs", false); if (temp == "not found") { temp = "0"; }
- convert(temp, tdiffs);
-
temp = validParameter.validFile(parameters, "bdiffs", false); if (temp == "not found") { temp = "0"; }
convert(temp, bdiffs);
temp = validParameter.validFile(parameters, "pdiffs", false); if (temp == "not found") { temp = "0"; }
convert(temp, pdiffs);
+ temp = validParameter.validFile(parameters, "tdiffs", false); if (temp == "not found") { int tempTotal = pdiffs + bdiffs; temp = toString(tempTotal); }
+ convert(temp, tdiffs);
+
if(tdiffs == 0){ tdiffs = bdiffs + pdiffs; }
temp = validParameter.validFile(parameters, "qfile", true);
temp = validParameter.validFile(parameters, "qthreshold", false); if (temp == "not found") { temp = "0"; }
convert(temp, qThreshold);
- temp = validParameter.validFile(parameters, "qtrim", false); if (temp == "not found") { temp = "F"; }
- qtrim = isTrue(temp);
+ temp = validParameter.validFile(parameters, "qtrim", false); if (temp == "not found") { temp = "t"; }
+ qtrim = m->isTrue(temp);
+
+ temp = validParameter.validFile(parameters, "rollaverage", false); if (temp == "not found") { temp = "0"; }
+ convert(temp, qRollAverage);
+
+ temp = validParameter.validFile(parameters, "qwindowaverage", false);if (temp == "not found") { temp = "0"; }
+ convert(temp, qWindowAverage);
+
+ temp = validParameter.validFile(parameters, "qwindowsize", false); if (temp == "not found") { temp = "50"; }
+ convert(temp, qWindowSize);
+
+ temp = validParameter.validFile(parameters, "qstepsize", false); if (temp == "not found") { temp = "1"; }
+ convert(temp, qWindowStep);
temp = validParameter.validFile(parameters, "qaverage", false); if (temp == "not found") { temp = "0"; }
convert(temp, qAverage);
+
+ temp = validParameter.validFile(parameters, "keepfirst", false); if (temp == "not found") { temp = "0"; }
+ convert(temp, keepFirst);
+
+ temp = validParameter.validFile(parameters, "removelast", false); if (temp == "not found") { temp = "0"; }
+ convert(temp, removeLast);
temp = validParameter.validFile(parameters, "allfiles", false); if (temp == "not found") { temp = "F"; }
- allFiles = isTrue(temp);
+ allFiles = m->isTrue(temp);
- temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = "1"; }
+ temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found") { temp = "1"; }
convert(temp, processors);
- if(allFiles && oligoFile == ""){
- m->mothurOut("You selected allfiles, but didn't enter an oligos file. Ignoring the allfiles request."); m->mothurOutEndLine();
+
+ if(allFiles && (oligoFile == "")){
+ m->mothurOut("You selected allfiles, but didn't enter an oligos. Ignoring the allfiles request."); m->mothurOutEndLine();
}
if((qAverage != 0 && qThreshold != 0) && qFileName == ""){
m->mothurOut("You didn't provide a quality file name, quality criteria will be ignored."); m->mothurOutEndLine();
exit(1);
}
}
+
//**********************************************************************************************************************
void TrimSeqsCommand::help(){
try {
- m->mothurOut("The trim.seqs command reads a fastaFile and creates .....\n");
- m->mothurOut("The trim.seqs command parameters are fasta, flip, oligos, maxambig, maxhomop, minlength, maxlength, qfile, qthreshold, qaverage, diffs, qtrim and allfiles.\n");
+ m->mothurOut("The trim.seqs command reads a fastaFile and creates 2 new fasta files, .trim.fasta and scrap.fasta, as well as group files if you provide and oligos file.\n");
+ m->mothurOut("The .trim.fasta contains sequences that meet your requirements, and the .scrap.fasta contains those which don't.\n");
+ m->mothurOut("The trim.seqs command parameters are fasta, flip, oligos, maxambig, maxhomop, minlength, maxlength, qfile, qthreshold, qaverage, diffs, qtrim, keepfirst, removelast and allfiles.\n");
m->mothurOut("The fasta parameter is required.\n");
- m->mothurOut("The flip parameter .... The default is 0.\n");
- m->mothurOut("The oligos parameter .... The default is "".\n");
- m->mothurOut("The maxambig parameter .... The default is -1.\n");
- m->mothurOut("The maxhomop parameter .... The default is 0.\n");
- m->mothurOut("The minlength parameter .... The default is 0.\n");
- m->mothurOut("The maxlength parameter .... The default is 0.\n");
- m->mothurOut("The tdiffs parameter is used to specify the total number of differences allowed in the sequence. The default is 0.\n");
+ m->mothurOut("The flip parameter will output the reverse compliment of your trimmed sequence. The default is false.\n");
+ m->mothurOut("The oligos parameter allows you to provide an oligos file.\n");
+ m->mothurOut("The maxambig parameter allows you to set the maximum number of ambigious bases allowed. The default is -1.\n");
+ m->mothurOut("The maxhomop parameter allows you to set a maximum homopolymer length. \n");
+ m->mothurOut("The minlength parameter allows you to set and minimum sequence length. \n");
+ m->mothurOut("The maxlength parameter allows you to set and maximum sequence length. \n");
+ m->mothurOut("The tdiffs parameter is used to specify the total number of differences allowed in the sequence. The default is pdiffs + bdiffs.\n");
m->mothurOut("The bdiffs parameter is used to specify the number of differences allowed in the barcode. The default is 0.\n");
m->mothurOut("The pdiffs parameter is used to specify the number of differences allowed in the primer. The default is 0.\n");
- m->mothurOut("The qfile parameter .....\n");
- m->mothurOut("The qthreshold parameter .... The default is 0.\n");
- m->mothurOut("The qaverage parameter .... The default is 0.\n");
- m->mothurOut("The allfiles parameter .... The default is F.\n");
- m->mothurOut("The qtrim parameter .... The default is F.\n");
+ m->mothurOut("The qfile parameter allows you to provide a quality file.\n");
+ m->mothurOut("The qthreshold parameter allows you to set a minimum quality score allowed. \n");
+ m->mothurOut("The qaverage parameter allows you to set a minimum average quality score allowed. \n");
+ m->mothurOut("The qwindowsize parameter allows you to set a number of bases in a window. Default=50.\n");
+ m->mothurOut("The qwindowaverage parameter allows you to set a minimum average quality score allowed over a window. \n");
+ m->mothurOut("The rollaverage parameter allows you to set a minimum rolling average quality score allowed over a window. \n");
+ m->mothurOut("The qstepsize parameter allows you to set a number of bases to move the window over. Default=1.\n");
+ m->mothurOut("The allfiles parameter will create separate group and fasta file for each grouping. The default is F.\n");
+ m->mothurOut("The qtrim parameter will trim sequence from the point that they fall below the qthreshold and put it in the .trim file if set to true. The default is T.\n");
+ m->mothurOut("The keepfirst parameter trims the sequence to the first keepfirst number of bases after the barcode or primers are removed, before the sequence is checked to see if it meets the other requirements. \n");
+ m->mothurOut("The removelast removes the last removelast number of bases after the barcode or primers are removed, before the sequence is checked to see if it meets the other requirements.\n");
m->mothurOut("The trim.seqs command should be in the following format: \n");
m->mothurOut("trim.seqs(fasta=yourFastaFile, flip=yourFlip, oligos=yourOligos, maxambig=yourMaxambig, \n");
m->mothurOut("maxhomop=yourMaxhomop, minlength=youMinlength, maxlength=yourMaxlength) \n");
int TrimSeqsCommand::execute(){
try{
- if (abort == true) { return 0; }
+ if (abort == true) { if (calledHelp) { return 0; } return 2; }
numFPrimers = 0; //this needs to be initialized
numRPrimers = 0;
+ vector<vector<string> > fastaFileNames;
+ vector<vector<string> > qualFileNames;
+
+ string trimSeqFile = outputDir + m->getRootName(m->getSimpleName(fastaFile)) + "trim.fasta";
+ outputNames.push_back(trimSeqFile); outputTypes["fasta"].push_back(trimSeqFile);
+
+ string scrapSeqFile = outputDir + m->getRootName(m->getSimpleName(fastaFile)) + "scrap.fasta";
+ outputNames.push_back(scrapSeqFile); outputTypes["fasta"].push_back(scrapSeqFile);
+
+ string trimQualFile = outputDir + m->getRootName(m->getSimpleName(fastaFile)) + "trim.qual";
+ string scrapQualFile = outputDir + m->getRootName(m->getSimpleName(fastaFile)) + "scrap.qual";
+ if (qFileName != "") {
+ outputNames.push_back(trimQualFile);
+ outputNames.push_back(scrapQualFile);
+ outputTypes["qfile"].push_back(trimQualFile);
+ outputTypes["qfile"].push_back(scrapQualFile);
+ }
- string trimSeqFile = outputDir + getRootName(getSimpleName(fastaFile)) + "trim.fasta";
- outputNames.push_back(trimSeqFile);
- string scrapSeqFile = outputDir + getRootName(getSimpleName(fastaFile)) + "scrap.fasta";
- outputNames.push_back(scrapSeqFile);
- string groupFile = outputDir + getRootName(getSimpleName(fastaFile)) + "groups";
-
- vector<string> fastaFileNames;
+ string outputGroupFileName;
if(oligoFile != ""){
- outputNames.push_back(groupFile);
- getOligos(fastaFileNames);
+ outputGroupFileName = outputDir + m->getRootName(m->getSimpleName(fastaFile)) + "groups";
+ outputNames.push_back(outputGroupFileName); outputTypes["group"].push_back(outputGroupFileName);
+ getOligos(fastaFileNames, qualFileNames);
}
-
- if(qFileName != "") { setLines(qFileName, qLines); }
+ vector<unsigned long int> fastaFilePos;
+ vector<unsigned long int> qFilePos;
+
+ setLines(fastaFile, qFileName, fastaFilePos, qFilePos);
+
+ for (int i = 0; i < (fastaFilePos.size()-1); i++) {
+ lines.push_back(new linePair(fastaFilePos[i], fastaFilePos[(i+1)]));
+ if (qFileName != "") { qLines.push_back(new linePair(qFilePos[i], qFilePos[(i+1)])); }
+ }
+ if(qFileName == "") { qLines = lines; } //files with duds
#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
if(processors == 1){
- ifstream inFASTA;
- openInputFile(fastaFile, inFASTA);
- int numSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
- inFASTA.close();
-
- lines.push_back(new linePair(0, numSeqs));
-
- driverCreateTrim(fastaFile, qFileName, trimSeqFile, scrapSeqFile, groupFile, fastaFileNames, lines[0], lines[0]);
-
- for (int j = 0; j < fastaFileNames.size(); j++) {
- rename((fastaFileNames[j] + toString(getpid()) + ".temp").c_str(), fastaFileNames[j].c_str());
- }
-
+ driverCreateTrim(fastaFile, qFileName, trimSeqFile, scrapSeqFile, trimQualFile, scrapQualFile, outputGroupFileName, fastaFileNames, qualFileNames, lines[0], qLines[0]);
}else{
- setLines(fastaFile, lines);
- if(qFileName == "") { qLines = lines; }
-
- createProcessesCreateTrim(fastaFile, qFileName, trimSeqFile, scrapSeqFile, groupFile, fastaFileNames);
-
- rename((trimSeqFile + toString(processIDS[0]) + ".temp").c_str(), trimSeqFile.c_str());
- rename((scrapSeqFile + toString(processIDS[0]) + ".temp").c_str(), scrapSeqFile.c_str());
- rename((groupFile + toString(processIDS[0]) + ".temp").c_str(), groupFile.c_str());
- for (int j = 0; j < fastaFileNames.size(); j++) {
- rename((fastaFileNames[j] + toString(processIDS[0]) + ".temp").c_str(), fastaFileNames[j].c_str());
- }
- //append files
- for(int i=1;i<processors;i++){
- appendFiles((trimSeqFile + toString(processIDS[i]) + ".temp"), trimSeqFile);
- remove((trimSeqFile + toString(processIDS[i]) + ".temp").c_str());
- appendFiles((scrapSeqFile + toString(processIDS[i]) + ".temp"), scrapSeqFile);
- remove((scrapSeqFile + toString(processIDS[i]) + ".temp").c_str());
- appendFiles((groupFile + toString(processIDS[i]) + ".temp"), groupFile);
- remove((groupFile + toString(processIDS[i]) + ".temp").c_str());
- for (int j = 0; j < fastaFileNames.size(); j++) {
- appendFiles((fastaFileNames[j] + toString(processIDS[i]) + ".temp"), fastaFileNames[j]);
- remove((fastaFileNames[j] + toString(processIDS[i]) + ".temp").c_str());
+ createProcessesCreateTrim(fastaFile, qFileName, trimSeqFile, scrapSeqFile, trimQualFile, scrapQualFile, outputGroupFileName, fastaFileNames, qualFileNames);
+ }
+ #else
+ driverCreateTrim(fastaFile, qFileName, trimSeqFile, scrapSeqFile, trimQualFile, scrapQualFile, outputGroupFileName, fastaFileNames, qualFileNames, lines[0], qLines[0]);
+ #endif
+
+ if (m->control_pressed) { return 0; }
+
+ if(allFiles){
+ map<string, string> uniqueFastaNames;// so we don't add the same groupfile multiple times
+ map<string, string>::iterator it;
+ set<string> namesToRemove;
+ for(int i=0;i<fastaFileNames.size();i++){
+ for(int j=0;j<fastaFileNames[0].size();j++){
+ if (fastaFileNames[i][j] != "") {
+ if(m->isBlank(fastaFileNames[i][j])){
+ remove(fastaFileNames[i][j].c_str());
+ namesToRemove.insert(fastaFileNames[i][j]);
+
+ if(qFileName != ""){
+ remove(qualFileNames[i][j].c_str());
+ namesToRemove.insert(qualFileNames[i][j]);
+ }
+ }else{
+ it = uniqueFastaNames.find(fastaFileNames[i][j]);
+ if (it == uniqueFastaNames.end()) {
+ uniqueFastaNames[fastaFileNames[i][j]] = barcodeNameVector[i];
+ }
}
}
}
+ }
+
+ //remove names for outputFileNames, just cleans up the output
+ vector<string> outputNames2;
+ for(int i = 0; i < outputNames.size(); i++) { if (namesToRemove.count(outputNames[i]) == 0) { outputNames2.push_back(outputNames[i]); } }
+ outputNames = outputNames2;
+
+ for (it = uniqueFastaNames.begin(); it != uniqueFastaNames.end(); it++) {
+ ifstream in;
+ m->openInputFile(it->first, in);
- if (m->control_pressed) { return 0; }
- #else
- ifstream inFASTA;
- openInputFile(fastafileNames[s], inFASTA);
- numSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
- inFASTA.close();
-
- lines.push_back(new linePair(0, numSeqs));
-
- driverCreateSummary(fastafile, qFileName, trimSeqFile, scrapSeqFile, groupFile, fastaFileNames, lines[0], lines[0]);
+ ofstream out;
+ string thisGroupName = outputDir + m->getRootName(m->getSimpleName(it->first)) + "groups";
+ outputNames.push_back(thisGroupName); outputTypes["group"].push_back(thisGroupName);
+ m->openOutputFile(thisGroupName, out);
- if (m->control_pressed) { return 0; }
- #endif
-
-
- for(int i=0;i<fastaFileNames.size();i++){
- ifstream inFASTA;
- string seqName;
- openInputFile(getRootName(fastaFile) + groupVector[i] + ".fasta", inFASTA);
- ofstream outGroups;
- openOutputFile(outputDir + getRootName(getSimpleName(fastaFile)) + groupVector[i] + ".groups", outGroups);
- outputNames.push_back(outputDir + getRootName(getSimpleName(fastaFile)) + groupVector[i] + ".groups");
-
- while(!inFASTA.eof()){
- if(inFASTA.get() == '>'){
- inFASTA >> seqName;
- outGroups << seqName << '\t' << groupVector[i] << endl;
+ while (!in.eof()){
+ if (m->control_pressed) { break; }
+
+ Sequence currSeq(in); m->gobble(in);
+ out << currSeq.getName() << '\t' << it->second << endl;
}
- while (!inFASTA.eof()) { char c = inFASTA.get(); if (c == 10 || c == 13){ break; } }
+ in.close();
+ out.close();
}
- outGroups.close();
- inFASTA.close();
}
- if (m->control_pressed) {
- for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); }
- return 0;
+ if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
+
+ //output group counts
+ m->mothurOutEndLine();
+ //int total = 0;
+// for (int i = 0; i < barcodeNameVector.size(); i++) {
+// if ((barcodeNameVector[i] != "") && (groupCounts[i] != 0)) { total += groupCounts[i]; m->mothurOut("Group " + barcodeNameVector[i] + " contains " + toString(groupCounts[i]) + " sequences."); m->mothurOutEndLine(); }
+// }
+// if (total != 0) { m->mothurOut("Total of all groups is " + toString(total)); m->mothurOutEndLine(); }
+
+ if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
+
+ //set fasta file as new current fastafile
+ string current = "";
+ itTypes = outputTypes.find("fasta");
+ if (itTypes != outputTypes.end()) {
+ if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setFastaFile(current); }
+ }
+
+ itTypes = outputTypes.find("qfile");
+ if (itTypes != outputTypes.end()) {
+ if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setQualFile(current); }
+ }
+
+ itTypes = outputTypes.find("group");
+ if (itTypes != outputTypes.end()) {
+ if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setGroupFile(current); }
}
m->mothurOutEndLine();
}
/**************************************************************************************/
-int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string trimFile, string scrapFile, string groupFile, vector<string> fastaNames, linePair* line, linePair* qline) {
+
+int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string trimFileName, string scrapFileName, string trimQFileName, string scrapQFileName, string groupFileName, vector<vector<string> > fastaFileNames, vector<vector<string> > qualFileNames, linePair* line, linePair* qline) {
+
try {
- ofstream outFASTA;
- int able = openOutputFile(trimFile, outFASTA);
+ ofstream trimFASTAFile;
+ m->openOutputFile(trimFileName, trimFASTAFile);
- ofstream scrapFASTA;
- openOutputFile(scrapFile, scrapFASTA);
+ ofstream scrapFASTAFile;
+ m->openOutputFile(scrapFileName, scrapFASTAFile);
- ofstream outGroups;
- vector<ofstream*> fastaFileNames;
- if (oligoFile != "") {
- openOutputFile(groupFile, outGroups);
- for (int i = 0; i < fastaNames.size(); i++) {
- fastaFileNames.push_back(new ofstream((fastaNames[i] + toString(getpid()) + ".temp").c_str(), ios::ate));
+ ofstream trimQualFile;
+ ofstream scrapQualFile;
+ if(qFileName != ""){
+ m->openOutputFile(trimQFileName, trimQualFile);
+ m->openOutputFile(scrapQFileName, scrapQualFile);
+ }
+
+ ofstream outGroupsFile;
+ if (oligoFile != ""){ m->openOutputFile(groupFileName, outGroupsFile); }
+ if(allFiles){
+ for (int i = 0; i < fastaFileNames.size(); i++) { //clears old file
+ for (int j = 0; j < fastaFileNames[i].size(); j++) { //clears old file
+ if (fastaFileNames[i][j] != "") {
+ ofstream temp;
+ m->openOutputFile(fastaFileNames[i][j], temp); temp.close();
+ if(qFileName != ""){
+ m->openOutputFile(qualFileNames[i][j], temp); temp.close();
+ }
+ }
+ }
}
}
ifstream inFASTA;
- openInputFile(filename, inFASTA);
+ m->openInputFile(filename, inFASTA);
+ inFASTA.seekg(line->start);
ifstream qFile;
- if(qFileName != "") { openInputFile(qFileName, qFile); }
-
- qFile.seekg(qline->start);
- inFASTA.seekg(line->start);
+ if(qFileName != "") {
+ m->openInputFile(qFileName, qFile);
+ qFile.seekg(qline->start);
+ }
- for(int i=0;i<line->num;i++){
+ int count = 0;
+ bool moreSeqs = 1;
+
+ while (moreSeqs) {
if (m->control_pressed) {
- inFASTA.close();
- outFASTA.close();
- scrapFASTA.close();
- if (oligoFile != "") { outGroups.close(); }
- if(qFileName != "") { qFile.close(); }
- for(int i=0;i<fastaFileNames.size();i++){
- fastaFileNames[i]->close();
- delete fastaFileNames[i];
- }
+ inFASTA.close(); trimFASTAFile.close(); scrapFASTAFile.close();
+ if (oligoFile != "") { outGroupsFile.close(); }
+
+ if(qFileName != ""){
+ qFile.close();
+ }
for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); }
+
return 0;
}
- bool success = 1;
-
- Sequence currSeq(inFASTA);
+ int success = 1;
+ string trashCode = "";
+ int currentSeqsDiffs = 0;
+
+ Sequence currSeq(inFASTA); m->gobble(inFASTA);
+
+ QualityScores currQual;
+ if(qFileName != ""){
+ currQual = QualityScores(qFile); m->gobble(qFile);
+ }
string origSeq = currSeq.getUnaligned();
if (origSeq != "") {
- int group;
- string trashCode = "";
- int currentSeqsDiffs = 0;
- currentSeqsTdiffs = 0;
- if(qFileName != ""){
- if(qThreshold != 0) { success = stripQualThreshold(currSeq, qFile); }
- else if(qAverage != 0) { success = cullQualAverage(currSeq, qFile); }
- if ((!qtrim) && (origSeq.length() != currSeq.getUnaligned().length())) {
- success = 0; //if you don't want to trim and the sequence does not meet quality requirements, move to scrap
- }
- if(!success) { trashCode += 'q'; }
- }
-
+ int barcodeIndex = 0;
+ int primerIndex = 0;
+
if(barcodes.size() != 0){
- success = stripBarcode(currSeq, group);
-// cout << "here: " << success << endl;
- if(!success){ trashCode += 'b'; }
- else{ currentSeqsDiffs += currentSeqsTdiffs; }
+ success = stripBarcode(currSeq, currQual, barcodeIndex);
+ if(success > bdiffs) { trashCode += 'b'; }
+ else{ currentSeqsDiffs += success; }
}
-
+
if(numFPrimers != 0){
- success = stripForward(currSeq);
- if(!success){ trashCode += 'f'; }
- else{ currentSeqsDiffs += currentSeqsTdiffs; }
+ success = stripForward(currSeq, currQual, primerIndex);
+ if(success > pdiffs) { trashCode += 'f'; }
+ else{ currentSeqsDiffs += success; }
}
- if (currentSeqsDiffs > tdiffs) { trashCode += 't'; }
-
+ if (currentSeqsDiffs > tdiffs) { trashCode += 't'; }
+
if(numRPrimers != 0){
- success = stripReverse(currSeq);
- if(!success){ trashCode += 'r'; }
+ success = stripReverse(currSeq, currQual);
+ if(!success) { trashCode += 'r'; }
}
+
+ if(keepFirst != 0){
+ success = keepFirstTrim(currSeq, currQual);
+ }
+
+ if(removeLast != 0){
+ success = removeLastTrim(currSeq, currQual);
+ if(!success) { trashCode += 'l'; }
+ }
+
+
+ if(qFileName != ""){
+ int origLength = currSeq.getNumBases();
+
+ if(qThreshold != 0) { success = currQual.stripQualThreshold(currSeq, qThreshold); }
+ else if(qAverage != 0) { success = currQual.cullQualAverage(currSeq, qAverage); }
+ else if(qRollAverage != 0) { success = currQual.stripQualRollingAverage(currSeq, qRollAverage); }
+ else if(qWindowAverage != 0){ success = currQual.stripQualWindowAverage(currSeq, qWindowStep, qWindowSize, qWindowAverage); }
+ else { success = 1; }
+
+ //you don't want to trim, if it fails above then scrap it
+ if ((!qtrim) && (origLength != currSeq.getNumBases())) { success = 0; }
+
+ if(!success) { trashCode += 'q'; }
+ }
if(minLength > 0 || maxLength > 0){
success = cullLength(currSeq);
- if(!success){ trashCode += 'l'; }
+ if(!success) { trashCode += 'l'; }
}
if(maxHomoP > 0){
success = cullHomoP(currSeq);
- if(!success){ trashCode += 'h'; }
+ if(!success) { trashCode += 'h'; }
}
if(maxAmbig != -1){
success = cullAmbigs(currSeq);
- if(!success){ trashCode += 'n'; }
+ if(!success) { trashCode += 'n'; }
}
- if(flip){ currSeq.reverseComplement(); } // should go last
+ if(flip){ // should go last
+ currSeq.reverseComplement();
+ if(qFileName != ""){
+ currQual.flipQScores();
+ }
+ }
if(trashCode.length() == 0){
- currSeq.setAligned(currSeq.getUnaligned()); //this is because of a modification we made to the sequence class to fix a bug. all seqs have an aligned version, which is the version that gets printed.
- currSeq.printSequence(outFASTA);
+ currSeq.setAligned(currSeq.getUnaligned());
+ currSeq.printSequence(trimFASTAFile);
+
+ if(qFileName != ""){
+ currQual.printQScores(trimQualFile);
+ }
+
if(barcodes.size() != 0){
- outGroups << currSeq.getName() << '\t' << groupVector[group] << endl;
+ outGroupsFile << currSeq.getName() << '\t' << barcodeNameVector[barcodeIndex] << endl;
+ groupCounts[barcodeIndex]++;
+ }
+
+
+ if(allFiles){
+ ofstream output;
+ m->openOutputFileAppend(fastaFileNames[barcodeIndex][primerIndex], output);
+ currSeq.printSequence(output);
+ output.close();
- if(allFiles){
- currSeq.printSequence(*fastaFileNames[group]);
+ if(qFileName != ""){
+ m->openOutputFileAppend(qualFileNames[barcodeIndex][primerIndex], output);
+ currQual.printQScores(output);
+ output.close();
}
}
}
else{
currSeq.setName(currSeq.getName() + '|' + trashCode);
currSeq.setUnaligned(origSeq);
- currSeq.printSequence(scrapFASTA);
+ currSeq.setAligned(origSeq);
+ currSeq.printSequence(scrapFASTAFile);
+ if(qFileName != ""){
+ currQual.printQScores(scrapQualFile);
+ }
}
+ count++;
}
- gobble(inFASTA);
+
+ #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+ unsigned long int pos = inFASTA.tellg();
+ if ((pos == -1) || (pos >= line->end)) { break; }
+ #else
+ if (inFASTA.eof()) { break; }
+ #endif
+
+ //report progress
+ if((count) % 1000 == 0){ m->mothurOut(toString(count)); m->mothurOutEndLine(); }
+
}
+ //report progress
+ if((count) % 1000 != 0){ m->mothurOut(toString(count)); m->mothurOutEndLine(); }
+
inFASTA.close();
- outFASTA.close();
- scrapFASTA.close();
- if (oligoFile != "") { outGroups.close(); }
- if(qFileName != "") { qFile.close(); }
-
- for(int i=0;i<fastaFileNames.size();i++){
- fastaFileNames[i]->close();
- delete fastaFileNames[i];
- }
+ trimFASTAFile.close();
+ scrapFASTAFile.close();
+ if (oligoFile != "") { outGroupsFile.close(); }
+ if(qFileName != "") { qFile.close(); scrapQualFile.close(); trimQualFile.close(); }
- return 0;
+ return count;
}
catch(exception& e) {
m->errorOut(e, "TrimSeqsCommand", "driverCreateTrim");
exit(1);
}
}
+
/**************************************************************************************************/
-int TrimSeqsCommand::createProcessesCreateTrim(string filename, string qFileName, string trimFile, string scrapFile, string groupFile, vector<string> fastaNames) {
+
+int TrimSeqsCommand::createProcessesCreateTrim(string filename, string qFileName, string trimFASTAFileName, string scrapFASTAFileName, string trimQualFileName, string scrapQualFileName, string groupFile, vector<vector<string> > fastaFileNames, vector<vector<string> > qualFileNames) {
try {
#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
- int process = 0;
+ int process = 1;
int exitCommand = 1;
processIDS.clear();
processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
process++;
}else if (pid == 0){
- driverCreateTrim(filename, qFileName, (trimFile + toString(getpid()) + ".temp"), (scrapFile + toString(getpid()) + ".temp"), (groupFile + toString(getpid()) + ".temp"), fastaNames, lines[process], qLines[process]);
+
+ vector<vector<string> > tempFASTAFileNames = fastaFileNames;
+ vector<vector<string> > tempPrimerQualFileNames = qualFileNames;
+
+ if(allFiles){
+ ofstream temp;
+
+ for(int i=0;i<tempFASTAFileNames.size();i++){
+ for(int j=0;j<tempFASTAFileNames[i].size();j++){
+ if (tempFASTAFileNames[i][j] != "") {
+ tempFASTAFileNames[i][j] += toString(getpid()) + ".temp";
+ m->openOutputFile(tempFASTAFileNames[i][j], temp); temp.close();
+
+ if(qFileName != ""){
+ tempPrimerQualFileNames[i][j] += toString(getpid()) + ".temp";
+ m->openOutputFile(tempPrimerQualFileNames[i][j], temp); temp.close();
+ }
+ }
+ }
+ }
+ }
+
+ driverCreateTrim(filename,
+ qFileName,
+ (trimFASTAFileName + toString(getpid()) + ".temp"),
+ (scrapFASTAFileName + toString(getpid()) + ".temp"),
+ (trimQualFileName + toString(getpid()) + ".temp"),
+ (scrapQualFileName + toString(getpid()) + ".temp"),
+ (groupFile + toString(getpid()) + ".temp"),
+ tempFASTAFileNames,
+ tempPrimerQualFileNames,
+ lines[process],
+ qLines[process]);
+
+ //pass groupCounts to parent
+ ofstream out;
+ string tempFile = filename + toString(getpid()) + ".num.temp";
+ m->openOutputFile(tempFile, out);
+ for(int i = 0; i < groupCounts.size(); i++) {
+ out << groupCounts[i] << endl;
+ }
+ out.close();
+
+ exit(0);
+ }else {
+ m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
+ for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
exit(0);
- }else { m->mothurOut("unable to spawn the necessary processes."); m->mothurOutEndLine(); exit(0); }
+ }
}
+ //parent do my part
+ ofstream temp;
+ m->openOutputFile(trimFASTAFileName, temp); temp.close();
+ m->openOutputFile(scrapFASTAFileName, temp); temp.close();
+ m->openOutputFile(trimQualFileName, temp); temp.close();
+ m->openOutputFile(scrapQualFileName, temp); temp.close();
+
+ driverCreateTrim(filename, qFileName, trimFASTAFileName, scrapFASTAFileName, trimQualFileName, scrapQualFileName, groupFile, fastaFileNames, qualFileNames, lines[0], qLines[0]);
+
//force parent to wait until all the processes are done
- for (int i=0;i<processors;i++) {
+ for (int i=0;i<processIDS.size();i++) {
int temp = processIDS[i];
wait(&temp);
}
+ //append files
+ for(int i=0;i<processIDS.size();i++){
+
+ m->mothurOut("Appending files from process " + toString(processIDS[i])); m->mothurOutEndLine();
+
+ m->appendFiles((trimFASTAFileName + toString(processIDS[i]) + ".temp"), trimFASTAFileName);
+ remove((trimFASTAFileName + toString(processIDS[i]) + ".temp").c_str());
+ m->appendFiles((scrapFASTAFileName + toString(processIDS[i]) + ".temp"), scrapFASTAFileName);
+ remove((scrapFASTAFileName + toString(processIDS[i]) + ".temp").c_str());
+
+ if(qFileName != ""){
+ m->appendFiles((trimQualFileName + toString(processIDS[i]) + ".temp"), trimQualFileName);
+ remove((trimQualFileName + toString(processIDS[i]) + ".temp").c_str());
+ m->appendFiles((scrapQualFileName + toString(processIDS[i]) + ".temp"), scrapQualFileName);
+ remove((scrapQualFileName + toString(processIDS[i]) + ".temp").c_str());
+ }
+
+ m->appendFiles((groupFile + toString(processIDS[i]) + ".temp"), groupFile);
+ remove((groupFile + toString(processIDS[i]) + ".temp").c_str());
+
+
+ if(allFiles){
+ for(int j=0;j<fastaFileNames.size();j++){
+ for(int k=0;k<fastaFileNames[j].size();k++){
+ if (fastaFileNames[j][k] != "") {
+ m->appendFiles((fastaFileNames[j][k] + toString(processIDS[i]) + ".temp"), fastaFileNames[j][k]);
+ remove((fastaFileNames[j][k] + toString(processIDS[i]) + ".temp").c_str());
+
+ if(qFileName != ""){
+ m->appendFiles((qualFileNames[j][k] + toString(processIDS[i]) + ".temp"), qualFileNames[j][k]);
+ remove((qualFileNames[j][k] + toString(processIDS[i]) + ".temp").c_str());
+ }
+ }
+ }
+ }
+ }
+
+ ifstream in;
+ string tempFile = filename + toString(processIDS[i]) + ".num.temp";
+ m->openInputFile(tempFile, in);
+ int count = 0;
+ int tempNum;
+ while (!in.eof()) {
+ in >> tempNum; m->gobble(in);
+ groupCounts[count] += tempNum;
+ count++;
+ }
+ in.close(); remove(tempFile.c_str());
+
+ }
+
return exitCommand;
#endif
}
exit(1);
}
}
+
/**************************************************************************************************/
-int TrimSeqsCommand::setLines(string filename, vector<linePair*>& lines) {
+int TrimSeqsCommand::setLines(string filename, string qfilename, vector<unsigned long int>& fastaFilePos, vector<unsigned long int>& qfileFilePos) {
try {
- lines.clear();
+ //set file positions for fasta file
+ fastaFilePos = m->divideFile(filename, processors);
- vector<long int> positions;
+ if (qfilename == "") { return processors; }
+
+ //get name of first sequence in each chunk
+ map<string, int> firstSeqNames;
+ for (int i = 0; i < (fastaFilePos.size()-1); i++) {
+ ifstream in;
+ m->openInputFile(filename, in);
+ in.seekg(fastaFilePos[i]);
+
+ Sequence temp(in);
+ firstSeqNames[temp.getName()] = i;
+
+ in.close();
+ }
+
+ //seach for filePos of each first name in the qfile and save in qfileFilePos
+ ifstream inQual;
+ m->openInputFile(qfilename, inQual);
- ifstream inFASTA;
- openInputFile(filename, inFASTA);
-
string input;
- while(!inFASTA.eof()){
- input = getline(inFASTA);
+ while(!inQual.eof()){
+ input = m->getline(inQual);
if (input.length() != 0) {
- if(input[0] == '>'){ long int pos = inFASTA.tellg(); positions.push_back(pos - input.length() - 1); }
+ if(input[0] == '>'){ //this is a sequence name line
+ istringstream nameStream(input);
+
+ string sname = ""; nameStream >> sname;
+ sname = sname.substr(1);
+
+ map<string, int>::iterator it = firstSeqNames.find(sname);
+
+ if(it != firstSeqNames.end()) { //this is the start of a new chunk
+ unsigned long int pos = inQual.tellg();
+ qfileFilePos.push_back(pos - input.length() - 1);
+ firstSeqNames.erase(it);
+ }
+ }
}
+
+ if (firstSeqNames.size() == 0) { break; }
}
- inFASTA.close();
+ inQual.close();
- int numFastaSeqs = positions.size();
-
+
+ if (firstSeqNames.size() != 0) {
+ for (map<string, int>::iterator it = firstSeqNames.begin(); it != firstSeqNames.end(); it++) {
+ m->mothurOut(it->first + " is in your fasta file and not in your quality file, not using quality file."); m->mothurOutEndLine();
+ }
+ qFileName = "";
+ return processors;
+ }
+
+ //get last file position of qfile
FILE * pFile;
- long size;
+ unsigned long int size;
//get num bytes in file
- pFile = fopen (filename.c_str(),"rb");
+ pFile = fopen (qfilename.c_str(),"rb");
if (pFile==NULL) perror ("Error opening file");
else{
fseek (pFile, 0, SEEK_END);
fclose (pFile);
}
- int numSeqsPerProcessor = numFastaSeqs / processors;
-
- for (int i = 0; i < processors; i++) {
-
- long int startPos = positions[ i * numSeqsPerProcessor ];
- if(i == processors - 1){
- numSeqsPerProcessor = numFastaSeqs - i * numSeqsPerProcessor;
- }else{
- long int myEnd = positions[ (i+1) * numSeqsPerProcessor ];
- }
- lines.push_back(new linePair(startPos, numSeqsPerProcessor));
- }
+ qfileFilePos.push_back(size);
- return numFastaSeqs;
+ return processors;
}
catch(exception& e) {
m->errorOut(e, "TrimSeqsCommand", "setLines");
exit(1);
}
}
+
//***************************************************************************************************************
-void TrimSeqsCommand::getOligos(vector<string>& outFASTAVec){ //vector<ofstream*>& outFASTAVec
+void TrimSeqsCommand::getOligos(vector<vector<string> >& fastaFileNames, vector<vector<string> >& qualFileNames){
try {
ifstream inOligos;
- openInputFile(oligoFile, inOligos);
+ m->openInputFile(oligoFile, inOligos);
ofstream test;
string type, oligo, group;
- int index=0;
+
+ int indexPrimer = 0;
+ int indexBarcode = 0;
while(!inOligos.eof()){
- inOligos >> type;
-
+
+ inOligos >> type; m->gobble(inOligos);
+
if(type[0] == '#'){
while (!inOligos.eof()) { char c = inOligos.get(); if (c == 10 || c == 13){ break; } } // get rest of line if there's any crap there
}
else{
+ //make type case insensitive
+ for(int i=0;i<type.length();i++){ type[i] = toupper(type[i]); }
+
inOligos >> oligo;
for(int i=0;i<oligo.length();i++){
if(oligo[i] == 'U') { oligo[i] = 'T'; }
}
- if(type == "forward"){
- forPrimer.push_back(oligo);
+ if(type == "FORWARD"){
+ group = "";
+
+ // get rest of line in case there is a primer name
+ while (!inOligos.eof()) {
+ char c = inOligos.get();
+ if (c == 10 || c == 13){ break; }
+ else if (c == 32 || c == 9){;} //space or tab
+ else { group += c; }
+ }
+
+ //check for repeat barcodes
+ map<string, int>::iterator itPrime = primers.find(oligo);
+ if (itPrime != primers.end()) { m->mothurOut("primer " + oligo + " is in your oligos file already."); m->mothurOutEndLine(); }
+
+ primers[oligo]=indexPrimer; indexPrimer++;
+ primerNameVector.push_back(group);
}
- else if(type == "reverse"){
+ else if(type == "REVERSE"){
Sequence oligoRC("reverse", oligo);
oligoRC.reverseComplement();
revPrimer.push_back(oligoRC.getUnaligned());
}
- else if(type == "barcode"){
+ else if(type == "BARCODE"){
inOligos >> group;
- barcodes[oligo]=index++;
- groupVector.push_back(group);
- if(allFiles){
- //outFASTAVec.push_back(new ofstream((outputDir + getRootName(getSimpleName(fastaFile)) + group + ".fasta").c_str(), ios::ate));
- outputNames.push_back((outputDir + getRootName(getSimpleName(fastaFile)) + group + ".fasta"));
- outFASTAVec.push_back((outputDir + getRootName(getSimpleName(fastaFile)) + group + ".fasta"));
- }
+ //check for repeat barcodes
+ map<string, int>::iterator itBar = barcodes.find(oligo);
+ if (itBar != barcodes.end()) { m->mothurOut("barcode " + oligo + " is in your oligos file already."); m->mothurOutEndLine(); }
+
+ barcodes[oligo]=indexBarcode; indexBarcode++;
+ barcodeNameVector.push_back(group);
}
+ else{ m->mothurOut(type + " is not recognized as a valid type. Choices are forward, reverse, and barcode. Ignoring " + oligo + "."); m->mothurOutEndLine(); }
}
+ m->gobble(inOligos);
+ }
+ inOligos.close();
+
+ if(barcodeNameVector.size() == 0 && primerNameVector[0] == ""){ allFiles = 0; }
+
+ //add in potential combos
+ if(barcodeNameVector.size() == 0){
+ barcodes[""] = 0;
+ barcodeNameVector.push_back("");
}
- inOligos.close();
+ if(primerNameVector.size() == 0){
+ primers[""] = 0;
+ primerNameVector.push_back("");
+ }
- numFPrimers = forPrimer.size();
- numRPrimers = revPrimer.size();
+ fastaFileNames.resize(barcodeNameVector.size());
+ for(int i=0;i<fastaFileNames.size();i++){
+ fastaFileNames[i].assign(primerNameVector.size(), "");
+ }
+ if(qFileName != ""){ qualFileNames = fastaFileNames; }
+ if(allFiles){
+ set<string> uniqueNames; //used to cleanup outputFileNames
+ for(map<string, int>::iterator itBar = barcodes.begin();itBar != barcodes.end();itBar++){
+ for(map<string, int>::iterator itPrimer = primers.begin();itPrimer != primers.end(); itPrimer++){
+
+ string primerName = primerNameVector[itPrimer->second];
+ string barcodeName = barcodeNameVector[itBar->second];
+
+ string comboGroupName = "";
+ string fastaFileName = "";
+ string qualFileName = "";
+
+ if(primerName == ""){
+ comboGroupName = barcodeNameVector[itBar->second];
+ }
+ else{
+ if(barcodeName == ""){
+ comboGroupName = primerNameVector[itPrimer->second];
+ }
+ else{
+ comboGroupName = barcodeNameVector[itBar->second] + "." + primerNameVector[itPrimer->second];
+ }
+ }
+
+ ofstream temp;
+ fastaFileName = outputDir + m->getRootName(m->getSimpleName(fastaFile)) + comboGroupName + ".fasta";
+ if (uniqueNames.count(fastaFileName) == 0) {
+ outputNames.push_back(fastaFileName);
+ outputTypes["fasta"].push_back(fastaFileName);
+ uniqueNames.insert(fastaFileName);
+ }
+
+ fastaFileNames[itBar->second][itPrimer->second] = fastaFileName;
+ m->openOutputFile(fastaFileName, temp); temp.close();
+
+ if(qFileName != ""){
+ qualFileName = outputDir + m->getRootName(m->getSimpleName(qFileName)) + comboGroupName + ".qual";
+ if (uniqueNames.count(fastaFileName) == 0) {
+ outputNames.push_back(qualFileName);
+ outputTypes["qfile"].push_back(qualFileName);
+ }
+
+ qualFileNames[itBar->second][itPrimer->second] = qualFileName;
+ m->openOutputFile(qualFileName, temp); temp.close();
+ }
+ }
+ }
+ }
+ numFPrimers = primers.size();
+ numRPrimers = revPrimer.size();
+ groupCounts.resize(barcodeNameVector.size(), 0);
+
}
catch(exception& e) {
m->errorOut(e, "TrimSeqsCommand", "getOligos");
exit(1);
}
}
+
//***************************************************************************************************************
-bool TrimSeqsCommand::stripBarcode(Sequence& seq, int& group){
+int TrimSeqsCommand::stripBarcode(Sequence& seq, QualityScores& qual, int& group){
try {
+
string rawSequence = seq.getUnaligned();
- bool success = 0; //guilty until proven innocent
+ int success = bdiffs + 1; //guilty until proven innocent
//can you find the barcode
for(map<string,int>::iterator it=barcodes.begin();it!=barcodes.end();it++){
string oligo = it->first;
if(rawSequence.length() < oligo.length()){ //let's just assume that the barcodes are the same length
- success = 0;
- break;
+ success = bdiffs + 10; //if the sequence is shorter than the barcode then bail out
+ break;
}
if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
group = it->second;
seq.setUnaligned(rawSequence.substr(oligo.length()));
- success = 1;
+
+ if(qual.getName() != ""){
+ qual.trimQScores(oligo.length(), -1);
+ }
+
+ success = 0;
break;
}
}
//if you found the barcode or if you don't want to allow for diffs
-// cout << success;
- if ((bdiffs == 0) || (success == 1)) { return success; }
+ if ((bdiffs == 0) || (success == 0)) { return success; }
else { //try aligning and see if you can find it
-// cout << endl;
int maxLength = 0;
maxLength = it->first.length();
}
}
- alignment = new NeedlemanOverlap(-2.0, 1.0, -1.0, (maxLength+bdiffs+1));
+ alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+bdiffs+1));
}else{ alignment = NULL; }
// int length = oligo.length();
if(rawSequence.length() < maxLength){ //let's just assume that the barcodes are the same length
- success = 0;
+ success = bdiffs + 10;
break;
}
}
oligo = oligo.substr(0,alnLength);
temp = temp.substr(0,alnLength);
-// cout << "barcode = " << oligo << " raw = " << rawSequence.substr(0,alnLength) << " raw aligned = " << temp << endl;
- cout << seq.getName() << endl;
- cout << temp << endl;
- cout << oligo << endl;
- cout << alnLength << endl;
- cout << endl;
- int newStart=0;
- int numDiff = countDiffs(oligo, temp);//, alnLength, newStart, bdiffs);
+ int numDiff = countDiffs(oligo, temp);
+
if(numDiff < minDiff){
minDiff = numDiff;
minCount = 1;
}
}
- if(minDiff > bdiffs){ success = 0; }
- else if(minCount > 1) { success = 0; }
- else{
+
+ if(minDiff > bdiffs) { success = minDiff; } //no good matches
+ else if(minCount > 1) { success = bdiffs + 100; } //can't tell the difference between multiple barcodes
+ else{ //use the best match
group = minGroup;
- seq.setUnaligned("*" + rawSequence.substr(minPos));
- success = 1;
+ seq.setUnaligned(rawSequence.substr(minPos));
+
+ if(qual.getName() != ""){
+ qual.trimQScores(minPos, -1);
+ }
+ success = minDiff;
}
if (alignment != NULL) { delete alignment; }
+
}
+
return success;
}
//***************************************************************************************************************
-bool TrimSeqsCommand::stripForward(Sequence& seq){
+int TrimSeqsCommand::stripForward(Sequence& seq, QualityScores& qual, int& group){
try {
string rawSequence = seq.getUnaligned();
- bool success = 0; //guilty until proven innocent
+ int success = pdiffs + 1; //guilty until proven innocent
- for(int i=0;i<numFPrimers;i++){
- string oligo = forPrimer[i];
-
- if(rawSequence.length() < oligo.length()){
- success = 0;
- break;
+ //can you find the primer
+ for(map<string,int>::iterator it=primers.begin();it!=primers.end();it++){
+ string oligo = it->first;
+ if(rawSequence.length() < oligo.length()){ //let's just assume that the primers are the same length
+ success = pdiffs + 10; //if the sequence is shorter than the barcode then bail out
+ break;
}
if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
+ group = it->second;
seq.setUnaligned(rawSequence.substr(oligo.length()));
- success = 1;
+ if(qual.getName() != ""){
+ qual.trimQScores(oligo.length(), -1);
+ }
+ success = 0;
break;
}
}
-
- //if you found the primer or if you don't want to allow for diffs
- if ((pdiffs == 0) || (success == 1)) { return success; }
+
+ //if you found the barcode or if you don't want to allow for diffs
+ if ((pdiffs == 0) || (success == 0)) { return success; }
else { //try aligning and see if you can find it
-
+
+ int maxLength = 0;
+
Alignment* alignment;
- if (numFPrimers > 0) { alignment = new NeedlemanOverlap(-2.0, 1.0, -1.0, (forPrimer[0].length()+pdiffs+1)); }
- else{ alignment = NULL; }
- //can you find the primer
- for(int i=0;i<numFPrimers;i++){
- string oligo = forPrimer[i];
- int length = oligo.length();
+ if (primers.size() > 0) {
+ map<string,int>::iterator it=primers.begin();
+
+ for(it;it!=primers.end();it++){
+ if(it->first.length() > maxLength){
+ maxLength = it->first.length();
+ }
+ }
+ alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+pdiffs+1));
+
+ }else{ alignment = NULL; }
- if(rawSequence.length() < oligo.length()){
- success = 0;
+ //can you find the barcode
+ int minDiff = 1e6;
+ int minCount = 1;
+ int minGroup = -1;
+ int minPos = 0;
+
+ for(map<string,int>::iterator it=primers.begin();it!=primers.end();it++){
+ string oligo = it->first;
+// int length = oligo.length();
+
+ if(rawSequence.length() < maxLength){
+ success = pdiffs + 100;
break;
}
-
- //resize if neccessary
- if ((length+pdiffs+1) > alignment->getnRows()) { alignment->resize(length+pdiffs+1); }
- //use needleman to align first primer.length()+numdiffs of sequence to each primer
- alignment->align(oligo, rawSequence.substr(0,length+pdiffs));
+ //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
+ alignment->align(oligo, rawSequence.substr(0,oligo.length()+pdiffs));
oligo = alignment->getSeqAAln();
string temp = alignment->getSeqBAln();
-
- int newStart = 0;
-// if(compareDNASeq(oligo, temp, length, newStart, pdiffs)){
-// seq.setUnaligned(rawSequence.substr(newStart));
-// success = 1;
-// break;
-// }
+
+ int alnLength = oligo.length();
+
+ for(int i=oligo.length()-1;i>=0;i--){
+ if(oligo[i] != '-'){ alnLength = i+1; break; }
+ }
+ oligo = oligo.substr(0,alnLength);
+ temp = temp.substr(0,alnLength);
+
+ int numDiff = countDiffs(oligo, temp);
+
+ if(numDiff < minDiff){
+ minDiff = numDiff;
+ minCount = 1;
+ minGroup = it->second;
+ minPos = 0;
+ for(int i=0;i<alnLength;i++){
+ if(temp[i] != '-'){
+ minPos++;
+ }
+ }
+ }
+ else if(numDiff == minDiff){
+ minCount++;
+ }
+
+ }
+
+ if(minDiff > pdiffs) { success = minDiff; } //no good matches
+ else if(minCount > 1) { success = pdiffs + 10; } //can't tell the difference between multiple primers
+ else{ //use the best match
+ group = minGroup;
+ seq.setUnaligned(rawSequence.substr(minPos));
+ if(qual.getName() != ""){
+ qual.trimQScores(minPos, -1);
+ }
+ success = minDiff;
}
if (alignment != NULL) { delete alignment; }
+
}
return success;
-
+
}
catch(exception& e) {
m->errorOut(e, "TrimSeqsCommand", "stripForward");
//***************************************************************************************************************
-bool TrimSeqsCommand::stripReverse(Sequence& seq){
+bool TrimSeqsCommand::stripReverse(Sequence& seq, QualityScores& qual){
try {
string rawSequence = seq.getUnaligned();
bool success = 0; //guilty until proven innocent
if(compareDNASeq(oligo, rawSequence.substr(rawSequence.length()-oligo.length(),oligo.length()))){
seq.setUnaligned(rawSequence.substr(0,rawSequence.length()-oligo.length()));
+ if(qual.getName() != ""){
+ qual.trimQScores(-1, rawSequence.length()-oligo.length());
+ }
success = 1;
break;
}
//***************************************************************************************************************
+bool TrimSeqsCommand::keepFirstTrim(Sequence& sequence, QualityScores& qscores){
+ try {
+ bool success = 1;
+ if(qscores.getName() != ""){
+ qscores.trimQScores(-1, keepFirst);
+ }
+ sequence.trim(keepFirst);
+ return success;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "keepFirstTrim", "countDiffs");
+ exit(1);
+ }
+
+}
+
+//***************************************************************************************************************
+
+bool TrimSeqsCommand::removeLastTrim(Sequence& sequence, QualityScores& qscores){
+ try {
+ bool success = 0;
+
+ int length = sequence.getNumBases() - removeLast;
+
+ if(length > 0){
+ if(qscores.getName() != ""){
+ qscores.trimQScores(-1, length);
+ }
+ sequence.trim(length);
+ success = 1;
+ }
+ else{
+ success = 0;
+ }
+
+ return success;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "removeLastTrim", "countDiffs");
+ exit(1);
+ }
+
+}
+
+//***************************************************************************************************************
+
bool TrimSeqsCommand::cullLength(Sequence& seq){
try {
}
}
+
//***************************************************************************************************************
-int TrimSeqsCommand::countDiffs(string oligo, string seq){//, int numBases, int& end, int diffs){
+int TrimSeqsCommand::countDiffs(string oligo, string seq){
try {
-// bool success = 1;
+
int length = oligo.length();
-// end = numBases;
-// int countBases = 0;
int countDiffs = 0;
-
-// if (length != 0) {
-// if ((oligo[0] == '-') || (oligo[0] == '.')) { return 1e6; } //no gaps allowed at beginning
-// }
-
for(int i=0;i<length;i++){
-
-// if ((oligo[i] != '-') && (oligo[i] != '.')) { countBases++; }
-
+
if(oligo[i] != seq[i]){
if(oligo[i] == 'A' || oligo[i] == 'T' || oligo[i] == 'G' || oligo[i] == 'C' || oligo[i] == '-' || oligo[i] == '.') { countDiffs++; }
else if((oligo[i] == 'N' || oligo[i] == 'I') && (seq[i] == 'N')) { countDiffs++; }
else if(oligo[i] == 'B' && (seq[i] != 'C' && seq[i] != 'T' && seq[i] != 'G')) { countDiffs++; }
else if(oligo[i] == 'D' && (seq[i] != 'A' && seq[i] != 'T' && seq[i] != 'G')) { countDiffs++; }
else if(oligo[i] == 'H' && (seq[i] != 'A' && seq[i] != 'T' && seq[i] != 'C')) { countDiffs++; }
- else if(oligo[i] == 'V' && (seq[i] != 'A' && seq[i] != 'C' && seq[i] != 'G')) { countDiffs++; }
-
-// if(countDiffs > diffs) { break; }
+ else if(oligo[i] == 'V' && (seq[i] != 'A' && seq[i] != 'C' && seq[i] != 'G')) { countDiffs++; }
}
-// else{
-// success = 1;
-// }
-// if (countBases >= numBases) { end = countBases; break; } //stop checking after end of barcode or primer
}
- //if it's a success we want to check for total diffs in driver, so save it.
-// if (success == 1) { currentSeqsTdiffs = countDiffs; }
return countDiffs;
}
}
}
-//***************************************************************************************************************
-
-bool TrimSeqsCommand::stripQualThreshold(Sequence& seq, ifstream& qFile){
- try {
- string rawSequence = seq.getUnaligned();
- int seqLength; // = rawSequence.length();
- string name, temp, temp2;
-
- qFile >> name >> temp;
-
- splitAtEquals(temp2, temp); //separates length=242, temp=length, temp2=242
- convert(temp, seqLength); //converts string to int
-
- if (name.length() != 0) { if(name.substr(1) != seq.getName()) { m->mothurOut("sequence name mismatch btwn fasta and qual file"); m->mothurOutEndLine(); } }
- while (!qFile.eof()) { char c = qFile.get(); if (c == 10 || c == 13){ break; } }
-
- int score;
- int end = seqLength;
-
- for(int i=0;i<seqLength;i++){
- qFile >> score;
-
- if(score <= qThreshold){
- end = i;
- break;
- }
- }
- for(int i=end+1;i<seqLength;i++){
- qFile >> score;
- }
-
- seq.setUnaligned(rawSequence.substr(0,end));
-
- return 1;
- }
- catch(exception& e) {
- m->errorOut(e, "TrimSeqsCommand", "stripQualThreshold");
- exit(1);
- }
-}
-
-//***************************************************************************************************************
-
-bool TrimSeqsCommand::cullQualAverage(Sequence& seq, ifstream& qFile){
- try {
- string rawSequence = seq.getUnaligned();
- int seqLength = seq.getNumBases();
- bool success = 0; //guilty until proven innocent
- string name;
-
- qFile >> name;
- if (name[0] == '>') { if(name.substr(1) != seq.getName()) { m->mothurOut("sequence name mismatch btwn fasta: " + seq.getName() + " and qual file: " + name); m->mothurOutEndLine(); } }
-
- while (!qFile.eof()) { char c = qFile.get(); if (c == 10 || c == 13){ break; } }
-
- float score;
- float average = 0;
-
- for(int i=0;i<seqLength;i++){
- qFile >> score;
- average += score;
- }
- average /= seqLength;
-
- if(average >= qAverage) { success = 1; }
- else { success = 0; }
-
- return success;
- }
- catch(exception& e) {
- m->errorOut(e, "TrimSeqsCommand", "cullQualAverage");
- exit(1);
- }
-}
//***************************************************************************************************************
projects / mothur.git / blobdiff