5 * Created by westcott on 9/8/09.
6 * Copyright 2009 Schloss Lab. All rights reserved.
10 #include "chimeracheckrdp.h"
12 //***************************************************************************************************************
13 ChimeraCheckRDP::ChimeraCheckRDP(string filename, string temp) { fastafile = filename; templateFile = temp; }
14 //***************************************************************************************************************
16 ChimeraCheckRDP::~ChimeraCheckRDP() {
18 for (int i = 0; i < querySeqs.size(); i++) { delete querySeqs[i]; }
23 errorOut(e, "ChimeraCheckRDP", "~AlignSim");
27 //***************************************************************************************************************
28 void ChimeraCheckRDP::print(ostream& out) {
33 //vector<bool> isChimeric; isChimeric.resize(querySeqs.size(), false);
35 for (int i = 0; i < querySeqs.size(); i++) {
37 out << querySeqs[i]->getName() << endl;
38 out << "IS scores: " << '\t';
40 //int lastChimericWindowFound = 0;
42 for (int k = 0; k < IS[i].size(); k++) {
43 out << IS[i][k].score << '\t';
44 //if (IS[i][k].score > chimeraCutoff) { isChimeric[i] = true; lastChimericWindowFound = k; }
47 //if (isChimeric[i]) {
48 //mothurOut(querySeqs[i]->getName() + "\tIS: " + toString(IS[i][lastChimericWindowFound].score) + "\tbreakpoint: " + toString(IS[i][lastChimericWindowFound].midpoint) + "\tleft parent: " + IS[i][lastChimericWindowFound].leftParent + "\tright parent: " + IS[i][lastChimericWindowFound].rightParent); mothurOutEndLine();
49 //out << endl << "chimera: YES" << endl;
51 //out << endl << "chimera: NO" << endl;
56 if (name != "") { //if user has specific names
57 map<string, string>::iterator it = names.find(querySeqs[i]->getName());
59 if (it != names.end()) { //user wants pic of this
60 makeSVGpic(IS[i], i); //zeros out negative results
62 }else{//output them all
63 makeSVGpic(IS[i], i); //zeros out negative results
70 errorOut(e, "ChimeraCheckRDP", "print");
75 //***************************************************************************************************************
76 void ChimeraCheckRDP::getChimeras() {
79 //read in query sequences and subject sequences
81 mothurOut("Reading query sequences... "); cout.flush();
82 querySeqs = readSeqs(fastafile);
84 //templateSeqs = readSeqs(templateFile);
85 templateDB = new AlignmentDB(templateFile, "kmer", kmerSize, 0.0,0.0,0.0,0.0);
88 int numSeqs = querySeqs.size();
91 closest.resize(numSeqs);
93 //break up file if needed
94 int linesPerProcess = numSeqs / processors ;
96 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
97 //find breakup of sequences for all times we will Parallelize
98 if (processors == 1) { lines.push_back(new linePair(0, numSeqs)); }
101 for (int i = 0; i < (processors-1); i++) {
102 lines.push_back(new linePair((i*linesPerProcess), ((i*linesPerProcess) + linesPerProcess)));
104 //this is necessary to get remainder of processors / numSeqs so you don't miss any lines at the end
105 int i = processors - 1;
106 lines.push_back(new linePair((i*linesPerProcess), numSeqs));
110 lines.push_back(new linePair(0, numSeqs));
113 kmer = new Kmer(kmerSize);
116 readName(name); //fills name map with names of seqs the user wants to have .svg for.
119 //find closest seq to each querySeq
120 for (int i = 0; i < querySeqs.size(); i++) {
121 closest[i] = templateDB->findClosestSequence(querySeqs[i]);
124 //for each query find IS value
125 if (processors == 1) {
126 for (int i = 0; i < querySeqs.size(); i++) {
129 }else { createProcessesIS(); }
131 //determine chimera report cutoff - window score above 95%
132 //getCutoff(); - not very acurate predictor
135 for (int i = 0; i < lines.size(); i++) { delete lines[i]; }
139 catch(exception& e) {
140 errorOut(e, "ChimeraCheckRDP", "getChimeras");
144 //***************************************************************************************************************
145 vector<sim> ChimeraCheckRDP::findIS(int query) {
149 vector< map<int, int> > queryKmerInfo; //vector of maps - each entry in the vector is a map of the kmers up to that spot in the unaligned seq
150 //example: seqKmerInfo[50] = map containing the kmers found in the first 50 + kmersize characters of ecoli.
151 //i chose to store the kmers numbers in a map so you wouldn't have to check for dupilcate entries and could easily find the
152 //kmers 2 seqs had in common. There may be a better way to do this thats why I am leaving so many comments...
153 vector< map<int, int> > subjectKmerInfo;
155 vector<sim> isValues;
156 string queryName = querySeqs[query]->getName();
157 string seq = querySeqs[query]->getUnaligned();
159 mothurOut("Finding IS values for sequence " + toString(query+1)); mothurOutEndLine();
161 queryKmerInfo = kmer->getKmerCounts(seq);
162 subjectKmerInfo = kmer->getKmerCounts(closest[query].getUnaligned());
164 //find total kmers you have in common with closest[query] by looking at the last entry in the vector of maps for each
165 int nTotal = calcKmers(queryKmerInfo[(queryKmerInfo.size()-1)], subjectKmerInfo[(subjectKmerInfo.size()-1)]);
167 //you don't want the starting point to be virtually at hte end so move it in 10%
168 int start = seq.length() / 10;
171 for (int m = start; m < (seq.length() - start); m+=increment) {
173 if ((m - kmerSize) < 0) { mothurOut("Your sequence is too short for your kmerSize."); mothurOutEndLine(); exit(1); }
177 string fragLeft = seq.substr(0, m); //left side of breakpoint
178 string fragRight = seq.substr(m); //right side of breakpoint
180 //make a sequence of the left side and right side
181 Sequence* left = new Sequence(queryName, fragLeft);
182 Sequence* right = new Sequence(queryName, fragRight);
184 //find seqs closest to each fragment
185 Sequence closestLeft = templateDB->findClosestSequence(left);
187 Sequence closestRight = templateDB->findClosestSequence(right);
189 //get kmerinfo for the closest left
190 vector< map<int, int> > closeLeftKmerInfo = kmer->getKmerCounts(closestLeft.getUnaligned());
192 //get kmerinfo for the closest right
193 vector< map<int, int> > closeRightKmerInfo = kmer->getKmerCounts(closestRight.getUnaligned());
195 //right side is tricky - since the counts grow on eachother to find the correct counts of only the right side you must subtract the counts of the left side
196 //iterate through left sides map to subtract the number of times you saw things before you got the the right side
197 map<int, int> rightside = queryKmerInfo[queryKmerInfo.size()-1];
198 for (map<int, int>::iterator itleft = queryKmerInfo[m-kmerSize].begin(); itleft != queryKmerInfo[m-kmerSize].end(); itleft++) {
199 int howManyTotal = queryKmerInfo[queryKmerInfo.size()-1][itleft->first]; //times that kmer was seen in total
201 //itleft->second is times it was seen in left side, so howmanytotal - leftside should give you right side
202 int howmanyright = howManyTotal - itleft->second;
204 //if any were seen just on the left erase
205 if (howmanyright == 0) {
206 rightside.erase(itleft->first);
210 map<int, int> closerightside = closeRightKmerInfo[closeRightKmerInfo.size()-1];
211 for (map<int, int>::iterator itright = closeRightKmerInfo[m-kmerSize].begin(); itright != closeRightKmerInfo[m-kmerSize].end(); itright++) {
212 int howManyTotal = closeRightKmerInfo[(closeRightKmerInfo.size()-1)][itright->first]; //times that kmer was seen in total
214 //itleft->second is times it was seen in left side, so howmanytotal - leftside should give you right side
215 int howmanyright = howManyTotal - itright->second;
217 //if any were seen just on the left erase
218 if (howmanyright == 0) {
219 closerightside.erase(itright->first);
224 int nLeft = calcKmers(closeLeftKmerInfo[m-kmerSize], queryKmerInfo[m-kmerSize]);
226 int nRight = calcKmers(closerightside, rightside);
228 int is = nLeft + nRight - nTotal;
230 //save IS, leftparent, rightparent, breakpoint
231 temp.leftParent = closestLeft.getName();
232 temp.rightParent = closestRight.getName();
236 isValues.push_back(temp);
245 catch(exception& e) {
246 errorOut(e, "ChimeraCheckRDP", "findIS");
250 //***************************************************************************************************************
251 void ChimeraCheckRDP::readName(string namefile) {
254 openInputFile(namefile, in);
267 catch(exception& e) {
268 errorOut(e, "ChimeraCheckRDP", "readName");
273 //***************************************************************************************************************
274 //find the smaller map and iterate through it and count kmers in common
275 int ChimeraCheckRDP::calcKmers(map<int, int> query, map<int, int> subject) {
279 map<int, int>::iterator small;
280 map<int, int>::iterator large;
282 if (query.size() < subject.size()) {
284 for (small = query.begin(); small != query.end(); small++) {
285 large = subject.find(small->first);
287 //if you found it they have that kmer in common
288 if (large != subject.end()) { common++; }
293 for (small = subject.begin(); small != subject.end(); small++) {
294 large = query.find(small->first);
296 //if you found it they have that kmer in common
297 if (large != query.end()) { common++; }
304 catch(exception& e) {
305 errorOut(e, "ChimeraCheckRDP", "calcKmers");
310 //***************************************************************************************************************
311 void ChimeraCheckRDP::getCutoff() {
316 //store all is scores for all windows
317 for (int i = 0; i < IS.size(); i++) {
318 for (int j = 0; j < IS[i].size(); j++) {
319 temp.push_back(IS[i][j].score);
324 sort(temp.begin(), temp.end());
327 chimeraCutoff = temp[int(temp.size() * 0.95)];
330 catch(exception& e) {
331 errorOut(e, "ChimeraCheckRDP", "getCutoff");
336 //***************************************************************************************************************
337 void ChimeraCheckRDP::makeSVGpic(vector<sim> info, int query) {
340 string file = querySeqs[query]->getName() + ".chimeracheck.svg";
342 openOutputFile(file, outsvg);
344 int width = (info.size()*5) + 150;
346 outsvg << "<svg xmlns:svg=\"http://www.w3.org/2000/svg\" xmlns=\"http://www.w3.org/2000/svg\" width=\"100%\" height=\"100%\" viewBox=\"0 0 700 " + toString(width) + "\">\n";
348 outsvg << "<text fill=\"black\" class=\"seri\" x=\"" + toString((width / 2) - 150) + "\" y=\"25\">Plotted IS values for " + querySeqs[query]->getName() + "</text>\n";
350 outsvg << "<line x1=\"75\" y1=\"600\" x2=\"" + toString((info.size()*5) + 75) + "\" y2=\"600\" stroke=\"black\" stroke-width=\"2\"/>\n";
351 outsvg << "<line x1=\"75\" y1=\"600\" x2=\"75\" y2=\"125\" stroke=\"black\" stroke-width=\"2\"/>\n";
353 outsvg << "<text fill=\"black\" class=\"seri\" x=\"80\" y=\"620\">" + toString(info[0].midpoint) + "</text>\n";
354 outsvg << "<text fill=\"black\" class=\"seri\" x=\"" + toString((info.size()*5) + 75) + "\" y=\"620\">" + toString(info[info.size()-1].midpoint) + "</text>\n";
355 outsvg << "<text fill=\"black\" class=\"seri\" x=\"" + toString((width / 2) - 150) + "\" y=\"650\">Base Positions</text>\n";
357 outsvg << "<text fill=\"black\" class=\"seri\" x=\"50\" y=\"580\">0</text>\n";
359 outsvg << "<text fill=\"black\" class=\"seri\" x=\"50\" y=\"350\">IS</text>\n";
364 for (int i = 0; i < info.size(); i++) {
365 if (info[i].score > biggest) {
366 biggest = info[i].score;
370 outsvg << "<text fill=\"black\" class=\"seri\" x=\"50\" y=\"135\">" + toString(biggest) + "</text>\n";
372 int scaler2 = 500 / biggest;
375 outsvg << "<polyline fill=\"none\" stroke=\"red\" stroke-width=\"2\" points=\"";
376 //160,200 180,230 200,210 234,220\"/> ";
377 for (int i = 0; i < info.size(); i++) {
378 if(info[i].score < 0) { info[i].score = 0; }
379 outsvg << ((i*5) + 75) << "," << (600 - (info[i].score * scaler2)) << " ";
383 outsvg << "</g>\n</svg>\n";
388 catch(exception& e) {
389 errorOut(e, "ChimeraCheckRDP", "makeSVGpic");
393 //***************************************************************************************************************
394 void ChimeraCheckRDP::createProcessesIS() {
396 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
398 vector<int> processIDS;
400 //loop through and create all the processes you want
401 while (process != processors) {
405 processIDS.push_back(pid);
409 for (int i = lines[process]->start; i < lines[process]->end; i++) {
413 //write out data to file so parent can read it
415 string s = toString(getpid()) + ".temp";
416 openOutputFile(s, out);
419 for (int i = lines[process]->start; i < lines[process]->end; i++) {
420 out << IS[i].size() << endl;
421 for (int j = 0; j < IS[i].size(); j++) {
422 out << IS[i][j].leftParent << '\t'<< IS[i][j].rightParent << '\t' << IS[i][j].midpoint << '\t' << IS[i][j].score << endl;
429 }else { mothurOut("unable to spawn the necessary processes."); mothurOutEndLine(); exit(0); }
432 //force parent to wait until all the processes are done
433 for (int i=0;i<processors;i++) {
434 int temp = processIDS[i];
438 //get data created by processes
439 for (int i=0;i<processors;i++) {
441 string s = toString(processIDS[i]) + ".temp";
442 openInputFile(s, in);
445 for (int k = lines[i]->start; k < lines[i]->end; k++) {
448 in >> size; gobble(in);
456 for (int j = 0; j < size; j++) {
457 in >> left >> right >> mid >> score; gobble(in);
460 temp.leftParent = left;
461 temp.rightParent = right;
465 IS[k].push_back(temp);
473 for (int i = 0; i < querySeqs.size(); i++) {
478 catch(exception& e) {
479 errorOut(e, "ChimeraCheckRDP", "createProcessesIS");
484 //***************************************************************************************************************