Last updated: 2020-12-15

Checks: 5 2

Knit directory: causal-TWAS/

This reproducible R Markdown analysis was created with workflowr (version 1.6.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

R Markdown file: uncommitted changes

The R Markdown file has unstaged changes. To know which version of the R Markdown file created these results, you’ll want to first commit it to the Git repo. If you’re still working on the analysis, you can ignore this warning. When you’re finished, you can run wflow_publish to commit the R Markdown file and build the HTML.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20191103)

The command set.seed(20191103) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: absolute

Using absolute paths to the files within your workflowr project makes it difficult for you and others to run your code on a different machine. Change the absolute path(s) below to the suggested relative path(s) to make your code more reproducible.

absolute relative
~/causalTWAS/causal-TWAS/code/gwas.R code/gwas.R
~/causalTWAS/causal-TWAS/code/ld.R code/ld.R
~/causalTWAS/causal-TWAS/code/qqplot.R code/qqplot.R

Repository version: dfb918e

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version dfb918e. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated: 


Ignored files:
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    code/workflow/.ipynb_checkpoints/
    Ignored:    data/

Untracked files:
    Untracked:  analysis/figure/

Unstaged changes:
    Modified:   analysis/simulation-susieI-ukbWG-gtex.adipose_s40.22.Rmd
    Modified:   analysis/simulation-susieI-ukbWG-gtex.adipose_visual.Rmd
    Modified:   code/run_ld.R

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/simulation-susieI-ukbWG-gtex.adipose_s40.22.Rmd) and HTML (docs/simulation-susieI-ukbWG-gtex.adipose_s40.22.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File Version Author Date Message
Rmd dfb918e simingz 2020-12-11 susieI-rss
html dfb918e simingz 2020-12-11 susieI-rss
Rmd 6e78356 simingz 2020-12-03 indididual level pipeline
html 6e78356 simingz 2020-12-03 indididual level pipeline
Rmd 07891d9 simingz 2020-11-25 filter and rerun
html 07891d9 simingz 2020-11-25 filter and rerun
Rmd 41a7f93 simingz 2020-11-24 credible set
html 41a7f93 simingz 2020-11-24 credible set
Rmd e81ebe1 simingz 2020-11-20 block feature
html e81ebe1 simingz 2020-11-20 block feature
Rmd 6fa5cfa simingz 2020-11-19 susie different prior variance
html 6fa5cfa simingz 2020-11-19 susie different prior variance
Rmd fd9b304 simingz 2020-11-18 index
Rmd 983dea2 simingz 2020-11-14 more PIP cali
html 983dea2 simingz 2020-11-14 more PIP cali
Rmd 447a401 simingz 2020-11-13 pip cali
html 447a401 simingz 2020-11-13 pip cali
Rmd 2e5d53d simingz 2020-11-04 SER
html 2e5d53d simingz 2020-11-04 SER
Rmd f057bb5 simingz 2020-11-03 SERI
html f057bb5 simingz 2020-11-03 SERI
Rmd f9c55de simingz 2020-10-30 log bayes factor
html f9c55de simingz 2020-10-30 log bayes factor
Rmd f7cd12a simingz 2020-10-24 qqplot
html f7cd12a simingz 2020-10-24 qqplot
Rmd ba58a99 simingz 2020-10-24 qqplot
html ba58a99 simingz 2020-10-24 qqplot
Rmd bf57f32 simingz 2020-10-24 qqplot
html bf57f32 simingz 2020-10-24 qqplot
Rmd 6526634 simingz 2020-10-23 filtered samples s40.22
html 6526634 simingz 2020-10-23 filtered samples s40.22 
library(mr.ash.alpha)
library(data.table)
suppressMessages({library(plotly)})
library(tidyr)
library(plyr)
library(stringr)
source("analysis/summarize_twas_plots.R")

Analysis description

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

The genotype data we used is from UKB biobank, randomly selecting 40000 samples. We then filtered samples based on relatedness, ethics and other qc metrics, that ended up with n = 22542 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 = 656321 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.

  1. Impute expression.

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

  1. Define and select regions

Next, the analysis is done at the “region” level, which is 500kb bins along the genome. We also used LDetect to define regions. 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.

  1. 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.

  2. 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/"
simtag <- '20201001-1-3'
source('~/causalTWAS/causal-TWAS/code/gwas.R')
source('~/causalTWAS/causal-TWAS/code/ld.R')
source('~/causalTWAS/causal-TWAS/code/qqplot.R')
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 ,"_s40.22.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)

Power estimation

We use gene.pve ~ 0.1, snp.pve ~ 0.5.

  • For SNPs, we use \(\pi_1 = 2.5e-3\), variance for effect size ~ \(0.03^2\), power at 5e-8 p value cutoff:
pow <- function(total, n, beta, cutp){
  rec <- rep(0, total)
  for (i in 1:total){
    x <- rnorm(n)
    y <- x * rnorm(1, sd = beta) + rnorm(n, sd = sqrt(2.5))
    lm.s <- lm(y~x)
    pv <- summary(lm.s)$coefficients[2,4]
    rec[i] <- pv
  }
  length(rec[rec < cutp])/length(rec)
}
load("data/power_s40.22.Rd")
total <- 1e3
n <- 22542
#p1 <- pow(total, n, 0.0276, 5e-8)
print(p1)
[1] 0.056
  • For genes, under low power setting, \(\pi_1 = 0.05\), variance for effect size ~ \(0.025^2\), power at 1e-5 cutoff:
#p2 <- pow(total, n, 0.025, 1e-5)
print(p2)
[1] 0.079

For genes, under high power setting, \(\pi_1 = 0.02\), variance for effect size ~ \(0.045^2\), power at 1e-5 cutoff:

#p3 <- pow(total, n, 0.045, 1e-5)
print(p3)
[1] 0.317
#save(p1,p2,p3, file = "data/power_s40.22.Rd")

p value distribution

  • TWAS p value of genes:
chrom <- 1
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", xlim= c(3,20), ylim =c(0,50)); grid()
h2 <- hist(-log10(a[a$ifcausal == 1, ]$PVALUE), breaks = h1$breaks, xlab = "-log10(p)", main = "P value distribution-causal", col = "salmon", xlim= c(3,20), ylim =c(0,50));grid()

cat("number of genes p < 1e-5:", nrow(a[a$PVALUE < 1e-5, ]))
number of genes p < 1e-5: 67
cat("number of causal genes p < 1e-5:", nrow(a[a$PVALUE < 1e-5 & a$ifcausal ==1, ]))
number of causal genes p < 1e-5: 33
plot(a[a$X.CHROM ==chrom, ]$BEGIN, -log10(a[a$X.CHROM ==chrom, ]$PVALUE), col = a[a$X.CHROM ==chrom, ]$ifcausal + 1, xlab = paste0("chr", chrom), ylab = "-log10(pvalue)")
points(a[a$X.CHROM ==chrom & a$ifcausal ==1, ]$BEGIN, -log10(a[a$X.CHROM ==chrom & a$ifcausal ==1, ]$PVALUE), col = "red", pch =19)
grid()

Version Author Date
f9c55de simingz 2020-10-30
bf57f32 simingz 2020-10-24
6526634 simingz 2020-10-23
  • qq plot for genes
gg_qqplot(a$PVALUE) +
  theme_bw(base_size = 24) +
  theme(
    axis.ticks = element_line(size = 0.5),
    panel.grid = element_blank()
    # panel.grid = element_line(size = 0.5, color = "grey80")
  )

Version Author Date
f9c55de simingz 2020-10-30
bf57f32 simingz 2020-10-24
  • p value of SNPs (GWAS):
snpgwasf <- paste0(simdatadir, simtag, ".snpgwas.txt.gz")

b <- fread(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", xlim= c(3,20), ylim =c(0,100)); grid()
h2 <- hist(-log10(b[b$ifcausal == 1, ]$PVALUE), breaks = h1$breaks, xlab = "-log10(p)", main = "P value distribution-causal", col = "salmon", xlim= c(3,20), ylim =c(0,100));grid()

cat("number of SNPs < 5e-8: ", nrow(b[b$PVALUE < 5e-8,]))
number of SNPs < 5e-8:  876
plot(b[b$X.CHROM ==chrom, ]$BEGIN, -log10(b[b$X.CHROM ==chrom, ]$PVALUE), col = b[b$X.CHROM ==chrom,]$ifcausal + 1, xlab = paste0("chr", chrom), ylab = "-log10(pvalue)")
points(b[b$X.CHROM ==chrom & b$ifcausal ==1, ]$BEGIN, -log10(b[b$X.CHROM ==chrom & b$ifcausal ==1, ]$PVALUE), col = "red", pch =19)
grid()

Version Author Date
f9c55de simingz 2020-10-30
bf57f32 simingz 2020-10-24
6526634 simingz 2020-10-23
  • qq plot for SNPs
gg_qqplot(b$PVALUE) +
  theme_bw(base_size = 24) +
  theme(
    axis.ticks = element_line(size = 0.5),
    panel.grid = element_blank()
    # panel.grid = element_line(size = 0.5, color = "grey80")
  )

Version Author Date
f9c55de simingz 2020-10-30
bf57f32 simingz 2020-10-24

Block features

a <- read.table(paste0(outputdir, "20201001-1-1.config9.gene.nofilter.r.txt"), header =T)
b <-  read.table(paste0(outputdir, "20201001-2-1.config9.gene.nofilter.r.txt"), header =T)
m1 <- read.table(paste0(outputdir, "20201001-1-1.config9.nofilter.r.txt"), header =T)
m2 <- read.table(paste0(outputdir, "20201001-2-1.config9.nofilter.r.txt"), header =T)
  
hist(a$p1-a$p0, main = NULL, xlab = "Block size(bp)", col = "salmon")

Version Author Date
e81ebe1 simingz 2020-11-20 
par(mfrow=c(2,2))
hist(a$nCausal, xlab = "No.causal genes", col = "salmon", breaks=100, xlim=c(0,15), main ="low power")
hist(b$nCausal, xlab = "No.causal genes", col = "salmon", breaks=100, xlim=c(0,15), main ="high power")
hist(m1$nCausal, xlab = "No.causal genes + SNPs", col = "salmon", breaks=100, xlim=c(0,15), main ="low power")
hist(m2$nCausal, xlab = "No.causal genes + SNPs", col = "salmon", breaks=100, xlim=c(0,15), main ="high power")

Version Author Date
e81ebe1 simingz 2020-11-20

TWAS FDP

library(ggplot2)
library(cowplot)
library(plotrix)

.obn <- function(pips, ifcausal, mode = c("PIP", "FDR")){
  a_bin <- cut(pips, breaks= seq(0, 1, by=0.1))
  if (mode == "PIP") {
    ob = c(by(ifcausal, a_bin, FUN = sum))
  } else if (mode == "FDR"){
    ob = c(by((1-ifcausal), a_bin, FUN = sum))
  }
  return(ob)
} 

nca_plot <- function(pips, ifcausal, runtag = NULL, mode = c("PIP", "FDR"), main = mode[1], ...){
    
  # ifcausal:0,1, runtag: for adding std.
  if (is.null(runtag)){
     se = 0 
  } else{
    dflist <- list()
    for (rt in unique(runtag)){
      pips.rt <- pips[runtag == rt]
      ifcausal.rt <- ifcausal[runtag == rt]
      dflist[[rt]] <- .obn(pips.rt, ifcausal.rt, mode = mode)
    }
    nca_mean <- colMeans(do.call(rbind, dflist))
    se <- apply(do.call(rbind, dflist), 2, plotrix::std.error)
  }
  
  df <- data.frame("ncausal" = nca_mean, "se" = se)
  
  fig <- ggplot(df) +
    geom_bar( aes(x=seq(0, 1, by=0.1)[1:10]+0.05, y=ncausal), color ="black", stat="identity", fill="salmon", alpha=0.7, width = 0.1) +
    geom_errorbar( aes(x= seq(0, 1, by=0.1)[1:10] + 0.05, ymin= ncausal-se, ymax=ncausal+se), width=0.05, colour="black", alpha=0.9, size=0.5) +
    ggtitle(main) +
    xlab("PIP") + ylab("No. causal genes")
    theme_cowplot()
  return(fig)
}

ncausal_plot <- function(pipfs, format = "susie", main = "PIP"){
  df <- NULL
  for (i in 1:length(pipfs)) {
    res <- fread(pipfs[i], header = T)
    res <- data.frame(res[res$type  =="gene", ])
    res$runtag <- i
    if (format == "susie"){
       res <- rename(res, c("susie_pip" = "pip") )
    } else if (format == "SER"){
       res <- rename(res, c("SERpip" = "pip"), )
    }
    res <- res[complete.cases(res),]
    df <- rbind(df, res)
  }

  fig <- nca_plot(df$pip, df$ifcausal, df$runtag, mode ="PIP", main = main)
  return(fig)
}

.exob <- function(pips, ifcausal, mode = c("PIP", "FDR")){
  a_bin <- cut(pips, breaks= seq(0, 1, by=0.1))
  if (mode == "PIP") {
    ex = c(by(pips, a_bin, FUN = mean))
    ob = c(by(ifcausal, a_bin, FUN = mean))
  } else if (mode == "FDR"){
    ex = c(by(pips, a_bin, FUN = mean))
    ob = 1 - c(by(ifcausal, a_bin, FUN = mean))
  }
  return(list("expected" = ex, "observed" = ob))
}


dot_plot = function(df, main) {
        ggplot(df, aes(x=mean_pip, y=observed_freq)) + 
          geom_errorbar(aes(ymin=observed_freq-se, ymax=observed_freq+se), colour="black", size = 0.5, width=.01) +
          geom_point(size=1.5, shape=21, fill="#002b36") + # 21 is filled circle
          xlab("Expected") +
          ylab("Observed") +
          coord_cartesian(ylim=c(0,1), xlim=c(0,1)) +
          geom_abline(slope=1,intercept=0,colour='red', size=0.2) +
          ggtitle(main) +
          expand_limits(y=0) +                        # Expand y range
          theme_cowplot()}

cp_plot <- function(pips, ifcausal, runtag = NULL, mode = c("PIP", "FDR"), main = mode[1], ...){
    
  # ifcausal:0,1, runtag: for adding std.
  if (is.null(runtag)){
     se = 0 
  } else{
    dflist <- list()
    for (rt in unique(runtag)){
      pips.rt <- pips[runtag == rt]
      ifcausal.rt <- ifcausal[runtag == rt]
      dflist[[rt]] <- .exob(pips.rt, ifcausal.rt, mode = mode)
    }
    mean_pip <- colMeans(do.call(rbind, lapply(dflist, '[[', "expected")))
    observed_freq <- colMeans(do.call(rbind, lapply(dflist, '[[', "observed")))
    se <- apply(do.call(rbind, lapply(dflist, '[[', "observed")), 2, plotrix::std.error)
  }
  
  df <- data.frame("mean_pip" = mean_pip, "observed_freq"= observed_freq, "se" = se)
  dot_plot(df, main = main)

  #plot(Expected, Observed, xlim= c(0,1), ylim=c(0,1), pch =19, main = main, ...)
  #lines(x = c(0,1), y = c(0,1), col ="grey", lty = 2)
}

caliPIP_plot <- function(pipfs, format = "susie", main = "PIP"){
  df <- NULL
  for (i in 1:length(pipfs)) {
    res <- fread(pipfs[i], header = T)
    res <- data.frame(res[res$type  =="gene", ])
    res$runtag <- i
    if (format == "susie"){
       res <- rename(res, c("susie_pip" = "pip"))
    } else if (format == "SER"){
       res <- rename(res, c("SERpip" = "pip"))
    }
    res <- res[complete.cases(res),]
    df <- rbind(df, res)
  }

  fig <- cp_plot(df$pip, df$ifcausal, df$runtag, mode ="PIP", main = main)
  return(fig)
}

# pipfs is just for ifcausal info
caliFDP_plot <- function(gwasfs, pipfs,  format= "susie",  main = "FDP"){
  df <- NULL
  for (i in 1:length(pipfs)) {
    pipres <- fread(pipfs[i], header = T)
    pipres <- data.frame(pipres[pipres$type  =="gene", ])
    pipres$runtag <- i
    if (format == "susie"){
       pipres <- rename(pipres, c("susie_pip" = "pip"))
    } else if (format == "SER"){
       pipres <- rename(pipres, c("SERpip" = "pip"))
    }
    gwasres <- read.table(gwasfs[i], header = T)
    gwasres <- rename(gwasres, c("MARKER_ID" = "name"))
    gwasres$FDR <- p.adjust(gwasres$PVALUE, method = "fdr")
    res <- merge(gwasres, pipres, by = "name", all = T)
    res <- res[complete.cases(res),]
    df <- rbind(df, res)
  }
  fig <- cp_plot(df$FDR, df$ifcausal, df$runtag, mode ="FDR", main = main)
  cat("FDP at bonferroni corrected p = 0.05: ", 1 - mean(df[df$PVALUE < 0.05 /dim(df)[1], "ifcausal"]))
  return(fig)
}


scatter_plot_PIP_p <- function(pipfs, gwasfs, pipformat = "susie", main ="PIP-p"){
  df <- NULL
  for (i in 1:length(pipfs)) {
    pipres <- fread(pipfs[i], header = T)
    pipres <- data.frame(pipres[pipres$type  =="gene", ])
    pipres$runtag <- i
    if (pipformat == "susie"){
       pipres <- rename(pipres, c("susie_pip" = "pip"))
    } else if (pipformat == "SER"){
       pipres <- rename(pipres, c("SERpip" = "pip"))
    }
    gwasres <- read.table(gwasfs[i], header = T)
    gwasres <- rename(gwasres, c("MARKER_ID" = "name"))
    res <- merge(gwasres, pipres, by = "name", all = T)
    res <- res[complete.cases(res),]
    df <- rbind(df, res)
  }
  
  df <- rename(df, c( "PVALUE" = "TWAS.p"))
  df[,"TWAS.p"] <- -log10(df[, "TWAS.p"])
  
  
                            
  df$ifcausal <- mapvalues(df$ifcausal, from=c(0,1),to=c("darkgreen", "salmon"))
  plot(df$TWAS.p, df$pip, col = df$ifcausal, main = main, xlab = "-log10(TWAS p value)", ylab = "PIP")

  # df$ifcausal <- mapvalues(df$ifcausal, from=c(0,1), to=c("Non causal", "Causal"))
  # fig <- plot_ly(data = df, x = ~ TWAS.p, y = ~ pip, color = ~ ifcausal,
  #                colors = c( "salmon", "darkgreen"), type ="scatter", text = ~ paste("Name: ", paste0(runtag,":",name),
  #                                                                  "\nChr: ", chr,  "\nPos:", pos))
                                                                    
}
tags <- paste(rep(1:2, each = 5), 1:5, sep = "-")
gwasfs <-  paste0(simdatadir, "20201001-", tags, ".exprgwas.txt.gz")
susieIfs2 <- paste0(outputdir, "20201001-", tags, ".config1.susieI.txt")

f1 <- caliFDP_plot(gwasfs[1:5], susieIfs2[1:5], format= "susie", main = "TWAS FDP (low power)")
FDP at bonferroni corrected p = 0.05:  0.5254237
f2 <- caliFDP_plot(gwasfs[6:10], susieIfs2[6:10], format= "susie", main ="TWAS FDP (high power)")
FDP at bonferroni corrected p = 0.05:  0.5250597
gridExtra::grid.arrange(f1, f2, ncol =2)

Version Author Date
447a401 simingz 2020-11-13

L’s effect on SuSiE

We run SuSiE on each block with true priors and true prior variances. We tested the difference between using L=1 and L=5. We only did this for the low power setting.

ncs_ncausal <- function(susief, main = "nCSvs.nCausal") {
  dt <- fread(susief, header = T)
  ncausal <- dt[, sum(ifcausal), by=list(b,rn)]
  nCS <- dt[, max(cs_index), by=list(b,rn)]
  plot(jitter(ncausal$V1), jitter(nCS$V1), xlab= "No.causal/block", ylab = "No. credible set/block", main= main)
  return(dt)
}

sumPIP_ncausal <- function(susief, main = "sumPIPvs.nCausal"){
  dt <- fread(susief, header = T)
  ncausal <- dt[, sum(ifcausal), by=list(b,rn)]
  sumpip <- dt[, sum(susie_pip), by=list(b,rn)]
  plot(jitter(ncausal$V1), sumpip$V1, xlab= "No.causal/block", ylab = "Regional sum PIP", main = main)
  return(dt)
}

L=1

Note: the data has been added some noise (both X and Y axes for credible set vs. no.causal, X axis for sum of PIP vs. no.causal) for visualization purposes.

tags <- paste(rep(1, each = 5), 1:5, sep = "-")
susieIfs2 <- paste0(outputdir, "20201001-", tags, ".config11.susieI.txt")
f1 <- caliPIP_plot(susieIfs2, main = "L=1")
f2 <- ncausal_plot(susieIfs2, format = "susie", main = "No. Causal Genes")
gridExtra::grid.arrange(f1, f2, ncol =2)

Version Author Date
41a7f93 simingz 2020-11-24 
par(mfrow=c(1,2))
dt <- sumPIP_ncausal(susieIfs2[2], main = "L=1")
dt <-ncs_ncausal(susieIfs2[2], main = "L=1, CS=0.95")

Version Author Date
07891d9 simingz 2020-11-25
41a7f93 simingz 2020-11-24
  • True parameters: gene \(\pi_1\), 0.0498753, SNP \(\pi_1\), 0.0025019. [enrichment = 19.9347551].
  • Average PIP: gene 0.0319922, SNP 0.0019147 . [enrichment = 16.7090917].
  • Proportion in credible set (0.95): gene 0.007606, SNP 0.0015938. [enrichment = 4.7722393].

L=5

Note: the data has been added some noise (both X and Y axes for credible set vs. no.causal, X axis for sum of PIP vs. no.causal) for visualization purposes.

tags <- paste(rep(1, each = 5), 1:5, sep = "-")
susieIfs2 <- paste0(outputdir, "20201001-", tags, ".config12.susieI.txt")

f1 <- caliPIP_plot(susieIfs2, main = "L=5")
f2 <- ncausal_plot(susieIfs2, format = "susie", main = "No. Causal Genes")
gridExtra::grid.arrange(f1, f2, ncol =2)

Version Author Date
41a7f93 simingz 2020-11-24 
par(mfrow=c(1,3))
dt <- sumPIP_ncausal(susieIfs2[2], main = "L=5")
dt <- ncs_ncausal(susieIfs2[2], main = "L=5, CS=0.95")

susieIfs2 <- paste0(outputdir, "20201001-", tags, ".config13.susieI.txt")
dt2 <- ncs_ncausal(susieIfs2[2], main = "L=5, CS=0.8")

Version Author Date
41a7f93 simingz 2020-11-24
  • True parameters: gene \(\pi_1\), 0.0498753, SNP \(\pi_1\), 0.0025019. [enrichment = 19.9347551].
  • Average PIP: gene 0.1409154, SNP 0.0091872 . [enrichment = 15.3381785].
  • Proportion in credible set (0.95): gene 0.0081047, SNP 0.0021012. [enrichment = 3.8572089].
  • Proportion in credible set (0.8): gene 0.0134663, SNP 0.0022246. [enrichment = 6.0533386].

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.SNP_susietruth"], 1:nrow(df), 2 + 1:nrow(df)/nrow(df)/3)
  s <- cbind(df[,"pi1.SNP_susieI"], 1:nrow(df), 3 + 1:nrow(df)/nrow(df)/3)
  t <- df[1,"pi1.SNP_truth"]
  dfp <- rbind(st,s)
  plot(dfp[,3], dfp[,1], col = dfp[,2], pch = 19, ylab = "SNP 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", ...)
}

susieI (1)

  • Regions: all regions, 500kb uniform regions.
  • Susie run parameters: L=1. We initialize with prior for genes and SNPs as uniform. gene.pve ~ 0.1, snp.pve ~ 0.5. Null weight is calculated based on prior of genes and SNPs ( 1 - sum of priors for snps and genes). Prior variance and residual variance were calculated by SUSIE for each region.
tags <- paste(rep(1:2, each = 5), 1:5, sep = "-")
phenofs <- paste0(simdatadir, "simu_20201001-", tags, "-pheno.Rd")
gwasfs <-  paste0(simdatadir, "20201001-", tags, ".exprgwas.txt.gz")
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,3))
plot_param(df[1:5,], main = "low power")
plot_param(df[6:10,], main = "high power")

