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
80 mothurOut("Reading sequences and template file... "); cout.flush();
81 querySeqs = readSeqs(fastafile);
82 //templateSeqs = readSeqs(templateFile);
83 templateDB = new KmerDB(templateFile, kmerSize);
84 mothurOut("Done."); mothurOutEndLine();
86 int numSeqs = querySeqs.size();
89 closest.resize(numSeqs);
91 //break up file if needed
92 int linesPerProcess = numSeqs / processors ;
94 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
95 //find breakup of sequences for all times we will Parallelize
96 if (processors == 1) { lines.push_back(new linePair(0, numSeqs)); }
99 for (int i = 0; i < (processors-1); i++) {
100 lines.push_back(new linePair((i*linesPerProcess), ((i*linesPerProcess) + linesPerProcess)));
102 //this is necessary to get remainder of processors / numSeqs so you don't miss any lines at the end
103 int i = processors - 1;
104 lines.push_back(new linePair((i*linesPerProcess), numSeqs));
108 lines.push_back(new linePair(0, numSeqs));
111 kmer = new Kmer(kmerSize);
114 readName(name); //fills name map with names of seqs the user wants to have .svg for.
117 //find closest seq to each querySeq
118 for (int i = 0; i < querySeqs.size(); i++) {
119 closest[i] = templateDB->findClosestSequence(querySeqs[i]);
122 //for each query find IS value
123 if (processors == 1) {
124 for (int i = 0; i < querySeqs.size(); i++) {
127 }else { createProcessesIS(); }
129 //determine chimera report cutoff - window score above 95%
130 //getCutoff(); - not very acurate predictor
133 for (int i = 0; i < lines.size(); i++) { delete lines[i]; }
137 catch(exception& e) {
138 errorOut(e, "ChimeraCheckRDP", "getChimeras");
142 //***************************************************************************************************************
143 vector<sim> ChimeraCheckRDP::findIS(int query) {
147 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
148 //example: seqKmerInfo[50] = map containing the kmers found in the first 50 + kmersize characters of ecoli.
149 //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
150 //kmers 2 seqs had in common. There may be a better way to do this thats why I am leaving so many comments...
151 vector< map<int, int> > subjectKmerInfo;
153 vector<sim> isValues;
154 string queryName = querySeqs[query]->getName();
155 string seq = querySeqs[query]->getUnaligned();
157 mothurOut("Finding IS values for sequence " + toString(query+1)); mothurOutEndLine();
159 queryKmerInfo = kmer->getKmerCounts(seq);
160 subjectKmerInfo = kmer->getKmerCounts(closest[query].getUnaligned());
162 //find total kmers you have in common with closest[query] by looking at the last entry in the vector of maps for each
163 int nTotal = calcKmers(queryKmerInfo[(queryKmerInfo.size()-1)], subjectKmerInfo[(subjectKmerInfo.size()-1)]);
165 //you don't want the starting point to be virtually at hte end so move it in 10%
166 int start = seq.length() / 10;
169 for (int m = start; m < (seq.length() - start); m+=increment) {
171 if ((m - kmerSize) < 0) { mothurOut("Your sequence is too short for your kmerSize."); mothurOutEndLine(); exit(1); }
175 string fragLeft = seq.substr(0, m); //left side of breakpoint
176 string fragRight = seq.substr(m); //right side of breakpoint
178 //make a sequence of the left side and right side
179 Sequence* left = new Sequence(queryName, fragLeft);
180 Sequence* right = new Sequence(queryName, fragRight);
182 //find seqs closest to each fragment
183 Sequence closestLeft = templateDB->findClosestSequence(left);
185 Sequence closestRight = templateDB->findClosestSequence(right);
187 //get kmerinfo for the closest left
188 vector< map<int, int> > closeLeftKmerInfo = kmer->getKmerCounts(closestLeft.getUnaligned());
190 //get kmerinfo for the closest right
191 vector< map<int, int> > closeRightKmerInfo = kmer->getKmerCounts(closestRight.getUnaligned());
193 //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
194 //iterate through left sides map to subtract the number of times you saw things before you got the the right side
195 map<int, int> rightside = queryKmerInfo[queryKmerInfo.size()-1];
196 for (map<int, int>::iterator itleft = queryKmerInfo[m-kmerSize].begin(); itleft != queryKmerInfo[m-kmerSize].end(); itleft++) {
197 int howManyTotal = queryKmerInfo[queryKmerInfo.size()-1][itleft->first]; //times that kmer was seen in total
199 //itleft->second is times it was seen in left side, so howmanytotal - leftside should give you right side
200 int howmanyright = howManyTotal - itleft->second;
202 //if any were seen just on the left erase
203 if (howmanyright == 0) {
204 rightside.erase(itleft->first);
208 map<int, int> closerightside = closeRightKmerInfo[closeRightKmerInfo.size()-1];
209 for (map<int, int>::iterator itright = closeRightKmerInfo[m-kmerSize].begin(); itright != closeRightKmerInfo[m-kmerSize].end(); itright++) {
210 int howManyTotal = closeRightKmerInfo[(closeRightKmerInfo.size()-1)][itright->first]; //times that kmer was seen in total
212 //itleft->second is times it was seen in left side, so howmanytotal - leftside should give you right side
213 int howmanyright = howManyTotal - itright->second;
215 //if any were seen just on the left erase
216 if (howmanyright == 0) {
217 closerightside.erase(itright->first);
222 int nLeft = calcKmers(closeLeftKmerInfo[m-kmerSize], queryKmerInfo[m-kmerSize]);
224 int nRight = calcKmers(closerightside, rightside);
226 int is = nLeft + nRight - nTotal;
228 //save IS, leftparent, rightparent, breakpoint
229 temp.leftParent = closestLeft.getName();
230 temp.rightParent = closestRight.getName();
234 isValues.push_back(temp);
243 catch(exception& e) {
244 errorOut(e, "ChimeraCheckRDP", "findIS");
248 //***************************************************************************************************************
249 void ChimeraCheckRDP::readName(string namefile) {
252 openInputFile(namefile, in);
265 catch(exception& e) {
266 errorOut(e, "ChimeraCheckRDP", "readName");
271 //***************************************************************************************************************
272 //find the smaller map and iterate through it and count kmers in common
273 int ChimeraCheckRDP::calcKmers(map<int, int> query, map<int, int> subject) {
277 map<int, int>::iterator small;
278 map<int, int>::iterator large;
280 if (query.size() < subject.size()) {
282 for (small = query.begin(); small != query.end(); small++) {
283 large = subject.find(small->first);
285 //if you found it they have that kmer in common
286 if (large != subject.end()) { common++; }
291 for (small = subject.begin(); small != subject.end(); small++) {
292 large = query.find(small->first);
294 //if you found it they have that kmer in common
295 if (large != query.end()) { common++; }
302 catch(exception& e) {
303 errorOut(e, "ChimeraCheckRDP", "calcKmers");
308 //***************************************************************************************************************
309 void ChimeraCheckRDP::getCutoff() {
314 //store all is scores for all windows
315 for (int i = 0; i < IS.size(); i++) {
316 for (int j = 0; j < IS[i].size(); j++) {
317 temp.push_back(IS[i][j].score);
322 sort(temp.begin(), temp.end());
325 chimeraCutoff = temp[int(temp.size() * 0.95)];
328 catch(exception& e) {
329 errorOut(e, "ChimeraCheckRDP", "getCutoff");
334 //***************************************************************************************************************
335 void ChimeraCheckRDP::makeSVGpic(vector<sim> info, int query) {
338 string file = querySeqs[query]->getName() + ".chimeracheck.svg";
340 openOutputFile(file, outsvg);
342 int width = (info.size()*5) + 150;
344 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";
346 outsvg << "<text fill=\"black\" class=\"seri\" x=\"" + toString((width / 2) - 150) + "\" y=\"25\">Plotted IS values for " + querySeqs[query]->getName() + "</text>\n";
348 outsvg << "<line x1=\"75\" y1=\"600\" x2=\"" + toString((info.size()*5) + 75) + "\" y2=\"600\" stroke=\"black\" stroke-width=\"2\"/>\n";
349 outsvg << "<line x1=\"75\" y1=\"600\" x2=\"75\" y2=\"125\" stroke=\"black\" stroke-width=\"2\"/>\n";
351 outsvg << "<text fill=\"black\" class=\"seri\" x=\"80\" y=\"620\">" + toString(info[0].midpoint) + "</text>\n";
352 outsvg << "<text fill=\"black\" class=\"seri\" x=\"" + toString((info.size()*5) + 75) + "\" y=\"620\">" + toString(info[info.size()-1].midpoint) + "</text>\n";
353 outsvg << "<text fill=\"black\" class=\"seri\" x=\"" + toString((width / 2) - 150) + "\" y=\"650\">Base Positions</text>\n";
355 outsvg << "<text fill=\"black\" class=\"seri\" x=\"50\" y=\"580\">0</text>\n";
357 outsvg << "<text fill=\"black\" class=\"seri\" x=\"50\" y=\"350\">IS</text>\n";
362 for (int i = 0; i < info.size(); i++) {
363 if (info[i].score > biggest) {
364 biggest = info[i].score;
368 outsvg << "<text fill=\"black\" class=\"seri\" x=\"50\" y=\"135\">" + toString(biggest) + "</text>\n";
370 int scaler2 = 500 / biggest;
373 outsvg << "<polyline fill=\"none\" stroke=\"red\" stroke-width=\"2\" points=\"";
374 //160,200 180,230 200,210 234,220\"/> ";
375 for (int i = 0; i < info.size(); i++) {
376 if(info[i].score < 0) { info[i].score = 0; }
377 outsvg << ((i*5) + 75) << "," << (600 - (info[i].score * scaler2)) << " ";
381 outsvg << "</g>\n</svg>\n";
386 catch(exception& e) {
387 errorOut(e, "ChimeraCheckRDP", "makeSVGpic");
391 //***************************************************************************************************************
392 void ChimeraCheckRDP::createProcessesIS() {
394 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
396 vector<int> processIDS;
398 //loop through and create all the processes you want
399 while (process != processors) {
403 processIDS.push_back(pid);
407 for (int i = lines[process]->start; i < lines[process]->end; i++) {
411 //write out data to file so parent can read it
413 string s = toString(getpid()) + ".temp";
414 openOutputFile(s, out);
417 for (int i = lines[process]->start; i < lines[process]->end; i++) {
418 out << IS[i].size() << endl;
419 for (int j = 0; j < IS[i].size(); j++) {
420 out << IS[i][j].leftParent << '\t'<< IS[i][j].rightParent << '\t' << IS[i][j].midpoint << '\t' << IS[i][j].score << endl;
427 }else { mothurOut("unable to spawn the necessary processes."); mothurOutEndLine(); exit(0); }
430 //force parent to wait until all the processes are done
431 for (int i=0;i<processors;i++) {
432 int temp = processIDS[i];
436 //get data created by processes
437 for (int i=0;i<processors;i++) {
439 string s = toString(processIDS[i]) + ".temp";
440 openInputFile(s, in);
443 for (int k = lines[i]->start; k < lines[i]->end; k++) {
446 in >> size; gobble(in);
454 for (int j = 0; j < size; j++) {
455 in >> left >> right >> mid >> score; gobble(in);
458 temp.leftParent = left;
459 temp.rightParent = right;
463 IS[k].push_back(temp);
471 for (int i = 0; i < querySeqs.size(); i++) {
476 catch(exception& e) {
477 errorOut(e, "ChimeraCheckRDP", "createProcessesIS");
482 //***************************************************************************************************************