--- /dev/null
+/*
+ * alignedsimilarity.cpp
+ * Mothur
+ *
+ * Created by westcott on 8/18/09.
+ * Copyright 2009 Schloss Lab. All rights reserved.
+ *
+ */
+
+#include "alignedsimilarity.h"
+#include "ignoregaps.h"
+
+
+//***************************************************************************************************************
+AlignSim::AlignSim(string filename, string temp) { fastafile = filename; templateFile = temp; }
+//***************************************************************************************************************
+
+AlignSim::~AlignSim() {
+ try {
+ for (int i = 0; i < querySeqs.size(); i++) { delete querySeqs[i]; }
+ for (int i = 0; i < templateSeqs.size(); i++) { delete templateSeqs[i]; }
+ delete distCalc;
+ }
+ catch(exception& e) {
+ errorOut(e, "AlignSim", "~AlignSim");
+ exit(1);
+ }
+}
+//***************************************************************************************************************
+void AlignSim::print(ostream& out) {
+ try {
+
+ mothurOutEndLine();
+
+ for (int i = 0; i < querySeqs.size(); i++) {
+
+ int j = 0; float largest = -10;
+ //find largest sim value
+ for (int k = 0; k < IS[i].size(); k++) {
+ //is this score larger
+ if (IS[i][k].score > largest) {
+ j = k;
+ largest = IS[i][k].score;
+ }
+ }
+
+ //find parental similarity
+ distCalc->calcDist(*(IS[i][j].leftParent), *(IS[i][j].rightParent));
+ float dist = distCalc->getDist();
+
+ //convert to similarity
+ dist = (1 - dist) * 100;
+
+ //warn about parental similarity - if its above 82% may not detect a chimera
+ if (dist >= 82) { mothurOut("When the chimeras parental similarity is above 82%, detection rates drop signifigantly."); mothurOutEndLine(); }
+
+ int index = ceil(dist);
+
+ if (index == 0) { index=1; }
+
+ //is your DE value higher than the 95%
+ string chimera;
+ if (IS[i][j].score > quantile[index-1][4]) { chimera = "Yes"; }
+ else { chimera = "No"; }
+
+ out << querySeqs[i]->getName() << "\tparental similarity: " << dist << "\tIS: " << IS[i][j].score << "\tbreakpoint: " << IS[i][j].midpoint << "\tchimera flag: " << chimera << endl;
+
+ if (chimera == "Yes") {
+ mothurOut(querySeqs[i]->getName() + "\tparental similarity: " + toString(dist) + "\tIS: " + toString(IS[i][j].score) + "\tbreakpoint: " + toString(IS[i][j].midpoint) + "\tchimera flag: " + chimera); mothurOutEndLine();
+ }
+ out << "Improvement Score\t";
+
+ for (int r = 0; r < IS[i].size(); r++) { out << IS[i][r].score << '\t'; }
+ out << endl;
+ }
+ }
+ catch(exception& e) {
+ errorOut(e, "AlignSim", "print");
+ exit(1);
+ }
+}
+
+//***************************************************************************************************************
+void AlignSim::getChimeras() {
+ try {
+
+ //read in query sequences and subject sequences
+ mothurOut("Reading sequences and template file... "); cout.flush();
+ querySeqs = readSeqs(fastafile);
+ templateSeqs = readSeqs(templateFile);
+ mothurOut("Done."); mothurOutEndLine();
+
+ int numSeqs = querySeqs.size();
+
+ IS.resize(numSeqs);
+
+ //break up file if needed
+ int linesPerProcess = numSeqs / processors ;
+
+ #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+ //find breakup of sequences for all times we will Parallelize
+ if (processors == 1) { lines.push_back(new linePair(0, numSeqs)); }
+ else {
+ //fill line pairs
+ for (int i = 0; i < (processors-1); i++) {
+ lines.push_back(new linePair((i*linesPerProcess), ((i*linesPerProcess) + linesPerProcess)));
+ }
+ //this is necessary to get remainder of processors / numSeqs so you don't miss any lines at the end
+ int i = processors - 1;
+ lines.push_back(new linePair((i*linesPerProcess), numSeqs));
+ }
+
+ //find breakup of templatefile for quantiles
+ if (processors == 1) { templateLines.push_back(new linePair(0, templateSeqs.size())); }
+ else {
+ for (int i = 0; i < processors; i++) {
+ templateLines.push_back(new linePair());
+ templateLines[i]->start = int (sqrt(float(i)/float(processors)) * templateSeqs.size());
+ templateLines[i]->end = int (sqrt(float(i+1)/float(processors)) * templateSeqs.size());
+ }
+ }
+ #else
+ lines.push_back(new linePair(0, numSeqs));
+ templateLines.push_back(new linePair(0, templateSeqs.size()));
+ #endif
+
+
+ distCalc = new ignoreGaps();
+
+ //find window breaks
+ windowBreak = findWindows();
+//for (int i = 0; i < windowBreak.size(); i++) { cout << windowBreak[i] << '\t'; }
+//cout << endl;
+
+ mothurOut("Finding the IS values for your sequences..."); cout.flush();
+ if (processors == 1) {
+ IS = findIS(lines[0]->start, lines[0]->end, querySeqs);
+ }else { IS = createProcessesIS(querySeqs, lines); }
+ mothurOut("Done."); mothurOutEndLine();
+
+ //quantiles are used to determine whether the de values found indicate a chimera
+ if (quanfile != "") { quantile = readQuantiles(); }
+ else {
+
+ mothurOut("Calculating quantiles for your template. This can take a while... I will output the quantiles to a .alignedsimilarity.quan file that you can input them using the quantiles parameter next time you run this command. Providing the .quan file will dramatically improve speed. "); cout.flush();
+ //get IS quantiles as a reference to these IS scores
+ //you are assuming the template is free of chimeras and therefore will give you a baseline as to the scores you would like to see
+ if (processors == 1) {
+ templateIS = findIS(templateLines[0]->start, templateLines[0]->end, templateSeqs);
+ }else { templateIS = createProcessesIS(templateSeqs, templateLines); }
+//cout << "here" << endl;
+ ofstream out4;
+ string o;
+
+ o = getRootName(templateFile) + "alignedsimilarity.quan";
+
+ openOutputFile(o, out4);
+
+ //adjust quantiles
+ //construct table
+ quantile.resize(100);
+
+ for (int i = 0; i < templateIS.size(); i++) {
+
+ if (templateIS[i].size() != 0) {
+ int j = 0; float score = -1000;
+ //find highest IS score
+ //cout << templateIS[i].size() << endl;
+ for (int k = 0; k < templateIS[i].size(); k++) {
+ if (templateIS[i][k].score > score) {
+ score = templateIS[i][k].score;
+ j = k;
+ }
+ }
+ //cout << j << endl;
+ //find similarity of parents
+ distCalc->calcDist(*(templateIS[i][j].leftParent), *(templateIS[i][j].rightParent));
+ float dist = distCalc->getDist();
+
+ //convert to similarity
+ dist = (1 - dist) * 100;
+ // cout << dist << endl;
+ int index = ceil(dist);
+ // cout << "index = " << index << endl;
+ if (index == 0) { index = 1; }
+ quantile[index-1].push_back(templateIS[i][j].score);
+ }
+ }
+
+
+ for (int i = 0; i < quantile.size(); i++) {
+ vector<float> temp;
+
+ if (quantile[i].size() == 0) {
+ //in case this is not a distance found in your template files
+ for (int g = 0; g < 6; g++) {
+ temp.push_back(0.0);
+ }
+ }else{
+
+ sort(quantile[i].begin(), quantile[i].end());
+
+ //save 10%
+ temp.push_back(quantile[i][int(quantile[i].size() * 0.10)]);
+ //save 25%
+ temp.push_back(quantile[i][int(quantile[i].size() * 0.25)]);
+ //save 50%
+ temp.push_back(quantile[i][int(quantile[i].size() * 0.5)]);
+ //save 75%
+ temp.push_back(quantile[i][int(quantile[i].size() * 0.75)]);
+ //save 95%
+ temp.push_back(quantile[i][int(quantile[i].size() * 0.95)]);
+ //save 99%
+ temp.push_back(quantile[i][int(quantile[i].size() * 0.99)]);
+
+ }
+
+ //output quan value
+ out4 << i+1 << '\t';
+ for (int u = 0; u < temp.size(); u++) { out4 << temp[u] << '\t'; }
+ out4 << endl;
+
+ quantile[i] = temp;
+
+ }
+
+ out4.close();
+
+ mothurOut("Done."); mothurOutEndLine();
+
+ }
+
+ //free memory
+ for (int i = 0; i < lines.size(); i++) { delete lines[i]; }
+ for (int i = 0; i < templateLines.size(); i++) { delete templateLines[i]; }
+
+ }
+ catch(exception& e) {
+ errorOut(e, "AlignSim", "getChimeras");
+ exit(1);
+ }
+}
+//***************************************************************************************************************
+vector<int> AlignSim::findWindows() {
+ try {
+
+ vector<int> win;
+
+ if (increment > querySeqs[0]->getAligned().length()) { mothurOut("You have selected an increment larger than the length of your sequences. I will use the default of 25."); increment = 25; }
+
+ for (int m = increment; m < (querySeqs[0]->getAligned().length() - increment); m+=increment) { win.push_back(m); }
+
+ return win;
+
+ }
+ catch(exception& e) {
+ errorOut(e, "AlignSim", "findWindows");
+ exit(1);
+ }
+}
+
+//***************************************************************************************************************
+vector< vector<sim> > AlignSim::findIS(int start, int end, vector<Sequence*> seqs) {
+ try {
+
+ vector< vector<sim> > isValues;
+ isValues.resize(seqs.size());
+
+ //for each sequence
+ for(int i = start; i < end; i++){
+
+ vector<sim> temp; temp.resize(windowBreak.size());
+
+ mothurOut("Finding IS value for sequence " + toString(i)); mothurOutEndLine();
+
+ //for each window
+ for (int m = 0; m < windowBreak.size(); m++) {
+
+ vector<Sequence*> closest; //left, right, overall
+ vector<int> similarity; //left+right and overall
+
+ //find closest left, closest right and closest overall
+ closest = findClosestSides(seqs[i], windowBreak[m], similarity, i);
+
+ int totalNumBases = seqs[i]->getUnaligned().length();
+
+ //IS = left+right-overall
+ float is = ((similarity[0]+similarity[1]) / (float) totalNumBases) - (similarity[2] / (float) totalNumBases);
+
+ ///save IS, leftparent, rightparent, breakpoint
+ temp[m].leftParent = closest[0];
+ temp[m].rightParent = closest[1];
+ temp[m].score = is;
+ temp[m].midpoint = windowBreak[m];
+//cout << is << '\t';
+ }
+//cout << endl;
+ isValues[i] = temp;
+
+ }
+
+ return isValues;
+
+ }
+ catch(exception& e) {
+ errorOut(e, "AlignSim", "findIS");
+ exit(1);
+ }
+}
+//***************************************************************************************************************
+vector<Sequence*> AlignSim::findClosestSides(Sequence* seq, int breakpoint, vector<int>& sim, int i) {
+ try{
+
+ vector<Sequence*> closest;
+
+ Sequence query, queryLeft, queryRight;
+ string frag = seq->getAligned();
+ string fragLeft = frag.substr(0, breakpoint);
+ string fragRight = frag.substr(breakpoint, frag.length());
+
+ //get pieces
+ query = *(seq);
+ queryLeft = *(seq);
+ queryRight = *(seq);
+ queryLeft.setAligned(fragLeft);
+ queryRight.setAligned(fragRight);
+
+ //initialize
+ Sequence* overall = templateSeqs[0];
+ Sequence* left = templateSeqs[0];
+ Sequence* right = templateSeqs[0];
+
+ float smallestOverall, smallestLeft, smallestRight;
+ smallestOverall = 1000; smallestLeft = 1000; smallestRight = 1000;
+
+ //go through the templateSeqs and search for the closest
+ for(int j = 0; j < templateSeqs.size(); j++){
+
+ //so you don't pick yourself
+ if (j != i) {
+ Sequence temp, tempLeft, tempRight;
+ string fragTemp = templateSeqs[j]->getAligned();
+ string fragTempLeft = fragTemp.substr(0, breakpoint);
+ string fragTempRight = fragTemp.substr(breakpoint, fragTemp.length());
+
+ //get pieces
+ temp = *(templateSeqs[j]);
+ tempLeft = *(templateSeqs[j]);
+ tempRight = *(templateSeqs[j]);
+ tempLeft.setAligned(fragTempLeft);
+ tempRight.setAligned(fragTempRight);
+
+ //find overall dist
+ distCalc->calcDist(query, temp);
+ float dist = distCalc->getDist();
+
+ if (dist < smallestOverall) {
+ overall = templateSeqs[j];
+ smallestOverall = dist;
+ }
+
+ //find left dist
+ distCalc->calcDist(queryLeft, tempLeft);
+ dist = distCalc->getDist();
+
+ if (dist < smallestLeft) {
+ left = templateSeqs[j];
+ smallestLeft = dist;
+ }
+
+ //find left dist
+ distCalc->calcDist(queryRight, tempRight);
+ dist = distCalc->getDist();
+
+ if (dist < smallestRight) {
+ right = templateSeqs[j];
+ smallestRight = dist;
+ }
+ }
+
+ }
+
+ closest.push_back(left);
+ closest.push_back(right);
+ closest.push_back(overall);
+
+ //fill sim with number of matched bases
+ Sequence tempL = *(left);
+ string tempLFrag = tempL.getAligned();
+ tempL.setAligned(tempLFrag.substr(0, breakpoint));
+
+ int bothMatches = findNumMatchedBases(queryLeft, tempL);
+ sim.push_back(bothMatches);
+
+ Sequence tempR = *(right);
+ string tempRFrag = tempR.getAligned();
+ tempR.setAligned(tempRFrag.substr(breakpoint, tempRFrag.length()));
+
+ bothMatches = findNumMatchedBases(queryRight, tempR);
+ sim.push_back(bothMatches);
+
+ bothMatches = findNumMatchedBases(query, *(overall));
+ sim.push_back(bothMatches);
+
+ return closest;
+
+
+ }
+ catch(exception& e) {
+ errorOut(e, "AlignSim", "findClosestSides");
+ exit(1);
+ }
+}
+//***************************************************************************************************************
+
+int AlignSim::findNumMatchedBases(Sequence seq, Sequence temp) {
+ try{
+
+ int num = 0;
+
+ string query = seq.getAligned();
+ string subject = temp.getAligned();
+//cout << seq.getName() << endl << endl;
+//cout << temp.getName() << endl << endl;
+
+ for (int i = 0; i < query.length(); i++) {
+ //count matches if query[i] is a base and subject is same bases or gap
+ if (isalpha(query[i])) {
+ if(!isalpha(subject[i])) { num++; }
+ else if (query[i] == subject[i]) { num++; }
+ else { }
+ }
+ }
+//cout << "num = " << num << endl;
+ return num;
+ }
+ catch(exception& e) {
+ errorOut(e, "AlignSim", "findNumMatchedBases");
+ exit(1);
+ }
+}
+
+/**************************************************************************************************/
+vector< vector<sim> > AlignSim::createProcessesIS(vector<Sequence*> seqs, vector<linePair*> linesToProcess) {
+ try {
+ vector< vector<sim> > localIs; localIs.resize(seqs.size());
+
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+ int process = 0;
+ vector<int> processIDS;
+
+ //loop through and create all the processes you want
+ while (process != processors) {
+ int pid = fork();
+
+ if (pid > 0) {
+ processIDS.push_back(pid);
+ process++;
+ }else if (pid == 0){
+
+ mothurOut("Finding IS values for sequences " + toString(linesToProcess[process]->start) + " to " + toString(linesToProcess[process]->end)); mothurOutEndLine();
+ localIs = findIS(linesToProcess[process]->start, linesToProcess[process]->end, seqs);
+ mothurOut("Done finding IS values for sequences " + toString(linesToProcess[process]->start) + " to " + toString(linesToProcess[process]->end)); mothurOutEndLine();
+
+ //write out data to file so parent can read it
+ ofstream out;
+ string s = toString(getpid()) + ".temp";
+ openOutputFile(s, out);
+
+ //output pairs
+ for (int i = linesToProcess[process]->start; i < linesToProcess[process]->end; i++) {
+ out << localIs[i].size() << endl;
+ for (int j = 0; j < localIs[i].size(); j++) {
+ localIs[i][j].leftParent->printSequence(out);
+ localIs[i][j].rightParent->printSequence(out);
+ out << ">" << '\t' << localIs[i][j].score << '\t' << localIs[i][j].midpoint << endl;
+ }
+ }
+ out.close();
+
+ exit(0);
+ }else { mothurOut("unable to spawn the necessary processes."); mothurOutEndLine(); exit(0); }
+ }
+
+ //force parent to wait until all the processes are done
+ for (int i=0;i<processors;i++) {
+ int temp = processIDS[i];
+ wait(&temp);
+ }
+
+ //get data created by processes
+ for (int i=0;i<processors;i++) {
+ ifstream in;
+ string s = toString(processIDS[i]) + ".temp";
+ openInputFile(s, in);
+
+ //get pairs
+ for (int k = linesToProcess[i]->start; k < linesToProcess[i]->end; k++) {
+ int size;
+ in >> size;
+ gobble(in);
+
+ vector<sim> tempVector;
+
+ for (int j = 0; j < size; j++) {
+
+ sim temp;
+
+ temp.leftParent = new Sequence(in);
+ gobble(in);
+ //temp.leftParent->printSequence(cout);
+ temp.rightParent = new Sequence(in);
+ gobble(in);
+//temp.rightParent->printSequence(cout);
+ string throwaway;
+ in >> throwaway >> temp.score >> temp.midpoint;
+ //cout << temp.score << '\t' << temp.midpoint << endl;
+ gobble(in);
+
+ tempVector.push_back(temp);
+ }
+
+ localIs[k] = tempVector;
+ }
+
+ in.close();
+ remove(s.c_str());
+ }
+
+
+#else
+ localIs = findIS(linesToProcess[0]->start, linesToProcess[0]->end, seqs);
+#endif
+
+ return localIs;
+ }
+ catch(exception& e) {
+ errorOut(e, "AlignSim", "createProcessesIS");
+ exit(1);
+ }
+}
+
+//***************************************************************************************************************
projects / mothur.git / blobdiff