Version Author Date
2e5d53d simingz 2020-11-04
f057bb5 simingz 2020-11-03 
f1 <- caliPIP_plot(susieIfs2[1:5], format ="susie", main = "PIP(low power)" )
f2 <- caliPIP_plot(susieIfs2[6:10], format ="susie", main = "PIP(high power)" )
f3 <- ncausal_plot(susieIfs2[1:5], format = "susie", main = "No. Causal Genes (low power)")
f4 <- ncausal_plot(susieIfs2[6:10], format = "susie", main = "No. Causal Genes (high power)")
gridExtra::grid.arrange(f1, f2, f3, f4, ncol =2)

Version Author Date
983dea2 simingz 2020-11-14
447a401 simingz 2020-11-13

susieI (2)

  • Regions: all regions, LD-defined regions.
  • Susie run parameters: L=1. We initialize with prior for genes and SNPs as uniform. gene.pve ~ 0.1, snp.pve ~ 0.5. Null weight is calculated based on prior of genes and SNPs ( 1 - sum of priors for snps and genes). Prior variance and residual variance were calculated by SUSIE for each region.
tags <- paste(rep(1:2, each = 5), 1:5, sep = "-")
phenofs <- paste0(simdatadir, "simu_20201001-", tags, "-pheno.Rd")
susieIfs <- paste0(outputdir, "20201001-", tags, ".config6.susieIres.Rd")
susieIfs2 <- paste0(outputdir, "20201001-", tags, ".config6.susieI.txt")

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

