Coexpression analysis

Introduction

While [BioFindr] is developed primarily for causal inference from genomics and transcriptomics data, coexpression analysis of transcriptomics data alone is also possible. In coexpression analysis, pairwise correlation is used as a simple measure for the probability of two genes being functionally related either through direct or indirect regulation, or through coregulation by a third factor. In BioFindr, significance of pairwise correlation is computed using a gene-specific background distribution, allowing for genes having different biological roles. For instance, it is known that many biological networks are scale-free, where a small number of so-called “hub genes” have a high number of interaction partners while most other genes only have few. In BioFindr, this is accomodated by modelling the distribution of correlation values between a given gene \(A\) and all other genes \(B\) as a 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.

We will illustrate how to run coexpression analysis with BioFindr using preprocessed data from the GEUVADIS study. See the installation instructions for the steps you need to take to reproduce this tutorial.

Set up the environment

We begin by setting up the environment and loading some necessary packages.

using DrWatson
using DataFrames
using Arrow
using Markdown

using BioFindr

Load expression data

BioFindr expects that expression data are stored in a DataFrame where columns correspond to variables (genes) and rows to samples. An expression DataFrame should not contain any other columns (e.g. gene annotation) than gene expression columns, and gene expression data should be stored as floating-point numbers. Internally, BioFindr operates on matrices, and if you have an expression DataFrame df, then Matrix(df) should return a matrix of floats.

At the moment, BioFindr does not support count-based expression data being provided as a DataFrame of integers. This is not an intrinsic limitation of the method, but simply to distinguish expression data from integer-valued genotype data. Future versions will remove this limitation by supporting scientific types.

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:

dt = DataFrame(Arrow.Table(datadir("exp_pro","findr-data-geuvadis", "dt.arrow")));
dm = DataFrame(Arrow.Table(datadir("exp_pro","findr-data-geuvadis", "dm.arrow")));

The mRNA data has expression data for 19947 genes in 360 samples, while miRNA data is available for 674 miRNAs in the same 360 samples.

We can confirm that the data frames are of the right type:

[typeof(Matrix(dt)) typeof(Matrix(dm))]
1×2 Matrix{DataType}:
 Matrix{Float64}  Matrix{Float64}

Run BioFindr

All-vs-all coexpression analysis

The simplest type of coexpression analysis tests for non-zero correlation among all possible pairs in a gene expression dataset. Let’s do this for the miRNA data:

dP_mir_all = findr(dm, FDR=0.05)
30827×4 DataFrame
30802 rows omitted
Row Source Target Probability qvalue
String String Float64 Float64
1 hsa-miR-574-3p hsa-miR-574-5p 1.0 0.0
2 hsa-miR-4482-1-3p hsa-miR-1304-3p 1.0 0.0
3 hsa-miR-3667-5p hsa-miR-3667-3p 1.0 0.0
4 hsa-miR-550a-2-5p hsa-miR-550a-2-3p 1.0 0.0
5 hsa-miR-550a-2-5p hsa-miR-550a-3-5p 1.0 0.0
6 hsa-miR-3130-2-5p hsa-miR-3130-1-5p 1.0 0.0
7 hsa-let-7b-5p hsa-let-7b-3p 1.0 0.0
8 hsa-miR-96-5p hsa-miR-182-5p 1.0 0.0
9 hsa-miR-96-5p hsa-miR-1271-5p 1.0 0.0
10 hsa-miR-3667-3p hsa-miR-3667-5p 1.0 0.0
11 hsa-miR-181c-5p hsa-miR-181c-3p 1.0 0.0
12 hsa-miR-181c-5p hsa-miR-181d-5p 1.0 0.0
13 hsa-miR-335-5p hsa-miR-335-3p 1.0 0.0
⋮ ⋮ ⋮ ⋮ ⋮
30816 hsa-miR-769-5p hsa-miR-25-5p 0.858747 0.0499658
30817 hsa-miR-3133-5p hsa-miR-1976-3p 0.858744 0.0499688
30818 hsa-miR-320c-2-3p hsa-miR-1260a-5p 0.858739 0.0499718
30819 hsa-miR-378a-3p hsa-miR-24-1-3p 0.858739 0.0499747
30820 hsa-miR-642a-5p hsa-miR-744-5p 0.858737 0.0499777
30821 hsa-miR-107-3p hsa-miR-152-3p 0.858735 0.0499807
30822 hsa-miR-17-3p hsa-miR-744-3p 0.858733 0.0499836
30823 hsa-miR-362-5p hsa-miR-616-5p 0.858731 0.0499866
30824 hsa-miR-223-3p hsa-miR-3688-1-3p 0.858729 0.0499895
30825 hsa-miR-21-3p hsa-miR-1260b-5p 0.858729 0.0499925
30826 hsa-miR-4517-5p hsa-miR-365a-3p 0.858728 0.0499955
30827 hsa-miR-342-3p hsa-miR-93-3p 0.858715 0.0499984

