slayer.cpp

   1 /*
   2  *  slayer.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 9/25/09.
   6  *  Copyright 2009 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "slayer.h"
  11
  12 /***********************************************************************/
  13 Slayer::Slayer(int win, int increment, int parentThreshold, float div, int i) :
  14                 windowSize(win), windowStep(increment), parentFragmentThreshold(parentThreshold), divRThreshold(div), iters(i){}
  15 /***********************************************************************/
  16 string Slayer::getResults(Sequence* query, vector<Sequence*> refSeqs) {
  17         try {
  18                 vector<data_struct> all; all.clear();
  19
  20                 for (int i = 0; i < refSeqs.size(); i++) {
  21
  22                         for (int j = i+1; j < refSeqs.size(); j++) {
  23
  24                                 //make copies of query and each parent because runBellerophon removes gaps and messes them up
  25                                 Sequence* q = new Sequence(query->getName(), query->getAligned());
  26                                 Sequence* leftParent = new Sequence(refSeqs[i]->getName(), refSeqs[i]->getAligned());
  27                                 Sequence* rightParent = new Sequence(refSeqs[j]->getName(), refSeqs[j]->getAligned());
  28
  29                                 vector<data_struct> divs = runBellerophon(q, leftParent, rightParent);
  30
  31                                 vector<data_struct> selectedDivs;
  32                                 for (int k = 0; k < divs.size(); k++) {
  33
  34                                         vector<snps> snpsLeft = getSNPS(divs[k].parentA.getAligned(), divs[k].querySeq.getAligned(), divs[k].parentB.getAligned(), divs[k].winLStart, divs[k].winLEnd);
  35                                         vector<snps> snpsRight = getSNPS(divs[k].parentA.getAligned(), divs[k].querySeq.getAligned(), divs[k].parentB.getAligned(), divs[k].winRStart, divs[k].winREnd);
  36
  37                                         int numSNPSLeft = snpsLeft.size();
  38                                         int numSNPSRight = snpsRight.size();
  39
  40                                         //require at least 3 SNPs on each side of the break
  41                                         if ((numSNPSLeft >= 3) && (numSNPSRight >= 3)) {
  42
  43                                                 int winSizeLeft = divs[k].winLEnd - divs[k].winLStart + 1;
  44                                                 int winSizeRight = divs[k].winREnd - divs[k].winRStart + 1;
  45
  46                                                 float snpRateLeft = numSNPSLeft / (float) winSizeLeft;
  47                                                 float snpRateRight = numSNPSRight / (float) winSizeRight;
  48                                                 float logR = log(snpRateLeft / snpRateRight) / log(2);
  49
  50                                                 // do not accept excess snp ratio on either side of the break
  51                                                 if (abs(logR) < 1 ) {
  52
  53                                                         float BS_A, BS_B;
  54                                                         bootstrapSNPS(snpsLeft, snpsRight, BS_A, BS_B);
  55
  56                                                         divs[k].bsa = BS_A;
  57                                                         divs[k].bsb = BS_B;
  58
  59                                                         divs[k].bsMax = max(BS_A, BS_B);
  60
  61                                                         divs[k].chimeraMax = max(divs[k].qla_qrb, divs[k].qlb_qra);
  62
  63                                                         selectedDivs.push_back(divs[k]);
  64                                                 }
  65                                         }
  66                                 }
  67
  68                                 //save selected
  69                                 for (int m = 0; m < selectedDivs.size(); m++) {  all.push_back(selectedDivs[m]);        }
  70
  71                                 delete q;
  72                                 delete leftParent;
  73                                 delete rightParent;
  74                         }
  75                 }
  76
  77
  78                 // compute bootstrap support
  79                 if (all.size() > 0) {
  80                         //sort them
  81                         sort(all.begin(), all.end(), compareDataStruct);
  82                         reverse(all.begin(), all.end());
  83
  84                         outputResults = all;
  85                         return "yes";
  86                 }
  87                 else {
  88                         outputResults = all;
  89                         return "no";
  90                 }
  91         }
  92         catch(exception& e) {
  93                 errorOut(e, "Slayer", "getResults");
  94                 exit(1);
  95         }
  96 }
  97 /***********************************************************************/
  98 vector<data_struct> Slayer::runBellerophon(Sequence* q, Sequence* pA, Sequence* pB) {
  99         try{
 100
 101                 vector<data_struct> data;
 102
 103                 //vertical filter
 104                 vector<Sequence*> temp;
 105                 temp.push_back(q); temp.push_back(pA); temp.push_back(pB);
 106                 map<int, int> spots = verticalFilter(temp);
 107
 108                 //get these to avoid numerous function calls
 109                 string query = q->getAligned();
 110                 string parentA = pA->getAligned();
 111                 string parentB = pB->getAligned();
 112                 int length = query.length();
 113
 114                 //check window size
 115                 if (length < (2*windowSize+windowStep)) {
 116                         mothurOut("Your window size is too large for " + q->getName() + ". I will make the window size " + toString(length/4) + " which is 1/4 the filtered length."); mothurOutEndLine();
 117                         windowSize = length / 4;
 118                 }
 119
 120                 for (int i = windowSize-1; i <= (length - windowSize); i += windowStep) {
 121
 122                         int breakpoint = i;
 123                         int leftLength = breakpoint + 1;
 124                         int rightLength = length - leftLength;
 125
 126                         float QLA = computePercentID(query, parentA, 0, breakpoint);
 127                         float QRB = computePercentID(query, parentB, breakpoint+1, length - 1);
 128
 129                         float QLB = computePercentID(query, parentB, 0, breakpoint);
 130                         float QRA = computePercentID(query, parentA, breakpoint+1, length - 1);
 131
 132                         float LAB = computePercentID(parentA, parentB, 0, breakpoint);
 133                         float RAB = computePercentID(parentA, parentB, breakpoint+1, length - 1);
 134
 135                         float AB = ((LAB*leftLength) + (RAB*rightLength)) / (float) length;
 136                         float QA = ((QLA*leftLength) + (QRA*rightLength)) / (float) length;
 137                         float QB = ((QLB*leftLength) + (QRB*rightLength)) / (float) length;
 138
 139                         float QLA_QRB = ((QLA*leftLength) + (QRB*rightLength)) / (float) length;
 140                         float QLB_QRA = ((QLB*leftLength) + (QRA*rightLength)) / (float) length;
 141
 142                         //in original and not used
 143                         //float avgQA_QB = ((QA*leftLength) + (QB*rightLength)) / (float) length;
 144
 145                         float divR_QLA_QRB = min((QLA_QRB/QA), (QLA_QRB/QB));
 146
 147                         float divR_QLB_QRA = min((QLB_QRA/QA), (QLB_QRA/QB));
 148
 149                         //is one of them above the
 150                         if (divR_QLA_QRB >= divRThreshold || divR_QLB_QRA >= divRThreshold) {
 151
 152                                 if (((QLA_QRB > QA) && (QLA_QRB > QB) && (QLA >= parentFragmentThreshold) && (QRB >= parentFragmentThreshold))  ||
 153                                         ((QLB_QRA > QA) && (QLB_QRA > QB) && (QLB >=parentFragmentThreshold) && (QRA >= parentFragmentThreshold)))  {
 154
 155                                         data_struct member;
 156
 157                                         member.divr_qla_qrb = divR_QLA_QRB;
 158                                         member.divr_qlb_qra = divR_QLB_QRA;
 159                                         member.qla_qrb = QLA_QRB;
 160                                         member.qlb_qra = QLB_QRA;
 161                                         member.qla = QLA;
 162                                         member.qrb = QRB;
 163                                         member.ab = AB;
 164                                         member.qa = QA;
 165                                         member.qb = QB;
 166                                         member.lab = LAB;
 167                                         member.rab = RAB;
 168                                         member.qra = QRA;
 169                                         member.qlb = QLB;
 170                                         member.winLStart = spots[0];
 171                                         member.winLEnd = spots[breakpoint];  //so breakpoint reflects spot in alignment before filter
 172                                         member.winRStart = spots[breakpoint+1];
 173                                         member.winREnd = spots[length-1];
 174                                         member.querySeq = *(q);
 175                                         member.parentA = *(pA);
 176                                         member.parentB = *(pB);
 177                                         member.bsa = 0;
 178                                         member.bsb = 0;
 179                                         member.bsMax = 0;
 180                                         member.chimeraMax = 0;
 181
 182                                         data.push_back(member);
 183
 184                                 }//if
 185                         }//if
 186                 }//for
 187
 188                 return data;
 189
 190         }
 191         catch(exception& e) {
 192                 errorOut(e, "Slayer", "runBellerophon");
 193                 exit(1);
 194         }
 195 }
 196 /***********************************************************************/
 197 vector<snps> Slayer::getSNPS(string parentA, string query, string parentB, int left, int right) {
 198         try {
 199
 200                 vector<snps> data;
 201
 202                 for (int i = left; i <= right; i++) {
 203
 204                         char A = parentA[i];
 205                         char Q = query[i];
 206                         char B = parentB[i];
 207
 208                         if ((A != Q) || (B != Q)) {
 209                                 snps member;
 210                                 member.queryChar = Q;
 211                                 member.parentAChar = A;
 212                                 member.parentBChar = B;
 213
 214                                 data.push_back(member);
 215                         }
 216                 }
 217
 218                 return data;
 219
 220         }
 221         catch(exception& e) {
 222                 errorOut(e, "Slayer", "getSNPS");
 223                 exit(1);
 224         }
 225 }
 226 /***********************************************************************/
 227 void Slayer::bootstrapSNPS(vector<snps> left, vector<snps> right, float& BSA, float& BSB) {
 228         try {
 229
 230                 srand((unsigned)time( NULL ));
 231
 232                 int count_A = 0; // sceneario QLA,QRB supported
 233                 int count_B = 0; // sceneario QLB,QRA supported
 234
 235                 int numLeft = max(1, int(left.size()/10 +0.5));
 236                 int numRight = max(1, int(right.size()/10 + 0.5));
 237
 238                 for (int i = 0; i < iters; i++) {
 239                         //random sampling with replacement.
 240
 241                         vector<snps> selectedLeft;
 242
 243                         for (int j = 0; j < numLeft; j++) {
 244                                 int index = int(rand() % left.size());
 245                                 selectedLeft.push_back(left[index]);
 246                         }
 247
 248                         vector<snps> selectedRight;
 249                         for (int j = 0; j < numRight; j++) {
 250                                 int index = int(rand() % right.size());
 251                                 selectedRight.push_back(right[index]);
 252                         }
 253
 254                         /* A  ------------------------------------------
 255                         #       QLA                     QRA
 256                         # Q  ------------------------------------------
 257                         #                      |
 258                         #                      |
 259                         # Q  ------------------------------------------
 260                         #       QLB                     QRB
 261                         # B  ------------------------------------------ */
 262
 263
 264                         float QLA = snpQA(selectedLeft);
 265                         float QRA = snpQA(selectedRight);
 266
 267                         float QLB = snpQB(selectedLeft);
 268                         float QRB = snpQB(selectedRight);
 269
 270                         //in original - not used - not sure why?
 271                         //float ALB = snpAB(selectedLeft);
 272                         //float ARB = snpAB(selectedRight);
 273
 274                         if ((QLA > QLB) && (QRB > QRA)) {
 275                                 count_A++;
 276                         }
 277
 278                         if ((QLB > QLA) && (QRA > QRB)) {
 279                                 count_B++;
 280                         }
 281
 282                 }
 283
 284                 BSA = ((float) count_A / (float) iters) * 100;
 285                 BSB = ((float) count_B / (float) iters) * 100;
 286
 287         }
 288         catch(exception& e) {
 289                 errorOut(e, "Slayer", "bootstrapSNPS");
 290                 exit(1);
 291         }
 292 }
 293 /***********************************************************************/
 294 float Slayer::snpQA(vector<snps> data) {
 295         try {
 296
 297                 int numIdentical = 0;
 298
 299                 for (int i = 0; i < data.size(); i++) {
 300                         if (data[i].parentAChar == data[i].queryChar) {
 301                                 numIdentical++;
 302                         }
 303                 }
 304
 305                 float percentID = (numIdentical / (float) data.size()) * 100;
 306
 307                 return percentID;
 308         }
 309         catch(exception& e) {
 310                 errorOut(e, "Slayer", "snpQA");
 311                 exit(1);
 312         }
 313 }
 314 /***********************************************************************/
 315 float Slayer::snpQB(vector<snps> data) {
 316         try {
 317
 318                 int numIdentical = 0;
 319
 320                 for (int i = 0; i < data.size(); i++) {
 321                         if (data[i].parentBChar == data[i].queryChar) {
 322                                 numIdentical++;
 323                         }
 324                 }
 325
 326                 float percentID = (numIdentical / (float) data.size()) * 100;
 327
 328                 return percentID;
 329
 330         }
 331         catch(exception& e) {
 332                 errorOut(e, "Slayer", "snpQB");
 333                 exit(1);
 334         }
 335 }
 336 /***********************************************************************/
 337 float Slayer::snpAB(vector<snps> data) {
 338         try {
 339                 int numIdentical = 0;
 340
 341                 for (int i = 0; i < data.size(); i++) {
 342                         if (data[i].parentAChar == data[i].parentBChar) {
 343                                 numIdentical++;
 344                         }
 345                 }
 346
 347                 float percentID = (numIdentical / (float) data.size()) * 100;
 348
 349                 return percentID;
 350
 351         }
 352         catch(exception& e) {
 353                 errorOut(e, "Slayer", "snpAB");
 354                 exit(1);
 355         }
 356 }
 357 /***********************************************************************/
 358 float Slayer::computePercentID(string queryFrag, string parent, int left, int right) {
 359         try {
 360                 int total = 0;
 361                 int matches = 0;
 362
 363                 for (int i = left; i <= right; i++) {
 364                         total++;
 365                         if (queryFrag[i] == parent[i]) {
 366                                 matches++;
 367                         }
 368                 }
 369
 370                 float percentID =( matches/(float)total) * 100;
 371
 372                 return percentID;
 373         }
 374         catch(exception& e) {
 375                 errorOut(e, "Slayer", "computePercentID");
 376                 exit(1);
 377         }
 378 }
 379 /***********************************************************************/
 380 //remove columns that contain any gaps
 381 map<int, int> Slayer::verticalFilter(vector<Sequence*> seqs) {
 382         try {
 383                 vector<int> gaps;       gaps.resize(seqs[0]->getAligned().length(), 0);
 384
 385                 string filterString = (string(seqs[0]->getAligned().length(), '1'));
 386
 387                 //for each sequence
 388                 for (int i = 0; i < seqs.size(); i++) {
 389
 390                         string seqAligned = seqs[i]->getAligned();
 391
 392                         for (int j = 0; j < seqAligned.length(); j++) {
 393                                 //if this spot is a gap
 394                                 if ((seqAligned[j] == '-') || (seqAligned[j] == '.') || (toupper(seqAligned[j]) == 'N'))        {       gaps[j]++;      }
 395                         }
 396                 }
 397
 398                 //zero out spot where all sequences have blanks
 399                 int numColRemoved = 0;
 400                 int count = 0;
 401                 map<int, int> maskMap; maskMap.clear();
 402
 403                 for(int i = 0; i < seqs[0]->getAligned().length(); i++){
 404                         if(gaps[i] != 0)        {       filterString[i] = '0';  numColRemoved++;  }
 405                         else {
 406                                 maskMap[count] = i;
 407                                 count++;
 408                         }
 409                 }
 410
 411                 //for each sequence
 412                 for (int i = 0; i < seqs.size(); i++) {
 413
 414                         string seqAligned = seqs[i]->getAligned();
 415                         string newAligned = "";
 416
 417                         for (int j = 0; j < seqAligned.length(); j++) {
 418                                 //if this spot is not a gap
 419                                 if (filterString[j] == '1') { newAligned += seqAligned[j]; }
 420                         }
 421
 422                         seqs[i]->setAligned(newAligned);
 423                 }
 424
 425                 return maskMap;
 426         }
 427         catch(exception& e) {
 428                 errorOut(e, "Slayer", "verticalFilter");
 429                 exit(1);
 430         }
 431 }
 432 /***********************************************************************/