Version Author Date
447a401 simingz 2020-11-13 
f1 <- caliPIP_plot(susieIfs2[1:5], format ="susie", main = "PIP(low power)" )
f2 <- caliPIP_plot(susieIfs2[6:10], format ="susie", main = "PIP(high power)" )
f3 <- ncausal_plot(susieIfs2[1:5], format = "susie", main = "No. Causal Genes (low power)")
f4 <- ncausal_plot(susieIfs2[6:10], format = "susie", main = "No. Causal Genes (high power)")
gridExtra::grid.arrange(f1, f2, f3, f4, ncol =2)

Version Author Date
983dea2 simingz 2020-11-14

susieI (3)

  • Regions: all regions, LD-defined regions.
  • Susie run parameters: L=10. We initialize with prior for genes and SNPs as uniform. gene.pve ~ 0.1, snp.pve ~ 0.5. Null weight is calculated based on prior of genes and SNPs ( 1 - sum of priors for snps and genes). Prior variance and residual variance were calculated by SUSIE for each region.
tags <- paste(rep(1:2, each = 5), 1:5, sep = "-")
phenofs <- paste0(simdatadir, "simu_20201001-", tags, "-pheno.Rd")
susieIfs <- paste0(outputdir, "20201001-", tags, ".config8.susieIres.Rd")
susieIfs2 <- paste0(outputdir, "20201001-", tags, ".config8.susieI.txt")

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

