OSCR

Convergent coexpression reveals shared biological mechanisms underlying common and rare variant risk in six neuropsychiatric disorders.

Code ↔ Paper

6 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 6 matches
  1. [1] § Methods › Heritability analysis for coexpression convergent genes ↔ example/run_BHR.R, lines 951–1069 · score 0.64 · baseline model, heritability explained, burden score, BHR, allele, variance
  2. [2] § Methods › Gene prioritization from GWAS and rare variant burden studies ↔ pops.py, lines 16–57 · score 0.57 · Polygenic Priority Score, PoPS, trained, gene annotation, chromosome, MAGMA
  3. [3] § Methods › Enrichment of evolutionarily constrained and functionally prioritized gene-sets ↔ R/gson.R, lines 289–426 · score 0.54 · access date, Targets Database, modules, matching, genes
  4. [4] § Methods › Gene prioritization from GWAS and rare variant burden studies ↔ R/gseAnalyzer.R, lines 1–32 · score 0.52 · pvalue cutoff, ranked genes, thresholds
  5. [5] § Methods › Functional annotation of convergent genes ↔ R/enricher.R, lines 83–115 · score 0.50 · cluster Profiler, Hochberg, FDR, BH, pvalue, enrichment
  6. [6] § Methods › Heritability analysis for coexpression convergent genes ↔ ldsc.py, lines 422–481 · score 0.50 · LD scores, kb, bin, window, partition, MAF

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R · 1,180 lines · 67 KB · MIT · 1 match

  1. library(bhr)
  2. source("correct_winners_curse.R")
  3. library(pracma)
  5. ########################################BHR basic estimates (h2, intercept) #############################
  7. path = "/Users/daniel/Desktop/rare_h2/ms/genebass_summary_statistics_null/"
  8. phenotype_file <- read.csv("~/rv_h2/reference_files/ms_phenotype_sheet.csv")
  9. traits = phenotype_file[phenotype_file$phenotype_rg == 1,"phenotype_key"]
  10. baseline_model <- read.table("~/rv_h2/reference_files/ms_baseline_oe5.txt")
  12. bhr_summary_statistic_names <- c("bhr_ms_gene_ss_400k_final_withnullburden_pLoF_low0_high1e-05_group1",
  13. "bhr_ms_gene_ss_400k_final_withnullburden_pLoF_low1e-05_high0.0001_group2",
  14. "bhr_ms_gene_ss_400k_final_withnullburden_pLoF_low0.0001_high0.001_group3",
  15. "bhr_ms_gene_ss_400k_final_withnullburden_synonymous_low0_high1e-05_group1",
  16. "bhr_ms_gene_ss_400k_final_withnullburden_synonymous_low1e-05_high0.0001_group2",
  17. "bhr_ms_gene_ss_400k_final_withnullburden_synonymous_low0.0001_high0.001_group3",
  18. "bhr_ms_gene_ss_400k_final_withnullburden_missense-benign_low0_high1e-05_group1",
  19. "bhr_ms_gene_ss_400k_final_withnullburden_missense-benign_low1e-05_high0.0001_group2",
  20. "bhr_ms_gene_ss_400k_final_withnullburden_missense-benign_low0.0001_high0.001_group3",
  21. "bhr_ms_gene_ss_400k_final_withnullburden_missense-notbenign_low0_high1e-05_group1",
  22. "bhr_ms_gene_ss_400k_final_withnullburden_missense-notbenign_low1e-05_high0.0001_group2",
  23. "bhr_ms_gene_ss_400k_final_withnullburden_missense-notbenign_low0.0001_high0.001_group3")
  25. results_holder <- matrix(data = NA, ncol = 23, nrow = length(bhr_summary_statistic_names) * length(traits))
  26. counter = 1
  27. for (ss in bhr_summary_statistic_names){
  28. summary_statistics <- readRDS(paste0(path,ss,".ms.munged.Rds"))
  29. for (trait in traits){
  30. n_bhr_trait <- head(summary_statistics[summary_statistics$phenotype_key == trait,"N"],1)
  31. print(c(trait, n_bhr_trait))
  32. output = BHR(mode = "univariate",
  33. trait1_sumstats = summary_statistics[summary_statistics$phenotype_key == trait,],
  34. annotations = list(baseline_model))
  36. results_holder[counter,] <- c(trait,
  37. ss,
  38. output$mixed_model$heritabilities[1,ncol(baseline_model)],
  39. output$mixed_model$heritabilities[2,ncol(baseline_model)],
  40. output$mixed_model$enrichments[1,1],
  41. output$mixed_model$enrichments[2,1],
  42. output$mixed_model$enrichments[1,2],
  43. output$mixed_model$enrichments[2,2],
  44. output$mixed_model$enrichments[1,3],
  45. output$mixed_model$enrichments[2,3],
  46. output$mixed_model$enrichments[1,4],
  47. output$mixed_model$enrichments[2,4],
  48. ((1-sum(output$mixed_model$fractions[1,]))/(1-sum(output$mixed_model$fraction_burden_score))),
  49. output$significant_genes$number_significant_genes,
  50. output$significant_genes$fraction_burdenh2_significant,
  51. output$significant_genes$fraction_burdenh2_significant_se,
  52. output$qc$intercept,
  53. output$qc$intercept_se,
  54. output$qc$attenuation_ratio,
  55. output$qc$attenuation_ratio_se,
  56. output$qc$lambda_gc,
  57. output$qc$lambda_gc_se,
  58. output$qc$mu_genome)
  60. print(paste0("Finished BHR estimate for ",trait," in summary statistic group ",ss))
  61. counter = counter + 1
  62. }
  63. }
  64. results_holder_df = as.data.frame(results_holder)
  65. results_holder_df[,3:23] <- sapply(results_holder_df[,3:23],as.numeric)
  66. colnames(results_holder_df) <- c("phenotype_key", "summary_statistic", "bhr_h2", "bhr_h2_se",
  67. "bhr_enrichment_oe1", "bhr_enrichment_oe1_se",
  68. "bhr_enrichment_oe2", "bhr_enrichment_oe2_se",
  69. "bhr_enrichment_oe3", "bhr_enrichment_oe3_se",
  70. "bhr_enrichment_oe4", "bhr_enrichment_oe4_se",
  71. "bhr_enrichment_oe5",
  72. "n_significant_genes","fraction_h2_significant_genes", "fraction_h2_significant_genes_se",
  73. "intercept", "intercept_se", "attenuation_ratio", "attenuation_ratio_se", "lambda_gc", "lambda_gc_se", "mu_genome")
  74. results_holder_df = merge(results_holder_df, phenotype_file[c("phenotype_key", "display_name", "phenocode", "n_bhr", "phenotype_core")], by = "phenotype_key")
  76. write.csv(results_holder_df, "~/rv_h2/outputs/bhr_output_h2.csv")
  79. ######################Basic BHR estimates, plus slope correction ##########################
  80. source("~/rv_h2/BHR.R")
  81. source("~/rv_h2/BHR_h2.R")
  82. source("~/rv_h2/randomeffects_jackknife.R")
  83. source("~/rv_h2/BHR_meta.R")
  85. path = "/Users/daniel/Desktop/rare_h2/ms/genebass_summary_statistics_null/"
  86. phenotype_file <- read.csv("~/rv_h2/reference_files/ms_phenotype_sheet.csv")
  87. traits = phenotype_file[phenotype_file$phenotype_rg == 1,"phenotype_key"]
  88. baseline_model <- read.table("~/rv_h2/reference_files/ms_baseline_oe5.txt")
  90. bhr_summary_statistic_names <- c("bhr_ms_gene_ss_400k_final_withnullburden_pLoF_low0_high1e-05_group1",
  91. "bhr_ms_gene_ss_400k_final_withnullburden_pLoF_low1e-05_high0.0001_group2",
  92. "bhr_ms_gene_ss_400k_final_withnullburden_pLoF_low0.0001_high0.001_group3",
  93. "bhr_ms_gene_ss_400k_final_withnullburden_synonymous_low0_high1e-05_group1",
  94. "bhr_ms_gene_ss_400k_final_withnullburden_synonymous_low1e-05_high0.0001_group2",
  95. "bhr_ms_gene_ss_400k_final_withnullburden_synonymous_low0.0001_high0.001_group3",
  96. "bhr_ms_gene_ss_400k_final_withnullburden_missense-benign_low0_high1e-05_group1",
  97. "bhr_ms_gene_ss_400k_final_withnullburden_missense-benign_low1e-05_high0.0001_group2",
  98. "bhr_ms_gene_ss_400k_final_withnullburden_missense-benign_low0.0001_high0.001_group3",
  99. "bhr_ms_gene_ss_400k_final_withnullburden_missense-notbenign_low0_high1e-05_group1",
  100. "bhr_ms_gene_ss_400k_final_withnullburden_missense-notbenign_low1e-05_high0.0001_group2",
  101. "bhr_ms_gene_ss_400k_final_withnullburden_missense-notbenign_low0.0001_high0.001_group3")
  103. results_holder <- matrix(data = NA, ncol = 23, nrow = length(bhr_summary_statistic_names) * length(traits))
  104. counter = 1
  105. for (ss in bhr_summary_statistic_names){
  106. summary_statistics <- readRDS(paste0(path,ss,".ms.munged.Rds"))
  107. for (trait in traits){
  108. n_bhr_trait <- head(summary_statistics[summary_statistics$phenotype_key == trait,"N"],1)
  109. print(c(trait, n_bhr_trait))
  110. output = BHR(mode = "univariate",
  111. trait1_sumstats = summary_statistics[summary_statistics$phenotype_key == trait,],
  112. annotations = list(baseline_model),
  113. slope_correction = 4.55087151/n_bhr_trait)
  115. results_holder[counter,] <- c(trait,
  116. ss,
  117. output$mixed_model$heritabilities[1,ncol(baseline_model)],
  118. output$mixed_model$heritabilities[2,ncol(baseline_model)],
  119. output$mixed_model$enrichments[1,1],
  120. output$mixed_model$enrichments[2,1],
  121. output$mixed_model$enrichments[1,2],
  122. output$mixed_model$enrichments[2,2],
  123. output$mixed_model$enrichments[1,3],
  124. output$mixed_model$enrichments[2,3],
  125. output$mixed_model$enrichments[1,4],
  126. output$mixed_model$enrichments[2,4],
  127. ((1-sum(output$mixed_model$fractions[1,]))/(1-sum(output$mixed_model$fraction_burden_score))),
  128. output$significant_genes$number_significant_genes,
  129. output$significant_genes$fraction_burdenh2_significant,
  130. output$significant_genes$fraction_burdenh2_significant_se,
  131. output$qc$intercept,
  132. output$qc$intercept_se,
  133. output$qc$attenuation_ratio,
  134. output$qc$attenuation_ratio_se,
  135. output$qc$lambda_gc,
  136. output$qc$lambda_gc_se,
  137. output$qc$mu_genome)
  139. print(paste0("Finished BHR estimate for ",trait," in summary statistic group ",ss))
  140. counter = counter + 1
  141. }
  142. }
  143. results_holder_df = as.data.frame(results_holder)
  144. results_holder_df[,3:23] <- sapply(results_holder_df[,3:23],as.numeric)
  145. colnames(results_holder_df) <- c("phenotype_key", "summary_statistic", "bhr_h2", "bhr_h2_se",
  146. "bhr_enrichment_oe1", "bhr_enrichment_oe1_se",
  147. "bhr_enrichment_oe2", "bhr_enrichment_oe2_se",
  148. "bhr_enrichment_oe3", "bhr_enrichment_oe3_se",
  149. "bhr_enrichment_oe4", "bhr_enrichment_oe4_se",
  150. "bhr_enrichment_oe5",
  151. "n_significant_genes","fraction_h2_significant_genes", "fraction_h2_significant_genes_se",
  152. "intercept", "intercept_se", "attenuation_ratio", "attenuation_ratio_se", "lambda_gc", "lambda_gc_se", "mu_genome")
  153. results_holder_df = merge(results_holder_df, phenotype_file[c("phenotype_key", "display_name", "phenocode", "n_bhr", "phenotype_core")], by = "phenotype_key")
  155. write.csv(results_holder_df, "~/rv_h2/outputs/bhr_output_h2_ldcorrection.csv")
  157. #############################################BHR gene set enrichments#####################
  158. source("~/rv_h2/BHR.R")
  159. source("~/rv_h2/BHR_h2.R")
  160. source("~/rv_h2/randomeffects_jackknife.R")
  161. source("~/rv_h2/BHR_meta.R")
  163. phenotype_file <- read.csv("~/rv_h2/reference_files/ms_phenotype_sheet.csv")
  164. gene_set_path = "/Users/daniel/Desktop/rare_h2/ms/gene_sets/gene_set_output/"
  165. bhr_gene_sets <- c("brain_GABAergic",
  166. "brain_Glutamatergic",
  167. "cosmic_all",
  168. "cosmic_oncogene",
  169. "cosmic_tsg",
  170. "ICA_cordblood_Erythroid",
  171. "ICA_cordblood_Megakaryocytes",
  172. "liver_Epithelial",
  173. "segblood",
  174. "segcortex",
  175. "segliver",
  176. "siggene_50NA",
  177. "siggene_2453NA",
  178. "siggene_3063NA",
  179. "siggene_21001NA",
  180. "siggene_30010NA",
  181. "siggene_30080NA",
  182. "siggene_30620NA",
  183. "siggene_30680NA",
  184. "siggene_30750NA",
  185. "siggene_30770NA",
  186. "siggene_30780NA",
  187. "siggene_50NA")
  189. summary_statistics <- readRDS("/Users/daniel/Desktop/rare_h2/ms/genebass_summary_statistics_null/bhr_ms_gene_ss_400k_final_withnullburden_pLoF_low0_high1e-05_group1.ms.munged.Rds")
  190. traits = phenotype_file[phenotype_file$phenotype_core == 1,"phenotype_key"]
  191. baseline_model <- read.table("~/rv_h2/reference_files/ms_baseline_oe5.txt")
  193. results_holder <- matrix(data = NA, ncol = 6, nrow = length(bhr_gene_sets) * length(traits))
  194. counter = 1
  195. for (bhr_gene_set in bhr_gene_sets){
  196. bhr_gene_set_annotation <- read.table(paste0(gene_set_path,"bhr_ms_" ,bhr_gene_set,".txt"), header = TRUE)
  197. for (trait in traits){
  198. output = BHR(mode = "univariate",
  199. trait1_sumstats = summary_statistics[summary_statistics$phenotype_key == trait,],
  200. annotations = list(baseline_model, bhr_gene_set_annotation))
  202. results_holder[counter,] <- c(trait,
  203. bhr_gene_set,
  204. output$mixed_model$fractions[1,ncol(baseline_model)],
  205. output$mixed_model$fractions[2,ncol(baseline_model)],
  206. output$mixed_model$enrichments[1,ncol(baseline_model)],
  207. output$mixed_model$enrichments[2,ncol(baseline_model)])
  208. counter = counter + 1
  209. print(paste0("Completed trait ",trait))
  210. }
  211. print(paste0("Completed gene set: ", bhr_gene_set))
  212. }
  213. results_holder_df = as.data.frame(results_holder)
  214. results_holder_df[,3:6] <- sapply(results_holder_df[,3:6],as.numeric)
  215. colnames(results_holder_df) <- c("phenotype_key", "gene_set", "fraction_h2", "fraction_h2_se", "enrichment", "enrichment_se")
  216. results_holder_df$enrichment_z <- (results_holder_df$enrichment - 1) / results_holder_df$enrichment_se
  217. results_holder_df = merge(results_holder_df, phenotype_file[c("phenotype_key", "display_name", "phenocode", "n_bhr", "phenotype_core")], by = "phenotype_key")
  218. results_holder_df_rename <- results_holder_df
  219. results_holder_df_rename$gene_set_display <- ifelse(results_holder_df_rename$gene_set == "brain_Glutamatergic", "Glutamatergic_neurons",
  220. ifelse(results_holder_df_rename$gene_set == "segliver", "Liver",
  221. ifelse(results_holder_df_rename$gene_set == "ICA_cordblood_Megakaryocytes", "Megakaryocytes",
  222. ifelse(results_holder_df_rename$gene_set == "cosmic_all", "Cancer_genes",
  223. ifelse(results_holder_df_rename$gene_set == "cosmic_oncogene", "Cancer_oncogenes",
  224. ifelse(results_holder_df_rename$gene_set == "segblood", "Whole_blood",
  225. ifelse(results_holder_df_rename$gene_set == "cosmic_tsg", "Cancer_TSG",
  226. ifelse(results_holder_df_rename$gene_set == "ICA_cordblood_Erythroid", "Erythrocyte",
  227. ifelse(results_holder_df_rename$gene_set == "brain_GABAergic", "GABAergic_neuron",
  228. ifelse(results_holder_df_rename$gene_set == "liver_Epithelial", "Hepatocyte",
  229. ifelse(results_holder_df_rename$gene_set == "segcortex", "Cortex", NA)))))))))))
  231. write.csv(results_holder_df_rename[c("gene_set", "gene_set_display", "phenotype_key", "display_name",
  232. "fraction_h2", "fraction_h2_se", "enrichment", "enrichment_se", "enrichment_z")], "~/rv_h2/outputs/bhr_output_geneset.csv")
  234. #################################################BHR aggregate for each trait, and for all traits############################################
  235. source("~/rv_h2/BHR.R")
  236. source("~/rv_h2/BHR_h2.R")
  237. source("~/rv_h2/randomeffects_jackknife.R")
  238. source("~/rv_h2/BHR_meta.R")
  240. path = "/Users/daniel/Desktop/rare_h2/ms/genebass_summary_statistics_null/"
  241. phenotype_file <- read.csv("~/rv_h2/reference_files/ms_phenotype_sheet.csv")
  242. traits = phenotype_file[phenotype_file$phenotype_rg == 1,"phenotype_key"]
  243. baseline_model <- read.table("~/rv_h2/reference_files/ms_baseline_oe5.txt")
  245. bhr_summary_statistic_names <- c("bhr_ms_gene_ss_400k_final_withnullburden_pLoF_low0_high1e-05_group1",
  246. "bhr_ms_gene_ss_400k_final_withnullburden_pLoF_low1e-05_high0.0001_group2",
  247. "bhr_ms_gene_ss_400k_final_withnullburden_pLoF_low0.0001_high0.001_group3",
  248. "bhr_ms_gene_ss_400k_final_withnullburden_missense-notbenign_low0_high1e-05_group1",
  249. "bhr_ms_gene_ss_400k_final_withnullburden_missense-notbenign_low1e-05_high0.0001_group2",
  250. "bhr_ms_gene_ss_400k_final_withnullburden_missense-notbenign_low0.0001_high0.001_group3")
  252. #Aggregate h2 for each trait for g1-g3 (ultra-rare + rare)
  253. ss_list <- list(readRDS(paste0(path, bhr_summary_statistic_names[1],".ms.munged.Rds")),
  254. readRDS(paste0(path, bhr_summary_statistic_names[2],".ms.munged.Rds")),
  255. readRDS(paste0(path, bhr_summary_statistic_names[3],".ms.munged.Rds")),
  256. readRDS(paste0(path, bhr_summary_statistic_names[4],".ms.munged.Rds")),
  257. readRDS(paste0(path, bhr_summary_statistic_names[5],".ms.munged.Rds")),
  258. readRDS(paste0(path, bhr_summary_statistic_names[6],".ms.munged.Rds")))
  260. results_holder <- matrix(data = NA, ncol = 3, nrow = length(traits))
  261. counter = 1
  262. for (trait in traits){
  263. output <- BHR(mode = 'aggregate',
  264. ss_list_trait1 = ss_list, trait_list = list(trait),
  265. annotations = baseline_model)
  266. results_holder[counter,] <- c(trait,
  267. output$aggregated_mixed_model_h2,
  268. output$aggregated_mixed_model_h2se)
  269. print(paste0("Completed ",trait))
  270. counter = counter + 1
  271. }
  272. results_holder_df = as.data.frame(results_holder)
  273. results_holder_df[,2:3] <- sapply(results_holder_df[,2:3],as.numeric)
  274. colnames(results_holder_df) <- c("phenotype_key", "aggregated_h2", "aggregate_h2_se")
  275. results_holder_df = merge(results_holder_df, phenotype_file[c("phenotype_key", "display_name")], by = "phenotype_key")
  276. results_holder_df$ss_group = "g13_lof_mis"
  278. write.csv(results_holder_df, "~/rv_h2/outputs/bhr_output_aggregated_per_trait_h2.csv")
  281. #Aggregate across all core traits
  282. results_holder <- matrix(data = NA, ncol = 3, nrow = 1)
  283. output <- BHR(mode = 'aggregate', ss_list_trait1 = ss_list,
  284. trait_list = phenotype_file[phenotype_file$phenotype_core == 1,"phenotype_key"],
  285. annotations = baseline_model)
  286. results_holder[1,] <- c("All_traits",
  287. output$aggregated_mixed_model_h2,
  288. output$aggregated_mixed_model_h2se)
  289. results_holder_df = as.data.frame(results_holder)
  290. results_holder_df[,2:3] <- sapply(results_holder_df[,2:3],as.numeric)
  291. colnames(results_holder_df) <- c("phenotype_key", "aggregated_h2", "aggregated_h2_se")
  292. write.csv(results_holder_df, "~/rv_h2/outputs/bhr_output_aggregated_all_traits_h2.csv")
  296. #################################################BHR aggregate for each trait, and for all traits, with slope correction ############################################
  297. source("~/rv_h2/BHR.R")
  298. source("~/rv_h2/BHR_h2.R")
  299. source("~/rv_h2/randomeffects_jackknife.R")
  300. source("~/rv_h2/BHR_meta.R")
  302. path = "/Users/daniel/Desktop/rare_h2/ms/genebass_summary_statistics_null/"
  303. phenotype_file <- read.csv("~/rv_h2/reference_files/ms_phenotype_sheet.csv")
  304. traits = phenotype_file[phenotype_file$phenotype_rg == 1,"phenotype_key"]
  305. baseline_model <- read.table("~/rv_h2/reference_files/ms_baseline_oe5.txt")
  307. bhr_summary_statistic_names <- c("bhr_ms_gene_ss_400k_final_withnullburden_pLoF_low0_high1e-05_group1",
  308. "bhr_ms_gene_ss_400k_final_withnullburden_pLoF_low1e-05_high0.0001_group2",
  309. "bhr_ms_gene_ss_400k_final_withnullburden_pLoF_low0.0001_high0.001_group3",
  310. "bhr_ms_gene_ss_400k_final_withnullburden_missense-notbenign_low0_high1e-05_group1",
  311. "bhr_ms_gene_ss_400k_final_withnullburden_missense-notbenign_low1e-05_high0.0001_group2",
  312. "bhr_ms_gene_ss_400k_final_withnullburden_missense-notbenign_low0.0001_high0.001_group3")
  314. #Aggregate h2 for each trait for g1-g3 (ultra-rare + rare)
  315. ss_list <- list(readRDS(paste0(path, bhr_summary_statistic_names[1],".ms.munged.Rds")),
  316. readRDS(paste0(path, bhr_summary_statistic_names[2],".ms.munged.Rds")),
  317. readRDS(paste0(path, bhr_summary_statistic_names[3],".ms.munged.Rds")),
  318. readRDS(paste0(path, bhr_summary_statistic_names[4],".ms.munged.Rds")),
  319. readRDS(paste0(path, bhr_summary_statistic_names[5],".ms.munged.Rds")),
  320. readRDS(paste0(path, bhr_summary_statistic_names[6],".ms.munged.Rds")))
  322. results_holder <- matrix(data = NA, ncol = 3, nrow = length(traits))
  323. counter = 1
  324. summary_statistics <- readRDS(paste0(path,"bhr_ms_gene_ss_400k_final_withnullburden_pLoF_low0_high1e-05_group1.ms.munged.Rds"))
  325. for (trait in traits){
  326. n_bhr_trait <- head(summary_statistics[summary_statistics$phenotype_key == trait,"N"],1)
  327. output <- BHR(mode = 'aggregate',
  328. ss_list_trait1 = ss_list, trait_list = list(trait),
  329. annotations = baseline_model,
  330. slope_correction = 4.55087151/n_bhr_trait)
  331. results_holder[counter,] <- c(trait,
  332. output$aggregated_mixed_model_h2,
  333. output$aggregated_mixed_model_h2se)
  334. print(paste0("Completed ",trait))
  335. counter = counter + 1
  336. }
  337. results_holder_df = as.data.frame(results_holder)
  338. results_holder_df[,2:3] <- sapply(results_holder_df[,2:3],as.numeric)
  339. colnames(results_holder_df) <- c("phenotype_key", "aggregated_h2", "aggregate_h2_se")
  340. results_holder_df = merge(results_holder_df, phenotype_file[c("phenotype_key", "display_name")], by = "phenotype_key")
  341. results_holder_df$ss_group = "g13_lof_mis"
  343. write.csv(results_holder_df, "~/rv_h2/outputs/bhr_output_aggregated_per_trait_h2_ldcorrection.csv")
  346. #Aggregate across all core traits
  347. source("~/rv_h2/BHR.R")
  348. source("~/rv_h2/BHR_h2.R")
  349. source("~/rv_h2/randomeffects_jackknife.R")
  350. source("~/rv_h2/BHR_meta.R")
  352. path = "/Users/daniel/Desktop/rare_h2/ms/genebass_summary_statistics_null/"
  353. phenotype_file <- read.csv("~/rv_h2/reference_files/ms_phenotype_sheet.csv")
  354. baseline_model <- read.table("~/rv_h2/reference_files/ms_baseline_oe5.txt")
  356. bhr_summary_statistic_names <- c("bhr_ms_gene_ss_400k_final_withnullburden_pLoF_low0_high1e-05_group1",
  357. "bhr_ms_gene_ss_400k_final_withnullburden_pLoF_low1e-05_high0.0001_group2",
  358. "bhr_ms_gene_ss_400k_final_withnullburden_pLoF_low0.0001_high0.001_group3",
  359. "bhr_ms_gene_ss_400k_final_withnullburden_missense-notbenign_low0_high1e-05_group1",
  360. "bhr_ms_gene_ss_400k_final_withnullburden_missense-notbenign_low1e-05_high0.0001_group2",
  361. "bhr_ms_gene_ss_400k_final_withnullburden_missense-notbenign_low0.0001_high0.001_group3")
  363. #Aggregate h2 for each trait for g1-g3 (ultra-rare + rare)
  364. ss_list <- list(readRDS(paste0(path, bhr_summary_statistic_names[1],".ms.munged.Rds")),
  365. readRDS(paste0(path, bhr_summary_statistic_names[2],".ms.munged.Rds")),
  366. readRDS(paste0(path, bhr_summary_statistic_names[3],".ms.munged.Rds")),
  367. readRDS(paste0(path, bhr_summary_statistic_names[4],".ms.munged.Rds")),
  368. readRDS(paste0(path, bhr_summary_statistic_names[5],".ms.munged.Rds")),
  369. readRDS(paste0(path, bhr_summary_statistic_names[6],".ms.munged.Rds")))
  370. #remember to add slope line to meta for this
  371. results_holder <- matrix(data = NA, ncol = 3, nrow = 1)
  372. output <- BHR(mode = 'aggregate', ss_list_trait1 = ss_list,
  373. trait_list = phenotype_file[phenotype_file$phenotype_core == 1,"phenotype_key"],
  374. annotations = baseline_model)
  375. results_holder[1,] <- c("All_traits",
  376. output$aggregated_mixed_model_h2,
  377. output$aggregated_mixed_model_h2se)
  378. results_holder_df = as.data.frame(results_holder)
  379. results_holder_df[,2:3] <- sapply(results_holder_df[,2:3],as.numeric)
  380. colnames(results_holder_df) <- c("phenotype_key", "aggregated_h2", "aggregated_h2_se")
  381. write.csv(results_holder_df, "~/rv_h2/outputs/bhr_output_aggregated_all_traits_h2_ldcorrection.csv")
  387. ####Significant gene fractions
  389. phenotype_file <- read.csv("~/rv_h2/reference_files/ms_phenotype_sheet.csv")
  390. sig_genes = read.csv("~/rv_h2/reference_files/siggenes_for_AMM.csv")
  391. traits = phenotype_file[phenotype_file$phenotype_core == 1,"phenotype_key"]
  392. summary_statistics_plof_grp1 <- readRDS("~/Documents/oconnor_rotation/rarevariantproject/bhr_ms_gene_ss_400k_withnullburden_pLoF_nvar449780_low0_high1e-05_group1.ms.munged.Rds")
  393. baseline_model <- read.table("~/rv_h2/reference_files/ms_baseline_oe5.txt")
  394. consensus_genes <- read.table("~/rv_h2/reference_files/bhr_ms_consensus_gene_list.txt")[c("gene_id", "gene")]
  396. get_frac_assns <- function(trait){
  397. print(trait)
  398. trait_sumstats = summary_statistics_plof_grp1[summary_statistics_plof_grp1$phenotype_key == trait,]
  400. sig_genes_trait = sig_genes$gene[sig_genes$id == trait]
  401. print(length(sig_genes_trait))
  402. output = BHR(mode = "univariate",
  403. trait1_sumstats = trait_sumstats,
  404. annotations = list(baseline_model),
  405. fixed_genes = sig_genes_trait,
  406. gwc_exclusion = FALSE)
  408. trait_sumstats_sig = trait_sumstats[trait_sumstats$gene %in% sig_genes_trait,]
  410. trait_sig_df <- data.frame(phenotype_key = trait,
  411. gene = trait_sumstats_sig$gene,
  412. varexplained = trait_sumstats_sig$w_t_beta^2/trait_sumstats_sig$burden_score,
  413. bhr_h2 = output$mixed_model$heritabilities[1,ncol(output$mixed_model$heritabilities)],
  414. frac_sig = output$significant_genes$fraction_burdenh2_significant,
  415. frac_sig_se = output$significant_genes$fraction_burdenh2_significant_se)
  416. trait_sig_df$chisq = trait_sumstats$N[1]*trait_sig_df$varexplained
  418. thresh = qchisq(p = 0.05/nrow(trait_sumstats),df = 1,lower.tail = FALSE)
  419. trait_sig_df$chisq_winnerscursecorr = correct_winners_curse(trait_sig_df$chisq,thresh)
  420. trait_sig_df$varexplained_winnerscursecorr = trait_sig_df$chisq_winnerscursecorr/trait_sumstats$N[1]
  421. trait_sig_df$frac_sig_winnerscursecorr = sum(trait_sig_df$varexplained_winnerscursecorr)/trait_sig_df$bhr_h2[1]
  422. return(trait_sig_df)
  423. }
  425. sig_results <- lapply(traits[traits %in% sig_genes$id],get_frac_assns)
  427. sig_df = sig_results[[1]]
  428. for (trait in 2:length(sig_results)){
  429. print(trait)
  430. sig_df = rbind(sig_df,sig_results[[trait]])
  431. }
  433. sig_df$display_name = phenotype_file$display_name[match(sig_df$phenotype_key,phenotype_file$phenotype_key)]
  434. sig_df$proportion_bhr_explained = sig_df$varexplained_winnerscursecorr/sig_df$bhr_h2
  435. sig_df$labelgene = consensus_genes$gene[match(sig_df$gene,consensus_genes$gene_id)]
  437. ordered_sigdf = sig_df[order(sig_df$phenotype_key,-sig_df$varexplained_winnerscursecorr),]
  438. write.csv(ordered_sigdf, "~/rv_h2/outputs/BHR_Significant_Genes_Info.csv", quote = FALSE, row.names = FALSE)
  441. ####Significant gene fractions, common variant space
  443. variant_table <- fread2("~/Documents/oconnor_rotation/rarevariantproject/variants_MAFfilt_subsetcols.tsv", header = TRUE)
  444. sumstat_lookup <- read.csv("reference_files/sumstat_lookup.csv", header = FALSE)
  446. sigclump_df = data.frame()
  447. for (trait in 1:nrow(sumstat_lookup)){
  448. print(trait)
  449. print(sumstat_lookup$V21[trait])
  451. print("loading files")
  452. clumpfile = paste0("~/Documents/oconnor_rotation/rarevariantproject/final_commonvar_manuscript/significant_clump/clump_output/", strsplit(sumstat_lookup$V21[trait],split = ".bgz|.gz")[[1]],".plink.tsv_clump.clumped")
  453. gwas_file = paste0("~/Documents/oconnor_rotation/rarevariantproject/final_commonvar_manuscript/significant_clump/gwas_sumstats/",strsplit(sumstat_lookup$V10[trait],split = ".bgz|.gz")[[1]],"_lean")
  454. ldsc_file = paste0("~/Documents/oconnor_rotation/rarevariantproject/final_commonvar_manuscript/significant_clump/input_sumstats/",strsplit(sumstat_lookup$V7[trait],split = ".bgz|.gz")[[1]])
  457. clumps = read.table(clumpfile, header = TRUE)
  458. clumps = clumps[clumps$P < 5e-8,]
  460. gwas = fread2(gwas_file, header = FALSE)
  461. ldsc = read.table(ldsc_file, header = TRUE)
  462. names(gwas) <- c("variant","minor_AF","beta")
  464. print("extracting significant clumps variances")
  465. snps = clumps$SNP
  466. chrs = variant_table$chr[match(snps,variant_table$rsid)]
  467. bps = variant_table$pos[match(snps,variant_table$rsid)]
  468. refs = variant_table$ref[match(snps,variant_table$rsid)]
  469. alts = variant_table$alt[match(snps,variant_table$rsid)]
  471. identifiers = paste(chrs,bps,refs,alts,sep = ":")
  473. betas = gwas$beta[match(identifiers,gwas$variant)]
  474. mafs = gwas$minor_AF[match(identifiers,gwas$variant)]
  475. varexplained = (betas*(sqrt(2*mafs*(1 - mafs))))^2
  477. traitclumpdf <- data.frame(trait = sumstat_lookup$V21[trait],
  478. SNP = snps,
  479. varexplained = varexplained,
  480. N = ldsc$N[1])
  481. traitclumpdf = traitclumpdf[!is.na(traitclumpdf$varexplained),]
  482. traitclumpdf$chisq = traitclumpdf$varexplained * ldsc$N[1]
  484. thresh = qchisq(5e-8,1,lower.tail = FALSE)
  486. traitclumpdf$chisq_winnerscursecorr = correct_winners_curse(traitclumpdf$chisq,thresh)
  487. traitclumpdf$varexplained_winnerscursecorr = traitclumpdf$chisq_winnerscursecorr/ldsc$N[1]
  489. sigclump_df <- rbind(sigclump_df,traitclumpdf)
  491. }
  493. #aside: get Ns
  495. get_N <- function(trait){
  496. print(trait)
  497. print(sumstat_lookup$V21[trait])
  499. ldsc_file = paste0("~/Documents/oconnor_rotation/rarevariantproject/final_commonvar_manuscript/significant_clump/input_sumstats/",strsplit(sumstat_lookup$V7[trait],split = ".bgz|.gz")[[1]])
  501. ldsc = read.table(ldsc_file, header = TRUE)
  503. return(ldsc$N[1])
  505. }
  506. ldsc_h2 = read.csv("~/Documents/oconnor_rotation/rarevariantproject/rv_h2/ldsc_h2.csv")
  507. ldsc_h2$N = sapply(1:nrow(ldsc_h2), get_N)
  508. write.csv(ldsc_h2, "outputs//ldsc_h2.csv", row.names = FALSE, quote = FALSE)
  509. sigclump_df <- sigclump_df[sigclump_df$trait %in% phenotype_file$phenotype_key[phenotype_file$phenotype_core ==1],]
  510. sigclump_df <- sigclump_df[!is.na(sigclump_df$varexplained),]
  511. get_HESSh2 <- function(trait){
  512. print(trait)
  513. HESS_h2_result = read.table(paste0("~/Documents/oconnor_rotation/rarevariantproject/final_commonvar_manuscript/HESS/step2_output/", trait,"_HESSformat.tsv_step2.log_HESSh2"), sep = " ")
  515. return(as.numeric(HESS_h2_result$V5))
  516. }
  518. get_HESSh2_se <- function(trait){
  519. HESS_h2_result = read.table(paste0("~/Documents/oconnor_rotation/rarevariantproject/final_commonvar_manuscript/HESS/step2_output/", trait,"_HESSformat.tsv_step2.log_HESSh2"), sep = " ")
  521. return(parse_number(HESS_h2_result$V6))
  522. }
  524. sigclump_df$HESS_h2 = sapply(sigclump_df$trait,get_HESSh2)
  525. sigclump_df$HESS_h2_se = sapply(sigclump_df$trait,get_HESSh2_se)
  526. sigclump_df$fraction_HESS = sigclump_df$varexplained_winnerscursecorr/sigclump_df$HESS_h2
  528. sigclump_df_ordered = sigclump_df[order(sigclump_df$trait,-sigclump_df$varexplained_winnerscursecorr),]
  529. names(sigclump_df_ordered)[1] <- "phenotype_key"
  531. sigclump_df_ordered$fraction_HESS_sig = sapply(sigclump_df_ordered$phenotype_key,
  532. function(x) {
  533. return(sum(sigclump_df_ordered$varexplained_winnerscursecorr[sigclump_df_ordered$phenotype_key == x])/(sigclump_df_ordered$HESS_h2[sigclump_df_ordered$phenotype_key == x][1]))
  534. })
  536. write.csv(sigclump_df_ordered,"outputs/common_significant_associations.csv", quote = FALSE, row.names = FALSE)
  538. ####Significant gene fractions, common variant space (HESS partitions)
  539. HESS_df <- data.frame()
  541. for (trait in phenotype_file$phenotype_key[phenotype_file$phenotype_core == 1]){
  542. #print(trait)
  543. trait_hess_results <- read.table(paste0("~/Documents/oconnor_rotation/rarevariantproject/final_commonvar_manuscript/HESS/step2_output/",
  544. trait,
  545. "_HESSformat.tsv_step2.txt"),
  546. header = TRUE)
  548. trait_hess_results$trait = trait
  549. HESS_df = rbind(HESS_df,trait_hess_results)
  551. }
  553. write.csv(HESS_df,"outputs/HESS_output.csv", row.names = FALSE, quote = FALSE)
  554. ######################################## Genetic Correlation #############################
  555. #Read in results from cluster
  557. ldsc_rg_guide <- read.table("~/Documents/oconnor_rotation/rarevariantproject/final_commonvar_manuscript/LDSC_rg/ldsc_rg_guide.tsv", sep = ":")
  559. get_ldsc_rg <- function(pair){
  560. input = read.table(paste0("~/Documents/oconnor_rotation/rarevariantproject/final_commonvar_manuscript/LDSC_rg/output/",pair,".txt"), sep = " ")
  561. return(as.numeric(input$V3[1]))
  562. }
  564. get_ldsc_rg_se <- function(pair){
  565. input = read.table(paste0("~/Documents/oconnor_rotation/rarevariantproject/final_commonvar_manuscript/LDSC_rg/output/",pair,".txt"), sep = " ")
  566. return(parse_number(input$V4[1]))
  567. }
  569. ldsc_rg_guide$ldsc_rg = sapply(1:nrow(ldsc_rg_guide), get_ldsc_rg)
  570. ldsc_rg_guide$ldsc_rg_se = sapply(1:nrow(ldsc_rg_guide), get_ldsc_rg_se)
  573. get_bhr_rg <- function(pair){
  574. print(pair)
  575. input = read.table(paste0("~/Documents/oconnor_rotation/rarevariantproject/BHR_rg_cluster/output/",pair,".txt"), header = TRUE)
  576. return(input$x[1])
  577. }
  579. get_bhr_rg_se <- function(pair){
  580. print(pair)
  581. input = read.table(paste0("~/Documents/oconnor_rotation/rarevariantproject/BHR_rg_cluster/output/",pair,".txt"), header = TRUE)
  582. return(input$x[2])
  583. }
  585. ldsc_rg_guide$bhr_rg = sapply(1:nrow(ldsc_rg_guide), get_bhr_rg)
  586. ldsc_rg_guide$bhr_rg_se = sapply(1:nrow(ldsc_rg_guide), get_bhr_rg_se)
  588. bhr_h2 = data.frame(googlesheets4::read_sheet("https://docs.google.com/spreadsheets/d/1AwvDRKEUJ6EtUbvnvSV34JkPV0Y7ozJHw-vDJOW2PpY/edit#gid=881218879", sheet = "bhr_h2"))
  589. bhr_h2_plofgrp1 =bhr_h2[bhr_h2$summary_statistic == "bhr_ms_gene_ss_400k_final_withnullburden_pLoF_low0_high1e-05_group1",]
  591. ldsc_rg_guide$bhr_h2_trait1 = bhr_h2_plofgrp1$bhr_h2[match(ldsc_rg_guide$V3,bhr_h2_plofgrp1$phenotype_key)]
  592. ldsc_rg_guide$bhr_h2_trait1_se = bhr_h2_plofgrp1$bhr_h2_se[match(ldsc_rg_guide$V3,bhr_h2_plofgrp1$phenotype_key)]
  593. ldsc_rg_guide$bhr_h2_trait2 = bhr_h2_plofgrp1$bhr_h2[match(ldsc_rg_guide$V6,bhr_h2_plofgrp1$phenotype_key)]
  594. ldsc_rg_guide$bhr_h2_trait2_se = bhr_h2_plofgrp1$bhr_h2_se[match(ldsc_rg_guide$V6,bhr_h2_plofgrp1$phenotype_key)]
  596. ldsc_rg_guide$trait1_bhrh2_z = ldsc_rg_guide$bhr_h2_trait1/ldsc_rg_guide$bhr_h2_trait1_se
  597. ldsc_rg_guide$trait2_bhrh2_z = ldsc_rg_guide$bhr_h2_trait2/ldsc_rg_guide$bhr_h2_trait2_se
  599. rg_output = ldsc_rg_guide[,c("V3","V6","ldsc_rg","ldsc_rg_se","bhr_rg","bhr_rg_se","trait1_bhrh2_z","trait2_bhrh2_z")]
  600. rg_output_sigbhrh2 = rg_output[abs(rg_output$trait1_bhrh2_z) > 1.96 & abs(rg_output$trait2_bhrh2_z) > 1.96,]
  601. names(rg_output_sigbhrh2) <- c("trait1","trait2","ldsc_rg","ldsc_rg_se","bhr_rg","bhr_rg_se")
  602. names(rg_output) <- c("trait1","trait2","ldsc_rg","ldsc_rg_se","bhr_rg","bhr_rg_se")
  605. write.csv(rg_output,"outputs/rg_output.csv", quote = FALSE, row.names = FALSE)
  607. #rg between missense and plof
  609. plof_sumstats = readRDS("~/Documents/oconnor_rotation/rarevariantproject/bhr_ms_gene_ss_400k_withnullburden_pLoF_nvar449780_low0_high1e-05_group1.ms.munged.Rds")
  610. missense_sumstats = readRDS("~/Documents/oconnor_rotation/rarevariantproject/bhr_ms_gene_ss_400k_final_withnullburden_missense-notbenign_nvar1590551_low0_high1e-05_group1.ms.munged.Rds")
  612. traits = phenotype_file$phenotype_key[phenotype_file$phenotype_core == 1]
  614. get_plof_missense_stats <- function(trait){
  615. plof_sumstats_trait = plof_sumstats[plof_sumstats$phenotype_key == trait,]
  616. missense_sumstats_trait = missense_sumstats[missense_sumstats$phenotype_key == trait,]
  618. bhr_plof = BHR_h2(plof_sumstats_trait,
  619. annotations = list(baseline_model),
  620. num_blocks = 100,
  621. genomewide_correction = FALSE,
  622. gwc_exclusion = NULL,
  623. overdispersion = FALSE,
  624. num_null_conditions = 0,
  625. output_jackknife_h2 = FALSE,
  626. fixed_genes = NULL,
  627. all_models = FALSE,
  628. slope_correction = FALSE )
  630. bhr_missense = BHR_h2(missense_sumstats_trait,
  631. annotations = list(baseline_model),
  632. num_blocks = 100,
  633. genomewide_correction = FALSE,
  634. gwc_exclusion = NULL,
  635. overdispersion = FALSE,
  636. num_null_conditions = 0,
  637. output_jackknife_h2 = FALSE,
  638. fixed_genes = NULL,
  639. all_models = FALSE,
  640. slope_correction = FALSE )
  641. if (bhr_plof$mixed_model$heritabilities[1,5] < 0 | bhr_missense$mixed_model$heritabilities[1,5] < 0 ){
  642. output = list(plof_h2 = bhr_plof$mixed_model$heritabilities[1,5],
  643. plof_h2_se = bhr_plof$mixed_model$heritabilities[2,5],
  644. missense_h2 = bhr_missense$mixed_model$heritabilities[1,5],
  645. missense_h2_se = bhr_missense$mixed_model$heritabilities[2,5],
  646. rg = NA,
  647. rg_se = NA)
  648. print(output)
  649. return(output)
  650. }
  652. bhr_rg_trait = BHR_rg(plof_sumstats_trait,
  653. missense_sumstats_trait,
  654. annotations = list(baseline_model),
  655. num_blocks = 100,
  656. genomewide_correction = FALSE,
  657. overdispersion = FALSE,
  658. num_null_conditions = 0,
  659. output_jackknife_rg = FALSE,
  660. fixed_genes = NULL)
  661. output = list(plof_h2 = bhr_plof$mixed_model$heritabilities[1,5],
  662. plof_h2_se = bhr_plof$mixed_model$heritabilities[2,5],
  663. missense_h2 = bhr_missense$mixed_model$heritabilities[1,5],
  664. missense_h2_se = bhr_missense$mixed_model$heritabilities[2,5],
  665. rg = bhr_rg_trait$rg$rg_mixed,
  666. rg_se = bhr_rg_trait$rg$rg_mixed_se)
  667. print(output)
  668. return(output)
  669. }
  671. plof_missense_compare_bhr <- sapply(traits, get_plof_missense_stats)
  672. names = rownames(plof_missense_compare_bhr)
  673. plof_missense_compare_bhr_df = as.data.frame(t(plof_missense_compare_bhr))
  674. plof_missense_compare_bhr_df = as.data.frame(sapply(1:ncol(plof_missense_compare_bhr_df), function(x) as.numeric(unlist(plof_missense_compare_bhr_df[,x]))))
  675. names(plof_missense_compare_bhr_df) <- names
  676. plof_missense_compare_bhr_df$display_name = phenotype_file$display_name[phenotype_file$phenotype_core ==1]
  678. write.csv(plof_missense_compare_bhr_df,"outputs/plof_missense_compare.csv", quote = FALSE, row.names = FALSE)
  681. ######################################## SCZ and BP #############################
  683. #BipEx
  684. #Read in the baseline model file
  685. library(pracma)
  686. library(tidyverse)
  687. library(readr)
  688. baseline_model <- read.table("~/Documents/oconnor_rotation/rarevariantproject/final_manuscript_repo/bhr/reference_files/ms_baseline_oe5.txt")
  690. #Read in the publicly available bipex variant table.
  691. bp_variantlevel_bipex <- bigreadr::fread2("~/Downloads/BipEx_variant_results.tsv")
  694. #Subset variants
  696. variant_filter = bp_variantlevel_bipex$group == "Bipolar Disorder" & #Filter to Bipolar Disorder counts
  697. bp_variantlevel_bipex$in_analysis == TRUE & #Use variant filter from Palmer et al, 2022
  698. !is.na(bp_variantlevel_bipex$gene_id) & #Remove variants with NA gene ID
  699. str_detect(bp_variantlevel_bipex$locus, "^chr\\d") #Subset to autosomal variants
  701. bp_variantlevel_bipex <- bp_variantlevel_bipex[variant_filter,
  702. c("gene_id",
  703. "consequence",
  704. "ac_case",
  705. "ac_ctrl",
  706. "locus",
  707. "mpc")]
  708. #Function to wrangle into BHR sumstats format
  709. wrangle_sumstats <- function(table,n_cases,n_controls, var_filter) {
  711. #Filter to variants of interest
  712. table = table[var_filter,]
  714. #Compute sample prevalence, will be used to compute per-sd beta
  715. prevalence = n_cases/(n_cases+n_controls)
  717. #Compute variant MAF in cases, will be used to compute per-sd beta
  718. table$AF_case = table$ac_case/(2*n_cases)
  720. #Compute variant MAFoverall, will be used to compute per-sd beta
  721. table$AF = (table$ac_case + table$ac_ctrl)/(2*(n_cases + n_controls))
  723. #calculate per-sd betas
  724. table$beta = (2 * (table$AF_case - table$AF) * prevalence)/sqrt(2 * table$AF * (1 - table$AF) * prevalence * (1 - prevalence))
  726. #calculate variant variances
  727. table$twopq = 2*table$AF * (1 - table$AF)
  729. #convert betas from per-sd (i.e. sqrt(variance explained)) to per-allele (i.e. in units of phenotype) betas.
  730. #per-allele are the usual betas reported by an exome wide association study.
  731. table$beta_perallele = table$beta/sqrt(table$twopq)
  733. #aggregate into gene-level table.
  734. sumstats = data.frame(gene = unique(table$gene_id))
  736. #each element of variant_variances is a list of variances for all variants in a gene.
  737. sumstats$variant_variances = lapply(sumstats$gene,
  738. function(x) table$twopq[table$gene_id == x])
  740. #each element of variant_variances is a list of per-allele betas for all variants in a gene.
  741. sumstats$betas = lapply(sumstats$gene,
  742. function(x) table$beta_perallele[table$gene_id == x])
  745. names(sumstats) <- c("gene", "variant_variances","betas")
  747. #add chromosome and position information
  748. #position doesn't need to be super precise, as it is only used to order genes for jackknife
  749. sumstats$gene_position <- parse_number(sapply(strsplit(table$locus[match(sumstats$gene,table$gene_id)], split = ":"), function(x) x[[2]]))
  750. sumstats$chromosome = parse_number(sapply(strsplit(table$locus[match(sumstats$gene,table$gene_id)], split = ":"), function(x) x[[1]]))
  752. #N = sum of case and control counts
  753. sumstats$N = n_cases + n_controls
  755. #we have found that in these smaller sample analyses, there are some genes with
  756. #large burden scores that are clearly outliers
  757. #we remove genes with burden scores more than 8 sd from the mean as a conservative filter.
  758. burdenscores = sapply(sumstats$variant_variances, function(x) sum(x))
  759. sumstats <- sumstats[abs(scale(burdenscores)) < 8,]
  762. return(sumstats)
  763. }
  765. bp_sumstats_ptv = wrangle_sumstats(bp_variantlevel_bipex,
  766. 14210, #N from Palmer et al, 2022
  767. 14422, #N from Palmer et al, 2022
  768. bp_variantlevel_bipex$consequence == "ptv")
  770. bp_sumstats_missenseMPC2 = wrangle_sumstats(bp_variantlevel_bipex[!is.na(bp_variantlevel_bipex$mpc),],
  771. 14210,
  772. 14422,
  773. bp_variantlevel_bipex$consequence[!is.na(bp_variantlevel_bipex$mpc)] %in% c("damaging_missense", "other_missense") &
  774. bp_variantlevel_bipex$mpc[!is.na(bp_variantlevel_bipex$mpc)] > 2)
  776. bp_sumstats_synonymous = wrangle_sumstats(bp_variantlevel_bipex,
  777. 14210,
  778. 14422,
  779. bp_variantlevel_bipex$consequence == "synonymous")
  781. #Run BHR
  782. bp_ptv_bhr <- bhr::BHR(bp_sumstats_ptv,
  783. annotations = list(baseline_model), #baseline model including constraint annotations
  784. num_blocks = 100, #number of blocks for jackknife
  785. num_null_conditions = 5, #5*num_genes null moment conditions
  786. mode = "univariate") #run in univariate mode to compute burden h2
  788. bp_missense_bhr <- bhr::BHR(bp_sumstats_missenseMPC2,
  789. annotations = list(baseline_model),
  790. num_blocks = 100,
  791. num_null_conditions = 5,
  792. mode = "univariate")
  794. bp_synonymous_bhr <- bhr::BHR(bp_sumstats_synonymous,
  795. annotations = list(baseline_model),
  796. num_blocks = 100,
  797. genomewide_correction = FALSE,
  798. num_null_conditions = 5,
  799. mode = "univariate")
  801. #convert observed scale to liability scale h2
  802. #calculate sample prevalence
  803. bp_sample_prevalence = 14210/(14210 + 14422)
  805. #this population prevalence estimate is from Ferrari et al, 2016 Bipolar Disorders
  806. population_prevalence_bp = 0.007
  809. obs2lia_factor <- function(K, P){
  810. X <- qnorm(K,lower.tail=FALSE)
  811. z <- (1/sqrt(2*pi))*(exp(-(X**2)/2))
  812. factor <- (K*(1-K)*K*(1-K))/(P*(1-P)*(z**2))
  813. return(factor)
  814. }
  816. bp_scalingfactor = obs2lia_factor(population_prevalence_bp,bp_sample_prevalence)
  818. #Burden heritability of schizophrenia
  819. #Read in table from SCHEMA website
  820. SCHEMA_variants <- bigreadr::fread2("~/Downloads/SCHEMA_variant_results.tsv",
  821. select = c("gene_id",
  822. "consequence",
  823. "ac_case",
  824. "ac_ctrl",
  825. "group",
  826. "locus",
  827. "in_analysis"))
  829. #Filter variants
  830. SCHEMA_variant_filter = (SCHEMA_variants$in_analysis | SCHEMA_variants$consequence == "synonymous_variant") & #Either filtered SCHEMA damaging variant, or synonymous
  831. SCHEMA_variants$group %in% c("EUR (exomes)","EUR (gnomAD exomes)","EUR-N (exomes)") & #In one of the primary EUR cohorts
  832. !is.na(SCHEMA_variants$gene_id) & #non-NA gene ID
  833. str_detect(SCHEMA_variants$locus, "X", negate = TRUE) & #autosomal
  834. str_detect(SCHEMA_variants$locus, "Y", negate = TRUE) #autosomal
  836. SCHEMA_variants = SCHEMA_variants[SCHEMA_variant_filter,]
  837. SCHEMA_variants = SCHEMA_variants[!is.na(SCHEMA_variants$ac_case),]
  839. #Aggregate variants into nextera and non-nextera variant tables
  840. EUR_Nextera <- SCHEMA_variants[SCHEMA_variants$group %in% c("EUR (exomes)","EUR (gnomAD exomes)"),]
  841. EUR_agg_Nextera <- aggregate(cbind(EUR_Nextera$ac_case,EUR_Nextera$ac_ctrl), by = list(EUR_Nextera$locus),FUN = sum)
  842. EUR_agg_Nextera <- EUR_agg_Nextera[EUR_agg_Nextera$V1 + EUR_agg_Nextera$V2 <= 5,]
  843. EUR_agg_Nextera$gene_id = EUR_Nextera$gene_id[match(EUR_agg_Nextera$Group.1,EUR_Nextera$locus)]
  844. EUR_agg_Nextera$consequence = EUR_Nextera$consequence[match(EUR_agg_Nextera$Group.1,EUR_Nextera$locus)]
  845. names(EUR_agg_Nextera) <- c("locus","ac_case","ac_ctrl","gene_id","consequence")
  847. EUR_NonNextera <- SCHEMA_variants[SCHEMA_variants$group %in% c("EUR-N (exomes)"),]
  848. EUR_agg_NonNextera <- aggregate(cbind(EUR_NonNextera$ac_case,EUR_NonNextera$ac_ctrl), by = list(EUR_NonNextera$locus),FUN = sum)
  849. EUR_agg_NonNextera <- EUR_agg_NonNextera[EUR_agg_NonNextera$V1 + EUR_agg_NonNextera$V2 <= 5,]
  850. EUR_agg_NonNextera$gene_id = EUR_NonNextera$gene_id[match(EUR_agg_NonNextera$Group.1,EUR_NonNextera$locus)]
  851. EUR_agg_NonNextera$consequence = EUR_NonNextera$consequence[match(EUR_agg_NonNextera$Group.1,EUR_NonNextera$locus)]
  852. names(EUR_agg_NonNextera) <- c("locus","ac_case","ac_ctrl","gene_id","consequence")
  854. #PTV wrangling
  855. scz_ptv_sumstats_Nextera <- wrangle_sumstats(EUR_agg_Nextera,
  856. 8874,
  857. 19074+23561,
  858. EUR_agg_Nextera$consequence %in% c("stop_gained",
  859. "frameshift_variant",
  860. "splice_acceptor_variant",
  861. "splice_donor_variant"))
  862. scz_ptv_sumstats_NonNextera <- wrangle_sumstats(EUR_agg_NonNextera,
  863. 7277,
  864. 11187,
  865. EUR_agg_NonNextera$consequence %in% c("stop_gained",
  866. "frameshift_variant",
  867. "splice_acceptor_variant",
  868. "splice_donor_variant"))
  870. #Missense wrangling
  871. scz_missenseMPC2_sumstats_Nextera <- wrangle_sumstats(EUR_agg_Nextera,
  872. 8874,
  873. 19074+23561,
  874. EUR_agg_Nextera$consequence %in% c("missense_variant_mpc_2-3",
  875. "missense_variant_mpc_>=3"))
  876. scz_missenseMPC2_sumstats_NonNextera <- wrangle_sumstats(EUR_agg_NonNextera,
  877. 7277,
  878. 11187,
  879. EUR_agg_NonNextera$consequence %in% c("missense_variant_mpc_2-3",
  880. "missense_variant_mpc_>=3"))
  882. #Synonymous wrangling
  883. scz_synonymous_sumstats_Nextera <- wrangle_sumstats(EUR_agg_Nextera,
  884. 8874,
  885. 19074+23561,
  886. EUR_agg_Nextera$consequence %in% c("synonymous_variant"))
  887. scz_synonymous_sumstats_NonNextera <- wrangle_sumstats(EUR_agg_NonNextera,
  888. 7277,
  889. 11187,
  890. EUR_agg_NonNextera$consequence %in% c("synonymous_variant"))
  894. #PTV BHR Models
  895. scz_ptv_bhr_Nextera <- bhr::BHR(scz_ptv_sumstats_Nextera,
  896. annotations = list(baseline_model),
  897. num_blocks = 100,
  898. mode = "univariate",
  899. output_jackknife_h2 = TRUE,
  900. all_models = TRUE)
  902. scz_ptv_bhr_NonNextera <- bhr::BHR(scz_ptv_sumstats_NonNextera,
  903. annotations = list(baseline_model),
  904. num_blocks = 100,
  905. mode = "univariate",
  906. output_jackknife_h2 = TRUE,
  907. all_models = TRUE)
  909. #Missense BHR Models
  910. scz_missenseMPC2_bhr_Nextera <- bhr::BHR(scz_missenseMPC2_sumstats_Nextera,
  911. annotations = list(baseline_model),
  912. num_blocks = 100,
  913. mode = "univariate",
  914. output_jackknife_h2 = TRUE,
  915. all_models = TRUE)
  917. scz_missenseMPC2_bhr_NonNextera <- bhr::BHR(scz_missenseMPC2_sumstats_NonNextera,
  918. annotations = list(baseline_model),
  919. num_blocks = 100,
  920. mode = "univariate",
  921. output_jackknife_h2 = TRUE,
  922. all_models = TRUE)
  924. #Synonymous BHR Models
  925. scz_synonymous_bhr_Nextera <- bhr::BHR(scz_synonymous_sumstats_Nextera,
  926. annotations = list(baseline_model),
  927. num_blocks = 100,
  928. mode = "univariate",
  929. output_jackknife_h2 = TRUE,
  930. all_models = TRUE)
  932. scz_synonymous_bhr_NonNextera <- bhr::BHR(scz_synonymous_sumstats_NonNextera,
  933. annotations = list(baseline_model),
  934. num_blocks = 100,
  935. mode = "univariate",
  936. output_jackknife_h2 = TRUE,
  937. all_models = TRUE)
  939. #Functions for meta-analysis
  940. get_metabeta <- function(betas, ses){
  941. we = 1 / (ses)^2
  942. return(sum(betas * we) / sum(we))
  943. }
  944. get_metase <- function(ses){
  945. we = 1 / (ses)^2
  946. return(sqrt(1/sum(we)))
  947. }
  949. #Function for meta-analyzing nextera and non-nextera samples, and extracting key parameters of interest
  951. SCZ_meta_analysis <- function(modelNextera,modelNonNextera,sumstatsNextera,sumstatsNonNextera, liability = TRUE){
  954. #calculate factor for conversion to liability scale.
  955. #Population prevalence is from Charlson et al, 2016, Schizophrenia Bulletin
  956. nextera_prevalence = 8874/(19074+23561+8874)
  957. nonnextera_prevalence = 7277/(11187+7277)
  958. population_prevalence = 0.0028
  960. obs2lia_factor <- function(K, P){
  961. X <- qnorm(K,lower.tail=FALSE)
  962. z <- (1/sqrt(2*pi))*(exp(-(X**2)/2))
  963. factor <- (K*(1-K)*K*(1-K))/(P*(1-P)*(z**2))
  964. return(factor)
  965. }
  967. nextera_scalingfactor = obs2lia_factor(population_prevalence,nextera_prevalence)
  968. nonnextera_scalingfactor = obs2lia_factor(population_prevalence,nonnextera_prevalence)
  970. #If observed scale desired, set scaling factor to 1
  971. if (liability == FALSE){
  972. nextera_scalingfactor = 1
  973. nonnextera_scalingfactor = 1
  975. }
  977. #extract heritability from constrained genes in nextera and non-nextera samples
  978. h2_constrained = c(modelNextera$mixed_model$heritabilities[1,1]*nextera_scalingfactor,
  979. modelNonNextera$mixed_model$heritabilities[1,1]*nonnextera_scalingfactor)
  980. h2_constrained_se = c(modelNextera$mixed_model$heritabilities[2,1]*nextera_scalingfactor,
  981. modelNonNextera$mixed_model$heritabilities[2,1]*nonnextera_scalingfactor)
  983. #extract overall heritability in nextera and non-nextera samples
  984. h2_all = c(modelNextera$mixed_model$heritabilities[1,5]*nextera_scalingfactor,
  985. modelNonNextera$mixed_model$heritabilities[1,5]*nonnextera_scalingfactor)
  986. h2_all_se = c(modelNextera$mixed_model$heritabilities[2,5]*nextera_scalingfactor,
  987. modelNonNextera$mixed_model$heritabilities[2,5]*nonnextera_scalingfactor)
  989. #meta-analyze constrained heritability ("numerator" of fraction explained by constrained genes)
  990. numerator = get_metabeta(h2_constrained,h2_constrained_se)
  991. numerator_se = get_metase(h2_constrained_se)
  993. #meta-analyze total heritability ("denominator" of fraction explained by constrained genes)
  994. denominator = get_metabeta(h2_all,h2_all_se)
  995. denominator_se = get_metase(h2_all_se)
  997. #get point estimate of fraction explained by constrained genes
  998. fraction = numerator/denominator
  1000. #compute fraction of alleles in constrained annotation
  1001. sumstatsNextera_annot = merge(sumstatsNextera, baseline_model, by.x = "gene", by.y = "gene")
  1002. sumstatsNonNextera_annot = merge(sumstatsNonNextera, baseline_model, by.x = "gene", by.y = "gene")
  1003. sumstatsNextera_annot$burden_score = sapply(sumstatsNextera_annot$variant_variances, function(x) sum(x))
  1004. sumstatsNonNextera_annot$burden_score = sapply(sumstatsNonNextera_annot$variant_variances, function(x) sum(x))
  1006. total_variants = sum(sumstatsNextera_annot$burden_score) + sum(sumstatsNonNextera_annot$burden_score)
  1007. constrained_variants = sum(sumstatsNextera_annot$burden_score * sumstatsNextera_annot$baseline_oe1_total5) + sum(sumstatsNonNextera_annot$burden_score*sumstatsNonNextera_annot$baseline_oe1_total5)
  1008. fraction_constrained_variants = constrained_variants/total_variants
  1010. #compute constraint enrichment point estimate
  1011. enrichment = fraction/fraction_constrained_variants
  1013. #get SE of fraction with delta method
  1015. #First, compute covariance matrix of (h2 total, h2 constrained), for nextera and non-nextera
  1016. sigma = matrix(data = NA, nrow = 2,ncol = 2)
  1017. sigma[1,1] = numerator_se^2
  1018. sigma[2,2] = denominator_se^2
  1020. #jackknife variances of h2 constrained and h2 total, and jackknfie covariance between h2 constrained and h2 total
  1021. Nextera_constrained_jackknife = modelNextera$subthreshold_genes$jackknife_h2[1,]*nextera_scalingfactor
  1022. Nextera_total_jackknife = modelNextera$subthreshold_genes$jackknife_h2[5,]*nextera_scalingfactor
  1024. variance_constrained_Nextera = ((length(Nextera_constrained_jackknife) -1)/length(Nextera_constrained_jackknife))*sum((Nextera_constrained_jackknife - mean(Nextera_constrained_jackknife))^2)
  1025. variance_total_Nextera = ((length(Nextera_total_jackknife) -1)/length(Nextera_total_jackknife))*sum((Nextera_total_jackknife - mean(Nextera_total_jackknife))^2)
  1026. covariance_Nextera_jackknife = ((length(Nextera_constrained_jackknife) -1)/length(Nextera_constrained_jackknife)) * sum((Nextera_constrained_jackknife - mean(Nextera_constrained_jackknife))*(Nextera_total_jackknife - mean(Nextera_total_jackknife)))
  1028. sigma_nextera = matrix(data = NA, nrow = 2,ncol = 2)
  1029. sigma_nextera[1,1] = variance_constrained_Nextera
  1030. sigma_nextera[2,2] = variance_total_Nextera
  1031. sigma_nextera[1,2] = covariance_Nextera_jackknife
  1032. sigma_nextera[2,1] = covariance_Nextera_jackknife
  1034. NonNextera_constrained_jackknife = modelNonNextera$subthreshold_genes$jackknife_h2[1,]*nonnextera_scalingfactor
  1035. NonNextera_total_jackknife = modelNonNextera$subthreshold_genes$jackknife_h2[5,]*nonnextera_scalingfactor
  1036. variance_constrained_NonNextera = ((length(NonNextera_constrained_jackknife) -1)/length(NonNextera_constrained_jackknife))*sum((NonNextera_constrained_jackknife - mean(NonNextera_constrained_jackknife))^2)
  1037. variance_total_NonNextera = ((length(NonNextera_total_jackknife) -1)/length(NonNextera_total_jackknife))*sum((NonNextera_total_jackknife - mean(NonNextera_total_jackknife))^2)
  1038. covariance_NonNextera_jackknife = ((length(NonNextera_constrained_jackknife) -1)/length(NonNextera_constrained_jackknife)) * sum((NonNextera_constrained_jackknife - mean(NonNextera_constrained_jackknife))*(NonNextera_total_jackknife - mean(NonNextera_total_jackknife)))
  1040. sigma_NonNextera = matrix(data = NA, nrow = 2,ncol = 2)
  1041. sigma_NonNextera[1,1] = variance_constrained_NonNextera
  1042. sigma_NonNextera[2,2] = variance_total_NonNextera
  1043. sigma_NonNextera[1,2] = covariance_NonNextera_jackknife
  1044. sigma_NonNextera[2,1] = covariance_NonNextera_jackknife
  1046. #meta-analyze covariance matrices
  1047. sigma_meta = pracma::inv(pracma::inv(sigma_nextera) + pracma::inv(sigma_NonNextera))
  1050. #get gradient
  1051. dg_dh2c = 1/denominator
  1052. dg_dh2all = (-numerator)/(denominator^2)
  1054. fraction_gradient = matrix(c(dg_dh2c,dg_dh2all), ncol = 1)
  1056. #compute variance of fraction of heritability explained by constrained genes
  1057. fraction_se = sqrt(t(fraction_gradient) %*% sigma_meta %*% fraction_gradient)
  1059. #convert fraction to enrichment
  1060. enrichment_se = fraction_se/fraction_constrained_variants
  1062. return(list(bhr_h2 = denominator,
  1063. bhr_h2_se = denominator_se,
  1064. fraction_constrained = fraction,
  1065. fraction_se = fraction_se,
  1066. enrichment_constrained = enrichment,
  1067. enrichment_constrained_se = enrichment_se))
  1069. }
  1071. scz_ptv_bhr_output = SCZ_meta_analysis(scz_ptv_bhr_Nextera,
  1072. scz_ptv_bhr_NonNextera,
  1073. scz_ptv_sumstats_Nextera,
  1074. scz_ptv_sumstats_NonNextera)
  1076. scz_missense_bhr_output = SCZ_meta_analysis(scz_missenseMPC2_bhr_Nextera,
  1077. scz_missenseMPC2_bhr_NonNextera,
  1078. scz_missenseMPC2_sumstats_Nextera,
  1079. scz_missenseMPC2_sumstats_NonNextera)
  1081. scz_synonymous_bhr_output = SCZ_meta_analysis(scz_synonymous_bhr_Nextera,
  1082. scz_synonymous_bhr_NonNextera,
  1083. scz_synonymous_sumstats_Nextera,
  1084. scz_synonymous_sumstats_NonNextera)
  1086. #Gather output
  1087. output_df_SCZBP <- data.frame(dx =c(rep("SCZ",3),rep("BP",3)),
  1088. class = c(c("pLoF","Missense (MPC >2)", "Syn"),c("pLoF","Missense (MPC >2)", "Syn")),
  1089. bhr_h2 = c(scz_ptv_bhr_output$bhr_h2,
  1090. scz_missense_bhr_output$bhr_h2,
  1091. scz_synonymous_bhr_output$bhr_h2,
  1092. bp_ptv_bhr$mixed_model$heritabilities[1,5]*bp_scalingfactor,
  1093. bp_missense_bhr$mixed_model$heritabilities[1,5]*bp_scalingfactor,
  1094. bp_synonymous_bhr$mixed_model$heritabilities[1,5]*bp_scalingfactor),
  1095. bhr_h2_se = c(scz_ptv_bhr_output$bhr_h2_se,
  1096. scz_missense_bhr_output$bhr_h2_se,
  1097. scz_synonymous_bhr_output$bhr_h2_se,
  1098. bp_ptv_bhr$mixed_model$heritabilities[2,5]*bp_scalingfactor,
  1099. bp_missense_bhr$mixed_model$heritabilities[2,5]*bp_scalingfactor,
  1100. bp_synonymous_bhr$mixed_model$heritabilities[2,5]*bp_scalingfactor),
  1101. fraction_constrained = c(scz_ptv_bhr_output$fraction_constrained,
  1102. scz_missense_bhr_output$fraction_constrained,
  1103. scz_synonymous_bhr_output$fraction_constrained,
  1104. bp_ptv_bhr$mixed_model$fractions[1,1],
  1105. bp_missense_bhr$mixed_model$fractions[1,1],
  1106. bp_synonymous_bhr$mixed_model$fractions[1,1]),
  1107. fraction_constrained_se = c(scz_ptv_bhr_output$fraction_constrained_se,
  1108. scz_missense_bhr_output$fraction_constrained_se,
  1109. scz_synonymous_bhr_output$fraction_constrained_se,
  1110. bp_ptv_bhr$mixed_model$fractions[2,1],
  1111. bp_missense_bhr$mixed_model$fractions[2,1],
  1112. bp_synonymous_bhr$mixed_model$fractions[2,1]),
  1113. enrichment_constrained = c(scz_ptv_bhr_output$enrichment_constrained,
  1114. scz_missense_bhr_output$enrichment_constrained,
  1115. scz_synonymous_bhr_output$enrichment_constrained,
  1116. bp_ptv_bhr$mixed_model$enrichments[1,1],
  1117. bp_missense_bhr$mixed_model$enrichments[1,1],
  1118. bp_synonymous_bhr$mixed_model$enrichments[1,1]),
  1119. enrichment_constrained_se = c(scz_ptv_bhr_output$enrichment_constrained_se,
  1120. scz_missense_bhr_output$enrichment_constrained_se,
  1121. scz_synonymous_bhr_output$enrichment_constrained_se,
  1122. bp_ptv_bhr$mixed_model$enrichments[2,1],
  1123. bp_missense_bhr$mixed_model$enrichments[2,1],
  1124. bp_synonymous_bhr$mixed_model$enrichments[2,1]))
  1126. #fractions of SCZ h2 explained by sig genes
  1127. #Load in gnomad information to match gene ids to gene names
  1128. gnomad_information <- data.frame(data.table::fread("~/Documents/oconnor_rotation/rarevariantproject/final_manuscript_repo/rv_h2/reference_files/gnomad.v2.1.1.lof_metrics.by_gene.txt",
  1129. select = c("gene","gene_id", "chromosome", "start_position", "end_position", "pLI")))
  1130. gnomad_information = gnomad_information[gnomad_information$chromosome != "X",]
  1131. gnomad_information = gnomad_information[gnomad_information$chromosome != "Y",]
  1132. gnomad_information_counts = data.frame(table(gnomad_information$gene))
  1133. gnomad_information_counts = gnomad_information_counts[gnomad_information_counts$Freq == 1,]
  1134. gene_information <- gnomad_information[gnomad_information$gene %in% gnomad_information_counts$Var1,]
  1135. gene_information$midpoint <- (gene_information$start_position + gene_information$end_position)/2
  1137. #9 autosomal SCHEMA genes
  1138. SCHEMAgenes = c("SETD1A", "CUL1", "XPO7","TRIO","CACNA1G","SP4","GRIA3","GRIN2A","HERC1","RB1CC1")
  1139. SCHEMAgenes_id = gene_information$gene_id[match(SCHEMAgenes,gene_information$gene)]
  1140. SCHEMAgenes_id = SCHEMAgenes_id[!is.na(SCHEMAgenes_id)]
  1142. #Run BHR with SCHEMA significant genes as fixed effects.
  1143. nextera_model_SCHEMA = bhr::BHR(scz_ptv_sumstats_Nextera,
  1144. annotations = list(baseline_model),
  1145. fixed_genes =SCHEMAgenes_id,
  1146. num_blocks = 100,
  1147. mode = "univariate")
  1149. nonnextera_model_SCHEMA = bhr::BHR(scz_ptv_sumstats_NonNextera,
  1150. annotations = list(baseline_model),
  1151. fixed_genes =SCHEMAgenes_id,
  1152. num_blocks = 100,
  1153. mode = "univariate")
  1155. #Meta analyze fraction of heritability explained by sig genes
  1156. meta_frac_SCHEMA = get_metabeta(c(nextera_model_SCHEMA$significant_genes$fraction_burdenh2_significant,
  1157. nonnextera_model_SCHEMA$significant_genes$fraction_burdenh2_significant),
  1158. c(nextera_model_SCHEMA$significant_genes$fraction_burdenh2_significant_se,
  1159. nonnextera_model_SCHEMA$significant_genes$fraction_burdenh2_significant_se))
  1161. meta_frac_SCHEMA_se = get_metase(c(nextera_model_SCHEMA$significant_genes$fraction_burdenh2_significant_se,
  1162. nonnextera_model_SCHEMA$significant_genes$fraction_burdenh2_significant_se))
  1164. #rg between scz and bp
  1165. nextera_bp_rg = bhr::BHR(mode = "bivariate",
  1166. trait1_sumstats = scz_ptv_sumstats_Nextera,
  1167. trait2_sumstats = bp_sumstats_ptv,
  1168. annotations = list(baseline_model),
  1169. num_blocks = 100)
  1171. nonnextera_bp_rg = bhr::BHR(mode = "bivariate",
  1172. trait1_sumstats = scz_ptv_sumstats_NonNextera,
  1173. trait2_sumstats = bp_sumstats_ptv,
  1174. annotations = list(baseline_model),
  1175. num_blocks = 100)
  1177. write.csv(output_df_SCZBP,"~/Documents/SCZ_BP_output.csv", quote = FALSE, row.names = FALSE)

