if (m->control_pressed) { return 0; }
#else
ifstream inFASTA;
- openInputFile(fastafileNames[s], inFASTA);
- numSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
+ openInputFile(fastaFile, inFASTA);
+ int 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]);
+ driverCreateTrim(fastaFile, qFileName, trimSeqFile, scrapSeqFile, groupFile, fastaFileNames, lines[0], lines[0]);
if (m->control_pressed) { return 0; }
#endif
if (oligoFile != "") {
openOutputFile(groupFile, outGroups);
for (int i = 0; i < fastaNames.size(); i++) {
+ #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
fastaFileNames.push_back(new ofstream((fastaNames[i] + toString(getpid()) + ".temp").c_str(), ios::ate));
+ #else
+ fastaFileNames.push_back(new ofstream((fastaNames[i] + toString(i) + ".temp").c_str(), ios::ate));
+ #endif
}
}
return 0;
}
- bool success = 1;
+ int success = 1;
Sequence currSeq(inFASTA);
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'; }
}
if(barcodes.size() != 0){
success = stripBarcode(currSeq, group);
-// cout << "here: " << success << endl;
if(success > bdiffs){ trashCode += 'b'; }
else{ currentSeqsDiffs += success; }
}
if(flip){ currSeq.reverseComplement(); } // should go last
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.setAligned(currSeq.getUnaligned());
currSeq.printSequence(outFASTA);
if(barcodes.size() != 0){
outGroups << currSeq.getName() << '\t' << groupVector[group] << endl;
else{
currSeq.setName(currSeq.getName() + '|' + trashCode);
currSeq.setUnaligned(origSeq);
+ currSeq.setAligned(origSeq);
currSeq.printSequence(scrapFASTA);
}
}
}
//***************************************************************************************************************
-bool TrimSeqsCommand::stripBarcode(Sequence& seq, int& group){
+int TrimSeqsCommand::stripBarcode(Sequence& seq, int& group){
try {
+
string rawSequence = seq.getUnaligned();
- bool success = bdiffs + 1; //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 = bdiffs + 1;
- break;
+ success = bdiffs + 10; //if the sequence is shorter than the barcode then bail out
+ break;
}
if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
// int length = oligo.length();
if(rawSequence.length() < maxLength){ //let's just assume that the barcodes are the same length
- success = bdiffs + 1;
+ 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);
}
}
- if(minDiff > bdiffs){ success = bdiffs + 1; }
- else if(minCount > 1) { success = bdiffs + 1; }
- 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));
+ seq.setUnaligned(rawSequence.substr(minPos));
success = minDiff;
}
if (alignment != NULL) { delete alignment; }
}
+// cout << success << endl;
+
return success;
}
//***************************************************************************************************************
-bool TrimSeqsCommand::stripForward(Sequence& seq){
+int TrimSeqsCommand::stripForward(Sequence& seq){
try {
string rawSequence = seq.getUnaligned();
- bool success = pdiffs + 1; //guilty until proven innocent
+ int success = pdiffs + 1; //guilty until proven innocent
//can you find the primer
for(int i=0;i<numFPrimers;i++){
string oligo = forPrimer[i];
if(rawSequence.length() < maxLength){
- success = pdiffs + 1;
+ success = pdiffs + 100;
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);
if(numDiff < minDiff){
}
}
- if(minDiff > pdiffs){ success = pdiffs + 1; }
- else if(minCount > 1) { success = pdiffs + 1; }
+ if(minDiff > pdiffs) { success = minDiff; }
+ else if(minCount > 1) { success = pdiffs + 10; }
else{
- seq.setUnaligned("*" + rawSequence.substr(minPos));
+ seq.setUnaligned(rawSequence.substr(minPos));
success = minDiff;
}
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(); } }
- 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;
for(int i=0;i<seqLength;i++){
qFile >> score;
- if(score <= qThreshold){
+ if(score < qThreshold){
end = i;
break;
}
projects / mothur.git / blobdiff