Version Author Date
41a7f93 simingz 2020-11-24
6fa5cfa simingz 2020-11-19
447a401 simingz 2020-11-13 
f1 <- caliPIP_plot(susieIfs2[1:5], format ="susie", main = "PIP(low power)" )
f2 <- caliPIP_plot(susieIfs2[6:10], format ="susie", main = "PIP(high power)" )
f3 <- ncausal_plot(susieIfs2[1:5], format = "susie", main = "No. Causal Genes (low power)")
f4 <- ncausal_plot(susieIfs2[6:10], format = "susie", main = "No. Causal Genes (high power)")
gridExtra::grid.arrange(f1, f2, f3, f4, ncol =2)

Version Author Date
41a7f93 simingz 2020-11-24
6fa5cfa simingz 2020-11-19
983dea2 simingz 2020-11-14

susieI (4)

  • Regions: all regions, LD-defined regions.
  • Susie run parameters: L=1. We use true prior variance and plug into SUSIE, so this should be the same as SER/EM (2). We initialize with prior for genes and SNPs as uniform. gene.pve ~ 0.1, snp.pve ~ 0.5. Null weight is calculated based on prior of genes and SNPs ( 1 - sum of priors for snps and genes).
tags <- paste(rep(1:2, each = 5), 1:5, sep = "-")
phenofs <- paste0(simdatadir, "simu_20201001-", tags, "-pheno.Rd")
susieIfs <- paste0(outputdir, "20201001-", tags, ".config9.susieIres.Rd")
susieIfs2 <- paste0(outputdir, "20201001-", tags, ".config9.susieI.txt")

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

