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) {
 784                         trimQuery = new Sequence(query->getName(), query->getAligned());
 785                         printResults.trimQuery = *trimQuery;
 786                 }
 787
 788                 chimeraFlags = "no";
 789                 printResults.flag = "no";
 790
 791                 //filter query
 792                 spotMap = runFilter(query);
 793                 printResults.spotMap = spotMap;
 794
 795                 querySeq = query;
 796
 797                 //you must create a template
 798                 vector<Sequence*> thisTemplate;
 799                 if (templateFileName != "self") { thisTemplate = templateSeqs; }
 800                 else { thisTemplate = getTemplate(query); } //fills this template and creates the databases
 801
 802                 if (m->control_pressed) {  return 0;  }
 803
 804                 if (thisTemplate.size() == 0) {  return 0; } //not chimeric
 805
 806                 //referenceSeqs, numWanted, matchScore, misMatchPenalty, divR, minSimilarity
 807                 Maligner maligner(thisTemplate, numWanted, match, misMatch, divR, minSim, minCov, searchMethod, databaseLeft, databaseRight);
 808                 Slayer slayer(window, increment, minSim, divR, iters, minSNP);
 809
 810                 if (templateFileName == "self") {
 811                         if (searchMethod == "kmer") {  delete databaseRight;  delete databaseLeft;  }
 812                         else if (searchMethod == "blast") {  delete databaseLeft; }
 813                 }
 814
 815                 if (m->control_pressed) {  return 0;  }
 816
 817                 string chimeraFlag = maligner.getResults(query, decalc);
 818
 819                 if (m->control_pressed) {  return 0;  }
 820
 821                 vector<results> Results = maligner.getOutput();
 822
 823                 if (realign) {
 824                         ChimeraReAligner realigner(thisTemplate, match, misMatch);
 825                         realigner.reAlign(query, Results);
 826                 }
 827
 828                 if (chimeraFlag == "yes") {
 829
 830                         //get sequence that were given from maligner results
 831                         vector<SeqDist> seqs;
 832                         map<string, float> removeDups;
 833                         map<string, float>::iterator itDup;
 834                         map<string, string> parentNameSeq;
 835                         map<string, string>::iterator itSeq;
 836                         for (int j = 0; j < Results.size(); j++) {
 837                                 float dist = (Results[j].regionEnd - Results[j].regionStart + 1) * Results[j].queryToParentLocal;
 838                                 //only add if you are not a duplicate
 839                                 itDup = removeDups.find(Results[j].parent);
 840                                 if (itDup == removeDups.end()) { //this is not duplicate
 841                                         removeDups[Results[j].parent] = dist;
 842                                         parentNameSeq[Results[j].parent] = Results[j].parentAligned;
 843                                 }else if (dist > itDup->second) { //is this a stronger number for this parent
 844                                         removeDups[Results[j].parent] = dist;
 845                                         parentNameSeq[Results[j].parent] = Results[j].parentAligned;
 846                                 }
 847                         }
 848
 849                         for (itDup = removeDups.begin(); itDup != removeDups.end(); itDup++) {
 850                                 itSeq = parentNameSeq.find(itDup->first);
 851                                 Sequence* seq = new Sequence(itDup->first, itSeq->second);
 852
 853                                 SeqDist member;
 854                                 member.seq = seq;
 855                                 member.dist = itDup->second;
 856
 857                                 seqs.push_back(member);
 858                         }
 859
 860                         //limit number of parents to explore - default 3
 861                         if (Results.size() > parents) {
 862                                 //sort by distance
 863                                 sort(seqs.begin(), seqs.end(), compareSeqDist);
 864                                 //prioritize larger more similiar sequence fragments
 865                                 reverse(seqs.begin(), seqs.end());
 866
 867                                 for (int k = seqs.size()-1; k > (parents-1); k--)  {
 868                                         delete seqs[k].seq;
 869                                         seqs.pop_back();
 870                                 }
 871                         }
 872
 873                         //put seqs into vector to send to slayer
 874                         vector<Sequence*> seqsForSlayer;
 875
 876                         for (int k = 0; k < seqs.size(); k++) {  seqsForSlayer.push_back(seqs[k].seq);  }
 877
 878                         //mask then send to slayer...
 879                         if (seqMask != "") {
 880                                 decalc->setMask(seqMask);
 881
 882                                 //mask querys
 883                                 decalc->runMask(query);
 884
 885                                 //mask parents
 886                                 for (int k = 0; k < seqsForSlayer.size(); k++) {
 887                                         decalc->runMask(seqsForSlayer[k]);
 888                                 }
 889
 890                                 spotMap = decalc->getMaskMap();
 891                         }
 892
 893                         if (m->control_pressed) {  for (int k = 0; k < seqs.size(); k++) {  delete seqs[k].seq;   }  return 0;  }
 894
 895                         //send to slayer
 896                         chimeraFlags = slayer.getResults(query, seqsForSlayer);
 897                         if (m->control_pressed) {  return 0;  }
 898                         chimeraResults = slayer.getOutput();
 899
 900                         //free memory
 901                         for (int k = 0; k < seqs.size(); k++) {  delete seqs[k].seq;   }
 902
 903                         printResults.spotMap = spotMap;
 904                         printResults.flag = chimeraFlags;
 905                         printResults.results = chimeraResults;
 906                 }
 907
 908                 return 0;
 909         }
 910         catch(exception& e) {
 911                 m->errorOut(e, "ChimeraSlayer", "getChimeras");
 912                 exit(1);
 913         }
 914 }
 915 //***************************************************************************************************************
 916 void ChimeraSlayer::printBlock(data_struct data, string flag, ostream& out){
 917         try {
 918                 out << querySeq->getName() << '\t';
 919                 out << data.parentA.getName() << "\t" << data.parentB.getName()  << '\t';
 920
 921                 out << data.divr_qla_qrb << '\t' << data.qla_qrb << '\t' << data.bsa << '\t';
 922                 out << data.divr_qlb_qra << '\t' << data.qlb_qra << '\t' << data.bsb << '\t';
 923
 924                 out << flag << '\t' << spotMap[data.winLStart] << "-" << spotMap[data.winLEnd] << '\t' << spotMap[data.winRStart] << "-" << spotMap[data.winREnd] << '\t';
 925
 926         }
 927         catch(exception& e) {
 928                 m->errorOut(e, "ChimeraSlayer", "printBlock");
 929                 exit(1);
 930         }
 931 }
 932 //***************************************************************************************************************
 933 void ChimeraSlayer::printBlock(data_results leftdata, data_results rightdata, bool leftChimeric, bool rightChimeric, string flag, ostream& out){
 934         try {
 935
 936                 if ((leftChimeric) && (!rightChimeric)) { //print left
 937                         out << querySeq->getName() << '\t';
 938                         out << leftdata.results[0].parentA.getName() << "\t" << leftdata.results[0].parentB.getName()  << '\t';
 939
 940                         out << leftdata.results[0].divr_qla_qrb << '\t' << leftdata.results[0].qla_qrb << '\t' << leftdata.results[0].bsa << '\t';
 941                         out << leftdata.results[0].divr_qlb_qra << '\t' << leftdata.results[0].qlb_qra << '\t' << leftdata.results[0].bsb << '\t';
 942
 943                         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';
 944
 945                 }else if ((!leftChimeric) && (rightChimeric)) {  //print right
 946                         out << querySeq->getName() << '\t';
 947                         out << rightdata.results[0].parentA.getName() << "\t" << rightdata.results[0].parentB.getName()  << '\t';
 948
 949                         out << rightdata.results[0].divr_qla_qrb << '\t' << rightdata.results[0].qla_qrb << '\t' << rightdata.results[0].bsa << '\t';
 950                         out << rightdata.results[0].divr_qlb_qra << '\t' << rightdata.results[0].qlb_qra << '\t' << rightdata.results[0].bsb << '\t';
 951
 952                         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';
 953
 954                 }else  { //print both results
 955                         if (leftdata.flag == "yes") {
 956                                 out << querySeq->getName() + "_LEFT" << '\t';
 957                                 out << leftdata.results[0].parentA.getName() << "\t" << leftdata.results[0].parentB.getName()  << '\t';
 958
 959                                 out << leftdata.results[0].divr_qla_qrb << '\t' << leftdata.results[0].qla_qrb << '\t' << leftdata.results[0].bsa << '\t';
 960                                 out << leftdata.results[0].divr_qlb_qra << '\t' << leftdata.results[0].qlb_qra << '\t' << leftdata.results[0].bsb << '\t';
 961
 962                                 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';
 963                         }
 964
 965                         if (rightdata.flag == "yes") {
 966                                 if (leftdata.flag == "yes") { out << endl; }
 967
 968                                 out << querySeq->getName() + "_RIGHT"<< '\t';
 969                                 out << rightdata.results[0].parentA.getName() << "\t" << rightdata.results[0].parentB.getName()  << '\t';
 970
 971                                 out << rightdata.results[0].divr_qla_qrb << '\t' << rightdata.results[0].qla_qrb << '\t' << rightdata.results[0].bsa << '\t';
 972                                 out << rightdata.results[0].divr_qlb_qra << '\t' << rightdata.results[0].qlb_qra << '\t' << rightdata.results[0].bsb << '\t';
 973
 974                                 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';
 975
 976                         }
 977                 }
 978         }
 979         catch(exception& e) {
 980                 m->errorOut(e, "ChimeraSlayer", "printBlock");
 981                 exit(1);
 982         }
 983 }
 984 //***************************************************************************************************************
 985 string ChimeraSlayer::getBlock(data_results leftdata, data_results rightdata, bool leftChimeric, bool rightChimeric, string flag){
 986         try {
 987
 988                 string out = "";
 989
 990                 if ((leftChimeric) && (!rightChimeric)) { //get left
 991                         out += querySeq->getName() + "\t";
 992                         out += leftdata.results[0].parentA.getName() + "\t" + leftdata.results[0].parentB.getName() + "\t";
 993
 994                         out += toString(leftdata.results[0].divr_qla_qrb) + "\t" + toString(leftdata.results[0].qla_qrb) + "\t" + toString(leftdata.results[0].bsa) + "\t";
 995                         out += toString(leftdata.results[0].divr_qlb_qra) + "\t" + toString(leftdata.results[0].qlb_qra) + "\t" + toString(leftdata.results[0].bsb) + "\t";
 996
 997                         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";
 998
 999                 }else if ((!leftChimeric) && (rightChimeric)) {  //print right
1000                         out += querySeq->getName() + "\t";
1001                         out += rightdata.results[0].parentA.getName() + "\t" + rightdata.results[0].parentB.getName()  + "\t";
1002
1003                         out += toString(rightdata.results[0].divr_qla_qrb) + "\t" + toString(rightdata.results[0].qla_qrb) + "\t" + toString(rightdata.results[0].bsa) + "\t";
1004                         out += toString(rightdata.results[0].divr_qlb_qra) + "\t" + toString(rightdata.results[0].qlb_qra) + "\t" + toString(rightdata.results[0].bsb) + "\t";
1005
1006                         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";
1007
1008                 }else  { //print both results
1009
1010                         if (leftdata.flag == "yes") {
1011                                 out += querySeq->getName() + "_LEFT\t";
1012                                 out += leftdata.results[0].parentA.getName() + "\t" + leftdata.results[0].parentB.getName() + "\t";
1013
1014                                 out += toString(leftdata.results[0].divr_qla_qrb) + "\t" + toString(leftdata.results[0].qla_qrb) + "\t" + toString(leftdata.results[0].bsa) + "\t";
1015                                 out += toString(leftdata.results[0].divr_qlb_qra) + "\t" + toString(leftdata.results[0].qlb_qra) + "\t" + toString(leftdata.results[0].bsb) + "\t";
1016
1017                                 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";
1018                         }
1019
1020                         if (rightdata.flag == "yes") {
1021                                 if (leftdata.flag == "yes") { out += "\n"; }
1022                                 out +=  querySeq->getName() + "_RIGHT\t";
1023                                 out += rightdata.results[0].parentA.getName() + "\t" + rightdata.results[0].parentB.getName()  + "\t";
1024
1025                                 out += toString(rightdata.results[0].divr_qla_qrb) + "\t" + toString(rightdata.results[0].qla_qrb) + "\t" + toString(rightdata.results[0].bsa) + "\t";
1026                                 out += toString(rightdata.results[0].divr_qlb_qra) + "\t" + toString(rightdata.results[0].qlb_qra) + "\t" + toString(rightdata.results[0].bsb) + "\t";
1027
1028                                 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";
1029                         }
1030                 }
1031
1032                 return out;
1033
1034         }
1035         catch(exception& e) {
1036                 m->errorOut(e, "ChimeraSlayer", "getBlock");
1037                 exit(1);
1038         }
1039 }
1040 //***************************************************************************************************************
1041 string ChimeraSlayer::getBlock(data_struct data, string flag){
1042         try {
1043
1044                 string outputString = "";
1045
1046                 outputString += querySeq->getName() + "\t";
1047                 outputString += data.parentA.getName() + "\t" + data.parentB.getName()  + "\t";
1048
1049                 outputString += toString(data.divr_qla_qrb) + "\t" + toString(data.qla_qrb) + "\t" + toString(data.bsa) + "\t";
1050                 outputString += toString(data.divr_qlb_qra) + "\t" + toString(data.qlb_qra) + "\t" + toString(data.bsb) + "\t";
1051
1052                 outputString += flag + "\t" + toString(spotMap[data.winLStart]) + "-" + toString(spotMap[data.winLEnd]) + "\t" + toString(spotMap[data.winRStart]) + "-" + toString(spotMap[data.winREnd]) + "\t";
1053
1054                 return outputString;
1055         }
1056         catch(exception& e) {
1057                 m->errorOut(e, "ChimeraSlayer", "getBlock");
1058                 exit(1);
1059         }
1060 }
1061 //***************************************************************************************************************/
1062