2578756179

Build Type

github

Committed by GitHub

Commit Message

<a href="https://github.com/sequana/sequana/commit/<a class=hub.com/sequana/sequana/commit/b4bbe2ce9558a99ceeecdc3677a2c492b31a3c9b">b4bbe2ce9<a href="https://github.com/sequana/sequana/commit/b4bbe2ce9558a99ceeecdc3677a2c492b31a3c9b">">Merge </a><a class="double-link" href="https://github.com/sequana/sequana/commit/<a class="double-link" href="https://github.com/sequana/sequana/commit/d4eb9f533b8baffc919554333b0c5d42b34f6987">d4eb9f533</a>">d4eb9f533</a><a href="https://github.com/sequana/sequana/commit/b4bbe2ce9558a99ceeecdc3677a2c492b31a3c9b"> into 01b52271f">01b52271f</a>

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

Run Details

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

11936 of 17082 relevant lines covered (69.87%)

2.1 hits per line

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

86.18

/sequana/enrichment/panther.py

#
#  This file is part of Sequana software
#
#  Copyright (c) 2016-2021 - 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 colorlog
from sequana.enrichment.ontology import Ontology
from sequana.enrichment.plot_go_terms import PlotGOTerms
from sequana.enrichment.quickgo import QuickGOGraph
from sequana.lazy import numpy as np
from sequana.lazy import pandas as pd
from sequana.lazy import pylab
from tqdm import tqdm

logger = colorlog.getLogger(__name__)


__all__ = ["PantherEnrichment"]


