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.
BiocHailBiocHail should be installed as follows:
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.
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.
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.
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).
The sample IDs:
Some methods return data immediately useful in R.
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)
ncol.hail.matrixtable.MatrixTable <- function(x) {
dim(x)[2]
}
nrow.hail.matrixtable.MatrixTable <- function(x) {
dim(x)[1]
}
nrow(mt)These can be useful on their own, or when calling python methods.
column fieldsWe combine the tab defined above, with the MatrixTable instance,
using python code reaching to R via r..
Aggregation methods can be used to obtain contingency tables or descriptive statistics.
First, we get the frequencies of superpopulation membership:
Then statistics on caffeine consumption:
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)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:
We’d also note that hail has a direct interface to ggplot2.
A high-level function adds quality metrics to the MatrixTable.
The call rate histogram is given by:
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.
Again we use the python code for filtering according to
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.
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()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)
lamAnd 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)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)
lamA manhattan plot for chr8:
The tutorial document proceeds with some illustrations of arbitrary aggregations. We will skip these for now.
Additional vignettes will address
## R Under development (unstable) (2024-10-21 r87258)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.21-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.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.5.1 BiocHail_1.7.1 BiocStyle_2.35.0
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.9 utf8_1.2.4 generics_0.1.3
## [4] RSQLite_2.3.8 lattice_0.22-6 digest_0.6.37
## [7] magrittr_2.0.3 evaluate_1.0.1 grid_4.5.0
## [10] bookdown_0.41 fastmap_1.2.0 blob_1.2.4
## [13] jsonlite_1.8.9 Matrix_1.7-1 DBI_1.2.3
## [16] BiocManager_1.30.25 httr_1.4.7 fansi_1.0.6
## [19] scales_1.3.0 jquerylib_0.1.4 cli_3.6.3
## [22] rlang_1.1.4 dbplyr_2.5.0 basilisk.utils_1.19.0
## [25] munsell_0.5.1 bit64_4.5.2 withr_3.0.2
## [28] cachem_1.1.0 yaml_2.3.10 tools_4.5.0
## [31] dir.expiry_1.15.0 parallel_4.5.0 memoise_2.0.1
## [34] dplyr_1.1.4 colorspace_2.1-1 filelock_1.0.3
## [37] basilisk_1.19.0 BiocGenerics_0.53.3 curl_6.0.1
## [40] reticulate_1.40.0 png_0.1-8 vctrs_0.6.5
## [43] R6_2.5.1 BiocFileCache_2.15.0 lifecycle_1.0.4
## [46] bit_4.5.0 pkgconfig_2.0.3 gtable_0.3.6
## [49] pillar_1.9.0 bslib_0.8.0 Rcpp_1.0.13-1
## [52] glue_1.8.0 xfun_0.49 tibble_3.2.1
## [55] tidyselect_1.2.1 knitr_1.49 htmltools_0.5.8.1
## [58] rmarkdown_2.29 compiler_4.5.0