1261710972

Build Type

github

Committed by GitHub

Commit Message

<a href="https://github.com/sequana/sequana/commit/<a class=hub.com/sequana/sequana/commit/60577b049802bcbf3dc6109db02723148d9ae166">60577b049<a href="https://github.com/sequana/sequana/commit/60577b049802bcbf3dc6109db02723148d9ae166">">Merge </a><a class="double-link" href="https://github.com/sequana/sequana/commit/<a class="double-link" href="https://github.com/sequana/sequana/commit/f0a7e5535543bf0d5f508c121bdd5778f4cade4c">f0a7e5535</a>">f0a7e5535</a><a href="https://github.com/sequana/sequana/commit/60577b049802bcbf3dc6109db02723148d9ae166"> into 8bdd94f22">8bdd94f22</a>

Pull Request Pull Request #678: Improvements related to enrichment

Run Details

83 of 83 new or added lines in 3 files covered. (100.0%)

11605 of 17202 relevant lines covered (67.46%)

1.34 hits per line

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

16.96

/sequana/modules_report/panther_enrichment.py

#
#  This file is part of Sequana software
#
#  Copyright (c) 2020 - 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
#
##############################################################################
"""Module to write enrichment report"""
import os
import sys

from sequana.lazy import pandas as pd
from sequana.lazy import pylab

from sequana.modules_report.base_module import SequanaBaseModule
from sequana.utils.datatables_js import DataTable
from sequana.enrichment.panther import PantherEnrichment
from sequana.utils import config


from easydev import Progress

import colorlog

logger = colorlog.getLogger(__name__)


