<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Calculating Annotation Coverage
This section shows how to calculate annotation coverage as described here
Step1:  2. Read associations
2a. You can read the associations one species at a time...
Step2:  2b. Or you can read 'gene2go' once and load all species...
Step3:  3. Import protein-coding information for human and fly
Step4:  4. Calculate Gene Ontology coverage
Store GO coverage information for human and fly in the list, cov_data.
Step5:  5 Report Gene Ontology coverage for human and fly
Print the human and fly GO coverage information that is stored in the list, cov_data.
<ASSISTANT_TASK:>
Python Code:
# Get ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/gene2go.gz
from goatools.base import download_ncbi_associations
gene2go = download_ncbi_associations()
Explanation: Calculating Annotation Coverage
This section shows how to calculate annotation coverage as described here:    
Annotation coverage of Gene Ontology (GO) terms to individual 
gene products is high for human or model organisms: 
  * 87% of ~20k human protein-coding genes have GO annotations
  * 76% of ~14k fly protein-coding genes have GO annotations
(Apr 27, 2016)
1. Download associations
NCBI's gene2go file contains annotations of GO terms to Entrez GeneIDs for over 35 different species. We are interested in human and fly which have the taxids 9606 and 7227 repectively.
End of explanation
from goatools.associations import read_ncbi_gene2go
geneid2gos_human = read_ncbi_gene2go(gene2go, taxids=[9606])
geneid2gos_fly = read_ncbi_gene2go(gene2go, taxids=[7227])
Explanation: 2. Read associations
2a. You can read the associations one species at a time...
End of explanation
from collections import defaultdict, namedtuple
taxid2asscs = defaultdict(lambda: defaultdict(lambda: defaultdict(set)))
geneid2gos_all = read_ncbi_gene2go(
        gene2go, 
        taxids=[9606, 7227], 
        taxid2asscs=taxid2asscs)
Explanation: 2b. Or you can read 'gene2go' once and load all species...
End of explanation
from goatools.test_data.genes_NCBI_9606_ProteinCoding import GeneID2nt as GeneID2nt_human
from goatools.test_data.genes_NCBI_7227_ProteinCoding import GeneID2nt as GeneID2nt_fly
lst = [
    (9606, GeneID2nt_human),
    (7227, GeneID2nt_fly)
]
Explanation: 3. Import protein-coding information for human and fly
End of explanation
cov_data = []
NtCov = namedtuple("NtCov", "taxid num_GOs num_covgenes coverage num_allgenes")
for taxid, pcGeneID2nt in lst:
    # Get GeneID2GOs association for current species
    geneid2gos = taxid2asscs[taxid]['GeneID2GOs']
    # Restrict GeneID2GOs to only protein-coding genes for this report
    pcgene_w_gos = set(geneid2gos.keys()).intersection(set(pcGeneID2nt.keys()))
    num_pcgene_w_gos = len(pcgene_w_gos)
    num_pc_genes = len(pcGeneID2nt)
    # Number of GO terms annotated to protein-coding genes
    gos_pcgenes = set()
    for geneid in pcgene_w_gos:
        gos_pcgenes |= geneid2gos[geneid]
    # Print report data
    cov_data.append(NtCov(
        taxid = taxid,
        num_GOs = len(gos_pcgenes),
        num_covgenes = num_pcgene_w_gos,
        coverage = 100.0*num_pcgene_w_gos/num_pc_genes,
        num_allgenes = num_pc_genes))
Explanation: 4. Calculate Gene Ontology coverage
Store GO coverage information for human and fly in the list, cov_data.
End of explanation
from __future__ import print_function
print(" taxid    GOs GeneIDs Coverage")
print("------ ------ ------- ----------------------")
fmtstr = "{TAXID:>6} {N:>6,} {M:>7,}  {COV:2.0f}% GO coverage of {TOT:,} protein-coding genes"
for nt in cov_data:
    print(fmtstr.format(
        TAXID = nt.taxid,
        N = nt.num_GOs,
        M = nt.num_covgenes,
        COV = nt.coverage,
        TOT = nt.num_allgenes))
Explanation: 5 Report Gene Ontology coverage for human and fly
Print the human and fly GO coverage information that is stored in the list, cov_data.
End of explanation
<END_TASK>