Exploring Bayesian approaches to eQTL mapping through probabilistic programming

TL;DR

This paper demonstrates how probabilistic programming with Stan can streamline Bayesian eQTL mapping by enabling automatic inference, thus accelerating model development and testing in genomic studies.

Contribution

It introduces a framework for automatic, black-box Bayesian inference in eQTL mapping models using Stan, simplifying the development process.

Findings

Facilitates rapid model prototyping and testing.

Enables automatic inference in complex Bayesian models.

Accelerates genomic data analysis workflows.

Abstract

The discovery of genomic polymorphisms influencing gene expression (also known as expression quantitative trait loci or eQTLs) can be formulated as a sparse Bayesian multivariate/multiple regression problem. An important aspect in the development of such models is the implementation of bespoke inference methodologies, a process which can become quite laborious, when multiple candidate models are being considered. We describe automatic, black-box inference in such models using Stan, a popular probabilistic programming language. The utilisation of systems like Stan can facilitate model prototyping and testing, thus accelerating the data modelling process. The code described in this chapter can be found at https://github.com/dvav/eQTLBookChapter.

Figures

Equations

$${\beta }_{ik}\sim (1-{\pi }_{ik}){\delta }_0+{\pi }_{ik}\mathcal{N}(0,{\sigma }_i^2)$$

$${\beta }_{ik}\sim \mathcal{N}\left(0,{\eta }^2{\zeta }_{ik}^2\right){\eta }^2\sim {\mathrm{p}}_{\eta }({\eta }^2){\zeta }_{ik}^2\sim {\mathrm{p}}_{\zeta }({\zeta }_{ik}^2)$$

$$y_{ij}\sim \mathcal{N}\left({\beta }_{0i}+\sum _k{\beta }_{ik}{\tilde{x}}_{jk},{\sigma }_i^2\right){\beta }_{0i}\sim 1\log {\sigma }_i^2\sim 1$$

$${\beta }_{ik}\sim \mathcal{N}\left(0,\frac{{\bar{\sigma }}^2}{NK}{\eta }^2{\zeta }_{ik}^2\right)\eta \sim {\mathcal{C}}^{+}(0,1){\zeta }_{ik}\sim {\mathcal{C}}^{+}(0,1)$$

$$z_{ij}\sim \mathcal{N}\mathcal{B}(m_{ij},{\varphi }_i)\log m_{ij}=\log c_{ij}+\log L_j+{\beta }_{0i}+\sum _k{\beta }_{ik}{\tilde{x}}_{jk}{\beta }_{0i}\sim 1\log {\varphi }_i\sim 1$$

$${\beta }_{ik}\sim \mathcal{N}\left(0,\frac{{\eta }^2{\zeta }_{ik}^2}{NK}\right)\eta \sim {\mathcal{C}}^{+}(0,1){\zeta }_{ik}\sim {\mathcal{C}}^{+}(0,1)$$

$$z_{ij}\sim \mathcal{P}(m_{ij})\log m_{ij}=\log c_{ij}+\log L_j+y_{ij}$$

$$z_{ij}\sim \frac{\Gamma (z_{ij}+{\alpha }_i)}{\Gamma (z_{ij}+1)\Gamma ({\alpha }_i)}{\left(\frac{c_j{\mu }_i}{c_j{\mu }_i+{\alpha }_i}\right)}^{z_{ij}}{\left(\frac{{\alpha }_i}{c_j{\mu }_i+{\alpha }_i}\right)}^{{\alpha }_i}$$

Full text

Keywords: Bayesian Variable Selection, Global-Local Shrinkage, Horseshoe Prior, RNA-seq, eQTL Mapping, Probabilistic Programming, Stan, R, Black-Box Bayesian Inference.

1 Introduction

The study of genomic variation and its association with gene expression is critical for elucidating the genetic basis of complex traits, including diseases. The advent of next-generation sequencing (NGS) made possible the detailed investigation of this relationship (also known as eQTL mapping) in large cohorts, but it also gave rise to novel statistical challenges[1, 2]. eQTL mapping can be examined in the context of sparse Bayesian multivariate/multiple regression, where we typically propose a number of candidate statistical models followed by benchmarking them against each other in terms of computational and statistical efficiency[3]. The most laborious aspect of this process is the development and software implementation of statistical inference algorithms for each model under consideration, a task that can be impeded by the potential absence of favourable mathematical properties (e.g. conjugacy) in any of these models.

The aim of this chapter is to demonstrate the utility of a popular probabilistic programming language (PPL), Stan, in the prototyping and testing phases of the data modelling process for eQTL mapping[4]. PPLs make it possible to describe and perform inference in hierarchical probabilistic models in a matter of minutes using a small amount of high-level code. Besides Stan, another popular PPL is PyMC, a Python-based software offering excellent performance, a wide range of inference algorithms, and the ability to mix these algorithms in the same inferential cycle for estimating different parts of a given model[5]. The reason for choosing Stan in this chapter is its simplicity, which stems from the fact that its syntax is very close to the mathematical notation already familiar to statisticians.

The practical part of this chapter covers: (a) the acquisition of genomic variation and gene expression data from online sources, (b) the use of these to simulate artificial eQTL datasets with known properties, (c) the implementation of statistical models in Stan for eQTL mapping, and (d) the estimation of unknown model parameters using the previously simulated datasets. In the remaining part of this introduction, we outline the statistical theory, which underpins the practical part of this chapter.

1.1 Theory

We assume an $N\times M$ matrix $Z=\{z_{ij}\}$ of read counts quantifying the expression of $N$ transcripts in $M$ samples, an $N\times M$ matrix $C=\{c_{ij}\}$ of transcript- and sample-specific normalization factors, and an $M\times K$ matrix $X=\{x_{jk}\}$ of genotypes indicating the number of minor alleles (0, 1 or 2) in each of $K$ bi-allelic genomic loci in $M$ samples. A matrix $\tilde{X}=\{\tilde{x}_{jk}\}$ is derived by standardizing each column of $X$.

We introduce an $N\times K$ sparse matrix $B$ of regression coefficients with elements $\beta_{ik}$, which measure the effect of variant $k$ on the expression of transcript $i$. Estimating $B$ is the main focus of subsequent analysis. Typically, sparsity is induced on $B$ by imposing appropriate priors. A common sparsity-inducing prior is a two-component mixture of the following form:

[TABLE]

where $\delta_{0}$ is a point-mass distribution centred at 0. This prior can set $\beta_{ik}$ exactly at 0, but posterior estimation requires stochastic exploration of a $2^{NK}$-dimensional discrete space. An alternative approach with obvious computational advantages would be to adopt a continuous shrinkage prior; for example, the important class of global-local shrinkage priors[6]:

[TABLE]

where $\zeta_{ik}$ are local (i.e. gene- and variant-specific) shrinkage scales, while $\eta$ is a global scale controlling the overall shrinkage of $B$. The shrinkage profile of $B$ depends on the form of $\mathrm{p}_{\eta}$ and $\mathrm{p}_{\zeta}$. Different choices of these distributions give rise to different shrinkage priors, but here we shall use the horseshoe prior[7] for which $\zeta_{ik}$ and $\eta$ follow standard and scaled half-Cauchy distributions, $\mathcal{C}^{+}(0,1)$ and $\mathcal{C}^{+}(0,\alpha)$, respectively.

1.1.1 Normal model

For the Normal model, we assume that read counts have been normalized, $\tilde{z}_{ij}=\frac{z_{ij}}{c_{ij}}$, and log-transformed, $y_{ij}=\log(\tilde{z}_{ij}+1)$. The model takes the following form:

[TABLE]

[TABLE]

where $\sigma^{2}_{i}$ is the variance of gene $i$ and $\sum_{k}\beta_{ik}\tilde{x}_{jk}$ is the effect of genomic variation on the baseline expression $\beta_{0i}$ of gene $i$. We assume that $M$ is large, so we can afford to impose a flat prior distribution on $\beta_{0i}$ and $\log\sigma_{i}$ over the interval ($-\infty,+\infty$). A more complicated model would include correlations between different $y_{ij}$ variables, additional (e.g. clinical, environmental and population) co-variates influencing the baseline gene expression $\beta_{0i}$, as well as hierarchical priors on $\beta_{0i}$ and $\sigma^{2}_{i}$. Notice that the variance of $\beta_{ik}$ is proportional to $\bar{\sigma}^{2}=\left(\frac{\sum_{i}\sigma_{i}}{N}\right){{}^{2}}$ and inversely proportional to the total number of coefficients. The implicit assumption under this formulation is that the true global scale parameter is $\xi\sim\mathcal{C}^{+}\left(0,\frac{\bar{\sigma}}{\sqrt{NK}}\right)$. Finally, in this and the subsequent models, we assume that we can ignore any gene-specific factors (e.g. length) affecting $c_{ij}$, thus $c_{ij}\equiv c_{j}$.

1.1.2 Negative Binomial model

The Negative Binomial distribution is immensely popular for modelling over-dispersed RNA-seq data[8, 9, 10], but the mathematical complexities associated with inference in this model might explain (at least partially) the popularity of transformation-based methods, such as voom[11]. Here, we examine the following model:

[TABLE]

[TABLE]

where $m_{ij}$ is the gene- and sample-specific mean of the Negative Binomial distribution, and $\phi_{i}$ is a gene-specific dispersion parameter, such that $\mathrm{Var}[z_{ij}]=m_{ij}+m^{2}_{ij}\phi_{i}$; $L_{j}=\sum_{i}z_{ij}$ is the total number of reads in sample $j$.

1.1.3 Poisson-LogNormal model

An alternative approach to work with non-transformed data is to impose a Poisson observational model on top of the Normal:

[TABLE]

where $y_{ij}$ serve as latent variables following the Normal model. The Poisson-LogNormal model is motivated by (a) the fact that the Negative Binomial model can be thought of as a Poisson-Gamma mixture, and (b) replacing the Gamma distribution in the previous mixture with LogNormal. By invoking the law of total expectation and the law of total variance, we can see that $E[z_{ij}]=c_{ij}L_{j}e^{\beta_{0i}+\sum_{k}\beta_{ik}\tilde{x}_{jk}+\sigma_{i}^{2}/2}$ and $\mathrm{Var}[z_{ij}]=E[z_{ij}]+E[z_{ij}]{{}^{2}}\phi_{i}$. Hence, the variance has the same form as in the Negative Binomial model, with dispersion parameter $\phi_{i}=e^{\sigma^{2}_{i}}-1$. When $\sigma^{2}_{i}=0$, the above model reduces to Poisson.

2 Materials

2.1 Operating system

A working UNIX environment (e.g. Linux or MacOS X) with a terminal emulator running the command shell bash. The work presented in this chapter was tested on Ubuntu Linux v18.04.

2.2 Software

A recent version of the command line tool and library curl for transferring data with URLs. 2. 2.

A recent version of vcftools, a set of tools written in Perl and C++ for working with VCF files[12]. 3. 3.

A recent version of R, the free software environment for statistical computing[13]. 4. 4.

A recent version of rstudio, an integrated development environment for working with R and the command line[14]. 5. 5.

A recent version of rstan, an R interface for stan. 6. 6.

A recent version of plyr, a set of R tools for splitting, modifying and combining data[15]. 7. 7.

A recent version of doMC, an R package providing a parallel backend for multicore computation. 8. 8.

A recent version of cowplot, an R package for plotting. 9. 9.

A recent version of tidyverse, a collection of R packages for data wrangling and plotting. 10. 10.

A recent version of reshape2, an R package for restructuring and aggregating data[16].

The above R packages can be installed either through the graphical interface provided by rstudio, or through the R console; for example install.packages(’tidyverse’)

3 Methods

3.1 Data acquisition

Start rstudio 2. 2.

Create a working directory tree by typing the following commands at the bash command prompt11endnote: 1 You can access the command prompt through a terminal emulator. From within rstudio, you can create a terminal emulator by going to Tools>Terminal>New Terminal. :

⬇

1mkdir -p eQTLchapter/{data,R,stan} 3. 3.

Make the root of the tree you just created your working directory by typing the following at the R console22endnote: 2 Instead of steps 2 and 3, you can just create a new rstudio project (File>New Project…) with equivalent results :

⬇

1setwd(’eQTLchapter’) 4. 4.