BioFindr computes a posterior probability of non-zero correlation for every Source and Target gene pair. By default the output is sorted by decreasing Probability and self-interactions are excluded. The optional parameter FDR can be used to limit the output to the set of pairs that has a global false discovery rate (FDR) less than a desired value. The qvalue column in the output can be used for further filtering of the output. Say you ran findr with an FRD threshold of 5% as above. If you now want to restrict the output to an FDR threshold of 1%, you can do:

filter!(row -> row.qvalue <= 0.01, dP_mir_all)
14682×4 DataFrame
14657 rows omitted
Row Source Target Probability qvalue
String String Float64 Float64
1 hsa-miR-574-3p hsa-miR-574-5p 1.0 0.0
2 hsa-miR-4482-1-3p hsa-miR-1304-3p 1.0 0.0
3 hsa-miR-3667-5p hsa-miR-3667-3p 1.0 0.0
4 hsa-miR-550a-2-5p hsa-miR-550a-2-3p 1.0 0.0
5 hsa-miR-550a-2-5p hsa-miR-550a-3-5p 1.0 0.0
6 hsa-miR-3130-2-5p hsa-miR-3130-1-5p 1.0 0.0
7 hsa-let-7b-5p hsa-let-7b-3p 1.0 0.0
8 hsa-miR-96-5p hsa-miR-182-5p 1.0 0.0
9 hsa-miR-96-5p hsa-miR-1271-5p 1.0 0.0
10 hsa-miR-3667-3p hsa-miR-3667-5p 1.0 0.0
11 hsa-miR-181c-5p hsa-miR-181c-3p 1.0 0.0
12 hsa-miR-181c-5p hsa-miR-181d-5p 1.0 0.0
13 hsa-miR-335-5p hsa-miR-335-3p 1.0 0.0
⋮ ⋮ ⋮ ⋮ ⋮
14671 hsa-miR-598-3p hsa-miR-505-3p 0.964104 0.00998035
14672 hsa-miR-99b-5p hsa-miR-1180-3p 0.964099 0.00998211
14673 hsa-miR-29b-2-3p hsa-miR-19b-1-3p 0.964094 0.00998388
14674 hsa-miR-150-3p hsa-miR-744-5p 0.964081 0.00998565
14675 hsa-miR-30a-3p hsa-miR-361-3p 0.964079 0.00998741
14676 hsa-miR-17-3p hsa-miR-3940-3p 0.964078 0.00998918
14677 hsa-miR-484-5p hsa-miR-92a-1-3p 0.964076 0.00999095
14678 hsa-miR-652-3p hsa-miR-100-5p 0.964074 0.00999272
14679 hsa-miR-378i-5p hsa-miR-425-3p 0.96407 0.00999448
14680 hsa-miR-212-3p hsa-miR-222-5p 0.964061 0.00999625
14681 hsa-miR-320b-1-3p hsa-miR-4797-3p 0.964053 0.00999802
14682 hsa-miR-23a-3p hsa-miR-7-3-5p 0.964052 0.00999979

Note that the filter! command modifies the input DataFrame in-place, that is, the rows not matching the selection criteria are deleted. Use filter to return a new DataFrame with the selected rows.

Finally, remember that the output of coexpression analysis in BioFindr is not symmetric, that is

\[ P(Source, Target) \neq P(Target, Source) \]

This is because the posterior probabilities are estimated using a Source-specific background distribution, accounting for the fact that different genes may have a different number of non-null interaction partners a priori. See the Findr paper for details.

Bipartite coexpression analysis

Since BioFindr’s posterior probabilities are Source gene-specific, we may be interested in computing probabilities only for a subset of Source genes, or using different Source and Target gene sets.

As an example of the first scenario, assume we are interested in finding microRNAs that are significantly correlated with microRNAs from the Mir-200 family. First find the Mir-200 microRNAs:

mir200 = names(dm)[startswith.(names(dm),"hsa-miR-200")]
3-element Vector{String}:
 "hsa-miR-200b-3p"
 "hsa-miR-200a-3p"
 "hsa-miR-200c-3p"

Then run

