Analysis description

n <- 40000 # number of samples
p <- 656321 # number of SNPs
J <- 8021 # number of genes

The genotype data we used is from UKB biobank, randomly selecting n = 40000 samples. We use SNP genotype data from chr 1 to chr 22 combined from UKB. SNPs are downsampled to 1/10 (randomly), eQTLs (see below for definition of eQTL) were added back. This ends up with p = p as.charater(p) SNPs.

Our analysis consists of the following steps:

1. Build expression predictors using another expression-genotype dataset.

The one we used in this analysis is GTEx Adipose tissue v7 dataset. This dataset contains ~ 380 samples. FUSION/TWAS were used to train expression model and we used their lasso results. SNPs included in eQTL anlaysis are restricted to cis-locus 500kb on either side of the gene boundary. eQTLs are defined as SNPs with abs(effectize) > 1e-8 in lasso results.

2. Impute expression.

We impute gene expression for our genotype data using expression models obtained from step 1. There are 8021 genes with expression model from chr17 to chr22. We imputed expression from genotypes using the expression predictors.

3. Define and select regions

Next, the analysis is done at the “region” level, which is 500kb bins along the genome. Each bin would contain all the SNPs, as well as all the genes in that bin. We are exploring several ways to select regions that contain true signals, e.g. based on regional sum of mr.ash PIP for genes/SNPs, region smallest TWAS p value for gene/SNPs, or regional bayes factors, etc.

4. Run susie iteratively We then run susie for each of these regions. So the features of SuSiE are: SNPs and “genes” (not cis-eQTLs of that gene). We use the same prior for all SNPs and another prior for all “genes” when running SUSIE. In some settings, we also run SUSIE with null weight, which is calculated as 1- prior.SNP * n.SNP - prior.gene * n.gene. We obtain the PIP for SNPs and gene in the region. After we run susie for all regions (one iteration), we take the average of all SNP PIPs as the prior of SNPs for the next iteration and similarly for the prior for genes.

5. We obtain PIP for genes from the last iteration as results.

simdatadir <- "~/causalTWAS/simulations/simulation_ashtest_20201001/"
outputdir <- "~/causalTWAS/simulations/simulation_susieI_20201001/"
susiedir <- "~/causalTWAS/simulations/simulation_susieI_20201001/"
tag <- '20201001-1-1'
Niter <- 20

p value distribution

• TWAS p value of genes:

simtag <- "20201001-1-1"
exprgwasf <- paste0(simdatadir, simtag, ".exprgwas.txt.gz")
load(paste0(simdatadir, "simu_", simtag, "-pheno.Rd"))
caulist <- list()
for (chrom in 1:22) {
  load(paste0("~/causalTWAS/ukbiobank/ukb_chr", chrom ,"_s40000.FBM.Rd"))
  load(paste0(simdatadir, "simu_s40000_GTEXadipose-B", chrom, "-cis-expr.Rd"))
  caulist[[chrom]]<- c(exprres$gnames[phenores$batch[[chrom]]$param$idx.cgene], dat$snp[phenores$batch[[chrom]]$param$idx.cSNP,])
}
cau <- unlist(caulist)

a <- read.table(exprgwasf, header = T)
a$ifcausal <- ifelse(a$MARKER_ID %in% cau, 1, 0)
ax <- pretty(0:max(-log10(a$PVALUE)), n = 30)

par(mfrow=c(3,1))
h1 <- hist(-log10(a$PVALUE), breaks = 100, xlab = "-log10(p)", main = "P value distribution-all", col = "grey"); grid()
h2 <- hist(-log10(a[a$ifcausal == 1, ]$PVALUE), breaks = h1$breaks, xlab = "-log10(p)", main = "P value distribution-causal", col = "salmon");grid()

plot(a[a$X.CHROM ==22, ]$BEGIN, -log10(a[a$X.CHROM ==22, ]$PVALUE), col = a[a$X.CHROM ==22, ]$ifcausal + 1, xlab = paste0("chr", chrom), ylab = "-log10(pvalue)")
points(a[a$X.CHROM ==22 & a$ifcausal ==1, ]$BEGIN, -log10(a[a$X.CHROM ==22 & a$ifcausal ==1, ]$PVALUE), col = "red", pch =19)

r <- do.call(cbind, lapply(1:22, function(x) phenores$batch[[x]]$Y.g))
GWAA(exprres$expr, matrix(rowSums(r[, c(1:12,14:21)]) + rnorm(N), ncol=1), snpname = exprres$gnames, anno = anno, "temp", family = gaussian, ncore = 1, nSplits = 1, compress = F)
m <- read.table("temp", header =T, comment.char = "")
m$ifcausal <- ifelse(m$MARKER_ID %in% cau, 1, 0)
plot(m[m$X.CHROM ==22, ]$BEGIN, -log10(m[m$X.CHROM ==22, ]$PVALUE), col = m$ifcausal + 1)


