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 fasta parameter is required.\n");
- m->mothurOut("The flip parameter .... 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 .... The default is "".\n");
m->mothurOut("The maxambig parameter .... The default is -1.\n");
m->mothurOut("The maxhomop parameter .... The default is 0.\n");
#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
if(processors == 1){
ifstream inFASTA;
+ int numSeqs;
openInputFile(fastaFile, inFASTA);
- int numSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
+ getNumSeqs(inFASTA, numSeqs);
inFASTA.close();
lines.push_back(new linePair(0, numSeqs));
if (m->control_pressed) { return 0; }
#else
ifstream inFASTA;
+ int numSeqs;
openInputFile(fastaFile, inFASTA);
- int numSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
+ getNumSeqs(inFASTA, numSeqs);
inFASTA.close();
lines.push_back(new linePair(0, numSeqs));
for(int i=0;i<fastaFileNames.size();i++){
ifstream inFASTA;
string seqName;
- openInputFile(getRootName(fastaFile) + groupVector[i] + ".fasta", inFASTA);
+ //openInputFile(getRootName(fastaFile) + groupVector[i] + ".fasta", inFASTA);
+ openInputFile(fastaFileNames[i], inFASTA);
ofstream outGroups;
- openOutputFile(outputDir + getRootName(getSimpleName(fastaFile)) + groupVector[i] + ".groups", outGroups);
- outputNames.push_back(outputDir + getRootName(getSimpleName(fastaFile)) + groupVector[i] + ".groups");
+ string outGroupFilename = outputDir + getRootName(getSimpleName(fastaFileNames[i])) + "groups";
+ openOutputFile(outGroupFilename, outGroups);
+ //openOutputFile(outputDir + getRootName(getSimpleName(fastaFile)) + groupVector[i] + ".groups", outGroups);
+ outputNames.push_back(outGroupFilename);
while(!inFASTA.eof()){
if(inFASTA.get() == '>'){
for(int i=0;i<line->num;i++){
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(); 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]; }
for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); }
return 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())) {
+
+ if (qtrim == 1 && (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'; }
}
}
else if(type == "barcode"){
inOligos >> group;
- barcodes[oligo]=index++;
+ barcodes[oligo]=index; 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"));
+ outputNames.push_back((outputDir + getRootName(getSimpleName(fastaFile)) + toString(index) + "." + group + ".fasta"));
+ outFASTAVec.push_back((outputDir + getRootName(getSimpleName(fastaFile)) + toString(index) + "." + group + ".fasta"));
}
}
}
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 newStart=0;
int numDiff = countDiffs(oligo, temp);
+
+// cout << oligo << '\t' << temp << '\t' << numDiff << endl;
+
if(numDiff < minDiff){
minDiff = numDiff;
minCount = 1;
maxLength = forPrimer[i].length();
}
}
- alignment = new NeedlemanOverlap(-2.0, 1.0, -1.0, (maxLength+pdiffs+1));
+ alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+pdiffs+1));
}else{ alignment = NULL; }
bool TrimSeqsCommand::stripQualThreshold(Sequence& seq, ifstream& qFile){
try {
+// string rawSequence = seq.getUnaligned();
+// int seqLength; // = rawSequence.length();
+// string name, temp, temp2;
+//
+// qFile >> name;
+//
+// //get rest of line
+// temp = "";
+// while (!qFile.eof()) {
+// char c = qFile.get();
+// if (c == 10 || c == 13){ break; }
+// else { temp += c; }
+// }
+//
+// int pos = temp.find("length");
+// if (pos == temp.npos) { m->mothurOut("Cannot find length in qfile for " + seq.getName()); m->mothurOutEndLine(); seqLength = 0; }
+// else {
+// string tempLength = temp.substr(pos);
+// istringstream iss (tempLength,istringstream::in);
+// iss >> 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(); } }
+
string rawSequence = seq.getUnaligned();
- int seqLength; // = rawSequence.length();
- string name, temp, temp2;
+ 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(); } }
- //get rest of line
- temp = "";
- while (!qFile.eof()) {
- char c = qFile.get();
- if (c == 10 || c == 13){ break; }
- else { temp += c; }
- }
-
- int pos = temp.find("length");
- if (pos == temp.npos) { m->mothurOut("Cannot find length in qfile for " + seq.getName()); m->mothurOutEndLine(); seqLength = 0; }
- else {
- string tempLength = temp.substr(pos);
- istringstream iss (tempLength,istringstream::in);
- iss >> 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){
+ if(score < qThreshold){
end = i;
break;
}
projects / mothur.git / blobdiff