Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Refactor universal enrichment #92

Merged
merged 2 commits into from
Aug 29, 2023
Merged

Refactor universal enrichment #92

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
06c6e03
Refactor universal enrichment logic following clusterProfiler's enricher
memgonzales Aug 29, 2023
acc554d
Merge branch 'main' of https://github.com/bioinfodlsu/rice-pilaf into…
memgonzales Aug 29, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 19 additions & 13 deletions callbacks/coexpression/util.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -190,27 +190,33 @@ def do_module_enrichment_analysis(implicated_gene_ids, genomic_intervals, addl_g
# provided by clusterProfiler
# ====================================================================================
with open(IMPLICATED_GENES_PATH) as implicated_genes_file, open(MODULES_PATH) as modules_file, open(f'{INPUT_GENES_DIR}/enriched_modules.tsv', 'w') as enriched_modules_file:
background_genes = set()
for line in modules_file:
background_genes = background_genes.union(
set(line.strip().split('\t')))

# There is only a single line, which lists all the implicated genes
for line in implicated_genes_file:
line = line.strip().split('\t')
implicated_genes = set(line)

modules = []
background_genes = set()
for idx, line in enumerate(modules_file):
module_genes = set(line.strip().split('\t'))
background_genes = background_genes.union(module_genes)
if implicated_genes.intersection(module_genes):
modules.append(idx)

p_values_indices = []
p_values = []
modules_file.seek(0)
for idx, line in enumerate(modules_file):
module = line.strip().split('\t')
module_genes = set(module)
table = construct_contigency_table(
background_genes, implicated_genes, module_genes)

p_value = fisher_exact(table, alternative='greater').pvalue
if not (0.999999999 < p_value and p_value < 1.000000001):
p_values.append(p_value)
if idx in modules:
module = line.strip().split('\t')
module_genes = set(module)
table = construct_contigency_table(
background_genes, implicated_genes, module_genes)

p_values.append(fisher_exact(
table, alternative='greater').pvalue)

# Add 1 since user-facing module number is one-based
p_values_indices.append(idx + 1)

_, adj_p_values, _, _ = sm.multipletests(
Expand Down