Version Author Date
6fa5cfa simingz 2020-11-19 
f1 <- caliPIP_plot(susieIfs2[1:5], format ="susie", main = "PIP(low power)" )
f2 <- caliPIP_plot(susieIfs2[6:10], format ="susie", main = "PIP(high power)" )
f3 <- ncausal_plot(susieIfs2[1:5], format = "susie", main = "No. Causal Genes (low power)")
f4 <- ncausal_plot(susieIfs2[6:10], format = "susie", main = "No. Causal Genes (high power)")
gridExtra::grid.arrange(f1, f2, f3, f4, ncol =2)

Version Author Date
6fa5cfa simingz 2020-11-19

susieI (5)

  • Regions: all regions, LD-defined regions.
  • Susie run parameters: L=5. We use true prior variance and plug into SUSIE. We initialize with prior for genes and SNPs as uniform. gene.pve ~ 0.1, snp.pve ~ 0.5. Null weight is calculated based on prior of genes and SNPs ( 1 - sum of priors for snps and genes).
tags <- paste(rep(1:2, each = 5), 1:5, sep = "-")
phenofs <- paste0(simdatadir, "simu_20201001-", tags, "-pheno.Rd")
susieIfs <- paste0(outputdir, "20201001-", tags, ".config10.susieIres.Rd")
susieIfs2 <- paste0(outputdir, "20201001-", tags, ".config10.susieI.txt")

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

Version Author Date
e81ebe1 simingz 2020-11-20
6fa5cfa simingz 2020-11-19 
f1 <- caliPIP_plot(susieIfs2[1:5], format ="susie", main = "PIP(low power)" )
f2 <- caliPIP_plot(susieIfs2[6:10], format ="susie", main = "PIP(high power)" )
f3 <- ncausal_plot(susieIfs2[1:5], format = "susie", main = "No. Causal Genes (low power)")
f4 <- ncausal_plot(susieIfs2[6:10], format = "susie", main = "No. Causal Genes (high poweunr)")
gridExtra::grid.arrange(f1, f2, f3, f4, ncol =2)

Version Author Date
e81ebe1 simingz 2020-11-20
6fa5cfa simingz 2020-11-19