class PantherEnrichment(Ontology, PlotGOTerms):
    """

    By default, we keep only the genes with a adjusted pvalue <= 0.05. The
    fold_change threshold is on a log2 scale and set to 0 (meaning no
    filtering). Only one value is requested and
    used to filter out positive and negative fold change (on the log2
    scale). In other word, a log2 fold change threshold of 2 means that we
    filter out values between -2 and 2.

    If you prefer to work in natural scale, just set the parameter
    fc_threshold, which overwrite the log2_fc_threshold parameter.

    ::

        pe = PantherEnrichment("input_file.tsv", taxon=10090, log2_fc_threshold=1)
        # compute enrichment for genes down and up
        pe.compute_enrichment()

        # Results for up case is stored in pe.enrichment
        # then, we plot the most mportat go terms
        df_up = pe.plot_go_terms("up")

        df = pe.plot_go_terms("up", pe.MF)
        pe.save_chart(df, "chart_MF_up.png")

        # all 3 main ontology
        df = pe.plot_go_terms("up")
        pe.save_chart(df, "chart_up.png")

    e.stats contains some statistics. One important is the list of unmapped
    genes. The results from the GO enrichment are stored in the attributes
    enrichment. There, we have again adjusted p-value and a fold enrichment,
    which can in turn be filtered or not.

    You can retrieve the cleaned data using the get_data method.

    You can also plot the GO terms that are significantly enriched using::

        e.plot_go_terms(['MF', 'CC', 'BP])

    This function returns the dataframe used during the plotting.

    If you want to look at the up regulated genes only::

        e.compute_enrichment(pe.mygenes["up"], 83333)
        df = e.plot_go_terms(['MF', 'CC', 'BP'],
                log=False, include_negative_enrichment=False,
                fontsize=8, sort_by='fold_enrichment',
                show_pvalues=True, fdr_threshold=0.05)


    The number of genes is limited to about 3100 depending (don't ask me
    why, this seem to be a hard-coded limitation on PantherDB website).
    In such case, you should add a filter e.g on padj or fold change


    """

    def __init__(
        self,
        gene_lists,
        taxon,
        requests_per_sec=10,
        padj_threshold=0.05,
        log2_fc_threshold=0,
        fc_threshold=None,
        enrichment_fdr=0.05,
        annot_col="Name",
    ):
        """


        rnadiff if provided, superseeds the input filename. This is useful for
        debugging
        """
        Ontology.__init__(self)
        PlotGOTerms.__init__(self)

        self.gene_lists = gene_lists
        self.enrichment_fdr = enrichment_fdr

        # users can set the fold change threshold in the log2 scale or normal
        # scale.
        assert log2_fc_threshold >= 0, "log2 fc_threshold must be >=0"
        if fc_threshold is not None:
            log2_fc_threshold = pylab.log2(fc_threshold)

        from bioservices import panther, quickgo

        self.quick_go_graph = QuickGOGraph()

        self.panther = panther.Panther(cache=True)
        self.valid_taxons = [x["taxon_id"] for x in self.panther.get_supported_genomes()]
        self.summary = {}

        self._taxon = None
        self.taxon = taxon

        self.quickgo = quickgo.QuickGO(cache=True)
        self.quickgo.requests_per_sec = requests_per_sec
        self.quickgo.settings.TIMEOUT = 120

        self._ancestors = {
            "MF": "GO:0003674",
            "CC": "GO:0005575",
            "BP": "GO:0008150",
            "SLIM_MF": "GO:0003674",
            "SLIM_CC": "GO:0005575",
            "SLIM_BP": "GO:0008150",
        }
        self.ontologies.extend(
            [
                "ANNOT_TYPE_ID_PANTHER_GO_SLIM_MF",
                "ANNOT_TYPE_ID_PANTHER_GO_SLIM_BP",
                "ANNOT_TYPE_ID_PANTHER_GO_SLIM_CC",
                "ANNOT_TYPE_ID_PANTHER_PC",
                "ANNOT_TYPE_ID_PANTHER_PATHWAY",
                "ANNOT_TYPE_ID_REACTOME_PATHWAY",
            ]
        )

        self.ontology_aliases.extend(
            [
                "SLIM_MF",
                "SLIM_BP",
                "SLIM_CC",
                "PROTEIN",
                "PANTHER_PATHWAY",
                "REACTOME_PATHWAY",
            ]
        )

        # panther accepts onyl ~2-3000 genes at max. Let us restrict the analysis
        # to the first 2000 genes based on their log2 fold change 2000 + and
        # 2000 negatives

        msg = "Ignoring DEGs with adjusted p-value > {} and fold change in [{}, {}]".format(
            padj_threshold, 1 / (2**log2_fc_threshold), 2**log2_fc_threshold
        )
        logger.info(msg)

        # used in report module
        self.summary["fold_change_range"] = [
            1 / (2**log2_fc_threshold),
            2**log2_fc_threshold,
        ]
        self.summary["padj_threshold"] = padj_threshold

        fc_threshold = log2_fc_threshold

        for x in sorted(gene_lists.keys()):

            N = len(gene_lists[x])
            logger.info(f"Starting with {N} genes from category '{x}'")

        self.summary["DGE_after_filtering"] = {k: len(v) for k, v in gene_lists.items()}

        self.enrichment = {}
        self.stats = {}
        self.obsolets = []

    def _set_taxon(self, taxon):
        if taxon not in self.valid_taxons:
            raise ValueError(f"taxon {taxon} not in pantherDB. please use one of {self.valid_taxons}")
        self.taxon_info = [x for x in self.panther.get_supported_genomes() if x["taxon_id"] == taxon]
        self.taxon_info = self.taxon_info[0]
        self._taxon_id = taxon

    def _get_taxon(self):
        return self._taxon_id

    taxon = property(_get_taxon, _set_taxon)

    def compute_enrichment(
        self,
        taxid=None,
        ontologies=None,
        enrichment_test="FISHER",
        correction="FDR",
    ):
        """
        :param enrichment_test: Fisher or Binomial
        :param correction: FDR or Bonferonni

        The field **number_in_reference** indicates from the reference, the number
        of genes that have a given ontolgy term. For instance, 998 genes have
        the term. This is stored in **number_in_reference**. If the reference
        contains 4391 genes, and you provided 49
        genes , the **expected** number of genes that have this ontology term is
        49*998/4391 that is 11.1369, which is stored in **"expected**.

        Now, if you actually find that 14 out of 49
        genes have the term, you need to compare the numbers 11.1369 and 14. Are
        they really different ? The ratio 14 / 11.1369 is stored in
        **fold_enrichment**. The pvalue and FDR are stored as well.

        Some genes may be missing If you provide 50 genes, you may end up with
        only 45 being mapped onto panther db database. This may explain some
        differences with the expected value.

        Fold enrichment is the number_in_list / expected ratio. Another close metric is the
        fractional difference: (observed - expected) / expected. This metric is
        slighlty less than the fold enrichment

        To get the number f genes from the database, simply take one output, and
        compute number_in_reference * number_in_list / expected

        The fold enrichment is also called odd-ratio.
        """

        for category, gene_list in self.gene_lists.items():
            logger.info(f"Computing enrichment for category {category}")
            self.enrichment[category], self.stats[category] = self._compute_enrichment(
                gene_list,
                taxid=taxid,
                ontologies=ontologies,
                enrichment_test=enrichment_test,
                correction=correction,
            )

    def get_mapping_stats(self):

        results = []
        for k in self.enrichment.keys():
            ontologies = self.enrichment[k].keys()
            for o in ontologies:
                M = self.enrichment[k][o]["input_list"]["mapped_count"]
                U = self.enrichment[k][o]["input_list"]["unmapped_count"]
                results.append([k, o, M, U])
        df = pd.DataFrame(results)
        df[4] = 100 * df[2] / (df[3] + df[2])
        df[5] = df[2] + df[3]
        df.columns = [
            "category",
            "ontology",
            "mapped",
            "unmapped",
            "mapped_percentage",
            "total",
        ]
        return df

    def _compute_enrichment(
        self,
        mygenes,
        taxid,
        ontologies=None,
        enrichment_test="FISHER",
        correction="FDR",
    ):

        # taxid=83333 # ecoli
        if taxid is None:
            taxid = self.taxon

        if isinstance(mygenes, list):
            mygenes = ",".join(mygenes)

        if mygenes.count(",") > 2500:  # pragma: no cover
            logger.warning("Please reduce the list input genes. may fail on pantherb otherwise")
        if len(mygenes) <= 2:  # pragma: no cover
            logger.error(f"Less than 2 genes are found for in the gene set: {mygenes}. No enrichment will be computed")
            return None, None

        if ontologies is None:
            ontologies = self.ontology_aliases
        else:
            for x in ontologies:
                assert x in self.ontology_aliases, f"{x} not valid ontology"

        # for each ontology, we will store one key/value item
        enrichment = {}

        def get_panther_ont(x):
            index = self.ontology_aliases.index(x)
            return self.ontologies[index]

        unclassified = {}
        for ontology in ontologies:
            unclassified[ontology] = 0
            logger.info(f" - Computing enrichment for {ontology}")

            try:
                del results
            except NameError:
                pass
            results = self.panther.get_enrichment(
                mygenes,
                taxid,
                get_panther_ont(ontology),
                enrichment_test=enrichment_test,
                correction=correction,
            )
            count = 0
            while count < 2 and results == 404:  # pragma: no cover
                logger.warning("Panther request failed Trying again")
                results = self.panther.get_enrichment(
                    mygenes,
                    taxid,
                    get_panther_ont(ontology),
                    enrichment_test=enrichment_test,
                    correction=correction,
                )
                count += 1

            if results == 404:  # pragma: no cover
                logger.warning("Invalid output from pantherdb (too many genes ?). skipping {}".format(ontology))
                enrichment[ontology] = None
                continue

            ##if isinstance(results["result"], dict):  # pragma: no cover
            #    results["result"] = [results["result"]]

            # extract the ID and label of the go term and save
            # as primary keys
            for i, k in enumerate(results["result"]):
                if results["result"][i]["term"]["label"] == "UNCLASSIFIED":
                    unclassified[ontology] += 1
                    # results["result"][i]['label'] = None
                else:
                    results["result"][i]["id"] = k["term"]["id"]
                    results["result"][i]["label"] = k["term"]["label"]

            # convert to dataframe and remove unclassified labels
            df = pd.DataFrame(results["result"])
            df = df[df["label"].isnull() == False]

            # compute pvalues adjusted
            pvalues = df["pValue"].copy()

            import statsmodels
            import statsmodels.stats.multitest

            if correction == "FDR":
                fdr = statsmodels.stats.multitest.multipletests(pvalues, method="fdr_bh")[1]
            elif correction.lower() == "bonferroni":
                fdr = statsmodels.stats.multitest.multipletests(pvalues, method="bonferroni")[1]

            df["fdr2"] = fdr
            if enrichment_test.lower() == "binomial":
                df["fdr"] = fdr

            # store all results
            enrichment[ontology] = results
            enrichment[ontology]["result"] = df

        stats = dict([(k, len(v["result"])) for k, v in enrichment.items()])
        stats["input_genes"] = len(mygenes.split(","))

        try:
            unmapped = enrichment[ontologies[0]]["input_list"]["unmapped_id"]
            stats["unmapped_genes"] = unmapped
            stats["N_unmapped_genes"] = len(unmapped)
        except Exception:
            stats["unmapped_genes"] = []
            stats["N_unmapped_genes"] = 0
        stats["unclassified_GO"] = unclassified

        # Here, looking at the FDr, it appears that when using bonferroni,
        # all FDR are set to zeros. Moreover, when using Fisher tests and
        # FDR (supposibly a FDR_BH, the results are noisy as compare to a
        # test from statsmodels. Moreover, when using binomial test, the FDR
        # is not computed... So, we will recompute the FDR ourself
        return enrichment, stats

    def get_data(self, category, ontologies, include_negative_enrichment=True, fdr=0.05):
        """

        From all input GO term that have been found and stored in
        enrichment[ONTOLOGY]['result'], we keep those with fdr<0.05. We also
        exclude UNCLASSIFIED entries. The final dataframe is returned

        ::

            pe.get_data(""up", MF")

        """

        df = self._get_data(category, ontologies)

        # some extra information for convenience
        df["pct_diff_expr"] = df["number_in_list"] * 100 / df["number_in_reference"]

        # could happen that all fold enrichment are set to 'NaN'
        df = df[df['fold_enrichment'] != 'NaN']
        try:
            df["log2_fold_enrichment"] = pylab.log2(df["fold_enrichment"])
        except:
            import pdb
            pdb.set_trace()
        df["abs_log2_fold_enrichment"] = abs(pylab.log2(df["fold_enrichment"]))
        df["expected"] = [int(x) for x in df.expected]

        # Some user may want to include GO terms with fold enrichment
        # significanyly below 1 or not.
        if include_negative_enrichment is False:
            df = df.query("fold_enrichment>=1").copy()
            logger.info("Found {} GO terms after keeping only positive enrichment".format(len(df)))

        # filter out FDR>0.05
        df = df.query("fdr<=@fdr").copy()
        logger.info("Found {} GO terms after keeping only FDR<{}".format(len(df), fdr))

        return df

    def get_functional_classification(self, mygenes, taxon):  # pragma: no cover ; too slow
        """Mapping information from pantherDB for the lisf of genes

        We also store uniprot persistent id

        """
        logger.warning("Very slow. Please wait")
        if isinstance(mygenes, list):
            mygenes = ",".join(mygenes)

        res = self.panther.get_mapping(mygenes, taxon)
        res = res["mapped"]

        for i, item in tqdm(enumerate(res)):
            accession = item["accession"]
            res[i]["persistent_id"] = self._get_name_given_accession(accession)
        return res

    def _get_name_given_accession(self, accession):  # pragma: no cover
        from bioservices import UniProt

        self.uniprot = UniProt(cache=True)
        self.uniprot.requests_per_sec = 10

        acc = [x for x in accession.split("|") if x.startswith("UniProtKB")]
        acc = acc[0].split("=")[1]
        res = self.uniprot.get_df(acc, limit=1)
        name = res["Gene names  (primary )"][0]
        return name

    def plot_piechart(self, df):
        # Here we show the GO terms that have number in list > 0
        # Note, that this is dangerous to look only at this picture without
        # the reference plot, which data is not available thourg the pathner API
        labels = []
        for this in df.query("number_in_list!=0").label.values:
            if len(this) > 50:
                labels.append(this)
            else:
                labels.append(this[0:50] + "...")
        pylab.pie(df.query("number_in_list!=0").number_in_list, labels=labels)
        pylab.tight_layout()

    def save_chart(self, df, filename="chart.png"):
        """

        pe = PantherEnrichment("B4052-V1.T1vsT0.complete.xls", fc_threshold=5,
            padj_threshold=0.05)
        df = pe.plot_go_terms("down", log=True, compute_levels=False)
        pe.save_chart(df, "chart.png")

        """
        self.quick_go_graph.save_chart(df, filename)

    def _get_graph(self, df, ontologies):
        return self.quick_go_graph._get_graph(df, ontologies=ontologies)

