5 * Created by westcott on 11/4/09.
6 * Copyright 2009 Schloss Lab. All rights reserved.
12 /**************************************************************************************************/
13 Knn::Knn(string tfile, string tempFile, string method, int kmerSize, float gapOpen, float gapExtend, float match, float misMatch, int n, int tid)
14 : Classify(), num(n), search(method) {
19 //create search database and names vector
20 generateDatabaseAndNames(tfile, tempFile, method, kmerSize, gapOpen, gapExtend, match, misMatch);
23 m->errorOut(e, "Knn", "Knn");
27 /**************************************************************************************************/
28 void Knn::setDistName(string s) {
32 m->openOutputFile(outDistName, outDistance);
33 outDistance << "Name\tBestMatch\tDistance" << endl;
37 m->errorOut(e, "Knn", "setDistName");
41 /**************************************************************************************************/
45 if (database != NULL) { delete database; }
48 m->errorOut(e, "Knn", "~Knn");
52 /**************************************************************************************************/
53 string Knn::getTaxonomy(Sequence* seq) {
57 //use database to find closest seq
58 vector<int> closest = database->findClosestSequences(seq, num);
60 if (search == "distance") { ofstream outDistance; m->openOutputFileAppend(outDistName, outDistance); outDistance << seq->getName() << '\t' << database->getName(closest[0]) << '\t' << database->getSearchScore() << endl; outDistance.close(); }
62 if (m->control_pressed) { return tax; }
64 vector<string> closestNames;
65 for (int i = 0; i < closest.size(); i++) {
66 //find that sequences taxonomy in map
67 it = taxonomy.find(names[closest[i]]);
69 //is this sequence in the taxonomy file
70 if (it == taxonomy.end()) { //error not in file
71 m->mothurOut("Error: sequence " + names[closest[i]] + " is not in the taxonomy file. It will be eliminated as a match to sequence " + seq->getName() + "."); m->mothurOutEndLine();
72 }else{ closestNames.push_back(it->first); }
75 if (closestNames.size() == 0) {
76 m->mothurOut("Error: All the matches for sequence " + seq->getName() + " have been eliminated. "); m->mothurOutEndLine();
79 tax = findCommonTaxonomy(closestNames);
80 if (tax == "") { m->mothurOut("There are no common levels for sequence " + seq->getName() + ". "); m->mothurOutEndLine(); tax = "unknown;"; }
87 m->errorOut(e, "Knn", "getTaxonomy");
91 /**************************************************************************************************/
92 string Knn::findCommonTaxonomy(vector<string> closest) {
94 /*vector< vector<string> > taxons; //taxon[0] = vector of taxonomy info for closest[0].
95 //so if closest[0] taxonomy is Bacteria;Alphaproteobacteria;Rhizobiales;Azorhizobium_et_rel.;Methylobacterium_et_rel.;Bosea;
96 //taxon[0][0] = Bacteria, taxon[0][1] = Alphaproteobacteria....
98 taxons.resize(closest.size());
101 for (int i = 0; i < closest.size(); i++) {
102 if (m->control_pressed) { return "control"; }
104 string tax = taxonomy[closest[i]]; //we know its there since we checked in getTaxonomy
107 taxons[i] = parseTax(tax);
109 //figure out who has the shortest taxonomy info. so you can start comparing there
110 if (taxons[i].size() < smallest) {
111 smallest = taxons[i].size();
115 //start at the highest level all the closest seqs have
117 for (int i = (smallest-1); i >= 0; i--) {
118 if (m->control_pressed) { return "control"; }
120 string thistax = taxons[0][i];
122 for (int j = 1; j < taxons.size(); j++) {
123 if (taxons[j][i] != thistax) { break; }
127 if (num == (taxons.size()-1)) { //they all match at this level
128 for (int k = 0; k <= i; k++) {
129 common += taxons[0][k] + ';';
137 //create a tree containing sequences from this bin
138 PhyloTree* p = new PhyloTree();
140 for (int i = 0; i < closest.size(); i++) {
141 p->addSeqToTree(closest[i], taxonomy[closest[i]]);
145 p->assignHeirarchyIDs(0);
147 TaxNode currentNode = p->get(0);
150 while (currentNode.children.size() != 0) { //you still have more to explore
153 int bestChildSize = 0;
155 //go through children
156 for (map<string, int>::iterator itChild = currentNode.children.begin(); itChild != currentNode.children.end(); itChild++) {
158 TaxNode temp = p->get(itChild->second);
160 //select child with largest accessions - most seqs assigned to it
161 if (temp.accessions.size() > bestChildSize) {
162 bestChild = p->get(itChild->second);
163 bestChildSize = temp.accessions.size();
168 if (bestChildSize == closest.size()) { //if yes, add it
169 conTax += bestChild.name + ";";
175 currentNode = bestChild;
182 catch(exception& e) {
183 m->errorOut(e, "Knn", "findCommonTaxonomy");
187 /**************************************************************************************************/