Association analysis

Introduction

While BioFindr is developed primarily for causal inference from genomics and transcriptomics data, association analysis between genomics and transcriptomics data is also possible. In association analysis, genetic effects on the transcriptome are measured by testing if genes are differentially expressed in different groups of samples defined by the genotype of a genetic variant of interest. In BioFindr, significance of association is computed using a categorical model and a variant-specific background distribution. Similar to what was done in the coexpression analysis tutorial, this is achieved by modelling the distribution of association values between a given variant \(A\) and all genes \(B\) as a mixture distribution of real and null (random) associations. The relative weight of each component then reflects the prior probability of finding a non-null \(B\) gene for a given variant \(A\), and is fitted for every \(A\) separately.

We will illustrate how to run association 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 BioFindr

Load data

Expression data

BioFindr expects that expression data are stored as floating-point numbers in a DataFrame where columns correspond to variables (genes) and rows to samples, see the coexpression analysis tutorial for more details.

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")));

Genotype data

BioFindr expects that genotype data are stored as integer numbers in a DataFrame where columns correspond to variables (genetic variants) and rows to samples. Since BioFindr uses a categorical association model, it does not matter how different genotypes (e.g. heterozygous vs. homozygous) are encoded as integers. Future versions will support scientific types for representing genotype data.

This tutorial uses two tables of genotype data from the same set of samples as the expression data, one with genotypes for mRNA eQTLs called dgt, and one for microRNA (miRNA) eQTLs called dgm:

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

Run BioFindr

Assume we are interested in identifying mRNA genes whose expression levels are associated to microRNA eQTLs. We run:

dP = findr(dt, dgm, FDR=0.05)

28×4 DataFrame

3 rows omitted

Row	Source	Target	Probability	qvalue
	String	String	Float64	Float64
1	rs768533	MDH1B	1.0	5.95578e-9
2	rs768533_1	MDH1B	1.0	5.95578e-9
3	rs73236618	TLR6	1.0	1.61902e-7
4	rs5743580	TLR6	0.999999	2.59411e-7
5	rs73236618	TLR1	0.999999	4.82133e-7
6	rs5743580	TLR1	0.999998	7.52003e-7
7	rs1059264	SLC25A29	0.999814	2.72175e-5
8	rs3744749	PRR34-AS1	0.999749	5.51865e-5
9	rs5743580	TLR10	0.999688	8.37096e-5
10	rs73236618	TLR10	0.999666	0.00010875
11	rs73236618	FAM114A1	0.999349	0.000158079
12	rs5743580	FAM114A1	0.999102	0.000219725
13	rs3744749	MIRLET7BHG	0.998505	0.000317786
⋮	⋮	⋮	⋮	⋮
17	rs768533	TMEM60	0.917629	0.0186698
18	rs768533_1	TMEM60	0.917629	0.0222088
19	rs768533	TOR3A	0.913258	0.0256053
20	rs768533_1	TOR3A	0.913258	0.0286621
21	rs768533	CD19	0.911178	0.0315268
22	rs768533_1	CD19	0.911178	0.0341311
23	rs4926170	ZFP64	0.900374	0.0369787
24	rs73054305	GIPR	0.89784	0.0396946
25	rs768533	GBF1	0.892085	0.0424234
26	rs768533_1	GBF1	0.892085	0.0449423
27	rs768533	ATP5F1A	0.891861	0.0472829
28	rs768533_1	ATP5F1A	0.891861	0.0494564

BioFindr computes a posterior probability of non-zero association for every Source variant (columns of dgm) and Target gene (columns of dt). 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) less than a desired value (here set to 5%). The qvalue column in the output can be used for further filtering of the output, see the coexpression analysis tutorial for further details.

Note the order of the arguments. The first argument dt is the Target DataFrame, and the second argument the Source DataFrame.