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(vector<Sequence*> t, int m, int mm) : match(m), misMatch(mm) {  templateSeqs = t;  }
  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                         //take query and break apart into pieces using breakpoints given by results of parents
  35                         for (int i = 0; i < parents.size(); i++) {
  36         cout << parents[i].parent << '\t' << parents[i].nastRegionStart <<  '\t' << parents[i].nastRegionEnd    << endl;
  37                                 int length = parents[i].nastRegionEnd - parents[i].nastRegionStart;
  38                                 string q = qAligned.substr(parents[i].nastRegionStart, length);
  39         cout << "query = " << q << endl;
  40                                 Sequence* queryFrag = new Sequence(query->getName(), q);
  41
  42                                 queryParts.push_back(queryFrag);
  43
  44                                 Sequence* parent = getSequence(parents[i].parent);
  45                                 string p = parent->getAligned();
  46
  47                                 p = p.substr(parents[i].nastRegionStart, length);
  48         cout << "parent = " << p << endl;
  49                                 parent->setAligned(p);
  50
  51                                 parentParts.push_back(parent);
  52
  53                                 if (q.length() > longest)       { longest = q.length(); }
  54                                 if (p.length() > longest)       { longest = p.length(); }
  55                         }
  56
  57                         //align each peice to correct parent from results
  58                         for (int i = 0; i < queryParts.size(); i++) {
  59                                 alignment = new NeedlemanOverlap(-2.0, match, misMatch, longest+1); //default gapopen, match, mismatch, longestbase
  60                                 Nast nast(alignment, queryParts[i], parentParts[i]);
  61                                 delete alignment;
  62                         }
  63
  64                         //recombine pieces to form new query sequence
  65                         for (int i = 0; i < queryParts.size(); i++) {
  66                                 newQuery += queryParts[i]->getAligned();
  67                         }
  68
  69
  70                         //make sure you don't cutoff end of query
  71                         if (parents[parents.size()-1].nastRegionEnd < qAligned.length()) {  newQuery += qAligned.substr(parents[parents.size()-1].nastRegionEnd);  }
  72
  73                         //set query to new aligned string
  74                         query->setAligned(newQuery);
  75
  76                         //free memory
  77                         for (int i = 0; i < queryParts.size(); i++) { delete queryParts[i];  }
  78                         for (int i = 0; i < parentParts.size(); i++) { delete parentParts[i];  }
  79
  80                 } //else leave query alone, you have bigger problems...
  81
  82         }
  83         catch(exception& e) {
  84                 errorOut(e, "ChimeraReAligner", "reAlign");
  85                 exit(1);
  86         }
  87 }
  88 //***************************************************************************************************************
  89 Sequence* ChimeraReAligner::getSequence(string name) {
  90         try{
  91                 Sequence* temp;
  92
  93                 //look through templateSeqs til you find it
  94                 int spot = -1;
  95                 for (int i = 0; i < templateSeqs.size(); i++) {
  96                         if (name == templateSeqs[i]->getName()) {
  97                                 spot = i;
  98                                 break;
  99                         }
 100                 }
 101
 102                 if(spot == -1) { mothurOut("Error: Could not find sequence."); mothurOutEndLine(); return NULL; }
 103
 104                 temp = new Sequence(templateSeqs[spot]->getName(), templateSeqs[spot]->getAligned());
 105
 106                 return temp;
 107         }
 108         catch(exception& e) {
 109                 errorOut(e, "ChimeraReAligner", "getSequence");
 110                 exit(1);
 111         }
 112 }
 113 //***************************************************************************************************************