Download genomic variation data from the 1000 Genomes project[17]. At the bash command prompt, type the following33endnote: 3 This chunk of code uses curl to stream data from the 1000 Genomes Project and pipe them (in gzipped VCF format) to vcftools for further processing. Only variants in chromosome 7 between positions 100K and 200K are retained. In addition, all variants with FILTER other than PASS are removed; indels are removed; variants in non-bi-allelic sites are removed; only variants with allele frequency higher than 5% are retained (i.e. 0.05<MAF<0.95). Filtered data are stored as a matrix with elements 0, 1 and 2 indicating the number of alleles in 496 variants across 2504 samples. If vcftools throws an error about not being able to open temporary files, you need to increase the maximum number of open files allowed by your system to more than double the number of samples in the VCF file. At the bash command prompt, type: ulimit -n 5100. You may need root permissions to run this command. :

⬇

1cd eQTLchapter

2curl -o - link \

3 | vcftools --gzvcf - \

4 --chr 7 --from-bp 100000 --to-bp 200000 \

5 --remove-filtered-all --remove-indels \

6 --min-alleles 2 --max-alleles 2 \

7 --maf 0.05 --max-maf 0.95 \

8 --012 --out data/chr7 5. 5.

Download RNA-seq read count data from the ReCount project[18]. At the bash prompt, type the following:

⬇

1curl -o data/montpick_count_table.txt link

2curl -o data/montpick_phenodata.txt link

3.2 Data importing

Create the following R scripts. At the bash prompt, type:44endnote: 4 Alternatively, you can create these files through rstudio (File>New File>R Script). :

⬇

1touch analysis.R R/{utils,viz}.R 2. 2.

Using the code editor provided by rstudio, add the line library(tidyverse) at the top of utils.R and viz.R, and the following lines at the top of analysis.R55endnote: 5Whenever the code in utils.R or viz.R changes, these files have to be sourced again to make the changes visible to analysis.R by executing lines 2 and 3 of this code chunk. A more structured approach would be to develop an R package, which is however beyond the scope of this chapter. The function doMC::registerDoMC(cores = 8) registers a parallel backend with 8 cores. This is required for multicore functionality in the plyr package.:

⬇

1source(’R/utils.R’)

2source(’R/viz.R’)

3

4doMC::registerDoMC(cores=8) 3. 3.

Create a function in utils.R for importing the genotypic data into R66endnote: 6 In this and all subsequent code listings, we make heavy use of the pipe operator %>% to form long linear chains of functions, where the output of each function becomes the input of the next. This allows for the generation of readable, easy-to-understand workflows. In this particular function, we read genotypic data from disk files in the form of an $M\times K$ matrix. Before returning, the matrix is filtered by removing all variants (i.e. columns) that may have zero variance across all samples (lines 24,25). Notice that file names are constructed in a device-independent manner using the function file.path. Finally, we prefer the operator = instead of the traditionally used <- for indicating assignment. :

⬇

1load_genotypes = function() {

2 # read sample names

3 samples =

4 file.path(’data’, ’chr7.012.indv’) %>%

5 readtsv(colnames = ’SAMPLE’) %>%

6 pull(SAMPLE)

7

8 # read loci

9 pos =

10 file.path(’data’, ’chr7.012.pos’) %>%

11 readtsv(colnames = c(’CHROM’, ’POS’)) %>%

12 unite(’POS’, CHROM, POS, sep = ’:’) %>%

13 pull(POS)

14

15 # read genotypes

16 geno =

17 file.path(’data’, ’chr7.012’) %>%

18 readtsv(na = ’-1’, colnames = c(’SAMPLE’, pos)) %>%

19 select(-1) %>%

20 as.matrix()

21 dimnames(geno) = list(samples = samples, variants = pos)

22

23 # remove loci with the same genotype across all samples

24 sds = apply(geno, 2, sd, na.rm = T)

25 geno[,sds > 0]

26} 4. 4.

Source utils.R, add the following line to analysis.R and execute it in order to import the genotypic data into R:

⬇

1obs_geno = load_genotypes() 5. 5.

Create the following function in viz.R for visualizing the correlation structure of the matrix of genotypes:

⬇

1plot_genotypes = function(geno) {

2 geno %>%

3 cor() %>%

4 reshape2::melt() %>%

5 ggplot() +

6 geom_raster(aes(x = Var1, y = Var2, fill = value)) +

7 theme(legend.position = ’none’,

8 axis.text.x = element_blank(),

9 axis.text.y = element_blank(),

10 axis.ticks.x = element_blank(),

11 axis.ticks.y = element_blank()) +

12 scale_fill_gradient(low = ’black’, high = ’white’) +

13 labs(x = ’genomic position’, y = ’genomic position’)

14} 6. 6.

Source viz.R, add the following line to analysis.R and execute it in order to visualize the genotypic data (Figure 1):

⬇

1plot_genotypes(obs_geno) 7. 7.

Create a function in utils.R for importing the read counts data into R77endnote: 7 Data are loaded in the form of an $N\times M$ matrix of read counts. By default, only data for the CEU population are read. If YRI data are preferred, set the function argument to pop=’YRI’, instead. As for the genotypic data, genes with constant expression across all samples are removed (lines 18,19). :

⬇

1load_counts = function(pop = ’CEU’) {

2 # load samples

3 samples =

4 read_delim(file.path(’data’, ’montpick_phenodata.txt’), delim = ’ ’) %>%

5 filter(population == pop) %>%

6 pull(sample.id)

7

8 # load read counts

9 counts =

10 read_tsv(file.path(’data’, ’montpick_count_table.txt’)) %>%

11 select(GENE = gene, samples) %>%

12 as.data.frame() %>%

13 column_to_rownames(’GENE’) %>%

14 as.matrix()

15 dimnames(counts) = list(genes = rownames(counts), samples = samples)

16

17 # remove genes with the same number of counts across all samples

18 sds = apply(counts, 1, sd, na.rm = T)

19 counts[sds > 0,]

20} 8. 8.

Source utils.R, add the following line to analysis.R and execute it in order to import the read counts data into R:

⬇

1obs_counts = load_counts()

3.3 Estimation of gene expression statistics

