//***************************************************************************************************************
-Bellerophon::Bellerophon(string name) {
+Bellerophon::Bellerophon(string name, string o) {
try {
fastafile = name;
+ outputDir = o;
}
catch(exception& e) {
errorOut(e, "Bellerophon", "Bellerophon");
out << pref[i].name << '\t' << setprecision(3) << pref[i].score[0] << '\t' << pref[i].leftParent[0] << '\t' << pref[i].rightParent[0] << endl;
//calc # of seqs with preference above 1.0
- if (pref[i].score[0] > 1.1) {
+ if (pref[i].score[0] > 1.0) {
above1++;
mothurOut(pref[i].name + " is a suspected chimera at breakpoint " + toString(pref[i].midpoint)); mothurOutEndLine();
mothurOut("It's score is " + toString(pref[i].score[0]) + " with suspected left parent " + pref[i].leftParent[0] + " and right parent " + pref[i].rightParent[0]); mothurOutEndLine();
//output results to screen
mothurOutEndLine();
- mothurOut("Sequence with preference score above 1.1: " + toString(above1)); mothurOutEndLine();
+ mothurOut("Sequence with preference score above 1.0: " + toString(above1)); mothurOutEndLine();
int spot;
spot = pref.size()-1;
mothurOut("Minimum:\t" + toString(pref[spot].score[0])); mothurOutEndLine();
}
//***************************************************************************************************************
-void Bellerophon::getChimeras() {
+int Bellerophon::getChimeras() {
try {
//do soft filter
if (filter) {
string optionString = "fasta=" + fastafile + ", soft=50";
+ if (outputDir != "") { optionString += ", outputdir=" + outputDir; }
+
filterSeqs = new FilterSeqsCommand(optionString);
filterSeqs->execute();
delete filterSeqs;
//reset fastafile to filtered file
- fastafile = getRootName(fastafile) + "filter.fasta";
+ if (outputDir == "") { fastafile = getRootName(fastafile) + "filter.fasta"; }
+ else { fastafile = outputDir + getRootName(getSimpleName(fastafile)) + "filter.fasta"; }
+
}
distCalculator = new eachGapDist();
//read in sequences
seqs = readSeqs(fastafile);
+ if (unaligned) { mothurOut("Your sequences need to be aligned when you use the bellerophon method."); mothurOutEndLine(); return 1; }
+
int numSeqs = seqs.size();
if (numSeqs == 0) { mothurOut("Error in reading you sequences."); mothurOutEndLine(); exit(1); }
vector<Sequence> left; vector<Sequence> right;
for (int i = 0; i < seqs.size(); i++) {
-//cout << "whole = " << seqs[i].getAligned() << endl;
+//cout << "midpoint = " << midpoint << "\twindow = " << window << endl;
+//cout << "whole = " << seqs[i]->getAligned().length() << endl;
//save left side
string seqLeft = seqs[i]->getAligned().substr(midpoint-window, window);
Sequence tempLeft;
tempLeft.setName(seqs[i]->getName());
tempLeft.setAligned(seqLeft);
left.push_back(tempLeft);
-//cout << "left = " << tempLeft.getAligned() << endl;
+//cout << "left = " << tempLeft.getAligned().length() << endl;
//save right side
string seqRight = seqs[i]->getAligned().substr(midpoint, window);
Sequence tempRight;
tempRight.setName(seqs[i]->getName());
tempRight.setAligned(seqRight);
right.push_back(tempRight);
-//cout << "right = " << seqRight << endl;
+//cout << "right = " << seqRight.length() << endl;
}
//adjust midpoint by increment
delete distCalculator;
//rank preference score to eachother
- /*float dme = 0.0;
+ float dme = 0.0;
float expectedPercent = 1 / (float) (pref.size());
for (int i = 0; i < pref.size(); i++) { dme += pref[i].score[0]; }
//how much higher or lower is this than expected
pref[i].score[0] = pref[i].score[0] / expectedPercent;
- }*/
+ }
//sort Preferences highest to lowest
sort(pref.begin(), pref.end(), comparePref);
+
+ return 0;
}
catch(exception& e) {
//calculate the dme
int count0 = 0;
- for (int i = 0; i < pref.size(); i++) { dme += pref[i].score[1]; if (pref[i].score[1] == 0.0) { count0++; } }
+ for (int i = 0; i < pref.size(); i++) { dme += pref[i].score[1]; if (pref[i].score[1] == 0.0) { count0++; } }
- //float expectedPercent = 1 / (float) (pref.size() - count0);
+ float expectedPercent = 1 / (float) (pref.size() - count0);
//cout << endl << "dme = " << dme << endl;
//recalculate prefernences based on dme
for (int i = 0; i < pref.size(); i++) {
-//cout << "unadjusted col[i] " << i << " = " << pref[i].score[1] << endl;
-//cout << i << '\t' << (dme / (dme - 2 * pref[i].score[1])) << endl;
-
+//cout << "unadjusted pref " << i << " = " << pref[i].score[1] << endl;
// gives the actual percentage of the dme this seq adds
- //pref[i].score[1] = pref[i].score[1] / dme;
+ pref[i].score[1] = pref[i].score[1] / dme;
//how much higher or lower is this than expected
- //pref[i].score[1] = pref[i].score[1] / expectedPercent;
+ pref[i].score[1] = pref[i].score[1] / expectedPercent;
+
+ //pref[i].score[1] = dme / (dme - 2 * pref[i].score[1]);
- //not 2 * pref[i].score[1] because i only calulate the pref scores once.
- pref[i].score[1] = dme / (dme - pref[i].score[1]);
-//cout << i << '\t' << pref[i].score[1] << endl;
//so a non chimeric sequence would be around 1, and a chimeric would be signifigantly higher.
//cout << "adjusted pref " << i << " = " << pref[i].score[1] << endl;
}