2597212350

Build Type

github

Committed by GitHub

Commit Message

<a href="https://github.com/sequana/sequana/commit/<a class=hub.com/sequana/sequana/commit/<a class="double-link" href="https://git"><a class=hub.com/sequana/sequana/commit/<a class="double-link" href="https://git"><a class=hub.com/sequana/sequana/commit/9a8f508f790608d15aa087c06a0073cb75328557">9a8f508f7<a href="https://github.com/sequana/sequana/commit/9a8f508f790608d15aa087c06a0073cb75328557">&quot;&gt;&amp;lt;a href=&amp;quot;https://github.com/sequana/sequana/commit/&lt;/a&gt;&lt;a class=&quot;double-link&quot; href=&quot;https://github.com/sequana/sequana/commit/&amp;lt;a class=&amp;quot;double-link&amp;quot; href=&amp;quot;https://git&quot;&gt;&amp;lt;a class=&lt;/a&gt;hub.com/sequana/sequana/commit/9a8f508f790608d15aa087c06a0073cb75328557&quot;&gt;9a8f508f7&lt;/a&gt;&lt;a href=&quot;https://github.com/sequana/sequana/commit/9a8f508f790608d15aa087c06a0073cb75328557&quot;&gt;&amp;lt;a href=&amp;quot;https://github.com/sequana/sequana/commit/9a8f508f790608d15aa087c06a0073cb&lt;/a&gt;75328557">&quot;&gt;Merge &lt;/a&gt;&lt;a class=&quot;double-link&quot; href=&quot;https://github.com/sequana/sequana/commit/&lt;a class=&quot;double-link&quot; href=&quot;https://github.com/sequana/sequana/commit/<a class="double-link" href="https://github.com/sequana/sequana/commit/6ea2104ef36649e07fcb5967ca0b12c1bd93d1ad">6ea2104ef</a><a href="https://github.com/sequana/sequana/commit/9a8f508f790608d15aa087c06a0073cb75328557">&amp;quot;&amp;gt;6ea2104ef&amp;lt;/a&amp;gt;&amp;quot;&amp;gt;6ea2104ef&amp;lt;/a&amp;gt;&amp;lt;a href=&amp;quot;https://github.com/sequana/sequana/commit/9a8f508f790608d15aa087c06a0073cb75328&lt;/a&gt;557&quot;&gt; into &lt;/a&gt;&lt;a class=&quot;double-link&quot; href=&quot;https://github.com/sequana/sequana/commit/<a class="double-link" href="https://github.com/sequana/sequana/commit/01b52271f">01b52271f">01b52271f</a>

Pull Request Pull Request #754: Fix incomplete annotation for kegg enrichment

Run Details

39 of 39 new or added lines in 8 files covered. (100.0%)

11997 of 17140 relevant lines covered (69.99%)

2.09 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

0.0

/sequana/scripts/main/enrichment_uniprot.py

#  This file is part of Sequana software
#
#  Copyright (c) 2016-2022 - Sequana Development Team
#
#  Distributed under the terms of the 3-clause BSD license.
#  The full license is in the LICENSE file, distributed with this software.
#
#  website: https://github.com/sequana/sequana
#  documentation: http://sequana.readthedocs.io
#
##############################################################################
import os
import sys
import json
import itertools

import click
import colorlog

from sequana.enrichment.uniprot_enrichment import UniprotEnrichment
from sequana.modules_report.uniprot_enrichment import ModuleUniprotEnrichment
from sequana.modules_report import ModulePantherEnrichment
from sequana.rnadiff import RNADiffResults
from sequana.utils import config

from .utils import CONTEXT_SETTINGS, common_logger, OptionEatAll


logger = colorlog.getLogger(__name__)


