Estimate mQTL in purified cells. Marginal effects obtained per cell type. Can estimate heterogeneity between of mQTL effect estimates across cell types
Estimate G x celltype proportion interaction for each celltype. Interaction is a deviation fron mean effect. Quite a complex model Can we just estimate heterogeneity between interaction terms?
library(ggplot2)library(dplyr)
Attaching package: 'dplyr'
The following objects are masked from 'package:stats':
filter, lag
The following objects are masked from 'package:base':
intersect, setdiff, setequal, union
sim <-function(nc, n, betas =runif(nc, -2, 2)) { g <-rbinom(n, 2, 0.4) m <-sapply(1:nc, function(i) { g * betas[i] +rnorm(n) })# for each individual sample cell type proportions cellprop <-sapply(1:n, function(x) {a <-runif(nc); a/sum(a)}) %>%t()# weighted sum M <- (scale(m) * cellprop) %>% rowSums res <-lapply(1:nc, function(i) { o <-summary(lm(M ~ g * cellprop[,i]))tibble(cell=i, b=o$coef[4,1], se=o$coef[4,2], pval=o$coef[4,4], method="interaction") }) %>%bind_rows() res2 <-lapply(1:nc, function(i) { o <-summary(lm(m[,i] ~ g))tibble(cell=i, b=o$coef[2,1], se=o$coef[2,2], pval=o$coef[2,4], method="marginal") }) %>%bind_rows() betas <-tibble(cell=1:nc, b=betas, method="true")# cellprop <- as_tibble(cellprop)# m <- as_tibble(m)# names(m) <- paste0("m", 1:nc)# names(cellprop) <- paste0("cp", 1:nc)# cellprop$M <- M# cellprop$g <- g# cellprop <- bind_cols(cellprop, m)# return(tibble(int=res, marg=res2, betas))# heterogeneity het <-bind_rows(fixed_effects_meta_analysis(res$b, res$se) %>%as_tibble() %>%slice_head(n=1) %>%mutate(method="interaction"),fixed_effects_meta_analysis(res2$b, res2$se) %>%as_tibble() %>%slice_head(n=1) %>%mutate(method="marginal"), ) effs <-bind_rows(res, res2, betas)return(list(effs=effs, het=het))}
# Provide vector of interaction betas and their standard errorsfixed_effects_meta_analysis <-function(beta_vec, se_vec) { w <-1/ se_vec^2 beta <-sum(beta_vec * w, na.rm=T) /sum(w, na.rm=T) se <-sqrt(1/sum(w, na.rm=T)) pval <-pnorm(abs(beta / se), lower.tail =FALSE) Qj <- w * (beta-beta_vec)^2 Q <-sum(Qj, na.rm=T) Qdf <-sum(!is.na(beta_vec))-1if(Qdf ==0) Q <-0 Qjpval <-pchisq(Qj, 1, lower.tail=FALSE) Qpval <-pchisq(Q, Qdf, lower.tail=FALSE) pv <-pnorm(abs(beta_vec)/se_vec, lower.tail=FALSE) min_pv <-min(pv, na.rm=T)return(list(beta=beta, se=se, Q=Q, Qdf=Qdf, Qpval=Qpval, Qj=Qj, Qjpval=Qjpval, min_pv=min_pv))}
sim2 <-function(nc, n, betas =runif(nc, -2, 2)) { g <-rbinom(n, 2, 0.4) m <-sapply(1:nc, function(i) { g * betas[i] +rnorm(n) })# for each individual sample cell type proportions cellprop <-sapply(1:n, function(x) {a <-runif(nc); a/sum(a)}) %>%t()# weighted sum M <- (scale(m) * cellprop) %>% rowSums res <-lapply(1:nc, function(i) { o <-summary(lm(M ~ g * cellprop[,i]))tibble(cell=i, b=o$coef[4,1], se=o$coef[4,2], pval=o$coef[4,4], method="interaction") }) %>%bind_rows() res2 <-lapply(1:nc, function(i) { o <-summary(lm(m[,i] ~ g))tibble(cell=i, b=o$coef[2,1], se=o$coef[2,2], pval=o$coef[2,4], method="marginal") }) %>%bind_rows() betas <-tibble(cell=1:nc, b=betas, method="true") cellprop <-as_tibble(cellprop) m <-as_tibble(m)names(m) <-paste0("m", 1:nc)names(cellprop) <-paste0("cp", 1:nc) cellprop$M <- M cellprop$g <- g cellprop <-bind_cols(cellprop, m)return(cellprop)}ana <-function(d) { nc <-length(grep("cp", names(d)))# per cell estimates o1 <-lapply(1:nc, \(x) { f <-paste0("M ~ g * cp", x) o <-summary(lm(f, data=d))$coef n <-rownames(o) o <-as_tibble(o)names(o) <-c("b", "se", "tval", "pval") o$term <- n o$cell <- x o$mod <-"interaction" o }) %>%bind_rows()# cell int estimates o2 <-lapply(1:nc, \(x) { f <-paste0("m", x, " ~ g") o <-summary(lm(f, data=d))$coef n <-rownames(o) o <-as_tibble(o)names(o) <-c("b", "se", "tval", "pval") o$term <- n o$cell <- x o$mod <-"marginal" o }) %>%bind_rows()bind_rows(o1, o2)}d <-sim2(2, 1000)
Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
`.name_repair` is omitted as of tibble 2.0.0.
ℹ Using compatibility `.name_repair`.
