5 * Created by westcott on 9/25/09.
6 * Copyright 2009 Pschloss Lab. All rights reserved.
10 #include "chimeraslayer.h"
11 #include "chimerarealigner.h"
13 //***************************************************************************************************************
14 ChimeraSlayer::ChimeraSlayer(string mode, bool r) : searchMethod(mode), realign(r) { decalc = new DeCalculator(); }
15 //***************************************************************************************************************
16 ChimeraSlayer::~ChimeraSlayer() { delete decalc; }
17 //***************************************************************************************************************
18 void ChimeraSlayer::printHeader(ostream& out) {
20 mothurOut("Only reporting sequence supported by 90% of bootstrapped results.");
23 //***************************************************************************************************************
24 void ChimeraSlayer::print(ostream& out) {
26 if (chimeraFlags == "yes") {
27 string chimeraFlag = "no";
28 if( (chimeraResults[0].bsa >= minBS && chimeraResults[0].divr_qla_qrb >= divR)
30 (chimeraResults[0].bsb >= minBS && chimeraResults[0].divr_qlb_qra >= divR) ) { chimeraFlag = "yes"; }
33 if (chimeraFlag == "yes") {
34 if ((chimeraResults[0].bsa >= minBS) || (chimeraResults[0].bsb >= minBS)) {
35 mothurOut(querySeq->getName() + "\tyes"); mothurOutEndLine();
38 out << querySeq->getName() << "\tyes" << endl;
39 printBlock(chimeraResults[0], out);
41 }else { out << querySeq->getName() << "\tno" << endl; }
45 errorOut(e, "ChimeraSlayer", "print");
49 //***************************************************************************************************************
50 int ChimeraSlayer::getChimeras(Sequence* query) {
55 for (int i = 0; i < query->getAligned().length(); i++) {
59 //referenceSeqs, numWanted, matchScore, misMatchPenalty, divR, minSimilarity
60 maligner = new Maligner(templateSeqs, numWanted, match, misMatch, divR, minSim, minCov, searchMethod);
61 slayer = new Slayer(window, increment, minSim, divR, iters, minSNP);
63 string chimeraFlag = maligner->getResults(query, decalc);
64 vector<results> Results = maligner->getOutput();
66 //realign query to parents to improve slayers detection rate???
68 ChimeraReAligner realigner(templateSeqs, match, misMatch);
69 realigner.reAlign(query, Results);
72 //if (chimeraFlag == "yes") {
74 //get sequence that were given from maligner results
76 map<string, float> removeDups;
77 map<string, float>::iterator itDup;
78 for (int j = 0; j < Results.size(); j++) {
79 float dist = (Results[j].regionEnd - Results[j].regionStart + 1) * Results[j].queryToParentLocal;
80 //only add if you are not a duplicate
81 itDup = removeDups.find(Results[j].parent);
82 if (itDup == removeDups.end()) { //this is not duplicate
83 removeDups[Results[j].parent] = dist;
84 }else if (dist > itDup->second) { //is this a stronger number for this parent
85 removeDups[Results[j].parent] = dist;
89 for (itDup = removeDups.begin(); itDup != removeDups.end(); itDup++) {
90 Sequence* seq = getSequence(itDup->first); //makes copy so you can filter and mask and not effect template
94 member.dist = itDup->second;
96 seqs.push_back(member);
99 //limit number of parents to explore - default 3
100 if (Results.size() > parents) {
102 sort(seqs.begin(), seqs.end(), compareSeqDist);
103 //prioritize larger more similiar sequence fragments
104 reverse(seqs.begin(), seqs.end());
106 for (int k = seqs.size()-1; k > (parents-1); k--) {
112 //put seqs into vector to send to slayer
113 vector<Sequence*> seqsForSlayer;
114 for (int k = 0; k < seqs.size(); k++) { seqsForSlayer.push_back(seqs[k].seq); }
115 //cout << i+1 << "num parents = " << seqsForSlayer.size() << '\t' << chimeraFlag << endl;
117 //string name = querySeqs[i]->getName();
118 //cout << name << endl;
119 //name = name.substr(name.find_first_of("|")+1);
120 //cout << name << endl;
121 //name = name.substr(name.find_first_of("|")+1);
122 //cout << name << endl;
123 //name = name.substr(0, name.find_last_of("|"));
124 //cout << name << endl;
125 //string filename = toString(i+1) + ".seqsforslayer";
126 //openOutputFile(filename, out);
127 //cout << querySeqs[i]->getName() << endl;
128 //for (int u = 0; u < seqsForSlayer.size(); u++) { cout << seqsForSlayer[u]->getName() << '\t'; seqsForSlayer[u]->printSequence(out); }
131 //filename = toString(i+1) + ".fasta";
132 //openOutputFile(filename, out);
133 //querySeqs[i]->printSequence(out);
138 //mask then send to slayer...
140 decalc->setMask(seqMask);
143 decalc->runMask(query);
146 for (int k = 0; k < seqsForSlayer.size(); k++) {
147 decalc->runMask(seqsForSlayer[k]);
150 spotMap = decalc->getMaskMap();
154 chimeraFlags = slayer->getResults(query, seqsForSlayer);
155 chimeraResults = slayer->getOutput();
158 for (int k = 0; k < seqs.size(); k++) { delete seqs[k].seq; }
163 catch(exception& e) {
164 errorOut(e, "ChimeraSlayer", "getChimeras");
168 //***************************************************************************************************************
169 void ChimeraSlayer::printBlock(data_struct data, ostream& out){
172 out << "parentA: " << data.parentA.getName() << " parentB: " << data.parentB.getName() << endl;
173 out << "Left Window: " << spotMap[data.winLStart] << " " << spotMap[data.winLEnd] << endl;
174 out << "Right Window: " << spotMap[data.winRStart] << " " << spotMap[data.winREnd] << endl;
176 out << "Divergence of Query to Leftside ParentA and Rightside ParentB: " << data.divr_qla_qrb << '\t' << "Bootstrap: " << data.bsa << endl;
177 out << "Divergence of Query to Rightside ParentA and Leftside ParentB: " << data.divr_qlb_qra << '\t' << "Bootstrap: " << data.bsb << endl;
179 out << "Similarity of parents: " << data.ab << endl;
180 out << "Similarity of query to parentA: " << data.qa << endl;
181 out << "Similarity of query to parentB: " << data.qb << endl;
183 out << "Percent_ID QLA_QRB: " << data.qla_qrb << endl;
184 out << "Percent_ID QLB_QRA: " << data.qlb_qra << endl;
185 //out << "Per_id(QL,A): " << data.qla << endl;
186 //out << "Per_id(QL,B): " << data.qlb << endl;
187 //out << "Per_id(QR,A): " << data.qra << endl;
188 //out << "Per_id(QR,B): " << data.qrb << endl;
191 out << "DeltaL: " << (data.qla - data.qlb) << endl;
192 out << "DeltaR: " << (data.qra - data.qrb) << endl;
195 catch(exception& e) {
196 errorOut(e, "ChimeraSlayer", "printBlock");
200 //***************************************************************************************************************