5 * Created by westcott on 2/12/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "chimerarealigner.h"
11 #include "needlemanoverlap.hpp"
14 //***************************************************************************************************************
15 ChimeraReAligner::ChimeraReAligner() { m = MothurOut::getInstance(); }
16 //***************************************************************************************************************
17 ChimeraReAligner::~ChimeraReAligner() {}
18 //***************************************************************************************************************
19 void ChimeraReAligner::reAlign(Sequence* query, vector<results> parents) {
21 if (parents.size() != 0) {
22 vector<Sequence*> queryParts;
23 vector<Sequence*> parentParts; //queryParts[0] relates to parentParts[0]
25 string qAligned = query->getAligned();
28 //sort parents by region start
29 sort(parents.begin(), parents.end(), compareRegionStart);
31 //make sure you don't cutoff beginning of query
32 if (parents[0].nastRegionStart > 0) { newQuery += qAligned.substr(0, parents[0].nastRegionStart); }
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);
40 Sequence* queryFrag = new Sequence(query->getName(), q);
41 queryParts.push_back(queryFrag);
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);
48 if (queryFrag->getUnaligned().length() > longest) { longest = queryFrag->getUnaligned().length(); }
49 if (parent->getUnaligned().length() > longest) { longest = parent->getUnaligned().length(); }
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() == "")) {;}
56 Alignment* alignment = new NeedlemanOverlap(-2.0, 1.0, -1.0, longest+1); //default gapopen, match, mismatch, longestbase
58 Nast nast(alignment, queryParts[i], parentParts[i]);
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
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);
75 newQuery += queryParts[i]->getAligned();
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); }
81 query->setAligned(newQuery);
84 for (int i = 0; i < queryParts.size(); i++) { delete queryParts[i]; }
85 for (int i = 0; i < parentParts.size(); i++) { delete parentParts[i]; }
87 } //else leave query alone, you have bigger problems...
91 m->errorOut(e, "ChimeraReAligner", "reAlign");
95 /***************************************************************************************************************
96 Sequence* ChimeraReAligner::getSequence(string name) {
100 //look through templateSeqs til you find it
102 for (int i = 0; i < templateSeqs.size(); i++) {
103 if (name == templateSeqs[i]->getName()) {
109 if(spot == -1) { m->mothurOut("Error: Could not find sequence."); m->mothurOutEndLine(); return NULL; }
111 temp = new Sequence(templateSeqs[spot]->getName(), templateSeqs[spot]->getAligned());
115 catch(exception& e) {
116 m->errorOut(e, "ChimeraReAligner", "getSequence");
120 //***************************************************************************************************************/