Create a function in utils.R for calculating sample-specific normalization factors given a matrix of count data88endnote: 8 This function implements the median-of-ratios method for normalizing a matrix of count data, as presented in [8]. Calculations are performed on the logarithmic scale. :

⬇

1calculate_norm_factors = function(counts) {

2 lcounts = log(counts)

3 lgm = rowSums(lcounts) / ncol(lcounts)

4 idxs = is.finite(lgm)

5 lratios = sweep(lcounts[idxs,], 1, lgm[idxs], ’-’)

6 apply(exp(lratios), 2, median)

7} 2. 2.

Create a function in utils.R for calculating the log-likelihood of a vector of counts, assuming each element is sampled from the Negative Binomial distribution99endnote: 9 The probability mass function of the Negative Binomial distribution assumed here is parameterized in terms of gene-specific mean $\mu_{i}$ and inverse dispersion $\alpha_{i}$ parameters, as follows:

$z_{ij}\sim\frac{\Gamma(z_{ij}+\alpha_{i})}{\Gamma(z_{ij}+1)\Gamma(\alpha_{i})}{\left(\frac{c_{j}\mu_{i}}{c_{j}\mu_{i}+\alpha_{i}}\right)}^{z_{ij}}{\left(\frac{\alpha_{i}}{c_{j}\mu_{i}+\alpha_{i}}\right)}^{\alpha_{i}}$

where $c_{j}$ are sample-specific normalisation factors. :

⬇

1lognbinom = function(pars, z, cc) {

2 mu = pars[1]

3 phi = pars[2]

4

5 m = cc * mu

6 alpha = 1 / phi

7 ll = lgamma(z + alpha) - lgamma(alpha) - lgamma(z + 1) +

8 z * log(m) + alpha * log(alpha) - (z + alpha) * log(m + alpha)

9

10 sum(ll)

11} 3. 3.

Create a function in utils.R for calculating gene-specific mean and dispersion maximum likelihood estimates for the Negative Binomial distribution given a matrix of count data1010endnote: 10 Given an $N\times M$ matrix of counts, we fit the Negative Binomial distribution to each gene/row, and we estimate gene-specific maximum likelihood estimates of mean and dispersion parameters. For parameter estimation, we use the function optim. The code possibly(optim, otherwise = NULL) on line 3 returns a version of the optim function, which returns NULL on error instead of halting execution, thus allowing processing of the whole set of available genes. The parameter fnscale=-1 on line 11 signals optim to perform maximization of the log-likelihood lognbinom, instead of minimization. Iterating over the rows of the data matrix is done using the plyr function alply, and it is performed in parallel (argument .parallel=T on line 15). After all rows have been processed, genes for which estimation failed or did not converge are discarded (lines 16,17), results are extracted and properly formatted (lines 18–21), and gene-specific variance values are calculated given the estimated values for mean and dispersion parameters (line 22). :

⬇

1calculate_gene_stats = function(counts) {

2 sizes = calculate_norm_factors(counts)

3 optim = possibly(optim, otherwise = NULL)

4 stats =

5 counts %>%

6 plyr::alply(1, function(cnts) {

7 optim(par = c(1, 1),

8 fn = lognbinom,

9 z = cnts,

10 cc = sizes,

11 control = list(fnscale = -1),

12 method = ’L-BFGS-B’,

13 hessian = T,

14 lower = 1e-12)

15 }, .parallel = T, .dims = T) %>%

16 compact() %>%

17 discard(~.x$convergence > 0) %>%

18 map(’par’) %>%

19 enframe() %>%

20 mutate(value = map(value, str_c, collapse = ’,’)) %>%

21 separate(value, into = c(’MEAN’, ’PHI’), sep = ’,’, convert = T) %>%

22 mutate(VAR = MEAN + PHI * MEAN^2) %>%

23 rename(genes = name) %>%

24 as.data.frame() %>%

25 column_to_rownames(’genes’) %>%

26 as.matrix()

27 dimnames(stats) = list(genes = rownames(stats), statistics = colnames(stats))

28 stats

29} 4. 4.

Source utils.R, add the following line to analysis.R and execute it in order to calculate basic gene-wise statistics from the read counts data:

⬇

1obs_counts_stats = calculate_gene_stats(obs_counts) 5. 5.

Create a function in viz.R for visualising empirical and estimated gene-wise statistics1111endnote: 11 An inner function plot_fn is defined at the top of plot_mean_variance, in order to avoid code repetition when generating plots gg1 and gg2 below. :

⬇

1plot_mean_variance = function(counts, stats) {

2 plot_fn = function(df, xlabel, ylabel) {

3 df %>%

4 ggplot() +

5 geom_point(aes(x = MEAN, y = VAR), size = 0.1) +

6 geom_smooth(aes(x = MEAN, y = VAR), size = 0.5,

7 color = ’red’, method = ’lm’, formula = y ~ poly(x, 2)) +

8 scale_x_continuous(trans = ’log10’) +

9 scale_y_continuous(trans = ’log10’) +

10 labs(x = xlabel, y = ylabel)

11 }

12

13 # normalise count data

14 sizes = calculate_norm_factors(counts)

15 counts = sweep(counts, 2, sizes, ’/’)

16

17 # plot of estimated stats

18 gg1 =

19 stats %>%

20 as.data.frame() %>%

21 plot_fn(xlabel = ’estimated mean’, ylabel = ’estimated variance’)

22

23 # plot of observed stats

24 gg2 =

25 data.frame(

26 MEAN = rowMeans(counts),

27 VAR = apply(counts, 1, var)

28 ) %>%

29 plot_fn(xlabel = ’observed mean’, ylabel = ’observed variance’)

30

31 # combine plots

32 cowplot::plot_grid(gg1, gg2, align = ’vh’, labels = ’AUTO’)

33} 6. 6.

Source viz.R, add the following line to analysis.R and execute it in order to visualise estimated and observed gene-wise statistics (Figure 2):

⬇

1plot_mean_variance(obs_counts, obs_counts_stats)

3.4 Data simulation