class ModulePantherEnrichment(SequanaBaseModule):
    """Write HTML report of variant calling. This class takes a csv file
    generated by sequana_variant_filter.
    """

    def __init__(
        self,
        gene_lists,
        taxon,
        enrichment_params={
            "padj": 0.05,
            "log2_fc": 3,
            "max_entries": 3000,
            "nmax": 50,
            "mapper": None,
            "plot_compute_levels": False,
            "plot_logx": True,
        },
        command="",
        ontologies=["MF", "BP", "CC"],
    ):
        """.. rubric:: constructor"""
        super().__init__()
        self.title = "Panther Enrichment"

        self.command = command
        self.gene_lists = gene_lists
        self.enrichment_params = enrichment_params
        self.nmax = enrichment_params.get("nmax", 50)

        # compute the enrichment here once for all, This may take time
        from sequana import logger

        logger.setLevel("INFO")

        logger.info(" === Module PanthEnrichment. === ")
        self.pe = PantherEnrichment(
            self.gene_lists,
            taxon,
            max_entries=self.enrichment_params["max_entries"],
            log2_fc_threshold=self.enrichment_params["log2_fc"],
        )

        self.ontologies = ontologies

        # Compute the enrichment
        self.pe.compute_enrichment(ontologies=self.ontologies)
        self.df_stats = self.pe.get_mapping_stats()[
            ["category", "mapped_percentage", "total"]
        ]

        self.create_report_content()
        self.create_html("enrichment.html")

    def create_report_content(self):
        self.sections = list()
        self.summary()
        self.add_go()
        self.sections.append(
            {"name": "5 - Info", "anchor": "command", "content": self.command}
        )

    def summary(self):
        """Add information."""

        total_up = len(self.gene_lists["up"])
        total_down = len(self.gene_lists["down"])
        total = total_up + total_down
        log2fc = self.enrichment_params["log2_fc"]

        # create html table for taxon information
        _taxon_id = self.pe.taxon_info["taxon_id"]
        _taxon_name = self.pe.taxon_info["long_name"]

        # table of undertermined IDs
        df_stats = self.df_stats.drop_duplicates()
        datatable = DataTable(df_stats, "unmapped")
        datatable.datatable.datatable_options = {
            "scrollX": "true",
            "pageLength": 10,
            "scrollCollapse": "true",
            "dom": "Bfrtip",
            "buttons": ["copy", "csv"],
        }
        js = datatable.create_javascript_function()
        html_table = datatable.create_datatable(float_format="%E")

        self.sections.append(
            {
                "name": "1 - Summary",
                "anchor": "filters_option",
                "content": f"""

<p>In the following sections, you will find the GO
terms enrichment. The input data for those analyis is the output of the RNADiff
analysis where adjusted p-values above 0.05 are excluded. Moreover, we removed 
candidates with log2 fold change below {log2fc}. Using these filters, the list of
differentially expressed genes is made of {total_up} up and {total_down} down genes (total {total})</p>
<p> In the following plots you can find the first GO terms that are enriched, keeping a 
maximum of {self.nmax} identifiers. </p>

<p>The taxon used is {_taxon_name} (ID {_taxon_id}).<br>

<p> Check the following table for the percentage of mapped genes on Panther
DB</p>

{js}{html_table}
""",
            }
        )

    def add_go(self):
        # somehow, logger used here and in add_kegg must be global. If you call
        # add_go and then add_kegg, logger becomes an unbound local variable.
        # https://stackoverflow.com/questions/10851906/python-3-unboundlocalerror-local-variable-referenced-before-assignment
        # global logger
        level = logger.level
        logger.setLevel(level)

        html_intro = """
<p>Here below is a set of plots showing the enriched GO terms using the down
regulated genes only, and then the up-regulated genes only. When possible a
graph of the found GO terms is provided. MF stands for molecular
function, CC for cellular components and BP for biological process.</p>
        </div>
        """.format(
            self.pe.summary["padj_threshold"],
            self.pe.summary["fold_change_range"][0],
            self.pe.summary["fold_change_range"][1],
            self.pe.summary["DGE_after_filtering"]["down"],
            self.pe.summary["DGE_after_filtering"]["up"],
            self.pe.summary["DGE_after_filtering"]["all"],
        )

        html = self._get_enrichment("down")
        self.sections.append(
            {
                "name": "2 - Enriched GO terms (Down cases)",
                "anchor": "go_down",
                "content": html,
            }
        )

        return
        html = self._get_enrichment("up")
        self.sections.append(
            {
                "name": "3 - Enriched GO terms (Up cases)",
                "anchor": "go_up",
                "content": html,
            }
        )

        html = self._get_enrichment("all")
        self.sections.append(
            {
                "name": "4 - Enriched GO terms (All cases)",
                "anchor": "go_all",
                "content": html,
            }
        )

    # a utility function to create the proper html table
    def get_html_table(self, this_df, identifier):
        df = this_df.copy()

        # depending on the Panther category, we add the link to the identifier
        links = []
        for x in df["id"]:
            if x.startswith("PC"):
                links.append(
                    f"http://www.pantherdb.org/panther/category.do?categoryAcc={x}"
                )
            elif x.startswith("R-"):
                links.append(f"https://reactome.org/PathwayBrowser/#/{x}")
            else:
                links.append(f"https://www.ebi.ac.uk/QuickGO/term/{x}")
        df["links"] = links

        # remove non-informative or redundant fields
        for x in ["term", "fdr2", "abs_log2_fold_enrichment", "pct_diff_expr"]:
            try:
                del df[x]
            except KeyError:
                pass

        first_col = df.pop("id")
        df.insert(0, "id", first_col)
        df = df.sort_values(by="fold_enrichment", ascending=False)

        datatable = DataTable(pd.DataFrame(df), identifier)
        datatable.datatable.set_links_to_column("links", "id")
        datatable.datatable.datatable_options = {
            "scrollX": "true",
            "pageLength": 10,
            "scrollCollapse": "true",
            "dom": "Bfrtip",
            "buttons": ["copy", "csv"],
        }
        js = datatable.create_javascript_function()
        html_table = datatable.create_datatable(float_format="%E")
        return js + html_table

    def _get_enrichment(self, category):
        # category is in down/up/all

        style = "width:95%"
        _temp_df = {}
        _minus = {}
        _plus = {}

        html = ""

        for ontology in self.ontologies:

            # get dat without plotting to store the entire set of GO terms
            df, subdf = self.pe._get_plot_go_terms_data(
                category, ontologies=ontology, compute_levels=False
            )

            # now plotting but showing only some restricted GO terms
            def plot_go_terms(filename, category, ontology):
                df = self.pe.plot_go_terms(
                    category,
                    ontologies=ontology,
                    compute_levels=self.enrichment_params["plot_compute_levels"],
                    log=self.enrichment_params["plot_logx"],
                    max_features=self.nmax,
                )
                _temp_df[ontology] = df.copy()
                _plus[ontology] = sum(df.plus_minus == "+")
                _minus[ontology] = sum(df.plus_minus == "-")
                pylab.savefig(f"{config.output_dir}/Panther_{category}_{ontology}.png")
                pylab.savefig(filename)
                pylab.close()

            if len(df):
                image = self.create_embedded_png(
                    plot_go_terms,
                    "filename",
                    style=style,
                    category=category,
                    ontology=ontology,
                )

                html += f"""
<h4>{category.title()} - {ontology}</h4>
<p>For {ontology}, we found {_plus[ontology]+_minus[ontology]} go terms.
Showing {self.nmax} here below (at most). The full list is downlodable from the CSV
 file hereafter.</p> {image} <br>"""
                html += self.get_html_table(
                    _temp_df[ontology], f"GO_table_{category}_{ontology}"
                )
            else:
                html += f"""
<h4>{category.title()} - {ontology}</h4><p>For {ontology} case, we found 0
enriched go terms. </p><br>"""

            filenames = []
            for ontology in self.ontologies:
                if "PROTEIN" in ontology or "PATHWAY" in ontology:
                    continue
                filename = f"{config.output_dir}/Chart_{category}_{ontology}.png"
                if ontology in _temp_df and len(_temp_df[ontology]):
                    logger.info(
                        f"Saving chart for case {ontology} ({category}) in {filename}"
                    )

                    self.pe.save_chart(_temp_df[ontology].iloc[0 : self.nmax], filename)
                    filenames.append(f"Chart_{category}_{ontology}.png")
            foto = self.add_fotorama(filenames, width=1000)
            html += f"<h4>Charts {category} -- </h4> {foto}"

        return html