Select block + susieI (1)

  • Regions: all regions 10 iterations. Then filter out regions with probability of having two or more effects > 0.2 to estimate paramters 20 iterations. then all regions to get all PIP. LD-defined regions.
  • Susie run parameters: L=1. We use true prior variance and plug into SUSIE. We initialize with prior for genes and SNPs as uniform. gene.pve ~ 0.1, snp.pve ~ 0.5. We use null weight when run susie. We run 10 iterations. Then filter out regions with probability of having two or more effects. Then rerun on selected regions to estimate parameters. Lastly, run on all regions with L=1 to get PIP.
tags <- paste(rep(1:2, each = 5), 1:5, sep = "-")
phenofs <- paste0(simdatadir, "simu_20201001-", tags, "-pheno.Rd")
susieIfs <- paste0(outputdir, "20201001-", tags, ".config14.fl.susieIres.Rd")
susieIfs2 <- paste0(outputdir, "20201001-", tags, ".config14.fl.susieI.txt")

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

Version Author Date
6e78356 simingz 2020-12-03
07891d9 simingz 2020-11-25 
a <- fread(susieIfs2[1], header =T)
cat("No.blocks selected for parameter estimation (low power):", nrow(unique(a[, c("b", "rn")])), "\n")
No.blocks selected for parameter estimation (low power): 1193 
a <- fread(susieIfs2[6], header =T)
cat("No.blocks selected for parameter estimation (high power):", nrow(unique(a[, c("b", "rn")])), "\n")
No.blocks selected for parameter estimation (high power): 1207

Select block + susieI (2)

  • Regions: all regions 10 iterations. Then filter out regions with probability of having two or more effects > 0.1 to estimate paramters 20 iterations. then all regions to get all PIP. LD-defined regions.
  • Susie run parameters: L=1. We use true prior variance and plug into SUSIE. We initialize with prior for genes and SNPs as uniform. gene.pve ~ 0.1, snp.pve ~ 0.5. We use null weight when run susie. We run 10 iterations. Then filter out regions with probability of having two or more effects. Then rerun on selected regions to estimate parameters. Lastly, run on all regions with L=1 to get PIP.
tags <- paste(rep(1:2, each = 5), 1:5, sep = "-")
phenofs <- paste0(simdatadir, "simu_20201001-", tags, "-pheno.Rd")
susieIfs <- paste0(outputdir, "20201001-", tags, ".config15.fl.susieIres.Rd")
susieIfs2 <- paste0(outputdir, "20201001-", tags, ".config15.fl.susieI.txt")

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

Version Author Date
6e78356 simingz 2020-12-03
07891d9 simingz 2020-11-25 
a <- fread(susieIfs2[1], header =T)
cat("No.blocks selected for parameter estimation (low power):", nrow(unique(a[, c("b", "rn")])), "\n")
No.blocks selected for parameter estimation (low power): 438 
a <- fread(susieIfs2[6], header =T)
cat("No.blocks selected for parameter estimation (high power):", nrow(unique(a[, c("b", "rn")])), "\n")
No.blocks selected for parameter estimation (high power): 447

Select block + susieI (3)

  • Regions: all regions 5 iterations. Then filter out regions with probability of having two or more effects > 0.2 to estimate paramters 20 iterations. then all regions to get all PIP. LD-defined regions.
  • Susie run parameters: L=1. We estimate prior variance based on EM and plug into SUSIE. We initialize with prior for genes and SNPs as uniform. gene.pve ~ 0.1, snp.pve ~ 0.5. We use null weight when run susie. We run 5 iterations. Then filter out regions with probability of having two or more effects. Then rerun on selected regions to estimate parameters (20 iterations). Lastly, run on all regions with L=5 to get PIP.
tags <- paste(rep(1:2, each = 5), 1:5, sep = "-")
phenofs <- paste0(simdatadir, "simu_20201001-", tags, "-pheno.Rd")
susieIfs <- paste0(outputdir, "20201001-", tags, ".config16.fl.susieIres.Rd")
susieIfs2 <- paste0(outputdir, "20201001-", tags, ".config16.fl.susieI.txt")

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

Version Author Date
6e78356 simingz 2020-12-03 
a <- fread(susieIfs2[1], header =T)
cat("No.blocks selected for parameter estimation (low power):", nrow(unique(a[, c("b", "rn")])), "\n")
No.blocks selected for parameter estimation (low power): 959 
a <- fread(susieIfs2[6], header =T)
cat("No.blocks selected for parameter estimation (high power):", nrow(unique(a[, c("b", "rn")])), "\n")
No.blocks selected for parameter estimation (high power): 960 
susieIfs2 <- paste0(outputdir, "20201001-", tags, ".config16.flrerun.susieI.txt")
f1 <- caliPIP_plot(susieIfs2[1:5], format ="susie", main = "PIP(low power)")
f2 <- caliPIP_plot(susieIfs2[6:10], format ="susie", main = "PIP(high power)")
f3 <- ncausal_plot(susieIfs2[1:5], format = "susie", main = "No. Causal Genes (low power)")
f4 <- ncausal_plot(susieIfs2[6:10], format = "susie", main = "No. Causal Genes (high power)")
gridExtra::grid.arrange(f1, f2, f3, f4, ncol =2)

Version Author Date
6e78356 simingz 2020-12-03 
par(mfrow=c(1,2))
f1 <- scatter_plot_PIP_p(susieIfs2[1:5], gwasfs[1:5], pipformat = "susie", main = "low power")
f2 <- scatter_plot_PIP_p(susieIfs2[6:10], gwasfs[6:10], pipformat = "susie", main = "high power")

Select block + susieI (4)

  • Regions: all regions 5 iterations. Then filter out regions with probability of having two or more effects > 0.1 to estimate paramters 20 iterations. then all regions to get all PIP. LD-defined regions.
  • Susie run parameters: L=1. We estimate prior variance based on EM and plug into SUSIE. We initialize with prior for genes and SNPs as uniform. gene.pve ~ 0.1, snp.pve ~ 0.5. We use null weight when run susie. We run 5 iterations. Then filter out regions with probability of having two or more effects. Then rerun on selected regions to estimate parameters (20 iterations). Lastly, run on all regions with L=5 to get PIP.