Create a function in utils.R for simulating genotypic and count gene expression data given observed data of the same kind and their estimated statistics1212endnote: 12 This function should be treated as a starting point for more sophisticated data simulation methods. Genotypic data are generated by randomly sampling without replacement nsamples rows and nvars columns from the observed matrix of genotypes. This approach breaks the spatial correlation between variants, so the variants/columns in the generated matrix $X$ can be treated as being mostly uncorrelated. The matrix $B$ of coefficients is initially generated as a matrix of zeros. In a second stage, nhits elements are randomly chosen, and set to non-zero values sampled randomly from a shifted exponential distribution. Half of the nhits elements are asigned positive values and the remaining half are assigned negative values. Given matrices $X$ and $B$ and the top ngenes estimated mean and dispersion parameters (where parameter values are first ranked in order of decreasing mean values; see line 20), a matrix $Z$ of read counts is generated with elements sampled from the Negative Binomial distribution. Alternatively, estimated mean and dispersion parameters can be randomly sampled by uncommenting line 19 (and commenting out line 20). :

⬇

1simulate_data = function(count_stats, genotypes, nsamples = 1000, ngenes = 100,

2 nvars = 50, nhits = 10, rate = 4) {

3 # fetch genotypes

4 X = genotypes[sample(1:nrow(genotypes), nsamples),

5 sample(1:ncol(genotypes), nvars)]

6 X = X[,apply(X, 2, sd, na.rm = T) > 0]

7 dimnames(X) = list(samples = str_c(’S’, 1:nrow(X)),

8 variants = str_c(’V’, 1:ncol(X)))

9

10 # simulate matrix of coefficients

11 B = matrix(0, nrow = ngenes, ncol = ncol(X),

12 dimnames = list(genes = str_c(’G’, 1:ngenes),

13 variants = str_c(’V’, 1:ncol(X))))

14 hits = c(1 + rexp(0.5 * nhits, rate = rate),

15 -1 - rexp(0.5 * nhits, rate = rate))

16 B[sample(length(B), length(hits))] = hits

17

18 # simulate counts

19 # df = count_stats[sample(1:nrow(count_stats), ngenes),]

20 df = count_stats[order(count_stats[,’MEAN’], decreasing = T),][1:ngenes,]

21

22 X0 = scale(X, center = T, scale = T)

23 m = sweep(exp(B %% t(X0)), 1, df[,’MEAN’], ’’)

24 alpha = 1 / df[,’PHI’]

25 Z = matrix(rnbinom(length(m), mu = m, size = alpha), nrow = nrow(m))

26 dimnames(Z) = list(genes = str_c(’G’, 1:nrow(Z)),

27 samples = str_c(’S’, 1:ncol(Z)))

28

29 # output

30 lst(B, X, Z, stats = df)

31} 2. 2.

Source utils.R to make simulate_data visible to the global environment.

3.5 Model implementation

Copy the files implementing the Normal, Poisson and Negative Binomial models to the stan/ directory (normal.stan, poissonln.stan and negbinom.stan, respectively; see end of chapter for code listings)1313endnote: 13 The stan code presented here is divided in a number of blocks. The data{…} block is read exactly once at the beginning of the inference procedure, and it declares required model data. The transformed data{…} block in negbinom.stan and poissonln.stan defines transformations of previously declared data, and it is also executed only once at initialisation. The parameters{…} block declares model parameters, an unconstrained version of which will be the target of inference procedures. The transformed parameters{…} block defines variables in terms of previously defined parameters. All variables that appear in these last two blocks will be returned at the output after inference completes. Finally, the model{…} block is where the joint log probability function is defined, using a syntax very close to the actual mathematical notation. Sampling statements (priors) for all variables in the parameters{…} block must appear here, otherwise a uniform, possibly unconstrained, prior is assumed. . 2. 2.

Compile the models by adding the following lines in analysis.R and executing them1414endnote: 14 The argument auto_write=T ensures that unecessary re-compilations of stan code are avoided. For this to work, the variables negbinom, poissonln and normal must reside in the global environment, as in this code listing. :

⬇

1normal =

2 file.path(’stan’, ’normal.stan’) %>%

3 rstan::stan_model(auto_write = T)

4poissonln =

5 file.path(’stan’, ’poissonln.stan’) %>%

6 rstan::stan_model(auto_write = T)

7negbinom =

8 file.path(’stan’, ’negbinom.stan’) %>%

9 rstan::stan_model(auto_write = T)

3.6 Model testing

Create an auxilliary function in utils.R for fitting a compiled model and extracting the estimated regression coefficients $B$1515endnote: 15 Stan provides three different inference methodologies[19]: (a) Full Bayesian inference using a self-adjusting Hamiltonian Monte Carlo approach, (b) approximate Bayesian inference using Automatic Differentiation Variational Inference (ADVI), and (c) point parameter estimates by maximizing the joint posterior. Here, for reasons of computational efficiency, we use the third approach (through function rstan::optimizing), which however does not provide any measure of uncertainty of the estimates. In principle, such estimates can be derived at a second stage by approximating locally the posterior distribution using the Laplace approximation. rstan::optimizing takes as input a compiled model and input data in the form of a list of variables, as they appear in the data{…} block of the model definition (variables not in this block are ignored). Function fit_model returns the matrices of true and estimated coefficients in vector format, and a vector of the corresponding indices. :

⬇