snpname = dat$snp[,1]
anno <- cbind(dat$chr, dat$pos, dat$pos) # EPACT format
  
GWAA(dat$G, r[,13, drop=F], snpname = snpname, anno = anno, "tempsnp", family = gaussian, ncore = 5, nSplits = 100, compress = F)

• TWAS p value of SNPs:

simtag <- "20201001-1-1"
snpgwasf <- paste0(simdatadir, simtag, ".snpgwas.txt.gz")

b <- read.table(snpgwasf, header = T)
b$ifcausal <- ifelse(b$MARKER_ID %in% cau, 1, 0)
ax <- pretty(0:max(-log10(b$PVALUE)), n = 30)

par(mfrow=c(3,1))
h1 <- hist(-log10(b$PVALUE), breaks = 100, xlab = "-log10(p)", main = "P value distribution-all", col = "grey"); grid()
h2 <- hist(-log10(b[b$ifcausal == 1, ]$PVALUE), breaks = h1$breaks, xlab = "-log10(p)", main = "P value distribution-causal", col = "salmon");grid()

plot(b[b$X.CHROM ==22, ]$BEGIN, -log10(b[b$X.CHROM ==22, ]$PVALUE), col = b$ifcausal + 1, xlab = paste0("chr", chrom), ylab = "-log10(pvalue)")
points(b[b$X.CHROM ==22 & b$ifcausal ==1, ]$BEGIN, -log10(b[b$X.CHROM ==22 & b$ifcausal ==1, ]$PVALUE), col = "red", pch =19)

Parameter estimation results

Results: Each row shows parameter estimation results from 5 simulation runs with similar settings (i.e. pi1 and PVE for genes and SNPs). each row has two plots, one for gene pi1 estimation, one for enrichment (gene pi1/snp pi1). Results from each run were represented by one dot, dots with the same color come from the same run. horizontal dash lines: simulation truth, susietruth, the truth in selected regions that were used to run susie iteractively (susieI).

show_param <- function(phenofs, susieIfs, susieIfs2){
  pars <- do.call(rbind, lapply(phenofs, function(x) {load(x); 
    c(phenores$param$pve.gene.truth,
      phenores$param$pve.snp.truth,
      length(phenores$batch[[1]]$param$idx.cgene)/phenores$batch[[1]]$param$J,
      length(phenores$batch[[1]]$param$idx.cSNP)/phenores$batch[[1]]$param$M)}))
 
  colnames(pars) <- c("PVE.gene_truth", "PVE.SNP_truth", "pi1.gene_truth", "pi1.SNP_truth")
    
  param.s <- do.call(rbind, lapply(susieIfs, function(x) {load(x); c(tail(prior.gene_rec[prior.gene_rec!=0], 1), tail(prior.SNP_rec[prior.SNP_rec!=0],1))}))
  
  param.s.truth <- do.call(rbind, lapply(susieIfs2, function(x) {
    a <- fread(x, header = T);
    c(nrow(a[a$ifcausal == 1 & a$type == "gene" ])/ nrow(a[a$type == "gene"]),
      nrow(a[a$ifcausal == 1 & a$type == "SNP"])/ nrow(a[a$type == "SNP"]))
    }))
  
  pars.s <- cbind(param.s.truth, param.s)[, c(1,3,2,4)]
  colnames(pars.s) <-  paste(rep(c("pi1.gene_", "pi1.SNP_"), each = 2), c("susietruth", "susieI"), sep = "")
  
  df <- cbind(tags, format(pars, digits = 4), format(pars.s, digits =4))
  rownames(df) <- NULL
  return(df)
  
  # df %>% 
  # kable("html", escape = F) %>%
  # kable_styling("striped", full_width = F) %>%
  #  row_spec(c(1:5, 11:15), background = "#FEF3B9") %>%
  # scroll_box(width = "100%", height = "600px", fixed_thead = T)
}

