Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def get_goeaobj(methods=None):
"""Test GOEA with method, fdr."""
obo_dag = GODag(ROOT + "goslim_generic.obo")
fin_assc = ROOT + "slim_association"
assoc = read_associations(fin_assc, 'id2gos', no_top=True)
popul_ids = [line.rstrip() for line in open(ROOT + "small_population")]
goeaobj = GOEnrichmentStudy(popul_ids, assoc, obo_dag, methods=methods)
return goeaobj
def _get_results(godag, propagate_counts, relationships, prt=sys.stdout):
"""Run a GOEA. Return results"""
taxid = 10090 # Mouse study
geneids_pop = set(GeneID2nt_mus.keys())
assoc_geneid2gos = get_assoc_ncbi_taxids([taxid], loading_bar=None)
geneids_study = get_geneid2symbol("nbt.3102-S4_GeneIDs.xlsx")
goeaobj = GOEnrichmentStudy(
geneids_pop,
assoc_geneid2gos,
godag,
propagate_counts=propagate_counts,
relationships=relationships,
alpha=0.05,
methods=['fdr_bh'])
return goeaobj.run_study(geneids_study, prt=prt)
"""Test to re-produce issue#96: Passes currently."""
# Trying to duplicate: ValueError("All values in table must be nonnegative.
# Get genes
print('CWD', os.getcwd())
study_ids = _get_geneids()
population_ids = GENEID2NT.keys()
# Get databases
print(os.getcwd())
fin = os.path.join(REPO, 'gene2go')
dnld_ncbi_gene_file(fin, loading_bar=None)
gene2go = read_ncbi_gene2go(fin, [9606])
fin_obo = os.path.join(REPO, "go-basic.obo")
godag = get_godag(fin_obo, loading_bar=None)
goeaobj = GOEnrichmentStudy(population_ids, gene2go, godag, methods=['fdr_bh'])
# Run GOEA Gene Ontology Enrichment Analysis
results_goeas = goeaobj.run_study(study_ids)
def test_i122():
"""Test to re-produce issue#122: Passes currently."""
obj = _Run(9606, 'gene2go', 'go-basic.obo')
study_ids, population_ids = obj.get_genes_study_n_bg()
# Result is the same whether fisher_scipy_stats of fisher
pvalcalc = 'fisher_scipy_stats'
goeaobj = GOEnrichmentStudy(population_ids, obj.gene2go, obj.godag, methods=['bonferroni', 'fdr_bh'], pvalcalc=pvalcalc)
# Run GOEA Gene Ontology Enrichment Analysis
results_goeas = goeaobj.run_study_nts(study_ids)
print('NS GO p stu_ratio pop_ratio p-uncorr bonferro fdr_bh stu ')
for ntd in results_goeas:
if ntd.study_count == 0:
doprt = False
if ntd.p_bonferroni < 0.05:
assert ntd.enrichment == 'p'
doprt = True
if ntd.p_fdr_bh < 0.05:
assert ntd.enrichment == 'p'
doprt = True
if doprt:
print(obj.str_nt(ntd))
# print(next(iter(results_goeas))._fields)
def _get_pvals(pvalfnc_names, prt=sys.stdout):
fisher2pvals = {}
taxid = 10090 # Mouse study
file_obo = os.path.join(os.getcwd(), "go-basic.obo")
obo_dag = get_godag(file_obo, prt, loading_bar=None)
geneids_pop = set(GeneID2nt_mus.keys())
assoc_geneid2gos = get_assoc_ncbi_taxids([taxid], loading_bar=None)
geneids_study = get_geneid2symbol("nbt.3102-S4_GeneIDs.xlsx")
for fisher in pvalfnc_names:
goeaobj = GOEnrichmentStudy(
geneids_pop,
assoc_geneid2gos,
obo_dag,
propagate_counts=False,
alpha=0.05,
methods=None,
pvalcalc=fisher)
fisher2pvals[fisher] = goeaobj.get_pval_uncorr(geneids_study, prt)
return fisher2pvals
def run_bonferroni():
"""Do Gene Ontology Enrichment Analysis w/Bonferroni multipletest. Print results 3 ways."""
# ---------------------------------------------------------------------
# Run Gene Ontology Analysis (GOEA)
#
# 1. Initialize
godag = get_godag(os.path.join(os.getcwd(), "go-basic.obo"), loading_bar=None)
fin_assc = os.path.join(REPO, "data/association")
assoc = read_associations(fin_assc, 'id2gos', no_top=True)
popul_ids = [line.rstrip() for line in open(os.path.join(REPO, "data/population"))]
study_ids = [line.rstrip() for line in open(os.path.join(REPO, "data/study"))]
# 2. Run enrichment analysis
goea = GOEnrichmentStudy(popul_ids, assoc, godag, alpha=0.05, methods=['bonferroni'])
results_nt = goea.run_study(study_ids)
return results_nt, goea
def get_goea_results(method="fdr_bh"):
"""Get GOEA results."""
root_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), "data")
obo_fin = os.path.join(root_dir, "goslim_generic.obo")
obo_dag = GODag(obo_fin)
fin_slim = os.path.join(root_dir, "slim_association")
assoc = read_associations(fin_slim, 'id2gos', no_top=True)
popul_ids = [line.rstrip() for line in open(os.path.join(root_dir, "small_population"))]
goeaobj = GOEnrichmentStudy(popul_ids, assoc, obo_dag, methods=[method])
study_ids = [line.rstrip() for line in open(os.path.join(root_dir, "small_study"))]
goea_results = goeaobj.run_study(study_ids, methods=[method])
return goea_results
def init_goea(**kws):
"""Initialize GODag and GOEnrichmentStudy."""
godag = get_godag(os.path.join(os.getcwd(), "go-basic.obo"), loading_bar=None)
fin_assc = ROOT + "association"
assoc = read_associations(fin_assc, 'id2gos', no_top=True)
popul_ids = [line.rstrip() for line in open(ROOT + "population")]
methods = kws['methods'] if 'methods' in kws else ['not_bonferroni']
study_ids = [line.rstrip() for line in open(ROOT + "study")]
return GOEnrichmentStudy(popul_ids, assoc, godag, methods=methods), study_ids
def _ns2o(pop, ns2assoc, godag, propagate_counts, alpha, methods, **kws):
return {
ns:GOEnrichmentStudy(pop, a, godag, propagate_counts, alpha, methods, name=ns, **kws) \
for ns, a in sorted(ns2assoc.items())}