1fit_model = function(data, model, fcts = calculate_norm_factors(data$Z), ...) {

2 Z_tilde = sweep(data$Z, 2, fcts, ’/’) # normalise count data

3 Y = log(Z_tilde + 1) # transform normalised count data

4 X0 = scale(data$X, center = T, scale = T) # standardise genotypes

5

6 # MAP estimation

7 fit = rstan::optimizing(object = model,

8 data = list(Z = data$Z,

9 Y = Y,

10 X = X0,

11 c = fcts,

12 s = colSums(data$Z),

13 N = nrow(Y),

14 M = ncol(Y),

15 K = ncol(X0)),

16 seed = 42,

17 ...)

18

19 # extract estimated matrix of regression coefficients B in vector format

20 tibble(EST =

21 fit %>%

22 pluck(’par’) %>%

23 enframe() %>%

24 filter(str_detect(name, ’^B\[’)) %>%

25 pull(value),

26 TRU = as.vector(data$B),

27 IDX = 1:length(TRU))

28} 2. 2.

Source utils.R to make fit_model visible to the global environment. 3. 3.

Make the function cpp_object_initializer visible to the global environment by adding the following line to analysis.R and executing it1616endnote: 16 This is necessary for successfully calling rstan::optimizing. Alternatively, we could have imported the whole of rstan using library(rstan), and called optimizing directly (i.e. without the :: operator). :

⬇

1cpp_object_initializer = rstan::cpp_object_initializer 4. 4.

Set the initial seed of R’s random number generator by adding the following line in analysis.R and executing it1717endnote: 17 This code facilitates reproducibility, since for the same seed, the same sequence of pseudo-random numbers is generated. :

⬇

1set.seed(42) 5. 5.

Test all three models on simulated data of various sizes (312, 625, 1250, 2500) by adding the following lines in analysis.R and executing them1818endnote: 18 This code chunk includes two nested loops: the outer loop iterates over different sample sizes, while the inner loop iterates over models, and it is executed in parallel. At each iteration of the outer loop, a dataset of appropriate size is simulated and then passed to each tested model in the inner loop. :

⬇

1fitted_models =

2 plyr::ldply(lst(312, 625, 1250, 2500), function(N) {

3 sim = simulate_data(obs_counts_stats, obs_geno, nsamples = N)

4 plyr::ldply(lst(normal, poissonln, negbinom), function(mdl) {

5 fit_model(sim, mdl)

6 }, .parallel = T, .id = ’MODEL’)

7 }, .progress = ’text’, .id = ’NSAMPLES’) %>%

8 as_tibble() 6. 6.

Create a function in viz.R for visualising the results1919endnote: 19 This function takes a thr argument: estimated coefficient values below this threshold are considered effectively 0 and set to NA (line 4). When plotting, NA values will be dropped with a warning. :

⬇

1plot_fitted_models = function(df, thr = 0.1 * max(abs(df$EST))) {

2 df %>%

3 mutate(TRU = na_if(TRU, 0),

4 EST = ifelse(abs(EST) < thr, NAreal_, EST)) %>%

5 ggplot() +

6 geom_ribbon(aes(x = IDX, ymin = -thr, ymax = thr), fill = ’grey85’) +

7 geom_hline(yintercept = 0, linetype = ’dashed’, size = 0.2) +

8 geom_point(aes(x = IDX, y = TRU), color = ’red’) +

9 geom_point(aes(x = IDX, y = EST), size = 0.5) +

10 facet_grid(NSAMPLES~MODEL) +

11 scale_x_continuous(expand = c(0.02, 0.02)) +

12 scale_y_continuous(expand = c(0.02, 0.02)) +

13 labs(x = ’index of regression coefficients (genes x variants)’,

14 y = ’value of regression coefficients’)

15} 7. 7.

Source viz.R to make plot_fitted_models.R visible to the global environment. 8. 8.

Add the following line in analysis.R and execute it to visualise the estimated coefficients of the fitted models (Figure 3):

⬇

1plot_fitted_models(fitted_models)

4 Conclusion

Although it is not prudent to make general statements based on a single set of simulations, we may observe from Figure 3 that all three models perform similarly. At small samples (N=312), the Normal, Poisson-LogNormal and Negative Binomial models correctly identify 7, 8 and 7 eQTLs, respectively, but the Negative Binomial model throws in addition a false positive. None of these models estimates sufficiently well the size of the eQTLs they identify correctly. At higher samples, all models estimate sufficiently well the true eQTLs in the data, although they do not avoid the ocassional false positive (see Normal and Poisson-LogNormal models at N=1250 and the Negative Binomial model at N=2500). The Normal model is the fastest to estimate at all sample sizes, followed by the Poisson-LogNormal model. Proper parameterisation can greatly facilitate inference in a hierarchical model, and reduce its associated computational cost (see the Stan User’s Guide for more information).

Subsequently, we may wish to experiment with additional sparse-inducing priors, and/or more systematic and extensive benchmarks, which would require implementing a custom inference algorithm for increased performance. This initial analysis using a fully automated inference system provides a first assessment of the efficiency of each model, and it can serve as a baseline for the subsequent development of novel statistical methods.

5 Acknowledgments

This work was supported by the National Institute for Health Research (NIHR) Oxford Biomedical Research Centre Program. The views expressed in this manuscript are not necessarily those of the Department of Health. We also acknowledge the Wellcome Trust Centre for Human Genetics Wellcome Trust Core Award Grant Number 090532/Z/09/Z. The funders of the study had no role in the study design, data collection, data analysis, data interpretation, or writing of the paper.

Notes

1. 1 You can access the command prompt through a terminal emulator. From within rstudio, you can create a terminal emulator by going to Tools>Terminal>New Terminal.
2. 2 Instead of steps 2 and 3, you can just create a new rstudio project (File>New Project…) with equivalent results
3. 3 This chunk of code uses curl to stream data from the 1000 Genomes Project and pipe them (in gzipped VCF format) to vcftools for further processing. Only variants in chromosome 7 between positions 100K and 200K are retained. In addition, all variants with FILTER other than PASS are removed; indels are removed; variants in non-bi-allelic sites are removed; only variants with allele frequency higher than 5% are retained (i.e. 0.05<MAF<0.95). Filtered data are stored as a matrix with elements 0, 1 and 2 indicating the number of alleles in 496 variants across 2504 samples. If vcftools throws an error about not being able to open temporary files, you need to increase the maximum number of open files allowed by your system to more than double the number of samples in the VCF file. At the bash command prompt, type: ulimit -n 5100. You may need root permissions to run this command.
4. 4 Alternatively, you can create these files through rstudio (File>New File>R Script).
5. 5 Whenever the code in utils.R or viz.R changes, these files have to be sourced again to make the changes visible to analysis.R by executing lines 2 and 3 of this code chunk. A more structured approach would be to develop an R package, which is however beyond the scope of this chapter. The function doMC::registerDoMC(cores = 8) registers a parallel backend with 8 cores. This is required for multicore functionality in the plyr package.
6. 6 In this and all subsequent code listings, we make heavy use of the pipe operator %>% to form long linear chains of functions, where the output of each function becomes the input of the next. This allows for the generation of readable, easy-to-understand workflows. In this particular function, we read genotypic data from disk files in the form of an $M\times K$ matrix. Before returning, the matrix is filtered by removing all variants (i.e. columns) that may have zero variance across all samples (lines 24,25). Notice that file names are constructed in a device-independent manner using the function file.path. Finally, we prefer the operator = instead of the traditionally used <- for indicating assignment.
7. 7 Data are loaded in the form of an $N\times M$ matrix of read counts. By default, only data for the CEU population are read. If YRI data are preferred, set the function argument to pop=’YRI’, instead. As for the genotypic data, genes with constant expression across all samples are removed (lines 18,19).
8. 8 This function implements the median-of-ratios method for normalizing a matrix of count data, as presented in [8]. Calculations are performed on the logarithmic scale.
9. 9 The probability mass function of the Negative Binomial distribution assumed here is parameterized in terms of gene-specific mean $\mu_{i}$ and inverse dispersion $\alpha_{i}$ parameters, as follows:

$z_{ij}\sim\frac{\Gamma(z_{ij}+\alpha_{i})}{\Gamma(z_{ij}+1)\Gamma(\alpha_{i})}{\left(\frac{c_{j}\mu_{i}}{c_{j}\mu_{i}+\alpha_{i}}\right)}^{z_{ij}}{\left(\frac{\alpha_{i}}{c_{j}\mu_{i}+\alpha_{i}}\right)}^{\alpha_{i}}$

where $c_{j}$ are sample-specific normalisation factors. 10. 10 Given an $N\times M$ matrix of counts, we fit the Negative Binomial distribution to each gene/row, and we estimate gene-specific maximum likelihood estimates of mean and dispersion parameters. For parameter estimation, we use the function optim. The code possibly(optim, otherwise = NULL) on line 3 returns a version of the optim function, which returns NULL on error instead of halting execution, thus allowing processing of the whole set of available genes. The parameter fnscale=-1 on line 11 signals optim to perform maximization of the log-likelihood lognbinom, instead of minimization. Iterating over the rows of the data matrix is done using the plyr function alply, and it is performed in parallel (argument .parallel=T on line 15). After all rows have been processed, genes for which estimation failed or did not converge are discarded (lines 16,17), results are extracted and properly formatted (lines 18–21), and gene-specific variance values are calculated given the estimated values for mean and dispersion parameters (line 22). 11. 11 An inner function plot_fn is defined at the top of plot_mean_variance, in order to avoid code repetition when generating plots gg1 and gg2 below. 12. 12 This function should be treated as a starting point for more sophisticated data simulation methods. Genotypic data are generated by randomly sampling without replacement nsamples rows and nvars columns from the observed matrix of genotypes. This approach breaks the spatial correlation between variants, so the variants/columns in the generated matrix $X$ can be treated as being mostly uncorrelated. The matrix $B$ of coefficients is initially generated as a matrix of zeros. In a second stage, nhits elements are randomly chosen, and set to non-zero values sampled randomly from a shifted exponential distribution. Half of the nhits elements are asigned positive values and the remaining half are assigned negative values. Given matrices $X$ and $B$ and the top ngenes estimated mean and dispersion parameters (where parameter values are first ranked in order of decreasing mean values; see line 20), a matrix $Z$ of read counts is generated with elements sampled from the Negative Binomial distribution. Alternatively, estimated mean and dispersion parameters can be randomly sampled by uncommenting line 19 (and commenting out line 20). 13. 13 The stan code presented here is divided in a number of blocks. The data{…} block is read exactly once at the beginning of the inference procedure, and it declares required model data. The transformed data{…} block in negbinom.stan and poissonln.stan defines transformations of previously declared data, and it is also executed only once at initialisation. The parameters{…} block declares model parameters, an unconstrained version of which will be the target of inference procedures. The transformed parameters{…} block defines variables in terms of previously defined parameters. All variables that appear in these last two blocks will be returned at the output after inference completes. Finally, the model{…} block is where the joint log probability function is defined, using a syntax very close to the actual mathematical notation. Sampling statements (priors) for all variables in the parameters{…} block must appear here, otherwise a uniform, possibly unconstrained, prior is assumed. 14. 14 The argument auto_write=T ensures that unecessary re-compilations of stan code are avoided. For this to work, the variables negbinom, poissonln and normal must reside in the global environment, as in this code listing. 15. 15 Stan provides three different inference methodologies[19]: (a) Full Bayesian inference using a self-adjusting Hamiltonian Monte Carlo approach, (b) approximate Bayesian inference using Automatic Differentiation Variational Inference (ADVI), and (c) point parameter estimates by maximizing the joint posterior. Here, for reasons of computational efficiency, we use the third approach (through function rstan::optimizing), which however does not provide any measure of uncertainty of the estimates. In principle, such estimates can be derived at a second stage by approximating locally the posterior distribution using the Laplace approximation. rstan::optimizing takes as input a compiled model and input data in the form of a list of variables, as they appear in the data{…} block of the model definition (variables not in this block are ignored). Function fit_model returns the matrices of true and estimated coefficients in vector format, and a vector of the corresponding indices. 16. 16 This is necessary for successfully calling rstan::optimizing. Alternatively, we could have imported the whole of rstan using library(rstan), and called optimizing directly (i.e. without the :: operator). 17. 17 This code facilitates reproducibility, since for the same seed, the same sequence of pseudo-random numbers is generated. 18. 18 This code chunk includes two nested loops: the outer loop iterates over different sample sizes, while the inner loop iterates over models, and it is executed in parallel. At each iteration of the outer loop, a dataset of appropriate size is simulated and then passed to each tested model in the inner loop. 19. 19 This function takes a thr argument: estimated coefficient values below this threshold are considered effectively 0 and set to NA (line 4). When plotting, NA values will be dropped with a warning.

Appendix A File normal.stan

⬇

1data {

2 int<lower=0> N; // number of genes

3 int<lower=0> M; // number of samples

4 int<lower=0> K; // number of variants

5 matrix[M, K] X; // matrix of genotypes

6 vector[M] Y[N]; // matrix of transformed read counts

7}

8

9parameters {

10 vector[N] b0; // baseline gene expression

11 vector[N] ls2; // log-variance of gene expression

12 real<lower=0> eta; // global scale parameter

13 vector<lower=0>[K] zeta[N]; // local scale parameters

14 vector[K] B[N]; // regression coefficients

15}

16

17transformed parameters {

18 vector<lower=0>[N] sigma = sqrt(exp(ls2)); // standard deviation of gene expression

19}

20

21model {

22 real sc = mean(sigma) / sqrt(N*K);

23 eta ~ cauchy(0, 1);

24 for(i in 1:N) {

25 zeta[i] ~ cauchy(0, 1);

26 B[i] ~ normal(0, eta * zeta[i] * sc);

27 Y[i] ~ normal(b0[i] + X * B[i], sigma[i]);

28 }

29}

Appendix B File poissonln.stan

⬇

1data {

2 int<lower=0> N; // number of genes

3 int<lower=0> M; // number of samples

4 int<lower=0> K; // number of variants

5 matrix[M, K] X; // matrix of genotypes

6 int<lower=0> Z[N, M]; // matrix of read counts

7 vector<lower=0>[M] c; // vector of normalisation factors

8 vector<lower=0>[M] s; // vector of library sizes

9}

10

11transformed data {

12 vector[M] lc = log(c); // log normalisation factors

13 vector[M] ls = log(s); // log library sizes

14}

15

16parameters {

17 vector[N] b0; // baseline gene expression

18 vector[N] ls2; // log-variance of gene expression

19 real<lower=0> eta; // global scale parameter

20 vector<lower=0>[K] zeta[N]; // local scale parameters

21 vector[K] B[N]; // regression coefficients

22 vector[M] Y[N]; // latent variables

23}

24

25transformed parameters {

26 vector<lower=0>[N] sigma = sqrt(exp(ls2)); // standard deviation of gene expression

27}

28

29model {

30 real sc = mean(sigma) / sqrt(N*K);

31 eta ~ cauchy(0, 1);

32 for(i in 1:N) {

33 zeta[i] ~ cauchy(0, 1);

34 B[i] ~ normal(0, eta * zeta[i] * sc);

35 Y[i] ~ normal(b0[i] + X * B[i], sigma[i]);

36 Z[i] ~ poisson_log(lc + ls + Y[i]);

37 }

38}

Appendix C File negbinom.stan

⬇

1data {

2 int<lower=0> N; // number of genes

3 int<lower=0> M; // number of samples

4 int<lower=0> K; // number of variants

5 matrix[M, K] X; // matrix of genotypes

6 int<lower=0> Z[N, M]; // matrix of read counts

7 vector<lower=0>[M] c; // vector of normalisation factors

8 vector<lower=0>[M] s; // vector of library sizes

9}

10

11transformed data {

12 vector[M] lc = log(c); // log normalisation factors

13 vector[M] ls = log(s); // log library sizes

14}

15

16parameters {

17 vector[N] b0; // baseline gene expression (log-scale)

18 vector[N] lphi; // log-dispersion

19 real<lower=0> eta; // global scale parameter

20 vector<lower=0>[K] zeta[N]; // local scale parameters

21 vector[K] B[N]; // regression coefficients

22}

23

24transformed parameters {

25 vector<lower=0>[N] phi = exp(lphi); // dispersion

26 vector<lower=0>[N] alpha = 1.0 ./ phi; // inverse-dispersion

27}

28

29model {

30 real sc = 1.0 / sqrt(N*K);

31 eta ~ cauchy(0, 1);

32 for(i in 1:N) {

33 zeta[i] ~ cauchy(0, 1);

34 B[i] ~ normal(0, eta * zeta[i] * sc);

35 Z[i] ~ neg_binomial_2_log(ls + lc + b0[i] + X * B[i], alpha[i]);

36 }

37}

Appendix D R sessionInfo()

• •

R version 3.5.1 (2018–07–02), x86_64-pc-linux-gnu

• •

Locale: LC_CTYPE=en_GB.UTF-8, LC_NUMERIC=C, LC_TIME=en_GB.UTF-8, LC_COLLATE=en_GB.UTF-8, LC_MONETARY=en_GB.UTF-8, LC_MESSAGES=en_GB.UTF-8, LC_PAPER=en_GB.UTF-8, LC_NAME=C, LC_ADDRESS=C, LC_TELEPHONE=C, LC_MEASUREMENT=en_GB.UTF-8, LC_IDENTIFICATION=C

• •

Running under: Ubuntu 18.04.1 LTS

• •

Matrix products: default

• •

BLAS: /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3

• •

LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so

• •

Base packages: base, datasets, graphics, grDevices, methods, stats, utils

• •

Other packages: bindrcpp 0.2.2, dplyr 0.7.8, forcats 0.3.0, ggplot2 3.1.0, purrr 0.2.5, readr 1.2.1, stringr 1.3.1, tibble 1.4.2, tidyr 0.8.2, tidyverse 1.2.1

• •

Loaded via a namespace (and not attached): assertthat 0.2.0, backports 1.1.2, base64enc 0.1–3, bindr 0.1.1, broom 0.5.1, callr 3.0.0, cellranger 1.1.0, cli 1.0.1, codetools 0.2–15, colorspace 1.3–2, compiler 3.5.1, cowplot 0.9.3, crayon 1.3.4, doMC 1.3.5, foreach 1.4.4, generics 0.0.2, glue 1.3.0, grid 3.5.1, gridExtra 2.3, gtable 0.2.0, haven 2.0.0, hms 0.4.2, httr 1.3.1, inline 0.3.15, iterators 1.0.10, jsonlite 1.5, labeling 0.3, lattice 0.20–38, lazyeval 0.2.1, loo 2.0.0, lubridate 1.7.4, magrittr 1.5, matrixStats 0.54.0, modelr 0.1.2, munsell 0.5.0, nlme 3.1–137, parallel 3.5.1, pillar 1.3.0, pkgbuild 1.0.2, pkgconfig 2.0.2, plyr 1.8.4, prettyunits 1.0.2, processx 3.2.1, ps 1.2.1, R6 2.3.0, Rcpp 1.0.0, readxl 1.1.0, reshape2 1.4.3, rlang 0.3.0.1, rstan 2.18.2, rstudioapi 0.8, rvest 0.3.2, scales 1.0.0, StanHeaders 2.18.0, stats4 3.5.1, stringi 1.2.4, tidyselect 0.2.5, tools 3.5.1, withr 2.1.2, xml2 1.2.0, yaml 2.2.0