chimeraslayer.cpp

   1 /*
   2  *  chimeraslayer.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 9/25/09.
   6  *  Copyright 2009 Pschloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "chimeraslayer.h"
  11 #include "chimerarealigner.h"
  12 #include "kmerdb.hpp"
  13 #include "blastdb.hpp"
  14
  15 //***************************************************************************************************************
  16 ChimeraSlayer::ChimeraSlayer(string file, string temp, string mode, int k, int ms, int mms, int win, float div,
  17 int minsim, int mincov, int minbs, int minsnp, int par, int it, int inc, int numw, bool r) : Chimera()  {
  18         try {
  19                 fastafile = file;
  20                 templateFileName = temp; templateSeqs = readSeqs(temp);
  21                 searchMethod = mode;
  22                 kmerSize = k;
  23                 match = ms;
  24                 misMatch = mms;
  25                 window = win;
  26                 divR = div;
  27                 minSim = minsim;
  28                 minCov = mincov;
  29                 minBS = minbs;
  30                 minSNP = minsnp;
  31                 parents = par;
  32                 iters = it;
  33                 increment = inc;
  34                 numWanted = numw;
  35                 realign = r;
  36
  37                 decalc = new DeCalculator();
  38
  39                 doPrep();
  40         }
  41         catch(exception& e) {
  42                 m->errorOut(e, "ChimeraSlayer", "ChimeraSlayer");
  43                 exit(1);
  44         }
  45 }
  46 //***************************************************************************************************************
  47 ChimeraSlayer::ChimeraSlayer(string file, string temp, string name, string mode, string abunds, int k, int ms, int mms, int win, float div,
  48                                                          int minsim, int mincov, int minbs, int minsnp, int par, int it, int inc, int numw, bool r) : Chimera()  {
  49         try {
  50                 fastafile = file; templateSeqs = readSeqs(fastafile);
  51                 templateFileName = temp;
  52                 searchMethod = mode;
  53                 kmerSize = k;
  54                 match = ms;
  55                 misMatch = mms;
  56                 window = win;
  57                 divR = div;
  58                 minSim = minsim;
  59                 minCov = mincov;
  60                 minBS = minbs;
  61                 minSNP = minsnp;
  62                 parents = par;
  63                 iters = it;
  64                 increment = inc;
  65                 numWanted = numw;
  66                 realign = r;
  67                 includeAbunds = abunds;
  68
  69                 //read name file and create nameMapRank
  70                 readNameFile(name);
  71
  72                 decalc = new DeCalculator();
  73
  74                 createFilter(templateSeqs, 0.0); //just removed columns where all seqs have a gap
  75
  76                 //run filter on template
  77                 for (int i = 0; i < templateSeqs.size(); i++) { runFilter(templateSeqs[i]);  }
  78
  79         }
  80         catch(exception& e) {
  81                 m->errorOut(e, "ChimeraSlayer", "ChimeraSlayer");
  82                 exit(1);
  83         }
  84 }
  85 //***************************************************************************************************************
  86 int ChimeraSlayer::readNameFile(string name) {
  87         try {
  88                 ifstream in;
  89                 m->openInputFile(name, in);
  90
  91                 int maxRank = 0;
  92                 int minRank = 10000000;
  93
  94                 while(!in.eof()){
  95
  96                         if (m->control_pressed) { in.close(); return 0; }
  97
  98                         string thisname, repnames;
  99
 100                         in >> thisname;         m->gobble(in);          //read from first column
 101                         in >> repnames;                 //read from second column
 102
 103                         map<string, vector<string> >::iterator it = nameMapRank.find(thisname);
 104                         if (it == nameMapRank.end()) {
 105
 106                                 vector<string> splitRepNames;
 107                                 m->splitAtComma(repnames, splitRepNames);
 108
 109                                 nameMapRank[thisname] = splitRepNames;
 110
 111                                 if (splitRepNames.size() > maxRank) { maxRank = splitRepNames.size(); }
 112                                 if (splitRepNames.size() < minRank) { minRank = splitRepNames.size(); }
 113
 114                         }else{  m->mothurOut(thisname + " is already in namesfile. I will use first definition."); m->mothurOutEndLine();  }
 115
 116                         m->gobble(in);
 117                 }
 118                 in.close();
 119
 120                 //sanity check to make sure files match
 121                 for (int i = 0; i < templateSeqs.size(); i++) {
 122                         map<string, vector<string> >::iterator it = nameMapRank.find(templateSeqs[i]->getName());
 123
 124                         if (it == nameMapRank.end()) { m->mothurOut("[ERROR]: " + templateSeqs[i]->getName() + " is not in namesfile, but is in fastafile. Every name in fasta file must be in first column of names file."); m->mothurOutEndLine(); m->control_pressed = true;  }
 125                 }
 126
 127                 if (maxRank == minRank) { m->mothurOut("[ERROR]: all sequences in namesfile have the same abundance, aborting."); m->mothurOutEndLine(); m->control_pressed = true;  }
 128
 129                 return 0;
 130
 131         }
 132         catch(exception& e) {
 133                 m->errorOut(e, "ChimeraSlayer", "readNameFile");
 134                 exit(1);
 135         }
 136 }
 137
 138 //***************************************************************************************************************
 139 int ChimeraSlayer::doPrep() {
 140         try {
 141
 142                 //read in all query seqs
 143                 vector<Sequence*> tempQuerySeqs = readSeqs(fastafile);
 144
 145                 vector<Sequence*> temp = templateSeqs;
 146                 for (int i = 0; i < tempQuerySeqs.size(); i++) {  temp.push_back(tempQuerySeqs[i]);  }
 147
 148                 createFilter(temp, 0.0); //just removed columns where all seqs have a gap
 149
 150                 for (int i = 0; i < tempQuerySeqs.size(); i++) { delete tempQuerySeqs[i];  }
 151
 152                 if (m->control_pressed) {  return 0; }
 153
 154                 //run filter on template
 155                 for (int i = 0; i < templateSeqs.size(); i++) {  if (m->control_pressed) {  return 0; }  runFilter(templateSeqs[i]);  }
 156
 157                 string  kmerDBNameLeft;
 158                 string  kmerDBNameRight;
 159
 160                 //generate the kmerdb to pass to maligner
 161                 if (searchMethod == "kmer") {
 162                         string templatePath = m->hasPath(templateFileName);
 163                         string rightTemplateFileName = templatePath + "right." + m->getRootName(m->getSimpleName(templateFileName));
 164                         databaseRight = new KmerDB(rightTemplateFileName, kmerSize);
 165
 166                         string leftTemplateFileName = templatePath + "left." + m->getRootName(m->getSimpleName(templateFileName));
 167                         databaseLeft = new KmerDB(leftTemplateFileName, kmerSize);
 168                 #ifdef USE_MPI
 169                         for (int i = 0; i < templateSeqs.size(); i++) {
 170
 171                                 if (m->control_pressed) { return 0; }
 172
 173                                 string leftFrag = templateSeqs[i]->getUnaligned();
 174                                 leftFrag = leftFrag.substr(0, int(leftFrag.length() * 0.33));
 175
 176                                 Sequence leftTemp(templateSeqs[i]->getName(), leftFrag);
 177                                 databaseLeft->addSequence(leftTemp);
 178                         }
 179                         databaseLeft->generateDB();
 180                         databaseLeft->setNumSeqs(templateSeqs.size());
 181
 182                         for (int i = 0; i < templateSeqs.size(); i++) {
 183                                 if (m->control_pressed) { return 0; }
 184
 185                                 string rightFrag = templateSeqs[i]->getUnaligned();
 186                                 rightFrag = rightFrag.substr(int(rightFrag.length() * 0.66));
 187
 188                                 Sequence rightTemp(templateSeqs[i]->getName(), rightFrag);
 189                                 databaseRight->addSequence(rightTemp);
 190                         }
 191                         databaseRight->generateDB();
 192                         databaseRight->setNumSeqs(templateSeqs.size());
 193
 194                 #else
 195                         //leftside
 196                         kmerDBNameLeft = leftTemplateFileName.substr(0,leftTemplateFileName.find_last_of(".")+1) + char('0'+ kmerSize) + "mer";
 197                         ifstream kmerFileTestLeft(kmerDBNameLeft.c_str());
 198                         bool needToGenerateLeft = true;
 199
 200                         if(kmerFileTestLeft){
 201                                 bool GoodFile = m->checkReleaseVersion(kmerFileTestLeft, m->getVersion());
 202                                 if (GoodFile) {  needToGenerateLeft = false;    }
 203                         }
 204
 205                         if(needToGenerateLeft){
 206
 207                                 for (int i = 0; i < templateSeqs.size(); i++) {
 208
 209                                         if (m->control_pressed) { return 0; }
 210
 211                                         string leftFrag = templateSeqs[i]->getUnaligned();
 212                                         leftFrag = leftFrag.substr(0, int(leftFrag.length() * 0.33));
 213
 214                                         Sequence leftTemp(templateSeqs[i]->getName(), leftFrag);
 215                                         databaseLeft->addSequence(leftTemp);
 216                                 }
 217                                 databaseLeft->generateDB();
 218
 219                         }else {
 220                                 databaseLeft->readKmerDB(kmerFileTestLeft);
 221                         }
 222                         kmerFileTestLeft.close();
 223
 224                         databaseLeft->setNumSeqs(templateSeqs.size());
 225
 226                         //rightside
 227                         kmerDBNameRight = rightTemplateFileName.substr(0,rightTemplateFileName.find_last_of(".")+1) + char('0'+ kmerSize) + "mer";
 228                         ifstream kmerFileTestRight(kmerDBNameRight.c_str());
 229                         bool needToGenerateRight = true;
 230
 231                         if(kmerFileTestRight){
 232                                 bool GoodFile = m->checkReleaseVersion(kmerFileTestRight, m->getVersion());
 233                                 if (GoodFile) {  needToGenerateRight = false;   }
 234                         }
 235
 236                         if(needToGenerateRight){
 237
 238                                 for (int i = 0; i < templateSeqs.size(); i++) {
 239                                         if (m->control_pressed) { return 0; }
 240
 241                                         string rightFrag = templateSeqs[i]->getUnaligned();
 242                                         rightFrag = rightFrag.substr(int(rightFrag.length() * 0.66));
 243
 244                                         Sequence rightTemp(templateSeqs[i]->getName(), rightFrag);
 245                                         databaseRight->addSequence(rightTemp);
 246                                 }
 247                                 databaseRight->generateDB();
 248
 249                         }else {
 250                                 databaseRight->readKmerDB(kmerFileTestRight);
 251                         }
 252                         kmerFileTestRight.close();
 253
 254                         databaseRight->setNumSeqs(templateSeqs.size());
 255                 #endif
 256                 }else if (searchMethod == "blast") {
 257
 258                         //generate blastdb
 259                         databaseLeft = new BlastDB(-2.0, -1.0, match, misMatch);
 260                         for (int i = 0; i < templateSeqs.size(); i++) {         databaseLeft->addSequence(*templateSeqs[i]);    }
 261                         databaseLeft->generateDB();
 262                         databaseLeft->setNumSeqs(templateSeqs.size());
 263                 }
 264
 265                 return 0;
 266
 267         }
 268         catch(exception& e) {
 269                 m->errorOut(e, "ChimeraSlayer", "doprep");
 270                 exit(1);
 271         }
 272 }
 273 //***************************************************************************************************************
 274 vector<Sequence*> ChimeraSlayer::getTemplate(Sequence* q) {
 275         try {
 276
 277                 vector<Sequence*> thisTemplate;
 278
 279                 int thisRank;
 280                 string thisName = q->getName();
 281                 map<string, vector<string> >::iterator itRank = nameMapRank.find(thisName); // you will find it because we already sanity checked
 282                 thisRank = (itRank->second).size();
 283
 284                 //create list of names we want to put into the template
 285                 set<string> namesToAdd;
 286                 for (itRank = nameMapRank.begin(); itRank != nameMapRank.end(); itRank++) {
 287                         if (itRank->first != thisName) {
 288                                 if (includeAbunds == "greaterequal") {
 289                                         if ((itRank->second).size() >= thisRank) {
 290                                                 //you are more abundant than me or equal to my abundance
 291                                                 for (int i = 0; i < (itRank->second).size(); i++) {
 292                                                         namesToAdd.insert((itRank->second)[i]);
 293                                                 }
 294                                         }
 295                                 }else if (includeAbunds == "greater") {
 296                                         if ((itRank->second).size() > thisRank) {
 297                                                 //you are more abundant than me
 298                                                 for (int i = 0; i < (itRank->second).size(); i++) {
 299                                                         namesToAdd.insert((itRank->second)[i]);
 300                                                 }
 301                                         }
 302                                 }else if (includeAbunds == "all") {
 303                                         //add everyone
 304                                         for (int i = 0; i < (itRank->second).size(); i++) {
 305                                                 namesToAdd.insert((itRank->second)[i]);
 306                                         }
 307                                 }
 308                         }
 309                 }
 310
 311                 for (int i = 0; i < templateSeqs.size(); i++) {
 312                         if (namesToAdd.count(templateSeqs[i]->getName()) != 0) {
 313                                 thisTemplate.push_back(templateSeqs[i]);
 314                         }
 315                 }
 316
 317                 string  kmerDBNameLeft;
 318                 string  kmerDBNameRight;
 319
 320                 //generate the kmerdb to pass to maligner
 321                 if (searchMethod == "kmer") {
 322                         string templatePath = m->hasPath(templateFileName);
 323                         string rightTemplateFileName = templatePath + "right." + m->getRootName(m->getSimpleName(templateFileName));
 324                         databaseRight = new KmerDB(rightTemplateFileName, kmerSize);
 325
 326                         string leftTemplateFileName = templatePath + "left." + m->getRootName(m->getSimpleName(templateFileName));
 327                         databaseLeft = new KmerDB(leftTemplateFileName, kmerSize);
 328 #ifdef USE_MPI
 329                         for (int i = 0; i < thisTemplate.size(); i++) {
 330
 331                                 if (m->control_pressed) { return thisTemplate; }
 332
 333                                 string leftFrag = thisTemplate[i]->getUnaligned();
 334                                 leftFrag = leftFrag.substr(0, int(leftFrag.length() * 0.33));
 335
 336                                 Sequence leftTemp(thisTemplate[i]->getName(), leftFrag);
 337                                 databaseLeft->addSequence(leftTemp);
 338                         }
 339                         databaseLeft->generateDB();
 340                         databaseLeft->setNumSeqs(thisTemplate.size());
 341
 342                         for (int i = 0; i < thisTemplate.size(); i++) {
 343                                 if (m->control_pressed) { return thisTemplate;  }
 344
 345                                 string rightFrag = thisTemplate[i]->getUnaligned();
 346                                 rightFrag = rightFrag.substr(int(rightFrag.length() * 0.66));
 347
 348                                 Sequence rightTemp(thisTemplate[i]->getName(), rightFrag);
 349                                 databaseRight->addSequence(rightTemp);
 350                         }
 351                         databaseRight->generateDB();
 352                         databaseRight->setNumSeqs(thisTemplate.size());
 353
 354 #else
 355
 356
 357                         for (int i = 0; i < thisTemplate.size(); i++) {
 358
 359                                 if (m->control_pressed) { return thisTemplate; }
 360
 361                                 string leftFrag = thisTemplate[i]->getUnaligned();
 362                                 leftFrag = leftFrag.substr(0, int(leftFrag.length() * 0.33));
 363
 364                                 Sequence leftTemp(thisTemplate[i]->getName(), leftFrag);
 365                                 databaseLeft->addSequence(leftTemp);
 366                         }
 367                         databaseLeft->generateDB();
 368                         databaseLeft->setNumSeqs(thisTemplate.size());
 369
 370                         for (int i = 0; i < thisTemplate.size(); i++) {
 371                                 if (m->control_pressed) { return thisTemplate; }
 372
 373                                 string rightFrag = thisTemplate[i]->getUnaligned();
 374                                 rightFrag = rightFrag.substr(int(rightFrag.length() * 0.66));
 375
 376                                 Sequence rightTemp(thisTemplate[i]->getName(), rightFrag);
 377                                 databaseRight->addSequence(rightTemp);
 378                         }
 379                         databaseRight->generateDB();
 380                         databaseRight->setNumSeqs(thisTemplate.size());
 381 #endif
 382                 }else if (searchMethod == "blast") {
 383
 384                         //generate blastdb
 385                         databaseLeft = new BlastDB(-2.0, -1.0, match, misMatch);
 386                         for (int i = 0; i < thisTemplate.size(); i++) { if (m->control_pressed) { return thisTemplate; }  databaseLeft->addSequence(*thisTemplate[i]);  }
 387                         databaseLeft->generateDB();
 388                         databaseLeft->setNumSeqs(thisTemplate.size());
 389                 }
 390
 391                 return thisTemplate;
 392
 393         }
 394         catch(exception& e) {
 395                 m->errorOut(e, "ChimeraSlayer", "getTemplate");
 396                 exit(1);
 397         }
 398 }
 399
 400 //***************************************************************************************************************
 401 ChimeraSlayer::~ChimeraSlayer() {
 402         delete decalc;
 403         if (templateFileName != "self") {
 404                 if (searchMethod == "kmer") {  delete databaseRight;  delete databaseLeft;  }
 405                 else if (searchMethod == "blast") {  delete databaseLeft; }
 406         }
 407 }
 408 //***************************************************************************************************************
 409 void ChimeraSlayer::printHeader(ostream& out) {
 410         m->mothurOutEndLine();
 411         m->mothurOut("Only reporting sequence supported by " + toString(minBS) + "% of bootstrapped results.");
 412         m->mothurOutEndLine();
 413
 414         out << "Name\tLeftParent\tRightParent\tDivQLAQRB\tPerIDQLAQRB\tBootStrapA\tDivQLBQRA\tPerIDQLBQRA\tBootStrapB\tFlag\tLeftWindow\tRightWindow\n";
 415 }
 416 //***************************************************************************************************************
 417 int ChimeraSlayer::print(ostream& out, ostream& outAcc) {
 418         try {
 419                 if (chimeraFlags == "yes") {
 420                         string chimeraFlag = "no";
 421                         if(  (chimeraResults[0].bsa >= minBS && chimeraResults[0].divr_qla_qrb >= divR)
 422                            ||
 423                            (chimeraResults[0].bsb >= minBS && chimeraResults[0].divr_qlb_qra >= divR) ) { chimeraFlag = "yes"; }
 424
 425
 426                         if (chimeraFlag == "yes") {
 427                                 if ((chimeraResults[0].bsa >= minBS) || (chimeraResults[0].bsb >= minBS)) {
 428                                         m->mothurOut(querySeq->getName() + "\tyes"); m->mothurOutEndLine();
 429                                         outAcc << querySeq->getName() << endl;
 430                                 }
 431                         }
 432
 433                         printBlock(chimeraResults[0], chimeraFlag, out);
 434                         out << endl;
 435                 }else {  out << querySeq->getName() << "\tno" << endl;  }
 436
 437                 return 0;
 438
 439         }
 440         catch(exception& e) {
 441                 m->errorOut(e, "ChimeraSlayer", "print");
 442                 exit(1);
 443         }
 444 }
 445 #ifdef USE_MPI
 446 //***************************************************************************************************************
 447 int ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc) {
 448         try {
 449                 MPI_Status status;
 450                 bool results = false;
 451                 string outAccString = "";
 452                 string outputString = "";
 453
 454                 if (chimeraFlags == "yes") {
 455                         string chimeraFlag = "no";
 456                         if(  (chimeraResults[0].bsa >= minBS && chimeraResults[0].divr_qla_qrb >= divR)
 457                            ||
 458                            (chimeraResults[0].bsb >= minBS && chimeraResults[0].divr_qlb_qra >= divR) ) { chimeraFlag = "yes"; }
 459
 460
 461                         if (chimeraFlag == "yes") {
 462                                 if ((chimeraResults[0].bsa >= minBS) || (chimeraResults[0].bsb >= minBS)) {
 463                                         cout << querySeq->getName() <<  "\tyes" << endl;
 464                                         outAccString += querySeq->getName() + "\n";
 465                                         results = true;
 466
 467                                         //write to accnos file
 468                                         int length = outAccString.length();
 469                                         char* buf2 = new char[length];
 470                                         memcpy(buf2, outAccString.c_str(), length);
 471
 472                                         MPI_File_write_shared(outAcc, buf2, length, MPI_CHAR, &status);
 473                                         delete buf2;
 474                                 }
 475                         }
 476
 477                         outputString = getBlock(chimeraResults[0], chimeraFlag);
 478                         outputString += "\n";
 479         //cout << outputString << endl;
 480                         //write to output file
 481                         int length = outputString.length();
 482                         char* buf = new char[length];
 483                         memcpy(buf, outputString.c_str(), length);
 484
 485                         MPI_File_write_shared(out, buf, length, MPI_CHAR, &status);
 486                         delete buf;
 487
 488                 }else {
 489                         outputString += querySeq->getName() + "\tno\n";
 490         //cout << outputString << endl;
 491                         //write to output file
 492                         int length = outputString.length();
 493                         char* buf = new char[length];
 494                         memcpy(buf, outputString.c_str(), length);
 495
 496                         MPI_File_write_shared(out, buf, length, MPI_CHAR, &status);
 497                         delete buf;
 498                 }
 499
 500
 501                 return results;
 502         }
 503         catch(exception& e) {
 504                 m->errorOut(e, "ChimeraSlayer", "print");
 505                 exit(1);
 506         }
 507 }
 508 #endif
 509
 510 //***************************************************************************************************************
 511 int ChimeraSlayer::getChimeras(Sequence* query) {
 512         try {
 513                 chimeraFlags = "no";
 514
 515                 //filter query
 516                 spotMap = runFilter(query);
 517
 518                 querySeq = query;
 519
 520                 //you must create a template
 521                 vector<Sequence*> thisTemplate;
 522                 if (templateFileName != "self") { thisTemplate = templateSeqs; }
 523                 else { thisTemplate = getTemplate(query); } //fills thistemplate and creates the databases
 524
 525                 if (m->control_pressed) {  return 0;  }
 526
 527                 if (thisTemplate.size() == 0) {  return 0; } //not chimeric
 528
 529                 //referenceSeqs, numWanted, matchScore, misMatchPenalty, divR, minSimilarity
 530                 Maligner maligner(thisTemplate, numWanted, match, misMatch, divR, minSim, minCov, searchMethod, databaseLeft, databaseRight);
 531                 Slayer slayer(window, increment, minSim, divR, iters, minSNP);
 532
 533                 if (templateFileName == "self") {
 534                         if (searchMethod == "kmer") {  delete databaseRight;  delete databaseLeft;  }
 535                         else if (searchMethod == "blast") {  delete databaseLeft; }
 536                 }
 537
 538                 if (m->control_pressed) {  return 0;  }
 539
 540                 string chimeraFlag = maligner.getResults(query, decalc);
 541                 if (m->control_pressed) {  return 0;  }
 542                 vector<results> Results = maligner.getOutput();
 543
 544                 //found in testing realigning only made things worse
 545                 if (realign) {
 546                         ChimeraReAligner realigner(thisTemplate, match, misMatch);
 547                         realigner.reAlign(query, Results);
 548                 }
 549
 550                 if (chimeraFlag == "yes") {
 551
 552                         //get sequence that were given from maligner results
 553                         vector<SeqDist> seqs;
 554                         map<string, float> removeDups;
 555                         map<string, float>::iterator itDup;
 556                         map<string, string> parentNameSeq;
 557                         map<string, string>::iterator itSeq;
 558                         for (int j = 0; j < Results.size(); j++) {
 559                                 float dist = (Results[j].regionEnd - Results[j].regionStart + 1) * Results[j].queryToParentLocal;
 560                                 //only add if you are not a duplicate
 561                                 itDup = removeDups.find(Results[j].parent);
 562                                 if (itDup == removeDups.end()) { //this is not duplicate
 563                                         removeDups[Results[j].parent] = dist;
 564                                         parentNameSeq[Results[j].parent] = Results[j].parentAligned;
 565                                 }else if (dist > itDup->second) { //is this a stronger number for this parent
 566                                         removeDups[Results[j].parent] = dist;
 567                                         parentNameSeq[Results[j].parent] = Results[j].parentAligned;
 568                                 }
 569                         }
 570
 571                         for (itDup = removeDups.begin(); itDup != removeDups.end(); itDup++) {
 572                                 //Sequence* seq = getSequence(itDup->first); //makes copy so you can filter and mask and not effect template
 573                                 itSeq = parentNameSeq.find(itDup->first);
 574 //cout << itDup->first << itSeq->second << endl;
 575                                 Sequence* seq = new Sequence(itDup->first, itSeq->second);
 576
 577                                 SeqDist member;
 578                                 member.seq = seq;
 579                                 member.dist = itDup->second;
 580
 581                                 seqs.push_back(member);
 582                         }
 583
 584                         //limit number of parents to explore - default 3
 585                         if (Results.size() > parents) {
 586                                 //sort by distance
 587                                 sort(seqs.begin(), seqs.end(), compareSeqDist);
 588                                 //prioritize larger more similiar sequence fragments
 589                                 reverse(seqs.begin(), seqs.end());
 590
 591                                 for (int k = seqs.size()-1; k > (parents-1); k--)  {
 592                                         delete seqs[k].seq;
 593                                         seqs.pop_back();
 594                                 }
 595                         }
 596
 597                         //put seqs into vector to send to slayer
 598                         vector<Sequence*> seqsForSlayer;
 599                         for (int k = 0; k < seqs.size(); k++) {  seqsForSlayer.push_back(seqs[k].seq);  }
 600
 601                         //mask then send to slayer...
 602                         if (seqMask != "") {
 603                                 decalc->setMask(seqMask);
 604
 605                                 //mask querys
 606                                 decalc->runMask(query);
 607
 608                                 //mask parents
 609                                 for (int k = 0; k < seqsForSlayer.size(); k++) {
 610                                         decalc->runMask(seqsForSlayer[k]);
 611                                 }
 612
 613                                 spotMap = decalc->getMaskMap();
 614                         }
 615
 616                         if (m->control_pressed) {  for (int k = 0; k < seqs.size(); k++) {  delete seqs[k].seq;   }  return 0;  }
 617
 618                         //send to slayer
 619                         chimeraFlags = slayer.getResults(query, seqsForSlayer);
 620                         if (m->control_pressed) {  return 0;  }
 621                         chimeraResults = slayer.getOutput();
 622
 623                         //free memory
 624                         for (int k = 0; k < seqs.size(); k++) {  delete seqs[k].seq;   }
 625                 }
 626
 627                 //delete maligner;
 628                 //delete slayer;
 629
 630                 return 0;
 631         }
 632         catch(exception& e) {
 633                 m->errorOut(e, "ChimeraSlayer", "getChimeras");
 634                 exit(1);
 635         }
 636 }
 637 //***************************************************************************************************************
 638 void ChimeraSlayer::printBlock(data_struct data, string flag, ostream& out){
 639         try {
 640         //out << ":)\n";
 641
 642                 out << querySeq->getName() << '\t';
 643                 out << data.parentA.getName() << "\t" << data.parentB.getName()  << '\t';
 644                 //out << "Left Window: " << spotMap[data.winLStart] << " " << spotMap[data.winLEnd] << endl;
 645                 //out << "Right Window: " << spotMap[data.winRStart] << " " << spotMap[data.winREnd] << endl;
 646
 647                 out << data.divr_qla_qrb << '\t' << data.qla_qrb << '\t' << data.bsa << '\t';
 648                 out << data.divr_qlb_qra << '\t' << data.qlb_qra << '\t' << data.bsb << '\t';
 649
 650                 out << flag << '\t' << spotMap[data.winLStart] << "-" << spotMap[data.winLEnd] << '\t' << spotMap[data.winRStart] << "-" << spotMap[data.winREnd] << '\t';
 651
 652                 //out << "Similarity of parents: " << data.ab << endl;
 653                 //out << "Similarity of query to parentA: " << data.qa << endl;
 654                 //out << "Similarity of query to parentB: " << data.qb << endl;
 655
 656
 657                 //out << "Per_id(QL,A): " << data.qla << endl;
 658                 //out << "Per_id(QL,B): " << data.qlb << endl;
 659                 //out << "Per_id(QR,A): " << data.qra << endl;
 660                 //out << "Per_id(QR,B): " << data.qrb << endl;
 661
 662
 663                 //out << "DeltaL: " << (data.qla - data.qlb) << endl;
 664                 //out << "DeltaR: " << (data.qra - data.qrb) << endl;
 665
 666         }
 667         catch(exception& e) {
 668                 m->errorOut(e, "ChimeraSlayer", "printBlock");
 669                 exit(1);
 670         }
 671 }
 672 //***************************************************************************************************************
 673 string ChimeraSlayer::getBlock(data_struct data, string flag){
 674         try {
 675
 676                 string outputString = "";
 677
 678                 outputString += querySeq->getName() + "\t";
 679                 outputString += data.parentA.getName() + "\t" + data.parentB.getName()  + "\t";
 680
 681                 outputString += toString(data.divr_qla_qrb) + "\t" + toString(data.qla_qrb) + "\t" + toString(data.bsa) + "\t";
 682                 outputString += toString(data.divr_qlb_qra) + "\t" + toString(data.qlb_qra) + "\t" + toString(data.bsb) + "\t";
 683
 684                 outputString += flag + "\t" + toString(spotMap[data.winLStart]) + "-" + toString(spotMap[data.winLEnd]) + "\t" + toString(spotMap[data.winRStart]) + "-" + toString(spotMap[data.winREnd]) + "\t";
 685
 686                 return outputString;
 687         }
 688         catch(exception& e) {
 689                 m->errorOut(e, "ChimeraSlayer", "getBlock");
 690                 exit(1);
 691         }
 692 }
 693 //***************************************************************************************************************/
 694