plot_param <- function(df, ...){
  df <- apply(df[ , 2:ncol(df)], 2, function(x) as.numeric(x))
  st <-  cbind(df[,"pi1.gene_susietruth"], 1:nrow(df), 2 + 1:nrow(df)/nrow(df)/3)
  s <- cbind(df[,"pi1.gene_susieI"], 1:nrow(df), 3 + 1:nrow(df)/nrow(df)/3)
  t <- df[1,"pi1.gene_truth"]
  dfp <- rbind(st,s)
  plot(dfp[,3], dfp[,1], col = dfp[,2], pch = 19, ylab = "gene pi1", xaxt = "n", xlab="", xlim = c(0.8, 3.5), frame.plot=FALSE, ylim = c(0, max(dfp[,1],t) *1.05), ...)
  axis(side=1, at=1:2, labels = FALSE, tick = F)
  text(x=2:3, 0, labels = c( "susieI_truth", "susieI"), xpd = T, pos =1)
  abline(h=t, lty = 2, col= "salmon", lwd=1.5)
  grid()
  
  st <-  cbind(df[,"pi1.gene_susietruth"]/df[,"pi1.SNP_susietruth"], 1:nrow(df), 2 + 1:nrow(df)/nrow(df)/3)
  s <- cbind(df[,"pi1.gene_susieI"]/df[,"pi1.SNP_susieI"], 1:nrow(df), 3 + 1:nrow(df)/nrow(df)/3)
  t <- df[1,"pi1.gene_truth"]/df[1,"pi1.SNP_truth"]
  dfp <- rbind(st,s)
  plot(dfp[,3], dfp[,1], col = dfp[,2], pch = 19, ylab = "Enrichment (gene/snp)", xaxt = "n", xlab="", xlim = c(0.8, 3.5),frame.plot=FALSE, ylim = c(0, min(max(dfp[,1],t) *1.05, 150)))
  axis(side=1, at=1:2, labels = FALSE, tick = F)
  text(x=2:3, 0, labels = c("susieI_truth", "susieI"), xpd = T, pos =1)
  abline(h= t, lty = 2, col= "darkgreen", lwd=1.5)
  grid()
}

gpip_dist <- function(susiefs, ...){
  dflist <- list()
  for (f in susiefs){
    dflist[[f]] <- read.table(f, header =T , stringsAsFactors = F)
  }
  df <- do.call(rbind, dflist)
  hist(df[df$type == "gene", "susie_pip"], xlab = "gene susie PIP", 
       breaks = 50, ylim = c(0,20), xlim=c(0,1), col = "salmon", ...)
}

tags <- paste(rep(1:2, each = 5), 1:5, sep = "-")
phenofs <- paste0(simdatadir, "simu_20201001-", tags, "-pheno.Rd")
susieIfs <- paste0(outputdir, "20201001-", tags, ".config1.susieIres.Rd")
susieIfs2 <- paste0(outputdir, "20201001-", tags, ".config1.susieI.txt")

df <- show_param(phenofs, susieIfs, susieIfs2)
par(mfrow = c(2,2))
plot_param(df[1:5,], main = "low power")
plot_param(df[6:10,], main = "high power")

tags <- paste(rep(1:2, each = 5), 1:5, sep = "-")
phenofs <- paste0(simdatadir, "simu_20201001-", tags, "-pheno.Rd")
susieIfs <- paste0(outputdir, "20201001-", tags, ".config5.susieIres.Rd")
susieIfs2 <- paste0(outputdir, "20201001-", tags, ".config5.susieI.txt")

df <- show_param(phenofs, susieIfs, susieIfs2)
par(mfrow = c(2,2))
plot_param(df[1:5,], main = "low power")
plot_param(df[6:10,], main = "high power")

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)

Matrix products: default
BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [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       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] kableExtra_1.2.1    stringr_1.4.0       plyr_1.8.6         
[4] tidyr_0.8.3         plotly_4.9.2.9000   ggplot2_3.3.1      
[7] data.table_1.12.7   mr.ash.alpha_0.1-34

loaded via a namespace (and not attached):
 [1] tidyselect_1.1.0   purrr_0.3.4        lattice_0.20-38   
 [4] bigassertr_0.1.3   colorspace_1.3-2   vctrs_0.3.1       
 [7] generics_0.0.2     htmltools_0.3.6    viridisLite_0.3.0 
[10] yaml_2.2.0         rlang_0.4.6        later_0.7.5       
[13] pillar_1.4.4       glue_1.4.1         withr_2.1.2       
[16] foreach_1.4.4      lifecycle_0.2.0    munsell_0.5.0     
[19] gtable_0.2.0       workflowr_1.6.2    rvest_0.3.2       
[22] htmlwidgets_1.3    codetools_0.2-15   evaluate_0.12     
[25] knitr_1.20         doParallel_1.0.15  httpuv_1.4.5      
[28] parallel_3.5.1     Rcpp_1.0.4.6       promises_1.0.1    
[31] scales_1.0.0       backports_1.1.2    webshot_0.5.1     
[34] jsonlite_1.6.1     fs_1.3.1           digest_0.6.25     
[37] bigparallelr_0.2.3 stringi_1.3.1      dplyr_1.0.0       
[40] cowplot_0.9.4      grid_3.5.1         rprojroot_1.3-2   
[43] here_0.1           tools_3.5.1        magrittr_1.5      
[46] lazyeval_0.2.1     tibble_3.0.1       crayon_1.3.4      
[49] pkgconfig_2.0.2    ellipsis_0.3.1     Matrix_1.2-15     
[52] xml2_1.2.0         bigstatsr_1.2.3    iterators_1.0.10  
[55] rmarkdown_1.10     httr_1.4.1         rstudioapi_0.11   
[58] R6_2.3.0           flock_0.7          git2r_0.26.1      
[61] compiler_3.5.1