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