//***************************************************************************************************************
-Bellerophon::Bellerophon(string name) {
+Bellerophon::Bellerophon(string name, string o) {
try {
fastafile = name;
+ outputDir = o;
}
catch(exception& e) {
- errorOut(e, "Bellerophon", "Bellerophon");
+ m->errorOut(e, "Bellerophon", "Bellerophon");
exit(1);
}
}
//***************************************************************************************************************
-void Bellerophon::print(ostream& out) {
+int Bellerophon::print(ostream& out, ostream& outAcc) {
try {
int above1 = 0;
out << "Name\tScore\tLeft\tRight\t" << endl;
//output prefenence structure to .chimeras file
for (int i = 0; i < pref.size(); i++) {
- out << pref[i].name << '\t' << pref[i].score[0] << '\t' << pref[i].leftParent[0] << '\t' << pref[i].rightParent[0] << endl;
+
+ if (m->control_pressed) { return 0; }
+
+ 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.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();
+ outAcc << pref[i].name << endl;
+ m->mothurOut(pref[i].name + " is a suspected chimera at breakpoint " + toString(pref[i].midpoint)); m->mothurOutEndLine();
+ m->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]); m->mothurOutEndLine();
}
}
//output results to screen
- mothurOutEndLine();
- mothurOut("Sequence with preference score above 1.0: " + toString(above1)); mothurOutEndLine();
+ m->mothurOutEndLine();
+ m->mothurOut("Sequence with preference score above 1.0: " + toString(above1)); m->mothurOutEndLine();
int spot;
spot = pref.size()-1;
- mothurOut("Minimum:\t" + toString(pref[spot].score[0])); mothurOutEndLine();
+ m->mothurOut("Minimum:\t" + toString(pref[spot].score[0])); m->mothurOutEndLine();
spot = pref.size() * 0.975;
- mothurOut("2.5%-tile:\t" + toString(pref[spot].score[0])); mothurOutEndLine();
+ m->mothurOut("2.5%-tile:\t" + toString(pref[spot].score[0])); m->mothurOutEndLine();
spot = pref.size() * 0.75;
- mothurOut("25%-tile:\t" + toString(pref[spot].score[0])); mothurOutEndLine();
+ m->mothurOut("25%-tile:\t" + toString(pref[spot].score[0])); m->mothurOutEndLine();
spot = pref.size() * 0.50;
- mothurOut("Median: \t" + toString(pref[spot].score[0])); mothurOutEndLine();
+ m->mothurOut("Median: \t" + toString(pref[spot].score[0])); m->mothurOutEndLine();
spot = pref.size() * 0.25;
- mothurOut("75%-tile:\t" + toString(pref[spot].score[0])); mothurOutEndLine();
+ m->mothurOut("75%-tile:\t" + toString(pref[spot].score[0])); m->mothurOutEndLine();
spot = pref.size() * 0.025;
- mothurOut("97.5%-tile:\t" + toString(pref[spot].score[0])); mothurOutEndLine();
+ m->mothurOut("97.5%-tile:\t" + toString(pref[spot].score[0])); m->mothurOutEndLine();
spot = 0;
- mothurOut("Maximum:\t" + toString(pref[spot].score[0])); mothurOutEndLine();
+ m->mothurOut("Maximum:\t" + toString(pref[spot].score[0])); m->mothurOutEndLine();
+
+ return 1;
}
catch(exception& e) {
- errorOut(e, "Bellerophon", "print");
+ m->errorOut(e, "Bellerophon", "print");
exit(1);
}
}
}
//***************************************************************************************************************
-void Bellerophon::getChimeras() {
+int Bellerophon::getChimeras() {
try {
//do soft filter
if (filter) {
- string optionString = "fasta=" + fastafile + ", soft=50, vertical=F";
+ string optionString = "fasta=" + fastafile + ", soft=50";
+ if (outputDir != "") { optionString += ", outputdir=" + outputDir; }
+
filterSeqs = new FilterSeqsCommand(optionString);
filterSeqs->execute();
delete filterSeqs;
+ if (m->control_pressed) { return 0; }
+
//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);
+ seqs = readSeqs(fastafile);
+
+ if (m->control_pressed) { return 0; }
+
+ if (unaligned) { m->mothurOut("Your sequences need to be aligned when you use the bellerophon method."); m->mothurOutEndLine(); return 1; }
int numSeqs = seqs.size();
- if (numSeqs == 0) { mothurOut("Error in reading you sequences."); mothurOutEndLine(); exit(1); }
+ if (numSeqs == 0) { m->mothurOut("Error in reading you sequences."); m->mothurOutEndLine(); exit(1); }
//set default window to 25% of sequence length
- string seq0 = seqs[0].getAligned();
+ string seq0 = seqs[0]->getAligned();
if (window == 0) { window = seq0.length() / 4; }
else if (window > (seq0.length() / 2)) {
- mothurOut("Your sequence length is = " + toString(seq0.length()) + ". You have selected a window size greater than the length of half your aligned sequence. I will run it with a window size of " + toString((seq0.length() / 2))); mothurOutEndLine();
+ m->mothurOut("Your sequence length is = " + toString(seq0.length()) + ". You have selected a window size greater than the length of half your aligned sequence. I will run it with a window size of " + toString((seq0.length() / 2))); m->mothurOutEndLine();
window = (seq0.length() / 2);
}
- if (increment > (seqs[0].getAlignLength() - (2*window))) {
+ if (increment > (seqs[0]->getAlignLength() - (2*window))) {
if (increment != 10) {
- mothurOut("You have selected a increment that is too large. I will use the default."); mothurOutEndLine();
+ m->mothurOut("You have selected a increment that is too large. I will use the default."); m->mothurOutEndLine();
increment = 10;
- if (increment > (seqs[0].getAlignLength() - (2*window))) { increment = 0; }
+ if (increment > (seqs[0]->getAlignLength() - (2*window))) { increment = 0; }
}else{ increment = 0; }
}
-cout << "increment = " << increment << endl;
+
if (increment == 0) { iters = 1; }
- else { iters = ((seqs[0].getAlignLength() - (2*window)) / increment); }
+ else { iters = ((seqs[0]->getAlignLength() - (2*window)) / increment); }
//initialize pref
pref.resize(numSeqs);
for (int i = 0; i < numSeqs; i++ ) {
pref[i].leftParent.resize(2); pref[i].rightParent.resize(2); pref[i].score.resize(2); pref[i].closestLeft.resize(2); pref[i].closestRight.resize(3);
- pref[i].name = seqs[i].getName();
+ pref[i].name = seqs[i]->getName();
pref[i].score[0] = 0.0; pref[i].score[1] = 0.0;
pref[i].closestLeft[0] = 100000.0; pref[i].closestLeft[1] = 100000.0;
pref[i].closestRight[0] = 100000.0; pref[i].closestRight[1] = 100000.0;
int count = 0;
while (count < iters) {
+ if (m->control_pressed) { return 0; }
+
//create 2 vectors of sequences, 1 for left side and one for right side
vector<Sequence> left; vector<Sequence> right;
for (int i = 0; i < seqs.size(); i++) {
-//cout << "whole = " << seqs[i].getAligned() << endl;
+
+ if (m->control_pressed) { return 0; }
+
+//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);
+ string seqLeft = seqs[i]->getAligned().substr(midpoint-window, window);
Sequence tempLeft;
- tempLeft.setName(seqs[i].getName());
+ 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);
+ string seqRight = seqs[i]->getAligned().substr(midpoint, window);
Sequence tempRight;
- tempRight.setName(seqs[i].getName());
+ 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
SparseMatrix* SparseRight = new SparseMatrix();
createSparseMatrix(0, left.size(), SparseLeft, left);
+
+ if (m->control_pressed) { delete SparseLeft; delete SparseRight; return 0; }
+
createSparseMatrix(0, right.size(), SparseRight, right);
+ if (m->control_pressed) { delete SparseLeft; delete SparseRight; return 0; }
+
vector<SeqMap> distMapRight;
vector<SeqMap> distMapLeft;
//cout << "left" << endl << endl;
for (MatData currentCell = SparseLeft->begin(); currentCell != SparseLeft->end(); currentCell++) {
distMapLeft[currentCell->row][currentCell->column] = currentCell->dist;
+ if (m->control_pressed) { delete SparseLeft; delete SparseRight; return 0; }
//cout << " i = " << currentCell->row << " j = " << currentCell->column << " dist = " << currentCell->dist << endl;
}
//cout << "right" << endl << endl;
for (MatData currentCell = SparseRight->begin(); currentCell != SparseRight->end(); currentCell++) {
distMapRight[currentCell->row][currentCell->column] = currentCell->dist;
+ if (m->control_pressed) { delete SparseLeft; delete SparseRight; return 0; }
//cout << " i = " << currentCell->row << " j = " << currentCell->column << " dist = " << currentCell->dist << endl;
}
delete SparseLeft;
delete SparseRight;
-
//fill preference structure
generatePreferences(distMapLeft, distMapRight, midpoint);
//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) {
- errorOut(e, "Bellerophon", "getChimeras");
+ m->errorOut(e, "Bellerophon", "getChimeras");
exit(1);
}
}
for(int i=startSeq; i<endSeq; i++){
for(int j=0;j<i;j++){
+
+ if (m->control_pressed) { return 0; }
distCalculator->calcDist(s[i], s[j]);
float dist = distCalculator->getDist();
return 1;
}
catch(exception& e) {
- errorOut(e, "Bellerophon", "createSparseMatrix");
+ m->errorOut(e, "Bellerophon", "createSparseMatrix");
exit(1);
}
}
/***************************************************************************************************************/
-void Bellerophon::generatePreferences(vector<SeqMap> left, vector<SeqMap> right, int mid){
+int Bellerophon::generatePreferences(vector<SeqMap> left, vector<SeqMap> right, int mid){
try {
float dme = 0.0;
SeqMap currentRight = right[i]; // same as left but with distances on the right side.
for (int j = 0; j < i; j++) {
+
+ if (m->control_pressed) { return 0; }
itL = currentLeft.find(j);
itR = currentRight.find(j);
if (itL->second < pref[i].closestLeft[1]) {
pref[i].closestLeft[1] = itL->second;
- pref[i].leftParent[1] = seqs[j].getName();
+ pref[i].leftParent[1] = seqs[j]->getName();
//cout << "updating closest left to " << pref[i].leftParent[1] << endl;
}
//cout << "pref[" << j << "].closestLeft[1] = " << pref[j].closestLeft[1] << " parent = " << pref[j].leftParent[1] << endl;
if (itL->second < pref[j].closestLeft[1]) {
pref[j].closestLeft[1] = itL->second;
- pref[j].leftParent[1] = seqs[i].getName();
+ pref[j].leftParent[1] = seqs[i]->getName();
//cout << "updating closest left to " << pref[j].leftParent[1] << endl;
}
//are you the closest right sequence
if (itR->second < pref[i].closestRight[1]) {
pref[i].closestRight[1] = itR->second;
- pref[i].rightParent[1] = seqs[j].getName();
+ pref[i].rightParent[1] = seqs[j]->getName();
}
if (itR->second < pref[j].closestRight[1]) {
pref[j].closestRight[1] = itR->second;
- pref[j].rightParent[1] = seqs[i].getName();
+ pref[j].rightParent[1] = seqs[i]->getName();
}
}
//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++; } }
//cout << endl << "dme = " << dme << endl;
//recalculate prefernences based on dme
for (int i = 0; i < pref.size(); i++) {
+
+ if (m->control_pressed) { return 0; }
//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;
//is this score bigger then the last score
for (int i = 0; i < pref.size(); i++) {
+ if (m->control_pressed) { return 0; }
+
//update biggest score
if (pref[i].score[1] > pref[i].score[0]) {
pref[i].score[0] = pref[i].score[1];
}
}
+
+ return 1;
}
catch(exception& e) {
- errorOut(e, "Bellerophon", "generatePreferences");
+ m->errorOut(e, "Bellerophon", "generatePreferences");
exit(1);
}
}