

@transform(buildLCA, suffix('.taxa.gz'), '.taxa.readcounts')
def countContributingReads(infile, outfile):
    '\n    count number of reads with a taxnomic assignment\n    '
    levels = ['phylum', 'class', 'order', 'family', 'genus', 'species']
    result = collections.OrderedDict()
    for level in levels:
        result[level] = 0
    inf = IOTools.openFile(infile)
    header = inf.readline().split('\t')
    indices = [3, 5, 7, 9, 11, 13]
    total = 0
    for line in inf.readlines():
        total += 1
        data = line[:(- 1)].split('\t')
        (phylum, _class, order, family, genus, species) = [data[i] for i in indices]
        if (phylum != 'NA'):
            result['phylum'] += 1
        if (_class != 'NA'):
            result['class'] += 1
        if (order != 'NA'):
            result['order'] += 1
        if (family != 'NA'):
            result['family'] += 1
        if (genus != 'NA'):
            result['genus'] += 1
        if (species != 'NA'):
            result['species'] += 1
    outf = open(outfile, 'w')
    outf.write('level\tn_reads\tpct_reads\n')
    for (level, count) in _class.iteritems():
        (nreads, prop) = (count, (float(count) / total))
        outf.write(('\t'.join([level, str(nreads), str((prop * 100))]) + '\n'))
    outf.close()

dataset/ETHPy150Open CGATOxford/cgat/obsolete/pipeline_metagenomecommunities.py countContributingReads/VarMisuse@290/407 `result`->`_class`