# This will be a complex command to provide HTML summary page for
# input files (e.g. bam), or results from pipelines. For each module,
# we should have corresponding option that starts with the module's name
# This can also takes as input various types of data (e.g. FastA)
@click.command(context_settings=CONTEXT_SETTINGS)
@click.argument("name", type=click.Path(exists=True), nargs=1)
@click.option(
    "--annotation-attribute",
    type=click.STRING,
    # required=True,
    default="index",
    show_default=True,
    help="a valid taxon identifiers",
)
@click.option(
    "--taxon",
    type=click.INT,
    required=True,
    help="a valid taxon identifiers",
)
@click.option(
    "--log2-foldchange-cutoff",
    type=click.FLOAT,
    default=1,
    show_default=True,
    help="remove events with absolute log2 fold change below this value",
)
@click.option(
    "--padj-cutoff",
    type=click.FLOAT,
    default=0.05,
    show_default=True,
    help="remove events with pvalue abobe this value default (0.05).",
)
@click.option(
    "--plot-linearx",
    type=click.BOOL,
    default=False,
    is_flag=True,
    show_default=True,
    help="""Default is log2 fold enrichment in the plots. use this to use linear scale""",
)
@click.option(
    "--compute-levels/--no-compute-levels",
    default=True,
    help="""Compute the levels of each go term, set --no-compute-levels to skip this step""",
)
@click.option(
    "--max-genes",
    type=click.INT,
    default=2500,
    show_default=True,
    help="""Maximum number of genes (up or down) to use in PantherDB.""",
)
@click.option(
    "--ontologies",
    default=("MF", "BP", "CC"),
    help="""Provide the ontologies to be included in the analysis and HTML report.
Valid choices are: from MF, BP, CC, SLIM_MF, SLIM_BP, SLIM_CC, PROTEIN,
PANTHER_PATHWAY, REACTOME_PATHWAY""",
    cls=OptionEatAll,
    show_default=True,
)
@click.option(
    "--max-enriched-go-terms",
    type=click.INT,
    default=40,
    show_default=True,
    help="""Max number of enriched go terms to show in the plots (most
enriched). All enriched GO terms are stored in tables""",
)
@click.option("--output-directory", show_default=True, default="enrichment_uniprot")
@common_logger
def enrichment_uniprot(**kwargs):
    """GO enrichment using Uniprot and HTML report creation

    \b
    The input should be the output of the RNADiff pipeline

    Example::

        sequana enrichment-uniprot rnadiff.csv --taxon 10090
            --log2-foldchange-cutoff 2 

        sequana enrichment rnadiff/rnadiff.csv
            --taxon 189518 \
            --log2-foldchange-cutoff 2 
            --ontologies MF 

    \b
    Valid ontologies are: MF, BP, CC


    """
    valid = {"MF", "BP",  "CC"}

    ontologies = eval(kwargs["ontologies"])
    for ontology in ontologies:
        if ontology not in valid:
            logger.error(f"Provided incorrect ontology ({ontology}). Must be in {valid}")
            sys.exit(1)

    logger.setLevel(kwargs["logger"])

    taxon = kwargs["taxon"]
    if taxon == 0:
        logger.error("You must provide a taxon with --taxon")
        sys.exit(1)

    params = {
        "padj": kwargs["padj_cutoff"],
        "log2_fc": kwargs["log2_foldchange_cutoff"],
        "max_entries": kwargs["max_genes"],
        "nmax": kwargs["max_enriched_go_terms"],
        "plot_logx": not kwargs["plot_linearx"],
        "plot_compute_levels": kwargs["compute_levels"],
    }

    logger.info(f"Reading RNAdiff results from {kwargs['name']}")
    dirpath = os.path.dirname(os.path.abspath(kwargs["name"]))
    rnadiff = RNADiffResults(dirpath, index_col=0, header=[0, 1])

    # now that we have loaded all results from a rnadiff analysis, let us
    # perform the enrichment for each comparison found in the file
    annot_col = kwargs.get("annotation_attribute", "index")

    logger.info(f"Using the annotation column '{annot_col}'")

    # setting these attributes set the gene list with log2fc and padj filter
    rnadiff._log2_fc = params["log2_fc"]
    rnadiff._alpha = params["padj"]
    gene_lists = rnadiff.get_gene_lists(annot_col=annot_col, Nmax=kwargs.get("max_genes", None), dropna=True)


    output_directory = kwargs["output_directory"]
    os.makedirs(output_directory, exist_ok=True)

    for compa, gene_dict in gene_lists.items():
        config.output_dir = f"{output_directory}/{compa}"
        os.makedirs(f"{output_directory}/{compa}", exist_ok=True)

        ue = UniprotEnrichment(gene_dict, 235443)
        ue.compute_enrichment()

        for x in itertools.product(['up', 'down', 'all'], ['CC', 'BP', 'MF']):
            mode, ont = x
            df = ue.plot_go_terms(mode, ontologies=[ont])
            if df is not None and len(df):
                from pylab import savefig
                savefig(f"{output_directory}/{compa}/plot_{mode}_{ont}.png", dpi=200)
                ue.save_chart(df, f'{output_directory}/{compa}/chart_{mode}_{ont}.png')
                df.to_csv(f'{output_directory}/{compa}/{mode}_{ont}.csv', sep=",")


        stats = {"taxon": taxon}
        stats["taxon_name"] = ue.taxon_name
        stats["up"] = len(gene_dict['up'])
        stats["down"] = len(gene_dict['down'])
        stats["all"] = len(gene_dict['all'])
        stats["name"] = compa
        with open(f"{output_directory}/{compa}/summary.json", "w") as fout:
            json.dump(stats, fout, indent=True)


        ModuleUniprotEnrichment(
            gene_dict,
            stats,
            enrichment_params=params,
            command=" ".join(["sequana"] + sys.argv[1:]),
             ontologies=["CC", "BP", "MF"],
        )