tags <- paste(rep(1:2, each = 5), 1:5, sep = "-")
phenofs <- paste0(simdatadir, "simu_20201001-", tags, "-pheno.Rd")
susieIfs <- paste0(outputdir, "20201001-", tags, ".config17.fl.susieIres.Rd")
susieIfs2 <- paste0(outputdir, "20201001-", tags, ".config17.fl.susieI.txt")

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

Version Author Date
6e78356 simingz 2020-12-03 
a <- fread(susieIfs2[1], header =T)
cat("No.blocks selected for parameter estimation (low power):", nrow(unique(a[, c("b", "rn")])), "\n")
No.blocks selected for parameter estimation (low power): 351 
a <- fread(susieIfs2[6], header =T)
cat("No.blocks selected for parameter estimation (high power):", nrow(unique(a[, c("b", "rn")])), "\n")
No.blocks selected for parameter estimation (high power): 351 
susieIfs2 <- paste0(outputdir, "20201001-", tags, ".config17.flrerun.susieI.txt")
f1 <- caliPIP_plot(susieIfs2[1:5], format ="susie", main = "PIP(low power)")
f2 <- caliPIP_plot(susieIfs2[6:10], format ="susie", main = "PIP(high power)")
f3 <- ncausal_plot(susieIfs2[1:5], format = "susie", main = "No. Causal Genes (low power)")
f4 <- ncausal_plot(susieIfs2[6:10], format = "susie", main = "No. Causal Genes (high power)")
gridExtra::grid.arrange(f1, f2, f3, f4, ncol =2)

Version Author Date
6e78356 simingz 2020-12-03 
par(mfrow=c(1,2))
f1 <- scatter_plot_PIP_p(susieIfs2[1:5], gwasfs[1:5], pipformat = "susie", main = "low power")
f2 <- scatter_plot_PIP_p(susieIfs2[6:10], gwasfs[6:10], pipformat = "susie", main = "high power")

Single effect model/EM (1)

  • Regions: all regions, 500kb uniform regions.
  • Run parameters: We initialize with prior for genes and SNPs as uniform. gene.pve ~ 0.1, snp.pve ~ 0.5. Null weight is calculated based on prior of genes and SNPs ( 1 - sum of priors for snps and genes). We use the true prior variance for genes and SNPs.
tags <- paste(rep(1:2, each = 5), 1:5, sep = "-")
phenofs <- paste0(simdatadir, "simu_20201001-", tags, "-pheno.Rd")
susieIfs <- paste0(outputdir, "20201001-", tags, ".config1.SERIres.Rd")
susieIfs2 <- paste0(outputdir, "20201001-", tags, ".config1.SERI.txt")

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

Version Author Date
2e5d53d simingz 2020-11-04
f057bb5 simingz 2020-11-03
f9c55de simingz 2020-10-30
6526634 simingz 2020-10-23 
f1 <- caliPIP_plot(susieIfs2[1:5], format ="SER", main = "PIP(low power)")
f2 <- caliPIP_plot(susieIfs2[6:10], format ="SER", main = "PIP(high power)")
f3 <- ncausal_plot(susieIfs2[1:5], format = "SER", main = "No. Causal Genes (low power)")
f4 <- ncausal_plot(susieIfs2[6:10], format = "SER", main = "No. Causal Genes (high power)")
gridExtra::grid.arrange(f1, f2, f3, f4, ncol =2)

Version Author Date
983dea2 simingz 2020-11-14
447a401 simingz 2020-11-13

Single effect model/EM (2)

  • Regions: all regions, defined based on LDetect.
  • Run parameters: We initialize with prior for genes and SNPs as uniform. gene.pve ~ 0.1, snp.pve ~ 0.5. Null weight is calculated based on prior of genes and SNPs ( 1 - sum of priors for snps and genes). We use the true prior variance for genes and SNPs.
tags <- paste(rep(1:2, each = 5), 1:5, sep = "-")
phenofs <- paste0(simdatadir, "simu_20201001-", tags, "-pheno.Rd")
susieIfs <- paste0(outputdir, "20201001-", tags, ".config6.SERIres.Rd")
susieIfs2 <- paste0(outputdir, "20201001-", tags, ".config6.SERI.txt")

df <- show_param(phenofs, susieIfs, susieIfs2)
Avoidable 0.601 seconds. This file is very unusual: it ends abruptly without a final newline, and also its size is a multiple of 4096 bytes. Please properly end the last row with a newline using for example 'echo >> file' to avoid this  time to copy.
Avoidable 4.501 seconds. This file is very unusual: it ends abruptly without a final newline, and also its size is a multiple of 4096 bytes. Please properly end the last row with a newline using for example 'echo >> file' to avoid this  time to copy.
par(mfrow = c(2,3))
plot_param(df[1:5,], main = "low power")
plot_param(df[6:10,], main = "high power")

Version Author Date
2e5d53d simingz 2020-11-04
f057bb5 simingz 2020-11-03 
f1 <- caliPIP_plot(susieIfs2[1:5], format ="SER", main = "PIP(low power)")
Avoidable 2.618 seconds. This file is very unusual: it ends abruptly without a final newline, and also its size is a multiple of 4096 bytes. Please properly end the last row with a newline using for example 'echo >> file' to avoid this  time to copy.
Avoidable 4.658 seconds. This file is very unusual: it ends abruptly without a final newline, and also its size is a multiple of 4096 bytes. Please properly end the last row with a newline using for example 'echo >> file' to avoid this  time to copy.
f2 <- caliPIP_plot(susieIfs2[6:10], format ="SER", main = "PIP(high power)")
f3 <- ncausal_plot(susieIfs2[1:5], format = "SER", main = "No. Causal Genes (low power)")
Avoidable 0.537 seconds. This file is very unusual: it ends abruptly without a final newline, and also its size is a multiple of 4096 bytes. Please properly end the last row with a newline using for example 'echo >> file' to avoid this  time to copy.
Avoidable 4.278 seconds. This file is very unusual: it ends abruptly without a final newline, and also its size is a multiple of 4096 bytes. Please properly end the last row with a newline using for example 'echo >> file' to avoid this  time to copy.
f4 <- ncausal_plot(susieIfs2[6:10], format = "SER", main = "No. Causal Genes (high power)")
gridExtra::grid.arrange(f1, f2, f3, f4, ncol =2)