dP_mir200_mir = findr(dm, colnames=mir200, FDR=0.01)
48×4 DataFrame
23 rows omitted
Row Source Target Probability qvalue
String String Float64 Float64
1 hsa-miR-200b-3p hsa-miR-200a-3p 1.0 2.0303e-10
2 hsa-miR-200a-3p hsa-miR-200b-3p 1.0 2.37622e-10
3 hsa-miR-200c-3p hsa-miR-93-5p 0.999999 2.26212e-7
4 hsa-miR-200c-3p hsa-miR-140-3p 0.999991 2.36372e-6
5 hsa-miR-200b-3p hsa-miR-429-3p 0.999927 1.65171e-5
6 hsa-miR-200c-3p hsa-miR-16-1-5p 0.999845 3.95838e-5
7 hsa-miR-200c-3p hsa-miR-664-3p 0.999839 5.69759e-5
8 hsa-miR-200c-3p hsa-miR-17-5p 0.999719 8.50297e-5
9 hsa-miR-200c-3p hsa-miR-769-5p 0.999366 0.000146048
10 hsa-miR-200b-3p hsa-miR-582-3p 0.999286 0.000202888
11 hsa-miR-200c-3p hsa-miR-130b-5p 0.999273 0.000250569
12 hsa-miR-200c-3p hsa-miR-425-5p 0.999193 0.000296912
13 hsa-miR-200c-3p hsa-miR-103a-2-3p 0.999146 0.00033977
⋮ ⋮ ⋮ ⋮ ⋮
37 hsa-miR-200a-3p hsa-miR-3667-3p 0.983748 0.00636739
38 hsa-miR-200b-3p hsa-miR-26b-5p 0.982875 0.00665048
39 hsa-miR-200c-3p hsa-miR-619-5p 0.982006 0.00694135
40 hsa-miR-200b-3p hsa-miR-328-3p 0.98103 0.00724207
41 hsa-miR-200c-3p hsa-miR-500a-5p 0.980897 0.00753137
42 hsa-miR-200b-3p hsa-miR-9-1-5p 0.980766 0.00781
43 hsa-miR-200c-3p hsa-miR-4677-3p 0.980721 0.00807672
44 hsa-miR-200c-3p hsa-miR-141-5p 0.979668 0.00835524
45 hsa-miR-200c-3p hsa-miR-30e-3p 0.979068 0.00863472
46 hsa-miR-200b-3p hsa-miR-93-5p 0.973906 0.00901427
47 hsa-miR-200c-3p hsa-miR-25-3p 0.972582 0.00940584
48 hsa-miR-200c-3p hsa-miR-4521-5p 0.972378 0.00978535

The parameter colnames must be a vector of strings containing a subset of variable names of the input DataFrame dm to be used as Source genes.

As an example of the second scenario, we may be interested in finding genes that are significantly correlated with all or a subset of microRNAs:

dP_mir_mrna = findr(dt, dm, FDR=0.01)
49499×4 DataFrame
49474 rows omitted
Row Source Target Probability qvalue
String String Float64 Float64
1 hsa-miR-3150b-3p MIR3150BHG 1.0 0.0
2 hsa-miR-3150b-3p PLEKHF2 1.0 0.0
3 hsa-let-7b-5p MIRLET7BHG 1.0 0.0
4 hsa-let-7b-5p PRR34 1.0 0.0
5 hsa-let-7b-5p PRR34-AS1 1.0 0.0
6 hsa-miR-335-5p CPA4 1.0 0.0
7 hsa-miR-335-5p MESTP1 1.0 0.0
8 hsa-miR-618-5p LIN7A 1.0 0.0
9 hsa-miR-335-3p CPA4 1.0 0.0
10 hsa-miR-335-3p MESTP1 1.0 0.0
11 hsa-let-7b-3p PRR34-AS1 1.0 0.0
12 hsa-miR-10a-5p HOXB3 1.0 0.0
13 hsa-miR-143-3p ALDH1L1 1.0 0.0
⋮ ⋮ ⋮ ⋮ ⋮
49488 hsa-miR-577-5p PEBP1 0.973475 0.00999603
49489 hsa-miR-193b-5p DTL 0.973474 0.00999637
49490 hsa-miR-944-3p RAPGEF2 0.973473 0.0099967
49491 hsa-miR-664-3p RPS7P10 0.973473 0.00999704
49492 hsa-miR-193a-5p RPL12P13 0.973471 0.00999737
49493 hsa-miR-21-3p ITPRID2-DT 0.973471 0.0099977
49494 hsa-miR-181c-3p EXOSC10 0.97347 0.00999804
49495 hsa-miR-151a-3p ZNF664 0.973469 0.00999837
49496 hsa-miR-671-3p GRINA 0.973468 0.00999871
49497 hsa-miR-25-5p VANGL2 0.973467 0.00999904
49498 hsa-miR-22-3p LILRB4 0.973466 0.00999937
49499 hsa-miR-125a-5p POLD4 0.973465 0.00999971

Note the order of the arguments: here we tested all microRNAs as \(A\) or Source genes (dm argument) against all mRNA transcripts as \(B\) or Target genes (dt argument), that is, background distributions are fitted for each microRNA (column of dm) from the log-likelihood ratios for all 23,722 mRNAs (columns of dt).