chimerarealigner.cpp

   1 /*
   2  *  chimerarealigner.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 2/12/10.
   6  *  Copyright 2010 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "chimerarealigner.h"
  11 #include "needlemanoverlap.hpp"
  12 #include "nast.hpp"
  13
  14 //***************************************************************************************************************
  15 ChimeraReAligner::ChimeraReAligner()  {  m = MothurOut::getInstance(); }
  16 //***************************************************************************************************************
  17 ChimeraReAligner::~ChimeraReAligner() {}
  18 //***************************************************************************************************************
  19 void ChimeraReAligner::reAlign(Sequence* query, vector<results> parents) {
  20         try {
  21                 if (parents.size() != 0) {
  22                         vector<Sequence*> queryParts;
  23                         vector<Sequence*> parentParts;  //queryParts[0] relates to parentParts[0]
  24
  25                         string qAligned = query->getAligned();
  26                         string newQuery = "";
  27
  28                         //sort parents by region start
  29                         sort(parents.begin(), parents.end(), compareRegionStart);
  30
  31                         //make sure you don't cutoff beginning of query
  32                         if (parents[0].nastRegionStart > 0) {  newQuery += qAligned.substr(0, parents[0].nastRegionStart);  }
  33                         int longest = 0;
  34
  35                         //take query and break apart into pieces using breakpoints given by results of parents
  36                         for (int i = 0; i < parents.size(); i++) {
  37                                 int length = parents[i].nastRegionEnd - parents[i].nastRegionStart+1;
  38                                 string q = qAligned.substr(parents[i].nastRegionStart, length);
  39
  40                                 Sequence* queryFrag = new Sequence(query->getName(), q);
  41                                 queryParts.push_back(queryFrag);
  42
  43                                 string p = parents[i].parentAligned;
  44                                 p = p.substr(parents[i].nastRegionStart, length);
  45                                 Sequence* parent = new Sequence(parents[i].parent, p);
  46                                 parentParts.push_back(parent);
  47
  48                                 if (queryFrag->getUnaligned().length() > longest)       { longest = queryFrag->getUnaligned().length(); }
  49                                 if (parent->getUnaligned().length() > longest)  { longest = parent->getUnaligned().length();    }
  50                         }
  51
  52                         //align each peice to correct parent from results
  53                         for (int i = 0; i < queryParts.size(); i++) {
  54                                 if ((queryParts[i]->getUnaligned() == "") || (parentParts[i]->getUnaligned() == "")) {;}
  55                                 else {
  56                                         Alignment* alignment = new NeedlemanOverlap(-2.0, 1.0, -1.0, longest+1); //default gapopen, match, mismatch, longestbase
  57
  58                                         Nast nast(alignment, queryParts[i], parentParts[i]);
  59                                         delete alignment;
  60                                 }
  61                         }
  62
  63                         //recombine pieces to form new query sequence
  64                         for (int i = 0; i < queryParts.size(); i++) {
  65                                 //sometimes the parent regions do not meet, for example region 1 may end at 1000 and region 2 starts at 1100.
  66                                 //we don't want to loose length so in this case we will leave query alone
  67                                 if (i != 0) {
  68                                         int space = parents[i].nastRegionStart - parents[i-1].nastRegionEnd - 1;
  69                                         if (space > 0) { //they don't meet and we need to add query piece
  70                                                 string q = qAligned.substr(parents[i-1].nastRegionEnd+1, space);
  71                                                 newQuery += q;
  72                                         }
  73                                 }
  74
  75                                 newQuery += queryParts[i]->getAligned();
  76                         }
  77
  78                         //make sure you don't cutoff end of query
  79                         if (parents[parents.size()-1].nastRegionEnd < (qAligned.length()-1)) {  newQuery += qAligned.substr(parents[parents.size()-1].nastRegionEnd+1);  }
  80
  81                         query->setAligned(newQuery);
  82
  83                         //free memory
  84                         for (int i = 0; i < queryParts.size(); i++) { delete queryParts[i];  }
  85                         for (int i = 0; i < parentParts.size(); i++) { delete parentParts[i];  }
  86
  87                 } //else leave query alone, you have bigger problems...
  88
  89         }
  90         catch(exception& e) {
  91                 m->errorOut(e, "ChimeraReAligner", "reAlign");
  92                 exit(1);
  93         }
  94 }
  95 /***************************************************************************************************************
  96 Sequence* ChimeraReAligner::getSequence(string name) {
  97         try{
  98                 Sequence* temp;
  99
 100                 //look through templateSeqs til you find it
 101                 int spot = -1;
 102                 for (int i = 0; i < templateSeqs.size(); i++) {
 103                         if (name == templateSeqs[i]->getName()) {
 104                                 spot = i;
 105                                 break;
 106                         }
 107                 }
 108
 109                 if(spot == -1) { m->mothurOut("Error: Could not find sequence."); m->mothurOutEndLine(); return NULL; }
 110
 111                 temp = new Sequence(templateSeqs[spot]->getName(), templateSeqs[spot]->getAligned());
 112
 113                 return temp;
 114         }
 115         catch(exception& e) {
 116                 m->errorOut(e, "ChimeraReAligner", "getSequence");
 117                 exit(1);
 118         }
 119 }
 120 //***************************************************************************************************************/