1	# This file is part of Sequana software
2	#
3	# Copyright (c) 2016-2022 - Sequana Development Team
4	#
5	# Distributed under the terms of the 3-clause BSD license.
6	# The full license is in the LICENSE file, distributed with this software.
7	#
8	# website: https://github.com/sequana/sequana
9	# documentation: http://sequana.readthedocs.io
10	#
11	##############################################################################
12	import os	×
13	import sys	×
14	import json	×
15	import itertools	×
16
17	import click	×
18	import colorlog	×
19
20	from sequana.enrichment.uniprot_enrichment import UniprotEnrichment	×
21	from sequana.modules_report.uniprot_enrichment import ModuleUniprotEnrichment	×
22	from sequana.modules_report import ModulePantherEnrichment	×
23	from sequana.rnadiff import RNADiffResults	×
24	from sequana.utils import config	×
25
26	from .utils import CONTEXT_SETTINGS, common_logger, OptionEatAll	×
27
28
29	logger = colorlog.getLogger(__name__)	×
30
31
32	# This will be a complex command to provide HTML summary page for
33	# input files (e.g. bam), or results from pipelines. For each module,
34	# we should have corresponding option that starts with the module's name
35	# This can also takes as input various types of data (e.g. FastA)
36	@click.command(context_settings=CONTEXT_SETTINGS)	×
37	@click.argument("name", type=click.Path(exists=True), nargs=1)	×
38	@click.option(	×
39	"--annotation-attribute",
40	type=click.STRING,
41	# required=True,
42	default="index",
43	show_default=True,
44	help="a valid taxon identifiers",
45	)
46	@click.option(	×
47	"--taxon",
48	type=click.INT,
49	required=True,
50	help="a valid taxon identifiers",
51	)
52	@click.option(	×
53	"--log2-foldchange-cutoff",
54	type=click.FLOAT,
55	default=1,
56	show_default=True,
57	help="remove events with absolute log2 fold change below this value",
58	)
59	@click.option(	×
60	"--padj-cutoff",
61	type=click.FLOAT,
62	default=0.05,
63	show_default=True,
64	help="remove events with pvalue abobe this value default (0.05).",
65	)
66	@click.option(	×
67	"--plot-linearx",
68	type=click.BOOL,
69	default=False,
70	is_flag=True,
71	show_default=True,
72	help="""Default is log2 fold enrichment in the plots. use this to use linear scale""",
73	)
74	@click.option(	×
75	"--compute-levels/--no-compute-levels",
76	default=True,
77	help="""Compute the levels of each go term, set --no-compute-levels to skip this step""",
78	)
79	@click.option(	×
80	"--max-genes",
81	type=click.INT,
82	default=2500,
83	show_default=True,
84	help="""Maximum number of genes (up or down) to use in PantherDB.""",
85	)
86	@click.option(	×
87	"--ontologies",
88	default=("MF", "BP", "CC"),
89	help="""Provide the ontologies to be included in the analysis and HTML report.
90	Valid choices are: from MF, BP, CC, SLIM_MF, SLIM_BP, SLIM_CC, PROTEIN,
91	PANTHER_PATHWAY, REACTOME_PATHWAY""",
92	cls=OptionEatAll,
93	show_default=True,
94	)
95	@click.option(	×
96	"--max-enriched-go-terms",
97	type=click.INT,
98	default=40,
99	show_default=True,
100	help="""Max number of enriched go terms to show in the plots (most
101	enriched). All enriched GO terms are stored in tables""",
102	)
103	@click.option("--output-directory", show_default=True, default="enrichment_uniprot")	×
104	@common_logger	×
105	def enrichment_uniprot(**kwargs):
106	"""GO enrichment using Uniprot and HTML report creation
107
108	\b
109	The input should be the output of the RNADiff pipeline
110
111	Example::
112
113	sequana enrichment-uniprot rnadiff.csv --taxon 10090
114	--log2-foldchange-cutoff 2
115
116	sequana enrichment rnadiff/rnadiff.csv
117	--taxon 189518 \
118	--log2-foldchange-cutoff 2
119	--ontologies MF
120
121	\b
122	Valid ontologies are: MF, BP, CC
123
124
125	"""
126	valid = {"MF", "BP", "CC"}	×
127
128	ontologies = eval(kwargs["ontologies"])	×
129	for ontology in ontologies:	×
130	if ontology not in valid:	×
131	logger.error(f"Provided incorrect ontology ({ontology}). Must be in {valid}")	×
132	sys.exit(1)	×
133
134	logger.setLevel(kwargs["logger"])	×
135
136	taxon = kwargs["taxon"]	×
137	if taxon == 0:	×
138	logger.error("You must provide a taxon with --taxon")	×
139	sys.exit(1)	×
140
141	params = {	×
142	"padj": kwargs["padj_cutoff"],
143	"log2_fc": kwargs["log2_foldchange_cutoff"],
144	"max_entries": kwargs["max_genes"],
145	"nmax": kwargs["max_enriched_go_terms"],
146	"plot_logx": not kwargs["plot_linearx"],
147	"plot_compute_levels": kwargs["compute_levels"],
148	}
149
150	logger.info(f"Reading RNAdiff results from {kwargs['name']}")	×
151	dirpath = os.path.dirname(os.path.abspath(kwargs["name"]))	×
152	rnadiff = RNADiffResults(dirpath, index_col=0, header=[0, 1])	×
153
154	# now that we have loaded all results from a rnadiff analysis, let us
155	# perform the enrichment for each comparison found in the file
156	annot_col = kwargs.get("annotation_attribute", "index")	×
157
158	logger.info(f"Using the annotation column '{annot_col}'")	×
159
160	# setting these attributes set the gene list with log2fc and padj filter
161	rnadiff._log2_fc = params["log2_fc"]	×
162	rnadiff._alpha = params["padj"]	×
163	gene_lists = rnadiff.get_gene_lists(annot_col=annot_col, Nmax=kwargs.get("max_genes", None), dropna=True)	×
164
165
166	output_directory = kwargs["output_directory"]	×
167	os.makedirs(output_directory, exist_ok=True)	×
168
169	for compa, gene_dict in gene_lists.items():	×
170	config.output_dir = f"{output_directory}/{compa}"	×
171	os.makedirs(f"{output_directory}/{compa}", exist_ok=True)	×
172
173	ue = UniprotEnrichment(gene_dict, 235443)	×
174	ue.compute_enrichment()	×
175
176	for x in itertools.product(['up', 'down', 'all'], ['CC', 'BP', 'MF']):	×
177	mode, ont = x	×
178	df = ue.plot_go_terms(mode, ontologies=[ont])	×
179	if df is not None and len(df):	×
180	from pylab import savefig	×
181	savefig(f"{output_directory}/{compa}/plot_{mode}_{ont}.png", dpi=200)	×
182	ue.save_chart(df, f'{output_directory}/{compa}/chart_{mode}_{ont}.png')	×
183	df.to_csv(f'{output_directory}/{compa}/{mode}_{ont}.csv', sep=",")	×
184
185
186	stats = {"taxon": taxon}	×
187	stats["taxon_name"] = ue.taxon_name	×
188	stats["up"] = len(gene_dict['up'])	×
189	stats["down"] = len(gene_dict['down'])	×
190	stats["all"] = len(gene_dict['all'])	×
191	stats["name"] = compa	×
192	with open(f"{output_directory}/{compa}/summary.json", "w") as fout:	×
193	json.dump(stats, fout, indent=True)	×
194
195
196	ModuleUniprotEnrichment(	×
197	gene_dict,
198	stats,
199	enrichment_params=params,
200	command=" ".join(["sequana"] + sys.argv[1:]),
201	ontologies=["CC", "BP", "MF"],
202	)

sequana / sequana / 2597212350

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous