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                 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();
  85
  86                 int numSeqs = querySeqs.size();
  87
  88                 IS.resize(numSeqs);
  89                 closest.resize(numSeqs);
  90
  91                 //break up file if needed
  92                 int linesPerProcess = numSeqs / processors ;
  93
  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));  }
  97                         else {
  98                                 //fill line pairs
  99                                 for (int i = 0; i < (processors-1); i++) {
 100                                         lines.push_back(new linePair((i*linesPerProcess), ((i*linesPerProcess) + linesPerProcess)));
 101                                 }
 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));
 105                         }
 106
 107                 #else
 108                         lines.push_back(new linePair(0, numSeqs));
 109                 #endif
 110
 111                 kmer = new Kmer(kmerSize);
 112
 113                 if (name != "") {
 114                         readName(name);  //fills name map with names of seqs the user wants to have .svg for.
 115                 }
 116
 117                 //find closest seq to each querySeq
 118                 for (int i = 0; i < querySeqs.size(); i++) {
 119                         closest[i] = templateDB->findClosestSequence(querySeqs[i]);
 120                 }
 121
 122                 //for each query find IS value
 123                 if (processors == 1) {
 124                         for (int i = 0; i < querySeqs.size(); i++) {
 125                                 IS[i] = findIS(i);
 126                         }
 127                 }else {         createProcessesIS();    }
 128
 129                 //determine chimera report cutoff - window score above 95%
 130                 //getCutoff();  - not very acurate predictor
 131
 132                 //free memory
 133                 for (int i = 0; i < lines.size(); i++)                                  {       delete lines[i];        }
 134
 135
 136         }
 137         catch(exception& e) {
 138                 errorOut(e, "ChimeraCheckRDP", "getChimeras");
 139                 exit(1);
 140         }
 141 }
 142 //***************************************************************************************************************
 143 vector<sim> ChimeraCheckRDP::findIS(int query) {
 144         try {
 145
 146
 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;
 152
 153                 vector<sim>  isValues;
 154                 string queryName = querySeqs[query]->getName();
 155                 string seq = querySeqs[query]->getUnaligned();
 156
 157                 mothurOut("Finding IS values for sequence " + toString(query+1)); mothurOutEndLine();
 158
 159                 queryKmerInfo = kmer->getKmerCounts(seq);
 160                 subjectKmerInfo = kmer->getKmerCounts(closest[query].getUnaligned());
 161
 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)]);
 164
 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;
 167
 168                 //for each window
 169                 for (int m = start; m < (seq.length() - start); m+=increment) {
 170
 171                         if ((m - kmerSize) < 0)  { mothurOut("Your sequence is too short for your kmerSize."); mothurOutEndLine(); exit(1); }
 172
 173                         sim temp;
 174
 175                         string fragLeft = seq.substr(0, m);  //left side of breakpoint
 176                         string fragRight = seq.substr(m);  //right side of breakpoint
 177
 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);
 181
 182                         //find seqs closest to each fragment
 183                         Sequence closestLeft = templateDB->findClosestSequence(left);
 184
 185                         Sequence closestRight = templateDB->findClosestSequence(right);
 186
 187                         //get kmerinfo for the closest left
 188                         vector< map<int, int> > closeLeftKmerInfo = kmer->getKmerCounts(closestLeft.getUnaligned());
 189
 190                         //get kmerinfo for the closest right
 191                         vector< map<int, int> > closeRightKmerInfo = kmer->getKmerCounts(closestRight.getUnaligned());
 192
 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
 198
 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;
 201
 202                                 //if any were seen just on the left erase
 203                                 if (howmanyright == 0) {
 204                                         rightside.erase(itleft->first);
 205                                 }
 206                         }
 207
 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
 211
 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;
 214
 215                                 //if any were seen just on the left erase
 216                                 if (howmanyright == 0) {
 217                                         closerightside.erase(itright->first);
 218                                 }
 219                         }
 220
 221
 222                         int nLeft = calcKmers(closeLeftKmerInfo[m-kmerSize], queryKmerInfo[m-kmerSize]);
 223
 224                         int nRight = calcKmers(closerightside, rightside);
 225
 226                         int is = nLeft + nRight - nTotal;
 227
 228                         //save IS, leftparent, rightparent, breakpoint
 229                         temp.leftParent = closestLeft.getName();
 230                         temp.rightParent = closestRight.getName();
 231                         temp.score = is;
 232                         temp.midpoint = m;
 233
 234                         isValues.push_back(temp);
 235
 236                         delete left;
 237                         delete right;
 238                 }
 239
 240                 return isValues;
 241
 242         }
 243         catch(exception& e) {
 244                 errorOut(e, "ChimeraCheckRDP", "findIS");
 245                 exit(1);
 246         }
 247 }
 248 //***************************************************************************************************************
 249 void ChimeraCheckRDP::readName(string namefile) {
 250         try{
 251                 ifstream in;
 252                 openInputFile(namefile, in);
 253                 string name;
 254
 255                 while (!in.eof()) {
 256
 257                         in >> name;
 258
 259                         names[name] = name;
 260
 261                         gobble(in);
 262                 }
 263
 264         }
 265         catch(exception& e) {
 266                 errorOut(e, "ChimeraCheckRDP", "readName");
 267                 exit(1);
 268         }
 269 }
 270
 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) {
 274         try{
 275
 276                 int common = 0;
 277                 map<int, int>::iterator small;
 278                 map<int, int>::iterator large;
 279
 280                 if (query.size() < subject.size()) {
 281
 282                         for (small = query.begin(); small != query.end(); small++) {
 283                                 large = subject.find(small->first);
 284
 285                                 //if you found it they have that kmer in common
 286                                 if (large != subject.end()) {   common++;       }
 287                         }
 288
 289                 }else {
 290
 291                         for (small = subject.begin(); small != subject.end(); small++) {
 292                                 large = query.find(small->first);
 293
 294                                 //if you found it they have that kmer in common
 295                                 if (large != query.end()) {             common++;        }
 296                         }
 297                 }
 298
 299                 return common;
 300
 301         }
 302         catch(exception& e) {
 303                 errorOut(e, "ChimeraCheckRDP", "calcKmers");
 304                 exit(1);
 305         }
 306 }
 307
 308 //***************************************************************************************************************
 309 void ChimeraCheckRDP::getCutoff() {
 310         try{
 311
 312                 vector<float> temp;
 313
 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);
 318                         }
 319                 }
 320
 321                 //sort them
 322                 sort(temp.begin(), temp.end());
 323
 324                 //get 95%
 325                 chimeraCutoff = temp[int(temp.size() * 0.95)];
 326
 327         }
 328         catch(exception& e) {
 329                 errorOut(e, "ChimeraCheckRDP", "getCutoff");
 330                 exit(1);
 331         }
 332 }
 333
 334 //***************************************************************************************************************
 335 void ChimeraCheckRDP::makeSVGpic(vector<sim> info, int query) {
 336         try{
 337
 338                 string file = querySeqs[query]->getName() + ".chimeracheck.svg";
 339                 ofstream outsvg;
 340                 openOutputFile(file, outsvg);
 341
 342                 int width = (info.size()*5) + 150;
 343
 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";
 345                 outsvg << "<g>\n";
 346                 outsvg << "<text fill=\"black\" class=\"seri\" x=\"" + toString((width / 2) - 150) + "\" y=\"25\">Plotted IS values for " + querySeqs[query]->getName() + "</text>\n";
 347
 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";
 350
 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";
 354
 355                 outsvg << "<text fill=\"black\" class=\"seri\" x=\"50\" y=\"580\">0</text>\n";
 356
 357                 outsvg << "<text fill=\"black\" class=\"seri\" x=\"50\" y=\"350\">IS</text>\n";
 358
 359
 360                 //find max is score
 361                 float biggest = 0.0;
 362                 for (int i = 0; i < info.size(); i++) {
 363                         if (info[i].score > biggest)  {
 364                                 biggest = info[i].score;
 365                         }
 366                 }
 367
 368                 outsvg << "<text fill=\"black\" class=\"seri\" x=\"50\" y=\"135\">" + toString(biggest) + "</text>\n";
 369
 370                 int scaler2 = 500 / biggest;
 371
 372
 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)) << " ";
 378                 }
 379
 380                 outsvg << "\"/> ";
 381                 outsvg << "</g>\n</svg>\n";
 382
 383                 outsvg.close();
 384
 385         }
 386         catch(exception& e) {
 387                 errorOut(e, "ChimeraCheckRDP", "makeSVGpic");
 388                 exit(1);
 389         }
 390 }
 391 //***************************************************************************************************************
 392 void ChimeraCheckRDP::createProcessesIS() {
 393         try {
 394         #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 395                 int process = 0;
 396                 vector<int> processIDS;
 397
 398                 //loop through and create all the processes you want
 399                 while (process != processors) {
 400                         int pid = fork();
 401
 402                         if (pid > 0) {
 403                                 processIDS.push_back(pid);
 404                                 process++;
 405                         }else if (pid == 0){
 406
 407                                 for (int i = lines[process]->start; i < lines[process]->end; i++) {
 408                                         IS[i] = findIS(i);
 409                                 }
 410
 411                                 //write out data to file so parent can read it
 412                                 ofstream out;
 413                                 string s = toString(getpid()) + ".temp";
 414                                 openOutputFile(s, out);
 415
 416                                 //output pairs
 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;
 421                                         }
 422                                 }
 423
 424                                 out.close();
 425
 426                                 exit(0);
 427                         }else { mothurOut("unable to spawn the necessary processes."); mothurOutEndLine(); exit(0); }
 428                 }
 429
 430                 //force parent to wait until all the processes are done
 431                 for (int i=0;i<processors;i++) {
 432                         int temp = processIDS[i];
 433                         wait(&temp);
 434                 }
 435
 436                 //get data created by processes
 437                 for (int i=0;i<processors;i++) {
 438                         ifstream in;
 439                         string s = toString(processIDS[i]) + ".temp";
 440                         openInputFile(s, in);
 441
 442                         //get pairs
 443                         for (int k = lines[i]->start; k < lines[i]->end; k++) {
 444
 445                                 int size;
 446                                 in >> size; gobble(in);
 447
 448                                 string left, right;
 449                                 int mid;
 450                                 float score;
 451
 452                                 IS[k].clear();
 453
 454                                 for (int j = 0; j < size; j++) {
 455                                         in >> left >> right >> mid >> score;  gobble(in);
 456
 457                                         sim temp;
 458                                         temp.leftParent = left;
 459                                         temp.rightParent = right;
 460                                         temp.midpoint = mid;
 461                                         temp.score = score;
 462
 463                                         IS[k].push_back(temp);
 464                                 }
 465                         }
 466
 467                         in.close();
 468                         remove(s.c_str());
 469                 }
 470 #else
 471                         for (int i = 0; i < querySeqs.size(); i++) {
 472                                 IS[i] = findIS(i);
 473                         }
 474 #endif
 475         }
 476         catch(exception& e) {
 477                 errorOut(e, "ChimeraCheckRDP", "createProcessesIS");
 478                 exit(1);
 479         }
 480 }
 481
 482 //***************************************************************************************************************
 483
 484