#Packages require(dplyr) require(tidyr) require(meta) require(data.table) require(TwoSampleMR) #Initial causal simulations #Input parameters n_humans <- 100000 n_genes <- 200 causal_genes <- 0.5 #Effect size matrix esm <- data.frame(gene_name = paste0('X', seq(1:n_genes))) #Effect of 1 SD unit increase in gene function on cellular exposure X esm <- esm %>% mutate(causal = rbinom(n_genes, 1, causal_genes)) %>% mutate(direction = rbinom(n_genes, 1, 0.5)) %>% mutate(direction = replace(direction, direction == 0, -1)) %>% mutate(e1 = 0) %>% mutate(e2 = runif(n_genes, 0.05, 0.15) * direction) %>% mutate(effect = case_when(causal == 0 ~ e1, causal == 1 ~ e2)) esm_out <- data.frame(effect = esm$effect) esm_out$gene <- 1:n_genes #Effect of experimental or LoF perturbation on gene function in SD units esm_out$perturb_rate_EXP <- runif(n_genes, -1.5, -1) esm_out$perturb_rate_LoF <- runif(n_genes, -1, -0.5) #MAF of burden esm_out <- esm_out %>% mutate(burden_freq = runif(nrow(esm_out), 0.005, 0.01)) #Effect of 1 SD unit increase in cellular exposure Y on organism-level outcome Z n_phen <- 100 #Loop over phenotypes for (k in 1:n_phen) { print(k) causal_effect <- runif(1, 0.2, 0.3) #Resimulate the effects of genes on cellular exposure esm_out <- esm_out %>% mutate(causal = rbinom(n_genes, 1, causal_genes)) %>% mutate(direction = rbinom(n_genes, 1, 0.5)) %>% mutate(direction = replace(direction, direction == 0, -1)) %>% mutate(e1 = 0) %>% mutate(e2 = runif(n_genes, 0.05, 0.15) * direction) %>% mutate(effect = case_when(causal == 0 ~ e1, causal == 1 ~ e2)) %>% mutate(balanced_pleiotropy_outcome = rnorm(nrow(esm_out), 0, 0.1)) %>% mutate(unbalanced_pleiotropy_outcome = rnorm(nrow(esm_out), 0.1, 0.1) * direction) %>% mutate(phenotypic_pleiotropy_outcome = runif(nrow(esm_out), 0.4, 0.6) * effect) %>% mutate(dir_incons = rbinom(nrow(esm_out), 1, 0.2)) %>% mutate(perturb_rate_LoF_dir = case_when(dir_incons == 0 ~ perturb_rate_LoF, dir_incons == 1 ~ perturb_rate_LoF * -1)) %>% select(effect, gene, perturb_rate_EXP, perturb_rate_LoF, perturb_rate_LoF_dir, burden_freq, balanced_pleiotropy_outcome, unbalanced_pleiotropy_outcome, phenotypic_pleiotropy_outcome) #Experimental simulation exposure_df <- data.frame() for(gene_use in 1:n_genes) { d1 <- esm_out %>% filter(gene == gene_use) %>% slice(rep(1:n(), each = 10000)) d2 <- d1 %>% mutate(Phen = rnorm(nrow(d1), 0, 1) + (effect * perturb_rate_EXP)) %>% select(gene, Phen) %>% mutate(perturb = 1) d3 <- d1 %>% mutate(Phen = rnorm(nrow(d1), 0, 1)) %>% select(gene, Phen) %>% mutate(perturb = 0) d4 <- rbind(d2, d3) lm1 <- lm(Phen ~ perturb, data = d4) temp <- c(gene_use, summary(lm1)$coefficients[2, 1], summary(lm1)$coefficients[2, 2], summary(lm1)$coefficients[2, 4]) exposure_df <- rbind(exposure_df, temp) } names(exposure_df) <- c("gene", "beta", "se", "pval") #Human phenotype simulation e1 <- data.frame(ID = rep(1:n_humans, n_genes), gene = rep(1:n_genes, 1, each = n_humans)) e1 <- e1 %>% inner_join(esm_out, by = "gene") %>% select(ID, gene, burden_freq) e1$gene_burden <- rbinom(n_humans * n_genes, 2, e1$burden_freq) e1 <- e1 %>% select(-burden_freq) e1_wide <- spread(e1, key = "gene", value = "gene_burden") e1_wide <- e1_wide %>% select(-ID) #Baseline model e1_wide_baseline <- as.data.frame(t(t(e1_wide) * (esm_out$effect * causal_effect * esm_out$perturb_rate_LoF))) e1_wide_baseline <- e1_wide_baseline %>% mutate(WeightedBurdenSum = rowSums(.[1:n_genes])) %>% mutate(Outcome = rnorm(nrow(e1_wide_baseline), 0, 1) + (WeightedBurdenSum)) #Baseline model with directional inconsistency (20% of pLoF in wrong direction) e1b_wide_baseline <- as.data.frame(t(t(e1_wide) * (esm_out$effect * causal_effect * esm_out$perturb_rate_LoF_dir))) e1b_wide_baseline <- e1b_wide_baseline %>% mutate(WeightedBurdenSum = rowSums(.[1:n_genes])) %>% mutate(Outcome_DC = rnorm(nrow(e1b_wide_baseline), 0, 1) + (WeightedBurdenSum)) #Balanced pleiotropy e1_wide_bp <- as.data.frame( t( t(e1_wide) * ((esm_out$effect * causal_effect) + esm_out$balanced_pleiotropy_outcome) * esm_out$perturb_rate_LoF ) ) e1_wide_bp <- e1_wide_bp %>% mutate(WeightedBurdenSum = rowSums(.[1:n_genes])) %>% mutate(Outcome_BP = rnorm(nrow(e1_wide_bp), 0, 1) + (WeightedBurdenSum)) #Unbalanced pleiotropy e1_wide_up <- as.data.frame( t( t(e1_wide) * ((esm_out$effect * causal_effect) + esm_out$unbalanced_pleiotropy_outcome) * esm_out$perturb_rate_LoF ) ) e1_wide_up <- e1_wide_up %>% mutate(WeightedBurdenSum = rowSums(.[1:n_genes])) %>% mutate(Outcome_UP = rnorm(nrow(e1_wide_up), 0, 1) + (WeightedBurdenSum)) #Phenotypic-level unbalanced pleiotropy e1_wide_pp <- as.data.frame( t( t(e1_wide) * ((esm_out$effect * causal_effect) + esm_out$phenotypic_pleiotropy_outcome) * esm_out$perturb_rate_LoF ) ) e1_wide_pp <- e1_wide_pp %>% mutate(WeightedBurdenSum = rowSums(.[1:n_genes])) %>% mutate(Outcome_PP = rnorm(nrow(e1_wide_pp), 0, 1) + (WeightedBurdenSum)) #Regr Dataset regr_dataset <- data.frame(Outcome = e1_wide_baseline$Outcome, Outcome_DC = e1b_wide_baseline$Outcome_DC, Outcome_BP = e1_wide_bp$Outcome_BP, Outcome_UP = e1_wide_up$Outcome_UP, Outcome_PP = e1_wide_pp$Outcome_PP, e1_wide) outcome_df <- data.frame() for (i in 1:n_genes) { temp <- regr_dataset %>% select(Outcome, Outcome_DC, Outcome_BP, Outcome_UP, Outcome_PP, i+5) names(temp) <- c("Outcome", "Outcome_DC", "Outcome_BP", "Outcome_UP", "Outcome_PP", "Gene") lm2 <- lm(Outcome ~ Gene, data = temp) lm3 <- lm(Outcome_BP ~ Gene, data = temp) lm4 <- lm(Outcome_UP ~ Gene, data = temp) lm5 <- lm(Outcome_PP ~ Gene, data = temp) lm6 <- lm(Outcome_DC ~ Gene, data = temp) temp <- c(i, summary(lm2)$coefficients[2, 1], summary(lm2)$coefficients[2, 2], summary(lm2)$coefficients[2, 4], summary(lm3)$coefficients[2, 1], summary(lm3)$coefficients[2, 2], summary(lm3)$coefficients[2, 4], summary(lm4)$coefficients[2, 1], summary(lm4)$coefficients[2, 2], summary(lm4)$coefficients[2, 4], summary(lm5)$coefficients[2, 1], summary(lm5)$coefficients[2, 2], summary(lm5)$coefficients[2, 4], summary(lm6)$coefficients[2, 1], summary(lm6)$coefficients[2, 2], summary(lm6)$coefficients[2, 4]) outcome_df <- rbind(outcome_df, temp) } names(outcome_df) <- c("gene", "beta_outcome_m1", "se_outcome_m1", "pval_outcome_m1", "beta_outcome_m2", "se_outcome_m2", "pval_outcome_m2", "beta_outcome_m3", "se_outcome_m3", "pval_outcome_m3", "beta_outcome_m4", "se_outcome_m4", "pval_outcome_m4", "beta_outcome_m5", "se_outcome_m5", "pval_outcome_m5") exposure_df$phenotype <- k outcome_df$phenotype <- k outcome_df$causal_effect <- causal_effect } #GPAT merge <- merge(exposure_df, outcome_df, by = c("gene", "phenotype")) output <- data.frame() #Loop over 100 simulations for (j in 1:100) { merge_use <- merge %>% filter(phenotype == j) %>% filter(pval < 0.005) %>% mutate(ratio_m1 = beta_outcome_m1 / beta) %>% mutate(ratio.se_m1 = se_outcome_m1/abs(beta)) %>% mutate(ratio_m2 = beta_outcome_m2 / beta) %>% mutate(ratio.se_m2 = se_outcome_m2/abs(beta)) %>% mutate(ratio_m3 = beta_outcome_m3 / beta) %>% mutate(ratio.se_m3 = se_outcome_m3/abs(beta)) %>% mutate(ratio_m4 = beta_outcome_m4 / beta) %>% mutate(ratio.se_m4 = se_outcome_m4/abs(beta)) %>% mutate(ratio_m5 = beta_outcome_m5 / beta) %>% mutate(ratio.se_m5 = se_outcome_m5 / abs(beta)) res_lof_m1 = meta::metagen(merge_use$ratio_m1, merge_use$ratio.se_m1) res_lof_m2 = meta::metagen(merge_use$ratio_m2, merge_use$ratio.se_m2) res_lof_m3 = meta::metagen(merge_use$ratio_m3, merge_use$ratio.se_m3) res_lof_m4 = meta::metagen(merge_use$ratio_m4, merge_use$ratio.se_m4) res_lof_m5 = meta::metagen(merge_use$ratio_m5, merge_use$ratio.se_m5) #MR Egger egger_m1 <- mr_egger_regression(b_exp=merge_use$beta, b_out=merge_use$beta_outcome_m1, se_exp=merge_use$se, se_out=merge_use$se_outcome_m1) egger_m2 <- mr_egger_regression(b_exp=merge_use$beta, b_out=merge_use$beta_outcome_m2, se_exp=merge_use$se, se_out=merge_use$se_outcome_m2) egger_m3 <- mr_egger_regression(b_exp=merge_use$beta, b_out=merge_use$beta_outcome_m3, se_exp=merge_use$se, se_out=merge_use$se_outcome_m3) egger_m4 <- mr_egger_regression(b_exp=merge_use$beta, b_out=merge_use$beta_outcome_m4, se_exp=merge_use$se, se_out=merge_use$se_outcome_m4) egger_m5 <- mr_egger_regression(b_exp=merge_use$beta, b_out=merge_use$beta_outcome_m5, se_exp=merge_use$se, se_out=merge_use$se_outcome_m5) causal_effect <- merge_use$causal_effect[1] temp_output <- cbind(j, causal_effect, res_lof_m1$TE.common, res_lof_m1$seTE.common, res_lof_m1$pval.common, res_lof_m2$TE.common, res_lof_m2$seTE.common, res_lof_m2$pval.common, res_lof_m3$TE.common, res_lof_m3$seTE.common, res_lof_m3$pval.common, res_lof_m4$TE.common, res_lof_m4$seTE.common, res_lof_m4$pval.common, res_lof_m5$TE.common, res_lof_m5$seTE.common, res_lof_m5$pval.common, egger_m1$b, egger_m1$se, egger_m1$pval, egger_m1$b_i, egger_m1$se_i, egger_m1$pval_i, egger_m2$b, egger_m2$se, egger_m2$pval, egger_m2$b_i, egger_m2$se_i, egger_m2$pval_i, egger_m3$b, egger_m3$se, egger_m3$pval, egger_m3$b_i, egger_m3$se_i, egger_m3$pval_i, egger_m4$b, egger_m4$se, egger_m4$pval, egger_m4$b_i, egger_m4$se_i, egger_m4$pval_i, egger_m5$b, egger_m5$se, egger_m5$pval, egger_m5$b_i, egger_m5$se_i, egger_m5$pval_i) output <- rbind(output, temp_output) } names(output) <- c("phenotype", "true_causal_effect", "ivw_beta_m1", "ivw_se_m1", "ivw_pval_m1", "ivw_beta_m2", "ivw_se_m2", "ivw_pval_m2", "ivw_beta_m3", "ivw_se_m3", "ivw_pval_m3", "ivw_beta_m4", "ivw_se_m4", "ivw_pval_m4", "ivw_beta_m5", "ivw_se_m5", "ivw_pval_m5", "egger_beta_m1", "egger_se_m1", "egger_pval_m1", "int_beta_m1", "int_se_m1", "int_pval_m1", "egger_beta_m2", "egger_se_m2", "egger_pval_m2", "int_beta_m2", "int_se_m2", "int_pval_m2", "egger_beta_m3", "egger_se_m3", "egger_pval_m3", "int_beta_m3", "int_se_m3", "int_pval_m3", "egger_beta_m4", "egger_se_m4", "egger_pval_m4", "int_beta_m4", "int_se_m4", "int_pval_m4", "egger_beta_m5", "egger_se_m5", "egger_pval_m5", "int_beta_m5", "int_se_m5", "int_pval_m5" )