chimeracheckrdp.cpp

   1 /*
   2  *  chimeracheckrdp.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 9/8/09.
   6  *  Copyright 2009 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "chimeracheckrdp.h"
  11
  12 //***************************************************************************************************************
  13 ChimeraCheckRDP::ChimeraCheckRDP(string filename, string temp) {  fastafile = filename;  templateFile = temp;  }
  14 //***************************************************************************************************************
  15
  16 ChimeraCheckRDP::~ChimeraCheckRDP() {
  17         try {
  18                 for (int i = 0; i < querySeqs.size(); i++)              {  delete querySeqs[i];         }
  19                 delete templateDB;
  20                 delete kmer;
  21         }
  22         catch(exception& e) {
  23                 errorOut(e, "ChimeraCheckRDP", "~AlignSim");
  24                 exit(1);
  25         }
  26 }
  27 //***************************************************************************************************************
  28 void ChimeraCheckRDP::print(ostream& out) {
  29         try {
  30
  31                 mothurOutEndLine();
  32
  33                 //vector<bool> isChimeric;  isChimeric.resize(querySeqs.size(), false);
  34
  35                 for (int i = 0; i < querySeqs.size(); i++) {
  36
  37                                 out << querySeqs[i]->getName() << endl;
  38                                 out << "IS scores: " << '\t';
  39
  40                                 //int lastChimericWindowFound = 0;
  41
  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;           }
  45                                 }
  46
  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;
  50                                 //}else{
  51                                         //out << endl << "chimera: NO" << endl;
  52                                 //}
  53
  54                                 if (svg) {
  55
  56                                         if (name != "") { //if user has specific names
  57                                                 map<string, string>::iterator it = names.find(querySeqs[i]->getName());
  58
  59                                                 if (it != names.end()) { //user wants pic of this
  60                                                         makeSVGpic(IS[i], i);  //zeros out negative results
  61                                                 }
  62                                         }else{//output them all
  63                                                 makeSVGpic(IS[i], i);  //zeros out negative results
  64                                         }
  65                                 }
  66                 }
  67
  68         }
  69         catch(exception& e) {
  70                 errorOut(e, "ChimeraCheckRDP", "print");
  71                 exit(1);
  72         }
  73 }
  74
  75 //***************************************************************************************************************
  76 void ChimeraCheckRDP::getChimeras() {
  77         try {
  78
  79                 //read in query sequences and subject sequences
  80                 mothurOutEndLine();
  81                 mothurOut("Reading query sequences... "); cout.flush();
  82                 querySeqs = readSeqs(fastafile);
  83                 mothurOut("Done.");
  84                 //templateSeqs = readSeqs(templateFile);
  85                 templateDB = new AlignmentDB(templateFile, "kmer", kmerSize, 0.0,0.0,0.0,0.0);
  86                 mothurOutEndLine();
  87
  88                 int numSeqs = querySeqs.size();
  89
  90                 IS.resize(numSeqs);
  91                 closest.resize(numSeqs);
  92
  93                 //break up file if needed
  94                 int linesPerProcess = numSeqs / processors ;
  95
  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));  }
  99                         else {
 100                                 //fill line pairs
 101                                 for (int i = 0; i < (processors-1); i++) {
 102                                         lines.push_back(new linePair((i*linesPerProcess), ((i*linesPerProcess) + linesPerProcess)));
 103                                 }
 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));
 107                         }
 108
 109                 #else
 110                         lines.push_back(new linePair(0, numSeqs));
 111                 #endif
 112
 113                 kmer = new Kmer(kmerSize);
 114
 115                 if (name != "") {
 116                         readName(name);  //fills name map with names of seqs the user wants to have .svg for.
 117                 }
 118
 119                 //find closest seq to each querySeq
 120                 for (int i = 0; i < querySeqs.size(); i++) {
 121                         closest[i] = templateDB->findClosestSequence(querySeqs[i]);
 122                 }
 123
 124                 //for each query find IS value
 125                 if (processors == 1) {
 126                         for (int i = 0; i < querySeqs.size(); i++) {
 127                                 IS[i] = findIS(i);
 128                         }
 129                 }else {         createProcessesIS();    }
 130
 131                 //determine chimera report cutoff - window score above 95%
 132                 //getCutoff();  - not very acurate predictor
 133
 134                 //free memory
 135                 for (int i = 0; i < lines.size(); i++)                                  {       delete lines[i];        }
 136
 137
 138         }
 139         catch(exception& e) {
 140                 errorOut(e, "ChimeraCheckRDP", "getChimeras");
 141                 exit(1);
 142         }
 143 }
 144 //***************************************************************************************************************
 145 vector<sim> ChimeraCheckRDP::findIS(int query) {
 146         try {
 147
 148
 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;
 154
 155                 vector<sim>  isValues;
 156                 string queryName = querySeqs[query]->getName();
 157                 string seq = querySeqs[query]->getUnaligned();
 158
 159                 mothurOut("Finding IS values for sequence " + toString(query+1)); mothurOutEndLine();
 160
 161                 queryKmerInfo = kmer->getKmerCounts(seq);
 162                 subjectKmerInfo = kmer->getKmerCounts(closest[query].getUnaligned());
 163
 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)]);
 166
 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;
 169
 170                 //for each window
 171                 for (int m = start; m < (seq.length() - start); m+=increment) {
 172
 173                         if ((m - kmerSize) < 0)  { mothurOut("Your sequence is too short for your kmerSize."); mothurOutEndLine(); exit(1); }
 174
 175                         sim temp;
 176
 177                         string fragLeft = seq.substr(0, m);  //left side of breakpoint
 178                         string fragRight = seq.substr(m);  //right side of breakpoint
 179
 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);
 183
 184                         //find seqs closest to each fragment
 185                         Sequence closestLeft = templateDB->findClosestSequence(left);
 186
 187                         Sequence closestRight = templateDB->findClosestSequence(right);
 188
 189                         //get kmerinfo for the closest left
 190                         vector< map<int, int> > closeLeftKmerInfo = kmer->getKmerCounts(closestLeft.getUnaligned());
 191
 192                         //get kmerinfo for the closest right
 193                         vector< map<int, int> > closeRightKmerInfo = kmer->getKmerCounts(closestRight.getUnaligned());
 194
 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
 200
 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;
 203
 204                                 //if any were seen just on the left erase
 205                                 if (howmanyright == 0) {
 206                                         rightside.erase(itleft->first);
 207                                 }
 208                         }
 209
 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
 213
 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;
 216
 217                                 //if any were seen just on the left erase
 218                                 if (howmanyright == 0) {
 219                                         closerightside.erase(itright->first);
 220                                 }
 221                         }
 222
 223
 224                         int nLeft = calcKmers(closeLeftKmerInfo[m-kmerSize], queryKmerInfo[m-kmerSize]);
 225
 226                         int nRight = calcKmers(closerightside, rightside);
 227
 228                         int is = nLeft + nRight - nTotal;
 229
 230                         //save IS, leftparent, rightparent, breakpoint
 231                         temp.leftParent = closestLeft.getName();
 232                         temp.rightParent = closestRight.getName();
 233                         temp.score = is;
 234                         temp.midpoint = m;
 235
 236                         isValues.push_back(temp);
 237
 238                         delete left;
 239                         delete right;
 240                 }
 241
 242                 return isValues;
 243
 244         }
 245         catch(exception& e) {
 246                 errorOut(e, "ChimeraCheckRDP", "findIS");
 247                 exit(1);
 248         }
 249 }
 250 //***************************************************************************************************************
 251 void ChimeraCheckRDP::readName(string namefile) {
 252         try{
 253                 ifstream in;
 254                 openInputFile(namefile, in);
 255                 string name;
 256
 257                 while (!in.eof()) {
 258
 259                         in >> name;
 260
 261                         names[name] = name;
 262
 263                         gobble(in);
 264                 }
 265
 266         }
 267         catch(exception& e) {
 268                 errorOut(e, "ChimeraCheckRDP", "readName");
 269                 exit(1);
 270         }
 271 }
 272
 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) {
 276         try{
 277
 278                 int common = 0;
 279                 map<int, int>::iterator small;
 280                 map<int, int>::iterator large;
 281
 282                 if (query.size() < subject.size()) {
 283
 284                         for (small = query.begin(); small != query.end(); small++) {
 285                                 large = subject.find(small->first);
 286
 287                                 //if you found it they have that kmer in common
 288                                 if (large != subject.end()) {   common++;       }
 289                         }
 290
 291                 }else {
 292
 293                         for (small = subject.begin(); small != subject.end(); small++) {
 294                                 large = query.find(small->first);
 295
 296                                 //if you found it they have that kmer in common
 297                                 if (large != query.end()) {             common++;        }
 298                         }
 299                 }
 300
 301                 return common;
 302
 303         }
 304         catch(exception& e) {
 305                 errorOut(e, "ChimeraCheckRDP", "calcKmers");
 306                 exit(1);
 307         }
 308 }
 309
 310 //***************************************************************************************************************
 311 void ChimeraCheckRDP::getCutoff() {
 312         try{
 313
 314                 vector<float> temp;
 315
 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);
 320                         }
 321                 }
 322
 323                 //sort them
 324                 sort(temp.begin(), temp.end());
 325
 326                 //get 95%
 327                 chimeraCutoff = temp[int(temp.size() * 0.95)];
 328
 329         }
 330         catch(exception& e) {
 331                 errorOut(e, "ChimeraCheckRDP", "getCutoff");
 332                 exit(1);
 333         }
 334 }
 335
 336 //***************************************************************************************************************
 337 void ChimeraCheckRDP::makeSVGpic(vector<sim> info, int query) {
 338         try{
 339
 340                 string file = querySeqs[query]->getName() + ".chimeracheck.svg";
 341                 ofstream outsvg;
 342                 openOutputFile(file, outsvg);
 343
 344                 int width = (info.size()*5) + 150;
 345
 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";
 347                 outsvg << "<g>\n";
 348                 outsvg << "<text fill=\"black\" class=\"seri\" x=\"" + toString((width / 2) - 150) + "\" y=\"25\">Plotted IS values for " + querySeqs[query]->getName() + "</text>\n";
 349
 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";
 352
 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";
 356
 357                 outsvg << "<text fill=\"black\" class=\"seri\" x=\"50\" y=\"580\">0</text>\n";
 358
 359                 outsvg << "<text fill=\"black\" class=\"seri\" x=\"50\" y=\"350\">IS</text>\n";
 360
 361
 362                 //find max is score
 363                 float biggest = 0.0;
 364                 for (int i = 0; i < info.size(); i++) {
 365                         if (info[i].score > biggest)  {
 366                                 biggest = info[i].score;
 367                         }
 368                 }
 369
 370                 outsvg << "<text fill=\"black\" class=\"seri\" x=\"50\" y=\"135\">" + toString(biggest) + "</text>\n";
 371
 372                 int scaler2 = 500 / biggest;
 373
 374
 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)) << " ";
 380                 }
 381
 382                 outsvg << "\"/> ";
 383                 outsvg << "</g>\n</svg>\n";
 384
 385                 outsvg.close();
 386
 387         }
 388         catch(exception& e) {
 389                 errorOut(e, "ChimeraCheckRDP", "makeSVGpic");
 390                 exit(1);
 391         }
 392 }
 393 //***************************************************************************************************************
 394 void ChimeraCheckRDP::createProcessesIS() {
 395         try {
 396         #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 397                 int process = 0;
 398                 vector<int> processIDS;
 399
 400                 //loop through and create all the processes you want
 401                 while (process != processors) {
 402                         int pid = fork();
 403
 404                         if (pid > 0) {
 405                                 processIDS.push_back(pid);
 406                                 process++;
 407                         }else if (pid == 0){
 408
 409                                 for (int i = lines[process]->start; i < lines[process]->end; i++) {
 410                                         IS[i] = findIS(i);
 411                                 }
 412
 413                                 //write out data to file so parent can read it
 414                                 ofstream out;
 415                                 string s = toString(getpid()) + ".temp";
 416                                 openOutputFile(s, out);
 417
 418                                 //output pairs
 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;
 423                                         }
 424                                 }
 425
 426                                 out.close();
 427
 428                                 exit(0);
 429                         }else { mothurOut("unable to spawn the necessary processes."); mothurOutEndLine(); exit(0); }
 430                 }
 431
 432                 //force parent to wait until all the processes are done
 433                 for (int i=0;i<processors;i++) {
 434                         int temp = processIDS[i];
 435                         wait(&temp);
 436                 }
 437
 438                 //get data created by processes
 439                 for (int i=0;i<processors;i++) {
 440                         ifstream in;
 441                         string s = toString(processIDS[i]) + ".temp";
 442                         openInputFile(s, in);
 443
 444                         //get pairs
 445                         for (int k = lines[i]->start; k < lines[i]->end; k++) {
 446
 447                                 int size;
 448                                 in >> size; gobble(in);
 449
 450                                 string left, right;
 451                                 int mid;
 452                                 float score;
 453
 454                                 IS[k].clear();
 455
 456                                 for (int j = 0; j < size; j++) {
 457                                         in >> left >> right >> mid >> score;  gobble(in);
 458
 459                                         sim temp;
 460                                         temp.leftParent = left;
 461                                         temp.rightParent = right;
 462                                         temp.midpoint = mid;
 463                                         temp.score = score;
 464
 465                                         IS[k].push_back(temp);
 466                                 }
 467                         }
 468
 469                         in.close();
 470                         remove(s.c_str());
 471                 }
 472 #else
 473                         for (int i = 0; i < querySeqs.size(); i++) {
 474                                 IS[i] = findIS(i);
 475                         }
 476 #endif
 477         }
 478         catch(exception& e) {
 479                 errorOut(e, "ChimeraCheckRDP", "createProcessesIS");
 480                 exit(1);
 481         }
 482 }
 483
 484 //***************************************************************************************************************
 485
 486