bellerophon.cpp

   1 /*
   2  *  bellerophon.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 7/9/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "bellerophon.h"
  11 #include "eachgapdist.h"
  12 #include "ignoregaps.h"
  13 #include "onegapdist.h"
  14
  15
  16 /***************************************************************************************************************/
  17
  18 Bellerophon::Bellerophon(string name, bool filterSeqs,  bool c, int win, int inc, int p, string o) : Chimera() {
  19         try {
  20                 fastafile = name;
  21                 correction = c;
  22                 outputDir = o;
  23                 window = win;
  24                 increment = inc;
  25                 processors = p;
  26
  27                 //read in sequences
  28                 seqs = readSeqs(fastafile);
  29                 numSeqs = seqs.size();
  30                 if (numSeqs == 0) { m->mothurOut("Error in reading you sequences."); m->mothurOutEndLine(); exit(1); }
  31
  32                 //do soft filter
  33                 if (filterSeqs)  {
  34                         createFilter(seqs, 0.5);
  35                         for (int i = 0; i < seqs.size(); i++) {  runFilter(seqs[i]);  }
  36                 }
  37
  38                 distCalculator = new eachGapDist();
  39
  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);
  46                 }
  47
  48                 if (increment > (seqs[0]->getAlignLength() - (2*window))) {
  49                         if (increment != 10) {
  50
  51                                 m->mothurOut("You have selected a increment that is too large. I will use the default."); m->mothurOutEndLine();
  52                                 increment = 10;
  53                                 if (increment > (seqs[0]->getAlignLength() - (2*window))) {  increment = 0;  }
  54
  55                         }else{ increment = 0; }
  56                 }
  57
  58                 if (increment == 0) { iters = 1; }
  59                 else { iters = ((seqs[0]->getAlignLength() - (2*window)) / increment); }
  60
  61                 //initialize pref
  62                 pref.resize(iters);
  63                 for (int i = 0; i < iters; i++) {
  64                         Preference temp;
  65                         for (int j = 0; j < numSeqs; j++) {
  66                                 pref[i].push_back(temp);
  67                         }
  68                 }
  69
  70         }
  71         catch(exception& e) {
  72                 m->errorOut(e, "Bellerophon", "Bellerophon");
  73                 exit(1);
  74         }
  75 }
  76
  77 //***************************************************************************************************************
  78 int Bellerophon::print(ostream& out, ostream& outAcc, string s) {
  79         try {
  80                 int above1 = 0;
  81
  82                 //sorted "best" preference scores for all seqs
  83                 vector<Preference> best = getBestPref();
  84
  85                 if (m->control_pressed) { return numSeqs; }
  86
  87                 out << "Name\tScore\tLeft\tRight\t" << endl;
  88                 //output prefenence structure to .chimeras file
  89                 for (int i = 0; i < best.size(); i++) {
  90
  91                         if (m->control_pressed) {  return numSeqs; }
  92
  93                         out << best[i].name << '\t' << setprecision(3) << best[i].score << '\t' << best[i].leftParent << '\t' << best[i].rightParent << endl;
  94
  95                         //calc # of seqs with preference above 1.0
  96                         if (best[i].score > 1.0) {
  97                                 above1++;
  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();
 101                         }
 102                 }
 103
 104                 //output results to screen
 105                 m->mothurOutEndLine();
 106                 m->mothurOut("Sequence with preference score above 1.0: " + toString(above1)); m->mothurOutEndLine();
 107                 int spot;
 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();
 120                 spot = 0;
 121                 m->mothurOut("Maximum:\t" + toString(best[spot].score)); m->mothurOutEndLine();
 122
 123                 return numSeqs;
 124
 125         }
 126         catch(exception& e) {
 127                 m->errorOut(e, "Bellerophon", "print");
 128                 exit(1);
 129         }
 130 }
 131 #ifdef USE_MPI
 132 //***************************************************************************************************************
 133 int Bellerophon::print(MPI_File& out, MPI_File& outAcc, string s) {
 134         try {
 135
 136                 int pid;
 137                 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
 138
 139                 if (pid == 0) {
 140                         string outString = "";
 141
 142                         //sorted "best" preference scores for all seqs
 143                         vector<Preference> best = getBestPref();
 144
 145                         int above1 = 0;
 146                         int ninetyfive = best.size() * 0.05;
 147                         float cutoffScore = best[ninetyfive].score;
 148
 149                         if (m->control_pressed) { return numSeqs; }
 150
 151                         outString = "Name\tScore\tLeft\tRight\n";
 152                         MPI_Status status;
 153                         int olength = outString.length();
 154                         char* buf5 = new char[olength];
 155                         memcpy(buf5, outString.c_str(), olength);
 156
 157                         MPI_File_write_shared(out, buf5, olength, MPI_CHAR, &status);
 158
 159                         delete buf5;
 160
 161                         //output prefenence structure to .chimeras file
 162                         for (int i = 0; i < best.size(); i++) {
 163
 164                                 if (m->control_pressed) {  return numSeqs; }
 165
 166                                 outString = best[i].name + "\t" +  toString(best[i].score) + "\t" + best[i].leftParent + "\t" + best[i].rightParent + "\n";
 167
 168                                 MPI_Status status;
 169                                 int length = outString.length();
 170                                 char* buf2 = new char[length];
 171                                 memcpy(buf2, outString.c_str(), length);
 172
 173                                 MPI_File_write_shared(out, buf2, length, MPI_CHAR, &status);
 174
 175                                 delete buf2;
 176
 177                                 //calc # of seqs with preference above 95%tile
 178                                 if (best[i].score >= cutoffScore) {
 179                                         above1++;
 180                                         string outAccString = "";
 181                                          outAccString += best[i].name + "\n";
 182
 183                                         MPI_Status statusAcc;
 184                                         length = outAccString.length();
 185                                         char* buf = new char[length];
 186                                         memcpy(buf, outAccString.c_str(), length);
 187
 188                                         MPI_File_write_shared(outAcc, buf, length, MPI_CHAR, &statusAcc);
 189
 190                                         delete buf;
 191
 192                                         cout << best[i].name << " is a suspected chimera at breakpoint " << toString(best[i].midpoint) << endl;
 193                                         cout << "It's score is " << toString(best[i].score) << " with suspected left parent " << best[i].leftParent << " and right parent " << best[i].rightParent << endl;
 194                                 }
 195                         }
 196
 197                         //output results to screen
 198                         m->mothurOutEndLine();
 199                         m->mothurOut("Sequence with preference score above " + toString(cutoffScore) +  ": " + toString(above1)); m->mothurOutEndLine();
 200                         int spot;
 201                         spot = best.size()-1;
 202                         m->mothurOut("Minimum:\t" + toString(best[spot].score)); m->mothurOutEndLine();
 203                         spot = best.size() * 0.975;
 204                         m->mothurOut("2.5%-tile:\t" + toString(best[spot].score)); m->mothurOutEndLine();
 205                         spot = best.size() * 0.75;
 206                         m->mothurOut("25%-tile:\t" + toString(best[spot].score)); m->mothurOutEndLine();
 207                         spot = best.size() * 0.50;
 208                         m->mothurOut("Median: \t" + toString(best[spot].score)); m->mothurOutEndLine();
 209                         spot = best.size() * 0.25;
 210                         m->mothurOut("75%-tile:\t" + toString(best[spot].score)); m->mothurOutEndLine();
 211                         spot = best.size() * 0.025;
 212                         m->mothurOut("97.5%-tile:\t" + toString(best[spot].score)); m->mothurOutEndLine();
 213                         spot = 0;
 214                         m->mothurOut("Maximum:\t" + toString(best[spot].score)); m->mothurOutEndLine();
 215
 216                 }
 217
 218                 return numSeqs;
 219
 220         }
 221         catch(exception& e) {
 222                 m->errorOut(e, "Bellerophon", "print");
 223                 exit(1);
 224         }
 225 }
 226 #endif
 227 //********************************************************************************************************************
 228 //sorts highest score to lowest
 229 inline bool comparePref(Preference left, Preference right){
 230         return (left.score > right.score);
 231 }
 232 //***************************************************************************************************************
 233 int Bellerophon::getChimeras() {
 234         try {
 235
 236                 //create breaking points
 237                 vector<int> midpoints;   midpoints.resize(iters, window);
 238                 for (int i = 1; i < iters; i++) {  midpoints[i] = midpoints[i-1] + increment;  }
 239
 240         #ifdef USE_MPI
 241                 int pid, numSeqsPerProcessor;
 242
 243                 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
 244                 MPI_Comm_size(MPI_COMM_WORLD, &processors);
 245
 246                 numSeqsPerProcessor = iters / processors;
 247
 248                 //each process hits this only once
 249                 unsigned long int startPos = pid * numSeqsPerProcessor;
 250                 if(pid == processors - 1){
 251                                 numSeqsPerProcessor = iters - pid * numSeqsPerProcessor;
 252                 }
 253                 lines.push_back(linePair(startPos, numSeqsPerProcessor));
 254
 255                 //fill pref with scores
 256                 driverChimeras(midpoints, lines[0]);
 257
 258                 if (m->control_pressed) { return 0; }
 259
 260                 //each process must send its parts back to pid 0
 261                 if (pid == 0) {
 262
 263                         //receive results
 264                         for (int j = 1; j < processors; j++) {
 265
 266                                 vector<string>  MPIBestSend;
 267                                 for (int i = 0; i < numSeqs; i++) {
 268
 269                                         if (m->control_pressed) { return 0; }
 270
 271                                         MPI_Status status;
 272                                         //receive string
 273                                         int length;
 274                                         MPI_Recv(&length, 1, MPI_INT, j, 2001, MPI_COMM_WORLD, &status);
 275
 276                                         char* buf = new char[length];
 277                                         MPI_Recv(&buf, length, MPI_CHAR, j, 2001, MPI_COMM_WORLD, &status);
 278
 279                                         string temp = buf;
 280                                         if (temp.length() > length) { temp = temp.substr(0, length); }
 281                                         delete buf;
 282
 283                                         MPIBestSend.push_back(temp);
 284                                 }
 285
 286                                 fillPref(j, MPIBestSend);
 287
 288                                 if (m->control_pressed) { return 0; }
 289                         }
 290
 291                 }else {
 292                         //takes best window for each sequence and turns Preference to string that can be parsed by pid 0.
 293                         //played with this a bit, but it may be better to try user-defined datatypes with set string lengths??
 294                         vector<string> MPIBestSend = getBestWindow(lines[0]);
 295                         pref.clear();
 296
 297                         //send your result to parent
 298                         for (int i = 0; i < numSeqs; i++) {
 299
 300                                 if (m->control_pressed) { return 0; }
 301
 302                                 int bestLength = MPIBestSend[i].length();
 303                                 char* buf = new char[bestLength];
 304                                 memcpy(buf, MPIBestSend[i].c_str(), bestLength);
 305
 306                                 MPI_Send(&bestLength, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD);
 307                                 MPI_Send(buf, bestLength, MPI_CHAR, 0, 2001, MPI_COMM_WORLD);
 308                                 delete buf;
 309                         }
 310
 311                         MPIBestSend.clear();
 312                 }
 313                 MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
 314         #else
 315
 316                 //divide breakpoints between processors
 317                 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 318                         if(processors == 1){
 319                                 lines.push_back(linePair(0, iters));
 320
 321                                 //fill pref with scores
 322                                 driverChimeras(midpoints, lines[0]);
 323
 324                         }else{
 325
 326                                 int numSeqsPerProcessor = iters / processors;
 327
 328                                 for (int i = 0; i < processors; i++) {
 329                                         unsigned long int startPos = i * numSeqsPerProcessor;
 330                                         if(i == processors - 1){
 331                                                 numSeqsPerProcessor = iters - i * numSeqsPerProcessor;
 332                                         }
 333                                         lines.push_back(linePair(startPos, numSeqsPerProcessor));
 334                                 }
 335
 336                                 createProcesses(midpoints);
 337                         }
 338                 #else
 339                         lines.push_back(linePair(0, iters));
 340
 341                         ///fill pref with scores
 342                         driverChimeras(midpoints, lines[0]);
 343                 #endif
 344
 345         #endif
 346
 347                 return 0;
 348
 349         }
 350         catch(exception& e) {
 351                 m->errorOut(e, "Bellerophon", "getChimeras");
 352                 exit(1);
 353         }
 354 }
 355 /**************************************************************************************************/
 356
 357 int Bellerophon::createProcesses(vector<int> mid) {
 358         try {
 359 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 360                 int process = 0;
 361                 int exitCommand = 1;
 362                 vector<int> processIDS;
 363
 364                 //loop through and create all the processes you want
 365                 while (process != processors) {
 366                         int pid = fork();
 367
 368                         if (pid > 0) {
 369                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
 370                                 process++;
 371                         }else if (pid == 0){
 372                                 exitCommand = driverChimeras(mid, lines[process]);
 373                                 string tempOut = outputDir + toString(getpid()) + ".temp";
 374                                 writePrefs(tempOut, lines[process]);
 375                                 exit(0);
 376                         }else {
 377                                 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
 378                                 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
 379                                 exit(0);
 380                         }
 381                 }
 382
 383                 //force parent to wait until all the processes are done
 384                 for (int i=0;i<processors;i++) {
 385                         int temp = processIDS[i];
 386                         wait(&temp);
 387                 }
 388
 389                 //get info that other processes created
 390                 for (int i = 0; i < processIDS.size(); i++) {
 391                         string tempIn = outputDir + toString(processIDS[i]) + ".temp";
 392                         readPrefs(tempIn);
 393                 }
 394
 395                 return exitCommand;
 396 #endif
 397         }
 398         catch(exception& e) {
 399                 m->errorOut(e, "AlignCommand", "createProcesses");
 400                 exit(1);
 401         }
 402 }
 403 //***************************************************************************************************************
 404 int Bellerophon::driverChimeras(vector<int> midpoints, linePair line) {
 405         try {
 406
 407                 for (int h = line.start; h < (line.start + line.num); h++) {
 408                         count = h;
 409                         int midpoint = midpoints[h];
 410
 411                         //initialize pref[count]
 412                         for (int i = 0; i < numSeqs; i++ ) {
 413                                 pref[count][i].name = seqs[i]->getName();
 414                                 pref[count][i].midpoint = midpoint;
 415                         }
 416
 417                         if (m->control_pressed) { return 0; }
 418
 419                         //create 2 vectors of sequences, 1 for left side and one for right side
 420                         vector<Sequence> left;  vector<Sequence> right;
 421
 422                         for (int i = 0; i < seqs.size(); i++) {
 423
 424                                 if (m->control_pressed) { return 0; }
 425
 426                                 //cout << "midpoint = " << midpoint << "\twindow = " << window << endl;
 427                                 //cout << "whole = " << seqs[i]->getAligned().length() << endl;
 428                                 //save left side
 429                                 string seqLeft = seqs[i]->getAligned().substr(midpoint-window, window);
 430                                 Sequence tempLeft;
 431                                 tempLeft.setName(seqs[i]->getName());
 432                                 tempLeft.setAligned(seqLeft);
 433                                 left.push_back(tempLeft);
 434                                 //cout << "left = " << tempLeft.getAligned().length() << endl;
 435                                 //save right side
 436                                 string seqRight = seqs[i]->getAligned().substr(midpoint, window);
 437                                 Sequence tempRight;
 438                                 tempRight.setName(seqs[i]->getName());
 439                                 tempRight.setAligned(seqRight);
 440                                 right.push_back(tempRight);
 441                                 //cout << "right = " << seqRight.length() << endl;
 442                         }
 443
 444                         //this should be parallelized
 445                         //perference = sum of (| distance of my left to sequence j's left - distance of my right to sequence j's right | )
 446                         //create a matrix containing the distance from left to left and right to right
 447                         //calculate distances
 448                         SparseMatrix* SparseLeft = new SparseMatrix();
 449                         SparseMatrix* SparseRight = new SparseMatrix();
 450
 451                         createSparseMatrix(0, left.size(), SparseLeft, left);
 452
 453                         if (m->control_pressed) { delete SparseLeft; delete SparseRight; return 0; }
 454
 455                         createSparseMatrix(0, right.size(), SparseRight, right);
 456
 457                         if (m->control_pressed) { delete SparseLeft; delete SparseRight; return 0; }
 458
 459                         left.clear(); right.clear();
 460                         vector<SeqMap> distMapRight;
 461                         vector<SeqMap> distMapLeft;
 462
 463                         // Create a data structure to quickly access the distance information.
 464                         //this is from thallingers reimplementation on get.oturep
 465                         // It consists of a vector of distance maps, where each map contains
 466                         // all distances of a certain sequence. Vector and maps are accessed
 467                         // via the index of a sequence in the distance matrix
 468                         distMapRight = vector<SeqMap>(numSeqs);
 469                         distMapLeft = vector<SeqMap>(numSeqs);
 470                         //cout << "left" << endl << endl;
 471                         for (MatData currentCell = SparseLeft->begin(); currentCell != SparseLeft->end(); currentCell++) {
 472                                 distMapLeft[currentCell->row][currentCell->column] = currentCell->dist;
 473                                 if (m->control_pressed) { delete SparseLeft; delete SparseRight; return 0; }
 474                                 //cout << " i = " << currentCell->row << " j = " << currentCell->column << " dist = " << currentCell->dist << endl;
 475                         }
 476                         //cout << "right" << endl << endl;
 477                         for (MatData currentCell = SparseRight->begin(); currentCell != SparseRight->end(); currentCell++) {
 478                                 distMapRight[currentCell->row][currentCell->column] = currentCell->dist;
 479                                 if (m->control_pressed) { delete SparseLeft; delete SparseRight; return 0; }
 480                                 //cout << " i = " << currentCell->row << " j = " << currentCell->column << " dist = " << currentCell->dist << endl;
 481                         }
 482
 483                         delete SparseLeft;
 484                         delete SparseRight;
 485
 486                         //fill preference structure
 487                         generatePreferences(distMapLeft, distMapRight, midpoint);
 488
 489                         if (m->control_pressed) { return 0; }
 490
 491                         //report progress
 492                         if((h+1) % 10 == 0){    cout << "Processing sliding window: " << toString(h+1) <<  "\n";  m->mothurOutJustToLog("Processing sliding window: " + toString(h+1) + "\n") ;         }
 493
 494                 }
 495
 496                 //report progress
 497                 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") ;             }
 498
 499                 return 0;
 500
 501         }
 502         catch(exception& e) {
 503                 m->errorOut(e, "Bellerophon", "driverChimeras");
 504                 exit(1);
 505         }
 506 }
 507
 508 /***************************************************************************************************************/
 509 int Bellerophon::createSparseMatrix(int startSeq, int endSeq, SparseMatrix* sparse, vector<Sequence> s){
 510         try {
 511
 512                 for(int i=startSeq; i<endSeq; i++){
 513
 514                         for(int j=0;j<i;j++){
 515
 516                                 if (m->control_pressed) { return 0; }
 517
 518                                 distCalculator->calcDist(s[i], s[j]);
 519                                 float dist = distCalculator->getDist();
 520
 521                                 PCell temp(i, j, dist);
 522                                 sparse->addCell(temp);
 523
 524                         }
 525                 }
 526
 527                 return 1;
 528         }
 529         catch(exception& e) {
 530                 m->errorOut(e, "Bellerophon", "createSparseMatrix");
 531                 exit(1);
 532         }
 533 }
 534 /***************************************************************************************************************/
 535 int Bellerophon::generatePreferences(vector<SeqMap> left, vector<SeqMap> right, int mid){
 536         try {
 537
 538                 SeqMap::iterator itR;
 539                 SeqMap::iterator itL;
 540
 541                 for (int i = 0; i < left.size(); i++) {
 542
 543                         SeqMap currentLeft = left[i];    //example i = 3;   currentLeft is a map of 0 to the distance of sequence 3 to sequence 0,
 544                                                                                                 //                                                                              1 to the distance of sequence 3 to sequence 1,
 545                                                                                                 //                                                                              2 to the distance of sequence 3 to sequence 2.
 546                         SeqMap currentRight = right[i];         // same as left but with distances on the right side.
 547
 548                         for (int j = 0; j < i; j++) {
 549
 550                                 if (m->control_pressed) {  return 0; }
 551
 552                                 itL = currentLeft.find(j);
 553                                 itR = currentRight.find(j);
 554 //cout << " i = " << i << " j = " << j << " distLeft = " << itL->second << endl;
 555 //cout << " i = " << i << " j = " << j << " distright = " << itR->second << endl;
 556
 557                                 //if you can find this entry update the preferences
 558                                 if ((itL != currentLeft.end()) && (itR != currentRight.end())) {
 559
 560                                         if (!correction) {
 561                                                 pref[count][i].score += abs((itL->second - itR->second));
 562                                                 pref[count][j].score += abs((itL->second - itR->second));
 563 //cout << "left " << i << " " << j << " = " << itL->second << " right " << i << " " << j << " = " << itR->second << endl;
 564 //cout << "abs = " << abs((itL->second - itR->second)) << endl;
 565 //cout << i << " score = " << pref[i].score[1] << endl;
 566 //cout << j << " score = " << pref[j].score[1] << endl;
 567                                         }else {
 568                                                 pref[count][i].score += abs((sqrt(itL->second) - sqrt(itR->second)));
 569                                                 pref[count][j].score += abs((sqrt(itL->second) - sqrt(itR->second)));
 570 //cout << "left " << i << " " << j << " = " << itL->second << " right " << i << " " << j << " = " << itR->second << endl;
 571 //cout << "abs = " << abs((sqrt(itL->second) - sqrt(itR->second))) << endl;
 572 //cout << i << " score = " << pref[i].score[1] << endl;
 573 //cout << j << " score = " << pref[j].score[1] << endl;
 574                                         }
 575 //cout << "pref[" << i << "].closestLeft[1] = " <<      pref[i].closestLeft[1] << " parent = " << pref[i].leftParent[1] << endl;
 576                                         //are you the closest left sequence
 577                                         if (itL->second < pref[count][i].closestLeft) {
 578
 579                                                 pref[count][i].closestLeft = itL->second;
 580                                                 pref[count][i].leftParent = seqs[j]->getName();
 581 //cout << "updating closest left to " << pref[i].leftParent[1] << endl;
 582                                         }
 583 //cout << "pref[" << j << "].closestLeft[1] = " <<      pref[j].closestLeft[1] << " parent = " << pref[j].leftParent[1] << endl;
 584                                         if (itL->second < pref[count][j].closestLeft) {
 585                                                 pref[count][j].closestLeft = itL->second;
 586                                                 pref[count][j].leftParent = seqs[i]->getName();
 587 //cout << "updating closest left to " << pref[j].leftParent[1] << endl;
 588                                         }
 589
 590                                         //are you the closest right sequence
 591                                         if (itR->second < pref[count][i].closestRight) {
 592                                                 pref[count][i].closestRight = itR->second;
 593                                                 pref[count][i].rightParent = seqs[j]->getName();
 594                                         }
 595                                         if (itR->second < pref[count][j].closestRight) {
 596                                                 pref[count][j].closestRight = itR->second;
 597                                                 pref[count][j].rightParent = seqs[i]->getName();
 598                                         }
 599
 600                                 }
 601                         }
 602
 603                 }
 604
 605
 606                 return 1;
 607
 608         }
 609         catch(exception& e) {
 610                 m->errorOut(e, "Bellerophon", "generatePreferences");
 611                 exit(1);
 612         }
 613 }
 614 /**************************************************************************************************/
 615 vector<Preference> Bellerophon::getBestPref() {
 616         try {
 617
 618                 vector<Preference> best;
 619
 620                 //for each sequence
 621                 for (int i = 0; i < numSeqs; i++) {
 622
 623                         //set best pref score to first one
 624                         Preference temp = pref[0][i];
 625
 626                         if (m->control_pressed) { return best;  }
 627
 628                         //for each window
 629                         for (int j = 1; j < pref.size(); j++) {
 630
 631                                 //is this a better score
 632                                 if (pref[j][i].score > temp.score) {    temp = pref[j][i];              }
 633                         }
 634
 635                         best.push_back(temp);
 636                 }
 637
 638                 //rank preference score to eachother
 639                 float dme = 0.0;
 640                 float expectedPercent = 1 / (float) (best.size());
 641
 642                 for (int i = 0; i < best.size(); i++) {  dme += best[i].score;  }
 643
 644                 for (int i = 0; i < best.size(); i++) {
 645
 646                         if (m->control_pressed) { return best; }
 647
 648                         //gives the actual percentage of the dme this seq adds
 649                         best[i].score = best[i].score / dme;
 650
 651                         //how much higher or lower is this than expected
 652                         best[i].score = best[i].score / expectedPercent;
 653
 654                 }
 655
 656                 //sort Preferences highest to lowest
 657                 sort(best.begin(), best.end(), comparePref);
 658
 659                 return best;
 660         }
 661         catch(exception& e) {
 662                 m->errorOut(e, "Bellerophon", "getBestPref");
 663                 exit(1);
 664         }
 665 }
 666 /**************************************************************************************************/
 667 int Bellerophon::writePrefs(string file, linePair tempLine) {
 668         try {
 669
 670                 ofstream outTemp;
 671                 m->openOutputFile(file, outTemp);
 672
 673                 //lets you know what part of the pref matrix you are writing
 674                 outTemp << tempLine.start << '\t' << tempLine.num << endl;
 675
 676                 for (int i = tempLine.start; i < (tempLine.start + tempLine.num); i++) {
 677
 678                         for (int j = 0; j < numSeqs; j++) {
 679
 680                                 if (m->control_pressed) { outTemp.close(); remove(file.c_str()); return 0; }
 681
 682                                 outTemp << pref[i][j].name << '\t' << pref[i][j].leftParent << '\t' << pref[i][j].rightParent << '\t';
 683                                 outTemp << pref[i][j].score << '\t' << pref[i][j].closestLeft << '\t' << pref[i][j].closestRight << '\t' << pref[i][j].midpoint <<  endl;
 684                         }
 685                 }
 686
 687                 outTemp.close();
 688
 689                 return 0;
 690         }
 691         catch(exception& e) {
 692                 m->errorOut(e, "Bellerophon", "writePrefs");
 693                 exit(1);
 694         }
 695 }
 696 /**************************************************************************************************/
 697 int Bellerophon::readPrefs(string file) {
 698         try {
 699
 700                 ifstream inTemp;
 701                 m->openInputFile(file, inTemp);
 702
 703                 int start, num;
 704
 705                 //lets you know what part of the pref matrix you are writing
 706                 inTemp >> start >> num;  m->gobble(inTemp);
 707
 708                 for (int i = start; i < num; i++) {
 709
 710                         for (int j = 0; j < numSeqs; j++) {
 711
 712                                 if (m->control_pressed) { inTemp.close(); remove(file.c_str()); return 0; }
 713
 714                                 inTemp >> pref[i][j].name >> pref[i][j].leftParent >> pref[i][j].rightParent;
 715                                 inTemp >> pref[i][j].score >> pref[i][j].closestLeft >> pref[i][j].closestRight >> pref[i][j].midpoint;
 716                                 m->gobble(inTemp);
 717                         }
 718                 }
 719
 720                 inTemp.close();
 721
 722                 remove(file.c_str());
 723
 724                 return 0;
 725         }
 726         catch(exception& e) {
 727                 m->errorOut(e, "Bellerophon", "writePrefs");
 728                 exit(1);
 729         }
 730 }
 731 /**************************************************************************************************/
 732 vector<string> Bellerophon::getBestWindow(linePair line) {
 733         try {
 734
 735                 vector<string> best;
 736
 737                 //for each sequence
 738                 for (int i = 0; i < numSeqs; i++) {
 739
 740                         //set best pref score to first one
 741                         Preference temp = pref[line.start][i];
 742
 743                         if (m->control_pressed) { return best;  }
 744
 745                         //for each window
 746                         for (int j = (line.start+1); j < (line.start+line.num); j++) {
 747
 748                                 //is this a better score
 749                                 if (pref[j][i].score > temp.score) {    temp = pref[j][i];              }
 750                         }
 751
 752                         string tempString = temp.name + '\t' + temp.leftParent + '\t' + temp.rightParent + '\t' + toString(temp.score);
 753                         best.push_back(tempString);
 754                 }
 755
 756                 return best;
 757
 758         }
 759         catch(exception& e) {
 760                 m->errorOut(e, "Bellerophon", "getBestWindow");
 761                 exit(1);
 762         }
 763 }
 764 /**************************************************************************************************/
 765 int Bellerophon::fillPref(int process, vector<string>& best) {
 766         try {
 767                 //figure out where you start so you can put the best scores there
 768                 int numSeqsPerProcessor = iters / processors;
 769                 int start = process * numSeqsPerProcessor;
 770
 771                 for (int i = 0; i < best.size(); i++) {
 772
 773                         if (m->control_pressed) { return 0;  }
 774
 775                         istringstream iss (best[i],istringstream::in);
 776
 777                         string tempScore;
 778                         iss >> pref[start][i].name >> pref[start][i].leftParent >> pref[start][i].rightParent >> tempScore;
 779                         convert(tempScore, pref[start][i].score);
 780                 }
 781
 782                 return 0;
 783         }
 784         catch(exception& e) {
 785                 m->errorOut(e, "Bellerophon", "fillPref");
 786                 exit(1);
 787         }
 788 }
 789
 790 /**************************************************************************************************/
 791