From 3abb236c602eb168ee112f080b563ebe2c705029 Mon Sep 17 00:00:00 2001 From: westcott Date: Mon, 20 Jul 2009 19:33:56 +0000 Subject: [PATCH] fixed bugs in venn and aligner --- aligncommand.cpp | 17 +- alignment.cpp | 123 ++++---- alignment.hpp | 3 +- bellerophon.h | 2 + blastalign.hpp | 2 + chao1.cpp | 2 +- chimera.h | 8 +- chimeraseqscommand.cpp | 22 +- chimeraseqscommand.h | 2 +- database.cpp | 16 +- database.hpp | 4 +- engine.cpp | 18 +- gotohoverlap.cpp | 99 +++--- gotohoverlap.hpp | 1 + kmerdb.cpp | 187 ++++++----- nast.cpp | 579 ++++++++++++++++++---------------- needlemanoverlap.cpp | 100 +++--- pintail.cpp | 699 +++++++++++++++++++++++++++-------------- pintail.h | 56 ++-- sharedchao1.cpp | 8 +- venn.cpp | 405 +++++++++++++++--------- 21 files changed, 1431 insertions(+), 922 deletions(-) diff --git a/aligncommand.cpp b/aligncommand.cpp index 2fda252..d672fc3 100644 --- a/aligncommand.cpp +++ b/aligncommand.cpp @@ -160,14 +160,16 @@ int AlignCommand::execute(){ templateDB = new KmerDB(templateFileName, kmerSize); } - if(align == "gotoh") { alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, 3000); } - else if(align == "needleman") { alignment = new NeedlemanOverlap(gapOpen, match, misMatch, 3000); } + int longestBase = templateDB->getLongestBase(); + + if(align == "gotoh") { alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase); } + else if(align == "needleman") { alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase); } else if(align == "blast") { alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch); } else if(align == "noalign") { alignment = new NoAlign(); } else { mothurOut(align + " is not a valid alignment option. I will run the command using needleman."); mothurOutEndLine(); - alignment = new NeedlemanOverlap(gapOpen, match, misMatch, 3000); + alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase); } string alignFileName = candidateFileName.substr(0,candidateFileName.find_last_of(".")+1) + "align"; @@ -184,7 +186,7 @@ int AlignCommand::execute(){ inFASTA.close(); lines.push_back(new linePair(0, numFastaSeqs)); - + driver(lines[0], alignFileName, reportFileName); } @@ -262,13 +264,14 @@ int AlignCommand::driver(linePair* line, string alignFName, string reportFName){ ifstream inFASTA; openInputFile(candidateFileName, inFASTA); + inFASTA.seekg(line->start); - + for(int i=0;inumSeqs;i++){ Sequence* candidateSeq = new Sequence(inFASTA); report.setCandidate(candidateSeq); - + Sequence temp = templateDB->findClosestSequence(candidateSeq); Sequence* templateSeq = &temp; @@ -276,7 +279,9 @@ int AlignCommand::driver(linePair* line, string alignFName, string reportFName){ report.setSearchParameters(search, templateDB->getSearchScore()); Nast nast(alignment, candidateSeq, templateSeq); + report.setAlignmentParameters(align, alignment); + report.setNastParameters(nast); alignmentFile << '>' << candidateSeq->getName() << '\n' << candidateSeq->getAligned() << endl; diff --git a/alignment.cpp b/alignment.cpp index 6014506..25b2760 100644 --- a/alignment.cpp +++ b/alignment.cpp @@ -20,72 +20,81 @@ Alignment::Alignment() { /* do nothing */ } /**************************************************************************************************/ Alignment::Alignment(int A) : nCols(A), nRows(A) { - - alignment.resize(nRows); // For the Gotoh and Needleman-Wunsch we initialize the dynamic programming - for(int i=0;i=0;i--){ + if(seqAaln[i] != '-' && seqBaln[i] == '-') { seqAend++; } + else if(seqAaln[i] == '-' && seqBaln[i] != '-') { seqBend++; } + else { break; } + } + pairwiseLength -= (seqAend + seqBend); + + seqAend = seqA.length() - seqAend - 1; + seqBend = seqB.length() - seqBend - 1; } - - pairwiseLength = seqAaln.length(); - seqAstart = 1; seqAend = 0; - seqBstart = 1; seqBend = 0; - - for(int i=0;i=0;i--){ - if(seqAaln[i] != '-' && seqBaln[i] == '-') { seqAend++; } - else if(seqAaln[i] == '-' && seqBaln[i] != '-') { seqBend++; } - else { break; } + catch(exception& e) { + errorOut(e, "Alignment", "traceBack"); + exit(1); } - pairwiseLength -= (seqAend + seqBend); - - seqAend = seqA.length() - seqAend - 1; - seqBend = seqB.length() - seqBend - 1; - } /**************************************************************************************************/ diff --git a/alignment.hpp b/alignment.hpp index 215e4b3..77b650c 100644 --- a/alignment.hpp +++ b/alignment.hpp @@ -26,6 +26,7 @@ public: virtual ~Alignment(); virtual void align(string, string) = 0; + // float getAlignmentScore(); string getSeqAAln(); string getSeqBAln(); @@ -33,7 +34,7 @@ public: int getCandidateEndPos(); int getTemplateStartPos(); int getTemplateEndPos(); - + int getPairwiseLength(); // int getLongestTemplateGap(); diff --git a/bellerophon.h b/bellerophon.h index 3f28e5b..05e93d7 100644 --- a/bellerophon.h +++ b/bellerophon.h @@ -39,6 +39,8 @@ class Bellerophon : public Chimera { void getChimeras(); void print(ostream&); + void setCons(string){}; + private: Dist* distCalculator; diff --git a/blastalign.hpp b/blastalign.hpp index 696c134..94bce82 100644 --- a/blastalign.hpp +++ b/blastalign.hpp @@ -20,6 +20,8 @@ public: BlastAlignment(float, float, float, float); ~BlastAlignment(); void align(string, string); + void setMatrix(int){}; + private: string candidateFileName; diff --git a/chao1.cpp b/chao1.cpp index b601b51..24521f5 100644 --- a/chao1.cpp +++ b/chao1.cpp @@ -19,7 +19,7 @@ EstOutput Chao1::getValues(SAbundVector* rank){ double singles = (double)rank->get(1); double doubles = (double)rank->get(2); double chaovar = 0.0000; - +//cout << "singles = " << singles << " doubles = " << doubles << " sobs = " << sobs << endl; double chao = sobs + singles*(singles-1)/(2*(doubles+1)); if(singles > 0 && doubles > 0){ diff --git a/chimera.h b/chimera.h index 07ad92f..e0e8944 100644 --- a/chimera.h +++ b/chimera.h @@ -26,13 +26,16 @@ class Chimera { Chimera(){}; Chimera(string); + Chimera(string, string); virtual ~Chimera(){}; virtual void setFilter(bool f) { filter = f; } virtual void setCorrection(bool c) { correction = c; } virtual void setProcessors(int p) { processors = p; } virtual void setWindow(int w) { window = w; } virtual void setIncrement(int i) { increment = i; } - virtual void setTemplate(string t) { templateFile = t; } + + virtual void setCons(string) {}; + //pure functions virtual void getChimeras() = 0; @@ -42,8 +45,7 @@ class Chimera { bool filter, correction; int processors, window, increment; - string templateFile; - + }; diff --git a/chimeraseqscommand.cpp b/chimeraseqscommand.cpp index 105fd11..28e3e6a 100644 --- a/chimeraseqscommand.cpp +++ b/chimeraseqscommand.cpp @@ -22,7 +22,7 @@ ChimeraSeqsCommand::ChimeraSeqsCommand(string option){ else { //valid paramters for this command - string Array[] = {"fasta", "filter", "correction", "processors", "method", "window", "increment", "template" }; + string Array[] = {"fasta", "filter", "correction", "processors", "method", "window", "increment", "template", "conservation" }; vector myArray (Array, Array+(sizeof(Array)/sizeof(string))); OptionParser parser(option); @@ -44,6 +44,10 @@ ChimeraSeqsCommand::ChimeraSeqsCommand(string option){ if (templatefile == "not open") { abort = true; } else if (templatefile == "not found") { templatefile = ""; } + consfile = validParameter.validFile(parameters, "conservation", true); + if (consfile == "not open") { abort = true; } + else if (consfile == "not found") { consfile = ""; } + string temp; temp = validParameter.validFile(parameters, "filter", false); if (temp == "not found") { temp = "T"; } @@ -58,12 +62,12 @@ ChimeraSeqsCommand::ChimeraSeqsCommand(string option){ temp = validParameter.validFile(parameters, "window", false); if (temp == "not found") { temp = "0"; } convert(temp, window); - temp = validParameter.validFile(parameters, "increment", false); if (temp == "not found") { temp = "10"; } + temp = validParameter.validFile(parameters, "increment", false); if (temp == "not found") { temp = "25"; } convert(temp, increment); method = validParameter.validFile(parameters, "method", false); if (method == "not found") { method = "pintail"; } - if ((method == "pintail") && (templatefile == "")) { mothurOut("You must provide a template file with the pintail method."); mothurOutEndLine(); abort = true; } + if ((method == "pintail") && (templatefile == "") && (consfile == "")) { mothurOut("You must provide a template or conservation file with the pintail method."); mothurOutEndLine(); abort = true; } } @@ -105,9 +109,11 @@ int ChimeraSeqsCommand::execute(){ if (abort == true) { return 0; } - if (method == "bellerophon") { chimera = new Bellerophon(fastafile); } - else if (method == "pintail") { chimera = new Pintail(fastafile); } - else { mothurOut("Not a valid method."); mothurOutEndLine(); return 0; } + if (method == "bellerophon") { chimera = new Bellerophon(fastafile); } + else if (method == "pintail") { chimera = new Pintail(fastafile, templatefile); + if (consfile != "") { chimera->setCons(consfile); } + else { chimera->setCons(""); } + }else { mothurOut("Not a valid method."); mothurOutEndLine(); return 0; } //set user options chimera->setFilter(filter); @@ -115,8 +121,7 @@ int ChimeraSeqsCommand::execute(){ chimera->setProcessors(processors); chimera->setWindow(window); chimera->setIncrement(increment); - chimera->setTemplate(templatefile); - + //find chimeras chimera->getChimeras(); @@ -139,6 +144,5 @@ int ChimeraSeqsCommand::execute(){ exit(1); } } - /**************************************************************************************************/ diff --git a/chimeraseqscommand.h b/chimeraseqscommand.h index dc7df8f..9bc8418 100644 --- a/chimeraseqscommand.h +++ b/chimeraseqscommand.h @@ -28,7 +28,7 @@ public: private: bool abort; - string method, fastafile, templatefile; + string method, fastafile, templatefile, consfile; bool filter, correction; int processors, midpoint, averageLeft, averageRight, window, iters, increment; Chimera* chimera; diff --git a/database.cpp b/database.cpp index 5fc9b67..534747a 100644 --- a/database.cpp +++ b/database.cpp @@ -15,7 +15,8 @@ Database::Database(string fastaFileName){ // This assumes that the template database is in fasta format, may // need to alter this in the future? - + longest = 0; + ifstream fastaFile; openInputFile(fastaFileName, fastaFile); @@ -43,6 +44,13 @@ Database::Database(string fastaFileName){ // This assumes that the template dat } } templateSequences[i] = Sequence(seqName, aligned); + + //necessary because Sequence constructor by default sets whatever it reads to unaligned + //the setUnaligned function remove gap char. + templateSequences[i].setUnaligned(templateSequences[i].getUnaligned()); + + if (templateSequences[i].getUnaligned().length() > longest) { longest = templateSequences[i].getUnaligned().length(); } + fastaFile.putback(letter); } @@ -64,4 +72,10 @@ Database::~Database(){ float Database::getSearchScore() { return searchScore; } // we're assuming that the search is already done + /**************************************************************************************************/ + +int Database::getLongestBase() { return longest+1; } + +/**************************************************************************************************/ + diff --git a/database.hpp b/database.hpp index 02cad03..01393c1 100644 --- a/database.hpp +++ b/database.hpp @@ -23,8 +23,10 @@ public: virtual ~Database(); virtual Sequence findClosestSequence(Sequence*) = 0; virtual float getSearchScore(); + virtual int getLongestBase(); + protected: - int numSeqs; + int numSeqs, longest; float searchScore; vector templateSequences; }; diff --git a/engine.cpp b/engine.cpp index 2215e08..95c0ff0 100644 --- a/engine.cpp +++ b/engine.cpp @@ -244,13 +244,21 @@ bool ScriptEngine::getInput(){ string ScriptEngine::getNextCommand(string& commandString) { try { string nextcommand = ""; + int count = 0; - nextcommand = commandString.substr(0,commandString.find_first_of(';')); - + //go through string until you reach ; or end + while (count < commandString.length()) { + + if (commandString[count] == ';') { break; } + else { nextcommand += commandString[count]; } + + count++; + } + + //if you are not at the end + if (count != commandString.length()) { commandString = commandString.substr(count+1, commandString.length()); } + else { commandString = ""; } - if ((commandString.find_first_of(';')+1) <= commandString.length()) { - commandString = commandString.substr(commandString.find_first_of(';')+1, commandString.length()); - }else { commandString = ""; } //you have reached the last command. //get rid of spaces in between commands if any if (commandString.length() > 0) { diff --git a/gotohoverlap.cpp b/gotohoverlap.cpp index e3be6c3..f6fef1d 100644 --- a/gotohoverlap.cpp +++ b/gotohoverlap.cpp @@ -28,62 +28,75 @@ GotohOverlap::GotohOverlap(float gO, float gE, float m, float mm, int r) : gapOpen(gO), gapExtend(gE), match(m), mismatch(mm), Alignment(r) { - for(int i=1;i alignment[i][j].dValue){ - if(alignment[i][j].iValue > diagonal){ - alignment[i][j].cValue = alignment[i][j].iValue; - alignment[i][j].prevCell = 'l'; + try { + seqA = ' ' + A; lA = seqA.length(); // the algorithm requires that the first character be a dummy value + seqB = ' ' + B; lB = seqB.length(); // the algorithm requires that the first character be a dummy value + + for(int i=1;i alignment[i][j].dValue){ + if(alignment[i][j].iValue > diagonal){ + alignment[i][j].cValue = alignment[i][j].iValue; + alignment[i][j].prevCell = 'l'; + } + else{ + alignment[i][j].cValue = diagonal; + alignment[i][j].prevCell = 'd'; + } } else{ - alignment[i][j].cValue = diagonal; - alignment[i][j].prevCell = 'd'; + if(alignment[i][j].dValue > diagonal){ + alignment[i][j].cValue = alignment[i][j].dValue; + alignment[i][j].prevCell = 'u'; + } + else{ + alignment[i][j].cValue = diagonal; + alignment[i][j].prevCell = 'd'; + } } + } - else{ - if(alignment[i][j].dValue > diagonal){ - alignment[i][j].cValue = alignment[i][j].dValue; - alignment[i][j].prevCell = 'u'; - } - else{ - alignment[i][j].cValue = diagonal; - alignment[i][j].prevCell = 'd'; - } - } - } + Overlap over; + over.setOverlap(alignment, lA, lB, 0); // Fix the gaps at the ends of the sequences + traceBack(); // Construct the alignment and set seqAaln and seqBaln + + } + catch(exception& e) { + errorOut(e, "GotohOverlap", "align"); + exit(1); } - Overlap over; - over.setOverlap(alignment, lA, lB, 0); // Fix the gaps at the ends of the sequences - traceBack(); // Construct the alignment and set seqAaln and seqBaln } /**************************************************************************************************/ diff --git a/gotohoverlap.hpp b/gotohoverlap.hpp index 85176ef..0b7d8ea 100644 --- a/gotohoverlap.hpp +++ b/gotohoverlap.hpp @@ -30,6 +30,7 @@ class GotohOverlap : public Alignment { public: GotohOverlap(float, float, float, float, int); void align(string, string); + ~GotohOverlap() {} private: diff --git a/kmerdb.cpp b/kmerdb.cpp index 0ec822f..86dba09 100644 --- a/kmerdb.cpp +++ b/kmerdb.cpp @@ -30,24 +30,31 @@ /**************************************************************************************************/ KmerDB::KmerDB(string fastaFileName, int kSize) : Database(fastaFileName), kmerSize(kSize) { - - string kmerDBName = fastaFileName.substr(0,fastaFileName.find_last_of(".")+1) + char('0'+ kmerSize) + "mer"; - ifstream kmerFileTest(kmerDBName.c_str()); - - int power4s[14] = { 1, 4, 16, 64, 256, 1024, 4096, 16384, 65536, 262144, 1048576, 4194304, 16777216, 67108864 }; - - maxKmer = power4s[kmerSize]; - kmerLocations.resize(maxKmer+1); + try { - if(!kmerFileTest){ // if we can open the kmer db file, then read it in... - mothurOut("Generating the " + kmerDBName + " database...\t"); cout.flush(); - generateKmerDB(kmerDBName); - } - else{ // ...otherwise generate it. - mothurOut("Reading in the " + kmerDBName + " database...\t"); cout.flush(); - readKmerDB(kmerDBName, kmerFileTest); + string kmerDBName = fastaFileName.substr(0,fastaFileName.find_last_of(".")+1) + char('0'+ kmerSize) + "mer"; + ifstream kmerFileTest(kmerDBName.c_str()); + + int power4s[14] = { 1, 4, 16, 64, 256, 1024, 4096, 16384, 65536, 262144, 1048576, 4194304, 16777216, 67108864 }; + + maxKmer = power4s[kmerSize]; + kmerLocations.resize(maxKmer+1); + + if(!kmerFileTest){ // if we can open the kmer db file, then read it in... + mothurOut("Generating the " + kmerDBName + " database...\t"); cout.flush(); + generateKmerDB(kmerDBName); + } + else{ // ...otherwise generate it. + mothurOut("Reading in the " + kmerDBName + " database...\t"); cout.flush(); + readKmerDB(kmerDBName, kmerFileTest); + } + mothurOut("DONE."); mothurOutEndLine(); mothurOutEndLine(); cout.flush(); + } - mothurOut("DONE."); mothurOutEndLine(); mothurOutEndLine(); cout.flush(); + catch(exception& e) { + errorOut(e, "KmerDB", "KmerDB"); + exit(1); + } } /**************************************************************************************************/ @@ -61,92 +68,114 @@ KmerDB::~KmerDB(){ /**************************************************************************************************/ Sequence KmerDB::findClosestSequence(Sequence* candidateSeq){ + try { - Kmer kmer(kmerSize); - - searchScore = 0; - int maxSequence = 0; - - vector matches(numSeqs, 0); // a record of the sequences with shared kmers - vector timesKmerFound(kmerLocations.size()+1, 0); // a record of the kmers that we have already found - - int numKmers = candidateSeq->getNumBases() - kmerSize + 1; - - for(int i=0;igetUnaligned(), i); // go through the query sequence and get a kmer number - if(timesKmerFound[kmerNumber] == 0){ // if we haven't seen it before... - for(int j=0;j matches(numSeqs, 0); // a record of the sequences with shared kmers + vector timesKmerFound(kmerLocations.size()+1, 0); // a record of the kmers that we have already found + + int numKmers = candidateSeq->getNumBases() - kmerSize + 1; + + for(int i=0;igetUnaligned(), i); // go through the query sequence and get a kmer number + if(timesKmerFound[kmerNumber] == 0){ // if we haven't seen it before... + for(int j=0;j searchScore){ // the query sequence - searchScore = matches[i]; - maxSequence = i; + + for(int i=0;i searchScore){ // the query sequence + searchScore = matches[i]; + maxSequence = i; + } } + + searchScore = 100 * searchScore / (float) numKmers; // return the Sequence object corresponding to the db + + return templateSequences[maxSequence]; // sequence with the most shared kmers. + } - - searchScore = 100 * searchScore / (float) numKmers; // return the Sequence object corresponding to the db - return templateSequences[maxSequence]; // sequence with the most shared kmers. + catch(exception& e) { + errorOut(e, "KmerDB", "findClosestSequence"); + exit(1); + } } /**************************************************************************************************/ void KmerDB::generateKmerDB(string kmerDBName){ + try { - Kmer kmer(kmerSize); - - for(int i=0;i seenBefore(maxKmer+1,0); - for(int j=0;j seenBefore(maxKmer+1,0); + for(int j=0;j> seqName >> numValues; + kmerDBFile.seekg(0); // start at the beginning of the file + + string seqName; + int seqNumber; - for(int j=0;j> seqNumber; // 1. number of sequences with the kmer number - kmerLocations[i].push_back(seqNumber); // 2. sequence indices + for(int i=0;i> seqName >> numValues; + + for(int j=0;j> seqNumber; // 1. number of sequences with the kmer number + kmerLocations[i].push_back(seqNumber); // 2. sequence indices + } } + kmerDBFile.close(); + } - kmerDBFile.close(); + catch(exception& e) { + errorOut(e, "KmerDB", "readKmerDB"); + exit(1); + } } /**************************************************************************************************/ diff --git a/nast.cpp b/nast.cpp index 0741ccd..4aed590 100644 --- a/nast.cpp +++ b/nast.cpp @@ -22,50 +22,70 @@ /**************************************************************************************************/ Nast::Nast(Alignment* method, Sequence* cand, Sequence* temp) : alignment(method), candidateSeq(cand), templateSeq(temp) { - maxInsertLength = 0; - pairwiseAlignSeqs(); // This is part A in Fig. 2 of DeSantis et al. - regapSequences(); // This is parts B-F in Fig. 2 of DeSantis et al. + try { + maxInsertLength = 0; + pairwiseAlignSeqs(); // This is part A in Fig. 2 of DeSantis et al. + regapSequences(); // This is parts B-F in Fig. 2 of DeSantis et al. + + } + catch(exception& e) { + errorOut(e, "Nast", "Nast"); + exit(1); + } + } /**************************************************************************************************/ void Nast::pairwiseAlignSeqs(){ // Here we call one of the pairwise alignment methods to align our unaligned candidate // and template sequences - alignment->align(candidateSeq->getUnaligned(), templateSeq->getUnaligned()); + try { + + alignment->align(candidateSeq->getUnaligned(), templateSeq->getUnaligned()); - string candAln = alignment->getSeqAAln(); - string tempAln = alignment->getSeqBAln(); + string candAln = alignment->getSeqAAln(); + string tempAln = alignment->getSeqBAln(); - if(candAln == ""){ - candidateSeq->setPairwise(""); - templateSeq->setPairwise(templateSeq->getUnaligned()); - } - else{ - if(tempAln[0] == '-'){ - int pairwiseAlignmentLength = tempAln.length(); // we need to make sure that the candidate sequence alignment - for(int i=0;isetPairwise(""); + templateSeq->setPairwise(templateSeq->getUnaligned()); + + } + else{ + + if(tempAln[0] == '-'){ + int pairwiseAlignmentLength = tempAln.length(); // we need to make sure that the candidate sequence alignment + for(int i=0;i=0; i--){// ends where the template sequence alignment ends, if it runs - if(isalpha(tempAln[i])){ // long, we nuke the end of the candidate sequence - candAln = candAln.substr(0,i+1); - tempAln = tempAln.substr(0,i+1); - break; - } + + int pairwiseAlignmentLength = tempAln.length(); + if(tempAln[pairwiseAlignmentLength-1] == '-'){ // we need to make sure that the candidate sequence alignment + for(int i=pairwiseAlignmentLength-1; i>=0; i--){// ends where the template sequence alignment ends, if it runs + if(isalpha(tempAln[i])){ // long, we nuke the end of the candidate sequence + candAln = candAln.substr(0,i+1); + tempAln = tempAln.substr(0,i+1); + break; + } + } } + } + + candidateSeq->setPairwise(candAln); // set the pairwise sequences in the Sequence objects for + templateSeq->setPairwise(tempAln); // the candidate and template sequences + } - candidateSeq->setPairwise(candAln); // set the pairwise sequences in the Sequence objects for - templateSeq->setPairwise(tempAln); // the candidate and template sequences - + catch(exception& e) { + errorOut(e, "Nast", "pairwiseAlignSeqs"); + exit(1); + } } /**************************************************************************************************/ @@ -73,274 +93,301 @@ void Nast::pairwiseAlignSeqs(){ // Here we call one of the pairwise alignment me void Nast::removeExtraGaps(string& candAln, string tempAln, string newTemplateAlign){ // here we do steps C-F of Fig. 2 from DeSantis et al. + try { - int longAlignmentLength = newTemplateAlign.length(); - - for(int i=0; i=0;leftIndex--){ // then we've got problems... - if(!isalpha(candAln[leftIndex])){ - leftRoom = 1; //count how far it is to the nearest gap on the LEFT side of the anomaly - while(leftIndex-leftRoom>=0 && !isalpha(candAln[leftIndex-leftRoom])) { leftRoom++; } - break; - } - } - - - for(rightIndex=i+1;rightIndex maxInsertLength){ maxInsertLength = insertLength; } + int longAlignmentLength = newTemplateAlign.length(); + + for(int i=0; i= insertLength){ - // Parts D & E from Fig. 2 of DeSantis et al. - if((i-leftIndex) <= (rightIndex-i)){ // the left gap is closer - > move stuff left there's - if(leftRoom >= insertLength){ // enough room to the left to move - string leftTemplateString = newTemplateAlign.substr(0,i); - string rightTemplateString = newTemplateAlign.substr(i+insertLength); - newTemplateAlign = leftTemplateString + rightTemplateString; - longAlignmentLength = newTemplateAlign.length(); - - string leftCandidateString = candAln.substr(0,leftIndex-insertLength+1); - string rightCandidateString = candAln.substr(leftIndex+1); - candAln = leftCandidateString + rightCandidateString; + // Part C of Fig. 2 from DeSantis et al. + if((isalpha(newTemplateAlign[i]) != isalpha(tempAln[i]))){ //if there is a discrepancy between the regapped + + rightRoom = 0; leftRoom = 0; + + // Part D of Fig. 2 from DeSantis et al. // template sequence and the official template sequence + for(leftIndex=i-1;leftIndex>0;leftIndex--){ // then we've got problems... + if(!isalpha(candAln[leftIndex])){ + leftRoom = 1; //count how far it is to the nearest gap on the LEFT side of the anomaly + while(leftIndex-leftRoom>=0 && !isalpha(candAln[leftIndex-leftRoom])) { leftRoom++; } + break; } - else{ // not enough room to the left, have to steal some space to - string leftTemplateString = newTemplateAlign.substr(0,i); // the right - string rightTemplateString = newTemplateAlign.substr(i+insertLength); - newTemplateAlign = leftTemplateString + rightTemplateString; - longAlignmentLength = newTemplateAlign.length(); - - string leftCandidateString = candAln.substr(0,leftIndex-leftRoom+1); - string insertString = candAln.substr(leftIndex+1,rightIndex-leftIndex-1); - string rightCandidateString = candAln.substr(rightIndex+(insertLength-leftRoom)); - candAln = leftCandidateString + insertString + rightCandidateString; + } + + + for(rightIndex=i+1;rightIndex move stuff right there's - if(rightRoom >= insertLength){ // enough room to the right to move - string leftTemplateString = newTemplateAlign.substr(0,i); - string rightTemplateString = newTemplateAlign.substr(i+insertLength); - newTemplateAlign = leftTemplateString + rightTemplateString; - longAlignmentLength = newTemplateAlign.length(); - - string leftCandidateString = candAln.substr(0,rightIndex); - string rightCandidateString = candAln.substr(rightIndex+insertLength); - candAln = leftCandidateString + rightCandidateString; - + + int insertLength = 0; // figure out how long the anomaly is + while(!isalpha(newTemplateAlign[i + insertLength])) { insertLength++; } + if(insertLength > maxInsertLength){ maxInsertLength = insertLength; } + + if((leftRoom + rightRoom) >= insertLength){ + + // Parts D & E from Fig. 2 of DeSantis et al. + if((i-leftIndex) <= (rightIndex-i)){ // the left gap is closer - > move stuff left there's + + if(leftRoom >= insertLength){ // enough room to the left to move + + string leftTemplateString = newTemplateAlign.substr(0,i); + string rightTemplateString = newTemplateAlign.substr(i+insertLength); + newTemplateAlign = leftTemplateString + rightTemplateString; + longAlignmentLength = newTemplateAlign.length(); + + string leftCandidateString = candAln.substr(0,leftIndex-insertLength+1); + string rightCandidateString = candAln.substr(leftIndex+1); + candAln = leftCandidateString + rightCandidateString; + + } + else{ // not enough room to the left, have to steal some space to + string leftTemplateString = newTemplateAlign.substr(0,i); // the right + string rightTemplateString = newTemplateAlign.substr(i+insertLength); + newTemplateAlign = leftTemplateString + rightTemplateString; + longAlignmentLength = newTemplateAlign.length(); + + string leftCandidateString = candAln.substr(0,leftIndex-leftRoom+1); + string insertString = candAln.substr(leftIndex+1,rightIndex-leftIndex-1); + string rightCandidateString = candAln.substr(rightIndex+(insertLength-leftRoom)); + candAln = leftCandidateString + insertString + rightCandidateString; + + } } - else{ // not enough room to the right, have to steal some - // space to the left lets move left and then right... - string leftTemplateString = newTemplateAlign.substr(0,i); - string rightTemplateString = newTemplateAlign.substr(i+insertLength); - newTemplateAlign = leftTemplateString + rightTemplateString; - longAlignmentLength = newTemplateAlign.length(); - - string leftCandidateString = candAln.substr(0,leftIndex-(insertLength-rightRoom)+1); - string insertString = candAln.substr(leftIndex+1,rightIndex-leftIndex-1); - string rightCandidateString = candAln.substr(rightIndex+rightRoom); - candAln = leftCandidateString + insertString + rightCandidateString; - + else{ // the right gap is closer - > move stuff right there's + if(rightRoom >= insertLength){ // enough room to the right to move + + string leftTemplateString = newTemplateAlign.substr(0,i); + string rightTemplateString = newTemplateAlign.substr(i+insertLength); + newTemplateAlign = leftTemplateString + rightTemplateString; + longAlignmentLength = newTemplateAlign.length(); + + string leftCandidateString = candAln.substr(0,rightIndex); + string rightCandidateString = candAln.substr(rightIndex+insertLength); + candAln = leftCandidateString + rightCandidateString; + + } + else{ // not enough room to the right, have to steal some + // space to the left lets move left and then right... + string leftTemplateString = newTemplateAlign.substr(0,i); + string rightTemplateString = newTemplateAlign.substr(i+insertLength); + newTemplateAlign = leftTemplateString + rightTemplateString; + longAlignmentLength = newTemplateAlign.length(); + + string leftCandidateString = candAln.substr(0,leftIndex-(insertLength-rightRoom)+1); + string insertString = candAln.substr(leftIndex+1,rightIndex-leftIndex-1); + string rightCandidateString = candAln.substr(rightIndex+rightRoom); + candAln = leftCandidateString + insertString + rightCandidateString; + + } } } - } - else{ // there could be a case where there isn't enough room in - string leftTemplateString = newTemplateAlign.substr(0,i); // either direction to move stuff - string rightTemplateString = newTemplateAlign.substr(i+leftRoom+rightRoom); - newTemplateAlign = leftTemplateString + rightTemplateString; - longAlignmentLength = newTemplateAlign.length(); + else{ + // there could be a case where there isn't enough room in + string leftTemplateString = newTemplateAlign.substr(0,i); // either direction to move stuff + string rightTemplateString = newTemplateAlign.substr(i+leftRoom+rightRoom); + newTemplateAlign = leftTemplateString + rightTemplateString; + longAlignmentLength = newTemplateAlign.length(); + + string leftCandidateString = candAln.substr(0,leftIndex-leftRoom+1); + string insertString = candAln.substr(leftIndex+1,rightIndex-leftIndex-1); + string rightCandidateString = candAln.substr(rightIndex+rightRoom); + candAln = leftCandidateString + insertString + rightCandidateString; + + } - string leftCandidateString = candAln.substr(0,leftIndex-leftRoom+1); - string insertString = candAln.substr(leftIndex+1,rightIndex-leftIndex-1); - string rightCandidateString = candAln.substr(rightIndex+rightRoom); - candAln = leftCandidateString + insertString + rightCandidateString; - } - - i -= insertLength; - } + i -= insertLength; + } + } } - + catch(exception& e) { + errorOut(e, "Nast", "removeExtraGaps"); + exit(1); + } + } /**************************************************************************************************/ void Nast::regapSequences(){ //This is essentially part B in Fig 2. of DeSantis et al. + try { - string candPair = candidateSeq->getPairwise(); - string candAln = ""; - - string tempPair = templateSeq->getPairwise(); - string tempAln = templateSeq->getAligned(); // we use the template aligned sequence as our guide - - int pairwiseLength = candPair.length(); - int fullAlignLength = tempAln.length(); - - if(candPair == ""){ - for(int i=0;isetAligned(candAln); - return; - } - - int fullAlignIndex = 0; - int pairwiseAlignIndex = 0; - string newTemplateAlign = ""; // this is going to be messy so we want a temporary template - // alignment string - while(tempAln[fullAlignIndex] == '.' || tempAln[fullAlignIndex] == '-'){ - candAln += '.'; // add the initial '-' and '.' to the candidate and template - newTemplateAlign += tempAln[fullAlignIndex];// pairwise sequences - fullAlignIndex++; - } - - string lastLoop = ""; - - while(pairwiseAlignIndexgetPairwise(); + string candAln = ""; + + string tempPair = templateSeq->getPairwise(); + string tempAln = templateSeq->getAligned(); // we use the template aligned sequence as our guide + + int pairwiseLength = candPair.length(); + int fullAlignLength = tempAln.length(); + + if(candPair == ""){ + for(int i=0;isetAligned(candAln); + return; } - else if(isalpha(tempPair[pairwiseAlignIndex]) && !isalpha(tempAln[fullAlignIndex]) - && isalpha(candPair[pairwiseAlignIndex])){ - // the template pairwise and candidate pairwise are characters and the template aligned is a gap - // need to insert gaps into the candidateSeq.aligned sequence - - candAln += '-'; - newTemplateAlign += '-';// + + int fullAlignIndex = 0; + int pairwiseAlignIndex = 0; + string newTemplateAlign = ""; // this is going to be messy so we want a temporary template + // alignment string + while(tempAln[fullAlignIndex] == '.' || tempAln[fullAlignIndex] == '-'){ + candAln += '.'; // add the initial '-' and '.' to the candidate and template + newTemplateAlign += tempAln[fullAlignIndex];// pairwise sequences fullAlignIndex++; } - else if(!isalpha(tempPair[pairwiseAlignIndex]) && isalpha(tempAln[fullAlignIndex]) - && isalpha(candPair[pairwiseAlignIndex])){ - // the template pairwise is a gap and the template aligned and pairwise sequences have characters - // this is the alpha scenario. add character to the candidateSeq.aligned sequence without progressing - // further through the tempAln sequence. - - candAln += candPair[pairwiseAlignIndex]; - newTemplateAlign += '-';// - pairwiseAlignIndex++; - } - else if(isalpha(tempPair[pairwiseAlignIndex]) && isalpha(tempAln[fullAlignIndex]) - && !isalpha(candPair[pairwiseAlignIndex])){ - // the template pairwise and full alignment are characters and the candidate sequence has a gap - // should not be a big deal, just add the gap position to the candidateSeq.aligned sequence; - - candAln += candPair[pairwiseAlignIndex]; - newTemplateAlign += tempAln[fullAlignIndex];// - fullAlignIndex++; - pairwiseAlignIndex++; + + string lastLoop = ""; + + while(pairwiseAlignIndexseems to be the opposite of the alpha scenario + + candAln += candPair[pairwiseAlignIndex]; + newTemplateAlign += tempAln[fullAlignIndex];// + pairwiseAlignIndex++; + fullAlignIndex++; + } + else if(isalpha(tempPair[pairwiseAlignIndex]) && !isalpha(tempAln[fullAlignIndex]) + && !isalpha(candPair[pairwiseAlignIndex])){ + // template pairwise has a character, but its full aligned sequence and candidate sequence have gaps + // this would happen like we need to add a gap. basically the opposite of the alpha situation + + newTemplateAlign += tempAln[fullAlignIndex];// + candAln += "-"; + fullAlignIndex++; + } + else if(!isalpha(tempPair[pairwiseAlignIndex]) && isalpha(tempAln[fullAlignIndex]) + && !isalpha(candPair[pairwiseAlignIndex])){ + // template and candidate pairwise are gaps and the template aligned is not a gap this should not be possible + // would skip the gaps and not progress through full alignment sequence + // not tested yet + + mothurOut("We're into D " + toString(fullAlignIndex) + " " + toString(pairwiseAlignIndex)); mothurOutEndLine(); + pairwiseAlignIndex++; + } + else{ + // everything has a gap - not possible + // not tested yet + + mothurOut("We're into F " + toString(fullAlignIndex) + " " + toString(pairwiseAlignIndex)); mothurOutEndLine(); + pairwiseAlignIndex++; + fullAlignIndex++; + } } - else if(!isalpha(tempPair[pairwiseAlignIndex]) && !isalpha(tempAln[fullAlignIndex]) - && isalpha(candPair[pairwiseAlignIndex])){ - // the template pairwise and aligned are gaps while the candidate pairwise has a character - // this would be an insertion, go ahead and add the character->seems to be the opposite of the alpha scenario - - candAln += candPair[pairwiseAlignIndex]; - newTemplateAlign += tempAln[fullAlignIndex];// - pairwiseAlignIndex++; - fullAlignIndex++; + + for(int i=fullAlignIndex;i=0;i--){ // ditto. + if(candAln[i] == 'Z' || !isalnum(candAln[i])) { candAln[i] = '.'; } + else{ end = i; break; } } - else if(!isalpha(tempPair[pairwiseAlignIndex]) && isalpha(tempAln[fullAlignIndex]) - && !isalpha(candPair[pairwiseAlignIndex])){ - // template and candidate pairwise are gaps and the template aligned is not a gap this should not be possible - // would skip the gaps and not progress through full alignment sequence - // not tested yet - - mothurOut("We're into D " + toString(fullAlignIndex) + " " + toString(pairwiseAlignIndex)); mothurOutEndLine(); - pairwiseAlignIndex++; + + for(int i=start;i<=end;i++){ // go through the candidate alignment sequence and make sure that + candAln[i] = toupper(candAln[i]); // everything is upper case } - else{ - // everything has a gap - not possible - // not tested yet - - mothurOut("We're into F " + toString(fullAlignIndex) + " " + toString(pairwiseAlignIndex)); mothurOutEndLine(); - pairwiseAlignIndex++; - fullAlignIndex++; - } - } - - for(int i=fullAlignIndex;i=0;i--){ // ditto. - if(candAln[i] == 'Z' || !isalnum(candAln[i])) { candAln[i] = '.'; } - else{ end = i; break; } - } - - for(int i=start;i<=end;i++){ // go through the candidate alignment sequence and make sure that - candAln[i] = toupper(candAln[i]); // everything is upper case + + + if(candAln.length() != tempAln.length()){ // if the regapped candidate sequence is longer than the official + removeExtraGaps(candAln, tempAln, newTemplateAlign);// template alignment then we need to do steps C-F in Fig. + } // 2 of Desantis et al. + + + candidateSeq->setAligned(candAln); } - - - if(candAln.length() != tempAln.length()){ // if the regapped candidate sequence is longer than the official - removeExtraGaps(candAln, tempAln, newTemplateAlign);// template alignment then we need to do steps C-F in Fig. - } // 2 of Desantis et al. - - - candidateSeq->setAligned(candAln); + catch(exception& e) { + errorOut(e, "Nast", "regapSequences"); + exit(1); + } } /**************************************************************************************************/ float Nast::getSimilarityScore(){ - - string cand = candidateSeq->getAligned(); - string temp = templateSeq->getAligned(); - int alignmentLength = temp.length(); - int mismatch = 0; - int denominator = 0; + try { - for(int i=0;igetAligned(); + string temp = templateSeq->getAligned(); + int alignmentLength = temp.length(); + int mismatch = 0; + int denominator = 0; + + for(int i=0;i= up){ - if(diagonal >= left){ - alignment[i][j].cValue = diagonal; - alignment[i][j].prevCell = 'd'; - } - else{ - alignment[i][j].cValue = left; - alignment[i][j].prevCell = 'l'; - } - } - else{ - if(up >= left){ - alignment[i][j].cValue = up; - alignment[i][j].prevCell = 'u'; + try { + seqA = ' ' + A; lA = seqA.length(); // algorithm requires a dummy space at the beginning of each string + seqB = ' ' + B; lB = seqB.length(); // algorithm requires a dummy space at the beginning of each string + + if (lA > nRows) { mothurOut("Your one of your candidate sequences is longer than you longest template sequence."); mothurOutEndLine(); } + + for(int i=1;i= up){ + if(diagonal >= left){ + alignment[i][j].cValue = diagonal; + alignment[i][j].prevCell = 'd'; + } + else{ + alignment[i][j].cValue = left; + alignment[i][j].prevCell = 'l'; + } } else{ - alignment[i][j].cValue = left; - alignment[i][j].prevCell = 'l'; + if(up >= left){ + alignment[i][j].cValue = up; + alignment[i][j].prevCell = 'u'; + } + else{ + alignment[i][j].cValue = left; + alignment[i][j].prevCell = 'l'; + } } } } + Overlap over; + over.setOverlap(alignment, lA, lB, 0); // Fix gaps at the beginning and end of the sequences + traceBack(); // Traceback the alignment to populate seqAaln and seqBaln + + } + catch(exception& e) { + errorOut(e, "NeedlemanOverlap", "align"); + exit(1); } - Overlap over; - over.setOverlap(alignment, lA, lB, 0); // Fix gaps at the beginning and end of the sequences - traceBack(); // Traceback the alignment to populate seqAaln and seqBaln + } /**************************************************************************************************/ diff --git a/pintail.cpp b/pintail.cpp index 4b812a7..5f4abd0 100644 --- a/pintail.cpp +++ b/pintail.cpp @@ -8,19 +8,11 @@ */ #include "pintail.h" -#include "ignoregaps.h" +#include "eachgapdist.h" //*************************************************************************************************************** -Pintail::Pintail(string name) { - try { - fastafile = name; - } - catch(exception& e) { - errorOut(e, "Pintail", "Pintail"); - exit(1); - } -} +Pintail::Pintail(string filename, string temp) { fastafile = filename; templateFile = temp; } //*************************************************************************************************************** Pintail::~Pintail() { @@ -37,24 +29,20 @@ Pintail::~Pintail() { void Pintail::print(ostream& out) { try { - for (itCoef = DE.begin(); itCoef != DE.end(); itCoef++) { + for (int i = 0; i < querySeqs.size(); i++) { - out << itCoef->first->getName() << '\t' << itCoef->second << endl; + out << querySeqs[i]->getName() << '\t' << "div: " << deviation[i] << "\tstDev: " << DE[i] << endl; out << "Observed\t"; - itObsDist = obsDistance.find(itCoef->first); - for (int i = 0; i < itObsDist->second.size(); i++) { out << itObsDist->second[i] << '\t'; } + for (int j = 0; j < obsDistance[i].size(); j++) { out << obsDistance[i][j] << '\t'; } out << endl; out << "Expected\t"; - itExpDist = expectedDistance.find(itCoef->first); - for (int i = 0; i < itExpDist->second.size(); i++) { out << itExpDist->second[i] << '\t'; } + for (int m = 0; m < expectedDistance[i].size(); m++) { out << expectedDistance[i][m] << '\t'; } out << endl; } - - } catch(exception& e) { errorOut(e, "Pintail", "print"); @@ -66,8 +54,6 @@ void Pintail::print(ostream& out) { void Pintail::getChimeras() { try { - distCalculator = new ignoreGaps(); - //read in query sequences and subject sequences mothurOut("Reading sequences and template file... "); cout.flush(); querySeqs = readSeqs(fastafile); @@ -76,71 +62,109 @@ void Pintail::getChimeras() { int numSeqs = querySeqs.size(); - //if window is set to default - if (window == 0) { if (querySeqs[0]->getAligned().length() > 800) { setWindow(200); } - else{ setWindow((querySeqs[0]->getAligned().length() / 4)); } } - else if (window > (querySeqs[0]->getAligned().length() / 4)) { - mothurOut("You have selected to large a window size for you sequences. I will choose a smaller window."); mothurOutEndLine(); - setWindow((querySeqs[0]->getAligned().length() / 4)); - } - - //calculate number of iters - iters = (querySeqs[0]->getAligned().length() - window + 1) / increment; -cout << "length = " << querySeqs[0]->getAligned().length() << " window = " << window << " increment = " << increment << " iters = " << iters << endl; + obsDistance.resize(numSeqs); + expectedDistance.resize(numSeqs); + seqCoef.resize(numSeqs); + DE.resize(numSeqs); + Qav.resize(numSeqs); + bestfit.resize(numSeqs); + trim.resize(numSeqs); + deviation.resize(numSeqs); + windowSizes.resize(numSeqs, window); + + //break up file if needed int linesPerProcess = processors / numSeqs; - //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))); + #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)); } - //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)); - } - - //map query sequences to their most similiar sequences in the template - Parallelized - mothurOut("Finding closest sequence in template to each sequence... "); cout.flush(); - if (processors == 1) { findPairs(lines[0]->start, lines[0]->end); } - else { createProcessesPairs(); } - mothurOut("Done."); mothurOutEndLine(); - - //find Oqs for each sequence - the observed distance in each window - Parallelized - mothurOut("Calculating observed percentage differences for each sequence... "); cout.flush(); - if (processors == 1) { calcObserved(lines[0]->start, lines[0]->end); } - else { createProcessesObserved(); } - mothurOut("Done."); mothurOutEndLine(); - + #else + lines.push_back(new linePair(0, numSeqs)); + #endif + + distcalculator = new eachGapDist(); + + if (processors == 1) { + mothurOut("Finding closest sequence in template to each sequence... "); cout.flush(); + bestfit = findPairs(lines[0]->start, lines[0]->end); +for (int m = 0; m < templateSeqs.size(); m++) { + if (templateSeqs[m]->getName() == "198806") { bestfit[0] = *(templateSeqs[m]); } + if (templateSeqs[m]->getName() == "198806") { bestfit[1] = *(templateSeqs[m]); } + if (templateSeqs[m]->getName() == "108139") { bestfit[2] = *(templateSeqs[m]); } +} + +for (int j = 0; j < bestfit.size(); j++) {//cout << querySeqs[j]->getName() << '\t' << "length = " << querySeqs[j]->getAligned().length() << '\t' << bestfit[j].getName() << " length = " << bestfit[j].getAligned().length() << endl; + //chops off beginning and end of sequences so they both start and end with a base + trimSeqs(querySeqs[j], bestfit[j], j); +//cout << "NEW SEQ PAIR" << querySeqs[j]->getAligned() << endl << "IN THE MIDDLE" << endl << bestfit[j].getAligned() << endl; + +} + + mothurOut("Done."); mothurOutEndLine(); + + windows = findWindows(lines[0]->start, lines[0]->end); + } else { createProcessesSpots(); } + //find P - mothurOut("Calculating expected percentage differences for each sequence... "); cout.flush(); - vector probabilityProfile = calcFreq(templateSeqs); + if (consfile == "") { probabilityProfile = calcFreq(templateSeqs); } + else { probabilityProfile = readFreq(); } //make P into Q for (int i = 0; i < probabilityProfile.size(); i++) { probabilityProfile[i] = 1 - probabilityProfile[i]; } - - //find Qav - averageProbability = findQav(probabilityProfile); - - //find Coefficient - maps a sequence to its coefficient - seqCoef = getCoef(averageProbability); - - //find Eqs for each sequence - the expected distance in each window - Parallelized - if (processors == 1) { calcExpected(lines[0]->start, lines[0]->end); } - else { createProcessesExpected(); } - mothurOut("Done."); mothurOutEndLine(); - //find deviation - Parallelized - mothurOut("Finding deviation from expected... "); cout.flush(); - if (processors == 1) { calcDE(lines[0]->start, lines[0]->end); } - else { createProcessesDE(); } - mothurOut("Done."); mothurOutEndLine(); + if (processors == 1) { + + mothurOut("Calculating observed distance... "); cout.flush(); + obsDistance = calcObserved(lines[0]->start, lines[0]->end); + mothurOut("Done."); mothurOutEndLine(); + + mothurOut("Finding variability... "); cout.flush(); + Qav = findQav(lines[0]->start, lines[0]->end); +for (int i = 0; i < Qav.size(); i++) { +cout << querySeqs[i]->getName() << " = "; +for (int u = 0; u < Qav[i].size();u++) { cout << Qav[i][u] << '\t'; } +cout << endl << endl; +} + + + mothurOut("Done."); mothurOutEndLine(); + + mothurOut("Calculating alpha... "); cout.flush(); + seqCoef = getCoef(lines[0]->start, lines[0]->end); +for (int i = 0; i < seqCoef.size(); i++) { +cout << querySeqs[i]->getName() << " coef = " << seqCoef[i] << endl; +} + + mothurOut("Done."); mothurOutEndLine(); + + mothurOut("Calculating expected distance... "); cout.flush(); + expectedDistance = calcExpected(lines[0]->start, lines[0]->end); + mothurOut("Done."); mothurOutEndLine(); + + mothurOut("Finding deviation... "); cout.flush(); + DE = calcDE(lines[0]->start, lines[0]->end); + deviation = calcDist(lines[0]->start, lines[0]->end); + mothurOut("Done."); mothurOutEndLine(); + + + + } + else { createProcesses(); } - + delete distcalculator; + //free memory for (int i = 0; i < lines.size(); i++) { delete lines[i]; } - delete distCalculator; + } catch(exception& e) { @@ -181,11 +205,83 @@ vector Pintail::readSeqs(string file) { } } +//*************************************************************************************************************** +//num is query's spot in querySeqs +void Pintail::trimSeqs(Sequence* query, Sequence& subject, int num) { + try { + + string q = query->getAligned(); + string s = subject.getAligned(); + + int front = 0; + for (int i = 0; i < q.length(); i++) { + if (isalpha(q[i]) && isalpha(s[i])) { front = i; break; } + } + + q = q.substr(front, q.length()); + s = s.substr(front, s.length()); + + int back = 0; + for (int i = q.length(); i >= 0; i--) { + if (isalpha(q[i]) && isalpha(s[i])) { back = i; break; } + } + + q = q.substr(0, back); + s = s.substr(0, back); + + trim[num][front] = back; + + //save + query->setAligned(q); + query->setUnaligned(q); + subject.setAligned(s); + subject.setUnaligned(s); + } + catch(exception& e) { + errorOut(e, "Pintail", "trimSeqs"); + exit(1); + } +} + +//*************************************************************************************************************** + +vector Pintail::readFreq() { + try { + + ifstream in; + openInputFile(consfile, in); + + vector prob; + + //read in probabilities and store in vector + int pos; float num; + + while(!in.eof()){ + + in >> pos >> num; + + prob.push_back(num); + + gobble(in); + } + + in.close(); + return prob; + + } + catch(exception& e) { + errorOut(e, "Pintail", "readFreq"); + exit(1); + } +} + //*************************************************************************************************************** //calculate the distances from each query sequence to all sequences in the template to find the closest sequence -void Pintail::findPairs(int start, int end) { +vector Pintail::findPairs(int start, int end) { try { + vector seqsMatches; seqsMatches.resize(end-start); + for(int i = start; i < end; i++){ float smallest = 10000.0; @@ -195,94 +291,206 @@ void Pintail::findPairs(int start, int end) { Sequence temp = *(templateSeqs[j]); - distCalculator->calcDist(query, temp); - float dist = distCalculator->getDist(); + distcalculator->calcDist(query, temp); + float dist = distcalculator->getDist(); if (dist < smallest) { - - bestfit[querySeqs[i]] = templateSeqs[j]; + seqsMatches[i] = *(templateSeqs[j]); smallest = dist; } } } + return seqsMatches; + } catch(exception& e) { errorOut(e, "Pintail", "findPairs"); exit(1); } } + //*************************************************************************************************************** -void Pintail::calcObserved(int start, int end) { +//find the window breaks for each sequence +vector< vector > Pintail::findWindows(int start, int end) { try { + + vector< vector > win; win.resize(end-start); + + //for each sequence + int count = 0; + for(int i = start; i < end; i++){ + + //if window is set to default + if (windowSizes[i] == 0) { if (querySeqs[i]->getAligned().length() > 1200) { windowSizes[i] = 300; } + else{ windowSizes[i] = (querySeqs[i]->getAligned().length() / 4); } } + else if (windowSizes[i] > (querySeqs[i]->getAligned().length() / 4)) { + mothurOut("You have selected to large a window size for sequence " + querySeqs[i]->getName() + ". I will choose an appropriate window size."); mothurOutEndLine(); + windowSizes[i] = (querySeqs[i]->getAligned().length() / 4); + } - + //cout << "length = " << querySeqs[i]->getAligned().length() << " window = " << windowSizes[i] << " increment = " << increment << endl; + + + string seq = querySeqs[i]->getAligned(); + int numBases = querySeqs[i]->getUnaligned().length(); + int spot = 0; + + //find location of first base + for (int j = 0; j < seq.length(); j++) { + if (isalpha(seq[j])) { spot = j; break; } + } + + //save start of seq + win[count].push_back(spot); + + + //move ahead increment bases at a time until all bases are in a window + int countBases = 0; + int totalBases = 0; //used to eliminate window of blanks at end of sequence + for (int m = spot; m < seq.length(); m++) { + + //count number of bases you see + if (isalpha(seq[m])) { countBases++; totalBases++; } + + //if you have seen enough bases to make a new window + if (countBases >= increment) { + win[count].push_back(m); //save spot in alignment + countBases = 0; //reset bases you've seen in this window + } + + //no need to continue if all your bases are in a window + if (totalBases == numBases) { break; } + } + + count++; + } + + + + return win; + + } + catch(exception& e) { + errorOut(e, "Pintail", "findWindows"); + exit(1); + } +} + +//*************************************************************************************************************** +vector< vector > Pintail::calcObserved(int start, int end) { + try { + + vector< vector > temp; + temp.resize(end-start); + + int count = 0; for(int i = start; i < end; i++){ - itBest = bestfit.find(querySeqs[i]); - Sequence* query; - Sequence* subject; + Sequence* query = querySeqs[i]; + Sequence subject = bestfit[i]; - if (itBest != bestfit.end()) { - query = itBest->first; - subject = itBest->second; - }else{ mothurOut("Error in calcObserved"); mothurOutEndLine(); } -//cout << query->getName() << '\t' << subject->getName() << endl; - int startpoint = 0; - for (int m = 0; m < iters; m++) { + for (int m = 0; m < windows[i].size(); m++) { - string seqFrag = query->getAligned().substr(startpoint, window); - string seqFragsub = subject->getAligned().substr(startpoint, window); + string seqFrag = query->getAligned().substr(windows[i][startpoint], windowSizes[i]); + string seqFragsub = subject.getAligned().substr(windows[i][startpoint], windowSizes[i]); int diff = 0; for (int b = 0; b < seqFrag.length(); b++) { - //if this is not a gap - if ((isalpha(seqFrag[b])) && (isalpha(seqFragsub[b]))) { + //if either the query or subject is not a gap + if ((isalpha(seqFrag[b])) || (isalpha(seqFragsub[b]))) { //and they are different - penalize if (seqFrag[b] != seqFragsub[b]) { diff++; } } } //percentage of mismatched bases - float dist = diff / (float)seqFrag.length(); + float dist; + dist = diff / (float) seqFrag.length() * 100; - obsDistance[query].push_back(dist); + temp[count].push_back(dist); - startpoint += increment; + startpoint++; } + + count++; } + return temp; } catch(exception& e) { errorOut(e, "Pintail", "calcObserved"); exit(1); } } +//*************************************************************************************************************** +vector Pintail::calcDist(int start, int end) { + try { + + vector temp; + + for(int i = start; i < end; i++){ + + Sequence* query = querySeqs[i]; + Sequence subject = bestfit[i]; + + string seqFrag = query->getAligned(); + string seqFragsub = subject.getAligned(); + + int diff = 0; + for (int b = 0; b < seqFrag.length(); b++) { + + //if either the query or subject is not a gap + if ((isalpha(seqFrag[b])) || (isalpha(seqFragsub[b]))) { + //and they are different - penalize + if (seqFrag[b] != seqFragsub[b]) { diff++; } + } + } + + //percentage of mismatched bases + float dist; + dist = diff / (float) seqFrag.length() * 100; + + temp.push_back(dist); + } + + return temp; + } + catch(exception& e) { + errorOut(e, "Pintail", "calcDist"); + exit(1); + } +} //*************************************************************************************************************** -void Pintail::calcExpected(int start, int end) { +vector< vector > Pintail::calcExpected(int start, int end) { try { - + vector< vector > temp; temp.resize(end-start); + //for each sequence + int count = 0; for(int i = start; i < end; i++){ - itCoef = seqCoef.find(querySeqs[i]); - float coef = itCoef->second; + float coef = seqCoef[i]; //for each window vector queryExpected; - for (int m = 0; m < iters; m++) { - float expected = averageProbability[m] * coef; + for (int m = 0; m < windows[i].size(); m++) { + float expected = Qav[i][m] * coef; queryExpected.push_back(expected); //cout << "average variabilty over window = " << averageProbability[m] << " coef = " << coef << " ei = " << expected << '\t' << "window = " << m << endl; } - expectedDistance[querySeqs[i]] = queryExpected; + temp[count] = queryExpected; + + count++; } + + return temp; } catch(exception& e) { @@ -291,28 +499,30 @@ void Pintail::calcExpected(int start, int end) { } } //*************************************************************************************************************** -void Pintail::calcDE(int start, int end) { +vector Pintail::calcDE(int start, int end) { try { + vector temp; temp.resize(end-start); //for each sequence + int count = 0; for(int i = start; i < end; i++){ - itObsDist = obsDistance.find(querySeqs[i]); - vector obs = itObsDist->second; + vector obs = obsDistance[i]; + vector exp = expectedDistance[i]; - itExpDist = expectedDistance.find(querySeqs[i]); - vector exp = itExpDist->second; // cout << "difference between obs and exp = " << abs(obs[m] - exp[m]) << endl; //for each window float sum = 0.0; //sum = sum from 1 to m of (oi-ei)^2 - for (int m = 0; m < iters; m++) { sum += ((obs[m] - exp[m]) * (obs[m] - exp[m])); } + for (int m = 0; m < windows[i].size(); m++) { sum += ((obs[m] - exp[m]) * (obs[m] - exp[m])); } - float de = sqrt((sum / (iters - 1))); + float de = sqrt((sum / (windows[i].size() - 1))); - DE[querySeqs[i]] = de; + temp[count] = de; + count++; } - + + return temp; } catch(exception& e) { errorOut(e, "Pintail", "calcDE"); @@ -326,6 +536,10 @@ vector Pintail::calcFreq(vector seqs) { try { vector prob; + string freqfile = getRootName(templateFile) + "probability"; + ofstream outFreq; + + openOutputFile(freqfile, outFreq); //at each position in the sequence for (int i = 0; i < seqs[0]->getAligned().length(); i++) { @@ -347,19 +561,25 @@ vector Pintail::calcFreq(vector seqs) { //find base with highest frequency int highest = 0; - for (int m = 0; m < freq.size(); m++) { if (freq[m] > highest) { highest = freq[m]; } } - - //add in gaps - so you can effectively "ignore them" - highest += gaps; + for (int m = 0; m < freq.size(); m++) { if (freq[m] > highest) { highest = freq[m]; } } - float highFreq = highest / (float) seqs.size(); + float highFreq; + //if ( (seqs.size() - gaps) == 0 ) { highFreq = 1.0; } + //else { highFreq = highest / (float) (seqs.size() - gaps); } + highFreq = highest / (float) seqs.size(); +cout << i << '\t' << highFreq << endl; float Pi; - Pi = (highFreq - 0.25) / 0.75; + Pi = (highFreq - 0.25) / 0.75; + + //saves this for later + outFreq << i << '\t' << Pi << endl; prob.push_back(Pi); } + outFreq.close(); + return prob; } @@ -369,23 +589,32 @@ vector Pintail::calcFreq(vector seqs) { } } //*************************************************************************************************************** -vector Pintail::findQav(vector prob) { +vector< vector > Pintail::findQav(int start, int end) { try { - vector averages; + vector< vector > averages; + map::iterator it; - //for each window find average - int startpoint = 0; - for (int m = 0; m < iters; m++) { - - float average = 0.0; - for (int i = startpoint; i < (startpoint+window); i++) { average += prob[i]; } - - average = average / window; -//cout << average << endl; - //save this windows average - averages.push_back(average); + for(int i = start; i < end; i++){ - startpoint += increment; + //for each window find average + vector temp; + for (int m = 0; m < windows[i].size(); m++) { + + float average = 0.0; + + it = trim[i].begin(); //trim[i] is a map of where this sequence was trimmed + + //while you are in the window for this sequence + for (int j = windows[i][m]+it->first; j < (windows[i][m]+windowSizes[i]); j++) { average += probabilityProfile[j]; } + + average = average / windowSizes[i]; + //cout << average << endl; + //save this windows average + temp.push_back(average); + } + + //save this qav + averages.push_back(temp); } return averages; @@ -396,31 +625,34 @@ vector Pintail::findQav(vector prob) { } } //*************************************************************************************************************** -map Pintail::getCoef(vector prob) { +vector Pintail::getCoef(int start, int end) { try { - map coefs; + vector coefs; + coefs.resize(end-start); - //find average prob - float probAverage = 0.0; - for (int i = 0; i < prob.size(); i++) { probAverage += prob[i]; } - probAverage = probAverage / (float) prob.size(); -cout << "(sum of ai) / m = " << probAverage << endl; //find a coef for each sequence - map >::iterator it; - for (it = obsDistance.begin(); it != obsDistance.end(); it++) { - - vector temp = it->second; - Sequence* tempSeq = it->first; + int count = 0; + for(int i = start; i < end; i++){ + + //find average prob for this seqs windows + float probAverage = 0.0; + for (int j = 0; j < Qav[i].size(); j++) { probAverage += Qav[i][j]; } + probAverage = probAverage / (float) Qav[i].size(); + cout << "(sum of ai) / m = " << probAverage << endl; + + vector temp = obsDistance[i]; //find observed average float obsAverage = 0.0; - for (int i = 0; i < temp.size(); i++) { obsAverage += temp[i]; } + for (int j = 0; j < temp.size(); j++) { obsAverage += temp[j]; } obsAverage = obsAverage / (float) temp.size(); -cout << tempSeq->getName() << '\t' << obsAverage << endl; +cout << "(sum of oi) / m = " << obsAverage << endl; float coef = obsAverage / probAverage; -cout << tempSeq->getName() << '\t' << "coef = " << coef << endl; + //save this sequences coefficient - coefs[tempSeq] = coef; + coefs[count] = coef; + + count++; } @@ -432,13 +664,16 @@ cout << tempSeq->getName() << '\t' << "coef = " << coef << endl; } } + /**************************************************************************************************/ -void Pintail::createProcessesPairs() { +void Pintail::createProcessesSpots() { try { #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) int process = 0; vector processIDS; + vector< vector > win; win.resize(querySeqs.size()); + vector< map > t; t.resize(querySeqs.size()); //loop through and create all the processes you want while (process != processors) { @@ -448,7 +683,31 @@ void Pintail::createProcessesPairs() { processIDS.push_back(pid); process++; }else if (pid == 0){ - findPairs(lines[process]->start, lines[process]->end); + + vector tempbest; + tempbest = findPairs(lines[process]->start, lines[process]->end); + int count = 0; + for (int i = lines[process]->start; i < lines[process]->end; i++) { + bestfit[i] = tempbest[count]; + + //chops off beginning and end of sequences so they both start and end with a base + trimSeqs(querySeqs[i], bestfit[i], i); + t[i] = trim[i]; + + count++; + } + + + + vector< vector > temp = findWindows(lines[process]->start, lines[process]->end); + + //move into best + count = 0; + for (int i = lines[process]->start; i < lines[process]->end; i++) { + win[i] = temp[count]; + count++; + } + exit(0); }else { mothurOut("unable to spawn the necessary processes."); mothurOutEndLine(); exit(0); } } @@ -459,62 +718,32 @@ void Pintail::createProcessesPairs() { wait(&temp); } + windows = win; + trim = t; #else - findPairs(lines[0]->start, lines[0]->end); + windows = findWindows(lines[0]->start, lines[0]->end); #endif } catch(exception& e) { - errorOut(e, "Pintail", "createProcessesPairs"); + errorOut(e, "Pintail", "createProcessesSpots"); exit(1); } } + /**************************************************************************************************/ -void Pintail::createProcessesObserved() { +void Pintail::createProcesses() { try { #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) int process = 0; vector 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){ - calcObserved(lines[process]->start, lines[process]->end); - exit(0); - }else { mothurOut("unable to spawn the necessary processes."); mothurOutEndLine(); exit(0); } - } + vector< vector > exp; exp.resize(querySeqs.size()); + vector de; de.resize(querySeqs.size()); + vector< vector > obs; obs.resize(querySeqs.size()); - //force parent to wait until all the processes are done - for (int i=0;istart, lines[0]->end); - -#endif - } - catch(exception& e) { - errorOut(e, "Pintail", "createProcessesObserved"); - exit(1); - } -} - -//*************************************************************************************************************** - -void Pintail::createProcessesExpected() { - try { -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) - int process = 0; - vector processIDS; //loop through and create all the processes you want while (process != processors) { @@ -524,45 +753,48 @@ void Pintail::createProcessesExpected() { processIDS.push_back(pid); process++; }else if (pid == 0){ - calcExpected(lines[process]->start, lines[process]->end); - 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;istart, lines[0]->end); + + vector< vector > temp; + vector tempde; + int count = 0; + + + temp = calcObserved(lines[process]->start, lines[process]->end); + count = 0; + for (int i = lines[process]->start; i < lines[process]->end; i++) { + obs[i] = temp[count]; + count++; + } -#endif - } - catch(exception& e) { - errorOut(e, "Pintail", "createProcessesExpected"); - exit(1); - } -} + temp = findQav(lines[process]->start, lines[process]->end); + count = 0; + for (int i = lines[process]->start; i < lines[process]->end; i++) { + Qav[i] = temp[count]; + count++; + } + + tempde = getCoef(lines[process]->start, lines[process]->end); + count = 0; + for (int i = lines[process]->start; i < lines[process]->end; i++) { + seqCoef[i] = tempde[count]; + count++; + } + + temp = calcExpected(lines[process]->start, lines[process]->end); + count = 0; + for (int i = lines[process]->start; i < lines[process]->end; i++) { + exp[i] = temp[count]; + count++; + } -/**************************************************************************************************/ + + tempde = calcDE(lines[process]->start, lines[process]->end); + count = 0; + for (int i = lines[process]->start; i < lines[process]->end; i++) { + de[i] = tempde[count]; + count++; + } -void Pintail::createProcessesDE() { - try { -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) - int process = 0; - vector 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){ - calcDE(lines[process]->start, lines[process]->end); exit(0); }else { mothurOut("unable to spawn the necessary processes."); mothurOutEndLine(); exit(0); } } @@ -573,13 +805,22 @@ void Pintail::createProcessesDE() { wait(&temp); } + obsDistance = obs; + expectedDistance = exp; + DE = de; + #else - calcDE(lines[0]->start, lines[0]->end); + bestfit = findPairs(lines[0]->start, lines[0]->end); + obsDistance = calcObserved(lines[0]->start, lines[0]->end); + Qav = findQav(lines[0]->start, lines[0]->end); + seqCoef = getCoef(lines[0]->start, lines[0]->end); + expectedDistance = calcExpected(lines[0]->start, lines[0]->end); + DE = calcDE(lines[0]->start, lines[0]->end); #endif } catch(exception& e) { - errorOut(e, "Pintail", "createProcessesDE"); + errorOut(e, "Pintail", "createProcesses"); exit(1); } } diff --git a/pintail.h b/pintail.h index 7c0f81b..797c99e 100644 --- a/pintail.h +++ b/pintail.h @@ -22,12 +22,14 @@ class Pintail : public Chimera { public: - Pintail(string); + Pintail(string, string); ~Pintail(); void getChimeras(); void print(ostream&); + void setCons(string c) { consfile = c; } + private: @@ -37,41 +39,47 @@ class Pintail : public Chimera { linePair(int i, int j) : start(i), end(j) {} }; - - Dist* distCalculator; - string fastafile; + Dist* distcalculator; int iters; + string fastafile, templateFile, consfile; vector lines; vector querySeqs; vector templateSeqs; - map bestfit; //maps a query sequence to its most similiar sequence in the template - map::iterator itBest; + vector bestfit; //bestfit[0] matches queryseqs[0]... - map > obsDistance; //maps a query sequence to its observed distance at each window - map > expectedDistance; //maps a query sequence to its expected distance at each window - map >::iterator itObsDist; - map >::iterator itExpDist; + vector< vector > obsDistance; //obsDistance[0] is the vector of observed distances for queryseqs[0]... + vector< vector > expectedDistance; //expectedDistance[0] is the vector of expected distances for queryseqs[0]... + vector deviation; //deviation[0] is the percentage of mismatched pairs over the whole seq between querySeqs[0] and its best match. + vector< vector > windows; // windows[0] is a vector containing the starting spot in queryseqs[0] aligned sequence for each window. + //this is needed so you can move by bases and not just spots in the alignment + vector< map > trim; //trim[0] is the start and end position of trimmed querySeqs[0]. Used to find the variability over each sequence window. + + vector windowSizes; //windowSizes[0] = window size of querySeqs[0] - vector averageProbability; //Qav - map seqCoef; //maps a sequence to its coefficient - map DE; //maps a sequence to its deviation - map::iterator itCoef; + vector< vector > Qav; //Qav[0] is the vector of average variablility for queryseqs[0]... + vector seqCoef; //seqCoef[0] is the coeff for queryseqs[0]... + vector DE; //DE[0] is the deviaation for queryseqs[0]... + vector probabilityProfile; vector readSeqs(string); - vector findQav(vector); + void trimSeqs(Sequence*, Sequence&, int); + vector readFreq(); + vector< vector > findQav(int, int); vector calcFreq(vector); - map getCoef(vector); + vector getCoef(int, int); - void findPairs(int, int); - void calcObserved(int, int); - void calcExpected(int, int); - void calcDE(int, int); + vector findPairs(int, int); + vector< vector > findWindows(int, int); + vector< vector > calcObserved(int, int); + vector< vector > calcExpected(int, int); + vector calcDE(int, int); + vector calcDist(int, int); - void createProcessesPairs(); - void createProcessesObserved(); - void createProcessesExpected(); - void createProcessesDE(); + void createProcessesSpots(); + void createProcesses(); + + }; diff --git a/sharedchao1.cpp b/sharedchao1.cpp index 402e0d9..803e953 100644 --- a/sharedchao1.cpp +++ b/sharedchao1.cpp @@ -29,7 +29,7 @@ EstOutput SharedChao1::getValues(vector shared){ //create and initialize trees to 0. initialTree(numGroups); - //loop through vectors calculating the f11, f1A, f2A, f1B, f2B, S12 values + for (int i = 0; i < shared[0]->size(); i++) { //get bin values and calc shared bool sharedByAll = true; @@ -50,7 +50,7 @@ EstOutput SharedChao1::getValues(vector shared){ //calculate chao1, (numleaves-1) because numleaves contains the ++ values. bool bias; for(int i=0;ilvalue == 0) { bias = true;}// break;} + if ((f2leaves[i]->lvalue == 0) || (f2leaves[i]->rvalue == 0)) { bias = true; }// break;} } if(bias){ @@ -60,6 +60,7 @@ EstOutput SharedChao1::getValues(vector shared){ if (i != (numLeaves-1)) { rightvalue = (float)(f1leaves[i]->rvalue * (f1leaves[i]->rvalue - 1)) / (float)((pow(2, (float)f2leaves[i]->rcoef)) * (f2leaves[i]->rvalue + 1)); }else{ + //add in sobs rightvalue = (float)(f1leaves[i]->rvalue); } Chao += leftvalue + rightvalue; @@ -73,6 +74,7 @@ EstOutput SharedChao1::getValues(vector shared){ if (i != (numLeaves-1)) { rightvalue = (float)(f1leaves[i]->rvalue * f1leaves[i]->rvalue) / (float)((pow(2, (float)f2leaves[i]->rcoef)) * f2leaves[i]->rvalue); }else{ + //add in sobs rightvalue = (float)(f1leaves[i]->rvalue); } Chao += leftvalue + rightvalue; @@ -84,7 +86,7 @@ EstOutput SharedChao1::getValues(vector shared){ delete f2leaves[i]; } - + data[0] = Chao; return data; } diff --git a/venn.cpp b/venn.cpp index 217dd35..7cb2134 100644 --- a/venn.cpp +++ b/venn.cpp @@ -133,9 +133,9 @@ void Venn::getPic(vector lookup, vector vCalcs singleCalc = new Sobs(); }else if (vCalcs[i]->getName() == "sharedchao") { singleCalc = new Chao1(); - }else if (vCalcs[i]->getName() == "sharedace") { - singleCalc = new Ace(10); - } + }//else if (vCalcs[i]->getName() == "sharedace") { + //singleCalc = new Ace(10); + //} //get estimates for numA vector numA = singleCalc->getValues(sabundA); @@ -190,154 +190,234 @@ void Venn::getPic(vector lookup, vector vCalcs //make a file for each calculator for(int i=0;igetName() == "sharedsobs") { - singleCalc = new Sobs(); - }else if (vCalcs[i]->getName() == "sharedchao") { - singleCalc = new Chao1(); - }else if (vCalcs[i]->getName() == "sharedace") { - singleCalc = new Ace(10); - } + + string filenamesvg = getRootName(globaldata->inputFileName) + lookup[0]->getLabel() + ".venn." + vCalcs[i]->getName() + ".svg"; + openOutputFile(filenamesvg, outsvg); - //get estimates for numA - vector numA = singleCalc->getValues(sabundA); + if (vCalcs[i]->getName() == "sharedace") { + + singleCalc = new Ace(10); + + //get estimates for numA + vector numA = singleCalc->getValues(sabundA); - //get estimates for numB - vector numB = singleCalc->getValues(sabundB); + //get estimates for numB + vector numB = singleCalc->getValues(sabundB); - //get estimates for numC - vector numC = singleCalc->getValues(sabundC); + //get estimates for numC + vector numC = singleCalc->getValues(sabundC); - - string filenamesvg = getRootName(globaldata->inputFileName) + lookup[0]->getLabel() + ".venn." + vCalcs[i]->getName() + ".svg"; - openOutputFile(filenamesvg, outsvg); - //get estimates for sharedAB, sharedAC and sharedBC - subset.clear(); - subset.push_back(lookup[0]); subset.push_back(lookup[1]); - vector sharedAB = vCalcs[i]->getValues(subset); - - subset.clear(); - subset.push_back(lookup[0]); subset.push_back(lookup[2]); - vector sharedAC = vCalcs[i]->getValues(subset); - - subset.clear(); - subset.push_back(lookup[1]); subset.push_back(lookup[2]); - vector sharedBC = vCalcs[i]->getValues(subset); - - vector sharedAwithBC; - vector sharedBwithAC; - vector sharedCwithAB; - - //find possible sharedABC values - float sharedABC1 = 0.0; float sharedABC2 = 0.0; float sharedABC3 = 0.0; float sharedABC = 0.0; + //get estimates for sharedAB, sharedAC and sharedBC + subset.clear(); + subset.push_back(lookup[0]); subset.push_back(lookup[1]); + vector sharedAB = vCalcs[i]->getValues(subset); + + subset.clear(); + subset.push_back(lookup[0]); subset.push_back(lookup[2]); + vector sharedAC = vCalcs[i]->getValues(subset); + + subset.clear(); + subset.push_back(lookup[1]); subset.push_back(lookup[2]); + vector sharedBC = vCalcs[i]->getValues(subset); + + vector sharedAwithBC; + vector sharedBwithAC; + vector sharedCwithAB; + + //find possible sharedABC values + float sharedABC1 = 0.0; float sharedABC2 = 0.0; float sharedABC3 = 0.0; float sharedABC = 0.0; - if (vCalcs[i]->getMultiple() == false) { - //merge BC and estimate with shared with A - SharedRAbundVector* merge = new SharedRAbundVector(); - for (int j = 0; j < lookup[1]->size(); j++) { - merge->push_back((lookup[1]->getAbundance(j) + lookup[2]->getAbundance(j)), j, ""); - } + if (vCalcs[i]->getMultiple() == false) { + //merge BC and estimate with shared with A + SharedRAbundVector* merge = new SharedRAbundVector(); + for (int j = 0; j < lookup[1]->size(); j++) { + merge->push_back((lookup[1]->getAbundance(j) + lookup[2]->getAbundance(j)), j, ""); + } + + subset.clear(); + subset.push_back(lookup[0]); subset.push_back(merge); + sharedAwithBC = vCalcs[i]->getValues(subset); - subset.clear(); - subset.push_back(lookup[0]); subset.push_back(merge); - sharedAwithBC = vCalcs[i]->getValues(subset); - - delete merge; - //merge AC and estimate with shared with B - merge = new SharedRAbundVector(); - for (int j = 0; j < lookup[0]->size(); j++) { - merge->push_back((lookup[0]->getAbundance(j) + lookup[2]->getAbundance(j)), j, ""); - } + delete merge; + //merge AC and estimate with shared with B + merge = new SharedRAbundVector(); + for (int j = 0; j < lookup[0]->size(); j++) { + merge->push_back((lookup[0]->getAbundance(j) + lookup[2]->getAbundance(j)), j, ""); + } + + subset.clear(); + subset.push_back(merge); subset.push_back(lookup[1]); + sharedBwithAC = vCalcs[i]->getValues(subset); - subset.clear(); - subset.push_back(merge); subset.push_back(lookup[1]); - sharedBwithAC = vCalcs[i]->getValues(subset); + delete merge; + //merge AB and estimate with shared with C + merge = new SharedRAbundVector(); + for (int j = 0; j < lookup[0]->size(); j++) { + merge->push_back((lookup[0]->getAbundance(j) + lookup[1]->getAbundance(j)), j, ""); + } + + subset.clear(); + subset.push_back(lookup[2]); subset.push_back(merge); + sharedCwithAB = vCalcs[i]->getValues(subset); + delete merge; + + sharedABC1 = sharedAB[0] + sharedAC[0] - sharedAwithBC[0]; + sharedABC2 = sharedAB[0] + sharedBC[0] - sharedBwithAC[0]; + sharedABC3 = sharedAC[0] + sharedBC[0] - sharedCwithAB[0]; + + //if any of the possible m's are - throw them out + if (sharedABC1 < 0.00001) { sharedABC1 = 0; } + if (sharedABC2 < 0.00001) { sharedABC2 = 0; } + if (sharedABC3 < 0.00001) { sharedABC3 = 0; } - delete merge; - //merge AB and estimate with shared with C - merge = new SharedRAbundVector(); - for (int j = 0; j < lookup[0]->size(); j++) { - merge->push_back((lookup[0]->getAbundance(j) + lookup[1]->getAbundance(j)), j, ""); + //sharedABC is the minimum of the 3 possibilities + if ((sharedABC1 < sharedABC2) && (sharedABC1 < sharedABC3)) { sharedABC = sharedABC1; } + else if ((sharedABC2 < sharedABC1) && (sharedABC2 < sharedABC3)) { sharedABC = sharedABC2; } + else if ((sharedABC3 < sharedABC1) && (sharedABC3 < sharedABC2)) { sharedABC = sharedABC3; } + }else{ + vector data = vCalcs[i]->getValues(lookup); + sharedABC = data[0]; + sharedAwithBC.push_back(sharedAB[0] + sharedAC[0] - sharedABC); + sharedBwithAC.push_back(sharedAB[0] + sharedBC[0] - sharedABC); + sharedCwithAB.push_back(sharedAC[0] + sharedBC[0] - sharedABC); } + + //image window + outsvg << "\n"; + outsvg << "\n"; + + //draw circles + outsvg << ""; + outsvg << "Venn Diagram at distance " + lookup[0]->getLabel() + "\n"; + outsvg << ""; + outsvg << ""; + outsvg << ""; + + //place labels within overlaps + outsvg << "" + toString(numA[0]-sharedAwithBC[0]) + "\n"; + outsvg << "getGroup().length() / 2)) + "\" y=\"150\">" + lookup[0]->getGroup() + "\n"; + outsvg << "" + toString(sharedAB[0] - sharedABC) + "\n"; + outsvg << "" + toString(numB[0]-sharedBwithAC[0]) + "\n"; + outsvg << "getGroup().length() / 2)) + "\" y=\"150\">" + lookup[1]->getGroup() + "\n"; + outsvg << "" + toString(sharedAC[0] - sharedABC) + "\n"; + outsvg << "" + toString(numC[0]-sharedCwithAB[0]) + "\n"; + outsvg << "getGroup().length() / 2)) + "\" y=\"410\">" + lookup[2]->getGroup() + "\n"; + outsvg << "" + toString(sharedBC[0] - sharedABC) + "\n"; + outsvg << "" + toString(sharedABC) + "\n"; + outsvg << "The number of sepecies shared between groups " + lookup[0]->getGroup() + " and " + lookup[1]->getGroup() + " is " + toString(sharedAB[0]) + "\n"; + outsvg << "The number of sepecies shared between groups " + lookup[0]->getGroup() + " and " + lookup[2]->getGroup() + " is " + toString(sharedAC[0]) + "\n"; + outsvg << "The number of sepecies shared between groups " + lookup[1]->getGroup() + " and " + lookup[2]->getGroup() + " is " + toString(sharedBC[0]) + "\n"; + outsvg << "The number of sepecies shared between groups " + lookup[0]->getGroup() + " and combined groups " + lookup[1]->getGroup() + lookup[2]->getGroup() + " is " + toString(sharedAwithBC[0]) + "\n"; + outsvg << "The number of sepecies shared between groups " + lookup[1]->getGroup() + " and combined groups " + lookup[0]->getGroup() + lookup[2]->getGroup() + " is " + toString(sharedBwithAC[0]) + "\n"; + outsvg << "The number of sepecies shared between groups " + lookup[2]->getGroup() + " and combined groups " + lookup[0]->getGroup() + lookup[1]->getGroup() + " is " + toString(sharedCwithAB[0]) + "\n"; + outsvg << "The number of species in group " + lookup[0]->getGroup() + " is " + toString(numA[0]); + if (numA.size() == 3) { + outsvg << " the lci is " + toString(numA[1]) + " and the hci is " + toString(numA[2]) + "\n"; + }else { outsvg << "\n"; } + + outsvg << "The number of species in group " + lookup[1]->getGroup() + " is " + toString(numB[0]); + if (numB.size() == 3) { + outsvg << " the lci is " + toString(numB[1]) + " and the hci is " + toString(numB[2]) + "\n"; + }else { outsvg << "\n"; } + + outsvg << "The number of species in group " + lookup[2]->getGroup() + " is " + toString(numC[0]); + if (numC.size() == 3) { + outsvg << " the lci is " + toString(numC[1]) + " and the hci is " + toString(numC[2]) + "\n"; + }else { outsvg << "\n"; } + + outsvg << "The total richness of all the groups is " + toString(numA[0] + numB[0] + numC[0] - sharedAB[0] - sharedAC[0] - sharedBC[0] + sharedABC) + "\n"; + outsvg << "The total shared richness is " + toString(sharedABC) + "\n"; + + delete singleCalc; + + }else { //sharedchao and sharedsobs are multigroup + + vector subset; + + //get estimates for numA + subset.push_back(lookup[0]); + vector numA = vCalcs[i]->getValues(subset); + + //get estimates for numB subset.clear(); - subset.push_back(lookup[2]); subset.push_back(merge); - sharedCwithAB = vCalcs[i]->getValues(subset); - delete merge; - - sharedABC1 = sharedAB[0] + sharedAC[0] - sharedAwithBC[0]; - sharedABC2 = sharedAB[0] + sharedBC[0] - sharedBwithAC[0]; - sharedABC3 = sharedAC[0] + sharedBC[0] - sharedCwithAB[0]; - - //if any of the possible m's are - throw them out - if (sharedABC1 < 0.00001) { sharedABC1 = 0; } - if (sharedABC2 < 0.00001) { sharedABC2 = 0; } - if (sharedABC3 < 0.00001) { sharedABC3 = 0; } - - //sharedABC is the minimum of the 3 possibilities - if ((sharedABC1 < sharedABC2) && (sharedABC1 < sharedABC3)) { sharedABC = sharedABC1; } - else if ((sharedABC2 < sharedABC1) && (sharedABC2 < sharedABC3)) { sharedABC = sharedABC2; } - else if ((sharedABC3 < sharedABC1) && (sharedABC3 < sharedABC2)) { sharedABC = sharedABC3; } - }else{ - vector data = vCalcs[i]->getValues(lookup); - sharedABC = data[0]; - sharedAwithBC.push_back(sharedAB[0] + sharedAC[0] - sharedABC); - sharedBwithAC.push_back(sharedAB[0] + sharedBC[0] - sharedABC); - sharedCwithAB.push_back(sharedAC[0] + sharedBC[0] - sharedABC); - } - - //image window - outsvg << "\n"; - outsvg << "\n"; + subset.push_back(lookup[1]); + vector numB = vCalcs[i]->getValues(subset); + + //get estimates for numC + subset.clear(); + subset.push_back(lookup[2]); + vector numC = vCalcs[i]->getValues(subset); - //draw circles - outsvg << ""; - outsvg << "Venn Diagram at distance " + lookup[0]->getLabel() + "\n"; - outsvg << ""; - outsvg << ""; - outsvg << ""; + subset.clear(); + subset.push_back(lookup[0]); subset.push_back(lookup[1]); + vector sharedab = vCalcs[i]->getValues(subset); + + subset.clear(); + subset.push_back(lookup[0]); subset.push_back(lookup[2]); + vector sharedac = vCalcs[i]->getValues(subset); + + subset.clear(); + subset.push_back(lookup[1]); subset.push_back(lookup[2]); + vector sharedbc = vCalcs[i]->getValues(subset); + + subset.clear(); + subset.push_back(lookup[0]); subset.push_back(lookup[1]); subset.push_back(lookup[2]); + vector sharedabc = vCalcs[i]->getValues(subset); + + //image window + outsvg << "\n"; + outsvg << "\n"; - //place labels within overlaps - outsvg << "" + toString(numA[0]-sharedAwithBC[0]) + "\n"; - outsvg << "getGroup().length() / 2)) + "\" y=\"150\">" + lookup[0]->getGroup() + "\n"; - outsvg << "" + toString(sharedAB[0] - sharedABC) + "\n"; - outsvg << "" + toString(numB[0]-sharedBwithAC[0]) + "\n"; - outsvg << "getGroup().length() / 2)) + "\" y=\"150\">" + lookup[1]->getGroup() + "\n"; - outsvg << "" + toString(sharedAC[0] - sharedABC) + "\n"; - outsvg << "" + toString(numC[0]-sharedCwithAB[0]) + "\n"; - outsvg << "getGroup().length() / 2)) + "\" y=\"410\">" + lookup[2]->getGroup() + "\n"; - outsvg << "" + toString(sharedBC[0] - sharedABC) + "\n"; - outsvg << "" + toString(sharedABC) + "\n"; - - outsvg << "The number of sepecies shared between groups " + lookup[0]->getGroup() + " and " + lookup[1]->getGroup() + " is " + toString(sharedAB[0]) + "\n"; - outsvg << "The number of sepecies shared between groups " + lookup[0]->getGroup() + " and " + lookup[2]->getGroup() + " is " + toString(sharedAC[0]) + "\n"; - outsvg << "The number of sepecies shared between groups " + lookup[1]->getGroup() + " and " + lookup[2]->getGroup() + " is " + toString(sharedBC[0]) + "\n"; - outsvg << "The number of sepecies shared between groups " + lookup[0]->getGroup() + " and combined groups " + lookup[1]->getGroup() + lookup[2]->getGroup() + " is " + toString(sharedAwithBC[0]) + "\n"; - outsvg << "The number of sepecies shared between groups " + lookup[1]->getGroup() + " and combined groups " + lookup[0]->getGroup() + lookup[2]->getGroup() + " is " + toString(sharedBwithAC[0]) + "\n"; - outsvg << "The number of sepecies shared between groups " + lookup[2]->getGroup() + " and combined groups " + lookup[0]->getGroup() + lookup[1]->getGroup() + " is " + toString(sharedCwithAB[0]) + "\n"; - outsvg << "The number of species in group " + lookup[0]->getGroup() + " is " + toString(numA[0]); - if (numA.size() == 3) { - outsvg << " the lci is " + toString(numA[1]) + " and the hci is " + toString(numA[2]) + "\n"; - }else { outsvg << "\n"; } - - outsvg << "The number of species in group " + lookup[1]->getGroup() + " is " + toString(numB[0]); - if (numB.size() == 3) { - outsvg << " the lci is " + toString(numB[1]) + " and the hci is " + toString(numB[2]) + "\n"; - }else { outsvg << "\n"; } - - outsvg << "The number of species in group " + lookup[2]->getGroup() + " is " + toString(numC[0]); - if (numC.size() == 3) { - outsvg << " the lci is " + toString(numC[1]) + " and the hci is " + toString(numC[2]) + "\n"; - }else { outsvg << "\n"; } + //draw circles + outsvg << ""; + outsvg << "Venn Diagram at distance " + lookup[0]->getLabel() + "\n"; + outsvg << ""; + outsvg << ""; + outsvg << ""; - outsvg << "The total richness of all the groups is " + toString(numA[0] + numB[0] + numC[0] - sharedAB[0] - sharedAC[0] - sharedBC[0] + sharedABC) + "\n"; - outsvg << "The total shared richness is " + toString(sharedABC) + "\n"; + //place labels within overlaps + outsvg << "" + toString(numA[0]-sharedab[0]-sharedac[0]+sharedabc[0]) + "\n"; + outsvg << "getGroup().length() / 2)) + "\" y=\"150\">" + lookup[0]->getGroup() + "\n"; + outsvg << "" + toString(sharedab[0] - sharedabc[0]) + "\n"; + outsvg << "" + toString(numB[0]-sharedab[0]-sharedbc[0]+sharedabc[0]) + "\n"; + outsvg << "getGroup().length() / 2)) + "\" y=\"150\">" + lookup[1]->getGroup() + "\n"; + outsvg << "" + toString(sharedac[0] - sharedabc[0]) + "\n"; + outsvg << "" + toString(numC[0]-sharedac[0]-sharedbc[0]+sharedabc[0]) + "\n"; + outsvg << "getGroup().length() / 2)) + "\" y=\"410\">" + lookup[2]->getGroup() + "\n"; + outsvg << "" + toString(sharedbc[0] - sharedabc[0]) + "\n"; + outsvg << "" + toString(sharedabc[0]) + "\n"; + outsvg << "The number of sepecies shared between groups " + lookup[0]->getGroup() + " and " + lookup[1]->getGroup() + " is " + toString(sharedab[0]) + "\n"; + outsvg << "The number of sepecies shared between groups " + lookup[0]->getGroup() + " and " + lookup[2]->getGroup() + " is " + toString(sharedac[0]) + "\n"; + outsvg << "The number of sepecies shared between groups " + lookup[1]->getGroup() + " and " + lookup[2]->getGroup() + " is " + toString(sharedbc[0]) + "\n"; + outsvg << "The number of species in group " + lookup[0]->getGroup() + " is " + toString(numA[0]); + if (numA.size() == 3) { + outsvg << " the lci is " + toString(numA[1]) + " and the hci is " + toString(numA[2]) + "\n"; + }else { outsvg << "\n"; } + + outsvg << "The number of species in group " + lookup[1]->getGroup() + " is " + toString(numB[0]); + if (numB.size() == 3) { + outsvg << " the lci is " + toString(numB[1]) + " and the hci is " + toString(numB[2]) + "\n"; + }else { outsvg << "\n"; } + + outsvg << "The number of species in group " + lookup[2]->getGroup() + " is " + toString(numC[0]); + if (numC.size() == 3) { + outsvg << " the lci is " + toString(numC[1]) + " and the hci is " + toString(numC[2]) + "\n"; + }else { outsvg << "\n"; } + + outsvg << "The total richness of all the groups is " + toString(numA[0] + numB[0] + numC[0] - sharedab[0] - sharedac[0] - sharedbc[0] + sharedabc[0]) + "\n"; + outsvg << "The total shared richness is " + toString(sharedabc[0]) + "\n"; + + + } + + //close file outsvg << "\n\n"; outsvg.close(); - delete singleCalc; + } @@ -451,7 +531,33 @@ void Venn::getPic(vector lookup, vector vCalcs //get estimate for all four data = vCalcs[i]->getValues(lookup); sharedABCD = data[0]; - //cout << "num abcd = " << sharedABCD << endl << endl; + //cout << "num abcd = " << sharedABCD << endl << endl; + + + + //image window + outsvg << "\n"; + outsvg << "\n"; + outsvg << ""; + outsvg << "Venn Diagram at distance " + lookup[0]->getLabel() + "\n"; + + outsvg << "The number of species in group " + lookup[0]->getGroup() + " is " + toString(numA) + "\n"; + outsvg << "The number of species in group " + lookup[1]->getGroup() + " is " + toString(numB) + "\n"; + outsvg << "The number of species in group " + lookup[2]->getGroup() + " is " + toString(numC) + "\n"; + outsvg << "The number of species in group " + lookup[3]->getGroup() + " is " + toString(numD) + "\n"; + + outsvg << "The number of sepecies shared between groups " + lookup[0]->getGroup() + " and " + lookup[1]->getGroup() + " is " + toString(sharedAB) + "\n"; + outsvg << "The number of sepecies shared between groups " + lookup[0]->getGroup() + " and " + lookup[2]->getGroup() + " is " + toString(sharedAC) + "\n"; + outsvg << "The number of sepecies shared between groups " + lookup[0]->getGroup() + " and " + lookup[3]->getGroup() + " is " + toString(sharedAD) + "\n"; + outsvg << "The number of sepecies shared between groups " + lookup[1]->getGroup() + " and " + lookup[2]->getGroup() + " is " + toString(sharedBC) + "\n"; + outsvg << "The number of sepecies shared between groups " + lookup[1]->getGroup() + " and " + lookup[3]->getGroup() + " is " + toString(sharedBD) + "\n"; + outsvg << "The number of sepecies shared between groups " + lookup[2]->getGroup() + " and " + lookup[3]->getGroup() + " is " + toString(sharedCD) + "\n"; + + outsvg << "The number of sepecies shared between groups " + lookup[0]->getGroup() + ", " + lookup[1]->getGroup() + " and " + lookup[2]->getGroup() + " is " + toString(sharedABC) + "\n"; + outsvg << "The number of sepecies shared between groups " + lookup[0]->getGroup() + ", " + lookup[1]->getGroup() + " and " + lookup[3]->getGroup() + " is " + toString(sharedABD) + "\n"; + outsvg << "The number of sepecies shared between groups " + lookup[0]->getGroup() + ", " + lookup[2]->getGroup() + " and " + lookup[3]->getGroup() + " is " + toString(sharedACD) + "\n"; + outsvg << "The number of sepecies shared between groups " + lookup[1]->getGroup() + ", " + lookup[2]->getGroup() + " and " + lookup[3]->getGroup() + " is " + toString(sharedBCD) + "\n"; + //make adjustments sharedABC = sharedABC - sharedABCD; //cout << "num abc = " << sharedABC << endl; @@ -462,13 +568,12 @@ void Venn::getPic(vector lookup, vector vCalcs sharedBCD = sharedBCD - sharedABCD; //cout << "num bcd = " << sharedBCD << endl; - sharedAB = sharedAB - sharedABC - sharedABCD - sharedABD; // cout << "num ab = " << sharedAB << endl; - sharedAC = sharedAC - sharedABC - sharedABCD - sharedACD; // cout << "num ac = " << sharedAC << endl; - sharedAD = sharedAD - sharedABD - sharedABCD - sharedACD; // cout << "num ad = " << sharedAD << endl; - - sharedBC = sharedBC - sharedABC - sharedABCD - sharedBCD; // cout << "num bc = " << sharedBC << endl; + sharedAB = sharedAB - sharedABC - sharedABCD - sharedABD; //cout << "num ab = " << sharedAB << endl; + sharedAC = sharedAC - sharedABC - sharedABCD - sharedACD; //cout << "num ac = " << sharedAC << endl; + sharedAD = sharedAD - sharedABD - sharedABCD - sharedACD; //cout << "num ad = " << sharedAD << endl; + sharedBC = sharedBC - sharedABC - sharedABCD - sharedBCD; //cout << "num bc = " << sharedBC << endl; sharedBD = sharedBD - sharedABD - sharedABCD - sharedBCD; // cout << "num bd = " << sharedBD << endl; - sharedCD = sharedCD - sharedACD - sharedABCD - sharedBCD; // cout << "num cd = " << sharedCD << endl; + sharedCD = sharedCD - sharedACD - sharedABCD - sharedBCD; //cout << "num cd = " << sharedCD << endl; numA = numA - sharedAB - sharedAC - sharedAD - sharedABCD - sharedABC - sharedACD - sharedABD; //cout << "num a = " << numA << endl; @@ -479,13 +584,7 @@ void Venn::getPic(vector lookup, vector vCalcs numD = numD - sharedAD - sharedBD - sharedCD - sharedABCD - sharedBCD - sharedACD - sharedABD; //cout << "num d = " << numD << endl; - //image window - outsvg << "\n"; - outsvg << "\n"; - //draw circles - outsvg << ""; - outsvg << "Venn Diagram at distance " + lookup[0]->getLabel() + "\n"; outsvg << "\n "; outsvg << "\n "; outsvg << "\n "; @@ -494,27 +593,31 @@ void Venn::getPic(vector lookup, vector vCalcs //A = red, B = green, C = blue, D = yellow //place labels within overlaps - outsvg << "" + toString(numA) + "\n"; + outsvg << "" + toString(numA) + "\n"; outsvg << "getGroup().length() / 2)) + "\" y=\"90\">" + lookup[0]->getGroup() + "\n"; outsvg << "" + toString(sharedAB) + "\n"; - outsvg << "" + toString(numB) + "\n"; + outsvg << "" + toString(numB) + "\n"; outsvg << "getGroup().length() / 2)) + "\" y=\"90\">" + lookup[1]->getGroup() + "\n"; outsvg << "" + toString(sharedAC) + "\n"; outsvg << "" + toString(numC) + "\n"; outsvg << "getGroup().length() / 2)) + "\" y=\"210\">" + lookup[2]->getGroup() + "\n"; - outsvg << "" + toString(sharedBC) + "\n"; + outsvg << "" + toString(sharedBD) + "\n"; outsvg << "" + toString(numD) + "\n"; outsvg << "getGroup().length() / 2)) + "\" y=\"210\">" + lookup[3]->getGroup() + "\n"; outsvg << "" + toString(sharedAD) + "\n"; - outsvg << "" + toString(sharedBD) + "\n"; + outsvg << "" + toString(sharedBC) + "\n"; outsvg << "" + toString(sharedCD) + "\n"; outsvg << "" + toString(sharedABD) + "\n"; outsvg << "" + toString(sharedBCD) + "\n"; outsvg << "" + toString(sharedACD) + "\n"; outsvg << "" + toString(sharedABC) + "\n"; outsvg << "" + toString(sharedABCD) + "\n"; - outsvg << "The total richness of all the groups is " + toString((float)(numA + numB + numC + numD + sharedAB + sharedAC + sharedAD + sharedBC + sharedBD + sharedCD + sharedABC + sharedABD + sharedACD + sharedBCD + sharedABCD)) + "\n"; - + + + + outsvg << "The total richness of all the groups is " + toString((float)(numA + numB + numC + numD + sharedAB + sharedAC + sharedAD + sharedBC + sharedBD + sharedCD + sharedABC + sharedABD + sharedACD + sharedBCD + sharedABCD)) + "\n"; + + outsvg << "\n\n"; outsvg.close(); delete singleCalc; -- 2.39.2