1	#
2	# This file is part of Sequana software
3	#
4	# Copyright (c) 2016-2021 - 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
14	import colorlog	3✔
15	from sequana.enrichment.ontology import Ontology	3✔
16	from sequana.enrichment.plot_go_terms import PlotGOTerms	3✔
17	from sequana.enrichment.quickgo import QuickGOGraph	3✔
18	from sequana.lazy import numpy as np	3✔
19	from sequana.lazy import pandas as pd	3✔
20	from sequana.lazy import pylab	3✔
21	from tqdm import tqdm	3✔
22
23	logger = colorlog.getLogger(__name__)	3✔
24
25
26	__all__ = ["PantherEnrichment"]	3✔
27
28
29	class PantherEnrichment(Ontology, PlotGOTerms):	3✔
30	"""
31
32	By default, we keep only the genes with a adjusted pvalue <= 0.05. The
33	fold_change threshold is on a log2 scale and set to 0 (meaning no
34	filtering). Only one value is requested and
35	used to filter out positive and negative fold change (on the log2
36	scale). In other word, a log2 fold change threshold of 2 means that we
37	filter out values between -2 and 2.
38
39	If you prefer to work in natural scale, just set the parameter
40	fc_threshold, which overwrite the log2_fc_threshold parameter.
41
42	::
43
44	pe = PantherEnrichment("input_file.tsv", taxon=10090, log2_fc_threshold=1)
45	# compute enrichment for genes down and up
46	pe.compute_enrichment()
47
48	# Results for up case is stored in pe.enrichment
49	# then, we plot the most mportat go terms
50	df_up = pe.plot_go_terms("up")
51
52	df = pe.plot_go_terms("up", pe.MF)
53	pe.save_chart(df, "chart_MF_up.png")
54
55	# all 3 main ontology
56	df = pe.plot_go_terms("up")
57	pe.save_chart(df, "chart_up.png")
58
59	e.stats contains some statistics. One important is the list of unmapped
60	genes. The results from the GO enrichment are stored in the attributes
61	enrichment. There, we have again adjusted p-value and a fold enrichment,
62	which can in turn be filtered or not.
63
64	You can retrieve the cleaned data using the get_data method.
65
66	You can also plot the GO terms that are significantly enriched using::
67
68	e.plot_go_terms(['MF', 'CC', 'BP])
69
70	This function returns the dataframe used during the plotting.
71
72	If you want to look at the up regulated genes only::
73
74	e.compute_enrichment(pe.mygenes["up"], 83333)
75	df = e.plot_go_terms(['MF', 'CC', 'BP'],
76	log=False, include_negative_enrichment=False,
77	fontsize=8, sort_by='fold_enrichment',
78	show_pvalues=True, fdr_threshold=0.05)
79
80
81	The number of genes is limited to about 3100 depending (don't ask me
82	why, this seem to be a hard-coded limitation on PantherDB website).
83	In such case, you should add a filter e.g on padj or fold change
84
85
86	"""
87
88	def __init__(	3✔
89	self,
90	gene_lists,
91	taxon,
92	requests_per_sec=10,
93	padj_threshold=0.05,
94	log2_fc_threshold=0,
95	fc_threshold=None,
96	enrichment_fdr=0.05,
97	annot_col="Name",
98	):
99	"""
100
101
102	rnadiff if provided, superseeds the input filename. This is useful for
103	debugging
104	"""
105	Ontology.__init__(self)	3✔
106	PlotGOTerms.__init__(self)	3✔
107
108	self.gene_lists = gene_lists	3✔
109	self.enrichment_fdr = enrichment_fdr	3✔
110
111	# users can set the fold change threshold in the log2 scale or normal
112	# scale.
113	assert log2_fc_threshold >= 0, "log2 fc_threshold must be >=0"	3✔
114	if fc_threshold is not None:	3✔
115	log2_fc_threshold = pylab.log2(fc_threshold)	3✔
116
117	from bioservices import panther, quickgo	3✔
118
119	self.quick_go_graph = QuickGOGraph()	3✔
120
121	self.panther = panther.Panther(cache=True)	3✔
122	self.valid_taxons = [x["taxon_id"] for x in self.panther.get_supported_genomes()]	3✔
123	self.summary = {}	3✔
124
125	self._taxon = None	3✔
126	self.taxon = taxon	3✔
127
128	self.quickgo = quickgo.QuickGO(cache=True)	3✔
129	self.quickgo.requests_per_sec = requests_per_sec	3✔
130	self.quickgo.settings.TIMEOUT = 120	3✔
131
132	self._ancestors = {	3✔
133	"MF": "GO:0003674",
134	"CC": "GO:0005575",
135	"BP": "GO:0008150",
136	"SLIM_MF": "GO:0003674",
137	"SLIM_CC": "GO:0005575",
138	"SLIM_BP": "GO:0008150",
139	}
140	self.ontologies.extend(	3✔
141	[
142	"ANNOT_TYPE_ID_PANTHER_GO_SLIM_MF",
143	"ANNOT_TYPE_ID_PANTHER_GO_SLIM_BP",
144	"ANNOT_TYPE_ID_PANTHER_GO_SLIM_CC",
145	"ANNOT_TYPE_ID_PANTHER_PC",
146	"ANNOT_TYPE_ID_PANTHER_PATHWAY",
147	"ANNOT_TYPE_ID_REACTOME_PATHWAY",
148	]
149	)
150
151	self.ontology_aliases.extend(	3✔
152	[
153	"SLIM_MF",
154	"SLIM_BP",
155	"SLIM_CC",
156	"PROTEIN",
157	"PANTHER_PATHWAY",
158	"REACTOME_PATHWAY",
159	]
160	)
161
162	# panther accepts onyl ~2-3000 genes at max. Let us restrict the analysis
163	# to the first 2000 genes based on their log2 fold change 2000 + and
164	# 2000 negatives
165
166	msg = "Ignoring DEGs with adjusted p-value > {} and fold change in [{}, {}]".format(	3✔
167	padj_threshold, 1 / (2log2_fc_threshold), 2log2_fc_threshold
168	)
169	logger.info(msg)	3✔
170
171	# used in report module
172	self.summary["fold_change_range"] = [	3✔
173	1 / (2**log2_fc_threshold),
174	2**log2_fc_threshold,
175	]
176	self.summary["padj_threshold"] = padj_threshold	3✔
177
178	fc_threshold = log2_fc_threshold	3✔
179
180	for x in sorted(gene_lists.keys()):	3✔
181
182	N = len(gene_lists[x])	3✔
183	logger.info(f"Starting with {N} genes from category '{x}'")	3✔
184
185	self.summary["DGE_after_filtering"] = {k: len(v) for k, v in gene_lists.items()}	3✔
186
187	self.enrichment = {}	3✔
188	self.stats = {}	3✔
189	self.obsolets = []	3✔
190
191	def _set_taxon(self, taxon):	3✔
192	if taxon not in self.valid_taxons:	3✔
193	raise ValueError(f"taxon {taxon} not in pantherDB. please use one of {self.valid_taxons}")	×
194	self.taxon_info = [x for x in self.panther.get_supported_genomes() if x["taxon_id"] == taxon]	3✔
195	self.taxon_info = self.taxon_info[0]	3✔
196	self._taxon_id = taxon	3✔
197
198	def _get_taxon(self):	3✔
199	return self._taxon_id	3✔
200
201	taxon = property(_get_taxon, _set_taxon)	3✔
202
203	def compute_enrichment(	3✔
204	self,
205	taxid=None,
206	ontologies=None,
207	enrichment_test="FISHER",
208	correction="FDR",
209	):
210	"""
211	:param enrichment_test: Fisher or Binomial
212	:param correction: FDR or Bonferonni
213
214	The field number_in_reference indicates from the reference, the number
215	of genes that have a given ontolgy term. For instance, 998 genes have
216	the term. This is stored in number_in_reference. If the reference
217	contains 4391 genes, and you provided 49
218	genes , the expected number of genes that have this ontology term is
219	49998/4391 that is 11.1369, which is stored in "expected*.
220
221	Now, if you actually find that 14 out of 49
222	genes have the term, you need to compare the numbers 11.1369 and 14. Are
223	they really different ? The ratio 14 / 11.1369 is stored in
224	fold_enrichment. The pvalue and FDR are stored as well.
225
226	Some genes may be missing If you provide 50 genes, you may end up with
227	only 45 being mapped onto panther db database. This may explain some
228	differences with the expected value.
229
230	Fold enrichment is the number_in_list / expected ratio. Another close metric is the
231	fractional difference: (observed - expected) / expected. This metric is
232	slighlty less than the fold enrichment
233
234	To get the number f genes from the database, simply take one output, and
235	compute number_in_reference * number_in_list / expected
236
237	The fold enrichment is also called odd-ratio.
238	"""
239
240	for category, gene_list in self.gene_lists.items():	3✔
241	logger.info(f"Computing enrichment for category {category}")	3✔
242	self.enrichment[category], self.stats[category] = self._compute_enrichment(	3✔
243	gene_list,
244	taxid=taxid,
245	ontologies=ontologies,
246	enrichment_test=enrichment_test,
247	correction=correction,
248	)
249
250	def get_mapping_stats(self):	3✔
251
252	results = []	×
253	for k in self.enrichment.keys():	×
254	ontologies = self.enrichment[k].keys()	×
255	for o in ontologies:	×
256	M = self.enrichment[k][o]["input_list"]["mapped_count"]	×
257	U = self.enrichment[k][o]["input_list"]["unmapped_count"]	×
258	results.append([k, o, M, U])	×
259	df = pd.DataFrame(results)	×
260	df[4] = 100 * df[2] / (df[3] + df[2])	×
261	df[5] = df[2] + df[3]	×
262	df.columns = [	×
263	"category",
264	"ontology",
265	"mapped",
266	"unmapped",
267	"mapped_percentage",
268	"total",
269	]
270	return df	×
271
272	def _compute_enrichment(	3✔
273	self,
274	mygenes,
275	taxid,
276	ontologies=None,
277	enrichment_test="FISHER",
278	correction="FDR",
279	):
280
281	# taxid=83333 # ecoli
282	if taxid is None:	3✔
283	taxid = self.taxon	3✔
284
285	if isinstance(mygenes, list):	3✔
286	mygenes = ",".join(mygenes)	3✔
287
288	if mygenes.count(",") > 2500: # pragma: no cover
289	logger.warning("Please reduce the list input genes. may fail on pantherb otherwise")
290	if len(mygenes) <= 2: # pragma: no cover
291	logger.error(f"Less than 2 genes are found for in the gene set: {mygenes}. No enrichment will be computed")
292	return None, None
293
294	if ontologies is None:	3✔
295	ontologies = self.ontology_aliases	×
296	else:
297	for x in ontologies:	3✔
298	assert x in self.ontology_aliases, f"{x} not valid ontology"	3✔
299
300	# for each ontology, we will store one key/value item
301	enrichment = {}	3✔
302
303	def get_panther_ont(x):	3✔
304	index = self.ontology_aliases.index(x)	3✔
305	return self.ontologies[index]	3✔
306
307	unclassified = {}	3✔
308	for ontology in ontologies:	3✔
309	unclassified[ontology] = 0	3✔
310	logger.info(f" - Computing enrichment for {ontology}")	3✔
311
312	try:	3✔
313	del results	3✔
314	except NameError:	3✔
315	pass	3✔
316	results = self.panther.get_enrichment(	3✔
317	mygenes,
318	taxid,
319	get_panther_ont(ontology),
320	enrichment_test=enrichment_test,
321	correction=correction,
322	)
323	count = 0	3✔
324	while count < 2 and results == 404: # pragma: no cover
325	logger.warning("Panther request failed Trying again")
326	results = self.panther.get_enrichment(
327	mygenes,
328	taxid,
329	get_panther_ont(ontology),
330	enrichment_test=enrichment_test,
331	correction=correction,
332	)
333	count += 1
334
335	if results == 404: # pragma: no cover
336	logger.warning("Invalid output from pantherdb (too many genes ?). skipping {}".format(ontology))
337	enrichment[ontology] = None
338	continue
339
340	##if isinstance(results["result"], dict): # pragma: no cover
341	# results["result"] = [results["result"]]
342
343	# extract the ID and label of the go term and save
344	# as primary keys
345	for i, k in enumerate(results["result"]):	3✔
346	if results["result"][i]["term"]["label"] == "UNCLASSIFIED":	3✔
347	unclassified[ontology] += 1	3✔
348	# results["result"][i]['label'] = None
349	else:
350	results["result"][i]["id"] = k["term"]["id"]	3✔
351	results["result"][i]["label"] = k["term"]["label"]	3✔
352
353	# convert to dataframe and remove unclassified labels
354	df = pd.DataFrame(results["result"])	3✔
355	df = df[df["label"].isnull() == False]	3✔
356
357	# compute pvalues adjusted
358	pvalues = df["pValue"].copy()	3✔
359
360	import statsmodels	3✔
361	import statsmodels.stats.multitest	3✔
362
363	if correction == "FDR":	3✔
364	fdr = statsmodels.stats.multitest.multipletests(pvalues, method="fdr_bh")[1]	3✔
365	elif correction.lower() == "bonferroni":	3✔
366	fdr = statsmodels.stats.multitest.multipletests(pvalues, method="bonferroni")[1]	3✔
367
368	df["fdr2"] = fdr	3✔
369	if enrichment_test.lower() == "binomial":	3✔
370	df["fdr"] = fdr	×
371
372	# store all results
373	enrichment[ontology] = results	3✔
374	enrichment[ontology]["result"] = df	3✔
375
376	stats = dict([(k, len(v["result"])) for k, v in enrichment.items()])	3✔
377	stats["input_genes"] = len(mygenes.split(","))	3✔
378
379	try:	3✔
380	unmapped = enrichment[ontologies[0]]["input_list"]["unmapped_id"]	3✔
381	stats["unmapped_genes"] = unmapped	3✔
382	stats["N_unmapped_genes"] = len(unmapped)	3✔
383	except Exception:	×
384	stats["unmapped_genes"] = []	×
385	stats["N_unmapped_genes"] = 0	×
386	stats["unclassified_GO"] = unclassified	3✔
387
388	# Here, looking at the FDr, it appears that when using bonferroni,
389	# all FDR are set to zeros. Moreover, when using Fisher tests and
390	# FDR (supposibly a FDR_BH, the results are noisy as compare to a
391	# test from statsmodels. Moreover, when using binomial test, the FDR
392	# is not computed... So, we will recompute the FDR ourself
393	return enrichment, stats	3✔
394
395	def get_data(self, category, ontologies, include_negative_enrichment=True, fdr=0.05):	3✔
396	"""
397
398	From all input GO term that have been found and stored in
399	enrichment[ONTOLOGY]['result'], we keep those with fdr<0.05. We also
400	exclude UNCLASSIFIED entries. The final dataframe is returned
401
402	::
403
404	pe.get_data(""up", MF")
405
406	"""
407
408	df = self._get_data(category, ontologies)	3✔
409
410	# some extra information for convenience
411	df["pct_diff_expr"] = df["number_in_list"] * 100 / df["number_in_reference"]	3✔
412
413	# could happen that all fold enrichment are set to 'NaN'
414	df = df[df['fold_enrichment'] != 'NaN']	3✔
415	try:	3✔
416	df["log2_fold_enrichment"] = pylab.log2(df["fold_enrichment"])	3✔
417	except:	×
418	import pdb	×
419	pdb.set_trace()	×
420	df["abs_log2_fold_enrichment"] = abs(pylab.log2(df["fold_enrichment"]))	3✔
421	df["expected"] = [int(x) for x in df.expected]	3✔
422
423	# Some user may want to include GO terms with fold enrichment
424	# significanyly below 1 or not.
425	if include_negative_enrichment is False:	3✔
426	df = df.query("fold_enrichment>=1").copy()	3✔
427	logger.info("Found {} GO terms after keeping only positive enrichment".format(len(df)))	3✔
428
429	# filter out FDR>0.05
430	df = df.query("fdr<=@fdr").copy()	3✔
431	logger.info("Found {} GO terms after keeping only FDR<{}".format(len(df), fdr))	3✔
432
433	return df	3✔
434
435	def get_functional_classification(self, mygenes, taxon): # pragma: no cover ; too slow
436	"""Mapping information from pantherDB for the lisf of genes
437
438	We also store uniprot persistent id
439
440	"""
441	logger.warning("Very slow. Please wait")
442	if isinstance(mygenes, list):
443	mygenes = ",".join(mygenes)
444
445	res = self.panther.get_mapping(mygenes, taxon)
446	res = res["mapped"]
447
448	for i, item in tqdm(enumerate(res)):
449	accession = item["accession"]
450	res[i]["persistent_id"] = self._get_name_given_accession(accession)
451	return res
452
453	def _get_name_given_accession(self, accession): # pragma: no cover
454	from bioservices import UniProt
455
456	self.uniprot = UniProt(cache=True)
457	self.uniprot.requests_per_sec = 10
458
459	acc = [x for x in accession.split("\|") if x.startswith("UniProtKB")]
460	acc = acc[0].split("=")[1]
461	res = self.uniprot.get_df(acc, limit=1)
462	name = res["Gene names (primary )"][0]
463	return name
464
465	def plot_piechart(self, df):	3✔
466	# Here we show the GO terms that have number in list > 0
467	# Note, that this is dangerous to look only at this picture without
468	# the reference plot, which data is not available thourg the pathner API
469	labels = []	3✔
470	for this in df.query("number_in_list!=0").label.values:	3✔
471	if len(this) > 50:	3✔
472	labels.append(this)	3✔
473	else:
474	labels.append(this[0:50] + "...")	3✔
475	pylab.pie(df.query("number_in_list!=0").number_in_list, labels=labels)	3✔
476	pylab.tight_layout()	3✔
477
478	def save_chart(self, df, filename="chart.png"):	3✔
479	"""
480
481	pe = PantherEnrichment("B4052-V1.T1vsT0.complete.xls", fc_threshold=5,
482	padj_threshold=0.05)
483	df = pe.plot_go_terms("down", log=True, compute_levels=False)
484	pe.save_chart(df, "chart.png")
485
486	"""
487	self.quick_go_graph.save_chart(df, filename)	3✔
488
489	def _get_graph(self, df, ontologies):	3✔
490	return self.quick_go_graph._get_graph(df, ontologies=ontologies)	3✔

sequana / sequana / 2578756179

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