Contents

1 Introduction

This document explores using Hail 0.2 with R via basilisk.

The computations follow the GWAS tutorial in the hail documentation. We won’t do all the computations there, and we add some material dealing with R-python interfacing. We’ll note that the actual computations on large data are done in Spark, but we don’t interact directly with Spark at all in this document.

Most of the computations are done via reticulate calls to python; the access to the hail environment is through basilisk. We also take advantage of R markdown’s capacity to execute python code directly. If an R chunk computes x, a python chunk can refer to it as r.x. If a python chunk computes r.x, an R chunk can refer to this value as x.

2 Installing BiocHail

BiocHail should be installed as follows:

if (!require("BiocManager"))
    install.packages("BiocManager")
BiocManager::install("BiocHail")

As of 1.0.0, a JDK for Java version <= 11.0 is necessary to use the version of Hail that is installed with the package. This package should be usable on MacOS with suitable java support. If Java version >= 8.x is used, warnings from Apache Spark may be observed. To the best of our knowledge the conditions to which the warnings pertain do not affect program performance.

3 Acquire a slice of the 1000 genomes genotypes and annotations

In this section we import the 1000 genomes VCF slice distributed by the hail project. hail_init uses basilisk, which ensures that a specific version of hail and its dependencies are available in an isolated virtual environment.

library(BiocHail)
library(ggplot2)

3.1 Initialization, data acquisition, rendering

Here is a curiosity of R-hail interaction. Note that the following chunk computes mt, a MatrixTable representation of 1000 genomes data, but our attempt to print it in markdown fails.

hl <- hail_init()
## + /home/biocbuild/.cache/R/basilisk/1.12.1/0/bin/conda 'create' '--yes' '--prefix' '/home/biocbuild/.cache/R/basilisk/1.12.1/BiocHail/1.0.1/bsklenv' 'python=3.8.13' '--quiet' '-c' 'conda-forge'
## + /home/biocbuild/.cache/R/basilisk/1.12.1/0/bin/conda 'install' '--yes' '--prefix' '/home/biocbuild/.cache/R/basilisk/1.12.1/BiocHail/1.0.1/bsklenv' 'python=3.8.13'
## + /home/biocbuild/.cache/R/basilisk/1.12.1/0/bin/conda 'install' '--yes' '--prefix' '/home/biocbuild/.cache/R/basilisk/1.12.1/BiocHail/1.0.1/bsklenv' '-c' 'conda-forge' 'python=3.8.13' 'pandas=1.3.5'
mt <- get_1kg(hl)
mt
## <hail.matrixtable.MatrixTable object at 0x7f18bb466310>
print(mt$rows()$select()$show(5L)) # limited info
## +---------------+------------+
## | locus         | alleles    |
## +---------------+------------+
## | locus<GRCh37> | array<str> |
## +---------------+------------+
## | 1:904165      | ["G","A"]  |
## | 1:909917      | ["G","A"]  |
## | 1:986963      | ["C","T"]  |
## | 1:1563691     | ["T","G"]  |
## | 1:1707740     | ["T","G"]  |
## +---------------+------------+
## showing top 5 rows
## 
## NULL

We can use the python syntax in a python R markdown chunk to see what we want. We use prefix r. to find references defined in our R session (compiling the vignette).

r.mt.rows().select().show(5) # python chunk!
## +---------------+------------+
## | locus         | alleles    |
## +---------------+------------+
## | locus<GRCh37> | array<str> |
## +---------------+------------+
## | 1:904165      | ["G","A"]  |
## | 1:909917      | ["G","A"]  |
## | 1:986963      | ["C","T"]  |
## | 1:1563691     | ["T","G"]  |
## | 1:1707740     | ["T","G"]  |
## +---------------+------------+
## showing top 5 rows

The sample IDs:

r.mt.s.show(5)  # python chunk!
## +-----------+
## | s         |
## +-----------+
## | str       |
## +-----------+
## | "HG00096" |
## | "HG00099" |
## | "HG00105" |
## | "HG00118" |
## | "HG00129" |
## +-----------+
## showing top 5 rows

3.2 Helper functions

Some methods return data immediately useful in R.

mt$count()
## [[1]]
## [1] 10961
## 
## [[2]]
## [1] 284

We can thus define a function dim to behave with hail MatrixTable instances in a familiar way, along with some others.

dim.hail.matrixtable.MatrixTable <- function(x) { 
  tmp <- x$count()
  c(tmp[[1]], tmp[[2]]) 
}
dim(mt)
## [1] 10961   284
ncol.hail.matrixtable.MatrixTable <- function(x) { 
 dim(x)[2]
}
nrow.hail.matrixtable.MatrixTable <- function(x) { 
 dim(x)[1]
}
nrow(mt)
## [1] 10961

These can be useful on their own, or when calling python methods.

3.3 Acquiring column fields

AS OF 21 AUGUST 2023 THE CODE BELOW SEGFAULTS IN VIGNETTE PROCESSING BUT NOT WHEN RUN IN REPL. MOST OF THE CODE HAS MIGRATED TO A TEST SCRIPT under tests/testthat.

