5 * Created by Sarah Westcott on 7/9/09.
6 * Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
10 #include "bellerophon.h"
11 #include "eachgapdist.h"
12 #include "ignoregaps.h"
13 #include "onegapdist.h"
16 /***************************************************************************************************************/
18 Bellerophon::Bellerophon(string name, bool filterSeqs, bool c, int win, int inc, int p, string o) : Chimera() {
28 seqs = readSeqs(fastafile);
29 numSeqs = seqs.size();
30 if (numSeqs == 0) { m->mothurOut("Error in reading you sequences."); m->mothurOutEndLine(); exit(1); }
34 createFilter(seqs, 0.5);
35 for (int i = 0; i < seqs.size(); i++) { runFilter(seqs[i]); }
38 distCalculator = new eachGapDist();
40 //set default window to 25% of sequence length
41 string seq0 = seqs[0]->getAligned();
42 if (window == 0) { window = seq0.length() / 4; }
43 else if (window > (seq0.length() / 2)) {
44 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();
45 window = (seq0.length() / 2);
48 if (increment > (seqs[0]->getAlignLength() - (2*window))) {
49 if (increment != 10) {
51 m->mothurOut("You have selected a increment that is too large. I will use the default."); m->mothurOutEndLine();
53 if (increment > (seqs[0]->getAlignLength() - (2*window))) { increment = 0; }
55 }else{ increment = 0; }
58 if (increment == 0) { iters = 1; }
59 else { iters = ((seqs[0]->getAlignLength() - (2*window)) / increment); }
63 for (int i = 0; i < iters; i++) {
65 for (int j = 0; j < numSeqs; j++) {
66 pref[i].push_back(temp);
72 m->errorOut(e, "Bellerophon", "Bellerophon");
77 //***************************************************************************************************************
78 int Bellerophon::print(ostream& out, ostream& outAcc) {
82 //sorted "best" preference scores for all seqs
83 vector<Preference> best = getBestPref();
85 if (m->control_pressed) { return numSeqs; }
87 out << "Name\tScore\tLeft\tRight\t" << endl;
88 //output prefenence structure to .chimeras file
89 for (int i = 0; i < best.size(); i++) {
91 if (m->control_pressed) { return numSeqs; }
93 out << best[i].name << '\t' << setprecision(3) << best[i].score << '\t' << best[i].leftParent << '\t' << best[i].rightParent << endl;
95 //calc # of seqs with preference above 1.0
96 if (best[i].score > 1.0) {
98 outAcc << best[i].name << endl;
99 m->mothurOut(best[i].name + " is a suspected chimera at breakpoint " + toString(best[i].midpoint)); m->mothurOutEndLine();
100 m->mothurOut("It's score is " + toString(best[i].score) + " with suspected left parent " + best[i].leftParent + " and right parent " + best[i].rightParent); m->mothurOutEndLine();
104 //output results to screen
105 m->mothurOutEndLine();
106 m->mothurOut("Sequence with preference score above 1.0: " + toString(above1)); m->mothurOutEndLine();
108 spot = best.size()-1;
109 m->mothurOut("Minimum:\t" + toString(best[spot].score)); m->mothurOutEndLine();
110 spot = best.size() * 0.975;
111 m->mothurOut("2.5%-tile:\t" + toString(best[spot].score)); m->mothurOutEndLine();
112 spot = best.size() * 0.75;
113 m->mothurOut("25%-tile:\t" + toString(best[spot].score)); m->mothurOutEndLine();
114 spot = best.size() * 0.50;
115 m->mothurOut("Median: \t" + toString(best[spot].score)); m->mothurOutEndLine();
116 spot = best.size() * 0.25;
117 m->mothurOut("75%-tile:\t" + toString(best[spot].score)); m->mothurOutEndLine();
118 spot = best.size() * 0.025;
119 m->mothurOut("97.5%-tile:\t" + toString(best[spot].score)); m->mothurOutEndLine();
121 m->mothurOut("Maximum:\t" + toString(best[spot].score)); m->mothurOutEndLine();
126 catch(exception& e) {
127 m->errorOut(e, "Bellerophon", "print");
132 //***************************************************************************************************************
133 int Bellerophon::print(MPI_File& out, MPI_File& outAcc) {
137 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
140 string outString = "";
142 //sorted "best" preference scores for all seqs
143 vector<Preference> best = getBestPref();
146 int ninetyfive = best.size() * 0.05;
147 float cutoffScore = best[ninetyfive].score;
149 if (m->control_pressed) { return numSeqs; }
151 outString += "Name\tScore\tLeft\tRight\n";
152 //output prefenence structure to .chimeras file
153 for (int i = 0; i < best.size(); i++) {
155 if (m->control_pressed) { return numSeqs; }
157 outString += best[i].name + "\t" + toString(best[i].score) + "\t" + best[i].leftParent + "\t" + best[i].rightParent + "\n";
160 int length = outString.length();
161 char* buf2 = new char[length];
\r
162 memcpy(buf2, outString.c_str(), length);
164 MPI_File_write_shared(out, buf2, length, MPI_CHAR, &status);
168 //calc # of seqs with preference above 95%tile
169 if (best[i].score >= cutoffScore) {
172 outAccString += best[i].name + "\n";
174 MPI_Status statusAcc;
175 length = outAccString.length();
176 char* buf = new char[length];
\r
177 memcpy(buf, outAccString.c_str(), length);
179 MPI_File_write_shared(outAcc, buf, length, MPI_CHAR, &statusAcc);
183 cout << best[i].name << " is a suspected chimera at breakpoint " << toString(best[i].midpoint) << endl;
184 cout << "It's score is " << toString(best[i].score) << " with suspected left parent " << best[i].leftParent << " and right parent " << best[i].rightParent << endl;
188 //output results to screen
189 m->mothurOutEndLine();
190 m->mothurOut("Sequence with preference score above " + toString(cutoffScore) + ": " + toString(above1)); m->mothurOutEndLine();
192 spot = best.size()-1;
193 m->mothurOut("Minimum:\t" + toString(best[spot].score)); m->mothurOutEndLine();
194 spot = best.size() * 0.975;
195 m->mothurOut("2.5%-tile:\t" + toString(best[spot].score)); m->mothurOutEndLine();
196 spot = best.size() * 0.75;
197 m->mothurOut("25%-tile:\t" + toString(best[spot].score)); m->mothurOutEndLine();
198 spot = best.size() * 0.50;
199 m->mothurOut("Median: \t" + toString(best[spot].score)); m->mothurOutEndLine();
200 spot = best.size() * 0.25;
201 m->mothurOut("75%-tile:\t" + toString(best[spot].score)); m->mothurOutEndLine();
202 spot = best.size() * 0.025;
203 m->mothurOut("97.5%-tile:\t" + toString(best[spot].score)); m->mothurOutEndLine();
205 m->mothurOut("Maximum:\t" + toString(best[spot].score)); m->mothurOutEndLine();
212 catch(exception& e) {
213 m->errorOut(e, "Bellerophon", "print");
218 //********************************************************************************************************************
219 //sorts highest score to lowest
220 inline bool comparePref(Preference left, Preference right){
221 return (left.score > right.score);
223 //***************************************************************************************************************
224 int Bellerophon::getChimeras() {
227 //create breaking points
228 vector<int> midpoints; midpoints.resize(iters, window);
229 for (int i = 1; i < iters; i++) { midpoints[i] = midpoints[i-1] + increment; }
232 int pid, numSeqsPerProcessor;
234 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
235 MPI_Comm_size(MPI_COMM_WORLD, &processors);
237 numSeqsPerProcessor = iters / processors;
239 //each process hits this only once
240 int startPos = pid * numSeqsPerProcessor;
241 if(pid == processors - 1){
242 numSeqsPerProcessor = iters - pid * numSeqsPerProcessor;
244 lines.push_back(linePair(startPos, numSeqsPerProcessor));
246 //fill pref with scores
247 driverChimeras(midpoints, lines[0]);
249 if (m->control_pressed) { return 0; }
251 //each process must send its parts back to pid 0
255 for (int j = 1; j < processors; j++) {
257 vector<string> MPIBestSend;
258 for (int i = 0; i < numSeqs; i++) {
260 if (m->control_pressed) { return 0; }
265 MPI_Recv(&length, 1, MPI_INT, j, 2001, MPI_COMM_WORLD, &status);
267 char* buf = new char[length];
268 MPI_Recv(&buf, length, MPI_CHAR, j, 2001, MPI_COMM_WORLD, &status);
271 if (temp.length() > length) { temp = temp.substr(0, length); }
274 MPIBestSend.push_back(temp);
277 fillPref(j, MPIBestSend);
279 if (m->control_pressed) { return 0; }
283 //takes best window for each sequence and turns Preference to string that can be parsed by pid 0.
284 //played with this a bit, but it may be better to try user-defined datatypes with set string lengths??
285 vector<string> MPIBestSend = getBestWindow(lines[0]);
288 //send your result to parent
289 for (int i = 0; i < numSeqs; i++) {
291 if (m->control_pressed) { return 0; }
293 int bestLength = MPIBestSend[i].length();
294 char* buf = new char[bestLength];
\r
295 memcpy(buf, MPIBestSend[i].c_str(), bestLength);
297 MPI_Send(&bestLength, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD);
298 MPI_Send(buf, bestLength, MPI_CHAR, 0, 2001, MPI_COMM_WORLD);
307 //divide breakpoints between processors
308 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
310 lines.push_back(linePair(0, iters));
312 //fill pref with scores
313 driverChimeras(midpoints, lines[0]);
317 int numSeqsPerProcessor = iters / processors;
319 for (int i = 0; i < processors; i++) {
320 int startPos = i * numSeqsPerProcessor;
321 if(i == processors - 1){
322 numSeqsPerProcessor = iters - i * numSeqsPerProcessor;
324 lines.push_back(linePair(startPos, numSeqsPerProcessor));
327 createProcesses(midpoints);
330 lines.push_back(linePair(0, iters));
332 ///fill pref with scores
333 driverChimeras(midpoints, lines[0]);
341 catch(exception& e) {
342 m->errorOut(e, "Bellerophon", "getChimeras");
346 /**************************************************************************************************/
348 int Bellerophon::createProcesses(vector<int> mid) {
350 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
353 vector<int> processIDS;
355 //loop through and create all the processes you want
356 while (process != processors) {
360 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
363 exitCommand = driverChimeras(mid, lines[process]);
364 string tempOut = outputDir + toString(getpid()) + ".temp";
365 writePrefs(tempOut, lines[process]);
367 }else { m->mothurOut("unable to spawn the necessary processes."); m->mothurOutEndLine(); exit(0); }
370 //force parent to wait until all the processes are done
371 for (int i=0;i<processors;i++) {
372 int temp = processIDS[i];
376 //get info that other processes created
377 for (int i = 0; i < processIDS.size(); i++) {
378 string tempIn = outputDir + toString(processIDS[i]) + ".temp";
385 catch(exception& e) {
386 m->errorOut(e, "AlignCommand", "createProcesses");
390 //***************************************************************************************************************
391 int Bellerophon::driverChimeras(vector<int> midpoints, linePair line) {
394 for (int h = line.start; h < (line.start + line.num); h++) {
396 int midpoint = midpoints[h];
398 //initialize pref[count]
399 for (int i = 0; i < numSeqs; i++ ) {
400 pref[count][i].name = seqs[i]->getName();
401 pref[count][i].midpoint = midpoint;
404 if (m->control_pressed) { return 0; }
406 //create 2 vectors of sequences, 1 for left side and one for right side
407 vector<Sequence> left; vector<Sequence> right;
409 for (int i = 0; i < seqs.size(); i++) {
411 if (m->control_pressed) { return 0; }
413 //cout << "midpoint = " << midpoint << "\twindow = " << window << endl;
414 //cout << "whole = " << seqs[i]->getAligned().length() << endl;
416 string seqLeft = seqs[i]->getAligned().substr(midpoint-window, window);
418 tempLeft.setName(seqs[i]->getName());
419 tempLeft.setAligned(seqLeft);
420 left.push_back(tempLeft);
421 //cout << "left = " << tempLeft.getAligned().length() << endl;
423 string seqRight = seqs[i]->getAligned().substr(midpoint, window);
425 tempRight.setName(seqs[i]->getName());
426 tempRight.setAligned(seqRight);
427 right.push_back(tempRight);
428 //cout << "right = " << seqRight.length() << endl;
431 //this should be parallelized
432 //perference = sum of (| distance of my left to sequence j's left - distance of my right to sequence j's right | )
433 //create a matrix containing the distance from left to left and right to right
434 //calculate distances
435 SparseMatrix* SparseLeft = new SparseMatrix();
436 SparseMatrix* SparseRight = new SparseMatrix();
438 createSparseMatrix(0, left.size(), SparseLeft, left);
440 if (m->control_pressed) { delete SparseLeft; delete SparseRight; return 0; }
442 createSparseMatrix(0, right.size(), SparseRight, right);
444 if (m->control_pressed) { delete SparseLeft; delete SparseRight; return 0; }
446 left.clear(); right.clear();
447 vector<SeqMap> distMapRight;
448 vector<SeqMap> distMapLeft;
450 // Create a data structure to quickly access the distance information.
451 //this is from thallingers reimplementation on get.oturep
452 // It consists of a vector of distance maps, where each map contains
453 // all distances of a certain sequence. Vector and maps are accessed
454 // via the index of a sequence in the distance matrix
455 distMapRight = vector<SeqMap>(numSeqs);
456 distMapLeft = vector<SeqMap>(numSeqs);
457 //cout << "left" << endl << endl;
458 for (MatData currentCell = SparseLeft->begin(); currentCell != SparseLeft->end(); currentCell++) {
459 distMapLeft[currentCell->row][currentCell->column] = currentCell->dist;
460 if (m->control_pressed) { delete SparseLeft; delete SparseRight; return 0; }
461 //cout << " i = " << currentCell->row << " j = " << currentCell->column << " dist = " << currentCell->dist << endl;
463 //cout << "right" << endl << endl;
464 for (MatData currentCell = SparseRight->begin(); currentCell != SparseRight->end(); currentCell++) {
465 distMapRight[currentCell->row][currentCell->column] = currentCell->dist;
466 if (m->control_pressed) { delete SparseLeft; delete SparseRight; return 0; }
467 //cout << " i = " << currentCell->row << " j = " << currentCell->column << " dist = " << currentCell->dist << endl;
473 //fill preference structure
474 generatePreferences(distMapLeft, distMapRight, midpoint);
476 if (m->control_pressed) { return 0; }
479 if((h+1) % 10 == 0){ cout << "Processing sliding window: " << toString(h+1) << "\n"; m->mothurOutJustToLog("Processing sliding window: " + toString(h+1) + "\n") ; }
484 if((line.start + line.num) % 10 != 0){ cout << "Processing sliding window: " << toString(line.start + line.num) << "\n"; m->mothurOutJustToLog("Processing sliding window: " + toString(line.start + line.num) + "\n") ; }
489 catch(exception& e) {
490 m->errorOut(e, "Bellerophon", "driverChimeras");
495 /***************************************************************************************************************/
496 int Bellerophon::createSparseMatrix(int startSeq, int endSeq, SparseMatrix* sparse, vector<Sequence> s){
499 for(int i=startSeq; i<endSeq; i++){
501 for(int j=0;j<i;j++){
503 if (m->control_pressed) { return 0; }
505 distCalculator->calcDist(s[i], s[j]);
506 float dist = distCalculator->getDist();
508 PCell temp(i, j, dist);
509 sparse->addCell(temp);
516 catch(exception& e) {
517 m->errorOut(e, "Bellerophon", "createSparseMatrix");
521 /***************************************************************************************************************/
522 int Bellerophon::generatePreferences(vector<SeqMap> left, vector<SeqMap> right, int mid){
526 SeqMap::iterator itR;
527 SeqMap::iterator itL;
529 for (int i = 0; i < left.size(); i++) {
531 SeqMap currentLeft = left[i]; //example i = 3; currentLeft is a map of 0 to the distance of sequence 3 to sequence 0,
532 // 1 to the distance of sequence 3 to sequence 1,
533 // 2 to the distance of sequence 3 to sequence 2.
534 SeqMap currentRight = right[i]; // same as left but with distances on the right side.
536 for (int j = 0; j < i; j++) {
538 if (m->control_pressed) { return 0; }
540 itL = currentLeft.find(j);
541 itR = currentRight.find(j);
542 //cout << " i = " << i << " j = " << j << " distLeft = " << itL->second << endl;
543 //cout << " i = " << i << " j = " << j << " distright = " << itR->second << endl;
545 //if you can find this entry update the preferences
546 if ((itL != currentLeft.end()) && (itR != currentRight.end())) {
549 pref[count][i].score += abs((itL->second - itR->second));
550 pref[count][j].score += abs((itL->second - itR->second));
551 //cout << "left " << i << " " << j << " = " << itL->second << " right " << i << " " << j << " = " << itR->second << endl;
552 //cout << "abs = " << abs((itL->second - itR->second)) << endl;
553 //cout << i << " score = " << pref[i].score[1] << endl;
554 //cout << j << " score = " << pref[j].score[1] << endl;
556 pref[count][i].score += abs((sqrt(itL->second) - sqrt(itR->second)));
557 pref[count][j].score += abs((sqrt(itL->second) - sqrt(itR->second)));
558 //cout << "left " << i << " " << j << " = " << itL->second << " right " << i << " " << j << " = " << itR->second << endl;
559 //cout << "abs = " << abs((sqrt(itL->second) - sqrt(itR->second))) << endl;
560 //cout << i << " score = " << pref[i].score[1] << endl;
561 //cout << j << " score = " << pref[j].score[1] << endl;
563 //cout << "pref[" << i << "].closestLeft[1] = " << pref[i].closestLeft[1] << " parent = " << pref[i].leftParent[1] << endl;
564 //are you the closest left sequence
565 if (itL->second < pref[count][i].closestLeft) {
567 pref[count][i].closestLeft = itL->second;
568 pref[count][i].leftParent = seqs[j]->getName();
569 //cout << "updating closest left to " << pref[i].leftParent[1] << endl;
571 //cout << "pref[" << j << "].closestLeft[1] = " << pref[j].closestLeft[1] << " parent = " << pref[j].leftParent[1] << endl;
572 if (itL->second < pref[count][j].closestLeft) {
573 pref[count][j].closestLeft = itL->second;
574 pref[count][j].leftParent = seqs[i]->getName();
575 //cout << "updating closest left to " << pref[j].leftParent[1] << endl;
578 //are you the closest right sequence
579 if (itR->second < pref[count][i].closestRight) {
580 pref[count][i].closestRight = itR->second;
581 pref[count][i].rightParent = seqs[j]->getName();
583 if (itR->second < pref[count][j].closestRight) {
584 pref[count][j].closestRight = itR->second;
585 pref[count][j].rightParent = seqs[i]->getName();
597 catch(exception& e) {
598 m->errorOut(e, "Bellerophon", "generatePreferences");
602 /**************************************************************************************************/
603 vector<Preference> Bellerophon::getBestPref() {
606 vector<Preference> best;
609 for (int i = 0; i < numSeqs; i++) {
611 //set best pref score to first one
612 Preference temp = pref[0][i];
614 if (m->control_pressed) { return best; }
617 for (int j = 1; j < pref.size(); j++) {
619 //is this a better score
620 if (pref[j][i].score > temp.score) { temp = pref[j][i]; }
623 best.push_back(temp);
626 //rank preference score to eachother
628 float expectedPercent = 1 / (float) (best.size());
630 for (int i = 0; i < best.size(); i++) { dme += best[i].score; }
632 for (int i = 0; i < best.size(); i++) {
634 if (m->control_pressed) { return best; }
636 //gives the actual percentage of the dme this seq adds
637 best[i].score = best[i].score / dme;
639 //how much higher or lower is this than expected
640 best[i].score = best[i].score / expectedPercent;
644 //sort Preferences highest to lowest
645 sort(best.begin(), best.end(), comparePref);
649 catch(exception& e) {
650 m->errorOut(e, "Bellerophon", "getBestPref");
654 /**************************************************************************************************/
655 int Bellerophon::writePrefs(string file, linePair tempLine) {
659 openOutputFile(file, outTemp);
661 //lets you know what part of the pref matrix you are writing
662 outTemp << tempLine.start << '\t' << tempLine.num << endl;
664 for (int i = tempLine.start; i < (tempLine.start + tempLine.num); i++) {
666 for (int j = 0; j < numSeqs; j++) {
668 if (m->control_pressed) { outTemp.close(); remove(file.c_str()); return 0; }
670 outTemp << pref[i][j].name << '\t' << pref[i][j].leftParent << '\t' << pref[i][j].rightParent << '\t';
671 outTemp << pref[i][j].score << '\t' << pref[i][j].closestLeft << '\t' << pref[i][j].closestRight << '\t' << pref[i][j].midpoint << endl;
679 catch(exception& e) {
680 m->errorOut(e, "Bellerophon", "writePrefs");
684 /**************************************************************************************************/
685 int Bellerophon::readPrefs(string file) {
689 openInputFile(file, inTemp);
693 //lets you know what part of the pref matrix you are writing
694 inTemp >> start >> num; gobble(inTemp);
696 for (int i = start; i < num; i++) {
698 for (int j = 0; j < numSeqs; j++) {
700 if (m->control_pressed) { inTemp.close(); remove(file.c_str()); return 0; }
702 inTemp >> pref[i][j].name >> pref[i][j].leftParent >> pref[i][j].rightParent;
703 inTemp >> pref[i][j].score >> pref[i][j].closestLeft >> pref[i][j].closestRight >> pref[i][j].midpoint;
710 remove(file.c_str());
714 catch(exception& e) {
715 m->errorOut(e, "Bellerophon", "writePrefs");
719 /**************************************************************************************************/
720 vector<string> Bellerophon::getBestWindow(linePair line) {
726 for (int i = 0; i < numSeqs; i++) {
728 //set best pref score to first one
729 Preference temp = pref[line.start][i];
731 if (m->control_pressed) { return best; }
734 for (int j = (line.start+1); j < (line.start+line.num); j++) {
736 //is this a better score
737 if (pref[j][i].score > temp.score) { temp = pref[j][i]; }
740 string tempString = temp.name + '\t' + temp.leftParent + '\t' + temp.rightParent + '\t' + toString(temp.score);
741 best.push_back(tempString);
747 catch(exception& e) {
748 m->errorOut(e, "Bellerophon", "getBestWindow");
752 /**************************************************************************************************/
753 int Bellerophon::fillPref(int process, vector<string>& best) {
755 //figure out where you start so you can put the best scores there
756 int numSeqsPerProcessor = iters / processors;
757 int start = process * numSeqsPerProcessor;
759 for (int i = 0; i < best.size(); i++) {
761 if (m->control_pressed) { return 0; }
763 istringstream iss (best[i],istringstream::in);
766 iss >> pref[start][i].name >> pref[start][i].leftParent >> pref[start][i].rightParent >> tempScore;
767 convert(tempScore, pref[start][i].score);
772 catch(exception& e) {
773 m->errorOut(e, "Bellerophon", "fillPref");
778 /**************************************************************************************************/