maligner.cpp

   1 /*
   2  *  maligner.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 9/23/09.
   6  *  Copyright 2009 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "maligner.h"
  11
  12 /***********************************************************************/
  13 Maligner::Maligner(vector<Sequence*> temp, int num, int match, int misMatch, float div, int minCov) :
  14                 db(temp), numWanted(num), matchScore(match), misMatchPenalty(misMatch), minDivR(div), minCoverage(minCov) {}
  15 /***********************************************************************/
  16 string Maligner::getResults(Sequence* q) {
  17         try {
  18
  19                 outputResults.clear();
  20
  21                 //make copy so trimming doesn't destroy query from calling class - remember to deallocate
  22                 query = new Sequence(q->getName(), q->getAligned());
  23
  24                 string chimera;
  25
  26                 decalc = new DeCalculator();
  27
  28                 //find closest seqs to query in template - returns copies of seqs so trim does not destroy - remember to deallocate
  29                 refSeqs = decalc->findClosest(query, db, numWanted);
  30
  31                 refSeqs = minCoverageFilter(refSeqs);
  32
  33                 if (refSeqs.size() < 2)  {
  34                         for (int i = 0; i < refSeqs.size(); i++) {  delete refSeqs[i];  }
  35                         percentIdenticalQueryChimera = 0.0;
  36                         return "unknown";
  37                 }
  38
  39                 int chimeraPenalty = computeChimeraPenalty();
  40
  41                 //trims seqs to first non gap char in all seqs and last non gap char in all seqs
  42                 decalc->trimSeqs(query, refSeqs);
  43
  44                 vector<Sequence*> temp = refSeqs;
  45                 temp.push_back(query);
  46
  47                 verticalFilter(temp);
  48
  49                 vector< vector<score_struct> > matrix = buildScoreMatrix(query->getAligned().length(), refSeqs.size()); //builds and initializes
  50
  51                 fillScoreMatrix(matrix, refSeqs, chimeraPenalty);
  52
  53                 vector<score_struct> path = extractHighestPath(matrix);
  54
  55                 vector<trace_struct> trace = mapTraceRegionsToAlignment(path, refSeqs);
  56
  57                 if (trace.size() > 1) {         chimera = "yes";        }
  58                 else { chimera = "no";  }
  59
  60                 int traceStart = path[0].col;
  61                 int traceEnd = path[path.size()-1].col;
  62
  63                 string queryInRange = query->getAligned();
  64                 queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart+1));
  65
  66                 string chimeraSeq = constructChimericSeq(trace, refSeqs);
  67
  68                 percentIdenticalQueryChimera = computePercentID(queryInRange, chimeraSeq);
  69
  70                 delete decalc;
  71
  72                 //save output results
  73                 for (int i = 0; i < trace.size(); i++) {
  74                         int regionStart = trace[i].col;
  75                         int regionEnd = trace[i].oldCol;
  76                         int seqIndex = trace[i].row;
  77
  78                         results temp;
  79
  80                         temp.parent = refSeqs[seqIndex]->getName();
  81                         temp.regionStart = regionStart;
  82                         temp.regionEnd = regionEnd;
  83
  84                         string parentInRange = refSeqs[seqIndex]->getAligned();
  85                         parentInRange = parentInRange.substr(traceStart, (traceEnd-traceStart+1));
  86
  87                         temp.queryToParent = computePercentID(queryInRange, parentInRange);
  88                         temp.divR = (percentIdenticalQueryChimera / temp.queryToParent);
  89
  90                         string queryInRegion = query->getAligned();
  91                         queryInRegion = queryInRegion.substr(regionStart, (regionEnd-regionStart+1));
  92
  93                         string parentInRegion = refSeqs[seqIndex]->getAligned();
  94                         parentInRegion = parentInRegion.substr(regionStart, (regionEnd-regionStart+1));
  95
  96                         temp.queryToParentLocal = computePercentID(queryInRegion, parentInRegion);
  97
  98                         outputResults.push_back(temp);
  99                 }
 100
 101                 //free memory
 102                 delete query;
 103                 for (int i = 0; i < refSeqs.size(); i++) {  delete refSeqs[i];  }
 104
 105                 return chimera;
 106         }
 107         catch(exception& e) {
 108                 errorOut(e, "Maligner", "getResults");
 109                 exit(1);
 110         }
 111 }
 112 /***********************************************************************/
 113 //removes top matches that do not have minimum coverage with query.
 114 vector<Sequence*> Maligner::minCoverageFilter(vector<Sequence*> ref){
 115         try {
 116                 vector<Sequence*> newRefs;
 117
 118                 string queryAligned = query->getAligned();
 119
 120                 for (int i = 0; i < ref.size(); i++) {
 121
 122                         string refAligned = ref[i]->getAligned();
 123
 124                         int numBases = 0;
 125                         int numCovered = 0;
 126
 127                         //calculate coverage
 128                         for (int j = 0; j < queryAligned.length(); j++) {
 129
 130                                 if (isalpha(queryAligned[j])) {
 131                                         numBases++;
 132
 133                                         if (isalpha(refAligned[j])) {
 134                                                 numCovered++;
 135                                         }
 136                                 }
 137                         }
 138
 139                         int coverage = ((numCovered/(float)numBases)*100);
 140
 141                         //if coverage above minimum
 142                         if (coverage > minCoverage) {
 143                                 newRefs.push_back(ref[i]);
 144                         }
 145                 }
 146
 147                 return newRefs;
 148         }
 149         catch(exception& e) {
 150                 errorOut(e, "Maligner", "minCoverageFilter");
 151                 exit(1);
 152         }
 153 }
 154 /***********************************************************************/
 155 // a breakpoint should yield fewer mismatches than this number with respect to the best parent sequence.
 156 int Maligner::computeChimeraPenalty() {
 157         try {
 158
 159                 int numAllowable = ((1.0 - (1.0/minDivR)) * query->getNumBases());
 160
 161                 int penalty = int(numAllowable + 1) * misMatchPenalty;
 162
 163                 return penalty;
 164
 165         }
 166         catch(exception& e) {
 167                 errorOut(e, "Maligner", "computeChimeraPenalty");
 168                 exit(1);
 169         }
 170 }
 171 /***********************************************************************/
 172 //this is a vertical filter
 173 void Maligner::verticalFilter(vector<Sequence*> seqs) {
 174         try {
 175                 vector<int> gaps;       gaps.resize(query->getAligned().length(), 0);
 176
 177                 string filterString = (string(query->getAligned().length(), '1'));
 178
 179                 //for each sequence
 180                 for (int i = 0; i < seqs.size(); i++) {
 181
 182                         string seqAligned = seqs[i]->getAligned();
 183
 184                         for (int j = 0; j < seqAligned.length(); j++) {
 185                                 //if this spot is a gap
 186                                 if ((seqAligned[j] == '-') || (seqAligned[j] == '.'))   {       gaps[j]++;      }
 187                         }
 188                 }
 189
 190                 //zero out spot where all sequences have blanks
 191                 int numColRemoved = 0;
 192                 for(int i = 0; i < seqs[0]->getAligned().length(); i++){
 193                         if(gaps[i] == seqs.size())      {       filterString[i] = '0';  numColRemoved++;  }
 194                 }
 195
 196                 //for each sequence
 197                 for (int i = 0; i < seqs.size(); i++) {
 198
 199                         string seqAligned = seqs[i]->getAligned();
 200                         string newAligned = "";
 201
 202                         for (int j = 0; j < seqAligned.length(); j++) {
 203                                 //if this spot is not a gap
 204                                 if (filterString[j] == '1') { newAligned += seqAligned[j]; }
 205                         }
 206
 207                         seqs[i]->setAligned(newAligned);
 208                 }
 209
 210
 211         }
 212         catch(exception& e) {
 213                 errorOut(e, "Maligner", "verticalFilter");
 214                 exit(1);
 215         }
 216 }
 217 //***************************************************************************************************************
 218 vector< vector<score_struct> > Maligner::buildScoreMatrix(int cols, int rows) {
 219         try{
 220
 221                 vector< vector<score_struct> > m; m.resize(rows);
 222
 223                 for (int i = 0; i < m.size(); i++) {
 224                         for (int j = 0; j < cols; j++) {
 225
 226                                 //initialize each cell
 227                                 score_struct temp;
 228                                 temp.prev = -1;
 229                                 temp.score = -9999999;
 230                                 temp.col = j;
 231                                 temp.row = i;
 232
 233                                 m[i].push_back(temp);
 234                         }
 235                 }
 236
 237                 return m;
 238         }
 239         catch(exception& e) {
 240                 errorOut(e, "Maligner", "buildScoreMatrix");
 241                 exit(1);
 242         }
 243 }
 244 //***************************************************************************************************************
 245 void Maligner::fillScoreMatrix(vector<vector<score_struct> >& m, vector<Sequence*> seqs, int penalty) {
 246         try{
 247
 248                 //get matrix dimensions
 249                 int numCols = query->getAligned().length();
 250                 int numRows = seqs.size();
 251
 252                 //initialize first col
 253                 string queryAligned = query->getAligned();
 254                 for (int i = 0; i < numRows; i++) {
 255                         string subjectAligned = seqs[i]->getAligned();
 256
 257                         //are you both gaps?
 258                         if ((!isalpha(queryAligned[0])) && (!isalpha(subjectAligned[0]))) {
 259                                 m[i][0].score = 0;
 260                         }else if (queryAligned[0] == subjectAligned[0]) {
 261                                 m[i][0].score = matchScore;
 262                         }else{
 263                                 m[i][0].score = 0;
 264                         }
 265                 }
 266
 267                 //fill rest of matrix
 268                 for (int j = 1; j < numCols; j++) {  //iterate through matrix columns
 269
 270                         for (int i = 0; i < numRows; i++) {  //iterate through matrix rows
 271
 272                                 string subjectAligned = seqs[i]->getAligned();
 273
 274                                 int matchMisMatchScore = 0;
 275                                 //are you both gaps?
 276                                 if ((!isalpha(queryAligned[j])) && (!isalpha(subjectAligned[j]))) {
 277                                         //leave the same
 278                                 }else if ((toupper(queryAligned[j]) == 'N') || (toupper(subjectAligned[j]) == 'N')) {
 279                                         //leave the same
 280                                 }else if (queryAligned[j] == subjectAligned[j]) {
 281                                         matchMisMatchScore = matchScore;
 282                                 }else if (queryAligned[j] != subjectAligned[j]) {
 283                                         matchMisMatchScore = misMatchPenalty;
 284                                 }
 285
 286                                 //compute score based on previous columns scores
 287                                 for (int prevIndex = 0; prevIndex < numRows; prevIndex++) { //iterate through rows
 288
 289                                         int sumScore = matchMisMatchScore + m[prevIndex][j-1].score;
 290
 291                                         //you are not at yourself
 292                                         if (prevIndex != i) {   sumScore += penalty;    }
 293                                         if (sumScore < 0)       {       sumScore = 0;                   }
 294
 295                                         if (sumScore > m[i][j].score) {
 296                                                 m[i][j].score = sumScore;
 297                                                 m[i][j].prev = prevIndex;
 298                                         }
 299                                 }
 300                         }
 301                 }
 302
 303         }
 304         catch(exception& e) {
 305                 errorOut(e, "Maligner", "fillScoreMatrix");
 306                 exit(1);
 307         }
 308 }
 309 //***************************************************************************************************************
 310 vector<score_struct> Maligner::extractHighestPath(vector<vector<score_struct> > m) {
 311         try {
 312
 313                 //get matrix dimensions
 314                 int numCols = query->getAligned().length();
 315                 int numRows = m.size();
 316
 317
 318                 //find highest score scoring matrix
 319                 score_struct highestStruct;
 320                 int highestScore = 0;
 321
 322                 for (int i = 0; i < numRows; i++) {
 323                         for (int j = 0; j < numCols; j++) {
 324                                 if (m[i][j].score > highestScore) {
 325                                         highestScore = m[i][j].score;
 326                                         highestStruct = m[i][j];
 327                                 }
 328                         }
 329                 }
 330
 331                 vector<score_struct> path;
 332
 333                 int rowIndex = highestStruct.row;
 334                 int pos = highestStruct.col;
 335                 int score = highestStruct.score;
 336
 337                 while (pos >= 0 && score > 0) {
 338                         score_struct temp = m[rowIndex][pos];
 339                         score = temp.score;
 340
 341                         if (score > 0) {        path.push_back(temp);   }
 342
 343                         rowIndex = temp.prev;
 344                         pos--;
 345                 }
 346
 347                 reverse(path.begin(), path.end());
 348
 349                 return path;
 350
 351         }
 352         catch(exception& e) {
 353                 errorOut(e, "Maligner", "extractHighestPath");
 354                 exit(1);
 355         }
 356 }
 357 //***************************************************************************************************************
 358 vector<trace_struct> Maligner::mapTraceRegionsToAlignment(vector<score_struct> path, vector<Sequence*> seqs) {
 359         try {
 360                 vector<trace_struct> trace;
 361
 362                 int region_index = path[0].row;
 363                 int region_start = path[0].col;
 364
 365                 for (int i = 1; i < path.size(); i++) {
 366
 367                         int next_region_index = path[i].row;
 368
 369                         if (next_region_index != region_index) {
 370
 371                                 // add trace region
 372                                 int col_index = path[i].col;
 373                                 trace_struct temp;
 374                                 temp.col = region_start;
 375                                 temp.oldCol = col_index-1;
 376                                 temp.row = region_index;
 377
 378                                 trace.push_back(temp);
 379
 380                                 region_index = path[i].row;
 381                                 region_start = col_index;
 382                         }
 383                 }
 384
 385                 // get last one
 386                 trace_struct temp;
 387                 temp.col = region_start;
 388                 temp.oldCol = path[path.size()-1].col;
 389                 temp.row = region_index;
 390                 trace.push_back(temp);
 391
 392                 return trace;
 393
 394         }
 395         catch(exception& e) {
 396                 errorOut(e, "Maligner", "mapTraceRegionsToAlignment");
 397                 exit(1);
 398         }
 399 }
 400 //***************************************************************************************************************
 401 string Maligner::constructChimericSeq(vector<trace_struct> trace, vector<Sequence*> seqs) {
 402         try {
 403                 string chimera = "";
 404
 405                 for (int i = 0; i < trace.size(); i++) {
 406                         string seqAlign = seqs[trace[i].row]->getAligned();
 407                         seqAlign = seqAlign.substr(trace[i].col, (trace[i].oldCol-trace[i].col+1));
 408                         chimera += seqAlign;
 409                 }
 410
 411                 return chimera;
 412         }
 413         catch(exception& e) {
 414                 errorOut(e, "Maligner", "constructChimericSeq");
 415                 exit(1);
 416         }
 417 }
 418 //***************************************************************************************************************
 419 float Maligner::computePercentID(string queryAlign, string chimera) {
 420         try {
 421
 422                 if (queryAlign.length() != chimera.length()) {
 423                         mothurOut("Error, alignment strings are of different lengths: "); mothurOutEndLine();
 424                         mothurOut(toString(queryAlign.length())); mothurOutEndLine(); mothurOutEndLine(); mothurOutEndLine(); mothurOutEndLine();
 425                         mothurOut(toString(chimera.length())); mothurOutEndLine();
 426                         return -1.0;
 427                 }
 428
 429
 430                 int numBases = 0;
 431                 int numIdentical = 0;
 432
 433                 for (int i = 0; i < queryAlign.length(); i++) {
 434                         if ((isalpha(queryAlign[i])) || (isalpha(chimera[i])))  {
 435                                 numBases++;
 436                                 if (queryAlign[i] == chimera[i]) {
 437                                         numIdentical++;
 438                                 }
 439                         }
 440                 }
 441
 442                 if (numBases == 0) { return 0; }
 443
 444                 float percentIdentical = (numIdentical/(float)numBases) * 100;
 445
 446                 return percentIdentical;
 447
 448         }
 449         catch(exception& e) {
 450                 errorOut(e, "Maligner", "computePercentID");
 451                 exit(1);
 452         }
 453 }
 454 //***************************************************************************************************************
 455