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, bool trim, 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                 trimChimera = trim;
  37
  38                 decalc = new DeCalculator();
  39
  40                 doPrep();
  41         }
  42         catch(exception& e) {
  43                 m->errorOut(e, "ChimeraSlayer", "ChimeraSlayer");
  44                 exit(1);
  45         }
  46 }
  47 //***************************************************************************************************************
  48 ChimeraSlayer::ChimeraSlayer(string file, string temp, bool trim, map<string, int>& prior, string mode, int k, int ms, int mms, int win, float div,
  49                                                          int minsim, int mincov, int minbs, int minsnp, int par, int it, int inc, int numw, bool r) : Chimera()  {
  50         try {
  51                 fastafile = file; templateSeqs = readSeqs(fastafile);
  52                 templateFileName = temp;
  53                 searchMethod = mode;
  54                 kmerSize = k;
  55                 match = ms;
  56                 misMatch = mms;
  57                 window = win;
  58                 divR = div;
  59                 minSim = minsim;
  60                 minCov = mincov;
  61                 minBS = minbs;
  62                 minSNP = minsnp;
  63                 parents = par;
  64                 iters = it;
  65                 increment = inc;
  66                 numWanted = numw;
  67                 realign = r;
  68                 trimChimera = trim;
  69                 priority = prior;
  70
  71                 decalc = new DeCalculator();
  72
  73                 createFilter(templateSeqs, 0.0); //just removed columns where all seqs have a gap
  74
  75                 //run filter on template
  76                 for (int i = 0; i < templateSeqs.size(); i++) {  if (m->control_pressed) {  break; }  runFilter(templateSeqs[i]);  }
  77
  78
  79         }
  80         catch(exception& e) {
  81                 m->errorOut(e, "ChimeraSlayer", "ChimeraSlayer");
  82                 exit(1);
  83         }
  84 }
  85 //***************************************************************************************************************
  86 int ChimeraSlayer::doPrep() {
  87         try {
  88                 if (searchMethod == "distance") {
  89                         //read in all query seqs
  90                         vector<Sequence*> tempQuerySeqs = readSeqs(fastafile);
  91
  92                         vector<Sequence*> temp = templateSeqs;
  93                         for (int i = 0; i < tempQuerySeqs.size(); i++) {  temp.push_back(tempQuerySeqs[i]);  }
  94
  95                         createFilter(temp, 0.0); //just removed columns where all seqs have a gap
  96
  97                         for (int i = 0; i < tempQuerySeqs.size(); i++) { delete tempQuerySeqs[i];  }
  98
  99                         if (m->control_pressed) {  return 0; }
 100
 101                         //run filter on template
 102                         for (int i = 0; i < templateSeqs.size(); i++) {  if (m->control_pressed) {  return 0; }  runFilter(templateSeqs[i]);  }
 103                 }
 104                 string  kmerDBNameLeft;
 105                 string  kmerDBNameRight;
 106
 107                 //generate the kmerdb to pass to maligner
 108                 if (searchMethod == "kmer") {
 109                         string templatePath = m->hasPath(templateFileName);
 110                         string rightTemplateFileName = templatePath + "right." + m->getRootName(m->getSimpleName(templateFileName));
 111                         databaseRight = new KmerDB(rightTemplateFileName, kmerSize);
 112
 113                         string leftTemplateFileName = templatePath + "left." + m->getRootName(m->getSimpleName(templateFileName));
 114                         databaseLeft = new KmerDB(leftTemplateFileName, kmerSize);
 115                 #ifdef USE_MPI
 116                         for (int i = 0; i < templateSeqs.size(); i++) {
 117
 118                                 if (m->control_pressed) { return 0; }
 119
 120                                 string leftFrag = templateSeqs[i]->getUnaligned();
 121                                 leftFrag = leftFrag.substr(0, int(leftFrag.length() * 0.33));
 122
 123                                 Sequence leftTemp(templateSeqs[i]->getName(), leftFrag);
 124                                 databaseLeft->addSequence(leftTemp);
 125                         }
 126                         databaseLeft->generateDB();
 127                         databaseLeft->setNumSeqs(templateSeqs.size());
 128
 129                         for (int i = 0; i < templateSeqs.size(); i++) {
 130                                 if (m->control_pressed) { return 0; }
 131
 132                                 string rightFrag = templateSeqs[i]->getUnaligned();
 133                                 rightFrag = rightFrag.substr(int(rightFrag.length() * 0.66));
 134
 135                                 Sequence rightTemp(templateSeqs[i]->getName(), rightFrag);
 136                                 databaseRight->addSequence(rightTemp);
 137                         }
 138                         databaseRight->generateDB();
 139                         databaseRight->setNumSeqs(templateSeqs.size());
 140
 141                 #else
 142                         //leftside
 143                         kmerDBNameLeft = leftTemplateFileName.substr(0,leftTemplateFileName.find_last_of(".")+1) + char('0'+ kmerSize) + "mer";
 144                         ifstream kmerFileTestLeft(kmerDBNameLeft.c_str());
 145                         bool needToGenerateLeft = true;
 146
 147                         if(kmerFileTestLeft){
 148                                 bool GoodFile = m->checkReleaseVersion(kmerFileTestLeft, m->getVersion());
 149                                 if (GoodFile) {  needToGenerateLeft = false;    }
 150                         }
 151
 152                         if(needToGenerateLeft){
 153
 154                                 for (int i = 0; i < templateSeqs.size(); i++) {
 155
 156                                         if (m->control_pressed) { return 0; }
 157
 158                                         string leftFrag = templateSeqs[i]->getUnaligned();
 159                                         leftFrag = leftFrag.substr(0, int(leftFrag.length() * 0.33));
 160
 161                                         Sequence leftTemp(templateSeqs[i]->getName(), leftFrag);
 162                                         databaseLeft->addSequence(leftTemp);
 163                                 }
 164                                 databaseLeft->generateDB();
 165
 166                         }else {
 167                                 databaseLeft->readKmerDB(kmerFileTestLeft);
 168                         }
 169                         kmerFileTestLeft.close();
 170
 171                         databaseLeft->setNumSeqs(templateSeqs.size());
 172
 173                         //rightside
 174                         kmerDBNameRight = rightTemplateFileName.substr(0,rightTemplateFileName.find_last_of(".")+1) + char('0'+ kmerSize) + "mer";
 175                         ifstream kmerFileTestRight(kmerDBNameRight.c_str());
 176                         bool needToGenerateRight = true;
 177
 178                         if(kmerFileTestRight){
 179                                 bool GoodFile = m->checkReleaseVersion(kmerFileTestRight, m->getVersion());
 180                                 if (GoodFile) {  needToGenerateRight = false;   }
 181                         }
 182
 183                         if(needToGenerateRight){
 184
 185                                 for (int i = 0; i < templateSeqs.size(); i++) {
 186                                         if (m->control_pressed) { return 0; }
 187
 188                                         string rightFrag = templateSeqs[i]->getUnaligned();
 189                                         rightFrag = rightFrag.substr(int(rightFrag.length() * 0.66));
 190
 191                                         Sequence rightTemp(templateSeqs[i]->getName(), rightFrag);
 192                                         databaseRight->addSequence(rightTemp);
 193                                 }
 194                                 databaseRight->generateDB();
 195
 196                         }else {
 197                                 databaseRight->readKmerDB(kmerFileTestRight);
 198                         }
 199                         kmerFileTestRight.close();
 200
 201                         databaseRight->setNumSeqs(templateSeqs.size());
 202                 #endif
 203                 }else if (searchMethod == "blast") {
 204
 205                         //generate blastdb
 206                         databaseLeft = new BlastDB(-1.0, -1.0, 1, -3);
 207
 208                         for (int i = 0; i < templateSeqs.size(); i++) {         databaseLeft->addSequence(*templateSeqs[i]);    }
 209                         databaseLeft->generateDB();
 210                         databaseLeft->setNumSeqs(templateSeqs.size());
 211                 }
 212
 213                 return 0;
 214
 215         }
 216         catch(exception& e) {
 217                 m->errorOut(e, "ChimeraSlayer", "doprep");
 218                 exit(1);
 219         }
 220 }
 221 //***************************************************************************************************************
 222 vector<Sequence*> ChimeraSlayer::getTemplate(Sequence* q) {
 223         try {
 224
 225                 //when template=self, the query file is sorted from most abundance to least abundant
 226                 //userTemplate grows as the query file is processed by adding sequences that are not chimeric and more abundant
 227                 vector<Sequence*> userTemplate;
 228
 229                 int myAbund = priority[q->getName()];
 230
 231                 for (int i = 0; i < templateSeqs.size(); i++) {
 232
 233                         if (m->control_pressed) { return userTemplate; }
 234
 235                         //have I reached a sequence with the same abundance as myself?
 236                         if (!(priority[templateSeqs[i]->getName()] > myAbund)) { break; }
 237
 238                         //if its am not chimeric add it
 239                         if (chimericSeqs.count(templateSeqs[i]->getName()) == 0) { userTemplate.push_back(templateSeqs[i]); }
 240                 }
 241
 242                 string  kmerDBNameLeft;
 243                 string  kmerDBNameRight;
 244
 245                 //generate the kmerdb to pass to maligner
 246                 if (searchMethod == "kmer") {
 247                         string templatePath = m->hasPath(templateFileName);
 248                         string rightTemplateFileName = templatePath + "right." + m->getRootName(m->getSimpleName(templateFileName));
 249                         databaseRight = new KmerDB(rightTemplateFileName, kmerSize);
 250
 251                         string leftTemplateFileName = templatePath + "left." + m->getRootName(m->getSimpleName(templateFileName));
 252                         databaseLeft = new KmerDB(leftTemplateFileName, kmerSize);
 253 #ifdef USE_MPI
 254                         for (int i = 0; i < userTemplate.size(); i++) {
 255
 256                                 if (m->control_pressed) { return userTemplate; }
 257
 258                                 string leftFrag = userTemplate[i]->getUnaligned();
 259                                 leftFrag = leftFrag.substr(0, int(leftFrag.length() * 0.33));
 260
 261                                 Sequence leftTemp(userTemplate[i]->getName(), leftFrag);
 262                                 databaseLeft->addSequence(leftTemp);
 263                         }
 264                         databaseLeft->generateDB();
 265                         databaseLeft->setNumSeqs(userTemplate.size());
 266
 267                         for (int i = 0; i < userTemplate.size(); i++) {
 268                                 if (m->control_pressed) { return userTemplate; }
 269
 270                                 string rightFrag = userTemplate[i]->getUnaligned();
 271                                 rightFrag = rightFrag.substr(int(rightFrag.length() * 0.66));
 272
 273                                 Sequence rightTemp(userTemplate[i]->getName(), rightFrag);
 274                                 databaseRight->addSequence(rightTemp);
 275                         }
 276                         databaseRight->generateDB();
 277                         databaseRight->setNumSeqs(userTemplate.size());
 278
 279 #else
 280
 281
 282                         for (int i = 0; i < userTemplate.size(); i++) {
 283
 284                                 if (m->control_pressed) { return userTemplate; }
 285
 286                                 string leftFrag = userTemplate[i]->getUnaligned();
 287                                 leftFrag = leftFrag.substr(0, int(leftFrag.length() * 0.33));
 288
 289                                 Sequence leftTemp(userTemplate[i]->getName(), leftFrag);
 290                                 databaseLeft->addSequence(leftTemp);
 291                         }
 292                         databaseLeft->generateDB();
 293                         databaseLeft->setNumSeqs(userTemplate.size());
 294
 295                         for (int i = 0; i < userTemplate.size(); i++) {
 296                                 if (m->control_pressed) { return userTemplate; }
 297
 298                                 string rightFrag = userTemplate[i]->getUnaligned();
 299                                 rightFrag = rightFrag.substr(int(rightFrag.length() * 0.66));
 300
 301                                 Sequence rightTemp(userTemplate[i]->getName(), rightFrag);
 302                                 databaseRight->addSequence(rightTemp);
 303                         }
 304                         databaseRight->generateDB();
 305                         databaseRight->setNumSeqs(userTemplate.size());
 306 #endif
 307                 }else if (searchMethod == "blast") {
 308
 309                         //generate blastdb
 310                         databaseLeft = new BlastDB(-1.0, -1.0, 1, -3);
 311
 312                         for (int i = 0; i < userTemplate.size(); i++) { if (m->control_pressed) { return userTemplate; }   databaseLeft->addSequence(*userTemplate[i]); }
 313                         databaseLeft->generateDB();
 314                         databaseLeft->setNumSeqs(userTemplate.size());
 315                 }
 316
 317                 return userTemplate;
 318
 319         }
 320         catch(exception& e) {
 321                 m->errorOut(e, "ChimeraSlayer", "getTemplate");
 322                 exit(1);
 323         }
 324 }
 325
 326 //***************************************************************************************************************
 327 ChimeraSlayer::~ChimeraSlayer() {
 328         delete decalc;
 329         if (templateFileName != "self") {
 330                 if (searchMethod == "kmer") {  delete databaseRight;  delete databaseLeft;  }
 331                 else if (searchMethod == "blast") {  delete databaseLeft; }
 332         }
 333 }
 334 //***************************************************************************************************************
 335 void ChimeraSlayer::printHeader(ostream& out) {
 336         m->mothurOutEndLine();
 337         m->mothurOut("Only reporting sequence supported by " + toString(minBS) + "% of bootstrapped results.");
 338         m->mothurOutEndLine();
 339
 340         out << "Name\tLeftParent\tRightParent\tDivQLAQRB\tPerIDQLAQRB\tBootStrapA\tDivQLBQRA\tPerIDQLBQRA\tBootStrapB\tFlag\tLeftWindow\tRightWindow\n";
 341 }
 342 //***************************************************************************************************************
 343 Sequence* ChimeraSlayer::print(ostream& out, ostream& outAcc) {
 344         try {
 345                 Sequence* trim = NULL;
 346                 if (trimChimera) { trim = new Sequence(trimQuery.getName(), trimQuery.getAligned()); }
 347
 348                 if (chimeraFlags == "yes") {
 349                         string chimeraFlag = "no";
 350                         if(  (chimeraResults[0].bsa >= minBS && chimeraResults[0].divr_qla_qrb >= divR)
 351                            ||
 352                            (chimeraResults[0].bsb >= minBS && chimeraResults[0].divr_qlb_qra >= divR) ) { chimeraFlag = "yes"; }
 353
 354
 355                         if (chimeraFlag == "yes") {
 356                                 if ((chimeraResults[0].bsa >= minBS) || (chimeraResults[0].bsb >= minBS)) {
 357                                         m->mothurOut(querySeq->getName() + "\tyes"); m->mothurOutEndLine();
 358                                         outAcc << querySeq->getName() << endl;
 359
 360                                         if (templateFileName == "self") {  chimericSeqs.insert(querySeq->getName()); }
 361
 362                                         if (trimChimera) {
 363                                                 int lengthLeft = spotMap[chimeraResults[0].winLEnd] - spotMap[chimeraResults[0].winLStart];
 364                                                 int lengthRight = spotMap[chimeraResults[0].winREnd] - spotMap[chimeraResults[0].winRStart];
 365
 366                                                 string newAligned = trim->getAligned();
 367
 368                                                 if (lengthLeft > lengthRight) { //trim right
 369                                                         for (int i = (spotMap[chimeraResults[0].winRStart]-1); i < newAligned.length(); i++) { newAligned[i] = '.'; }
 370                                                 }else { //trim left
 371                                                         for (int i = 0; i < spotMap[chimeraResults[0].winLEnd]; i++) { newAligned[i] = '.'; }
 372                                                 }
 373                                                 trim->setAligned(newAligned);
 374                                         }
 375                                 }
 376                         }
 377
 378                         printBlock(chimeraResults[0], chimeraFlag, out);
 379                         out << endl;
 380                 }else {
 381                         out << querySeq->getName() << "\tno" << endl;
 382                 }
 383
 384                 return trim;
 385
 386         }
 387         catch(exception& e) {
 388                 m->errorOut(e, "ChimeraSlayer", "print");
 389                 exit(1);
 390         }
 391 }
 392 //***************************************************************************************************************
 393 Sequence* ChimeraSlayer::print(ostream& out, ostream& outAcc, data_results leftPiece, data_results rightPiece) {
 394         try {
 395                 Sequence* trim = NULL;
 396
 397                 if (trimChimera) {
 398                         string aligned = leftPiece.trimQuery.getAligned() + rightPiece.trimQuery.getAligned();
 399                         trim = new Sequence(leftPiece.trimQuery.getName(), aligned);
 400                 }
 401
 402                 if ((leftPiece.flag == "yes") || (rightPiece.flag == "yes")) {
 403
 404                         string chimeraFlag = "no";
 405                         if (leftPiece.flag == "yes") {
 406
 407                                 if(  (leftPiece.results[0].bsa >= minBS && leftPiece.results[0].divr_qla_qrb >= divR)
 408                                         ||
 409                                         (leftPiece.results[0].bsb >= minBS && leftPiece.results[0].divr_qlb_qra >= divR) ) { chimeraFlag = "yes"; }
 410                         }
 411
 412                         if (rightPiece.flag == "yes") {
 413                                 if ( (rightPiece.results[0].bsa >= minBS && rightPiece.results[0].divr_qla_qrb >= divR)
 414                                  ||
 415                                  (rightPiece.results[0].bsb >= minBS && rightPiece.results[0].divr_qlb_qra >= divR) ) { chimeraFlag = "yes"; }
 416                         }
 417
 418                         bool rightChimeric = false;
 419                         bool leftChimeric = false;
 420
 421                         if (chimeraFlag == "yes") {
 422                                 //which peice is chimeric or are both
 423                                 if (rightPiece.flag == "yes") { if ((rightPiece.results[0].bsa >= minBS) || (rightPiece.results[0].bsb >= minBS)) { rightChimeric = true; } }
 424                                 if (leftPiece.flag == "yes") { if ((leftPiece.results[0].bsa >= minBS) || (leftPiece.results[0].bsb >= minBS))  { leftChimeric = true;  } }
 425
 426                                 if (rightChimeric || leftChimeric) {
 427                                         m->mothurOut(querySeq->getName() + "\tyes"); m->mothurOutEndLine();
 428                                         outAcc << querySeq->getName() << endl;
 429
 430                                         if (templateFileName == "self") {  chimericSeqs.insert(querySeq->getName()); }
 431
 432                                         if (trimChimera) {
 433                                                 string newAligned = trim->getAligned();
 434
 435                                                 //right side is fine so keep that
 436                                                 if ((leftChimeric) && (!rightChimeric)) {
 437                                                         for (int i = 0; i < leftPiece.spotMap[leftPiece.results[0].winREnd]; i++) { newAligned[i] = '.'; }
 438                                                 }else if ((!leftChimeric) && (rightChimeric)) { //leftside is fine so keep that
 439                                                         for (int i = (rightPiece.spotMap[rightPiece.results[0].winLStart]-1); i < newAligned.length(); i++) { newAligned[i] = '.'; }
 440                                                 }else { //both sides are chimeric, keep longest piece
 441
 442                                                         int lengthLeftLeft = leftPiece.spotMap[leftPiece.results[0].winLEnd] - leftPiece.spotMap[leftPiece.results[0].winLStart];
 443                                                         int lengthLeftRight = leftPiece.spotMap[leftPiece.results[0].winREnd] - leftPiece.spotMap[leftPiece.results[0].winRStart];
 444
 445                                                         int longest = 1; // leftleft = 1, leftright = 2, rightleft = 3 rightright = 4
 446                                                         int length = lengthLeftLeft;
 447                                                         if (lengthLeftLeft < lengthLeftRight) { longest = 2;  length = lengthLeftRight; }
 448
 449                                                         int lengthRightLeft = rightPiece.spotMap[rightPiece.results[0].winLEnd] - rightPiece.spotMap[rightPiece.results[0].winLStart];
 450                                                         int lengthRightRight = rightPiece.spotMap[rightPiece.results[0].winREnd] - rightPiece.spotMap[rightPiece.results[0].winRStart];
 451
 452                                                         if (lengthRightLeft > length) { longest = 3; length = lengthRightLeft;  }
 453                                                         if (lengthRightRight > length) { longest = 4; }
 454
 455                                                         if (longest == 1) { //leftleft
 456                                                                 for (int i = (leftPiece.spotMap[leftPiece.results[0].winRStart]-1); i < newAligned.length(); i++) { newAligned[i] = '.'; }
 457                                                         }else if (longest == 2) { //leftright
 458                                                                 //get rid of leftleft
 459                                                                 for (int i = (leftPiece.spotMap[leftPiece.results[0].winLStart]-1); i < (leftPiece.spotMap[leftPiece.results[0].winLEnd]-1); i++) { newAligned[i] = '.'; }
 460                                                                 //get rid of right
 461                                                                 for (int i = (rightPiece.spotMap[rightPiece.results[0].winLStart]-1); i < newAligned.length(); i++) { newAligned[i] = '.'; }
 462                                                         }else if (longest == 3) { //rightleft
 463                                                                 //get rid of left
 464                                                                 for (int i = 0; i < leftPiece.spotMap[leftPiece.results[0].winREnd]; i++) { newAligned[i] = '.'; }
 465                                                                 //get rid of rightright
 466                                                                 for (int i = (rightPiece.spotMap[rightPiece.results[0].winRStart]-1); i < newAligned.length(); i++) { newAligned[i] = '.'; }
 467                                                         }else { //rightright
 468                                                                 //get rid of left
 469                                                                 for (int i = 0; i < leftPiece.spotMap[leftPiece.results[0].winREnd]; i++) { newAligned[i] = '.'; }
 470                                                                 //get rid of rightleft
 471                                                                 for (int i = (rightPiece.spotMap[rightPiece.results[0].winLStart]-1); i < (rightPiece.spotMap[rightPiece.results[0].winLEnd]-1); i++) { newAligned[i] = '.'; }
 472                                                         }
 473                                                 }
 474
 475                                                 trim->setAligned(newAligned);
 476                                         }
 477
 478                                 }
 479                         }
 480
 481                         printBlock(leftPiece, rightPiece, leftChimeric, rightChimeric, chimeraFlag, out);
 482                         out << endl;
 483                 }else {
 484                         out << querySeq->getName() << "\tno" << endl;
 485                 }
 486
 487                 return trim;
 488
 489         }
 490         catch(exception& e) {
 491                 m->errorOut(e, "ChimeraSlayer", "print");
 492                 exit(1);
 493         }
 494 }
 495
 496 #ifdef USE_MPI
 497 //***************************************************************************************************************
 498 Sequence* ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc, data_results leftPiece, data_results rightPiece) {
 499         try {
 500                 MPI_Status status;
 501                 bool results = false;
 502                 string outAccString = "";
 503                 string outputString = "";
 504
 505                 Sequence* trim = NULL;
 506
 507                 if (trimChimera) {
 508                         string aligned = leftPiece.trimQuery.getAligned() + rightPiece.trimQuery.getAligned();
 509                         trim = new Sequence(leftPiece.trimQuery.getName(), aligned);
 510                 }
 511
 512
 513                 if ((leftPiece.flag == "yes") || (rightPiece.flag == "yes")) {
 514
 515                         string chimeraFlag = "no";
 516                         if (leftPiece.flag == "yes") {
 517
 518                                 if(  (leftPiece.results[0].bsa >= minBS && leftPiece.results[0].divr_qla_qrb >= divR)
 519                                    ||
 520                                    (leftPiece.results[0].bsb >= minBS && leftPiece.results[0].divr_qlb_qra >= divR) ) { chimeraFlag = "yes"; }
 521                         }
 522
 523                         if (rightPiece.flag == "yes") {
 524                                 if ( (rightPiece.results[0].bsa >= minBS && rightPiece.results[0].divr_qla_qrb >= divR)
 525                                         ||
 526                                         (rightPiece.results[0].bsb >= minBS && rightPiece.results[0].divr_qlb_qra >= divR) ) { chimeraFlag = "yes"; }
 527                         }
 528
 529                         bool rightChimeric = false;
 530                         bool leftChimeric = false;
 531
 532                         if (chimeraFlag == "yes") {
 533                                 //which peice is chimeric or are both
 534                                 if (rightPiece.flag == "yes") { if ((rightPiece.results[0].bsa >= minBS) || (rightPiece.results[0].bsb >= minBS)) { rightChimeric = true; } }
 535                                 if (leftPiece.flag == "yes") { if ((leftPiece.results[0].bsa >= minBS) || (leftPiece.results[0].bsb >= minBS))  { leftChimeric = true;  } }
 536
 537                                 if (rightChimeric || leftChimeric) {
 538                                         cout << querySeq->getName() <<  "\tyes" << endl;
 539                                         outAccString += querySeq->getName() + "\n";
 540                                         results = true;
 541
 542                                         if (templateFileName == "self") {  chimericSeqs.insert(querySeq->getName()); }
 543
 544                                         //write to accnos file
 545                                         int length = outAccString.length();
 546                                         char* buf2 = new char[length];
 547                                         memcpy(buf2, outAccString.c_str(), length);
 548
 549                                         MPI_File_write_shared(outAcc, buf2, length, MPI_CHAR, &status);
 550                                         delete buf2;
 551
 552                                         if (trimChimera) {
 553                                                 string newAligned = trim->getAligned();
 554
 555                                                 //right side is fine so keep that
 556                                                 if ((leftChimeric) && (!rightChimeric)) {
 557                                                         for (int i = 0; i < leftPiece.spotMap[leftPiece.results[0].winREnd]; i++) { newAligned[i] = '.'; }
 558                                                 }else if ((!leftChimeric) && (rightChimeric)) { //leftside is fine so keep that
 559                                                         for (int i = (rightPiece.spotMap[rightPiece.results[0].winLStart]-1); i < newAligned.length(); i++) { newAligned[i] = '.'; }
 560                                                 }else { //both sides are chimeric, keep longest piece
 561
 562                                                         int lengthLeftLeft = leftPiece.spotMap[leftPiece.results[0].winLEnd] - leftPiece.spotMap[leftPiece.results[0].winLStart];
 563                                                         int lengthLeftRight = leftPiece.spotMap[leftPiece.results[0].winREnd] - leftPiece.spotMap[leftPiece.results[0].winRStart];
 564
 565                                                         int longest = 1; // leftleft = 1, leftright = 2, rightleft = 3 rightright = 4
 566                                                         int length = lengthLeftLeft;
 567                                                         if (lengthLeftLeft < lengthLeftRight) { longest = 2;  length = lengthLeftRight; }
 568
 569                                                         int lengthRightLeft = rightPiece.spotMap[rightPiece.results[0].winLEnd] - rightPiece.spotMap[rightPiece.results[0].winLStart];
 570                                                         int lengthRightRight = rightPiece.spotMap[rightPiece.results[0].winREnd] - rightPiece.spotMap[rightPiece.results[0].winRStart];
 571
 572                                                         if (lengthRightLeft > length) { longest = 3; length = lengthRightLeft;  }
 573                                                         if (lengthRightRight > length) { longest = 4; }
 574
 575                                                         if (longest == 1) { //leftleft
 576                                                                 for (int i = (leftPiece.spotMap[leftPiece.results[0].winRStart]-1); i < newAligned.length(); i++) { newAligned[i] = '.'; }
 577                                                         }else if (longest == 2) { //leftright
 578                                                                 //get rid of leftleft
 579                                                                 for (int i = (leftPiece.spotMap[leftPiece.results[0].winLStart]-1); i < (leftPiece.spotMap[leftPiece.results[0].winLEnd]-1); i++) { newAligned[i] = '.'; }
 580                                                                 //get rid of right
 581                                                                 for (int i = (rightPiece.spotMap[rightPiece.results[0].winLStart]-1); i < newAligned.length(); i++) { newAligned[i] = '.'; }
 582                                                         }else if (longest == 3) { //rightleft
 583                                                                 //get rid of left
 584                                                                 for (int i = 0; i < leftPiece.spotMap[leftPiece.results[0].winREnd]; i++) { newAligned[i] = '.'; }
 585                                                                 //get rid of rightright
 586                                                                 for (int i = (rightPiece.spotMap[rightPiece.results[0].winRStart]-1); i < newAligned.length(); i++) { newAligned[i] = '.'; }
 587                                                         }else { //rightright
 588                                                                 //get rid of left
 589                                                                 for (int i = 0; i < leftPiece.spotMap[leftPiece.results[0].winREnd]; i++) { newAligned[i] = '.'; }
 590                                                                 //get rid of rightleft
 591                                                                 for (int i = (rightPiece.spotMap[rightPiece.results[0].winLStart]-1); i < (rightPiece.spotMap[rightPiece.results[0].winLEnd]-1); i++) { newAligned[i] = '.'; }
 592                                                         }
 593                                                 }
 594
 595                                                 trim->setAligned(newAligned);
 596                                         }
 597
 598                                 }
 599                         }
 600
 601                         outputString = getBlock(leftPiece, rightPiece, leftChimeric, rightChimeric, chimeraFlag);
 602                         outputString += "\n";
 603
 604                         //write to output file
 605                         int length = outputString.length();
 606                         char* buf = new char[length];
 607                         memcpy(buf, outputString.c_str(), length);
 608
 609                         MPI_File_write_shared(out, buf, length, MPI_CHAR, &status);
 610                         delete buf;
 611
 612                 }else {
 613                         outputString += querySeq->getName() + "\tno\n";
 614
 615                         //write to output file
 616                         int length = outputString.length();
 617                         char* buf = new char[length];
 618                         memcpy(buf, outputString.c_str(), length);
 619
 620                         MPI_File_write_shared(out, buf, length, MPI_CHAR, &status);
 621                         delete buf;
 622                 }
 623
 624
 625                 return trim;
 626         }
 627         catch(exception& e) {
 628                 m->errorOut(e, "ChimeraSlayer", "print");
 629                 exit(1);
 630         }
 631 }
 632 //***************************************************************************************************************
 633 Sequence* ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc) {
 634         try {
 635                 MPI_Status status;
 636                 bool results = false;
 637                 string outAccString = "";
 638                 string outputString = "";
 639
 640                 Sequence* trim = NULL;
 641                 if (trimChimera) { trim = new Sequence(trimQuery.getName(), trimQuery.getAligned()); }
 642
 643                 if (chimeraFlags == "yes") {
 644                         string chimeraFlag = "no";
 645                         if(  (chimeraResults[0].bsa >= minBS && chimeraResults[0].divr_qla_qrb >= divR)
 646                            ||
 647                            (chimeraResults[0].bsb >= minBS && chimeraResults[0].divr_qlb_qra >= divR) ) { chimeraFlag = "yes"; }
 648
 649
 650                         if (chimeraFlag == "yes") {
 651                                 if ((chimeraResults[0].bsa >= minBS) || (chimeraResults[0].bsb >= minBS)) {
 652                                         cout << querySeq->getName() <<  "\tyes" << endl;
 653                                         outAccString += querySeq->getName() + "\n";
 654                                         results = true;
 655
 656                                         if (templateFileName == "self") {  chimericSeqs.insert(querySeq->getName()); }
 657
 658                                         //write to accnos file
 659                                         int length = outAccString.length();
 660                                         char* buf2 = new char[length];
 661                                         memcpy(buf2, outAccString.c_str(), length);
 662
 663                                         MPI_File_write_shared(outAcc, buf2, length, MPI_CHAR, &status);
 664                                         delete buf2;
 665
 666                                         if (trimChimera) {
 667                                                 int lengthLeft = spotMap[chimeraResults[0].winLEnd] - spotMap[chimeraResults[0].winLStart];
 668                                                 int lengthRight = spotMap[chimeraResults[0].winREnd] - spotMap[chimeraResults[0].winRStart];
 669
 670                                                 string newAligned = trim->getAligned();
 671                                                 if (lengthLeft > lengthRight) { //trim right
 672                                                         for (int i = (spotMap[chimeraResults[0].winRStart]-1); i < newAligned.length(); i++) { newAligned[i] = '.'; }
 673                                                 }else { //trim left
 674                                                         for (int i = 0; i < (spotMap[chimeraResults[0].winLEnd]-1); i++) { newAligned[i] = '.'; }
 675                                                 }
 676                                                 trim->setAligned(newAligned);
 677                                         }
 678                                 }
 679                         }
 680
 681                         outputString = getBlock(chimeraResults[0], chimeraFlag);
 682                         outputString += "\n";
 683
 684                         //write to output file
 685                         int length = outputString.length();
 686                         char* buf = new char[length];
 687                         memcpy(buf, outputString.c_str(), length);
 688
 689                         MPI_File_write_shared(out, buf, length, MPI_CHAR, &status);
 690                         delete buf;
 691
 692                 }else {
 693                         outputString += querySeq->getName() + "\tno\n";
 694
 695                         //write to output file
 696                         int length = outputString.length();
 697                         char* buf = new char[length];
 698                         memcpy(buf, outputString.c_str(), length);
 699
 700                         MPI_File_write_shared(out, buf, length, MPI_CHAR, &status);
 701                         delete buf;
 702                 }
 703
 704                 return trim;
 705         }
 706         catch(exception& e) {
 707                 m->errorOut(e, "ChimeraSlayer", "print");
 708                 exit(1);
 709         }
 710 }
 711 #endif
 712
 713 //***************************************************************************************************************
 714 int ChimeraSlayer::getChimeras(Sequence* query) {
 715         try {
 716
 717                 trimQuery.setName(query->getName()); trimQuery.setAligned(query->getAligned());
 718                 printResults.trimQuery = trimQuery;
 719
 720                 chimeraFlags = "no";
 721                 printResults.flag = "no";
 722
 723                 //filter query
 724                 spotMap = runFilter(query);
 725                 printResults.spotMap = spotMap;
 726
 727                 querySeq = query;
 728
 729                 //you must create a template
 730                 vector<Sequence*> thisTemplate;
 731                 if (templateFileName != "self") { thisTemplate = templateSeqs; }
 732                 else {  thisTemplate = getTemplate(query);  } //fills this template and creates the databases
 733
 734                 if (m->control_pressed) {  return 0;  }
 735
 736                 if (thisTemplate.size() == 0) {  return 0; } //not chimeric
 737
 738                 //referenceSeqs, numWanted, matchScore, misMatchPenalty, divR, minSimilarity
 739                 Maligner maligner(thisTemplate, numWanted, match, misMatch, divR, minSim, minCov, searchMethod, databaseLeft, databaseRight);
 740                 Slayer slayer(window, increment, minSim, divR, iters, minSNP);
 741
 742                 if (templateFileName == "self") {
 743                         if (searchMethod == "kmer") {  delete databaseRight;  delete databaseLeft;  }
 744                         else if (searchMethod == "blast") {  delete databaseLeft; }
 745                 }
 746
 747                 if (m->control_pressed) {  return 0;  }
 748
 749                 string chimeraFlag = maligner.getResults(query, decalc);
 750
 751                 if (m->control_pressed) {  return 0;  }
 752
 753                 vector<results> Results = maligner.getOutput();
 754
 755                 //cout << query->getName() << endl;
 756                 //for (int i = 0; i < Results.size(); i++) {
 757                         //cout << Results[i].parent << '\t' << Results[i].regionStart << '\t' << Results[i].regionEnd << '\t' << Results[i].nastRegionStart << '\t' << Results[i].nastRegionEnd << '\t' << Results[i].queryToParent << '\t' << Results[i].queryToParentLocal << endl;
 758                 //}
 759                 //cout << "done\n" << endl;
 760                 if (chimeraFlag == "yes") {
 761
 762                         if (realign) {
 763                                 ChimeraReAligner realigner(thisTemplate, match, misMatch);
 764                                 realigner.reAlign(query, Results);
 765                         }
 766
 767                         //get sequence that were given from maligner results
 768                         vector<SeqDist> seqs;
 769                         map<string, float> removeDups;
 770                         map<string, float>::iterator itDup;
 771                         map<string, string> parentNameSeq;
 772                         map<string, string>::iterator itSeq;
 773                         for (int j = 0; j < Results.size(); j++) {
 774                                 float dist = (Results[j].regionEnd - Results[j].regionStart + 1) * Results[j].queryToParentLocal;
 775                                 //only add if you are not a duplicate
 776                                 itDup = removeDups.find(Results[j].parent);
 777                                 if (itDup == removeDups.end()) { //this is not duplicate
 778                                         removeDups[Results[j].parent] = dist;
 779                                         parentNameSeq[Results[j].parent] = Results[j].parentAligned;
 780                                 }else if (dist > itDup->second) { //is this a stronger number for this parent
 781                                         removeDups[Results[j].parent] = dist;
 782                                         parentNameSeq[Results[j].parent] = Results[j].parentAligned;
 783                                 }
 784                         }
 785
 786                         for (itDup = removeDups.begin(); itDup != removeDups.end(); itDup++) {
 787                                 itSeq = parentNameSeq.find(itDup->first);
 788                                 Sequence* seq = new Sequence(itDup->first, itSeq->second);
 789
 790                                 SeqDist member;
 791                                 member.seq = seq;
 792                                 member.dist = itDup->second;
 793
 794                                 seqs.push_back(member);
 795                         }
 796
 797                         //limit number of parents to explore - default 3
 798                         if (Results.size() > parents) {
 799                                 //sort by distance
 800                                 sort(seqs.begin(), seqs.end(), compareSeqDist);
 801                                 //prioritize larger more similiar sequence fragments
 802                                 reverse(seqs.begin(), seqs.end());
 803
 804                                 for (int k = seqs.size()-1; k > (parents-1); k--)  {
 805                                         delete seqs[k].seq;
 806                                         seqs.pop_back();
 807                                 }
 808                         }
 809
 810                         //put seqs into vector to send to slayer
 811                         vector<Sequence*> seqsForSlayer;
 812
 813                         for (int k = 0; k < seqs.size(); k++) {  seqsForSlayer.push_back(seqs[k].seq);  }
 814
 815                         //mask then send to slayer...
 816                         if (seqMask != "") {
 817                                 decalc->setMask(seqMask);
 818
 819                                 //mask querys
 820                                 decalc->runMask(query);
 821
 822                                 //mask parents
 823                                 for (int k = 0; k < seqsForSlayer.size(); k++) {
 824                                         decalc->runMask(seqsForSlayer[k]);
 825                                 }
 826
 827                                 spotMap = decalc->getMaskMap();
 828                         }
 829
 830                         if (m->control_pressed) {  for (int k = 0; k < seqs.size(); k++) {  delete seqs[k].seq;   }  return 0;  }
 831
 832                         //send to slayer
 833                         chimeraFlags = slayer.getResults(query, seqsForSlayer);
 834                         if (m->control_pressed) {  return 0;  }
 835                         chimeraResults = slayer.getOutput();
 836
 837                         //free memory
 838                         for (int k = 0; k < seqs.size(); k++) {  delete seqs[k].seq;   }
 839
 840                         printResults.spotMap = spotMap;
 841                         printResults.flag = chimeraFlags;
 842                         printResults.results = chimeraResults;
 843                 }
 844
 845                 return 0;
 846         }
 847         catch(exception& e) {
 848                 m->errorOut(e, "ChimeraSlayer", "getChimeras");
 849                 exit(1);
 850         }
 851 }
 852 //***************************************************************************************************************
 853 void ChimeraSlayer::printBlock(data_struct data, string flag, ostream& out){
 854         try {
 855                 out << querySeq->getName() << '\t';
 856                 out << data.parentA.getName() << "\t" << data.parentB.getName()  << '\t';
 857
 858                 out << data.divr_qla_qrb << '\t' << data.qla_qrb << '\t' << data.bsa << '\t';
 859                 out << data.divr_qlb_qra << '\t' << data.qlb_qra << '\t' << data.bsb << '\t';
 860
 861                 out << flag << '\t' << spotMap[data.winLStart] << "-" << spotMap[data.winLEnd] << '\t' << spotMap[data.winRStart] << "-" << spotMap[data.winREnd] << '\t';
 862
 863         }
 864         catch(exception& e) {
 865                 m->errorOut(e, "ChimeraSlayer", "printBlock");
 866                 exit(1);
 867         }
 868 }
 869 //***************************************************************************************************************
 870 void ChimeraSlayer::printBlock(data_results leftdata, data_results rightdata, bool leftChimeric, bool rightChimeric, string flag, ostream& out){
 871         try {
 872
 873                 if ((leftChimeric) && (!rightChimeric)) { //print left
 874                         out << querySeq->getName() << '\t';
 875                         out << leftdata.results[0].parentA.getName() << "\t" << leftdata.results[0].parentB.getName()  << '\t';
 876
 877                         out << leftdata.results[0].divr_qla_qrb << '\t' << leftdata.results[0].qla_qrb << '\t' << leftdata.results[0].bsa << '\t';
 878                         out << leftdata.results[0].divr_qlb_qra << '\t' << leftdata.results[0].qlb_qra << '\t' << leftdata.results[0].bsb << '\t';
 879
 880                         out << flag << '\t' << leftdata.spotMap[leftdata.results[0].winLStart] << "-" << leftdata.spotMap[leftdata.results[0].winLEnd] << '\t' << leftdata.spotMap[leftdata.results[0].winRStart] << "-" << leftdata.spotMap[leftdata.results[0].winREnd] << '\t';
 881
 882                 }else if ((!leftChimeric) && (rightChimeric)) {  //print right
 883                         out << querySeq->getName() << '\t';
 884                         out << rightdata.results[0].parentA.getName() << "\t" << rightdata.results[0].parentB.getName()  << '\t';
 885
 886                         out << rightdata.results[0].divr_qla_qrb << '\t' << rightdata.results[0].qla_qrb << '\t' << rightdata.results[0].bsa << '\t';
 887                         out << rightdata.results[0].divr_qlb_qra << '\t' << rightdata.results[0].qlb_qra << '\t' << rightdata.results[0].bsb << '\t';
 888
 889                         out << flag << '\t' << rightdata.spotMap[rightdata.results[0].winLStart] << "-" << rightdata.spotMap[rightdata.results[0].winLEnd] << '\t' << rightdata.spotMap[rightdata.results[0].winRStart] << "-" << rightdata.spotMap[rightdata.results[0].winREnd] << '\t';
 890
 891                 }else  { //print both results
 892                         if (leftdata.flag == "yes") {
 893                                 out << querySeq->getName() + "_LEFT" << '\t';
 894                                 out << leftdata.results[0].parentA.getName() << "\t" << leftdata.results[0].parentB.getName()  << '\t';
 895
 896                                 out << leftdata.results[0].divr_qla_qrb << '\t' << leftdata.results[0].qla_qrb << '\t' << leftdata.results[0].bsa << '\t';
 897                                 out << leftdata.results[0].divr_qlb_qra << '\t' << leftdata.results[0].qlb_qra << '\t' << leftdata.results[0].bsb << '\t';
 898
 899                                 out << flag << '\t' << leftdata.spotMap[leftdata.results[0].winLStart] << "-" << leftdata.spotMap[leftdata.results[0].winLEnd] << '\t' << leftdata.spotMap[leftdata.results[0].winRStart] << "-" << leftdata.spotMap[leftdata.results[0].winREnd] << '\t';
 900                         }
 901
 902                         if (rightdata.flag == "yes") {
 903                                 if (leftdata.flag == "yes") { out << endl; }
 904
 905                                 out << querySeq->getName() + "_RIGHT"<< '\t';
 906                                 out << rightdata.results[0].parentA.getName() << "\t" << rightdata.results[0].parentB.getName()  << '\t';
 907
 908                                 out << rightdata.results[0].divr_qla_qrb << '\t' << rightdata.results[0].qla_qrb << '\t' << rightdata.results[0].bsa << '\t';
 909                                 out << rightdata.results[0].divr_qlb_qra << '\t' << rightdata.results[0].qlb_qra << '\t' << rightdata.results[0].bsb << '\t';
 910
 911                                 out << flag << '\t' << rightdata.spotMap[rightdata.results[0].winLStart] << "-" << rightdata.spotMap[rightdata.results[0].winLEnd] << '\t' << rightdata.spotMap[rightdata.results[0].winRStart] << "-" << rightdata.spotMap[rightdata.results[0].winREnd] << '\t';
 912
 913                         }
 914                 }
 915         }
 916         catch(exception& e) {
 917                 m->errorOut(e, "ChimeraSlayer", "printBlock");
 918                 exit(1);
 919         }
 920 }
 921 //***************************************************************************************************************
 922 string ChimeraSlayer::getBlock(data_results leftdata, data_results rightdata, bool leftChimeric, bool rightChimeric, string flag){
 923         try {
 924
 925                 string out = "";
 926
 927                 if ((leftChimeric) && (!rightChimeric)) { //get left
 928                         out += querySeq->getName() + "\t";
 929                         out += leftdata.results[0].parentA.getName() + "\t" + leftdata.results[0].parentB.getName() + "\t";
 930
 931                         out += toString(leftdata.results[0].divr_qla_qrb) + "\t" + toString(leftdata.results[0].qla_qrb) + "\t" + toString(leftdata.results[0].bsa) + "\t";
 932                         out += toString(leftdata.results[0].divr_qlb_qra) + "\t" + toString(leftdata.results[0].qlb_qra) + "\t" + toString(leftdata.results[0].bsb) + "\t";
 933
 934                         out += flag + "\t" + toString(leftdata.spotMap[leftdata.results[0].winLStart]) + "-" + toString(leftdata.spotMap[leftdata.results[0].winLEnd]) + "\t" + toString(leftdata.spotMap[leftdata.results[0].winRStart]) + "-" + toString(leftdata.spotMap[leftdata.results[0].winREnd]) + "\t";
 935
 936                 }else if ((!leftChimeric) && (rightChimeric)) {  //print right
 937                         out += querySeq->getName() + "\t";
 938                         out += rightdata.results[0].parentA.getName() + "\t" + rightdata.results[0].parentB.getName()  + "\t";
 939
 940                         out += toString(rightdata.results[0].divr_qla_qrb) + "\t" + toString(rightdata.results[0].qla_qrb) + "\t" + toString(rightdata.results[0].bsa) + "\t";
 941                         out += toString(rightdata.results[0].divr_qlb_qra) + "\t" + toString(rightdata.results[0].qlb_qra) + "\t" + toString(rightdata.results[0].bsb) + "\t";
 942
 943                         out += flag + "\t" + toString(rightdata.spotMap[rightdata.results[0].winLStart]) + "-" + toString(rightdata.spotMap[rightdata.results[0].winLEnd]) + "\t" + toString(rightdata.spotMap[rightdata.results[0].winRStart]) + "-" + toString(rightdata.spotMap[rightdata.results[0].winREnd]) + "\t";
 944
 945                 }else  { //print both results
 946
 947                         if (leftdata.flag == "yes") {
 948                                 out += querySeq->getName() + "_LEFT\t";
 949                                 out += leftdata.results[0].parentA.getName() + "\t" + leftdata.results[0].parentB.getName() + "\t";
 950
 951                                 out += toString(leftdata.results[0].divr_qla_qrb) + "\t" + toString(leftdata.results[0].qla_qrb) + "\t" + toString(leftdata.results[0].bsa) + "\t";
 952                                 out += toString(leftdata.results[0].divr_qlb_qra) + "\t" + toString(leftdata.results[0].qlb_qra) + "\t" + toString(leftdata.results[0].bsb) + "\t";
 953
 954                                 out += flag + "\t" + toString(leftdata.spotMap[leftdata.results[0].winLStart]) + "-" + toString(leftdata.spotMap[leftdata.results[0].winLEnd]) + "\t" + toString(leftdata.spotMap[leftdata.results[0].winRStart]) + "-" + toString(leftdata.spotMap[leftdata.results[0].winREnd]) + "\t";
 955                         }
 956
 957                         if (rightdata.flag == "yes") {
 958                                 if (leftdata.flag == "yes") { out += "\n"; }
 959                                 out +=  querySeq->getName() + "_RIGHT\t";
 960                                 out += rightdata.results[0].parentA.getName() + "\t" + rightdata.results[0].parentB.getName()  + "\t";
 961
 962                                 out += toString(rightdata.results[0].divr_qla_qrb) + "\t" + toString(rightdata.results[0].qla_qrb) + "\t" + toString(rightdata.results[0].bsa) + "\t";
 963                                 out += toString(rightdata.results[0].divr_qlb_qra) + "\t" + toString(rightdata.results[0].qlb_qra) + "\t" + toString(rightdata.results[0].bsb) + "\t";
 964
 965                                 out += flag + "\t" + toString(rightdata.spotMap[rightdata.results[0].winLStart]) + "-" + toString(rightdata.spotMap[rightdata.results[0].winLEnd]) + "\t" + toString(rightdata.spotMap[rightdata.results[0].winRStart]) + "-" + toString(rightdata.spotMap[rightdata.results[0].winREnd]) + "\t";
 966                         }
 967                 }
 968
 969                 return out;
 970
 971         }
 972         catch(exception& e) {
 973                 m->errorOut(e, "ChimeraSlayer", "getBlock");
 974                 exit(1);
 975         }
 976 }
 977 //***************************************************************************************************************
 978 string ChimeraSlayer::getBlock(data_struct data, string flag){
 979         try {
 980
 981                 string outputString = "";
 982
 983                 outputString += querySeq->getName() + "\t";
 984                 outputString += data.parentA.getName() + "\t" + data.parentB.getName()  + "\t";
 985
 986                 outputString += toString(data.divr_qla_qrb) + "\t" + toString(data.qla_qrb) + "\t" + toString(data.bsa) + "\t";
 987                 outputString += toString(data.divr_qlb_qra) + "\t" + toString(data.qlb_qra) + "\t" + toString(data.bsb) + "\t";
 988
 989                 outputString += flag + "\t" + toString(spotMap[data.winLStart]) + "-" + toString(spotMap[data.winLEnd]) + "\t" + toString(spotMap[data.winRStart]) + "-" + toString(spotMap[data.winREnd]) + "\t";
 990
 991                 return outputString;
 992         }
 993         catch(exception& e) {
 994                 m->errorOut(e, "ChimeraSlayer", "getBlock");
 995                 exit(1);
 996         }
 997 }
 998 //***************************************************************************************************************/
 999