Commit 81479249 by Eric Coissac

--no commit message

parent 63bf2cfe
......@@ -4,25 +4,18 @@ import math
import sys
from obitools.ecopcr import taxonomy
from obitools.ecopcr import EcoPCRFile
from obitools.graph import Graph
from obitools.utils import progressBar
from obitools.align import LCS
from obitools.align import isLCSReachable
from obitools.format.options import addInOutputOption, sequenceWriterGenerator
from obitools.options import getOptionManager
from obitools.graph.algorithms.component import componentIterator
from obitools.ecopcr.options import addTaxonomyDBOptions, loadTaxonomyDatabase
def addTaxonomyOptions(optionManager):
optionManager.add_option('-d','--ecopcrdb',
action="store", dest="db",
metavar="<FILENAME>",
type="string",
help="ecoPCR Database "
"name")
def addSpecificityOptions(optionManager):
optionManager.add_option('-e','--errors',
action="store", dest="dist",
metavar="###",
......@@ -38,15 +31,15 @@ def addTaxonomyOptions(optionManager):
if __name__=='__main__':
optionParser = getOptionManager([addTaxonomyOptions],
)
optionParser = getOptionManager([addTaxonomyDBOptions,addSpecificityOptions])
(options, entries) = optionParser()
tax = taxonomy.EcoTaxonomyDB(options.db)
loadTaxonomyDatabase(options)
tax =options.taxonomy
ranks = set(x for x in tax.rankIterator())
results = [seq for seq in entries]
# dist = 0
graph = Graph("error",directed=False)
xx = 0
......@@ -103,7 +96,7 @@ if __name__=='__main__':
yy = 0
for c in componentIterator(graph):
sub = graph.subgraph(c)
sub = graph.subgraph(c)
first = True
s = ""
for node in sub: #all nodes of a component should go with same key (taken as first sequence in comp)
......@@ -151,31 +144,31 @@ if __name__=='__main__':
if tag in indexbykey:
key = indexbykey[tag] #get component key for this seq
if options.quorum > 0.0:
for x in indexbyseq[key]:
txn = tax.getTaxonAtRank(x['taxid'],rankid)
if txn not in taxonindividuals:
taxonindividuals[txn] = set([x['taxid']])
else:
taxonindividuals[txn].add(x['taxid'])
for x in indexbyseq[key]:
txn = tax.getTaxonAtRank(x['taxid'],rankid)
if txn not in taxonindividuals:
taxonindividuals[txn] = set([x['taxid']])
else:
taxonindividuals[txn].add(x['taxid'])
taxlist |=set(tax.getTaxonAtRank(x['taxid'],rankid) for x in indexbyseq[key])
else:
if tag not in lostSeqs:
lostSeqs.append(tag)
taxoncount+=1
if options.quorum > 0.0:
max = 0
sum = 0
for k in taxonindividuals.keys ():
if len(taxonindividuals[k]) > max:
max = len(taxonindividuals[k])
sum = sum + len(taxonindividuals[k])
if max >= (sum-sum*options.quorum):
taxonok += 1
max = 0
sum = 0
for k in taxonindividuals.keys ():
if len(taxonindividuals[k]) > max:
max = len(taxonindividuals[k])
sum = sum + len(taxonindividuals[k])
if max >= (sum-sum*options.quorum):
taxonok += 1
else:
if len(taxlist)==1:
taxonok+=1
if len(taxlist)==1:
taxonok+=1
if taxoncount:
print '%-20s\t%10d\t%10d\t%8.2f' % (rank,taxonok,taxoncount,float(taxonok)/taxoncount*100)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment