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