@transform(buildLCA, suffix(".taxa.gz"), ".taxa.readcounts")
def countContributingReads(infile, outfile):
    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")

    # column indices
    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 result.iteritems():
        nreads, prop = count, float(count) / total
        outf.write("\t".join([level, str(nreads), str(prop * 100)]) + "\n")
    outf.close()

###################################################################
###################################################################
###################################################################
count number of reads with a taxnomic assignment
dataset/ETHPy150Open CGATOxford/cgat/obsolete/pipeline_metagenomecommunities.py/countContributingReads