1	#
2	# This file is part of Sequana software
3	#
4	# Copyright (c) 2020 - Sequana Development Team
5	#
6	# Distributed under the terms of the 3-clause BSD license.
7	# The full license is in the LICENSE file, distributed with this software.
8	#
9	# website: https://github.com/sequana/sequana
10	# documentation: http://sequana.readthedocs.io
11	#
12	##############################################################################
13	"""Module to write enrichment report"""	2✔
14	import os	2✔
15	import sys	2✔
16
17	from sequana.lazy import pandas as pd	2✔
18	from sequana.lazy import pylab	2✔
19
20	from sequana.modules_report.base_module import SequanaBaseModule	2✔
21	from sequana.utils.datatables_js import DataTable	2✔
22	from sequana.enrichment.panther import PantherEnrichment	2✔
23	from sequana.utils import config	2✔
24
25
26	from easydev import Progress	2✔
27
28	import colorlog	2✔
29
30	logger = colorlog.getLogger(__name__)	2✔
31
32
33	class ModulePantherEnrichment(SequanaBaseModule):	2✔
34	"""Write HTML report of variant calling. This class takes a csv file
35	generated by sequana_variant_filter.
36	"""
37
38	def __init__(	2✔
39	self,
40	gene_lists,
41	taxon,
42	enrichment_params={
43	"padj": 0.05,
44	"log2_fc": 3,
45	"max_entries": 3000,
46	"nmax": 50,
47	"mapper": None,
48	"plot_compute_levels": False,
49	"plot_logx": True,
50	},
51	command="",
52	ontologies=["MF", "BP", "CC"],
53	):
54	""".. rubric:: constructor"""
55	super().__init__()	×
56	self.title = "Panther Enrichment"	×
57
58	self.command = command	×
59	self.gene_lists = gene_lists	×
60	self.enrichment_params = enrichment_params	×
61	self.nmax = enrichment_params.get("nmax", 50)	×
62
63	# compute the enrichment here once for all, This may take time
64	from sequana import logger	×
65
66	logger.setLevel("INFO")	×
67
68	logger.info(" === Module PanthEnrichment. === ")	×
69	self.pe = PantherEnrichment(	×
70	self.gene_lists,
71	taxon,
72	max_entries=self.enrichment_params["max_entries"],
73	log2_fc_threshold=self.enrichment_params["log2_fc"],
74	)
75
76	self.ontologies = ontologies	×
77
78	# Compute the enrichment
79	self.pe.compute_enrichment(ontologies=self.ontologies)	×
80	self.df_stats = self.pe.get_mapping_stats()[	×
81	["category", "mapped_percentage", "total"]
82	]
83
84	self.create_report_content()	×
85	self.create_html("enrichment.html")	×
86
87	def create_report_content(self):	2✔
88	self.sections = list()	×
89	self.summary()	×
90	self.add_go()	×
91	self.sections.append(	×
92	{"name": "5 - Info", "anchor": "command", "content": self.command}
93	)
94
95	def summary(self):	2✔
96	"""Add information."""
97
98	total_up = len(self.gene_lists["up"])	×
99	total_down = len(self.gene_lists["down"])	×
100	total = total_up + total_down	×
101	log2fc = self.enrichment_params["log2_fc"]	×
102
103	# create html table for taxon information
104	_taxon_id = self.pe.taxon_info["taxon_id"]	×
105	_taxon_name = self.pe.taxon_info["long_name"]	×
106
107	# table of undertermined IDs
108	df_stats = self.df_stats.drop_duplicates()	×
109	datatable = DataTable(df_stats, "unmapped")	×
110	datatable.datatable.datatable_options = {	×
111	"scrollX": "true",
112	"pageLength": 10,
113	"scrollCollapse": "true",
114	"dom": "Bfrtip",
115	"buttons": ["copy", "csv"],
116	}
117	js = datatable.create_javascript_function()	×
118	html_table = datatable.create_datatable(float_format="%E")	×
119
120	self.sections.append(	×
121	{
122	"name": "1 - Summary",
123	"anchor": "filters_option",
124	"content": f"""
125
126	<p>In the following sections, you will find the GO
127	terms enrichment. The input data for those analyis is the output of the RNADiff
128	analysis where adjusted p-values above 0.05 are excluded. Moreover, we removed
129	candidates with log2 fold change below {log2fc}. Using these filters, the list of
130	differentially expressed genes is made of {total_up} up and {total_down} down genes (total {total})</p>
131	<p> In the following plots you can find the first GO terms that are enriched, keeping a
132	maximum of {self.nmax} identifiers. </p>
133
134	<p>The taxon used is {_taxon_name} (ID {_taxon_id}).<br>
135
136	<p> Check the following table for the percentage of mapped genes on Panther
137	DB</p>
138
139	{js}{html_table}
140	""",
141	}
142	)
143
144	def add_go(self):	2✔
145	# somehow, logger used here and in add_kegg must be global. If you call
146	# add_go and then add_kegg, logger becomes an unbound local variable.
147	# https://stackoverflow.com/questions/10851906/python-3-unboundlocalerror-local-variable-referenced-before-assignment
148	# global logger
149	level = logger.level	×
150	logger.setLevel(level)	×
151
152	html_intro = """	×
153	<p>Here below is a set of plots showing the enriched GO terms using the down
154	regulated genes only, and then the up-regulated genes only. When possible a
155	graph of the found GO terms is provided. MF stands for molecular
156	function, CC for cellular components and BP for biological process.</p>
157	</div>
158	""".format(
159	self.pe.summary["padj_threshold"],
160	self.pe.summary["fold_change_range"][0],
161	self.pe.summary["fold_change_range"][1],
162	self.pe.summary["DGE_after_filtering"]["down"],
163	self.pe.summary["DGE_after_filtering"]["up"],
164	self.pe.summary["DGE_after_filtering"]["all"],
165	)
166
167	html = self._get_enrichment("down")	×
168	self.sections.append(	×
169	{
170	"name": "2 - Enriched GO terms (Down cases)",
171	"anchor": "go_down",
172	"content": html,
173	}
174	)
175
176	return	×
177	html = self._get_enrichment("up")
178	self.sections.append(
179	{
180	"name": "3 - Enriched GO terms (Up cases)",
181	"anchor": "go_up",
182	"content": html,
183	}
184	)
185
186	html = self._get_enrichment("all")
187	self.sections.append(	×
188	{
189	"name": "4 - Enriched GO terms (All cases)",
190	"anchor": "go_all",
191	"content": html,
192	}
193	)
194
195	# a utility function to create the proper html table
196	def get_html_table(self, this_df, identifier):	2✔
197	df = this_df.copy()	×
198
199	# depending on the Panther category, we add the link to the identifier
200	links = []	×
201	for x in df["id"]:	×
202	if x.startswith("PC"):	×
203	links.append(	×
204	f"http://www.pantherdb.org/panther/category.do?categoryAcc={x}"
205	)
206	elif x.startswith("R-"):	×
207	links.append(f"https://reactome.org/PathwayBrowser/#/{x}")	×
208	else:
209	links.append(f"https://www.ebi.ac.uk/QuickGO/term/{x}")	×
210	df["links"] = links	×
211
212	# remove non-informative or redundant fields
213	for x in ["term", "fdr2", "abs_log2_fold_enrichment", "pct_diff_expr"]:	×
214	try:	×
215	del df[x]	×
216	except KeyError:	×
217	pass	×
218
219	first_col = df.pop("id")	×
220	df.insert(0, "id", first_col)	×
221	df = df.sort_values(by="fold_enrichment", ascending=False)	×
222
223	datatable = DataTable(pd.DataFrame(df), identifier)	×
224	datatable.datatable.set_links_to_column("links", "id")	×
225	datatable.datatable.datatable_options = {	×
226	"scrollX": "true",
227	"pageLength": 10,
228	"scrollCollapse": "true",
229	"dom": "Bfrtip",
230	"buttons": ["copy", "csv"],
231	}
232	js = datatable.create_javascript_function()	×
233	html_table = datatable.create_datatable(float_format="%E")	×
234	return js + html_table	×
235
236	def _get_enrichment(self, category):	2✔
237	# category is in down/up/all
238
239	style = "width:95%"	×
240	_temp_df = {}	×
241	_minus = {}	×
242	_plus = {}	×
243
244	html = ""	×
245
246	for ontology in self.ontologies:	×
247
248	# get dat without plotting to store the entire set of GO terms
249	df, subdf = self.pe._get_plot_go_terms_data(	×
250	category, ontologies=ontology, compute_levels=False
251	)
252
253	# now plotting but showing only some restricted GO terms
254	def plot_go_terms(filename, category, ontology):	×
255	df = self.pe.plot_go_terms(	×
256	category,
257	ontologies=ontology,
258	compute_levels=self.enrichment_params["plot_compute_levels"],
259	log=self.enrichment_params["plot_logx"],
260	max_features=self.nmax,
261	)
262	_temp_df[ontology] = df.copy()	×
263	_plus[ontology] = sum(df.plus_minus == "+")	×
264	_minus[ontology] = sum(df.plus_minus == "-")	×
265	pylab.savefig(f"{config.output_dir}/Panther_{category}_{ontology}.png")	×
266	pylab.savefig(filename)	×
267	pylab.close()	×
268
269	if len(df):	×
270	image = self.create_embedded_png(	×
271	plot_go_terms,
272	"filename",
273	style=style,
274	category=category,
275	ontology=ontology,
276	)
277
278	html += f"""	×
279	<h4>{category.title()} - {ontology}</h4>
280	<p>For {ontology}, we found {_plus[ontology]+_minus[ontology]} go terms.
281	Showing {self.nmax} here below (at most). The full list is downlodable from the CSV
282	file hereafter.</p> {image} <br>"""
283	html += self.get_html_table(	×
284	_temp_df[ontology], f"GO_table_{category}_{ontology}"
285	)
286	else:
287	html += f"""	×
288	<h4>{category.title()} - {ontology}</h4><p>For {ontology} case, we found 0
289	enriched go terms. </p><br>"""
290
291	filenames = []	×
292	for ontology in self.ontologies:	×
293	if "PROTEIN" in ontology or "PATHWAY" in ontology:	×
294	continue	×
295	filename = f"{config.output_dir}/Chart_{category}_{ontology}.png"	×
296	if ontology in _temp_df and len(_temp_df[ontology]):	×
297	logger.info(	×
298	f"Saving chart for case {ontology} ({category}) in {filename}"
299	)
300
301	self.pe.save_chart(_temp_df[ontology].iloc[0 : self.nmax], filename)	×
302	filenames.append(f"Chart_{category}_{ontology}.png")	×
303	foto = self.add_fotorama(filenames, width=1000)	×
304	html += f"<h4>Charts {category} -- </h4> {foto}"	×
305
306	return html	×

sequana / sequana / 1261710972

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