In variani/bigcov: Covariance/Correlation for big data

opts_chunk$set(comment = NA, results = 'markup', tidy = F, message = F, warning = F, echo = T,
  fig.width = 6, fig.height = 6)

Preliminaries

library(utils) # for downloads
library(BEDMatrix) # to read plink data efficiently

library(RSpectra) # to perform PCA with a few components, e.g. 2

library(magrittr) 
library(dplyr) 

#library(bigcov)
library(devtools)
load_all("~/git/variani/bigcov/")

Function to plot PCA

plot_pop <- function(mod, labs, ...)
{
  plot(mod$vectors[, 1], mod$vectors[, 2], type = "n", 
    xlab = "PC1", ylab = "PC2", ...)
  text(mod$vectors[, 1], mod$vectors[, 2], label = labs, col = as.numeric(labs))
}

Download and extract data

http://zzz.bwh.harvard.edu/plink/res.shtml#hapmap

The HapMap genotype data (the latest is release 23) are available here as PLINK binary filesets. The SNPs are currently coded according NCBI build 36 coordinates on the forward strand. Several versions are available here: the entire dataset (a single, very large fileset: you will need a computer with at least 2Gb of RAM to load this file).

The filtered SNP set refers to a list of SNPs that have MAF greater than 0.01 and genotyping rate greater than 0.95 in the 60 CEU founders. This fileset is probably a good starting place for imputation in samples of European descent. Filtered versions of the other HapMap panels will be made available shortly.

dir_data <- "data-hapmap" 
dir.create(dir_data, showWarnings = FALSE)

# description file
url_pop <- "http://zzz.bwh.harvard.edu/plink/dist/hapmap.pop"
fn_pop <- basename(url_pop)
download.file(url_pop, file.path(dir_data, fn_pop))

# genotype file
url_gen <- "http://zzz.bwh.harvard.edu/plink/dist/hapmap_JPT_CHB_r23a.zip"
fn_gen <- basename(url_gen)
download.file(url_gen, file.path(dir_data, fn_gen))

unzip(file.path(dir_data, fn_gen), exdir = dir_data)

list.files(dir_data)

Read files by BEDMatrix

bed_gen <- paste0(strsplit(fn_gen, '[.]')[[1]][1], ".bed")
system.time({
  bmat <- BEDMatrix(file.path(dir_data, bed_gen))
})

bmat

dim(bmat)

Extract IDs of individuals

ids <- bmat@dnames[[1]]
head(ids)

ids <- strsplit(ids, "_") %>% sapply(function(x) x[1])
head(ids)

pop <- read.table(file.path(dir_data, fn_pop), header = FALSE)
names(pop) <- c("name", "id", "pop")

pop <- left_join(data_frame(id = ids), pop, by = "id")
stopifnot(nrow(pop) == length(ids))

Summary on data set

| Characteristic | Value | |--|--| | File | r fn_gen | | Description | JPT+CHB (release 23, 90 individuals, 3.99 million SNPs) | | Orignial file | r bed_gen | | File format in R | r class(bmat) | | Number of individuals | r nrow(bmat) | | Number of markers | r ncol(bmat) %>% format(big.mark = ",") |

GRM on all markers

system.time({
  grm <- bigdat_grm(bmat, batch_size = 1e5, verbose = 1)
})

mod <- eigs(grm, k = 2)

labs <- pop$pop
plot_pop(mod, labs, main = "GRM on all markers")

GRM on common markers (> 5%)

system.time({
  grm <- bigdat_grm(bmat, batch_size = 1e5, maf_min = 0.05, verbose = 1)
})

mod <- eigs(grm, k = 2)

labs <- pop$pop
plot_pop(mod, labs, main = "GRM on common markers")

Jacard on rare markers (<5%)

The genotype data does not have enough rare variants, e.g. < 1%. Thus, the Jacard matrix is the identity matrix if the analysis is run on genotypes with MAF below 1%.

system.time({
  grm <- bigdat_grm(bmat, "Jacard", batch_size = 1e5, maf_max = 0.05, verbose = 1) 
})

mod <- eigs(grm, k = 2)

labs <- pop$pop
plot_pop(mod, labs, main = "Jacard on rare markers")

GRM on rare markers (<5%)

system.time({
  grm <- bigdat_grm(bmat, batch_size = 1e5, maf_max = 0.05, verbose = 1) 
})

mod <- eigs(grm, k = 2)

labs <- pop$pop
plot_pop(mod, labs, main = "GRM on rare markers")

Parallel computing

We compute the previous example on a single or two CPUs. The following code does not work on Mac with BLAS supporting multi threading.

system.time(bigdat_grm(bmat, batch_size = 1e5, maf_max = 0.05))

system.time(bigdat_grm(bmat, batch_size = 1e5, maf_max = 0.05, cores = 2))

R session info

sessionInfo()

License

This document is licensed under the Creative Commons Attribution 4.0 International Public License.

References

variani/bigcov documentation built on May 3, 2019, 4:34 p.m.