annopath <- path_1kg_annotations()
tab <- hl$import_table(annopath, impute=TRUE)$key_by("Sample")
r.tab.describe()# python chunk!
r.tab.show(width=100)

3.4 Adding the sample annotation to the MatrixTable; aggregation

We combine the tab defined above, with the MatrixTable instance, using python code reaching to R via r..

r.mt = r.mt.annotate_cols(pheno = r.tab[r.mt.s])   # python chunk! 
r.mt.col.describe()

Aggregation methods can be used to obtain contingency tables or descriptive statistics.

First, we get the frequencies of superpopulation membership:

mt$aggregate_cols(hl$agg$counter(mt$pheno$SuperPopulation))

Then statistics on caffeine consumption:

uu <- mt$aggregate_cols(hl$agg$stats(mt$pheno$CaffeineConsumption))
names(uu)
uu$mean
uu$stdev

4 Working with variants; quality assessment

The significance of the aggregation functions is that the computations are performed by Spark, on potentially huge distributed data structures.

Now we aggregate over rows (SNPs). We’ll use python directly:

from pprint import pprint  # python chunk!
snp_counts = r.mt.aggregate_rows(r.hl.agg.counter(r.hl.Struct(ref=r.mt.alleles[0], alt=r.mt.alleles[1])))
pprint(snp_counts)

4.1 A histogram of read depths

Hail uses the concept of ‘entries’ for matrix elements, and each ‘entry’ is a ‘struct’ with potentially many fields.

Here we’ll use R to compute a histogram of sequencing depths over all samples and variants.

p_hist <- mt$aggregate_entries(
     hl$expr$aggregators$hist(mt$DP, 0L, 30L, 30L))
names(p_hist)
length(p_hist$bin_edges)
length(p_hist$bin_freq)
midpts <- function(x) diff(x)/2+x[-length(x)]
dpdf <- data.frame(x=midpts(p_hist$bin_edges), y=p_hist$bin_freq)
ggplot(dpdf, aes(x=x,y=y)) + geom_col() + ggtitle("DP") + ylab("Frequency")

An exercise: produce a function mt_hist that produces a histogram of measures from any of the relevant VCF components of a MatrixTable instance.

Note also all the aggregators available:

names(hl$expr$aggregators)

We’d also note that hail has a direct interface to ggplot2.

4.2 Quality summaries

A high-level function adds quality metrics to the MatrixTable.

mt <- hl$sample_qc(mt)
r.mt.col.describe()  # python!

The call rate histogram is given by:

callrate <- mt$sample_qc$call_rate$collect()
graphics::hist(callrate)

4.3 Filtering

4.3.1 Sample quality

We’ll use the python code given for filtering, in which per-sample mean read depth and call rate are must exceed (arbitrarily chosen) thresholds.

# python chunk!
r.mt = r.mt.filter_cols((r.mt.sample_qc.dp_stats.mean >= 4) & (r.mt.sample_qc.call_rate >= 0.97))
print('After filter, %d/284 samples remain.' % r.mt.count_cols())

4.3.2 Genotype quality

Again we use the python code for filtering according to

  • relative purity of reads underlying homozygous reference or alt calls
  • good balance of ref/alt counts for het calls
ab = r.mt.AD[1] / r.hl.sum(r.mt.AD)

filter_condition_ab = ((r.mt.GT.is_hom_ref() & (ab <= 0.1)) |
                        (r.mt.GT.is_het() & (ab >= 0.25) & (ab <= 0.75)) |
                        (r.mt.GT.is_hom_var() & (ab >= 0.9)))

fraction_filtered = r.mt.aggregate_entries(r.hl.agg.fraction(~filter_condition_ab))
print(f'Filtering {fraction_filtered * 100:.2f}% entries out of downstream analysis.')
r.mt = r.mt.filter_entries(filter_condition_ab)

Note that filtering entries does not change MatrixTable shape.

dim(mt)

4.3.3 Variant characteristics

Allele frequencies, Hardy-Weinberg equilibrium test results and other summaries are obtained using the variant_qc function.

mt = hl$variant_qc(mt)
r.mt.row.describe()  #! python

5 GWAS execution

A built-in procedure for testing for association between the (simulated) caffeine consumption measure and genotype will be used.

The following commands eliminate variants with minor allele frequency less than 0.01, along with those with small \(p\)-values in tests of Hardy-Weinberg equilibrium.

r.mt = r.mt.filter_rows(r.mt.variant_qc.AF[1] > 0.01)
r.mt = r.mt.filter_rows(r.mt.variant_qc.p_value_hwe > 1e-6)
r.mt.count()

5.1 Association test for quantitative response

Now we perform a naive test of association. The Manhattan plot generated by hail can be displayed for interaction using bokeh. We comment this out for now; it should be possible to embed the bokeh display in this document but the details are not ready-to-hand.

r.gwas = r.hl.linear_regression_rows(y=r.mt.pheno.CaffeineConsumption,
                                 x=r.mt.GT.n_alt_alleles(),
                                 covariates=[1.0])
# r.pl = r.hl.plot.manhattan(r.gwas.p_value)
# import bokeh
# bokeh.plotting.show(r.pl)

