================================

Differential Expression Analysis

================================

InMoose offers a Python port of the well-known R Bioconductor packages:

- `DESeq2 package <https://bioconductor.org/packages/release/bioc/html/DESeq2.html>`_ [Love2014]_ in module :doc:`deseq`.

- `edgeR package <https://bioconductor.org/packages/release/bioc/html/edgeR.html>`_ [Chen2016]_ in module :doc:`edgepy`.

- `limma package <https://bioconductor.org/packages/release/bioc/html/limma.html>`_ [Ritchie2015]_. in module :doc:`limma`.

Note that not all features of the R packages are necessarily ported. Extending

the functionality of these modules will be based on user requests, so do not

hesitate to open an issue if your favorite feature is missing.

In addition, InMoose provides a meta-analysis feature to combine the results

from different differential expression analysis tools.

Please refer to [Colange2024]_ for a detailed comparison of InMoose

implementation with the original R implementations.

Differential Expression Meta-Analysis

=====================================

.. currentmodule:: inmoose.diffexp

We illustrate the differential expression meta-analysis capabilities of InMoose

along two approaches:

- the Aggregate Data (AD) approach consists in running classical differential

  expression tools on individual cohorts then combining the results through

  *e.g.* random-effect models.

- the Individual Sample Data (ISD) consists in merging individual cohorts into a

  large meta-cohort, accounting for batch effects to eliminate inter-cohort

  biases, then running a classical differential expression analysis on the

  resulting meta-cohort.

We start by simulating RNA-Seq data, using the :mod:`sim` module of InMoose.

.. repl::

   import numpy as np

   import pandas as pd

   from inmoose.sim import sim_rnaseq

   # number of genes

   N = 10000

   # number of samples

   M = 2000

   assert M % 10 == 0

   P = M // 10  # 10% of M, helper variable

   # 3 batches: 20% 30% 50% of the samples

   batch = (2 * P) * [0] + (3 * P) * [1] + (5 * P) * [2]

   batch = np.array([f"batch{b}" for b in batch])

   batch0 = batch == "batch0"

   batch1 = batch == "batch1"

   batch2 = batch == "batch2"

   # 2 condition groups

   #   - group 1: 50% batch 1, 33% batch 2, 60% batch 3

   #   - group 2: 50% batch 1, 67% batch 2, 40% batch 3

   group = P * [0] + P * [1] + P * [0] + (2 * P) * [1] + (2 * P) * [0] + (3 * P) * [1]

   group = np.array([f"group{g}" for g in group])

   assert len(batch) == M and len(group) == M

   # store clinical metadata (i.e. batch and group) as a DataFrame

   clinical = pd.DataFrame({"batch": batch, "group": group})

   clinical.index = [f"sample{i}" for i in range(M)]

   # simulate data

   # random_state passes a seed to the PRNG for reproducibility

   counts = sim_rnaseq(N, M, batch=batch, group=group, random_state=42).T

We then run the two meta-analysis approaches on the obtained data.

.. repl::

   from inmoose.deseq2 import DESeq, DESeqDataSet

   from inmoose.diffexp import meta_de

   from inmoose.pycombat import pycombat_seq

   # run AD meta-analysis on batches 1, 2 and 3

   # first run the differential expression analysis on each batch individually

   cohorts = [DESeqDataSet(counts[b], clinical.loc[b], design="~ group")

              for b in [batch0, batch1, batch2]]

   individual_de = [DESeq(c).results() for c in cohorts]

   # then aggregate the results

   ad = meta_de([de for de in individual_de])

   # run ISD meta-analysis on merged batches

   # first correct batch effects to properly merge batches

   # transpositions account for the difference of formats for pycombat_seq and deseq2

   harmonized_counts = pycombat_seq(counts.T, batch).T

   # then run the differential expression analysis on the merged batches

   isd = DESeq(DESeqDataSet(harmonized_counts, clinical, design="~ group")).results()

We can now compare the results obtained by the two approaches.

.. repl-quiet::

   from matplotlib import rcParams

   # repl default config replaces '.' by '-' in the savefig.directory :/

   rcParams['savefig.directory'] = rcParams['savefig.directory'].replace("readthedocs-org", "readthedocs.org")

.. repl::

   import matplotlib.pyplot as plt

   from scipy.stats import pearsonr

   ax = plt.gca()

   ax.axline((0,0), slope=1, ls="--", lw=0.3, c=".4")

   ax.scatter(isd["log2FoldChange"], ad["combined logFC"], s=3)

   corr = pearsonr(isd["log2FoldChange"], ad["combined logFC"]).statistic

   ax.annotate(f"Pearson correlation: {100*corr:.1f}%", (2, -4), size="small")

   ax.set_title("ISD", loc="center")

   ax.set_ylabel("AD", loc="center")

   plt.show()

It is possible to combine results obtained from different tools, as long as the

results of the differential expression analysis are stored as

:class:`DEResults`.  All three modules :doc:`limma`, :doc:`edgepy` and

:doc:`deseq` return sub-classes of :class:`DEResults`, thus allowing users to

perform cross-technology meta-analysis (*e.g.* by combining results from

:doc:`limma` with results from :doc:`deseq`).

References

==========

.. [Colange2024] M. Colange, G. Appé, L. Meunier, S. Weill, A. Nordor, A.

   Behdenna. 2024. Differential Expression Analysis with InMoose, the Integrated

   Multi-Omic Open-Source Environment in Python. *Bioarxiv*.

   :doi:`10.1101/2024.11.14.623578`

Code documentation

==================

.. autosummary::

   :toctree: generated/

   ~DEResults.DEResults

   meta_de