run_BHR.R at commit 108e9c8, under MIT · at the source

Overview

Authors: Hanna Abe1, Calwing Liao2,3,4, Lide Han5, Theodore Morley5, Michael E Talkowski2,3,4, Kristen J Brennand6, Douglas M Ruderfer5,7
ORCID iDs: Hanna Abe, Lide Han, Kristen J Brennand
  1. Vanderbilt University, Vanderbilt Genetics Institute, Nashville, TN USA
  2. Analytic and Translational Genetics Unit, Department of Medicine, Massachusetts General Hospital, Boston, MA USA
  3. Stanley Center for Psychiatric Research, Broad Institute of MIT and Harvard, Cambridge, MA USA
  4. Center for Genomic Medicine, Massachusetts General Hospital, Boston, MA USA
  5. Division of Genetic Medicine, Department of Medicine, Vanderbilt University Medical Center, Nashville, TN USA
  6. Departments of Psychiatry and Genetics, Division of Molecular Psychiatry, Department of Genetics, Wu Tsai Institute, Yale University School of Medicine, New Haven, CT USA
  7. Department of Biomedical Informatics and Psychiatry and Behavioral Sciences, Vanderbilt University Medical Center, Nashville, TN USA
Institutions: Vanderbilt University (United States); Broad Institute (United States); Massachusetts General Hospital (United States); Stanley Center for Psychiatric Research; Vanderbilt University Medical Center (United States); Yale University (United States)
Journal: Molecular psychiatry, volume 31, issue 8, pages 4787-4798
Dates: received 26 August 2025; accepted 16 March 2026; published online 7 April 2026; in print 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41380-026-03571-x · PMID 41946833 · PMCID PMC13364661 · OpenAlex W4413379657
Open access: hybrid, a free copy (OpenAlex)
Preprint: doi:10.1101/2025.08.13.25333593
Status: code verified
Categories: genetics / omics (modality), human (organism), cellular / molecular (subfield)
Methods: Statistics, Machine learning, Preprocessing, Connectivity
Keywords: Genetics, Psychiatric disorders
MeSH: Mental Disorders*, Brain, Gene Expression Profiling, Genetic Predisposition to Disease, Genetic Variation, Genome-Wide Association Study, Humans, Neurodegenerative Diseases, Polymorphism, Single Nucleotide, Transcriptome (* major topic)
Topic: Genetic Associations and Epidemiology (Genetics, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: U.S. Department of Health & Human Services | NIH | National Institute of Mental Health (NIMH) (R01MH123155); U.S. Department of Health &amp; Human Services | NIH | National Institute of Mental Health (R01MH123155); NIMH NIH HHS (R01 MH123155)
Citations: not cited yet (Europe PMC); 53 references in the paper

Abstract

Genome-wide association studies (GWAS) and large-scale rare variant burden analyses have identified both common and rare loss-of-function variants associated with neuropsychiatric and neurodegenerative disorders. Yet, the shared biological processes influenced by both classes of variation remain poorly characterized. In this study, we utilized transcriptomic data from 933 post-mortem brain samples to identify genes that show convergent coexpression with GWAS and rare variant burden risk genes across six brain disorders. Despite largely distinct sets of significant risk genes from GWAS and rare variant burden studies, we found a significant overlap in their convergently coexpressed genes. These convergent genes showed enrichment for common and rare variant heritability and highlighted key biological pathways and cell-type markers impacted by both types of genetic variation. Compared to genes coexpressed with one variant class, shared convergent genes exhibited stronger evolutionary constraint and greater enrichment for known drug targets, underscoring their potential therapeutic relevance. Collectively, our results establish a systematic and generalizable framework for integrating coexpression data with genetic risk to reveal transcriptional programs supported by both common and rare variant evidence, offering mechanistic insights into neuropsychiatric diseases.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 6 matches between paragraphs and lines of code.

FinucaneLab/gene_features

License: GPL-3.0
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: 41d006d5c0e9cbdb8dfafdf0ee821f165e5e708c, 19 May 2022
Languages: R (80), Shell (1)
Size: 1,381 files, 81 scripts
Software Heritage: not archived
Found in: the text, “Gene prioritization from GWAS and rare variant b”
Holds: README, license file, environment (code/install.R)
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: tidyverse (80 files), data.table (79 files), Seurat (79 files), reticulate (77 files), Harmony (2 files), Monocle 3 (1 file), reshape2 (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
83 files
At the source: github.com/FinucaneLab/gene_features

ajaynadig/bhr

License: MIT
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: 108e9c897d1b7da2a6b0bace9f63e54da55e25a0, 19 June 2026
Languages: MATLAB (13), R (10), Python (1)
Size: 37 files, 24 scripts
Software Heritage: archived
Found in: “Code availability”
Holds: README, license file, environment (DESCRIPTION), tests, 1 notebook
Not found: CITATION.cff, continuous integration, documentation
Tools: tidyverse (4 files), data.table (2 files), Statistics and Machine Learning Toolbox (2 files), ggplot2 (1 file), patchwork (1 file), reshape2 (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
26 files
At the source: github.com/ajaynadig/bhr

FinucaneLab/pops

License: GPL-3.0
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: 76eb86cba10254490003c8f4dc7ff5ce492d3667, 21 July 2025
Languages: Python (2)
Size: 25 files, 2 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, environment (requirements.txt)
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (2 files), pandas (2 files), scikit-learn (1 file), SciPy (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
4 files
At the source: github.com/FinucaneLab/pops

hakyimlab/MetaXcan

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: e069063a10539fe92adadb645fc667a40c2cf885, 8 September 2026
Languages: Python (118), R (2), Shell (2)
Size: 318 files, 122 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, environment (software/conda_env.yaml, software/requirements.txt, software/setup.cfg, software/setup.py), tests, documentation, 1 notebook
Not found: CITATION.cff, continuous integration
Tools: NumPy (52 files), pandas (40 files), SciPy (7 files), h5py (4 files), statsmodels (4 files)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
124 files
At the source: github.com/hakyimlab/MetaXcan

RuderferLab/coexpression_convergence

License: none: the authors keep all their rights
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: 8e847ffc072415891ecd6f5178c3f982836c65b5, 21 July 2022
Languages: R (2)
Size: 5 files, 2 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: data.table (1 file), tidyverse (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
3 files
At the source: github.com/RuderferLab/coexpression_convergence

bulik/ldsc

License: GPL-3.0
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Commit: 2fdeeb3b44379408794154993dbd6101b8946b7e, 16 January 2026
Languages: Python (19), MATLAB (4), Perl (1), R (1)
Size: 1,093 files, 25 scripts
Software Heritage: archived
Found in: “Code availability”
Holds: README, license file, environment (environment.yml, requirements.txt, setup.py), tests
Not found: CITATION.cff, continuous integration, documentation
Tools: NumPy (18 files), pandas (11 files), SciPy (5 files), BEDTools (2 files)
Availability: 1 check, the latest on 26 September 2026: the link answers
  • 26 September 2026: the link answers
27 files
At the source: github.com/bulik/ldsc

YuLab-SMU/clusterProfiler

License: none: the authors keep all their rights
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: d10e74853722ca5f3fb3a0c3466c4649e0fedada, 26 September 2026
Languages: R (59), Quarto (1)
Size: 163 files, 60 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, environment (DESCRIPTION), tests, continuous integration, documentation, 3 notebooks
Not found: license file, CITATION.cff
Tools: clusterProfiler (9 files), tidyverse (8 files), igraph (3 files), ggplot2 (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
61 files
At the source: github.com/YuLab-SMU/clusterProfiler

Code availability

For gene prioritization: MAGMA (https://cncr.nl/research/magma), Polygenic priority score (PoPs, https://github.com/FinucaneLab/pops), TWAS-SPrediXcan (https://github.com/hakyimlab/MetaXcan) For convergent coexpression analysis (https://github.com/RuderferLab/coexpression_convergence). For SNP and rare variant heritability estimation: SLDSC (https://github.com/bulik/ldsc), BHR (https://github.com/ajaynadig/bhr). For gene ontology and other functional enrichment analysis: enrichr (https://maayanlab.cloud/Enrichr), ClusterProfiler (https://github.com/YuLab-SMU/clusterProfiler), rrvgo (10.18129/B9.bioc.rrvgo), and SMR (https://yanglab.westlake.edu.cn/software/smr).

Reproduced under the paper's license (CC BY), from the paper cited above.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 7 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 316 scripts, each with its path and the digest of its content;
  • 6 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data Availability Statement

GWAS summary statistics data from Psychiatric genomics consortium (PGC) : (https://pgc.unc.edu/for-researchers/download-results). Exome burden summary: (https://schema.broadinstitute.org). Postmortem transcriptome data can be accessed from: CMC (10.7303/syn22000731.1), GTEx brain tissues (https://www.gtexportal.org/). The Open Targets Platform (https://platform.opentargets.org). SMR database (https://cnsgenomics.com/software/smr/#DataResource).

For gene prioritization: MAGMA (https://cncr.nl/research/magma), Polygenic priority score (PoPs, https://github.com/FinucaneLab/pops), TWAS-SPrediXcan (https://github.com/hakyimlab/MetaXcan) For convergent coexpression analysis (https://github.com/RuderferLab/coexpression_convergence). For SNP and rare variant heritability estimation: SLDSC (https://github.com/bulik/ldsc), BHR (https://github.com/ajaynadig/bhr). For gene ontology and other functional enrichment analysis: enrichr (https://maayanlab.cloud/Enrichr), ClusterProfiler (https://github.com/YuLab-SMU/clusterProfiler), rrvgo (10.18129/B9.bioc.rrvgo), and SMR (https://yanglab.westlake.edu.cn/software/smr).

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 29 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 7 authors, 2 keywords, 10 MeSH terms, 3 funders, 53 references.

Cite

This paper

Abe, H., Liao, C., Han, L., Morley, T., Talkowski, M. E., Brennand, K. J., & Ruderfer, D. M. (2026). Convergent coexpression reveals shared biological mechanisms underlying common and rare variant risk in six neuropsychiatric disorders. Molecular psychiatry, 31(8), 4787-4798. https://doi.org/10.1038/s41380-026-03571-x

BibTeX

@article{abe2026convergent,
author = {Abe, Hanna and Liao, Calwing and Han, Lide and Morley, Theodore and Talkowski, Michael E and Brennand, Kristen J and Ruderfer, Douglas M},
title = {{Convergent coexpression reveals shared biological mechanisms underlying common and rare variant risk in six neuropsychiatric disorders}},
journal = {Molecular psychiatry},
year = {2026},
month = apr,
volume = {31},
number = {8},
pages = {4787--4798},
publisher = {Springer Nature},
issn = {1359-4184},
doi = {10.1038/s41380-026-03571-x},
url = {https://doi.org/10.1038/s41380-026-03571-x},
pmid = {41946833},
pmcid = {PMC13364661}
}

RIS

TY - JOUR
AU - Abe, Hanna
AU - Liao, Calwing
AU - Han, Lide
AU - Morley, Theodore
AU - Talkowski, Michael E
AU - Brennand, Kristen J
AU - Ruderfer, Douglas M
TI - Convergent coexpression reveals shared biological mechanisms underlying common and rare variant risk in six neuropsychiatric disorders
T2 - Molecular psychiatry
J2 - Mol Psychiatry
PY - 2026
DA - 2026/04/07
VL - 31
IS - 8
SP - 4787
EP - 4798
SN - 1359-4184
PB - Springer Nature
DO - 10.1038/s41380-026-03571-x
UR - https://doi.org/10.1038/s41380-026-03571-x
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41380-026-03571-x",
"type": "article-journal",
"title": "Convergent coexpression reveals shared biological mechanisms underlying common and rare variant risk in six neuropsychiatric disorders",
"container-title": "Molecular psychiatry",
"author": [
{
"family": "Abe",
"given": "Hanna"
},
{
"family": "Liao",
"given": "Calwing"
},
{
"family": "Han",
"given": "Lide"
},
{
"family": "Morley",
"given": "Theodore"
},
{
"family": "Talkowski",
"given": "Michael E"
},
{
"family": "Brennand",
"given": "Kristen J"
},
{
"family": "Ruderfer",
"given": "Douglas M"
}
],
"container-title-short": "Mol Psychiatry",
"volume": "31",
"issue": "8",
"page": "4787-4798",
"DOI": "10.1038/s41380-026-03571-x",
"PMID": "41946833",
"PMCID": "PMC13364661",
"ISSN": "1359-4184",
"publisher": "Springer Nature",
"URL": "https://doi.org/10.1038/s41380-026-03571-x",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
7
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41467-026-75193-4 [code]
Multi-ancestry gene expression models amplify transcriptome-wide association study discovery and validation.
Journal: Nature communications
In common: Monocle 3, igraph, Seurat, 10 other tools, genetics / omics, cellular / molecular, 5 references
[2] doi:10.1038/s41593-026-02247-7 [code]
Transcriptomic and phenotypic convergence of neurodevelopmental disorder risk genes in vitro and in vivo.
Journal: Nature neuroscience
In common: Statistics and Machine Learning Toolbox, genetics / omics, cellular / molecular, 12 references, author Kristen J. Brennand
[3] doi:10.1038/s41588-026-02646-3 [code]
Co-expression-based models improve eQTL predictions for transcriptome-wide association studies and highlight new schizophrenia-associated genes.
Journal: Nature genetics
In common: reshape2, h5py, data.table, 5 other tools, genetics / omics, cellular / molecular, 11 references
[4] doi:10.1186/s12967-026-08266-z [code]
Single-cell multi-omic integration analysis prioritizes druggable genes and reveals cell-type-specific causal effects in glioblastomagenesis.
Journal: Journal of translational medicine
In common: Monocle 3, reticulate, igraph, 7 other tools, genetics / omics, cellular / molecular, 5 references
[5] doi:10.1038/s44318-026-00818-9 [code]
FAM134B-mediated ER-phagy degrades APP and suppresses Alzheimer's disease pathology.
Journal: The EMBO journal
In common: Monocle 3, BEDTools, Harmony, 13 other tools, cellular / molecular
[6] doi:10.1038/s41467-026-71542-5 [code]
Astrocyte fatty acid metabolism as a driver of risk for major depressive disorder.
Journal: Nature communications
In common: igraph, Seurat, reshape2, 4 other tools, genetics / omics, cellular / molecular, 10 references
[7] doi:10.1038/s41467-026-76675-1 [code]
Long-read proteogenomic atlas of human neuronal differentiation reveals isoform diversity informing neurodevelopmental risk mechanisms.
Journal: Nature communications
In common: Monocle 3, BEDTools, reticulate, 12 other tools, genetics / omics, 1 reference
[8] doi:10.1016/j.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: Monocle 3, Harmony, igraph, 12 other tools, genetics / omics, cellular / molecular, 1 reference
[9] doi:10.1038/s41467-026-76676-0 [code]
Determinants of functional burden pleiotropy and gene dosage responses across human traits.
Journal: Nature communications
In common: igraph, reshape2, data.table, 8 other tools, genetics / omics, cellular / molecular, 5 references
[10] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: Monocle 3, reticulate, igraph, 13 other tools, cellular / molecular

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.