The “QQ plot” that helps evaluate adequacy of the analysis can be formed using hl.plot.qq for very large applications; here we collect the results for plotting in R.

First we estimate \(\lambda_{GC}\)

pv = gwas$p_value$collect()
x2 = stats::qchisq(1-pv,1)
lam = stats::median(x2, na.rm=TRUE)/stats::qchisq(.5,1)
lam

And the qqplot:

qqplot(-log10(ppoints(length(pv))), -log10(pv), xlim=c(0,8), ylim=c(0,8),
  ylab="-log10 p", xlab="expected")
abline(0,1)

There is hardly any point to examining a Manhattan plot in this situation. But let’s see how it might be done. We’ll use igvR to get an interactive and extensible display.

locs <- gwas$locus$collect()
conts <- sapply(locs, function(x) x$contig)
pos <- sapply(locs, function(x) x$position)
library(igvR)
mytab <- data.frame(chr=as.character(conts), pos=pos, pval=pv)
gt <- GWASTrack("simp", mytab, chrom.col=1, pos.col=2, pval.col=3)
igv <- igvR()
setGenome(igv, "hg19")
displayTrack(igv, gt)

5.2 Evaluating population stratification

An approach to assessing population stratification is provided as hwe_normalized_pca. See the hail methods docs for details.

We are avoiding a tuple assignment in the tutorial document.

r.pcastuff = r.hl.hwe_normalized_pca(r.mt.GT)
r.mt = r.mt.annotate_cols(scores=r.pcastuff[1][r.mt.s].scores)

We’ll collect the key information and plot.

sc <- pcastuff[[2]]$scores$collect()
pc1 = sapply(sc, "[", 1)
pc2 = sapply(sc, "[", 2)
fac = mt$pheno$SuperPopulation$collect()
myd = data.frame(pc1, pc2, pop=fac)
library(ggplot2)
ggplot(myd, aes(x=pc1, y=pc2, colour=factor(pop))) + geom_point()

Now repeat the association test with adjustments for population of origin and gender.

r.gwas2 = r.hl.linear_regression_rows(
    y=r.mt.pheno.CaffeineConsumption,
    x=r.mt.GT.n_alt_alleles(),
    covariates=[1.0,r.mt.pheno.isFemale,r.mt.scores[0],
        r.mt.scores[1], r.mt.scores[2]])

New value of \(\lambda_{GC}\):

pv = gwas2$p_value$collect()
x2 = stats::qchisq(1-pv,1)
lam = stats::median(x2, na.rm=TRUE)/stats::qchisq(.5,1)
lam

A manhattan plot for chr8:

locs <- gwas2$locus$collect()
conts <- sapply(locs, function(x) x$contig)
pos <- sapply(locs, function(x) x$position)
mytab <- data.frame(chr=as.character(conts), pos=pos, pval=pv)
ggplot(mytab[mytab$chr=="8",], aes(x=pos, y=-log10(pval))) + geom_point()

6 Conclusions

The tutorial document proceeds with some illustrations of arbitrary aggregations. We will skip these for now.

Additional vignettes will address

SessionInfo

sessionInfo()
## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.17-bioc/R/lib/libRblas.so 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB              LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: America/New_York
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] ggplot2_3.4.3    BiocHail_1.0.1   BiocStyle_2.28.0
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.7            utf8_1.2.3            generics_0.1.3       
##  [4] RSQLite_2.3.1         lattice_0.21-8        digest_0.6.33        
##  [7] magrittr_2.0.3        evaluate_0.21         grid_4.3.1           
## [10] bookdown_0.35         fastmap_1.1.1         blob_1.2.4           
## [13] jsonlite_1.8.7        Matrix_1.6-1          DBI_1.1.3            
## [16] BiocManager_1.30.22   httr_1.4.7            purrr_1.0.2          
## [19] fansi_1.0.4           scales_1.2.1          jquerylib_0.1.4      
## [22] cli_3.6.1             rlang_1.1.1           dbplyr_2.3.3         
## [25] basilisk.utils_1.12.1 munsell_0.5.0         bit64_4.0.5          
## [28] withr_2.5.0           cachem_1.0.8          yaml_2.3.7           
## [31] tools_4.3.1           dir.expiry_1.8.0      parallel_4.3.1       
## [34] memoise_2.0.1         dplyr_1.1.2           colorspace_2.1-0     
## [37] filelock_1.0.2        basilisk_1.12.1       BiocGenerics_0.46.0  
## [40] curl_5.0.2            reticulate_1.31       png_0.1-8            
## [43] vctrs_0.6.3           R6_2.5.1              BiocFileCache_2.8.0  
## [46] lifecycle_1.0.3       bit_4.0.5             pkgconfig_2.0.3      
## [49] gtable_0.3.4          pillar_1.9.0          bslib_0.5.1          
## [52] Rcpp_1.0.11           glue_1.6.2            xfun_0.40            
## [55] tibble_3.2.1          tidyselect_1.2.0      knitr_1.43           
## [58] htmltools_0.5.6       rmarkdown_2.24        compiler_4.3.1