Version Author Date
983dea2 simingz 2020-11-14
447a401 simingz 2020-11-13 
par(mfrow=c(1,2))
f1 <- scatter_plot_PIP_p(susieIfs2[1:5], gwasfs[1:5], pipformat = "SER", main = "low power")
Avoidable 0.571 seconds. This file is very unusual: it ends abruptly without a final newline, and also its size is a multiple of 4096 bytes. Please properly end the last row with a newline using for example 'echo >> file' to avoid this  time to copy.
Avoidable 4.547 seconds. This file is very unusual: it ends abruptly without a final newline, and also its size is a multiple of 4096 bytes. Please properly end the last row with a newline using for example 'echo >> file' to avoid this  time to copy.
f2 <- scatter_plot_PIP_p(susieIfs2[6:10], gwasfs[6:10], pipformat = "SER", main = "high power")

Version Author Date
983dea2 simingz 2020-11-14
447a401 simingz 2020-11-13

Single effect model/EM (3)

  • Regions: regions with at most 1 causal signal, regions are defined based on LDetect.
  • Run parameters: We initialize with prior for genes and SNPs as uniform. gene.pve ~ 0.1, snp.pve ~ 0.5. Null weight is calculated based on prior of genes and SNPs ( 1 - sum of priors for snps and genes). We use the true prior variance for genes and SNPs.
tags <- paste(rep(1:2, each = 5), 1:5, sep = "-")
phenofs <- paste0(simdatadir, "simu_20201001-", tags, "-pheno.Rd")
susieIfs <- paste0(outputdir, "20201001-", tags, ".config7.SERIres.Rd")
susieIfs2 <- paste0(outputdir, "20201001-", tags, ".config7.SERI.txt")

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

Version Author Date
2e5d53d simingz 2020-11-04
f057bb5 simingz 2020-11-03 
f1 <- caliPIP_plot(susieIfs2[1:5], format ="SER", main = "PIP(low power)")
f2 <- caliPIP_plot(susieIfs2[6:10], format ="SER", main = "PIP(high power)")
f3 <- ncausal_plot(susieIfs2[1:5], format = "SER", main = "No. Causal Genes (low power)")
f4 <- ncausal_plot(susieIfs2[6:10], format = "SER", main = "No. Causal Genes (high power)")
gridExtra::grid.arrange(f1, f2, f3, f4, ncol =2)

Version Author Date
983dea2 simingz 2020-11-14
447a401 simingz 2020-11-13

Summary stats version

  • The run setting is the same as in “Select block + susieI (3)” except that we used summary statistics and in-sample LD as input.
simdatadir <- "~/causalTWAS/simulations/simulation_ashtest_20201001/"
outputdir <- "~/causalTWAS/simulations/simulation_susieI_rss_20201001/"
susiedir <- "~/causalTWAS/simulations/simulation_susieI_rss_20201001/"
tags <- paste(rep(1:2, each = 5), 1:5, sep = "-")
phenofs <- paste0(simdatadir, "simu_20201001-", tags, "-pheno.Rd")
susieIfs <- paste0(outputdir, "20201001-", tags, ".config1.fl.susieIrssres.Rd")
susieIfs2 <- paste0(outputdir, "20201001-", tags, ".config1.fl.susieIrss.txt")

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

a <- fread(susieIfs2[1], header =T)
cat("No.blocks selected for parameter estimation (low power):", nrow(unique(a[, c("b", "rn")])), "\n")
No.blocks selected for parameter estimation (low power): 946 
a <- fread(susieIfs2[6], header =T)
cat("No.blocks selected for parameter estimation (high power):", nrow(unique(a[, c("b", "rn")])), "\n")
No.blocks selected for parameter estimation (high power): 946 
susieIfs2 <- paste0(outputdir, "20201001-", tags, ".config1.flrerun.susieIrss.txt")
f1 <- caliPIP_plot(susieIfs2[1:5], format ="susie", main = "PIP(low power)")
f2 <- caliPIP_plot(susieIfs2[6:10], format ="susie", main = "PIP(high power)")
f3 <- ncausal_plot(susieIfs2[1:5], format = "susie", main = "No. Causal Genes (low power)")
f4 <- ncausal_plot(susieIfs2[6:10], format = "susie", main = "No. Causal Genes (high power)")
gridExtra::grid.arrange(f1, f2, f3, f4, ncol =2)


sessionInfo()
R version 3.6.1 (2019-07-05)
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] parallel  stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] plotrix_3.7-6       cowplot_1.0.0       snpStats_1.34.0    
 [4] Matrix_1.2-18       survival_2.44-1.1   doParallel_1.0.14  
 [7] iterators_1.0.10    foreach_1.4.4       stringr_1.4.0      
[10] plyr_1.8.4          tidyr_1.1.0         plotly_4.9.0       
[13] ggplot2_3.2.1       data.table_1.13.2   mr.ash.alpha_0.1-34

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.5          lattice_0.20-38     assertthat_0.2.1   
 [4] rprojroot_1.3-2     digest_0.6.20       R6_2.4.0           
 [7] backports_1.1.4     evaluate_0.14       httr_1.4.1         
[10] highr_0.8           pillar_1.4.2        zlibbioc_1.30.0    
[13] rlang_0.4.6         lazyeval_0.2.2      whisker_0.3-2      
[16] R.utils_2.9.0       R.oo_1.22.0         rmarkdown_1.13     
[19] labeling_0.3        splines_3.6.1       htmlwidgets_1.3    
[22] munsell_0.5.0       compiler_3.6.1      httpuv_1.5.1       
[25] xfun_0.8            pkgconfig_2.0.2     BiocGenerics_0.30.0
[28] htmltools_0.3.6     tidyselect_1.1.0    gridExtra_2.3      
[31] tibble_2.1.3        bigstatsr_0.9.9     workflowr_1.6.2    
[34] codetools_0.2-16    viridisLite_0.3.0   crayon_1.3.4       
[37] dplyr_0.8.3         withr_2.1.2         later_0.8.0        
[40] R.methodsS3_1.7.1   grid_3.6.1          jsonlite_1.6       
[43] gtable_0.3.0        lifecycle_0.1.0     git2r_0.26.1       
[46] magrittr_1.5        scales_1.1.0        stringi_1.4.3      
[49] farver_2.0.1        fs_1.3.1            promises_1.0.1     
[52] vctrs_0.3.1         tools_3.6.1         glue_1.3.1         
[55] purrr_0.3.4         yaml_2.2.0          colorspace_1.4-1   
[58] knitr_1.23