Fine-mapping Height ZBTB38 – SuSiE check
Yuxin Zou
10/07/2019
We perform some check for the SuSiE result on region around ZBTB38.
Load pacakges:
library(readr)
library(dplyr)
library(gridExtra)
library(susieR)Load plotting functions:
knitr::read_chunk("../scripts/plots.R")library(ggplot2)
#' plot gene name annotations
#' @param dat a matrix of gene names with 'start' and 'end' base-pair position
#' @param xrange range of x axis, base-pair position
plot_geneName = function(dat, xrange, chr){
ngene = 2:nrow(dat)
line = 1
dat$lines = NA
dat$lines[1] = 1
gene.end = dat[1, 'end']
while(length(ngene) != 0){
id = which(dat[ngene, 'start'] > gene.end + 0.02)[1]
if(!is.na(id)){
dat$lines[ngene[id]] = line
gene.end = dat[ngene[id],'end']
ngene = ngene[-id]
}else{
line = line + 1
dat$lines[ngene[1]] = line
gene.end = dat[ngene[1],'end']
ngene = ngene[-1]
}
}
dat$start = pmax(dat$start, xrange[1])
dat$end = pmin(dat$end, xrange[2])
dat$mean = rowMeans(dat[,c('start', 'end')])
pl = ggplot(dat, aes(xmin = xrange[1], xmax = xrange[2])) + xlim(xrange[1], xrange[2]) + ylim(min(-dat$lines-0.6), -0.8) +
geom_rect(aes(xmin = start, xmax = end, ymin = -lines-0.05, ymax = -lines+0.05), fill='blue') +
geom_text(aes(x = mean, y=-lines-0.4, label=geneName), size=4) +
xlab(paste0('base-pair position (Mb) on chromosome ', chr)) + ylab('Gene') +
theme_bw() + theme(axis.text.x=element_blank(),
axis.ticks = element_blank(),
axis.text.y=element_blank(),
axis.title = element_text(size=15),
plot.title=element_text(size=11),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
pl
}
discrete_gradient_pal <- function(colours, bins = 5) {
ramp <- scales::colour_ramp(colours)
function(x) {
if (length(x) == 0) return(character())
i <- floor(x * bins)
i <- ifelse(i > bins-1, bins-1, i)
ramp(i/(bins-1))
}
}
scale_colour_discrete_gradient <- function(..., colours, bins = 5, na.value = "grey50", guide = "colourbar", aesthetics = "colour", colors) {
colours <- if (missing(colours))
colors
else colours
continuous_scale(
aesthetics,
"discrete_gradient",
discrete_gradient_pal(colours, bins),
na.value = na.value,
guide = guide,
...
)
}
#' Locuszoom plot
#' @param z a vector of z scores with SNP names
#' @param pos base-pair positions
#' @param gene.pos.map a matrix of gene names with 'start' and 'end' base-pair position
#' @param z.ref.name the reference SNP
#' @param ld correlations between teh reference SNP and the rests
#' @param title title of the plot
#' @param title.size the size of the title
#' @param true the true value
#' @param y.height height of -log10(p) plot and height of the gene name annotation plot
#' @param y.lim range of y axis
locus.zoom = function(z, pos, chr, gene.pos.map=NULL, z.ref.name=NULL, ld=NULL,
title = NULL, title.size = 10, true = NULL,
y.height=c(5,1.5), y.lim=NULL, xrange=NULL){
if(is.null(xrange)){
xrange = c(min(pos), max(pos))
}
tmp = data.frame(POS = pos, p = -(pnorm(-abs(z), log.p = T) + log(2))/log(10))
pl_zoom = ggplot(tmp, aes(x = POS, y = p)) + geom_point(color = 'darkblue') +
ylab("-log10(p value)") + ggtitle(title) + xlim(xrange[1], xrange[2]) +
theme_bw() + theme(axis.title.x=element_blank(),
plot.title = element_text(size=title.size))
if(!is.null(ld) && !is.null(z.ref.name)){
tmp$ref = names(z) == z.ref.name
tmp$r2 = ld^2
pl_zoom = ggplot(tmp, aes(x = POS, y = p, shape = ref, size=ref, color=r2)) + geom_point() +
ylab("-log10(p value)") + ggtitle(title) + xlim(xrange[1], xrange[2]) +
scale_color_gradientn(colors = c("darkblue", "deepskyblue", "lightgreen", "orange", "red"),
values = seq(0,1,0.2), breaks=seq(0,1,0.2)) +
# scale_colour_discrete_gradient(
# colours = c("darkblue", "deepskyblue", "lightgreen", "orange", "red"),
# limits = c(0, 1.01),
# breaks = c(0,0.2,0.4,0.6,0.8,1),
# guide = guide_colourbar(nbin = 100, raster = FALSE, frame.colour = "black", ticks.colour = NA)
# ) +
scale_shape_manual(values=c(20, 18), guide=FALSE) + scale_size_manual(values=c(2,5), guide=FALSE) +
theme_bw() + theme(axis.title.x=element_blank(),
axis.title.y = element_text(size=15),axis.text = element_text(size=15),
plot.title = element_text(size=title.size))
}
if(!is.null(y.lim)){
pl_zoom = pl_zoom + ylim(y.lim[1], y.lim[2])
}
# pl_zoom = pl_zoom + geom_hline(yintercept=-log10(5e-08), linetype='dashed', color = 'red')
if(!is.null(true)){
tmp.true = data.frame(POS = which(true!=0), p = tmp$p[which(true!=0)],
ref = (names(z) == z.ref.name)[which(true!=0)],
label = paste0('SNP',1:length(which(true!=0))))
pl_zoom = pl_zoom + geom_point(data=tmp.true, aes(x=POS, y=p),
color='red', show.legend = FALSE, shape=1, stroke = 1) +
geom_text(data=tmp.true, aes(x = POS-30, y=p+1, label=label), size=3, color='red')
}
if(!is.null(gene.pos.map)){
pl_gene = plot_geneName(gene.pos.map, xrange = xrange, chr=chr)
g = egg::ggarrange(pl_zoom, pl_gene, nrow=2, heights = y.height, draw=FALSE)
}else{
g = pl_zoom
}
g
}
#' SuSiE plot with Locuszoom plot
#' @param z a vector of z scores with SNP names
#' @param model the fitted SuSiE model
#' @param pos base-pair positions
#' @param gene.pos.map a matrix of gene names with 'start' and 'end' base-pair position
#' @param z.ref.name the reference SNP
#' @param ld correlations between teh reference SNP and the rests
#' @param title title of the plot
#' @param title.size the size of the title
#' @param true the true value
#' @param plotz whether to plot locuszoom plot
#' @param y.lim range of y axis
#' @param y.susie the y axis of the SuSiE plot, 'PIP' or 'p', 'p' refers to -log10(p)
susie_plot_locuszoom = function(z, model, pos, chr, gene.pos.map = NULL, z.ref.name, ld,
title = NULL, title.size = 10, true = NULL,
plotz = TRUE, y.lim=NULL, y.susie='PIP', xrange=NULL){
if(is.null(xrange)){
xrange = c(min(pos), max(pos))
}
if(plotz){
pl_zoom = locus.zoom(z, pos = pos, chr = chr, ld=ld, z.ref.name = z.ref.name, title = title, title.size = title.size, y.lim=y.lim, xrange=xrange)
}
pip = model$pip
tmp = data.frame(POS = pos, PIP = pip, p = -(pnorm(-abs(z), log.p = T) + log(2))/log(10))
if(y.susie == 'PIP'){
if(plotz){
pl_susie = ggplot(tmp, aes(x = POS, y = PIP)) + geom_point(show.legend = FALSE, size=3) +
xlim(xrange[1], xrange[2]) +
theme_bw() + theme(axis.title.x=element_blank(), axis.text.x=element_blank())
}else{
pl_susie = ggplot(tmp, aes(x = POS, y = PIP)) + geom_point(show.legend = FALSE, size=3) +
xlim(xrange[1], xrange[2]) + ggtitle(title) +
theme_bw() + theme(axis.title.x=element_blank(),
axis.text = element_text(size=15),
axis.title.y = element_text(size=15),
plot.title = element_text(size=title.size))
}
}else if(y.susie == 'p'){
if(plotz){
pl_susie = ggplot(tmp, aes(x = POS, y = p)) + geom_point(show.legend = FALSE, size=3) +
ylab("-log10(p value)") + xlim(xrange[1], xrange[2]) +
theme_bw() + theme(axis.title.x=element_blank(), axis.text.x=element_blank())
# pl_susie = pl_susie + geom_hline(yintercept=-log10(5e-08), linetype='dashed', color = 'red')
}else{
pl_susie = ggplot(tmp, aes(x = POS, y = p)) + geom_point(show.legend = FALSE, size=3) +
ylab("-log10(p value)") + xlim(xrange[1], xrange[2]) + ggtitle(title) +
theme_bw() + theme(axis.title.x=element_blank(),
axis.text = element_text(size=15),
axis.title.y = element_text(size=15),
plot.title = element_text(size=title.size))
}
}
if(!is.null(true)){
tmp.true = data.frame(POS = pos[which(true!=0)], PIP = pip[which(true!=0)])
pl_susie = pl_susie + geom_point(data=tmp.true, aes(x=POS, y=PIP),
color='red', size=3, show.legend = FALSE)
}
model.cs = model$sets$cs
if(!is.null(model.cs)){
tmp$CS = numeric(length(z))
for(i in 1:length(model.cs)){
tmp$CS[model.cs[[i]]] = gsub('L', 'CS', names(model.cs)[i])
}
tmp.cs = tmp[unlist(model.cs),]
tmp.cs$CS = factor(tmp.cs$CS)
levels(tmp.cs$CS) = paste0('CS', 1:length(model.cs))
colors = c('red', 'cyan', 'green', 'orange', 'dodgerblue', 'violet', 'gold',
'#FF00FF', 'forestgreen', '#7A68A6')
if(y.susie == 'PIP'){
pl_susie = pl_susie + geom_point(data=tmp.cs, aes(x=POS, y=PIP, color=CS),
size=3, shape=1, stroke = 2) +
scale_color_manual(values=colors)
}else if(y.susie == 'p'){
pl_susie = pl_susie + geom_point(data=tmp.cs, aes(x=POS, y=p, color=CS),
shape=1, size=3, stroke=1.5) +
scale_color_manual(values=colors)
}
}
if(!is.null(gene.pos.map)){
pl_gene = plot_geneName(gene.pos.map, xrange = xrange, chr=chr)
if(plotz){
g = egg::ggarrange(pl_zoom, pl_susie, pl_gene, nrow=3, heights = c(4,4,1.5), draw=FALSE)
}else{
g = egg::ggarrange(pl_susie, pl_gene, nrow=2, heights = c(5.5,1.5), draw=FALSE)
}
}else{
if(plotz){
g = egg::ggarrange(pl_zoom, pl_susie, nrow=2, heights = c(4,4), draw=FALSE)
}else{
g = pl_susie
}
}
g
}locus.zoom.cs = function(z, cs, pos, chr, gene.pos.map=NULL, z.ref.name=NULL, ld=NULL, title = NULL, title.size = 10, xrange = NULL, y.lab=NULL){
if(is.null(y.lab)){
y.lab = '-log10(p value)'
}
tmp = data.frame(POS = pos, log10p = -(pnorm(-abs(z), log.p = T) + log(2))/log(10))
tmp$ref = names(z) == z.ref.name
tmp$r2 = ld^2
tmp$CS = rep(4, length(z))
tmp$CS[cs] = 16
tmp$CS = as.factor(tmp$CS)
pl_zoom = ggplot(tmp, aes(x = POS, y = log10p, shape = CS, color = r2, size=CS)) +
geom_point() + ylab(y.lab) +
scale_color_gradientn(colors = c("darkblue", "deepskyblue", "lightgreen", "orange", "red"),
values = seq(0,1,0.2), breaks=seq(0,1,0.2)) +
# scale_colour_discrete_gradient(
# colours = c("darkblue", "deepskyblue", "lightgreen", "orange", "red"),
# limits = c(0, 1.01),
# breaks = c(0,0.2,0.4,0.6,0.8,1),
# guide = guide_colourbar(nbin = 100, raster = FALSE, frame.colour = "black", ticks.colour = NA)) +
scale_shape_manual(values = c(4, 19), guide=FALSE) +
scale_size_manual(values=c(1.5,4), guide=FALSE) +
ggtitle(title) +
theme_bw() + theme(axis.title.x=element_blank(), axis.text=element_text(size=15),
axis.title.y=element_text(size=12),
plot.title = element_text(size=title.size))
tmp.sub = tmp[cs,]
pl_zoom = pl_zoom + geom_point(data = tmp.sub, aes(x=POS, y=log10p), shape=1, size=4, color='black', stroke=0.1)
if(!is.null(xrange)){
pl_zoom = pl_zoom + xlim(xrange[1], xrange[2])
}
pl_gene = plot_geneName(gene.pos.map, xrange = xrange, chr=chr)
g = egg::ggarrange(pl_zoom, pl_gene, nrow=2, heights = c(5.5,1.5), draw=FALSE)
g
}Get summary statistics:
ss.dat = readRDS('../data/height.ZBTB38.XtX.Xty.rds')
betas = as.vector(ss.dat$Xty/diag(ss.dat$XtX))
rss = c(ss.dat$yty) - betas * as.vector(ss.dat$Xty)
se = sqrt(rss/((ss.dat$n-1)*diag(ss.dat$XtX)))
z = betas/se
pval = 2*pnorm(-abs(z))
R = as.matrix(t(ss.dat$XtX * (1/sqrt(diag(ss.dat$XtX)))) * (1/ sqrt(diag(ss.dat$XtX))))
names(z) = rownames(R)Get gene data:
genes <- read_delim("../data/seq_gene.md.gz",delim = "\t",quote = "")
class(genes) <- "data.frame"
genes <- subset(genes,
group_label == "GRCh37.p5-Primary Assembly" &
feature_type == "GENE")
start.pos <- min(ss.dat$pos$POS)
stop.pos <- max(ss.dat$pos$POS)
plot.genes <- subset(genes,
chromosome == 3 &
((chr_start > start.pos & chr_start < stop.pos) |
(chr_stop > start.pos & chr_start < stop.pos)) & feature_type == 'GENE')
gene.pos.map = plot.genes %>% select(feature_name, chr_start, chr_stop)
colnames(gene.pos.map) = c('geneName', 'start', 'end')
gene.pos.map = as.data.frame(gene.pos.map)
gene.pos.map = gene.pos.map %>% mutate(start = start/1e6, end = end/1e6)
gene.pos.map = gene.pos.map[-c(4, 9, 10, 15),]SuSiE bhat result:
mod_bhat = susie_bhat(bhat = betas, shat = se, R = R, n = ss.dat$n, var_y = as.numeric(ss.dat$yty/(ss.dat$n-1)), track_fit=TRUE)
z.max = which.max(abs(z))
p1 = susie_plot_locuszoom(z, mod_bhat, pos = ss.dat$pos$POS/1e6, chr=3, gene.pos.map = gene.pos.map, ld = R[z.max,], z.ref.name = 'rs2871960_C')
p2 = susie_plot_locuszoom(z, mod_bhat, pos = ss.dat$pos$POS/1e6, chr=3, gene.pos.map = gene.pos.map, ld = R[z.max,], z.ref.name = 'rs2871960_C', y.susie ='p')
grid.arrange(p1, p2, ncol=2)
For 1Mb region about ZBTB38, SuSiE found 2 CSs. The SNP with the strongest marginal p value in CS1 is rs2871960 (p = 1.8670e-207). For CS2, the SNP with the strongest marginal p value is rs11919556 (p = 0.0211).
The correlation between rs2871960 and rs11919556 is 0.221051. The average correlation between SNPs in CS1 and CS2 is
round(mean(abs(R[mod_bhat$sets$cs$L1, mod_bhat$sets$cs$L2])), 4)## [1] 0.1849Zoom in plot:
p1 = susie_plot_locuszoom(z, mod_bhat, pos = ss.dat$pos$POS/1e6, chr=3, gene.pos.map = gene.pos.map, z.ref.name = names(z), ld = R[z.max,], title='ZBTB38 Credible Sets', plotz = FALSE, y.susie = 'p', xrange=c(140.9,141.5), title.size = 20)## Warning: Removed 1231 rows containing missing values (geom_point).## Warning: Removed 5 rows containing missing values (geom_rect).## Warning: Removed 5 rows containing missing values (geom_text).p1
# pdf('GitHub/finemap-uk-biobank/ZBTB38_CSs.pdf', height = 5, width=8)
# p1
# dev.off()CS 1
- Removing effect of top SNP in CS 2
library(data.table)
library(readr)
library(Matrix)
geno.file = 'height.ZBTB38.raw.gz'
cat("Reading genotype data.\n")
geno <- fread(geno.file,sep = "\t",header = TRUE,stringsAsFactors = FALSE)
class(geno) <- "data.frame"
# Extract the genotypes.
X <- as(as.matrix(geno[-(1:6)]),'dgCMatrix')
pheno.file <- "/gpfs/data/stephens-lab/finemap-uk-biobank/data/raw/height.csv.gz"
out.pheno.file <- "pheno.height.txt"
out.covar.file <- "covar.removeCS2.height.txt"
# LOAD PHENOTYPE and COVARIATES DATA
# -------------------
# Read the phenotype data from the CSV file.
cat("Reading phenotype data.\n")
pheno <- suppressMessages(read_csv(pheno.file))
class(pheno) <- "data.frame"
pheno$sex = factor(pheno$sex)
pheno$assessment_centre = factor(pheno$assessment_centre)
pheno$genotype_measurement_batch = factor(pheno$genotype_measurement_batch)
pheno$age2 = pheno$age^2
# match individual order with genotype file
ind = fread('height.ZBTB38.0.01.psam')
match.idx = match(ind$IID, pheno$id)
pheno = pheno[match.idx,]
Z = model.matrix(~ sex + age + age2 + assessment_centre + genotype_measurement_batch +
pc_genetic1 + pc_genetic2 + pc_genetic3 + pc_genetic4 + pc_genetic5 +
pc_genetic6 + pc_genetic7 + pc_genetic8 + pc_genetic9 + pc_genetic10 +
pc_genetic11 + pc_genetic12 + pc_genetic13 + pc_genetic14 + pc_genetic15 +
pc_genetic16 + pc_genetic17 + pc_genetic18 + pc_genetic19 + pc_genetic20 + X[,896], data = pheno)
# Remove intercept
Z = Z[,-1]
colnames(Z)[150] = 'X896'
Z = scale(Z, center=T, scale=F)
# standardize quantitative columns
cols = which(colnames(Z) %in% c("age","pc_genetic1","pc_genetic2","pc_genetic3","pc_genetic4",
"pc_genetic5","pc_genetic6","pc_genetic7","pc_genetic8","pc_genetic9",
"pc_genetic10","pc_genetic11","pc_genetic12","pc_genetic13","pc_genetic14",
"pc_genetic15","pc_genetic16","pc_genetic17","pc_genetic18","pc_genetic19","pc_genetic20", "X896"))
Z[,cols] = scale(Z[,cols])
Z[,'age2'] = Z[,'age']^2
# Compute XtX and Xty
y = pheno$height
names(y) = pheno$id
# Center y
y = y - mean(y)
# Center scale X
cm = Matrix::colMeans(X, na.rm = TRUE)
csd = susieR:::compute_colSds(X)
csd[csd == 0] = 1
X = as.matrix(t((t(X) - cm) / csd))
A <- crossprod(Z) # Z'Z
# chol decomposition for (Z'Z)^(-1)
R = chol(solve(A)) # R'R = (Z'Z)^(-1)
W = R %*% crossprod(Z, X) # RZ'X
S = R %*% crossprod(Z, y) # RZ'y
# Load LD matrix from raw genotype
ld.matrix = as.matrix(fread(paste0('height.ZBTB38.matrix')))
# X'X
XtX = ld.matrix*(nrow(X)-1) - crossprod(W) # W'W = X'ZR'RZ'X = X'Z(Z'Z)^{-1}Z'X
rownames(XtX) = colnames(XtX) = colnames(X)
# X'y
Xty = as.vector(y %*% X)
Xty = Xty - crossprod(W, S) # W'S = X'ZR'RZ'y = X'Z(Z'Z)^{-1}Z'y
saveRDS(list(XtX = XtX, Xty = Xty, yty = sum(y^2) - crossprod(S), n = length(y), pos=pos),
paste0('height.ZBTB38.removeCS2.XtX.Xty.rds'))
# # Save phenotype file
# pheno.table = data.frame(FID = pheno$id, IID = pheno$id, height = pheno$height)
# write.table(pheno.table, file = out.pheno.file, quote = FALSE, row.names = FALSE)
# # Save covariates file
# IID = pheno$id
# FID = pheno$id
# cov.table = cbind(FID, IID, Z)
# write.table(cov.table, file= out.covar.file, quote=FALSE, row.names = FALSE)plink2 --linear hide-covar no-x-sex omit-ref --pfile height.ZBTB38.0.01 --covar covar.removeCS2.height.txt --pheno pheno.height.txt --vif 100 --out height.ZBTB38.0.01.removeCS2.plink2After removing the effect of rs11919556 (p = 0.0213),
ss.dat = readRDS('../data/height.ZBTB38.removeCS2.XtX.Xty.rds')
betas = as.vector(ss.dat$Xty/diag(ss.dat$XtX))
rss = c(ss.dat$yty) - betas * ss.dat$Xty
se = as.vector(sqrt(rss/((ss.dat$n-1)*diag(ss.dat$XtX))))
z = betas/se
pval = 2*pnorm(-abs(z))
R = as.matrix(t(ss.dat$XtX * (1/sqrt(diag(ss.dat$XtX)))) * (1/ sqrt(diag(ss.dat$XtX))))
names(z) = rownames(R)
z.max = which.max(abs(z))
locus.zoom.cs(z, cs = mod_bhat$sets$cs$L1, pos=ss.dat$pos$POS/1e6, chr=3, gene.pos.map = gene.pos.map, z.ref.name = names(z)[z.max], ld = R[z.max,], xrange=c(140.9,141.5), title='ZBTB38 CS1', y.lab = '-log10(p value condition on CS2)')
# lm.res = read.delim('../data/height.ZBTB38.0.01.removeCS2.plink2.height.glm.linear')
# lm.res.t = lm.res$T_STAT
# names(lm.res.t) = lm.res$ID
# z.max = which.max(abs(lm.res.t))
# locus.zoom.cs(lm.res.t, cs = mod_bhat$sets$cs$L1, pos=lm.res$POS/1e6, chr=3, gene.pos.map = gene.pos.map, z.ref.name = lm.res$ID[z.max], ld = R[z.max,], xrange=c(140.9,141.5), title='ZBTB38 CS1', y.lab = '-log10(p value condition on CS2)')- Remove effect of all SNPs based on weights in other CSs
XtXr = mod_bhat$XtXr - ss.dat$XtX %*% (mod_bhat$alpha[1,] * mod_bhat$mu[1,])
Xtr = ss.dat$Xty - XtXr
betas = Xtr/diag(ss.dat$XtX)
b_2 = colSums(mod_bhat$alpha[-1,]*mod_bhat$mu[-1,])/mod_bhat$X_column_scale_factors
rss = c(ss.dat$yty - 2*sum(ss.dat$Xty * b_2) + sum(XtXr * b_2)) - betas * Xtr
se = as.vector(sqrt(rss/((ss.dat$n-1)*diag(ss.dat$XtX))))## Warning in sqrt(rss/((ss.dat$n - 1) * diag(ss.dat$XtX))): NaNs producedz.CS1 = betas/se
names(z.CS1) = colnames(R)
z.cs1.max = which.max(abs(z.CS1))
p2 = locus.zoom.cs(z.CS1, mod_bhat$sets$cs$L1, pos=ss.dat$pos$POS/1e6, chr = 3, gene.pos.map = gene.pos.map, z.ref.name = names(z.cs1.max), ld = R[z.cs1.max,], xrange=c(140.9,141.5), title='ZBTB38 Credible Set 1', y.lab = '-log10(p value condition on other CSs)', title.size = 20)## Warning: Removed 1232 rows containing missing values (geom_point).## Warning: Removed 5 rows containing missing values (geom_rect).## Warning: Removed 5 rows containing missing values (geom_text).p2
# pdf('GitHub/finemap-uk-biobank/ZBTB38_CS1.pdf', height=5, width=8)
# p2
# dev.off()tmp = data.frame(POS = ss.dat$pos$POS/1e6, logPO = log(mod_bhat$alpha[1,]/(1-mod_bhat$alpha[1,])))
tmp$CS = rep(4, length(z))
tmp$CS[mod_bhat$sets$cs$L1] = 16
tmp$CS = as.factor(tmp$CS)
pl_zoom = ggplot(tmp, aes(x = POS, y = logPO, shape = CS, size=CS)) +
geom_point() + ylab('log posterior odd') +
scale_shape_manual(values = c(4, 19), guide=FALSE) +
scale_size_manual(values=c(1.5,4), guide=FALSE) +
ggtitle('ZBTB38 Credible Set 1') +
theme_bw() + theme(axis.title.x=element_blank(), plot.title = element_text(size=12))
tmp.sub = tmp[mod_bhat$sets$cs$L1,]
pl_zoom = pl_zoom + geom_point(data = tmp.sub, aes(x=POS, y=logPO), shape=1, size=4, color='black', stroke=0.1)
pl_zoom = pl_zoom + xlim(140.9,141.5)
pl_gene = plot_geneName(gene.pos.map, xrange = c(140.9,141.5), chr=3)
g = egg::ggarrange(pl_zoom, pl_gene, nrow=2, heights = c(5,1), draw=FALSE)## Warning: Removed 1231 rows containing missing values (geom_point).## Warning: Removed 5 rows containing missing values (geom_rect).## Warning: Removed 5 rows containing missing values (geom_text).g
CS 2
- Removing effect of top SNP in CS 1
library(data.table)
library(readr)
library(Matrix)
geno.file = 'height.ZBTB38.raw.gz'
cat("Reading genotype data.\n")
geno <- fread(geno.file,sep = "\t",header = TRUE,stringsAsFactors = FALSE)
class(geno) <- "data.frame"
# Extract the genotypes.
X <- as(as.matrix(geno[-(1:6)]),'dgCMatrix')
pheno.file <- "/gpfs/data/stephens-lab/finemap-uk-biobank/data/raw/height.csv.gz"
out.pheno.file <- "pheno.height.txt"
out.covar.file <- "covar.removeCS1.height.txt"
# LOAD PHENOTYPE and COVARIATES DATA
# -------------------
# Read the phenotype data from the CSV file.
cat("Reading phenotype data.\n")
pheno <- suppressMessages(read_csv(pheno.file))
class(pheno) <- "data.frame"
pheno$sex = factor(pheno$sex)
pheno$assessment_centre = factor(pheno$assessment_centre)
pheno$genotype_measurement_batch = factor(pheno$genotype_measurement_batch)
pheno$age2 = pheno$age^2
# match individual order with genotype file
ind = fread('height.ZBTB38.psam')
match.idx = match(ind$IID, pheno$id)
pheno = pheno[match.idx,]
Z = model.matrix(~ sex + age + age2 + assessment_centre + genotype_measurement_batch +
pc_genetic1 + pc_genetic2 + pc_genetic3 + pc_genetic4 + pc_genetic5 +
pc_genetic6 + pc_genetic7 + pc_genetic8 + pc_genetic9 + pc_genetic10 +
pc_genetic11 + pc_genetic12 + pc_genetic13 + pc_genetic14 + pc_genetic15 +
pc_genetic16 + pc_genetic17 + pc_genetic18 + pc_genetic19 + pc_genetic20 + X[,865], data = pheno)
# Remove intercept
Z = Z[,-1]
colnames(Z)[150] = 'X865'
Z = scale(Z, center=T, scale = F)
# standardize quantitative columns
cols = which(colnames(Z) %in% c("age","pc_genetic1","pc_genetic2","pc_genetic3","pc_genetic4",
"pc_genetic5","pc_genetic6","pc_genetic7","pc_genetic8","pc_genetic9",
"pc_genetic10","pc_genetic11","pc_genetic12","pc_genetic13","pc_genetic14",
"pc_genetic15","pc_genetic16","pc_genetic17","pc_genetic18","pc_genetic19","pc_genetic20", "X865"))
Z[,cols] = scale(Z[,cols])
Z[,'age2'] = Z[,'age']^2
# Compute XtX and Xty
y = pheno$height
names(y) = pheno$id
# Center y
y = y - mean(y)
# Center scale X
cm = Matrix::colMeans(X, na.rm = TRUE)
csd = susieR:::compute_colSds(X)
csd[csd == 0] = 1
X = as.matrix(t((t(X) - cm) / csd))
A <- crossprod(Z) # Z'Z
# chol decomposition for (Z'Z)^(-1)
R = chol(solve(A)) # R'R = (Z'Z)^(-1)
W = R %*% crossprod(Z, X) # RZ'X
S = R %*% crossprod(Z, y) # RZ'y
# Load LD matrix from raw genotype
ld.matrix = as.matrix(fread(paste0('height.ZBTB38.matrix')))
# X'X
XtX = ld.matrix*(nrow(X)-1) - crossprod(W) # W'W = X'ZR'RZ'X = X'Z(Z'Z)^{-1}Z'X
rownames(XtX) = colnames(XtX) = colnames(X)
# X'y
Xty = as.vector(y %*% X)
Xty = Xty - crossprod(W, S) # W'S = X'ZR'RZ'y = X'Z(Z'Z)^{-1}Z'y
saveRDS(list(XtX = XtX, Xty = Xty, yty = sum(y^2) - crossprod(S), n = length(y), pos=pos),
paste0('height.ZBTB38.removeCS1.XtX.Xty.rds'))
# # Save phenotype file
# pheno.table = data.frame(FID = pheno$id, IID = pheno$id, height = pheno$height)
# write.table(pheno.table, file = out.pheno.file, quote = FALSE, row.names = FALSE)
# # Save covariates file
# IID = pheno$id
# FID = pheno$id
# cov.table = cbind(FID, IID, Z)
# write.table(cov.table, file= out.covar.file, quote=FALSE, row.names = FALSE)plink2 --linear hide-covar no-x-sex omit-ref --pfile height.ZBTB38.0.01 --covar covar.removeCS1.height.txt --pheno pheno.height.txt --vif 100 --out height.ZBTB38.0.01.removeCS1.plink2After removing the effect of rs2871960 (p = 7.200e-207), the conditional p value for rs11919556 becomes 4.279e-06, and it becomes the strongest one among all SNPs!
ss.dat = readRDS('../data/height.ZBTB38.removeCS1.XtX.Xty.rds')
betas = as.vector(ss.dat$Xty/diag(ss.dat$XtX))
rss = c(ss.dat$yty) - betas * ss.dat$Xty
se = as.vector(sqrt(rss/((ss.dat$n-1)*diag(ss.dat$XtX))))
z = betas/se
pval = 2*pnorm(-abs(z))
R = as.matrix(t(ss.dat$XtX * (1/sqrt(diag(ss.dat$XtX)))) * (1/ sqrt(diag(ss.dat$XtX))))
names(z) = rownames(R)
z.max = which.max(abs(z))
locus.zoom.cs(z, cs = mod_bhat$sets$cs$L2, pos=ss.dat$pos$POS/1e6, chr=3, gene.pos.map = gene.pos.map, z.ref.name = names(z)[z.max], ld = R[z.max,], xrange=c(140.9,141.5), title='ZBTB38 CS2', y.lab = '-log10(p value condition on CS1)')
# lm.res = read.delim('../data/height.ZBTB38.removeCS1.plink2.height.glm.linear')
# lm.res.t = lm.res$T_STAT
# names(lm.res.t) = lm.res$ID
# z.max = which.max(abs(lm.res.t))
# locus.zoom.cs(lm.res.t, cs = mod_bhat$sets$cs$L2, pos=lm.res$POS/1e6, chr=3, gene.pos.map = gene.pos.map, z.ref.name = lm.res$ID[z.max], ld = R[z.max,], xrange=c(140.9,141.5), title='ZBTB38 CS2', y.lab = '-log10(p value condition on CS1)')- Remove effect of all SNPs based on weights in other CSs
XtXr = mod_bhat$XtXr - ss.dat$XtX %*% (mod_bhat$alpha[2,] * mod_bhat$mu[2,])
Xtr = ss.dat$Xty - XtXr
betas = Xtr/diag(ss.dat$XtX)
b_2 = colSums(mod_bhat$alpha[-2,]*mod_bhat$mu[-2,])/mod_bhat$X_column_scale_factors
rss = c(ss.dat$yty - 2*sum(ss.dat$Xty * b_2) + sum(XtXr * b_2)) - betas * Xtr
se = as.vector(sqrt(rss/((ss.dat$n-1)*diag(ss.dat$XtX))))## Warning in sqrt(rss/((ss.dat$n - 1) * diag(ss.dat$XtX))): NaNs producedz.CS2 = betas/se
names(z.CS2) = colnames(R)
z.cs2.max = which.max(abs(z.CS2))
p3 = locus.zoom.cs(z.CS2, mod_bhat$sets$cs$L2, pos=ss.dat$pos$POS/1e6, chr=3, gene.pos.map = gene.pos.map, z.ref.name = names(z.cs2.max), ld = R[z.cs2.max,], xrange=c(140.9,141.5), title='ZBTB38 Credible Set 2', y.lab = '-log10(p value condition on other CSs)', title.size = 20)## Warning: Removed 1232 rows containing missing values (geom_point).## Warning: Removed 5 rows containing missing values (geom_rect).## Warning: Removed 5 rows containing missing values (geom_text).p3
# pdf('GitHub/finemap-uk-biobank/ZBTB38_CS2.pdf', height=5, width=8)
# p3
# dev.off()tmp = data.frame(POS = ss.dat$pos$POS/1e6, logPO = log(mod_bhat$alpha[2,]/(1-mod_bhat$alpha[2,])))
tmp$CS = rep(4, length(z))
tmp$CS[mod_bhat$sets$cs$L2] = 16
tmp$CS = as.factor(tmp$CS)
pl_zoom = ggplot(tmp, aes(x = POS, y = logPO, shape = CS, size=CS)) +
geom_point() + ylab('log posterior odd') +
scale_shape_manual(values = c(4, 19), guide=FALSE) +
scale_size_manual(values=c(1.5,4), guide=FALSE) +
ggtitle('ZBTB38 Credible Set 2') +
theme_bw() + theme(axis.title.x=element_blank(), plot.title = element_text(size=12))
tmp.sub = tmp[mod_bhat$sets$cs$L2,]
pl_zoom = pl_zoom + geom_point(data = tmp.sub, aes(x=POS, y=logPO), shape=1, size=4, color='black', stroke=0.1)
pl_zoom = pl_zoom + xlim(140.9,141.5)
pl_gene = plot_geneName(gene.pos.map, xrange = c(140.9,141.5), chr=3)
g = egg::ggarrange(pl_zoom, pl_gene, nrow=2, heights = c(5,1), draw=FALSE)## Warning: Removed 1231 rows containing missing values (geom_point).## Warning: Removed 5 rows containing missing values (geom_rect).## Warning: Removed 5 rows containing missing values (geom_text).g
Simulation under the estimated model
In the following simulation, we treat rs2871960 and rs11919556 as true signals. The effect sizes are from estimated model. We use the fitted residual variance in the simulation. The response y is simulated from \[ y \sim N_n(Xb, \sigma^2 I) \] , where X the genotype matrix that column centered, scaled, and removed the effect of covariates.
# ON CRI
library(data.table)
library(readr)
library(Matrix)
library(susieR)
geno.file = 'height.ZBTB38.0.01.raw.gz'
cat("Reading genotype data.\n")
geno <- fread(geno.file,sep = "\t",header = TRUE,stringsAsFactors = FALSE)
class(geno) <- "data.frame"
# Extract the genotypes.
X <- as(as.matrix(geno[-(1:6)]),'dgCMatrix')
pheno.file <- "/gpfs/data/stephens-lab/finemap-uk-biobank/data/raw/height.csv.gz"
# LOAD PHENOTYPE and COVARIATES DATA
# -------------------
# Read the phenotype data from the CSV file.
cat("Reading phenotype data.\n")
pheno <- suppressMessages(read_csv(pheno.file))
class(pheno) <- "data.frame"
pheno$sex = factor(pheno$sex)
pheno$assessment_centre = factor(pheno$assessment_centre)
pheno$genotype_measurement_batch = factor(pheno$genotype_measurement_batch)
pheno$age2 = pheno$age^2
# match individual order with genotype file
ind = fread('height.ZBTB38.0.01.psam')
match.idx = match(ind$IID, pheno$id)
pheno = pheno[match.idx,]
Z = model.matrix(~ sex + age + age2 + assessment_centre + genotype_measurement_batch +
pc_genetic1 + pc_genetic2 + pc_genetic3 + pc_genetic4 + pc_genetic5 +
pc_genetic6 + pc_genetic7 + pc_genetic8 + pc_genetic9 + pc_genetic10 +
pc_genetic11 + pc_genetic12 + pc_genetic13 + pc_genetic14 + pc_genetic15 +
pc_genetic16 + pc_genetic17 + pc_genetic18 + pc_genetic19 + pc_genetic20, data = pheno)
# Remove intercept
Z = Z[,-1]
# standardize quantitative columns
cols = which(colnames(Z) %in% c("age","pc_genetic1","pc_genetic2","pc_genetic3","pc_genetic4",
"pc_genetic5","pc_genetic6","pc_genetic7","pc_genetic8","pc_genetic9",
"pc_genetic10","pc_genetic11","pc_genetic12","pc_genetic13","pc_genetic14",
"pc_genetic15","pc_genetic16","pc_genetic17","pc_genetic18","pc_genetic19","pc_genetic20"))
Z[,cols] = scale(Z[,cols])
Z[,'age2'] = Z[,'age']^2
# Center scale X
cm = Matrix::colMeans(X, na.rm = TRUE)
csd = susieR:::compute_colSds(X)
csd[csd == 0] = 1
X = t((t(X) - cm) / csd)
A <- crossprod(Z)
# chol decomposition for (Z'Z)^(-1)
R = chol(solve(A))
W = R %*% t(Z) %*% X
# Remove Covariates from X
X <- as.matrix(X - Z %*% crossprod(R,W))
# Get estimated parameters
mod = readRDS('height.chr3.140800000.141800000.0.01.susie.model.rds')
ss.dat = readRDS('height.ZBTB38.0.01.XtX.Xty.rds')
R = as.matrix(t(ss.dat$XtX * (1/sqrt(diag(ss.dat$XtX)))) * (1/ sqrt(diag(ss.dat$XtX))))
sigma2 = mod$sigma2
V = mod$V[1:2]
n = nrow(X)
b2 = c(mod$mu[1, 945], mod$mu[2, 978])
Xb = as.vector(X[, c(945,978)] %*% b2)
result = vector("list", 1000)
set.seed(201910)
for(i in 1:1000){
## generate y
y = Xb + rnorm(n, 0, sqrt(sigma2))
## compute Xty, yty
Xty = as.vector(y %*% X)
yty = sum(y^2)
## compute summary stats
betas = as.vector(Xty/diag(ss.dat$XtX))
rss = yty - betas * Xty
se = as.vector(sqrt(rss/((n-1)*diag(ss.dat$XtX))))
z = betas/se
result[[i]] = susie_bhat(betas, se, R=R, n=n, var_y = as.numeric(yty/(n-1)))
}
saveRDS(result, 'height.ZBTB38.0.01.simulation1000.rds')Load simulation results:
result = readRDS('../data/height.ZBTB38.0.01.simulation1000.rds')In 1000 simulations, 655 runs have only 1 CS, 344 runs have 2 CSs.
table(sapply(result, function(mod) length(mod$sets$cs)))##
## 1 2 3
## 655 344 1323 runs contain both true signals.
contain.true = sapply(result, function(mod) all(c(945, 978) %in% unlist(mod$sets$cs)))
sum(contain.true)## [1] 323959 runs contain rs2871960.
contain.true.1 = sapply(result, function(mod) all( 945%in% unlist(mod$sets$cs)))
sum(contain.true.1)## [1] 959333 runs contain rs11919556.
contain.true.2 = sapply(result, function(mod) all( 978 %in% unlist(mod$sets$cs)))
sum(contain.true.2)## [1] 333