causal-inference.ipynb

{
  "cells": [
    {
      "cell_type": "raw",
      "metadata": {},
      "source": [
        "---\n",
        "title: Causal inference\n",
        "format:\n",
        "  html:\n",
        "    code-fold: false\n",
        "execute:\n",
        "  cache: true\n",
        "  keep-ipynb: true\n",
        "---"
      ],
      "id": "e3a40232"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Introduction\n",
        "\n",
        "The primary use of [BioFindr][1] is for causal inference from genomics and transcriptomics data. In causal inference, cis-acting eQTLs are used as causal anchors or instrumental variables to orient the direction of causality between coexpressed genes, see the [Findr paper][5] for more details. As in [coexpression analysis](coexpression.qmd) and [association analysis](association.qmd), the significance of a causal effect is computed using a gene-specific background distribution. This is again achieved by modelling the distribution of test statistics between a given gene $A$ and all other genes $B$ as a [mixture distribution](https://en.wikipedia.org/wiki/Mixture_distribution) of real and null (random) correlations. The relative weight of each component reflects the prior probability of finding a non-null $B$ gene for a given $A$ gene, and is fitted for every $A$ gene separately.\n",
        "\n",
        "Unlike in [coexpression analysis](coexpression.qmd) and [association analysis](association.qmd), causal inference cannot be performed using a single statistical test, but requires the combination of multiple tests. In [BioFindr][1], [six likelihood ratio tests](https://tmichoel.github.io/BioFindr.jl/dev/realLLR/) are implemented, which can be combined in multiple ways for causal inference. Tests are combined by addition or mulitplication of the [posterior probabilities](https://tmichoel.github.io/BioFindr.jl/dev/posteriorprobs/) of individial tests, an approach first proposed in [this paper](https://doi.org/10.1186/gb-2007-8-10-r219).\n",
        "\n",
        "\n",
        "We will illustrate how to run causal inference with [BioFindr][1] using [preprocessed data][2] from the [GEUVADIS study][3]. See the [installation instructions](installation.qmd) for the steps you need to take to reproduce this tutorial.\n",
        "\n",
        "## Set up the environment\n"
      ],
      "id": "1cc8f3b6"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "using DrWatson\n",
        "quickactivate(@__DIR__)\n",
        "\n",
        "using DataFrames\n",
        "using Arrow\n",
        "\n",
        "using BioFindr"
      ],
      "id": "1c9b5546",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Load data\n",
        "\n",
        "### Expression data\n",
        "\n",
        "[BioFindr][1] expects that expression data are stored as [floating-point numbers](https://docs.julialang.org/en/v1/manual/integers-and-floating-point-numbers/) in a [DataFrame][4] where columns correspond to variables (genes) and rows to samples, see the [coexpression analysis tutorial](coexpression.qmd) for more details.\n",
        "\n",
        "This tutorial uses two tables of expression data from the same set of samples, one for mRNA expression data called `dt`, and one for microRNA (miRNA) expression data called `dm`:\n"
      ],
      "id": "eaae4e84"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "dt = DataFrame(Arrow.Table(datadir(\"exp_pro\",\"findr-data-geuvadis\", \"dt.arrow\")));\n",
        "dm = DataFrame(Arrow.Table(datadir(\"exp_pro\",\"findr-data-geuvadis\", \"dm.arrow\")));"
      ],
      "id": "9fa3d624",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Genotype data\n",
        "\n",
        "[BioFindr][1] expects that genotype data are stored as [integer numbers](https://docs.julialang.org/en/v1/manual/integers-and-floating-point-numbers/) in a [DataFrame][4] where columns correspond to variables (genetic variants) and rows to samples. [BioFindr][1] uses a categorical model to associate genetic variation to variation in gene expression, hence how different genotypes (e.g. heterozygous vs. homozygous) are encoded as integers does not matter. Future versions will support [scientific types](https://juliaai.github.io/ScientificTypes.jl/dev/) for representing genotype data.\n",
        "\n",
        "This tutorial uses two tables of expression data from the same set of samples, one with genotypes for mRNA eQTLs called `dgt`, and one for microRNA (miRNA) eQTLs called `dgm`:\n"
      ],
      "id": "caa4c9e4"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "dgt = DataFrame(Arrow.Table(datadir(\"exp_pro\",\"findr-data-geuvadis\", \"dgt.arrow\")));\n",
        "dgm = DataFrame(Arrow.Table(datadir(\"exp_pro\",\"findr-data-geuvadis\", \"dgm.arrow\")));"
      ],
      "id": "b3646576",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### eQTL mapping data\n",
        "\n",
        "The [preprocessed GEUVADIS data][2] has been organized such that each column of the genotype data is the strongest eQTLs for the corresponding column in the matching expression data. Usually however, eQTL mapping data will be available in the form of a table with variant IDs, gene IDs, and various eQTL associaion statistics (see the [original GEUVADIS file](https://www.ebi.ac.uk/biostudies/files/E-GEUV-1/E-GEUV-1/analysis_results/EUR373.gene.cis.FDR5.all.rs137.txt.gz) for an example). [BioFindr][1] expects that such a table is read into a [DataFrame][4], and that only the best (most strongly associated) eQTL is kept for each gene, that is, genes appear only once in the eQTL mapping table. Let's artificially create such tables for our data:\n"
      ],
      "id": "c4876037"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "dpt = DataFrame(SNP_ID = names(dgt), GENE_ID=names(dt)[1:ncol(dgt)]);\n",
        "dpm = DataFrame(SNP_ID = names(dgm), GENE_ID=names(dm)[1:ncol(dgm)]);"
      ],
      "id": "abc30770",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Run BioFindr\n",
        "\n",
        "### Subset-to-all causal inference\n",
        "\n",
        "In the most common scenario, we have a dataset of gene expression values, with significant cis-eQTL instruments for a subset of genes (\"eGenes\"). We can then infer causal interactions from the eGenes to all other genes (including other eGenes). For instance, to infer causal microRNA $\\to$ microRNA interaction, run:\n"
      ],
      "id": "4bebe51f"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "dP = findr(dm, dgm, dpm; FDR=0.25)"
      ],
      "id": "3cdfa059",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "BioFindr computes a [posterior probability](https://tmichoel.github.io/BioFindr.jl/dev/posteriorprobs/) of a non-zero causal interaction for every pair pf **Source** and **Target** genes (columns of `dm`). The possible **Source** genes are the eGenes, the subset of genes with cis-eQTLs (columns of `dgm`), as defined by the eQTL-to-gene mapping in `dpm`.\n",
        "\n",
        "By default the output is sorted by decreasing **Probability**. The optional parameter **FDR** can be used to limit the output to the set of pairs that has a [global false discovery rate (FDR)](https://en.wikipedia.org/wiki/False_discovery_rate#Storey-Tibshirani_procedure) less than a desired value (here set to 10%). The **qvalue** column in the output can be used for further filtering of the output, see the [coexpression analysis tutorial](coexpression.qmd) for more details.\n",
        "\n",
        "Note the order of input arguments: first `dm`, the expression data, then `dgm`, the genotype data, and then `dpm`, the eQTL mapping of variants to eGenes.\n",
        "\n",
        "By default, BioFindr assumes that the first column of `dpm` is a list of variant names that can be found in the column names of `dgm`, and that the second column of `dpm` is a list of gene names that can be found in the column names of `dm`:\n"
      ],
      "id": "c50813bc"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "dpm"
      ],
      "id": "3a33a6dc",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "If your eQTL mapping DataFrame has the relevant columns in a different place, you need to use the optional arguments `colX` (for the gene expression IDs) and `colG` (for the eQTL genotype IDs) to specify either the relevant columns index or name:\n"
      ],
      "id": "5d5b62d0"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "dP = findr(dm, dgm, dpm; colX=2, colG=1, FDR=0.1);"
      ],
      "id": "aa24636d",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "or\n"
      ],
      "id": "dd2004e1"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "dP = findr(dm, dgm, dpm; colX=\"GENE_ID\", colG=\"SNP_ID\", FDR=0.1);"
      ],
      "id": "f64ccd68",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Any other columns in the eQTL mapping DataFrame (such as association summary statistics) are ignored.\n",
        "\n",
        "### Bipartite causal inference\n",
        "\n",
        "In some applications, we have multiple omics datasets from the same samples, and may be interested in causal inference from one set of variables to another one. In the GEUVADIS study for instance, we may be interested in inferring causal interactions between microRNAs and target genes. To infer interactions *from* microRNAs *to* mRNAs, run:\n"
      ],
      "id": "6e9f1eb8"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "dP = findr(dt, dm, dgm, dpm; FDR=0.1)"
      ],
      "id": "29006ba5",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Causal test combinations\n",
        "\n",
        "As indicated above, causal inference requires the combination of multiple statistical tests, and [BioFindr][1] supports [multiple such combinations](https://tmichoel.github.io/BioFindr.jl/dev/testLLR/#Causal-inference), each with their own strengths and weaknesses. Comparisons between the different methods can be found in the [original BioFindr paper][5] and in a [follow-up paper][6]. In brief, the following combinations are available:\n",
        "\n",
        "#### Instrumental variables (default)\n",
        "\n",
        "The *instrumental variables* (IV) test combination infers a causal $A\\to B$ interaction if there is a genetic association $E\\to B$, where $E$ is the causal anchor (best cis-eQTL for $A$), and if $A$ and $B$ are not independently associated with $E$ (no pleiotropic effect). This test combination has been shown to provide a good balance between false negative and false positive rate and has a clear genetic interpretation. It is therefore recommended as the default combination. To use the instrumental variable test combination, you don't need to do anything (it's the default): \n"
      ],
      "id": "e0f1797b"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "dP = findr(dm, dgm, dpm; FDR=0.1);"
      ],
      "id": "d5fe34f4",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "or, you can explicitly set the optional parameter `combination=\"IV\"`:\n"
      ],
      "id": "83bf22dc"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "dP = findr(dm, dgm, dpm; FDR=0.1, combination=\"IV\");"
      ],
      "id": "0bb8cc0d",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "#### Mediation\n",
        "\n",
        "The *mediation test* infers a causal $A\\to B$ interaction if there is a genetic association $E\\to B$ and if this association disappears after conditioning on $A$. If these two tests are true, the causal graph must be $E\\to A\\to B$. As a result, this test combination has a very low false positive rate, but in most applications it is hampered by a high false negative rate: if $A$ and $B$ are jointly regulated by a third factor $C$, the mediation test will fail even in the presence of a true $A\\to B$ interaction. To use the mediation test combination, set the optional parameter `combination=\"mediation\"`:\n"
      ],
      "id": "a6da7cd5"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "dP_med = findr(dm, dgm, dpm; FDR=0.1, combination=\"mediation\");"
      ],
      "id": "e7a094cc",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "#### Original BioFindr combination\n",
        "\n",
        "In the [original BioFindr paper][5], a new test combination was proposed that combines the instrumental variable test combination with a new [relevance test](https://tmichoel.github.io/BioFindr.jl/dev/realLLR/#Relevance-test). On simulated data, this new combination further reduced the false negative rate compared to the IV combination. at the expense of an increased false positive rate, for an overall improved performance. However, the relevance test does not have a clear genetic or causal interpretation, and on real data, performance is generally equivalent to the IV combination, and therefore the IV combination is now recommended as the default. To use the original BioFindr test combination, set the optional parameter `combination=\"orig\"`:\n"
      ],
      "id": "c8f38f56"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "dP_orig = findr(dm, dgm, dpm; FDR=0.1, combination=\"orig\");"
      ],
      "id": "52d636cd",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "[1]: https://github.com/tmichoel/BioFindr.jl\n",
        "[2]: https://github.com/lingfeiwang/findr-data-geuvadis\n",
        "[3]: https://doi.org/10.1038/nature12531\n",
        "[4]: https://dataframes.juliadata.org/stable/\n",
        "[5]: https://doi.org/10.1371/journal.pcbi.1005703\n",
        "[6]: https://doi.org/10.1039/D0MO00140F"
      ],
      "id": "3c477424"
    }
  ],
  "metadata": {
    "kernelspec": {
      "name": "julia-1.10",
      "language": "julia",
      "display_name": "Julia 1.10.0"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}