Multi-ancestry gene expression models amplify transcriptome-wide association study discovery and validation.

1. ---
2. title: "PGC_multi_ancestry_analysis"
3. author: "Xavier Bledsoe"
4. date: "4/30/2025"
5. output:
6. html_document:
7. code_folding: show
8. toc: true
9. toc_float: true
10. number_sections: true
11. ---
13. ```{r setup, include=FALSE}
14. knitr::opts_chunk$set(echo = TRUE, collapse = TRUE)
15. library(data.table)
16. library(stringr)
17. library(ggplot2)
18. library(ggrepel)
19. library(viridis)
20. library(pheatmap)
21. library(RColorBrewer)
22. library(ggpubr)
23. #library(kableExtra)
25. library(data.table)
26. library(dplyr)
27. library(ggrepel)
28. library(stringr)
29. library(colorspace)
31. library(data.table)
32. library(ggplot2)
33. library(pals)
34. library(stringr)
35. library(gridExtra)
36. library(grid)
37. library(ggstance)
38. library(viridis)
39. library(ggpointdensity)
40. library(ggh4x)
41. library(forcats)
44. ```
48. # PGC TWAS analysis
50. ## TWAS descriptive statistics
52. ### Fig 1a: Nashville Plots of TWAS Data
54. For all 6 GWAS from the PGC consortium, we want to visualize the distribution of
55. TWAS results across the GALA II/SAGE (G2S) and GTEx whole blood results.
58. ```{r fig.height=8, fig.width = 15}
59. dt_ <- fread('Data/all_exc_assocs_PGCancestryTWAS.txt', header = TRUE)
61. # exclude replication/comparison data from initial analysis.
62. gwas <- c('PGC_AA_MDD_2023', 'PGC_AA_PTSD_2024','PGC_AA_SCZ_2022','PGC_ASIAN_SCZ_2022','PGC_EUR_PTSD_2024','PGC_HNA_PTSD_2024', 'PGC_EURO_SCZ_2022', 'PGC_LATINO_SCZ_2022')
63. mod <- c('PrediXcan_Brain_Cortex')
65. dt <- dt_[!(gwas_phenotype %in% gwas | training_model %in% mod),]
68. dt[training_model == 'All_Whole_Blood', `gene expression model` := 'G2S All WB']
69. dt[training_model == 'AA_Whole_Blood', `gene expression model` := 'G2S African American WB']
70. dt[training_model == 'MX_Whole_Blood', `gene expression model` := 'G2S Mexican American WB']
71. dt[training_model == 'PR_Whole_Blood', `gene expression model` := 'G2S Puerto Rican WB']
72. dt[training_model == 'PrediXcan_Whole_Blood', `gene expression model` := 'GTEx European WB']
75. tms <- rev(c('G2S All WB', 'G2S African American WB', 'G2S Mexican American WB', 'G2S Puerto Rican WB', 'GTEx European WB'))
77. dt$`gene expression model` <- factor(dt$`gene expression model`, levels = tms)
79. adtwas2 <- dt
81. # print manhattan plots for both the brain specific associations and the raw associations across all tissues
83. chrom_sizes <- data.table(CHR = c(1:22), chr_length = c(
84. 248956422, #1
85. 242193529,
86. 198295559,
87. 190214555,
88. 181538259, #5
89. 170805979,
90. 159345973,
91. 145138636,
92. 138394717,
93. 133797422, #10
94. 135086622,
95. 133275309,
96. 114364328,
97. 114364328,
98. 101991189, #15
99. 90338345,
100. 83257441,
101. 80373285,
102. 58617616,
103. 64444167, #20
104. 46709983,
105. 50818468
106. ))
107. # Get the cumulative start locations for each chromosome in terms of base pairs.
108. chrom_sizes[, chr_start := (cumsum(chr_length)-248956422)]
110. # Set the x-axis label locations for each chromosome to be right in the middle of the chromosome block
111. chrom_sizes[, chr_label_loc := (chr_start + chr_length/2)]
116. all_genes <- fread('Data/all_ensembl.txt.gz', header = TRUE, stringsAsFactors=FALSE)
117. # This file was downloaded from biomart website on 5/25/2022. The parameters were as follows:
118. # Dataset: Human genes (GRCh38.p13)
119. # Filters: Chromosome/scaffold: 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 , 10 , 11 , 12 , 13 , 14 , 15 , 16 , 17 , 18 , 19 , 20 , 21 , 22
120. # Attributes: Gene stable ID, Gene start (bp), Chromosome/scaffold name
121. # Export all files to: TSV
122. #
123. # I then uploaded the file to my personal directory and compressed it via the gzip commandline command.
125. all_genes <- setNames(all_genes, c('ensembl_gene_id', 'start_position', 'chromosome_name'))
126. all_genes_loc <- as.data.table(merge(all_genes, chrom_sizes[, c('CHR','chr_start')], by.x = 'chromosome_name', by.y = 'CHR'))
127. all_genes_loc[, gn_start := rowSums(.SD), .SDcols = c("start_position", "chr_start")]
128. all_genes_loc <- all_genes_loc[, c('ensembl_gene_id', 'chromosome_name', 'gn_start')]
129. #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
131. dt=adtwas2
132. dt$gene <- str_remove_all(dt$gene, '\\..+')
133. # Iterate through the tissue models and NIDPs from the TWAS results
136. # Get the top 15 p-values for labeling purposes
137. top15 <- dt[, .(pvalue = min(pvalue)), by = c('gene','gene_name')][order(pvalue)][1:8, c('gene','gene_name','pvalue')]
138. top15x = setDT(merge(top15, all_genes_loc, by.x = 'gene', by.y = 'ensembl_gene_id'))
139. top15x$short_name = 'top_15'
141. top25 <- dt[, .(pvalue = min(pvalue)), by = c('gene','gene_name')][order(pvalue)][1:25, c('gene','gene_name','pvalue')]
142. top25x = setDT(merge(top25, all_genes_loc, by.x = 'gene', by.y = 'ensembl_gene_id'))
143. top25x$short_name = 'top_25'
145. twas2 <- merge(dt, all_genes_loc, by.x = 'gene', by.y = 'ensembl_gene_id')
147. twas2$`gene expression model` <- factor(twas2$`gene expression model`, levels = tms)
150. fig2 = ggplot(twas2, aes(x=gn_start, y=-log10(pvalue))) +
151. # Show all points
152. geom_point( aes(color=`gene expression model`), alpha=0.8, size=1.5) +
153. geom_hline(yintercept=-log10(0.05/(dim(adtwas2)[1])), color = "red", size=0.5) +
154. facet_wrap(vars(gwas_phenotype), scales = 'free', nrow = 3) +
155. xlab('chromosome') +
156. ylab('-log 10 pvalue')+
157. #scale_color_manual(values = rep(c("grey", "skyblue"), 22 )) +
158. scale_x_continuous(label = chrom_sizes[c(1:18,20,22),]$CHR, breaks= chrom_sizes[c(1:18,20,22),]$chr_label_loc ) +
159. #scale_y_continuous(expand = c(0, 0),limits = c(0, (max(-log10(twas2$pvalue)) + 3)) ) + # remove space between plot area and x axis
160. scale_color_discrete_sequential(palette = 'Batlow') +
161. #scale_color_manual(values=as.vector(watlington(15)[c(3,4,1,5,12)])) +
164. # Custom the theme:
165. #guides(color = guide_legend(title = "Training Model", override.aes = list(size = 3), nrow=2, byrow=FALSE)) +
166. guides(color = guide_legend(nrow=3, byrow=FALSE)) +
167. theme_minimal()+
168. theme(
169. legend.position="bottom",
170. panel.border = element_blank(),
171. panel.grid.major.x = element_blank(),
172. panel.grid.minor.x = element_blank(),
173. axis.title.x = element_blank(),
174. plot.title = element_text(hjust = 0.5),
175. legend.text=element_text()
176. )
178. name = 'Output/all_assocs_nashplot_pgc.pdf'
179. pdf(file = name, width = 15, height = 8, pointsize = 12, bg = "white")
180. print(fig2)
181. dev.off()
184. print(fig2)
185. ```
188. ### Fig 1b: Gene sharing in models
190. Certain genes are only assessed in G2S or GTEx but not both. We noted in computation
191. analyses that there are far more statistically significant results from the G2S
192. models than the GTEx. We want to see if differential gene characterization is
193. unerlying this feature. To answer this question, we go to the original predictDB
194. models for G2S and GTEx and examine the genes characterized in all the models.
196. Then we can take the genes that were sgnificantly associated with a disease and
197. characterize those genes according to if they are captured in a G2S model, GTEx
198. GTEx model, or both.
200. ```{r fig.height=7, fig.width = 4.5}
201. load('Data/dbdata.r')
202. cotested <- dbdata$cotested
203. g2s_only <- dbdata$g2s_only
204. gtx_only <- dbdata$gtx_only
209. dt_ <- fread('Data/all_exc_BHsig_assocs_PGCancestryTWAS.txt', header = TRUE)
211. # exclude replication/comparison data from initial analysis.
212. gwas <- c('PGC_AA_MDD_2023', 'PGC_AA_PTSD_2024','PGC_AA_SCZ_2022','PGC_ASIAN_SCZ_2022','PGC_EUR_PTSD_2024','PGC_HNA_PTSD_2024', 'PGC_EURO_SCZ_2022')
213. mod <- c('PrediXcan_Brain_Cortex')
215. dt <- dt_[!(gwas_phenotype %in% gwas | training_model %in% mod),]
218. dt[training_model == 'All_Whole_Blood', `gene expression model` := 'G2S All WB']
219. dt[training_model == 'AA_Whole_Blood', `gene expression model` := 'G2S African American WB']
220. dt[training_model == 'MX_Whole_Blood', `gene expression model` := 'G2S Mexican American WB']
221. dt[training_model == 'PR_Whole_Blood', `gene expression model` := 'G2S Puerto Rican WB']
222. dt[training_model == 'PrediXcan_Whole_Blood', `gene expression model` := 'GTEx European WB']
225. tms <- rev(c('G2S All WB', 'G2S African American WB', 'G2S Mexican American WB', 'G2S Puerto Rican WB', 'GTEx European WB'))
227. dt$`gene expression model` <- factor(dt$`gene expression model`, levels = tms)
230. dt[`gene expression model` == 'GTEx European WB', model := 'GTEX WB']
231. dt[`gene expression model` != 'GTEx European WB', model := 'G2S']
234. z1 <- unique(dt[, c('gene', 'gwas_phenotype', 'model')])
235. z2 <- z1[, .(.N), by = c('gene', 'gwas_phenotype')]
236. z2[N == 2, typ := 'co-sig']
237. z2[N==1 & gene %in% cotested, typ := 'co-tested']
238. z2[N==1 & !(gene %in% cotested), typ := 'only']
240. z3 <- setDT(merge(z1, z2, by = c('gene', 'gwas_phenotype')))
242. z3[typ != 'co-sig', category := paste(model, typ, sep = ' ')]
243. z3[typ == 'co-sig', category := 'co-significant']
245. z4 <- z3[, .(.N), by = c('gwas_phenotype', 'model', 'category')]
248. z4[gwas_phenotype == 'PGC_ALL_PTSD_2024', dz := 'PTSD']
249. z4[gwas_phenotype == 'PGC_all_SCZ_2022', dz := 'SCZ']
250. z4[gwas_phenotype == 'PGC_allADHD_2022', dz := 'ADHD']
251. z4[gwas_phenotype == 'PGC_allMDD_2023', dz := 'MDD']
252. z4[gwas_phenotype == 'PGC_BD1_2021', dz := 'BD1']
253. z4[gwas_phenotype == 'SUD_alc_2019', dz := 'AUD']
255. z4$category = factor(z4$category, levels = c("G2S only", "G2S co-tested", "GTEX WB only", "GTEX WB co-tested", "co-significant"))
257. tm_gp_plot6 <- ggplot(z4, aes(x = model, y = N, fill = category)) +
258. geom_bar(stat = 'identity') +
259. facet_wrap(vars(dz), scales = 'free', nrow = 3) +
260. #scale_fill_manual(values=as.vector(watlington(15)[c(3,4,5,12,1)])) +
261. scale_fill_manual(values=c('#cc6576ff', 'deeppink4', 'darkblue', 'cyan3', 'grey30')) +
262. theme_minimal() +
263. theme(legend.position = 'bottom', legend.title=element_blank()) +
264. guides(fill = guide_legend(nrow = 2))
266. name = 'Output/all_pgc_ancestry_TMxGP_bar_formatted2g.pdf'
267. pdf(file = name, width = 4.5, height = 7, pointsize = 12, bg = "white")
268. print(tm_gp_plot6)
269. dev.off()
272. print(tm_gp_plot6)
274. ```
279. ### Table of gene overlap across models
281. Uncomment for results displayed as table
282. ```{r}
284. #kable(z4, caption = "Overlap data for genes included in GTEx vs G2S models", align = "ccc", digits = 2) %>%
285. # kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = F) %>%
286. # column_spec(1, bold = T) %>%
287. # row_spec(0, bold = T, color = "white", background = "#D7261E")%>%
288. # scroll_box(width = "100%", height = "400px") # Adjust height as needed
290. ```
294. ## SNP level analyses
296. ### Fig 2e, S12: Correlation of SNP weights
298. Each model has SNP weights which serve as predictors of gene expression. We can take the
299. same genes and get all the shared SNP predictors across G2S and GTEx and characterize
300. the effect size of the weight on the gene. This will tell us if the models expect the same
301. SNPs to have different effects on the same genes.
303. ```{r fig.height=11, fig.width = 9}
308. all_dat <- list.files('Data/rdata', pattern = '.+\\.txt')
310. duo_snps <- data.table()
311. for(fil in str_subset(all_dat, '^snpfx_dat_all.+\\.txt')){
312. newdat <- fread(paste0('Data/rdata/', fil), header = T)
313. duo_snps <- rbind(duo_snps, newdat)
314. }
316. r2_values <- duo_snps[, .(r2 = summary(lm(weight.GTEx ~ weight.G2S))$adj.r.squared), by = .(gwas_phenotype, model)]
317. r_values <- duo_snps[, .(r_val = cor.test(weight.GTEx, weight.G2S)$estimate), by = .(gwas_phenotype, model)]
319. fx_shared <- ggplot(data = duo_snps, aes(x = weight.G2S, y = weight.GTEx)) +
320. geom_point() +
321. geom_smooth(method = 'lm', se=F) +
322. geom_text(data = r_values, aes(x = Inf, y = -Inf, label = paste0("R = ", round(r_val, 3))),
323. hjust = 1.1, vjust = -0.5, check_overlap = TRUE) +
324. xlab(paste0('SNP weight'))+
325. ylab(paste0('SNP weight GTEx')) +
326. facet_grid(gwas_phenotype~model) +
327. theme(legend.position = 'none') +
328. theme_bw()
330. name = 'Output/fxcorr_grid24_cor.pdf'
331. pdf(file = name, width = 9, height = 11, pointsize = 12, bg = "white")
332. print(fx_shared)
333. dev.off()
336. print(fx_shared)
339. ```
341. ### Correlation of gene zscores
343. While the previous analysis examined the model specific effects of SNP weights on gene expression
344. this analysis examines the relationship between GReX on disease as characterized by GTEx vs G2S.
345. We look at the zscores for the relationship of Gene X on disease Y as modeled by the different TWAS.
346. This tells us if the models predict that the effect of genes on disease differs across models.
348. ```{r}
350. dsgene <- data.table()
351. for(fil in str_subset(all_dat, '^dsgene_dat_all.+\\.txt')){
352. newdat <- fread(paste0('Data/rdata/', fil), header = T)
353. dsgene <- rbind(dsgene, newdat)
354. }
356. ```
358. ```{r fig.height=12.5, fig.width = 10}
360. cor_vals <- dsgene[, .(Zcor = cor(zscore.G2S, zscore.GTEx, use="complete.obs"), Pcor = cor(pvalue.G2S, pvalue.GTEx, use="complete.obs")), by = c( 'gwas_phenotype', 'model')]
361. dg2 <- setDT(merge(dsgene, cor_vals, by = c('gwas_phenotype', 'model')))
364. zsc_plot <- ggplot(data = dg2, aes(x = zscore.G2S, y = zscore.GTEx)) +
365. #geom_point() +
366. geom_pointdensity(adjust = 4) +
367. scale_color_viridis() +
368. geom_smooth(method = 'lm', se=F, color = 'black') +
369. geom_text(aes(x = Inf, y = Inf, label = paste0("R = ", round(Zcor, 3))),
370. hjust = 1.1, vjust = -0.3, check_overlap = TRUE, color = 'black', size = 3) +
371. facet_grid2(gwas_phenotype~model, axes = 'all') +
372. #facet_grid(gwas_phenotype~model) +
373. xlab('zscore.G2S') +
374. ylab('zscore.GTEx') +
375. coord_cartesian( clip = "off")+
376. theme_minimal()+
377. theme(panel.spacing = unit(1, "cm", data = NULL), strip.text.x = element_text(margin = margin(5, 5, 12, 5, "pt")))
380. name = 'Output/zsc_corr_grid24.pdf'
381. pdf(file = name, width = 10, height = 12.5, pointsize = 12, bg = "white")
382. print(zsc_plot)
383. dev.off()
385. print(zsc_plot)
386. ```
388. ### Correlation of gene zscores 2 (Fig 2b, S6-7)
390. Repeat the analysis but instead of a density plot, use black points with red trendline.
391. ```{r fig.height=12.5, fig.width = 10}
393. zsc_plotblk <- ggplot(data = dg2, aes(x = zscore.G2S, y = zscore.GTEx)) +
394. geom_point() +
395. #geom_pointdensity(adjust = 4) +
396. #scale_color_viridis() +
397. geom_smooth(method = 'lm', se=F, color = 'red') +
398. geom_text(aes(x = -Inf, y = Inf, label = paste0("R = ", round(Zcor, 3))),
399. hjust = 0.1, vjust = -0.3, check_overlap = TRUE, color = 'red', size = 3) +
400. facet_grid2(gwas_phenotype~model, axes = 'all') +
401. #facet_grid(gwas_phenotype~model) +
402. xlab('zscore.G2S') +
403. ylab('zscore.GTEx') +
404. coord_cartesian( clip = "off")+
405. theme_minimal()+
406. theme(panel.spacing = unit(1, "cm", data = NULL), strip.text.x = element_text(margin = margin(5, 5, 12, 5, "pt")))
409. name = 'Output/zsc_corr_grid24_black.pdf'
410. pdf(file = name, width = 10, height = 12.5, pointsize = 12, bg = "white")
411. print(zsc_plotblk)
412. dev.off()
414. print(zsc_plotblk)
416. ```
418. ### Correlation of pvalues
420. Similar to the correlation between zscores, we can also examine the relationship between p-values.
421. This will tell us if the different models call the same gene-disease association with differing
422. degrees of statistical significance.
424. ```{r fig.height=12.5, fig.width = 10}
425. pval_plot <- ggplot(data = dg2, aes(x = pvalue.G2S, y = pvalue.GTEx)) +
426. #geom_point() +
427. geom_pointdensity(adjust = 4) +
428. scale_color_viridis() +
429. geom_smooth(method = 'lm', se=F, color = 'black') +
430. geom_text(aes(x = Inf, y = Inf, label = paste0("R = ", round(Pcor, 3))),
431. hjust = 1.1, vjust = -0.3, check_overlap = TRUE, color = 'black', size = 3) +
432. scale_fill_manual(values=as.vector(kelly()[8:11])) +
433. facet_grid2(gwas_phenotype~model, axes = 'all') +
434. #facet_grid(gwas_phenotype~model) +
435. xlab('pvalue.G2S') +
436. ylab('pvalue.GTEx') +
437. coord_cartesian( clip = "off")+
438. theme_minimal() +
439. theme(panel.spacing = unit(1, "cm", data = NULL), strip.text.x = element_text(margin = margin(5, 5, 12, 5, "pt")))
442. name = 'Output/pval_corr_grid24.pdf'
443. pdf(file = name, width = 10, height = 12.5, pointsize = 12, bg = "white")
444. print(pval_plot)
445. dev.off()
447. print(pval_plot)
448. ```
450. ### Correlation of p-values 2 (Figure 2c, S8-9)
452. And again, repeat with black points and red trendline.
454. ```{r fig.height=12.5, fig.width = 10}
455. pval_plotblk <- ggplot(data = dg2, aes(x = pvalue.G2S, y = pvalue.GTEx)) +
456. geom_point() +
457. #geom_pointdensity(adjust = 4) +
458. #scale_color_viridis() +
459. geom_smooth(method = 'lm', se=F, color = 'red') +
460. geom_text(aes(x = -Inf, y = Inf, label = paste0("R = ", round(Pcor, 3))),
461. hjust = 0.1, vjust = -0.3, check_overlap = TRUE, color = 'red', size = 3) +
462. scale_fill_manual(values=as.vector(kelly()[8:11])) +
463. facet_grid2(gwas_phenotype~model, axes = 'all') +
464. #facet_grid(gwas_phenotype~model) +
465. xlab('pvalue.G2S') +
466. ylab('pvalue.GTEx') +
467. coord_cartesian( clip = "off")+
468. theme_minimal() +
469. theme(panel.spacing = unit(1, "cm", data = NULL), strip.text.x = element_text(margin = margin(5, 5, 12, 5, "pt")))
472. name = 'Output/pval_corr_grid24_black.pdf'
473. pdf(file = name, width = 10, height = 12.5, pointsize = 12, bg = "white")
474. print(pval_plotblk)
475. dev.off()
477. print(pval_plotblk)
480. ```
483. ### Fig S12:Comparison of r2 values for GTeX vs G2S
485. Each gene that is analyzed in the TWAS relies on a number of SNP weights that are used as predictors.
486. The PredictDB pipeline that is used to make the model includes information about the correlation
487. between the predicted gene expression using its SNP weights and the measured gene expression from
488. RNA seq in the original dataset. This R2 is thus a measure of the predictive performance of the model
489. regarding that specific gene. The difference in TWAS results may be related to the accuracy with which
490. the different model are predicting gene expression. As a result, the code below compares the R2 values
491. for genes across the different models that emerged from the TWAS.
493. ```{r fig.height=12.5, fig.width = 10}
494. r2corr <- ggplot(data = dsgene, aes(x = pred_perf_r2.G2S, y = pred_perf_r2.GTEx)) +
495. geom_point() +
496. geom_abline(intercept = 0, slope = 1, color = 'grey80', linetype = "dashed")+
497. geom_smooth(method = 'lm', se=F) +
498. xlab(paste0('R2 G2S'))+
499. ylab(paste0('R2 GTEx')) +
500. facet_grid(gwas_phenotype~model) +
501. theme(legend.position = 'none') +
502. theme_bw()
504. name = 'Output/r2_mod_corr_grid24.pdf'
505. pdf(file = name, width = 9, height = 11, pointsize = 12, bg = "white")
506. print(r2corr)
507. dev.off()
509. print(r2corr)
511. ```
513. ### Correlation of zscores by gene significance
515. Earlier we examined the total zscore correlations across models. It is possible to examine these
516. in more granularity, splitting categories according to if the genes of interest were deemed significant
517. in one, both, or none of the models.
519. ```{r fig.height=5, fig.width = 5}
520. fdga <- dsgene[, .(Ngroups = sum(!is.na(zscore.G2S) & !is.na(zscore.GTEx))), by = c('eqtl_source', 'twas_specificity', 'gwas_phenotype', 'model')]
522. fdsgene2 <- setDT(merge(dsgene[!is.na(zscore.G2S) & !is.na(zscore.GTEx),], fdga, by = c('eqtl_source','twas_specificity', 'gwas_phenotype', 'model')))
523. ffx_corr <- fdsgene2[Ngroups >2,
524. .(fx_r2 = cor.test(zscore.G2S, zscore.GTEx)$estimate,
525. fx_r2_ciL = cor.test(zscore.G2S, zscore.GTEx)$conf.int[1],
526. fx_r2_ciH = cor.test(zscore.G2S, zscore.GTEx)$conf.int[2],
527. N = .N), by = c('eqtl_source','twas_specificity', 'gwas_phenotype', 'model')]
528. modnm = c('G2S_aa', 'G2S_all', 'G2S_mx', 'G2S_pr')
529. ffx_corr$twas_specificity <- factor(ffx_corr$twas_specificity, levels = c('none', 'GTEx_WB', modnm, 'shared'))
532. ffx_corr[twas_specificity == 'shared', spec := 'G2S & GTEx']
533. ffx_corr[twas_specificity == 'none', spec := 'none']
534. ffx_corr[!(twas_specificity %in% c('shared','none')), spec := 'G2S or GTEx']
535. fx2 <- ffx_corr[, .(mfxr = mean(fx_r2)), by = c('eqtl_source', 'spec')]
538. fx2_cord <- ggplot(data = fx2, aes(x = spec, y = mfxr, color = eqtl_source, group = eqtl_source)) +
539. geom_hline(yintercept = 0, color = 'red') +
540. #geom_errorbar(aes(ymin = fx_r2_ciL, ymax = fx_r2_ciH), width = 0.1, color = 'black', position=position_dodge(width=0.3)) +
541. geom_point(position=position_dodge(width=0.3), size = 2) +
542. scale_color_manual(values=c('green4', 'slateblue3')) +
543. ylab('effect size correlation') +
544. xlab('twas model significance') +
545. ylim(c(-1, 1))+
546. #facet_grid(gwas_phenotype~model, scales = 'free') +
547. theme_bw() +
548. theme(axis.text.x = element_text(angle = 60, hjust = 1))
552. name = 'Output/genefx_corplot_grid24d_all.pdf'
553. pdf(file = name, width = 5, height = 5, pointsize = 12, bg = "white")
554. print(fx2_cord)
555. dev.off()
557. print(fx2_cord)
559. ```
561. ### Fig S16: correlation of zscores by gene significance 2
563. We extend the analyses from above but split across models and clinical phenotypes
565. ```{r fig.height=11, fig.width = 13}
568. fdga <- dsgene[, .(Ngroups = sum(!is.na(zscore.G2S) & !is.na(zscore.GTEx))), by = c('eqtl_source', 'twas_specificity', 'gwas_phenotype', 'model')]
570. fdsgene2 <- setDT(merge(dsgene[!is.na(zscore.G2S) & !is.na(zscore.GTEx),], fdga, by = c('eqtl_source','twas_specificity', 'gwas_phenotype', 'model')))
571. ffx_corr <- fdsgene2[Ngroups >2,
572. .(fx_r2 = cor.test(zscore.G2S, zscore.GTEx)$estimate,
573. fx_r2_ciL = cor.test(zscore.G2S, zscore.GTEx)$conf.int[1],
574. fx_r2_ciH = cor.test(zscore.G2S, zscore.GTEx)$conf.int[2],
575. N = .N), by = c('eqtl_source','twas_specificity', 'gwas_phenotype', 'model')]
576. modnm = c('G2S_aa', 'G2S_all', 'G2S_mx', 'G2S_pr')
577. ffx_corr$twas_specificity <- factor(ffx_corr$twas_specificity, levels = c('none', 'GTEx_WB', modnm, 'shared'))
580. fx_cord <- ggplot(data = ffx_corr, aes(x = twas_specificity, y = fx_r2, color = eqtl_source, group = eqtl_source)) +
581. geom_hline(yintercept = 0, color = 'red') +
582. geom_errorbar(aes(ymin = fx_r2_ciL, ymax = fx_r2_ciH), width = 0.1, color = 'black', position=position_dodge(width=0.3)) +
583. geom_point(position=position_dodge(width=0.3), size = 2) +
584. scale_color_manual(values=c('green4', 'slateblue3')) +
585. ylab('effect size correlation') +
586. xlab('twas model significance') +
587. ylim(c(-1, 1))+
588. facet_grid(gwas_phenotype~model, scales = 'free') +
589. theme_bw() +
590. theme(axis.text.x = element_text(angle = 60, hjust = 1))
594. name = 'Output/genefx_corplot_grid24d.pdf'
595. pdf(file = name, width = 13, height = 11, pointsize = 12, bg = "white")
596. print(fx_cord)
597. dev.off()
602. print(fx_cord)
605. ```
609. ### Fig S17: correlation of pvalues by gene significance
611. And again, we can examine the same analysis according to pvalue to get a visual representation
612. of how the models handle statistical significance.
614. ```{r fig.height=11, fig.width = 13}
616. pdga <- dsgene[, .(Ngroups = sum(!is.na(pvalue.G2S) & !is.na(pvalue.GTEx))), by = c('eqtl_source', 'twas_specificity', 'gwas_phenotype', 'model')]
618. pdsgene2 <- setDT(merge(dsgene[!is.na(pvalue.G2S) & !is.na(pvalue.GTEx),], pdga, by = c('eqtl_source','twas_specificity', 'gwas_phenotype', 'model')))
619. pfx_corr <- pdsgene2[Ngroups >2,
620. .(pv_r = cor.test(pvalue.G2S, pvalue.GTEx)$estimate,
621. pv_r_ciL = cor.test(pvalue.G2S, pvalue.GTEx)$conf.int[1],
622. pv_r_ciH = cor.test(pvalue.G2S, pvalue.GTEx)$conf.int[2],
623. N = .N), by = c('eqtl_source','twas_specificity', 'gwas_phenotype', 'model')]
624. modnm = c('G2S_aa', 'G2S_all', 'G2S_mx', 'G2S_pr')
625. pfx_corr$twas_specificity <- factor(pfx_corr$twas_specificity, levels = c('none', 'GTEx_WB', modnm, 'shared'))
628. pv_cor <- ggplot(data = pfx_corr, aes(x = twas_specificity, y = pv_r, color = eqtl_source, group = eqtl_source)) +
629. geom_hline(yintercept = 0, color = 'red') +
630. geom_errorbar(aes(ymin = pv_r_ciL, ymax = pv_r_ciH), width = 0.1, color = 'black', position=position_dodge(width=0.3)) +
631. geom_point(position=position_dodge(width=0.3), size = 2) +
632. scale_color_manual(values=c('coral1', 'deepskyblue2')) +
633. #geom_label(aes(y = -0.75 + 0.25*sum('shared_distinct' %in% eqtl_source), label = paste('N =', N))) +
634. #geom_label(aes(y = -1 , label = paste(N), color = eqtl_source), position = position_dodgev(height = -0.7)) +
635. ylab('pvalue correlation') +
636. xlab('twas model significance') +
637. ylim(c(-1, 1))+
638. facet_grid(gwas_phenotype~model, scales = 'free') +
639. theme_bw() +
640. theme(axis.text.x = element_text(angle = 60, hjust = 1))
644. name = 'Output/genepv_corplot_grid24.pdf'
645. pdf(file = name, width = 13, height = 11, pointsize = 12, bg = "white")
646. print(pv_cor)
647. dev.off()
649. print(pv_cor)
651. ```
654. ### Fig 3b,S15: Median SNP weight magnitudes
656. The SNPs used for the analysis may be different across the models but overall, I'm interested in
657. seeing if the different models are able to assign SNPs of different effect sizes to the genes. To
658. capture this, the code below characterizes the median magnitude of the SNP weights across G2S and
659. GTEx.
661. ```{r fig.height=11, fig.width = 9}
665. all_snps <- data.table()
666. for(fil in str_subset(all_dat, '^snp_dat_all.+\\.txt')){
667. newdat <- fread(paste0('Data/rdata/', fil), header = T)
668. all_snps <- rbind(all_snps, newdat)
669. }
671. wts <- setDT(melt(all_snps[, c('model', 'gwas_phenotype', 'weight.G2S', 'weight.GTEx')], id.vars = c('model', 'gwas_phenotype')))
673. eqtl_medwtplot <- ggplot(data = wts, aes(x = variable, y = abs(value), fill = variable)) +
674. geom_boxplot(outlier.shape = NA) +
675. xlab('model')+
676. ylab('median weight magnitude') +
677. theme_minimal()+
678. facet_grid(gwas_phenotype~model, scales = 'free') +
679. coord_cartesian(ylim= c(0,0.15)) +
680. theme(legend.position = 'none') +
681. theme_bw()
683. name = 'Output/eqtl_medwtplot_grid24.pdf'
684. pdf(file = name, width = 9, height = 11, pointsize = 12, bg = "white")
685. print(eqtl_medwtplot)
686. dev.off()
688. print(eqtl_medwtplot)
690. ```
693. ### Fig 2f: SNP, Gene, and Sig Gene correlations
695. Let's group the zscore and SNP weight correlations into a single figure for easier comparitive
696. analysis.
698. ```{r fig.width = 5.5, fig.height = 11}
700. wt_r2_values <- duo_snps[, .(r2_SNPpredictors = summary(lm(weight.GTEx ~ weight.G2S))$adj.r.squared, wt_rse = summary(lm(weight.GTEx ~ weight.G2S))$sigma, N_r2_SNP = .N), by = .(gwas_phenotype, model)]
701. #wt_r2_valuessig <- duo_sig[, .(r2_SNPpredictors = summary(lm(weight.GTEx ~ weight.G2S))$adj.r.squared, wt_rse = summary(lm(weight.GTEx ~ weight.G2S))$sigma, N_r2_SNP = .N), by = .(gwas_phenotype, model)]
706. wt_r_values <- duo_snps[, .(r_SNPpredictors = cor.test(weight.GTEx,weight.G2S)$estimate,
707. wt_ciL = cor.test(weight.GTEx,weight.G2S)$conf.int[1],
708. wt_ciH = cor.test(weight.GTEx,weight.G2S)$conf.int[2],
709. N_wt_SNP = .N), by = .(gwas_phenotype, model)]
713. gnfx_r_values <- dsgene[, .(r_gnfx = cor.test(zscore.GTEx,zscore.G2S)$estimate,
714. gnfx_ciL = cor.test(zscore.GTEx,zscore.G2S)$conf.int[1],
715. gnfx_ciH = cor.test(zscore.GTEx,zscore.G2S)$conf.int[2],
716. N_gnfx = .N), by = .(gwas_phenotype, model)]
719. gnfx_r_valsig <- dsgene[twas_specificity != 'none', .(r_gnfx_sig = cor.test(zscore.GTEx,zscore.G2S)$estimate,
720. gnfx_sig_ciL = cor.test(zscore.GTEx,zscore.G2S)$conf.int[1],
721. gnfx_sig_ciH = cor.test(zscore.GTEx,zscore.G2S)$conf.int[2],
722. N_sig_gnfx = .N), by = .(gwas_phenotype, model)]
727. cordat <- setDT(merge(wt_r_values, merge(gnfx_r_values, gnfx_r_valsig, by = c('gwas_phenotype', 'model')), by = c('gwas_phenotype', 'model')))
731. wt_r_values <- duo_snps[, .(r = cor.test(weight.GTEx,weight.G2S)$estimate,
732. ciL = cor.test(weight.GTEx,weight.G2S)$conf.int[1],
733. ciH = cor.test(weight.GTEx,weight.G2S)$conf.int[2],
734. N = .N,
735. measure = 'SNP_predictor'), by = .(gwas_phenotype, model)]
739. gnfx_r_values <- dsgene[, .(r = cor.test(zscore.GTEx,zscore.G2S)$estimate,
740. ciL = cor.test(zscore.GTEx,zscore.G2S)$conf.int[1],
741. ciH = cor.test(zscore.GTEx,zscore.G2S)$conf.int[2],
742. N = .N,
743. measure = 'effect size'), by = .(gwas_phenotype, model)]
746. gnfx_r_valsig <- dsgene[twas_specificity != 'none', .(r = cor.test(zscore.GTEx,zscore.G2S)$estimate,
747. ciL = cor.test(zscore.GTEx,zscore.G2S)$conf.int[1],
748. ciH = cor.test(zscore.GTEx,zscore.G2S)$conf.int[2],
749. N = .N,
750. measure = 'effect size (sig)'), by = .(gwas_phenotype, model)]
754. cordat <- rbind(wt_r_values, gnfx_r_values, gnfx_r_valsig)
756. cordat <- rbind(wt_r_values, gnfx_r_values, gnfx_r_valsig)
758. cordat$measure <- factor(cordat$measure, levels = c('SNP_predictor', 'effect size', 'effect size (sig)'))
760. cordat[gwas_phenotype == 'PGC_BD1_2021', disease := 'BD1']
761. cordat[gwas_phenotype == 'PGC_all_SCZ_2022', disease := 'SCZ']
762. cordat[gwas_phenotype == 'SUD_alc_2019', disease := 'SUD-A']
763. cordat[gwas_phenotype == 'PGC_allADHD_2022', disease := 'ADHD']
764. cordat[gwas_phenotype == 'PGC_ALL_PTSD_2024', disease := 'PTSD']
765. cordat[gwas_phenotype == 'PGC_allMDD_2023', disease := 'MDD']
772. allcor_plot <- ggplot(data = cordat, aes(x = measure, y = r, group = model, color = model)) +
773. #geom_line(position=position_dodge(width=0.5))+
774. geom_errorbar(aes(ymin = ciL, ymax = ciH), width = 0.3, color = 'black', position=position_dodge(width=0.5)) +
775. geom_point(position=position_dodge(width=0.5)) +
776. #facet_grid2(disease~., axes = 'all') +
777. facet_grid(disease~.) +
778. guides(color = guide_legend(nrow = 2)) +
779. xlab('measurement') +
780. ylab('r correlation coefficient') +
781. ylim(c(0,1)) +
782. theme_bw() +
783. theme(axis.text = element_text(size = 14),
784. axis.title = element_text(size = 14),
785. strip.text = element_text(size = 14),
786. legend.text = element_text(size=14),
787. legend.title = element_text(size=14),
788. legend.position = 'bottom')
791. name = 'Output/allcordat_grid.pdf'
792. pdf(file = name, width = 5.5, height = 11, pointsize = 12, bg = "white")
793. print(allcor_plot)
794. dev.off()
796. print(allcor_plot)
798. ```
800. ### Fig 3e,S22: r2 by model specificity
802. As before, we split the zscore correlations by model specificity. Here we replicate this analytical design
803. with the r2 but instead of dots and bars for range, we'll visualize this with violin plots and imbedded boxplots.
804. ```{r fig.width = 9, fig.height = 11}
806. r2_dat <- data.table()
807. for(fil in str_subset(all_dat, '^r2_dat_all.+\\.txt')){
808. newdat <- fread(paste0('Data/rdata/', fil), header = T)
809. r2_dat <- rbind(r2_dat, newdat)
810. }
813. r2_vplot <- ggplot(data = r2_dat, aes(x = variable, y = value, fill = eqtl_source)) +
814. geom_violin(scale = 'width', position=position_dodge(1), drop = F) +
815. geom_boxplot(width = 0.1, color = 'black', position=position_dodge(1), outlier.size = 0.5) +
816. xlab('eQTL model specificity') +
817. ylab('R2 per gene') +
818. facet_grid(gwas_phenotype~model, scales = 'free') +
819. theme_minimal() +
820. theme(axis.text.x = element_text(angle = 60, hjust = 1))
822. name = 'Output/r2_vplot_grid24.pdf'
823. pdf(file = name, width = 9, height = 11, pointsize = 12, bg = "white")
824. print(r2_vplot)
825. dev.off()
827. print(r2_vplot)
829. ```
831. ### r2 plot agregated across all factors
833. ```{r}
835. duo_sig <- data.table()
836. for(fil in str_subset(all_dat, '^fxsig_dat_all.+\\.txt')){
837. newdat <- fread(paste0('Data/rdata/', fil), header = T)
838. duo_sig <- rbind(duo_sig, newdat)
839. }
841. r2_values <- duo_sig[, .(r2 = summary(lm(weight.GTEx ~ weight.G2S))$adj.r.squared), by = .(gwas_phenotype, model)]
842. r_values <- duo_sig[, .(r_val = cor.test(weight.GTEx, weight.G2S)$estimate), by = .(gwas_phenotype, model)]
844. fx_sharedsig <- ggplot(data = duo_sig, aes(x = weight.G2S, y = weight.GTEx, color = model_specificity)) +
845. geom_point() +
846. geom_smooth(method = 'lm', se=F, weight = 1, fullrange = T) +
847. geom_text(inherit.aes = F, data = r_values, aes(x = Inf, y = -Inf, label = paste0("R = ", round(r_val, 3))),
848. hjust = 1.1, vjust = -0.5, check_overlap = TRUE) +
849. xlab(paste0('SNP weight'))+
850. ylab(paste0('SNP weight GTEx')) +
851. facet_grid(gwas_phenotype~model, scales = 'free') +
852. theme(legend.position = 'none') +
853. theme_bw()
855. name = 'Output/r2_aggregate_1.pdf'
856. pdf(file = name, width = 4, height = 4, pointsize = 12, bg = "white")
857. print(fx_sharedsig)
858. dev.off()
860. print(fx_sharedsig)
862. ```
866. ### Fig 2d, S11: Number of used eQTL's per gene
869. How many eQTLs were used for each gene that was tested in the different TWAS analyses?
871. ```{r fig.height=12.5, fig.width = 10}
874. eqtl_dat <- data.table()
875. for(fil in str_subset(all_dat, '^eqtl_dat_all.+\\.txt')){
876. newdat <- fread(paste0('Data/rdata/', fil), header = T)
877. eqtl_dat <- rbind(eqtl_dat, newdat)
878. }
881. eqtl_vplot <- ggplot(data = eqtl_dat, aes(x = variable, y = value, fill = variable)) +
882. geom_violin(scale = 'width') +
883. geom_boxplot(width = 0.1, color = 'black', fill = 'white',outlier.shape = NA) +
884. xlab('eQTL model specificity')+
885. ylab('Number of eQTLs per gene') +
886. theme_minimal()+
887. facet_grid(gwas_phenotype~model, scales = 'free') +
888. theme_bw() +
889. theme(legend.position = 'none', axis.text.x = element_text(angle = 60, hjust = 1))
891. name = 'Output/eqtl_vplot_grid24.pdf'
892. pdf(file = name, width = 9, height = 11, pointsize = 12, bg = "white")
893. print(eqtl_vplot)
894. dev.off()
896. print(eqtl_vplot)
900. ```
902. ### Fig 3c, S13: r2 by model
904. split r2 by model instead of significance or correlation.
906. ```{r fig.height=12.5, fig.width = 10}
909. preds <- melt(dsgene[, c('gwas_phenotype', 'gene', 'model', 'pred_perf_r2.G2S', 'pred_perf_r2.GTEx')], id.vars = c('gwas_phenotype', 'gene', 'model' ))
910. r2_vplot <- ggplot(data = preds, aes(x = variable, y = value, fill = variable)) +
911. geom_violin(scale = 'width') +
912. geom_boxplot(width = 0.1, color = 'black', fill = 'white',outlier.shape = NA) +
913. xlab('model')+
914. ylab('Prediction performance r2') +
915. theme_minimal()+
916. facet_grid(gwas_phenotype~model, scales = 'free') +
917. theme_bw() +
918. theme(legend.position = 'none', axis.text.x = element_text(angle = 60, hjust = 1))
921. name = 'Output/predperf_vplot_grid24.pdf'
922. pdf(file = name, width = 9, height = 11, pointsize = 12, bg = "white")
923. print(r2_vplot)
924. dev.off()
926. print(r2_vplot)
930. ```
932. ## MDD Ancestry-concordant TWAS
934. We wanted to see if the results we got were the specific to a study design where we applied two different models
935. to the same GWAS summary statistics. We pulled GWAS of MDD done in predominantly African Americans vs Europeans
936. and then applied more ancestry concordant models to each for the TWAS. (G2S_AA + African American MDD GWAS and then
937. GTEx + Eur GWAS). For this we repeated a number of the descriptive analyzes performed above.
939. ### zscore corrrelation all
941. Correlation of Zscores. I performed this multiple times considering different statistical significance thresholds.
943. ```{r}
944. library(data.table)
945. library(ggplot2)
946. library(pals)
947. library(gridExtra)
948. library(ggh4x)
949. library(DBI)
952. # get the correlations needed between the different models and twas analyses
954. mda <- fread('Data/dual_data/all_exc_assocs_PGCancestryTWAS_AA.txt', header = T)[gwas_phenotype == 'PGC_allMDD_2023',]
955. mde <- fread('Data/dual_data/all_exc_assocs_PGCancestryTWAS_EUR.txt', header = T)
957. mdeg <- unique(mde$gene_name)
958. mdag <- unique(mda$gene_name)
960. uniq_mde <- setdiff(mdeg, mdag)
961. uniq_mda <- setdiff(mdag, mdeg)
964. jointmdd <- setDT(merge(mda[training_model != 'PrediXcan_Brain_Cortex',], mde, by = c('gene_name', 'training_model'), suffixes = c('.multi', '.eur')))
966. cor.test(jointmdd$zscore.multi, jointmdd$zscore.eur)
971. ```
973. ### zscore cor all shared
975. ```{r fig.width = 8, fig.height = 2.6}
977. # correlation between zscores across the two cohorts is 0.803 across all models.
978. mdcor <- jointmdd[, .(zsc_cor = cor.test(zscore.multi, zscore.eur)$estimate), by = c('training_model')]
979. mdcorsig <- jointmdd[bhpval.eur <= 0.05 & bhpval.multi <= 0.05, .(zsc_cor = cor.test(zscore.multi, zscore.eur)$estimate), by = c('training_model')]
980. mdcorsig2 <- jointmdd[bhpval.eur <= 0.05 | bhpval.multi <= 0.05, .(zsc_cor = cor.test(zscore.multi, zscore.eur)$estimate), by = c('training_model')]
984. crossmdd <- setDT(merge(mda[training_model != 'PrediXcan_Brain_Cortex' & training_model != 'PrediXcan_Whole_Blood',], mde[training_model == 'PrediXcan_Whole_Blood',], by = c('gene_name'), suffixes = c('.multi_g2S', '.eur_gtex')))
986. cmdcor <- crossmdd[, .(zsc_cor = cor.test(zscore.multi_g2S, zscore.eur_gtex)$estimate,
987. zsc_ciL = cor.test(zscore.multi_g2S, zscore.eur_gtex)$conf.int[1],
988. zsc_ciH = cor.test(zscore.multi_g2S, zscore.eur_gtex)$conf.int[2],
989. .N), by = c('training_model.multi_g2S')]
990. cmdcorsig <- crossmdd[bhpval.eur_gtex <= 0.05 & bhpval.multi_g2S <= 0.05, .(zsc_cor = cor.test(zscore.multi_g2S, zscore.eur_gtex)$estimate,
991. zsc_ciL = cor.test(zscore.multi_g2S, zscore.eur_gtex)$conf.int[1],
992. zsc_ciH = cor.test(zscore.multi_g2S, zscore.eur_gtex)$conf.int[2],
993. .N), by = c('training_model.multi_g2S')]
994. cmdcorsig2 <- crossmdd[bhpval.eur_gtex <= 0.05 | bhpval.multi_g2S <= 0.05, .(zsc_cor = cor.test(zscore.multi_g2S, zscore.eur_gtex)$estimate,
995. zsc_ciL = cor.test(zscore.multi_g2S, zscore.eur_gtex)$conf.int[1],
996. zsc_ciH = cor.test(zscore.multi_g2S, zscore.eur_gtex)$conf.int[2],
997. .N), by = c('training_model.multi_g2S')]
1001. ## Make plots
1003. pdat1 <- setDT(merge(crossmdd, cmdcor, by = 'training_model.multi_g2S'))
1005. cmd_genecorr <- ggplot(data = pdat1, aes(x = zscore.multi_g2S, y=zscore.eur_gtex)) +
1006. geom_point() +
1007. #geom_pointdensity(adjust = 4) +
1008. #scale_color_viridis() +
1009. geom_smooth(method = 'lm', se=F, color = 'green4') +
1010. geom_text(aes(x = -Inf, y = Inf, label = paste0("R = ", round(zsc_cor, 3))),
1011. hjust = 0.1, vjust = -0.3, check_overlap = TRUE, color = 'green4', size = 3) +
1012. facet_grid2(.~training_model.multi_g2S, axes = 'all') +
1013. ggtitle('All shared genes') +
1014. xlab('zscore.G2S') +
1015. ylab('zscore.GTEx') +
1016. coord_cartesian( clip = "off")+
1017. theme_minimal()+
1018. theme(panel.spacing = unit(1, "cm", data = NULL),
1019. strip.text.x = element_text(margin = margin(5, 5, 12, 5, "pt")))
1022. name='Output/dual_zsc_corr_grid_all.pdf'
1023. pdf(file = name, width = 8, height = 2.6, pointsize = 12, bg = "white")
1024. print(cmd_genecorr)
1025. dev.off()
1027. print(cmd_genecorr)
1028. ```
1030. ### zscore cor all co-sig
1033. ```{r fig.width = 8, fig.height = 2.6}
1036. pdat2 <- setDT(merge(crossmdd[bhpval.eur_gtex <= 0.05 & bhpval.multi_g2S <= 0.05,], cmdcorsig, by = 'training_model.multi_g2S'))
1038. cmd_genecorr_sig <- ggplot(data = pdat2, aes(x = zscore.multi_g2S, y=zscore.eur_gtex)) +
1039. geom_point() +
1040. #geom_pointdensity(adjust = 4) +
1041. #scale_color_viridis() +
1042. geom_smooth(method = 'lm', se=F, color = 'springgreen3') +
1043. geom_text(aes(x = -Inf, y = Inf, label = paste0("R = ", round(zsc_cor, 3))),
1044. hjust = 0.1, vjust = -0.3, check_overlap = TRUE, color = 'springgreen3', size = 3) +
1045. facet_grid2(.~training_model.multi_g2S, axes = 'all') +
1046. ggtitle('All co-significant genes') +
1047. xlab('zscore.G2S') +
1048. ylab('zscore.GTEx') +
1049. coord_cartesian( clip = "off")+
1050. theme_minimal()+
1051. theme(panel.spacing = unit(1, "cm", data = NULL),
1052. strip.text.x = element_text(margin = margin(5, 5, 12, 5, "pt")))
1055. name='Output/dual_zsc_corr_grid_sig.pdf'
1056. pdf(file = name, width = 8, height = 2.6, pointsize = 12, bg = "white")
1057. print(cmd_genecorr_sig)
1058. dev.off()
1061. print(cmd_genecorr_sig)
1063. ```
1065. ### zscore cor all sig
1067. ```{r fig.width = 8, fig.height = 2.6}
1070. pdat3 <- setDT(merge(crossmdd[bhpval.eur_gtex <= 0.05 | bhpval.multi_g2S <= 0.05,], cmdcorsig2, by = 'training_model.multi_g2S'))
1072. cmd_genecorr_sig2 <- ggplot(data = pdat3, aes(x = zscore.multi_g2S, y=zscore.eur_gtex)) +
1073. geom_point() +
1074. #geom_pointdensity(adjust = 4) +
1075. #scale_color_viridis() +
1076. geom_smooth(method = 'lm', se=F, color = 'lightgreen') +
1077. geom_text(aes(x = -Inf, y = Inf, label = paste0("R = ", round(zsc_cor, 3))),
1078. hjust = 0.1, vjust = -0.3, check_overlap = TRUE, color = 'lightgreen', size = 3) +
1079. facet_grid2(.~training_model.multi_g2S, axes = 'all') +
1080. ggtitle('All significant genes') +
1081. xlab('zscore.G2S') +
1082. ylab('zscore.GTEx') +
1083. coord_cartesian( clip = "off")+
1084. theme_minimal()+
1085. theme(panel.spacing = unit(1, "cm", data = NULL),
1086. strip.text.x = element_text(margin = margin(5, 5, 12, 5, "pt")))
1089. name='Output/dual_zsc_corr_grid2_sig2.pdf'
1090. pdf(file = name, width = 8, height = 2.6, pointsize = 12, bg = "white")
1091. print(cmd_genecorr_sig2)
1092. dev.off()
1095. cmd_genecorr_sig2
1097. ```
1099. ### pval cor all shared 1
1100. ```{r fig.width = 8, fig.height = 2.6}
1103. cmdcor <- crossmdd[, .(zsc_cor = cor.test(pvalue.multi_g2S, pvalue.eur_gtex)$estimate,
1104. zsc_ciL = cor.test(pvalue.multi_g2S, pvalue.eur_gtex)$conf.int[1],
1105. zsc_ciH = cor.test(pvalue.multi_g2S, pvalue.eur_gtex)$conf.int[2],
1106. .N), by = c('training_model.multi_g2S')]
1107. cmdcorsig <- crossmdd[bhpval.eur_gtex <= 0.05 & bhpval.multi_g2S <= 0.05, .(zsc_cor = cor.test(pvalue.multi_g2S, pvalue.eur_gtex)$estimate,
1108. zsc_ciL = cor.test(pvalue.multi_g2S, pvalue.eur_gtex)$conf.int[1],
1109. zsc_ciH = cor.test(pvalue.multi_g2S, pvalue.eur_gtex)$conf.int[2],
1110. .N), by = c('training_model.multi_g2S')]
1111. cmdcorsig2 <- crossmdd[bhpval.eur_gtex <= 0.05 | bhpval.multi_g2S <= 0.05, .(zsc_cor = cor.test(pvalue.multi_g2S, pvalue.eur_gtex)$estimate,
1112. zsc_ciL = cor.test(pvalue.multi_g2S, pvalue.eur_gtex)$conf.int[1],
1113. zsc_ciH = cor.test(pvalue.multi_g2S, pvalue.eur_gtex)$conf.int[2],
1114. .N), by = c('training_model.multi_g2S')]
1118. ## Make plots
1120. pdat1 <- setDT(merge(crossmdd, cmdcor, by = 'training_model.multi_g2S'))
1122. cmd_genecorr <- ggplot(data = pdat1, aes(x = zscore.multi_g2S, y=zscore.eur_gtex)) +
1123. geom_point() +
1124. #geom_pointdensity(adjust = 4) +
1125. #scale_color_viridis() +
1126. geom_smooth(method = 'lm', se=F, color = 'gold') +
1127. geom_text(aes(x = -Inf, y = Inf, label = paste0("R = ", round(zsc_cor, 3))),
1128. hjust = 0.1, vjust = -0.3, check_overlap = TRUE, color = 'gold', size = 3) +
1129. facet_grid2(.~training_model.multi_g2S, axes = 'all') +
1130. ggtitle('All shared genes') +
1131. xlab('pvalue.G2S') +
1132. ylab('pvalue.GTEx') +
1133. coord_cartesian( clip = "off")+
1134. theme_minimal()+
1135. theme(panel.spacing = unit(1, "cm", data = NULL),
1136. strip.text.x = element_text(margin = margin(5, 5, 12, 5, "pt")))
1139. name='Output/dual_pval_corr_grid_all.pdf'
1140. pdf(file = name, width = 8, height = 2.6, pointsize = 12, bg = "white")
1141. print(cmd_genecorr)
1142. dev.off()
1144. print(cmd_genecorr)
1147. ```
1149. ### pval cor all co-sig
1151. And then again with pvalues.
1153. ```{r fig.width = 8, fig.height = 2.6}
1155. pdat2 <- setDT(merge(crossmdd[bhpval.eur_gtex <= 0.05 & bhpval.multi_g2S <= 0.05,], cmdcorsig, by = 'training_model.multi_g2S'))
1157. cmd_genecorr_sig <- ggplot(data = pdat2, aes(x = zscore.multi_g2S, y=zscore.eur_gtex)) +
1158. geom_point() +
1159. #geom_pointdensity(adjust = 4) +
1160. #scale_color_viridis() +
1161. geom_smooth(method = 'lm', se=F, color = 'yellow3') +
1162. geom_text(aes(x = -Inf, y = Inf, label = paste0("R = ", round(zsc_cor, 3))),
1163. hjust = 0.1, vjust = -0.3, check_overlap = TRUE, color = 'yellow3', size = 3) +
1164. facet_grid2(.~training_model.multi_g2S, axes = 'all') +
1165. ggtitle('All co-significant genes') +
1166. xlab('pvalue.G2S') +
1167. ylab('pvalue.GTEx') +
1168. coord_cartesian( clip = "off")+
1169. theme_minimal()+
1170. theme(panel.spacing = unit(1, "cm", data = NULL),
1171. strip.text.x = element_text(margin = margin(5, 5, 12, 5, "pt")))
1174. name='Output/dual_pval_corr_grid2_sig.pdf'
1175. pdf(file = name, width = 8, height = 2.6, pointsize = 12, bg = "white")
1176. print(cmd_genecorr_sig)
1177. dev.off()
1179. print(cmd_genecorr_sig)
1181. ```
1183. ### pval cor all sig
1185. ```{r fig.width = 8, fig.height = 2.6}
1188. pdat3 <- setDT(merge(crossmdd[bhpval.eur_gtex <= 0.05 | bhpval.multi_g2S <= 0.05,], cmdcorsig2, by = 'training_model.multi_g2S'))
1190. cmd_genecorr_sig2 <- ggplot(data = pdat3, aes(x = zscore.multi_g2S, y=zscore.eur_gtex)) +
1191. geom_point() +
1192. #geom_pointdensity(adjust = 4) +
1193. #scale_color_viridis() +
1194. geom_smooth(method = 'lm', se=F, color = 'lightgreen') +
1195. geom_text(aes(x = -Inf, y = Inf, label = paste0("R = ", round(zsc_cor, 3))),
1196. hjust = 0.1, vjust = -0.3, check_overlap = TRUE, color = 'goldenrod', size = 3) +
1197. facet_grid2(.~training_model.multi_g2S, axes = 'all') +
1198. ggtitle('All significant genes') +
1199. xlab('pvalue.G2S') +
1200. ylab('pvalue.GTEx') +
1201. coord_cartesian( clip = "off")+
1202. theme_minimal()+
1203. theme(panel.spacing = unit(1, "cm", data = NULL),
1204. strip.text.x = element_text(margin = margin(5, 5, 12, 5, "pt")))
1207. name='Output/dual_pval_corr_grid2_sig2.pdf'
1208. pdf(file = name, width = 8, height = 2.6, pointsize = 12, bg = "white")
1209. print(cmd_genecorr_sig2)
1210. dev.off()
1212. print(cmd_genecorr_sig2)
1214. ```
1216. ### Pval cor All shared
1217. ```{r fig.width = 8, fig.height = 2.6}
1219. snp_dat_all <- fread('Data/dual_data/dual_duo_snps_used.txt', header = T)
1221. cmd_snpcor <- snp_dat_all[, .(wt_cor = cor.test(weight.G2S, weight.GTEx)$estimate,
1222. zsc_ciL = cor.test(weight.G2S, weight.GTEx)$conf.int[1],
1223. zsc_ciH = cor.test(weight.G2S, weight.GTEx)$conf.int[2],
1224. .N), by = c('model')]
1228. cordat <- setDT(merge(snp_dat_all, cmd_snpcor, by = 'model'))
1230. cmd_snpcorr <- ggplot(data = cordat, aes(x = weight.G2S, y=weight.GTEx)) +
1231. geom_point() +
1232. #geom_pointdensity(adjust = 4) +
1233. #scale_color_viridis() +
1234. geom_smooth(method = 'lm', se=F, color = 'purple') +
1235. geom_text(aes(x = -Inf, y = Inf, label = paste0("R = ", round(wt_cor, 3))),
1236. hjust = 0.1, vjust = -0.3, check_overlap = TRUE, color = 'purple', size = 3) +
1237. facet_grid2(.~model, axes = 'all') +
1238. ggtitle('All shared genes') +
1239. xlab('weight.G2S') +
1240. ylab('weight.GTEx') +
1241. coord_cartesian( clip = "off")+
1242. theme_minimal()+
1243. theme(panel.spacing = unit(1, "cm", data = NULL),
1244. strip.text.x = element_text(margin = margin(5, 5, 12, 5, "pt")))
1247. name='Output/dual_weight_corr_grid_all.pdf'
1248. pdf(file = name, width = 8, height = 2.6, pointsize = 12, bg = "white")
1249. print(cmd_snpcorr)
1250. dev.off()
1252. print(cmd_snpcorr)
1254. ```
1256. ### Effect size correlation
1258. normalized zscore correlation across models
1259. ```{r fig.width = 12, fig.height = 4.25}
1262. duo_snps <- snp_dat_all
1264. ds_gene <- duo_snps[, .(N_g2s_SNPs = sum(!(is.na(weight.G2S)) & shared == F),
1265. N_GTeX_SNPs = sum(!(is.na(weight.GTEx)) & shared == F),
1266. N_shared_SNPs = sum(shared == TRUE),
1267. N_total_snps = length(unique(rsid)),
1268. prop_g2s = sum(!(is.na(weight.G2S)))/length(unique(rsid))), by = 'gene']
1270. #group genes by model presence
1271. ds_gene[, eqtl_source := fcase(N_g2s_SNPs == 0 & N_shared_SNPs == 0, "GTEx",
1272. N_GTeX_SNPs == 0 & N_shared_SNPs == 0, "G2S",
1273. N_shared_SNPs != 0, "shared_overlapping",
1274. N_shared_SNPs == 0 & (N_g2s_SNPs != 0 & N_GTeX_SNPs != 0), "shared_distinct")]
1279. dubs <- setDT(merge(crossmdd, ds_gene, by.x = c('gene.multi_g2S'), by.y = c('gene')))
1282. cmdcor_source <- dubs[, .(zsc_cor = cor.test(zscore.multi_g2S, zscore.eur_gtex)$estimate,
1283. zsc_ciL = cor.test(zscore.multi_g2S, zscore.eur_gtex)$conf.int[1],
1284. zsc_ciH = cor.test(zscore.multi_g2S, zscore.eur_gtex)$conf.int[2],
1285. .N), by = c('training_model.multi_g2S', 'eqtl_source')]
1289. fx_cor <- ggplot(data = cmdcor_source, aes(x = eqtl_source, y = zsc_cor, color = training_model.multi_g2S)) +
1290. geom_hline(yintercept = 0, color = 'red') +
1291. geom_errorbar(aes(ymin = zsc_ciL, ymax = zsc_ciH), width = 0.1, color = 'black') +
1292. geom_point(size = 2) +
1293. geom_label(aes(y = -0.75, label = paste('N =', N)), size = 3) +
1294. #scale_color_manual(values=c('green4', 'slateblue3')) +
1295. ylab('normalized effect size correlation') +
1296. xlab('SNP predictor model specificity') +
1297. ylim(c(-1, 1))+
1298. facet_grid(.~training_model.multi_g2S, scales = 'free') +
1299. theme_bw() +
1300. theme(axis.text.x = element_text(angle = 60, hjust = 1))
1304. name = 'Output/dual_genefx_corplot_grid.pdf'
1305. pdf(file = name, width = 12, height = 4.25, pointsize = 12, bg = "white")
1306. print(fx_cor)
1307. dev.off()
1309. print(fx_cor)
1311. ```
1313. ### Fig 4b: N SNPs per gene
1315. N Snp weights used per gene across the two models.
1317. ```{r fig.width = 12, fig.height = 4.25}
1320. cmdcor_source <- dubs[eqtl_source == 'shared_overlapping' & (bhpval.eur_gtex <= 0.05 | bhpval.multi_g2S <= 0.05), .(zsc_cor = cor.test(zscore.multi_g2S, zscore.eur_gtex)$estimate,
1321. zsc_ciL = cor.test(zscore.multi_g2S, zscore.eur_gtex)$conf.int[1],
1322. zsc_ciH = cor.test(zscore.multi_g2S, zscore.eur_gtex)$conf.int[2],
1323. .N), by = c('training_model.multi_g2S', 'eqtl_source')]
1326. cmd_sigs <- dubs[bhpval.eur_gtex <= 0.05 | bhpval.multi_g2S <= 0.05, .(.N), by = c('training_model.multi_g2S', 'eqtl_source')]
1328. snp_dat_all_L <- setDT(melt(snp_dat_all[, -c('shared')], id.vars = c('gene', 'rsid','model')))
1329. snp_dat_N <- snp_dat_all_L[!is.na(value) & variable == 'weight.G2S', .(N = length(unique(rsid))), by = c('gene', 'model')]
1331. snp_dat_all_L2 <- snp_dat_all_L
1332. snp_dat_all_L2[,gr := paste0(gene, rsid)]
1333. snp_dat_N2 <- snp_dat_all_L2[!is.na(value) & variable == 'weight.GTEx', .(N = length(unique(gr))), by = 'gene']
1334. snp_dat_N2[, model := 'GTEx_Whole_Blood']
1336. snp_dat_Ns <- rbind(snp_dat_N, snp_dat_N2)
1337. # median weight across all SNPs used
1339. eqtl_vplot <- ggplot(data = snp_dat_Ns, aes(x = model, y = N, fill = model)) +
1340. geom_violin(scale = 'width', position=position_dodge(1), drop = F) +
1341. geom_boxplot(width = 0.1, color = 'black', position=position_dodge(1), outlier.size = 0.5) +
1342. xlab('SNP predictor model ')+
1343. ylab('SNP predictors per gene') +
1344. theme_minimal()+
1345. theme(legend.position = 'none') +
1346. theme_bw() +
1347. theme(text = element_text(size = 13), axis.text = element_text(size = 13))
1351. name = 'Output/dual_snpPredNplot_gridvplot.pdf'
1352. pdf(file = name, width = 11, height = 3, pointsize = 12, bg = "white")
1353. print(eqtl_vplot)
1354. dev.off()
1356. print(eqtl_vplot)
1357. ```
1359. ### Fig 4c: Median SNP magnitude
1361. Median magnitude of the used SNP weights across the models.
1363. ```{r fig.width = 12, fig.height = 4.25}
1366. snp_dat_all_L1 <- snp_dat_all_L[!is.na(value) & variable == 'weight.G2S', -c('variable')]
1367. snp_dat_all_L2 <- unique(snp_dat_all_L[!is.na(value) & variable == 'weight.GTEx', -c('variable', 'model')])
1368. snp_dat_all_L2[, model := 'GTEx_Whole_Blood']
1369. snp_dat_wts <- rbind(snp_dat_all_L1, snp_dat_all_L2)
1372. eqtl_wtplot <- ggplot(data = snp_dat_wts, aes(x = model, y = abs(value), fill = model)) +
1373. geom_boxplot(color = 'black', position=position_dodge(1), outlier.shape = NA) +
1374. xlab('SNP predictor model ')+
1375. ylab('median SNP weight magnitude') +
1376. coord_cartesian(ylim= c(0,0.125)) +
1377. theme_minimal()+
1378. theme(legend.position = 'none') +
1379. theme_bw() +
1380. theme(text = element_text(size = 13), axis.text = element_text(size = 13))
1384. name = 'Output/dual_snpPred_wt_plot_gridvplot.pdf'
1385. pdf(file = name, width = 11, height = 3, pointsize = 12, bg = "white")
1386. print(eqtl_wtplot)
1387. dev.off()
1389. print(eqtl_wtplot)
1392. ```
1394. ### Fig 4d: All correlation Comparison
1396. Compile into single figure
1398. ```{r fig.width = 12, fig.height = 4.25}
1400. crossmdd2 <- crossmdd
1402. crossmdd2[bhpval.eur_gtex <= 0.05 | bhpval.multi_g2S <= 0.05, sig := 'sig']
1403. crossmdd2[!(bhpval.eur_gtex <= 0.05 | bhpval.multi_g2S <= 0.05), sig := 'nonsig']
1406. cmd_snpcor <- snp_dat_all[, .(wt_cor = cor.test(weight.G2S, weight.GTEx)$estimate,
1407. zsc_ciL = cor.test(weight.G2S, weight.GTEx)$conf.int[1],
1408. zsc_ciH = cor.test(weight.G2S, weight.GTEx)$conf.int[2],
1409. .N), by = c('model')]
1412. cmdcor <- crossmdd[, .(zsc_cor = cor.test(zscore.multi_g2S, zscore.eur_gtex)$estimate,
1413. zsc_ciL = cor.test(zscore.multi_g2S, zscore.eur_gtex)$conf.int[1],
1414. zsc_ciH = cor.test(zscore.multi_g2S, zscore.eur_gtex)$conf.int[2],
1415. .N), by = c('training_model.multi_g2S')]
1416. cmdcorsig <- crossmdd[bhpval.eur_gtex <= 0.05 & bhpval.multi_g2S <= 0.05, .(zsc_cor = cor.test(zscore.multi_g2S, zscore.eur_gtex)$estimate,
1417. zsc_ciL = cor.test(zscore.multi_g2S, zscore.eur_gtex)$conf.int[1],
1418. zsc_ciH = cor.test(zscore.multi_g2S, zscore.eur_gtex)$conf.int[2],
1419. .N), by = c('training_model.multi_g2S')]
1420. cmdcorsig2 <- crossmdd[bhpval.eur_gtex <= 0.05 | bhpval.multi_g2S <= 0.05, .(zsc_cor = cor.test(zscore.multi_g2S, zscore.eur_gtex)$estimate,
1421. zsc_ciL = cor.test(zscore.multi_g2S, zscore.eur_gtex)$conf.int[1],
1422. zsc_ciH = cor.test(zscore.multi_g2S, zscore.eur_gtex)$conf.int[2],
1423. .N), by = c('training_model.multi_g2S')]
1426. snp1 <- setnames(cmd_snpcor, c('wt_cor', 'zsc_ciL', 'zsc_ciH'), c('cor', 'ciL', 'ciH'))
1427. #snp1 <- cmd_snpcor
1428. snp1$group <- 'SNP weights'
1430. gene1 <- setnames(cmdcor, c('training_model.multi_g2S', 'zsc_cor', 'zsc_ciL', 'zsc_ciH'), c('model','cor', 'ciL', 'ciH'))
1431. #gene1 <- cmdcor
1432. gene1$group <- 'all genes'
1434. genesig <- setnames(cmdcorsig2, c('training_model.multi_g2S', 'zsc_cor', 'zsc_ciL', 'zsc_ciH'), c('model','cor', 'ciL', 'ciH'))
1435. #genesig <- cmdcorsig2
1436. genesig$group <- 'sig genes'
1438. compdata <- setDT(rbind(snp1, gene1, genesig))
1440. compdata$group <- factor(compdata$group, levels=c('SNP weights', 'all genes', 'sig genes'))
1443. allcor_plot <- ggplot(data = compdata, aes(x = group, y = cor, color = model, group = model)) +
1444. geom_errorbar(aes(ymin = ciL, ymax = ciH), width = 0.3, color = 'black', position=position_dodge(width=0.9)) +
1445. geom_point( position=position_dodge(width = 0.9), size = 3) +
1446. xlab('measurement') +
1447. ylab('r correlation coefficient') +
1448. ylim(c(0,1)) +
1449. theme_bw() +
1450. theme(text = element_text(size = 13), axis.text = element_text(size = 13))
1454. name = 'Output/dual_allcordat_mddcomp.pdf'
1455. pdf(file = name, width = 11, height = 3, pointsize = 12, bg = "white")
1456. print(allcor_plot)
1457. dev.off()
1459. print(allcor_plot)
1462. ```
1465. ## NeuroimaGene analysis
1467. Use the NeuroimaGene resource to get the neuroimaging features associated with the
1468. TWAS results for each disease.
1470. ### Fig 6a: Neuroimaging visualizations
1472. ```{r neuroimaGene1}
1473. library(neuroimaGene)
1475. ########## NeuroimaGene analysis
1477. all <- dt
1479. bdg <- unique(all[gwas_phenotype == 'PGC_BD1_2021',]$gene)
1480. szg <- unique(all[gwas_phenotype == 'PGC_all_SCZ_2022']$gene)
1481. adg <- unique(all[gwas_phenotype == 'PGC_allADHD_2022']$gene)
1482. mdg <- unique(all[gwas_phenotype == 'PGC_allMDD_2023']$gene)
1483. ptg <- unique(all[gwas_phenotype == 'PGC_ALL_PTSD_2024']$gene)
1484. aug <- unique(all[gwas_phenotype == 'SUD_alc_2019']$gene)
1487. bdng <- neuroimaGene(bdg)
1488. neuro_vis(bdng)
1489. ```
1493. ```{r}
1494. szng <- neuroimaGene(szg)
1495. neuro_vis(szng)
1496. ```
1501. ```{r}
1502. adng <- neuroimaGene(adg)
1503. neuro_vis(adng)
1504. ```
1507. ```{r}
1508. mdng <- neuroimaGene(mdg)
1509. neuro_vis(mdng)
1510. ```
1513. ```{r}
1514. ptng <- neuroimaGene(ptg)
1515. neuro_vis(ptng)
1516. ```
1520. ```{r}
1521. aung <- neuroimaGene(aug)
1522. neuro_vis(aung)
1523. ```
1526. ```{r}
1527. bdng2 <- neuroimaGene(bdg, atlas = 'aseg_volume')
1528. neuro_vis(bdng2, atlas = 'Subcortex')
1529. ```
1532. ```{r}
1533. szng2 <- neuroimaGene(szg, atlas = 'aseg_volume')
1534. neuro_vis(szng2, atlas = 'Subcortex')
1535. ```
1538. ```{r}
1539. adng2 <- neuroimaGene(adg, atlas = 'aseg_volume')
1540. neuro_vis(adng2, atlas = 'Subcortex')
1541. ```
1543. ```{r}
1544. mdng2 <- neuroimaGene(mdg, atlas = 'aseg_volume')
1545. neuro_vis(mdng2, atlas = 'Subcortex')
1546. ```
1550. ```{r}
1551. ptng2 <- neuroimaGene(ptg, atlas = 'aseg_volume')
1552. neuro_vis(ptng2, atlas = 'Subcortex')
1553. ```
1555. ```{r}
1556. aung2 <- neuroimaGene(aug, atlas = 'aseg_volume')
1557. neuro_vis(aung2, atlas = 'Subcortex')
1558. ```
1560. Compile data
1562. ```{r}
1563. adng$group <- 'ADHD'
1564. adng2$group <- 'ADHD'
1565. aung2$group <- 'AUD'
1566. aung$group <- 'AUD'
1567. bdng$group <- 'BD1'
1568. bdng2$group <- 'BD1'
1569. mdng$group <- 'MDD'
1570. mdng2$group <- 'MDD'
1571. ptng$group <- 'PTSD'
1572. ptng2$group <- 'PTSD'
1573. szng$group <- 'SCZ'
1574. szng2$group <- 'SCZ'
1575. ng_data_ <- rbind(bdng, szng, adng, mdng, ptng, aung, bdng2, szng2, adng2, mdng2, ptng2, aung2)
1577. ng_data <- setDT(merge(ng_data_, neuroimaGene::anno, by = 'gwas_phenotype'))
1578. ```
1581. Comparison with no G2S
1583. ```{r}
1585. eur <- dt[training_model == 'PrediXcan_Whole_Blood']
1587. bdg_e <- unique(eur[gwas_phenotype == 'PGC_BD1_2021',]$gene)
1588. szg_e <- unique(eur[gwas_phenotype == 'PGC_all_SCZ_2022']$gene)
1589. adg_e <- unique(eur[gwas_phenotype == 'PGC_allADHD_2022']$gene)
1590. mdg_e <- unique(eur[gwas_phenotype == 'PGC_allMDD_2023']$gene)
1591. ptg_e <- unique(eur[gwas_phenotype == 'PGC_ALL_PTSD_2024']$gene)
1592. aug_e <- unique(eur[gwas_phenotype == 'SUD_alc_2019']$gene)
1595. bdng_e <- neuroimaGene(bdg_e)
1596. szng_e <- neuroimaGene(szg_e)
1597. adng_e <- neuroimaGene(adg_e)
1598. mdng_e <- neuroimaGene(mdg_e)
1599. ptng_e <- neuroimaGene(ptg_e)
1600. aung_e <- neuroimaGene(aug_e)
1601. bdng2_e <- neuroimaGene(bdg_e, atlas = 'aseg_volume')
1602. szng2_e <- neuroimaGene(szg_e, atlas = 'aseg_volume')
1603. adng2_e <- neuroimaGene(adg_e, atlas = 'aseg_volume')
1604. mdng2_e <- neuroimaGene(mdg_e, atlas = 'aseg_volume')
1605. ptng2_e <- neuroimaGene(ptg_e, atlas = 'aseg_volume')
1606. aung2_e <- neuroimaGene(aug_e, atlas = 'aseg_volume')
1609. adng_e$group <- 'ADHD'
1610. adng2_e$group <- 'ADHD'
1611. aung2_e$group <- 'AUD'
1612. aung_e$group <- 'AUD'
1613. bdng_e$group <- 'BD1'
1614. bdng2_e$group <- 'BD1'
1615. mdng_e$group <- 'MDD'
1616. mdng2_e$group <- 'MDD'
1617. ptng_e$group <- 'PTSD'
1618. ptng2_e$group <- 'PTSD'
1619. szng_e$group <- 'SCZ'
1620. szng2_e$group <- 'SCZ'
1621. ng_data_e_ <- rbind(bdng_e, szng_e, adng_e, mdng_e, ptng_e, aung_e, bdng2_e, szng2_e, adng2_e, mdng2_e, ptng2_e, aung2_e)
1623. ng_data_e <- setDT(merge(ng_data_e_, neuroimaGene::anno, by = 'gwas_phenotype'))[, -c('zscore', 'atl_BHpval', 'fMRI_node_1', 'fMRI_node_2')]
1624. ```
1627. now g2s
1629. ```{r}
1630. g2s <- dt[training_model != 'PrediXcan_Whole_Blood']
1632. bdg_g <- unique(g2s[gwas_phenotype == 'PGC_BD1_2021',]$gene)
1633. szg_g <- unique(g2s[gwas_phenotype == 'PGC_all_SCZ_2022']$gene)
1634. adg_g <- unique(g2s[gwas_phenotype == 'PGC_allADHD_2022']$gene)
1635. mdg_g <- unique(g2s[gwas_phenotype == 'PGC_allMDD_2023']$gene)
1636. ptg_g <- unique(g2s[gwas_phenotype == 'PGC_ALL_PTSD_2024']$gene)
1637. aug_g <- unique(g2s[gwas_phenotype == 'SUD_alc_2019']$gene)
1640. bdng_g <- neuroimaGene(bdg_g)
1641. szng_g <- neuroimaGene(szg_g)
1642. adng_g <- neuroimaGene(adg_g)
1643. mdng_g <- neuroimaGene(mdg_g)
1644. ptng_g <- neuroimaGene(ptg_g)
1645. aung_g <- neuroimaGene(aug_g)
1646. bdng2_g <- neuroimaGene(bdg_g, atlas = 'aseg_volume')
1647. szng2_g <- neuroimaGene(szg_g, atlas = 'aseg_volume')
1648. adng2_g <- neuroimaGene(adg_g, atlas = 'aseg_volume')
1649. mdng2_g <- neuroimaGene(mdg_g, atlas = 'aseg_volume')
1650. ptng2_g <- neuroimaGene(ptg_g, atlas = 'aseg_volume')
1651. aung2_g <- neuroimaGene(aug_g, atlas = 'aseg_volume')
1654. adng_g$group <- 'ADHD'
1655. adng2_g$group <- 'ADHD'
1656. aung2_g$group <- 'AUD'
1657. aung_g$group <- 'AUD'
1658. bdng_g$group <- 'BD1'
1659. bdng2_g$group <- 'BD1'
1660. mdng_g$group <- 'MDD'
1661. mdng2_g$group <- 'MDD'
1662. ptng_g$group <- 'PTSD'
1663. ptng2_g$group <- 'PTSD'
1664. szng_g$group <- 'SCZ'
1665. szng2_g$group <- 'SCZ'
1666. ng_data_g_ <- rbind(bdng_g, szng_g, adng_g, mdng_g, ptng_g, aung_g, bdng2_g, szng2_g, adng2_g, mdng2_g, ptng2_g, aung2_g)
1668. ng_data_g <- setDT(merge(ng_data_g_, neuroimaGene::anno, by = 'gwas_phenotype'))[, -c('zscore', 'atl_BHpval', 'fMRI_node_1', 'fMRI_node_2')]
1669. ```
1671. Separate out and print the neuroimagene findings specific to G2S models
1673. ```{r}
1674. ng_data_g$id <- paste0(ng_data_g$gwas_phenotype, ng_data_g$group)
1675. ng_data_e$id <- paste0(ng_data_e$gwas_phenotype, ng_data_e$group)
1677. ng_data_gonly <- ng_data_g[!(id %in% ng_data_e$id),][measurement != 'QC',]
1679. ng_data_gonly[, xlabel := paste0(atlas, '\n', measurement)]
1681. ng_plotdata <- ng_data_gonly[, .(N = length(unique(gwas_phenotype))), by = c('group', 'xlabel')]
1683. ng_plotdata$group <- factor(ng_plotdata$group, levels = c("PTSD", "SCZ", "ADHD", "MDD", "BD1", "AUD"))
1686. ng_barplot <- ggplot(ng_plotdata, aes(y = N, x = xlabel, fill = xlabel)) +
1687. geom_bar(stat = 'identity') +
1688. facet_grid(group~.) +
1689. scale_fill_brewer(palette = 'Dark2') +
1690. ylab('G2S specific NIDPs') +
1691. xlab('measurement') +
1692. theme_minimal() +
1693. theme(axis.text.x = element_text(angle = 270), legend.position = 'none')
1695. name = 'Output/ng_g2s_barplot.pdf'
1696. pdf(file = name, width = 2.5, height = 7, pointsize = 12, bg = "white")
1697. print(ng_barplot)
1698. dev.off()
1700. print(ng_barplot)
1701. ```
1703. ### Fig S1: Descriptive Summary Data
1705. Show the number of results for each model/disease pair and split by shared vs individually significant.
1707. **Special thanks to Nathan Watkins and Tavian Bowen-Moore**
1709. ```{r, echo=FALSE}
1711. library(data.table)
1712. library(ggplot2)
1713. library(gridExtra)
1716. dt[training_model == 'All_Whole_Blood', `gene expression model` := 'G2S All WB']
1717. dt[training_model == 'AA_Whole_Blood', `gene expression model` := 'G2S African American WB']
1718. dt[training_model == 'MX_Whole_Blood', `gene expression model` := 'G2S Mexican American WB']
1719. dt[training_model == 'PR_Whole_Blood', `gene expression model` := 'G2S Puerto Rican WB']
1720. dt[training_model == 'PrediXcan_Whole_Blood', `gene expression model` := 'GTEx European WB']
1722. dt1 <- dt[!training_model == 'PrediXcan_Whole_Blood']
1723. dt2 <- dt[training_model == 'PrediXcan_Whole_Blood']
1725. # Function to clean gene names
1726. clean_gene_name <- function(gene_name) {
1727. sub("\\..*", "", gene_name)
1728. }
1730. # Function to create and store plots
1731. create_plots <- function(MAtwas, MUtwas) {
1732. # Clean gene names in both datasets
1733. MAtwas[, gene_name := clean_gene_name(gene_name)]
1734. MUtwas[, gene_name := clean_gene_name(gene_name)]
1736. # List to store plots
1737. plots <- list()
1739. for (blood_type in unique(MUtwas$training_model)) {
1741. # Loop through each unique blood type in MUtwas
1742. for (phenotype in unique(MAtwas$gwas_phenotype)) {
1744. # Subset the data based on the current phenotype and blood type
1745. tut <- MUtwas[training_model == blood_type & gwas_phenotype == phenotype,]
1746. tot <- MAtwas[training_model == 'PrediXcan_Whole_Blood' & gwas_phenotype == phenotype,]
1748. # Check if all entries in tut are 'GTEX_whole_blood'
1749. if (all(tut$training_model == 'PrediXcan_Whole_Blood')) {
1750. next
1751. }
1753. # Split the tut data.table into 2 data.tables according to training model
1754. aa <- tut[training_model == blood_type,]
1755. eur <- tot
1758. # Get vectors of all gene associations with the disease in each set
1759. aa_gns <- unique(aa$gene_name)
1760. eur_gns <- unique(eur$gene_name)
1762. # Create vectors of gene names that are unique to AA, GTEx, and shared
1763. unique_eur <- setdiff(eur_gns, aa_gns)
1764. unique_aa <- setdiff(aa_gns, eur_gns)
1765. shared <- intersect(aa_gns, eur_gns)
1767. # Custom names and color mapping for plotting
1768. custom_names <- c('AA_Whole_Blood' = 'AA',
1769. 'All_Whole_Blood' = 'ALL',
1770. 'MX_Whole_Blood' = 'MX',
1771. 'PR_Whole_Blood' = 'PR',
1772. 'PrediXcan_Whole_Blood' = 'GTEx')
1774. color_mapping <- c('AA' = 'orange',
1775. 'ALL' = 'red',
1776. 'MX' = 'pink',
1777. 'PR' = 'green',
1778. 'GTEx' = 'blue',
1779. 'shared' = 'cyan')
1781. dz_names <- c('PGC_allADHD_2022' = 'ADHD',
1782. 'PGC_all_SCZ_2022' = 'SCZ',
1783. 'PGC_allMDD_2023' = 'MDD',
1784. 'PGC_BD1_2021' = 'BD1',
1785. 'SUD_alc_2019' = 'SUD-A',
1786. 'PGC_ALL_PTSD_2024' = 'PTSD')
1788. # Set the disease name
1789. dz_name <- dz_names[unique(aa$gwas_phenotype)]
1791. # Create data.table with data for plotting
1792. gene_cts <- data.table(study_specificity = c('GTEx', custom_names[blood_type], 'shared'),
1793. N = c(length(unique_eur), length(unique_aa), length(shared)))
1795. # Set the factor levels for study_specificity
1796. gene_cts$study_specificity <- factor(gene_cts$study_specificity, levels = names(color_mapping))
1798. # Create and plot the ggplot object
1799. p2 <- ggplot(data = gene_cts,
1800. aes(x = study_specificity, y = N, fill = study_specificity)) +
1801. geom_bar(stat = 'identity') +
1802. ggtitle(dz_name) +
1803. xlab("Study Specificity") +
1804. ylab("Number of Genes") +
1805. scale_fill_manual(values = color_mapping) +
1806. theme_minimal() +
1807. theme(
1808. plot.title = element_text(size = 10), # Reduce title size
1809. axis.title = element_text(size = 9), # Reduce axis titles size
1810. axis.text = element_text(size = 9), # Reduce axis text size
1811. plot.margin = margin(10, 10, 10, 10),
1812. legend.position = 'none'
1813. )
1815. # Store the plot in the list
1816. plots[[paste(phenotype, blood_type, sep = "_")]] <- p2
1817. }
1818. }
1820. return(plots)
1821. }
1823. # Generate and save plots to PDF
1824. save_plots_to_pdf <- function(plots, file_name) {
1825. pdf(file_name, width = 13, height = 8)
1826. do.call(grid.arrange, c(plots, ncol = 6, padding = unit(1, "lines")))
1827. dev.off()
1828. }
1830. # usage
1831. plots <- create_plots(dt2, dt1)
1832. save_plots_to_pdf(plots, "Output/gene_unique_comp2.pdf")
1834. knitr::include_graphics("Output/gene_unique_comp2.pdf")
1836. ```
1839. # Corrplots
1841. ### Fig 5a: correlation plots
1842. Get correlation of disase TWAS results across models to see if the models are predicting the same
1843. inter-disease transcriptomic relationships.
1845. ```{r}
1848. library(data.table)
1849. library(corrplot)
1851. dt_ <- fread('Data/all_exc_assocs_PGCancestryTWAS.txt', header = TRUE)
1853. # exclude replication/comparison data from initial analysis.
1854. gwas <- c('PGC_AA_MDD_2023', 'PGC_AA_PTSD_2024','PGC_AA_SCZ_2022','PGC_ASIAN_SCZ_2022','PGC_EUR_PTSD_2024','PGC_HNA_PTSD_2024', 'PGC_EURO_SCZ_2022', 'PGC_LATINO_SCZ_2022')
1855. mod <- c('PrediXcan_Brain_Cortex')
1857. dt <- dt_[!(gwas_phenotype %in% gwas | training_model %in% mod),]
1860. dt[training_model == 'All_Whole_Blood', `gene expression model` := 'G2S All WB']
1861. dt[training_model == 'AA_Whole_Blood', `gene expression model` := 'G2S African American WB']
1862. dt[training_model == 'MX_Whole_Blood', `gene expression model` := 'G2S Mexican American WB']
1863. dt[training_model == 'PR_Whole_Blood', `gene expression model` := 'G2S Puerto Rican WB']
1864. dt[training_model == 'PrediXcan_Whole_Blood', `gene expression model` := 'GTEx European WB']
1866. dt1a <- dt[!training_model == 'PrediXcan_Whole_Blood']
1867. dt2a <- dt[training_model == 'PrediXcan_Whole_Blood']
1869. # Function to clean gene names
1870. clean_gene_name <- function(gene_name) {
1871. sub("\\..*", "", gene_name)
1872. }
1874. # Function to create and save correlation plots
1875. create_and_save_corr_plots <- function(MAtwas, MUtwas) {
1876. # Clean gene names in both datasets
1877. MAtwas[, gene_name := clean_gene_name(gene_name)]
1878. MUtwas[, gene_name := clean_gene_name(gene_name)]
1880. # Get unique training models from MUtwas and add GTEX from MAtwas
1881. training_models <- unique(MUtwas$training_model)
1882. training_models <- c(training_models, 'PrediXcan_Whole_Blood')
1884. # Define a custom color palette
1885. custom_colors <- colorRampPalette(c("red", "grey60", "blue"))(1000)
1887. # Open a PDF device with Letter dimensions
1888. pdf("Output/correlograms.pdf", width = 8.5, height = 11)
1890. for (model in training_models) {
1891. modl <- copy(model)
1892. # Filter data for the current training model
1893. if (model == 'PrediXcan_Whole_Blood') {
1894. filtered_data <- MAtwas[training_model == 'PrediXcan_Whole_Blood', c('gwas_phenotype', 'gene', 'zscore')]
1895. } else {
1896. filtered_data <- MUtwas[training_model == modl, c('gwas_phenotype', 'gene', 'zscore')]
1897. }
1899. unique_genes <- unique(filtered_data$gene)
1900. unique_phenotypes <- unique(filtered_data$gwas_phenotype)
1902. matrix_data <- matrix(NA, nrow = length(unique_genes), ncol = length(unique_phenotypes))
1903. rownames(matrix_data) <- unique_genes
1904. colnames(matrix_data) <- unique_phenotypes
1906. for (i in 1:nrow(filtered_data)) {
1907. row <- filtered_data[i, ]
1908. matrix_data[row$gene, row$gwas_phenotype] <- row$zscore
1909. }
1911. cor_matrix <- cor(matrix_data, use = "pairwise.complete.obs")
1913. # Adjust plot margins to accommodate longer titles
1914. par(mar = c(1, 1, 5, 1))
1916. corrplot(corr = cor_matrix,
1917. method = "circle",
1918. type = "lower",
1919. tl.pos = "tl",
1920. order = "original",
1921. col = custom_colors, # Use custom colors
1922. title = paste("Correlation plot for training model:", model),
1923. mar = c(0, 0, 2, 0)) # Adjust margins to fit the title
1924. }
1926. # Close the PDF device
1927. dev.off()
1928. }
1930. # Example usage
1931. create_and_save_corr_plots(dt2a, dt1a)
1934. ```
1937. ## Cell type analysis
1939. What are the cell types most implicated in the brain analyses based on the TWAS genes that emerged as significant.
1941. ```{r}
1942. pg <- fread('../../PanglaoDB_markers_27_Mar_2020.tsv', header = T)
1945. overlap <- merge(
1946. dt[, .(`gene_name`, gwas_phenotype)],
1947. pg[, .(`official gene symbol`, `cell type`)],
1948. by.x = "gene_name",
1949. by.y = "official gene symbol",
1950. allow.cartesian = TRUE
1951. )
1953. # Count overlaps for each disease × cell_type pair
1954. result <- overlap[, .(n_genes = uniqueN(gene_name), genes = paste(unique(gene_name), collapse = '|')), by = c('gwas_phenotype', 'cell type')]
1956. # Optional: ensure long format (already long here)
1957. setorder(result, gwas_phenotype, `n_genes`)
1960. write.table(result[n_genes > 1, ], 'Output/disease_cell_types.txt', col.names = T, row.names = F, quote = F, sep = '\t')
1961. ```
1963. ## TWAS Power from GWAS ancestries
1965. ### Fig S24: Ancestry specific TWAS
1966. We hypothesized that the derivation of statistically significant twas results from different GWAS was affected by the N for each ancestry specific GWAS. We next examine how similar the TWAS results from the GWAS subsets are to the general data across multiple significance cutoffs.
1968. ```{r}
1969. allnom <- fread('Data/all_exc_assocs_PGCancestryTWAS.txt', header = T)
1970. allnom <- allnom[pvalue <= 0.05,]
1974. # exclude replication/comparison data from initial analysis.
1975. gwas <- c('PGC_AA_MDD_2023', 'PGC_AA_PTSD_2024','PGC_AA_SCZ_2022','PGC_ASIAN_SCZ_2022','PGC_EUR_PTSD_2024','PGC_HNA_PTSD_2024')
1976. mod <- c('PrediXcan_Brain_Cortex')
1978. anc <- allnom[gwas_phenotype %in% gwas & !(training_model %in% mod),]
1981. gpdict1 <- data.table(gwas_phenotype = gwas, diagnosis = c('MDD', 'PTSD', 'SCZ', 'SCZ', 'PTSD', 'PTSD'), gwas_ancestry = c('AA', 'AA', 'AA', 'EA', 'EUR', 'HNA'))
1982. gpdict2 <- data.table(gwas_phenotype = c('PGC_allADHD_2022', 'PGC_allMDD_2023', 'PGC_ALL_PTSD_2024', 'PGC_all_SCZ_2022', 'PGC_BD1_2021', 'SUD_alc_2019'),
1983. diagnosis = c('ADHD', 'MDD', 'PTSD', 'SCZ', 'BD1', 'AUD'),
1984. gwas_ancestry = c('multi', 'multi', 'multi', 'multi', 'EUR', 'multi'))
1991. dt_ <- fread('Data/all_exc_BHsig_assocs_PGCancestryTWAS.txt', header = TRUE)
1993. dt <- dt_[!(gwas_phenotype %in% gwas | training_model %in% mod),]
1995. dt[training_model == 'All_Whole_Blood', `gene expression model` := 'G2S All WB']
1996. dt[training_model == 'AA_Whole_Blood', `gene expression model` := 'G2S African American WB']
1997. dt[training_model == 'MX_Whole_Blood', `gene expression model` := 'G2S Mexican American WB']
1998. dt[training_model == 'PR_Whole_Blood', `gene expression model` := 'G2S Puerto Rican WB']
1999. dt[training_model == 'PrediXcan_Whole_Blood', `gene expression model` := 'GTEx European WB']
2003. anc2 <- setDT(merge(anc, gpdict1, by = 'gwas_phenotype'))
2004. dt2 <- setDT(merge(dt, gpdict2, by = c('gwas_phenotype')))
2006. duo <- setDT(merge(dt2, anc2, by = c('gene', 'diagnosis', 'training_model'), suffixes = c('.multi', '.sa')))
2011. pdata <- data.table()
2012. for (dz in unique(duo$diagnosis)) {
2013. for (g_anc in unique(duo$gwas_ancestry.sa)) {
2014. for(pthresh in c(5*10^seq(-2, -5, -1), 1*10^seq(-2, -5, -1))) {
2015. if(nrow(duo[pvalue.sa <= pthresh & (diagnosis == dz & gwas_ancestry.sa == g_anc),]) >= 5) {
2016. mod <- lm(data = duo[pvalue.sa <= pthresh & (diagnosis == dz & gwas_ancestry.sa == g_anc),], zscore.multi~zscore.sa)
2017. newdat <- data.table(pvalue = pthresh, diagnosis = dz, gwas_ancestry = g_anc, N = nrow(duo[pvalue.sa <= pthresh & (diagnosis == dz & gwas_ancestry.sa == g_anc),]), r2 = summary(mod)$adj.r.squared, regression_p = summary(mod)$coefficients[2,4])
2018. pdata <- rbind(pdata, newdat)
2019. }
2020. }
2021. }
2022. }
2025. pz <- ggplot(pdata, aes(x = -log10(pvalue), y = r2, color = -log10(regression_p))) +
2026. geom_hline(yintercept = 0) +
2027. geom_point() +
2028. facet_wrap(diagnosis~gwas_ancestry, scales = 'free', nrow = 6) +
2029. scale_y_continuous(limits = c(-0.5, 1)) +
2030. scale_x_continuous(limits = c(0, 5)) +
2031. theme_minimal() +
2032. theme(legend.position = 'bottom')
2036. name = 'Output/anc_specific_r2.pdf'
2037. pdf(file = name, width = 3.5, height = 10, pointsize = 12, bg = "white")
2038. print(pz)
2039. dev.off()
2041. print(pz)
2042. ```
2043. ### Fig S25: produce TWAS plots for multi-ancestry GWAS
2045. ```{r}
2046. library(data.table)
2047. library(ggplot2)
2048. library(ggrepel)
2049. library(stringr)
2050. library(colorspace)
2052. dt_ <- fread('Data/all_exc_assocs_PGCancestryTWAS.txt', header = TRUE)
2054. # include replication/comparison data from initial analysis.
2055. gwas <- c('PGC_AA_MDD_2023', 'PGC_AA_PTSD_2024','PGC_AA_SCZ_2022','PGC_ASIAN_SCZ_2022','PGC_EUR_PTSD_2024','PGC_HNA_PTSD_2024', 'PGC_EURO_SCZ_2022','PGC_LATINO_SCZ_2022')
2056. mod <- c('PrediXcan_Brain_Cortex')
2058. dt <- dt_[gwas_phenotype %in% gwas & !(training_model %in% mod),]
2061. dt[training_model == 'All_Whole_Blood', `gene expression model` := 'G2S All WB']
2062. dt[training_model == 'AA_Whole_Blood', `gene expression model` := 'G2S African American WB']
2063. dt[training_model == 'MX_Whole_Blood', `gene expression model` := 'G2S Mexican American WB']
2064. dt[training_model == 'PR_Whole_Blood', `gene expression model` := 'G2S Puerto Rican WB']
2065. dt[training_model == 'PrediXcan_Whole_Blood', `gene expression model` := 'GTEx European WB']
2068. tms <- rev(c('G2S All WB', 'G2S African American WB', 'G2S Mexican American WB', 'G2S Puerto Rican WB', 'GTEx European WB'))
2070. dt$`gene expression model` <- factor(dt$`gene expression model`, levels = tms)
2073. dt$gene <- str_remove_all(dt$gene, '\\..+')
2075. # Iterate through the tissue models and NIDPs from the TWAS results
2078. # Get the top 15 p-values for labeling purposes
2079. top15 <- dt[, .(pvalue = min(pvalue)), by = c('gene','gene_name')][order(pvalue)][1:8, c('gene','gene_name','pvalue')]
2080. top15x = setDT(merge(top15, all_genes_loc, by.x = 'gene', by.y = 'ensembl_gene_id'))
2081. top15x$short_name = 'top_15'
2083. top25 <- dt[, .(pvalue = min(pvalue)), by = c('gene','gene_name')][order(pvalue)][1:25, c('gene','gene_name','pvalue')]
2084. top25x = setDT(merge(top25, all_genes_loc, by.x = 'gene', by.y = 'ensembl_gene_id'))
2085. top25x$short_name = 'top_25'
2087. twas2 <- merge(dt, all_genes_loc, by.x = 'gene', by.y = 'ensembl_gene_id')
2089. twas2$`gene expression model` <- factor(twas2$`gene expression model`, levels = tms)
2092. fig2 = ggplot(twas2, aes(x=gn_start, y=-log10(pvalue))) +
2093. # Show all points
2094. geom_point( aes(color=`gene expression model`), alpha=0.8, size=1.5) +
2095. geom_hline(yintercept=-log10(0.05/(dim(adtwas2)[1])), color = "red", size=0.5) +
2096. facet_wrap(vars(gwas_phenotype), scales = 'free', nrow = 6) +
2097. xlab('chromosome') +
2098. ylab('-log 10 pvalue')+
2099. #scale_color_manual(values = rep(c("grey", "skyblue"), 22 )) +
2100. scale_x_continuous(label = chrom_sizes[c(1:18,20,22),]$CHR, breaks= chrom_sizes[c(1:18,20,22),]$chr_label_loc ) +
2101. scale_color_discrete_diverging(palette = 'Berlin') +
2102. guides(color = guide_legend(nrow=3, byrow=FALSE)) +
2103. theme_minimal()+
2104. theme(
2105. legend.position="bottom",
2106. panel.border = element_blank(),
2107. panel.grid.major.x = element_blank(),
2108. panel.grid.minor.x = element_blank(),
2109. axis.title.x = element_blank(),
2110. plot.title = element_text(hjust = 0.5),
2111. legend.text=element_text()
2112. )
2114. name = 'Output/all_assocs_nashplot_pgc_ancs2.pdf'
2115. pdf(file = name, width = 14, height = 15, pointsize = 12, bg = "white")
2116. print(fig2)
2117. dev.off()
2122. ```
2124. ```{r}
2125. print(fig2)
2126. ```
2128. # Revisions Round 1
2130. ```{r}
2131. dir.create('Output/revisions')
2132. ```
2133. ### Fig S2: Statistical contributors to model power
2135. ```{r}
2136. # Do power comparison for N individuals in training set vs number of associations.
2138. # Read in all significant data and subset GWAS findings of interest
2140. sig <- fread('Data/all_exc_BHsig_assocs_PGCancestryTWAS.txt', header = T)
2143. dt <- sig[!(gwas_phenotype %in% c('PGC_AA_MDD_2023', 'PGC_AA_PTSD_2024','PGC_AA_SCZ_2022','PGC_ASIAN_SCZ_2022','PGC_EUR_PTSD_2024','PGC_HNA_PTSD_2024', 'PGC_EURO_SCZ_2022')),]
2146. dt[training_model == 'All_Whole_Blood', `gene expression model` := 'G2S All WB']
2147. dt[training_model == 'AA_Whole_Blood', `gene expression model` := 'G2S African American WB']
2148. dt[training_model == 'MX_Whole_Blood', `gene expression model` := 'G2S Mexican American WB']
2149. dt[training_model == 'PR_Whole_Blood', `gene expression model` := 'G2S Puerto Rican WB']
2150. dt[training_model == 'PrediXcan_Whole_Blood', `gene expression model` := 'GTEx European WB']
2151. dt[training_model == 'PrediXcan_Brain_Cortex', `gene expression model` := 'GTEx European Brain']
2155. dtn <- dt[, .(N_assocs =.N), by = c('gene expression model')]
2159. dtn[ `gene expression model` == 'G2S All WB', N_samples := 2733]
2160. dtn[`gene expression model` == 'G2S African American WB', N_samples := 757]
2161. dtn[`gene expression model` == 'G2S Mexican American WB', N_samples := 784]
2162. dtn[ `gene expression model` == 'G2S Puerto Rican WB', N_samples := 893]
2163. dtn[ `gene expression model` == 'GTEx European WB', N_samples := 670]
2164. dtn[ `gene expression model` == 'GTEx European Brain', N_samples := 205]
2166. dtn[ `gene expression model` == 'G2S All WB', Mean_Afr := 0.32]
2167. dtn[`gene expression model` == 'G2S African American WB', Mean_Afr := 0.8]
2168. dtn[`gene expression model` == 'G2S Mexican American WB', Mean_Afr := 0.04]
2169. dtn[ `gene expression model` == 'G2S Puerto Rican WB', Mean_Afr := 0.22]
2172. dtn[ `gene expression model` == 'G2S All WB', Mean_IAM := 0.24]
2173. dtn[`gene expression model` == 'G2S African American WB', Mean_IAM := 0.01]
2174. dtn[`gene expression model` == 'G2S Mexican American WB', Mean_IAM := 0.57]
2175. dtn[ `gene expression model` == 'G2S Puerto Rican WB', Mean_IAM := 0.10]
2177. library(ggrepel)
2179. lmod <- lm(data = dtn, N_assocs ~ N_samples)
2182. lmod <- lm(data = dtn, N_assocs ~ N_samples)
2183. lmodall <- lm(data = dtn[c(1,2,3,6),], N_assocs ~ Mean_Afr)
2184. summary(lmodall)
2186. lmod_g2s <- lm(data = dtn[!(`gene expression model` %in% c('GTEx European WB', 'GTEx European Brain')),], N_assocs ~ N_samples)
2188. dtn_plot_g2s <- ggplot(dtn[!(`gene expression model` %in% c('GTEx European WB', 'GTEx European Brain')),], aes(x = N_samples, y = N_assocs, color = `gene expression model`)) +
2189. geom_point() +
2190. geom_smooth(method = 'lm', color = 'black', se = F) +
2191. #geom_text(aes(x = -Inf, y = -Inf, label = summary(lmod_g2s)$pvalue)) +
2192. annotate(geom="text", x=Inf, y=-Inf, label=paste('r2=',round(summary(lmod_g2s)$adj.r.squared, 3)), color="red3", hjust = 1.2, vjust = -0.25) +
2193. geom_label_repel(aes(label = `gene expression model`)) +
2194. theme_minimal() +
2195. theme(legend.position = 'none')
2198. name = 'Output/revisions/assoc_by_modelN_g2s.pdf'
2199. pdf(file = name, width = 5, height = 4, pointsize = 12, bg = "white")
2200. print(dtn_plot_g2s)
2201. dev.off()
2202. print(dtn_plot_g2s)
2203. ```
2205. ```{r}
2206. lmod_g2s2 <- lm(data = dtn[!(`gene expression model` %in% c('GTEx European WB', 'GTEx European Brain', 'G2S All WB')),], N_assocs ~ N_samples)
2208. summary(lmod_g2s2)
2210. dtn_plot_g2s2 <- ggplot(dtn[!(`gene expression model` %in% c('GTEx European WB', 'GTEx European Brain', 'G2S All WB')),], aes(x = N_samples, y = N_assocs, color = `gene expression model`)) +
2211. geom_point() +
2212. geom_smooth(method = 'lm', color = 'black', se = F) +
2213. annotate(geom="text", x=Inf, y=-Inf, label=paste('r2=',round(summary(lmod_g2s2)$adj.r.squared, 3)), color="red3", hjust = 1.2, vjust = -0.25) +
2214. geom_label_repel(aes(label = `gene expression model`)) +
2215. theme(legend.position = 'none')
2218. name = 'Output/revisions/assoc_by_modelN_g2s_single.pdf'
2219. pdf(file = name, width = 5, height = 4, pointsize = 12, bg = "white")
2220. print(dtn_plot_g2s2)
2221. dev.off()
2222. print(dtn_plot_g2s2)
2224. ```
2226. ```{r}
2228. lmod_g2s2b <- lm(data = dtn[!(`gene expression model` %in% c('GTEx European WB', 'GTEx European Brain', 'G2S All WB')),], N_assocs ~ Mean_Afr)
2230. summary(lmod_g2s2b)
2232. dtn_plot_g2s2b <- ggplot(dtn[!(`gene expression model` %in% c('GTEx European WB', 'GTEx European Brain', 'G2S All WB')),], aes(x = Mean_Afr, y = N_assocs, color = `gene expression model`)) +
2233. geom_point() +
2234. geom_smooth(method = 'lm', color = 'black', se = F) +
2235. scale_color_brewer(palette = 'Accent') +
2236. annotate(geom="text", x=Inf, y=-Inf, label=paste('r2=',round(summary(lmod_g2s2b)$adj.r.squared, 3)), color="red3", hjust = 1.2, vjust = -0.25) +
2237. geom_label_repel(aes(label = `gene expression model`)) +
2238. theme(legend.position = 'none')
2241. name = 'Output/revisions/assoc_by_Afr_g2s.pdf'
2242. pdf(file = name, width = 5, height = 4, pointsize = 12, bg = "white")
2243. print(dtn_plot_g2s2b)
2244. dev.off()
2245. print(dtn_plot_g2s2b)
2247. ```
2249. ```{r}
2250. lmod_g2s2c <- lm(data = dtn[!(`gene expression model` %in% c('GTEx European WB', 'GTEx European Brain', 'G2S All WB')),], N_assocs ~ Mean_IAM)
2252. dtn_plot_g2s2c <- ggplot(dtn[!(`gene expression model` %in% c('GTEx European WB', 'GTEx European Brain', 'G2S All WB')),], aes(x = Mean_IAM, y = N_assocs, color = `gene expression model`)) +
2253. geom_point() +
2254. geom_smooth(method = 'lm', color = 'black', se = F) +
2255. scale_color_brewer(palette = 'Accent') +
2256. annotate(geom="text", x=Inf, y=-Inf, label=paste('r2=',round(summary(lmod_g2s2c)$adj.r.squared, 3)), color="red3", hjust = 1.2, vjust = -0.25) +
2257. geom_label_repel(aes(label = `gene expression model`)) +
2258. theme(legend.position = 'none')
2261. name = 'Output/revisions/assoc_by_IAM_g2s.pdf'
2262. pdf(file = name, width = 5, height = 4, pointsize = 12, bg = "white")
2263. print(dtn_plot_g2s2c)
2264. dev.off()
2266. print(dtn_plot_g2s2c)
2267. ```
2271. ## Werth Replication Analysis
2273. ### Figure S5: Comparisons of TWAS vs prior discoveries
2274. ```{r}
2276. library(data.table)
2277. library(ggplot2)
2278. library(reshape2)
2279. library(pals)
2280. library(stringr)
2283. dt_ <- fread('Data/all_exc_BHsig_assocs_PGCancestryTWAS.txt', header = TRUE)
2285. # exclude replication/comparison data from initial analysis.
2286. gwas <- c('PGC_AA_MDD_2023', 'PGC_AA_PTSD_2024','PGC_AA_SCZ_2022','PGC_ASIAN_SCZ_2022','PGC_EUR_PTSD_2024','PGC_HNA_PTSD_2024', 'PGC_EURO_SCZ_2022')
2287. mod <- c('PrediXcan_Brain_Cortex')
2289. dt <- dt_[!(gwas_phenotype %in% gwas | training_model %in% mod),]
2292. euro <- dt[training_model == 'PrediXcan_Whole_Blood',]
2293. anc <- dt[training_model != 'PrediXcan_Whole_Blood',]
2294. dup <- setDT(merge(anc, euro, by = c('gene_name', 'gwas_phenotype')))
2296. euro$gptm <- paste0(euro$gene_name, '_', euro$gwas_phenotype)
2297. anc$gptm <- paste0(anc$gene_name, '_', anc$gwas_phenotype)
2299. anc2 <- anc[!(gptm %in% euro$gptm),]
2300. euro2 <- euro[!(gptm %in% anc$gptm),]
2304. #Run ADHD replication analysis.
2306. adhd_w <- fread('Data/twas_atlas_data/replication_datasets/adhd_22_werth.txt', header = T)
2308. adhd_w_wb <- adhd_w[, c('gene_name', 'Whole_Blood_unconditioned')][!is.na(Whole_Blood_unconditioned),]
2309. adhd_w_wb$fdr_p <- p.adjust(adhd_w_wb$Whole_Blood_unconditioned, 'BH')
2311. adhd_ta <- adhd_w_wb[Whole_Blood_unconditioned <= 0.05,]
2313. adhd_anc_specific <- anc2[gwas_phenotype == 'PGC_allADHD_2022' & !(gene_name %in% adhd_ta$gene_name),]
2314. a1 <- nrow(adhd_anc_specific[, .(.N), by = 'gene_name'])
2315. adhd_euro_specific <- euro2[gwas_phenotype == 'PGC_allADHD_2022' & !(gene_name %in% adhd_ta$gene_name),]
2316. a2 <- nrow(adhd_euro_specific[, .(.N), by = 'gene_name'])
2317. a3 <- adhd_ta_specific <- length(unique(adhd_ta[!(gene_name %in% dt$gene_name),]$gene_name))
2318. a4 <- length(unique(anc2[gwas_phenotype == 'PGC_allADHD_2022' & (gene_name %in% adhd_ta$gene_name),]$gene_name))
2319. a5 <- length(unique(euro2[gwas_phenotype == 'PGC_allADHD_2022' & (gene_name %in% adhd_ta$gene_name),]$gene_name))
2323. # then BD
2325. bd_w <- fread('Data/twas_atlas_data/replication_datasets/bd_22_werth.txt', header = T)
2327. bd_w_wb <- bd_w[, c('gene_name', 'Whole_Blood_unconditioned')][!is.na(Whole_Blood_unconditioned),]
2328. bd_w_wb$fdr_p <- p.adjust(bd_w_wb$Whole_Blood_unconditioned, 'BH')
2330. bd_ta <- bd_w_wb[Whole_Blood_unconditioned <= 0.05,]
2332. bd_anc_specific <- anc2[gwas_phenotype == 'PGC_BD1_2021' & !(gene_name %in% bd_ta$gene_name),]
2333. b1 <- nrow(bd_anc_specific[, .(.N), by = 'gene_name'])
2334. bd_euro_specific <- euro2[gwas_phenotype == 'PGC_BD1_2021' & !(gene_name %in% bd_ta$gene_name),]
2335. b2 <- nrow(bd_euro_specific[, .(.N), by = 'gene_name'])
2336. b3 <- bd_ta_specific <- length(unique(bd_ta[!(gene_name %in% dt$gene_name),]$gene_name))
2337. b4 <- length(unique(anc2[gwas_phenotype == 'PGC_BD1_2021' & (gene_name %in% bd_ta$gene_name),]$gene_name))
2338. b5 <- length(unique(euro2[gwas_phenotype == 'PGC_BD1_2021' & (gene_name %in% bd_ta$gene_name),]$gene_name))
2343. # and MDD
2345. mdd_w <- fread('Data/twas_atlas_data/replication_datasets/mdd_22_werth.txt', header = T)
2347. mdd_w_wb <- mdd_w[, c('gene_name', 'Whole_Blood_unconditioned')][!is.na(Whole_Blood_unconditioned),]
2348. mdd_w_wb$fdr_p <- p.adjust(mdd_w_wb$Whole_Blood_unconditioned, 'BH')
2350. mdd_ta <- mdd_w_wb[Whole_Blood_unconditioned <= 0.05,]
2352. mdd_anc_specific <- anc2[gwas_phenotype == 'PGC_allMDD_2023' & !(gene_name %in% mdd_ta$gene_name),]
2353. m1 <- nrow(mdd_anc_specific[, .(.N), by = 'gene_name'])
2354. mdd_euro_specific <- euro2[gwas_phenotype == 'PGC_allMDD_2023' & !(gene_name %in% mdd_ta$gene_name),]
2355. m2 <- nrow(mdd_euro_specific[, .(.N), by = 'gene_name'])
2356. m3 <- mdd_ta_specific <- length(unique(mdd_ta[!(gene_name %in% dt$gene_name),]$gene_name))
2357. m4 <- length(unique(anc2[gwas_phenotype == 'PGC_allMDD_2023' & (gene_name %in% mdd_ta$gene_name),]$gene_name))
2358. m5 <- length(unique(euro2[gwas_phenotype == 'PGC_allMDD_2023' & (gene_name %in% mdd_ta$gene_name),]$gene_name))
2360. # and lastly scz
2363. scz_w <- fread('Data/twas_atlas_data/replication_datasets/scz_22_werth.txt', header = T)
2365. scz_w_wb <- scz_w[, c('gene_name', 'Whole_Blood_unconditioned')][!is.na(Whole_Blood_unconditioned),]
2366. scz_w_wb$fdr_p <- p.adjust(scz_w_wb$Whole_Blood_unconditioned, 'BH')
2368. scz_ta <- scz_w_wb[Whole_Blood_unconditioned <= 0.05,]
2370. scz_anc_specific <- anc2[gwas_phenotype == 'PGC_all_SCZ_2022' & !(gene_name %in% scz_ta$gene_name),]
2371. s1 <- nrow(scz_anc_specific[, .(.N), by = 'gene_name'])
2372. scz_euro_specific <- euro2[gwas_phenotype == 'PGC_all_SCZ_2022' & !(gene_name %in% scz_ta$gene_name),]
2373. s2 <- nrow(scz_euro_specific[, .(.N), by = 'gene_name'])
2374. s3 <- scz_ta_specific <- length(unique(scz_ta[!(gene_name %in% dt$gene_name),]$gene_name))
2375. s4 <- length(unique(anc2[gwas_phenotype == 'PGC_all_SCZ_2022' & (gene_name %in% scz_ta$gene_name),]$gene_name))
2376. s5 <- length(unique(euro2[gwas_phenotype == 'PGC_all_SCZ_2022' & (gene_name %in% scz_ta$gene_name),]$gene_name))
2379. #compile data into single object
2380. ant_twas_dat <- data.table(study = c('ADHD', 'MDD', 'BD1', 'SCZ'),
2381. g2s_new = c(a1, m1, b1, s1),
2382. euro_new = c(a2, m2, b2, s2),
2383. Werth_unique = c(a3,m3,b3,s3),
2384. g2s_rep = c(a4,m4,b4,s4),
2385. euro_rep = c(a5, m5, b5, s5))
2388. atd <- setDT(melt(ant_twas_dat, id.vars = 'study'))
2389. atd[variable %in% c('g2s_new', 'g2s_rep'), cohort := 'G2S']
2390. atd[variable %in% c('euro_new', 'euro_rep'), cohort := 'GTEx']
2391. atd[variable %in% c('Werth_unique'), cohort := 'Werth']
2393. tots <- atd[, .(value = sum(value)), by = c('variable', 'cohort')]
2394. tots$study = 'All Conditions'
2396. atd2 <- rbind(atd, tots)
2399. #visualize data
2401. p10 <- ggplot(atd[variable %in% c('g2s_rep', 'euro_rep'),], aes(x = cohort, y = value, color = cohort, fill = cohort)) +
2402. geom_bar(stat = 'identity', position = 'dodge') +
2403. geom_label(aes(label = value, y = value +2), fill = "white", color = 'black', label.size = NA) +
2404. facet_grid(.~study) +
2405. ylab('Replicated Genes') +
2406. scale_fill_brewer(palette = 'Dark2') +
2407. scale_color_brewer(palette = 'Dark2') +
2408. theme_minimal() +
2409. theme(axis.text.x = element_text(angle = 60, hjust = 1),
2410. axis.title.x = element_blank(),
2411. legend.position = 'none')
2413. name = 'Output/revisions/werth_rep.pdf'
2414. pdf(file = name, width = 7, height = 3, pointsize = 12, bg = "white")
2415. print(p10)
2416. dev.off()
2418. print(p10)
2420. p11 <- ggplot(tots[variable %in% c('g2s_rep', 'euro_rep'),], aes(x = cohort, y = value, color = cohort, fill = cohort)) +
2421. geom_bar(stat = 'identity', position = 'dodge') +
2422. geom_label(aes(label = value, y = value +4), fill = "white", color = 'black', label.size = NA) +
2423. facet_grid(.~study) +
2424. ylab('Replicated Genes') +
2425. scale_fill_brewer(palette = 'Dark2') +
2426. scale_color_brewer(palette = 'Dark2') +
2427. theme_minimal() +
2428. theme(axis.text.x = element_text(angle = 60, hjust = 1),
2429. axis.title = element_blank(),
2430. legend.position = 'none')
2432. name = 'Output/revisions/werth_total_rep.pdf'
2433. pdf(file = name, width = 1.8, height = 3, pointsize = 12, bg = "white")
2434. print(p11)
2435. dev.off()
2437. ```
2439. ## Pvalue vs Zscore monotonicity
2441. ### Fig S10: Pvalue change relative to zscores
2442. ```{r}
2444. dt <- fread("Data/all_exc_assocs_PGCancestryTWAS.txt", header = T)
2448. plotdata <- dt[gwas_phenotype == 'PGC_allADHD_2022' & training_model %in% c('PrediXcan_Whole_Blood', 'AA_Whole_Blood')]
2451. pz_comp <- setDT(merge(plotdata[training_model == 'PrediXcan_Whole_Blood', c('gene_name', 'gwas_phenotype', 'zscore', 'pvalue')],
2452. plotdata[training_model == 'AA_Whole_Blood', c('gene_name', 'gwas_phenotype', 'zscore', 'pvalue')],
2453. suffixes = c('.gtex', '.g2s_AA'), by = c('gene_name', 'gwas_phenotype')))
2457. pz_comp[, zdif := abs(zscore.gtex - zscore.g2s_AA)]
2458. pz_comp[, pdif := abs(pvalue.gtex - pvalue.g2s_AA)]
2460. zp_gtg2 <- ggplot(pz_comp, aes(x = zscore.gtex, y = pvalue.g2s_AA)) +
2461. geom_point(color = 'seagreen')
2464. zp_g2g2 <- ggplot(pz_comp, aes(x = zscore.g2s_AA, y = pvalue.g2s_AA)) +
2465. geom_point(color = 'magenta')
2468. zp_g2gt <- ggplot(pz_comp, aes(x = zscore.g2s_AA, y = pvalue.gtex)) +
2469. geom_point(color = 'goldenrod2')
2472. zp_gtgt <- ggplot(pz_comp, aes(x = zscore.gtex, y = pvalue.gtex)) +
2473. geom_point(color = 'slateblue')
2477. z_vs_p_plot <- ggarrange(zp_gtg2, zp_g2g2, zp_gtgt, zp_g2gt, ncol = 2, nrow = 2)
2479. pdf(file = 'Output/revisions/z_vs_p_analysis.pdf', height = 9, width = 9)
2480. print(z_vs_p_plot)
2481. dev.off()
2483. print(z_vs_p_plot)
2485. #
2488. zp_gtg2 <- ggplot(pz_comp, aes(x = zscore.gtex, y = pvalue.g2s_AA)) +
2489. geom_point(color = 'seagreen')
2492. zp_g2g2 <- ggplot(pz_comp, aes(x = zscore.g2s_AA, y = pvalue.g2s_AA)) +
2493. geom_point(color = 'magenta')
2496. zp_g2gt <- ggplot(pz_comp, aes(x = zscore.g2s_AA, y = pvalue.gtex)) +
2497. geom_point(color = 'goldenrod2')
2500. zp_gtgt <- ggplot(pz_comp, aes(x = zscore.gtex, y = pvalue.gtex)) +
2501. geom_point(color = 'slateblue')
2506. z_vs_p_plot <- ggarrange(zp_gtg2, zp_g2g2, zp_gtgt, zp_g2gt, ncol = 2, nrow = 2)
2508. pdf(file = 'Output/revisions/z_vs_p_analysis.pdf', height = 9, width = 9)
2509. print(z_vs_p_plot)
2510. dev.off()
2516. pp_gtg2 <- ggplot(pz_comp, aes(x = pvalue.gtex, y = pvalue.g2s_AA)) +
2517. geom_point(color = 'seagreen')
2520. pp_g2g2 <- ggplot(pz_comp, aes(x = pvalue.g2s_AA, y = pvalue.g2s_AA)) +
2521. geom_point(color = 'magenta')
2524. pp_g2gt <- ggplot(pz_comp, aes(x = pvalue.g2s_AA, y = pvalue.gtex)) +
2525. geom_point(color = 'goldenrod2')
2528. pp_gtgt <- ggplot(pz_comp, aes(x = pvalue.gtex, y = pvalue.gtex)) +
2529. geom_point(color = 'slateblue')
2534. p_vs_p_plot <- ggarrange(pp_gtg2, pp_g2g2, pp_gtgt, pp_g2gt, ncol = 2, nrow = 2)
2536. pdf(file = 'Output/revisions/p_vs_p_analysis.pdf', height = 9, width = 9)
2537. print(p_vs_p_plot)
2538. dev.off()
2540. print(p_vs_p_plot)
2544. zz_gtg2 <- ggplot(pz_comp, aes(x = zscore.gtex, y = zscore.g2s_AA)) +
2545. geom_point(color = 'seagreen')
2548. zz_g2g2 <- ggplot(pz_comp, aes(x = zscore.g2s_AA, y = zscore.g2s_AA)) +
2549. geom_point(color = 'magenta')
2552. zz_g2gt <- ggplot(pz_comp, aes(x = zscore.g2s_AA, y = zscore.gtex)) +
2553. geom_point(color = 'goldenrod2')
2556. zz_gtgt <- ggplot(pz_comp, aes(x = zscore.gtex, y = zscore.gtex)) +
2557. geom_point(color = 'slateblue')
2562. z_vs_z_plot <- ggarrange(zz_gtg2, zz_g2g2, zz_gtgt, zz_g2gt, ncol = 2, nrow = 2)
2564. pdf(file = 'Output/revisions/z_vs_z_analysis.pdf', height = 9, width = 9)
2565. print(z_vs_z_plot)
2566. dev.off()
2568. print(z_vs_z_plot)
2571. #####
2573. pzc <- copy(pz_comp)
2574. pzc[, pdif := pvalue.gtex - pvalue.g2s_AA]
2575. pzc[, zdif := abs(zscore.gtex) - abs(zscore.g2s_AA)]
2577. pzc[, zpsign := sign(zdif) * sign(pdif)]
2580. zmagpplot <- ggplot(pzc, aes(x = zdif, y = pdif)) +
2581. geom_point()
2583. pdf(file = 'Output/revisions/deltazmagvsp.pdf', height = 9, width = 9)
2584. print(zmagpplot)
2585. dev.off()
2587. print(zmagpplot)
2592. azp_gtg2 <- ggplot(pz_comp, aes(x = abs(zscore.gtex), y = pvalue.g2s_AA)) +
2593. geom_point(color = 'seagreen')
2596. azp_g2g2 <- ggplot(pz_comp, aes(x = abs(zscore.g2s_AA), y = pvalue.g2s_AA)) +
2597. geom_point(color = 'magenta')
2600. azp_g2gt <- ggplot(pz_comp, aes(x = abs(zscore.g2s_AA), y = pvalue.gtex)) +
2601. geom_point(color = 'goldenrod2')
2604. azp_gtgt <- ggplot(pz_comp, aes(x = abs(zscore.gtex), y = pvalue.gtex)) +
2605. geom_point(color = 'slateblue')
2609. zm_vs_p_plot <- ggarrange(azp_gtg2, azp_g2g2, azp_gtgt, azp_g2gt, ncol = 2, nrow = 2)
2611. pdf(file = 'Output/revisions/zmag_vs_p_analysis.pdf', height = 9, width = 9)
2612. print(zm_vs_p_plot)
2613. dev.off()
2615. print(zm_vs_p_plot)
2618. ```
2621. ## Fig S18-20: Linkage Disequilibrium analyses
2623. How does linkage disequilibrium affect the relationship of genes that are called significant for a phenotype?
2625. ```{r}
2628. lddir <- 'Data/LD_genelists_shared_distinct/'
2632. assocs <- c ('ENSG00000172247_MX_Whole_Blood_SUD_alc_2019_snps_20',
2633. # 'ENSG00000237513_AA_Whole_Blood_PGC_ALL_PTSD_2024_snps_15',
2634. 'ENSG00000237513_AA_Whole_Blood_PGC_all_SCZ_2022_snps_15',
2635. 'ENSG00000255284_MX_Whole_Blood_PGC_allADHD_2022_snps_17')
2638. assoc <- c ('ENSG00000172247_MX_Whole_Blood_SUD_alc_2019_snps_20')
2641. pops <- c('CEU', 'GBR', 'YRI', 'MXL', 'PUR', 'ASW')
2643. for(assoc in assocs){
2644. tab_all <- fread(paste0(lddir, 'LD_results2/',assoc,'_LD_all.txt'), header = F)
2645. tab_all$population <- str_remove(tab_all$V11, '.+:')
2647. pop_charts <- c()
2648. for( pop in pops){
2649. tab <- tab_all[population == pop,]
2650. head(pop)
2651. setnames(tab, c('V1', 'V5', 'V9', 'V10'), c('SNP1', 'SNP2', 'R2', 'Dprime'))
2652. init <- fread(paste0(lddir, assoc, '.txt'), header = T)
2654. g2s_rsids <- init[!is.na(weight.G2S),]$rsid
2655. gtex_rsids <- init[!is.na(weight.GTEx),]$rsid
2658. comp <- tab[, c('SNP1', 'SNP2', 'R2', 'Dprime')]
2660. comp2 <- comp[(SNP1 %in% g2s_rsids & SNP2 %in% gtex_rsids),]
2661. comp3 <- comp[(SNP2 %in% g2s_rsids & SNP1 %in% gtex_rsids) , ]
2662. comp3 <- setnames(comp3, c('SNP2' , 'SNP1'), c('SNP1' , 'SNP2'))
2663. comp4 <- setDT(rbind(comp2, comp3))
2664. setnames(comp4, c('SNP1' , 'SNP2'), c('G2S', 'GTEx'))
2666. mattabl <- matrix(, nrow = length(g2s_rsids), ncol = length(gtex_rsids))
2667. rownames(mattabl) <- g2s_rsids
2668. colnames(mattabl) <- gtex_rsids
2670. for(i in 1:nrow(comp4)){
2671. g2 <- comp4[i,]$G2S
2672. gt <- comp4[i,]$GTEx
2673. r2 <- comp4[i,]$R2
2674. mattabl[g2, gt] <- r2
2675. }
2678. snp_locs1 <- tab[, c(1,2)]
2679. snp_locs1$pos <- str_remove(snp_locs1$V2, '.+:')
2680. snp_locs1 <- unique(snp_locs1)
2681. setnames(snp_locs1, 'SNP1', 'SNP')
2683. snp_locs2 <- tab[, c(5,6)]
2684. snp_locs2$pos <- str_remove(snp_locs2$V6, '.+:')
2685. snp_locs2 <- unique(snp_locs2)
2686. setnames(snp_locs2, 'SNP2', 'SNP')
2688. allsnps <- setDT(rbind(snp_locs1[, c('SNP', 'pos')], snp_locs2[, c('SNP', 'pos')]))
2690. g2_ord <- unique(allsnps[SNP %in% g2s_rsids,])[order(pos)]$SNP
2691. gt_ord <- unique(allsnps[SNP %in% gtex_rsids,])[order(pos)]$SNP
2694. comp4$G2S <- factor(comp4$G2S, levels = g2_ord)
2695. comp4$GTEx <- factor(comp4$GTEx, levels = gt_ord)
2697. LDplot1 <- ggplot(comp4, aes( x = G2S, y = GTEx, fill = R2)) +
2698. geom_tile() +
2699. theme_minimal() +
2700. ggtitle(pop) +
2701. scale_fill_distiller(palette = "RdYlGn", direction = 1) +
2702. theme(axis.text.x = element_text(angle = 40, hjust = 1))
2704. pop_charts[[pop]] = LDplot1
2705. }
2706. name <- paste0('Output/revisions/ALL_', assoc,'_LD_corplots.pdf')
2708. pdf(file = name, width = 9, height = 10, pointsize = 12, bg = "white")
2709. do.call(grid.arrange, c(pop_charts, ncol = 2, padding = unit(1, "lines")))
2710. dev.off()
2713. }
2715. #images are output to output directory, not in std out.
2716. ```
2718. ### Fig S21: Max R2 to SNP weight comparison
2720. For each gene, how does the maximum LD between SNPs used in the model compare to the weight of the SNP? We want to assess whether SNPs that are more frequently in LD have greater effects on gene expression.
2721. ```{r}
2722. ######## get relationship between max r2 and SNP weight
2725. library(data.table)
2726. library(ggplot2)
2727. library(gridExtra)
2728. library(stringr)
2731. lddir <- 'Data/LD_genelists_shared_distinct/'
2735. assocs <- c('ENSG00000172247_MX_Whole_Blood_SUD_alc_2019_snps_20',
2736. #'ENSG00000237513_AA_Whole_Blood_PGC_ALL_PTSD_2024_snps_15',
2737. 'ENSG00000237513_AA_Whole_Blood_PGC_all_SCZ_2022_snps_15',
2738. 'ENSG00000255284_MX_Whole_Blood_PGC_allADHD_2022_snps_17')
2741. pops <- c('CEU', 'GBR', 'YRI', 'MXL', 'PUR', 'ASW')
2745. all_dat <- list.files('Data/rdata', pattern = '.+\\.txt')
2747. snp_dat <- data.table()
2748. for(fil in str_subset(all_dat, '^snp_dat_all.+\\.txt')){
2749. newdat <- fread(paste0('Data/rdata/', fil), header = T)
2750. snp_dat <- rbind(snp_dat , newdat)
2751. }
2754. #initiate data.table to hold r2 and pvalue results
2755. cordat <- data.table(association = c(), gene = c(), training_model = c(), gwas_phenotype = c(), population = c(), group = c(), adj.r.squared = c(), lm_pvalue = c())
2756. #initiate data.table to hold all points from plots
2757. fulldata <- data.table()
2758. for(assoc in assocs){
2759. tab_all <- fread(paste0(lddir, 'LD_results2/',assoc,'_LD_all.txt'), header = F)
2760. tab_all$population <- str_remove(tab_all$V11, '.+:')
2762. gene_nm <- str_remove(assoc, '_.+')
2763. mod1 <- str_remove(assoc, paste0(gene_nm , '_'))
2764. mod <- str_remove(mod1, '_PGC.+')
2765. mod <- str_remove(mod, '_SUD.+')
2766. phen1 <- str_remove(assoc, paste0(gene_nm , '_', mod, '_'))
2767. phen <- str_remove(phen1, '_snps.+')
2769. sdat <- snp_dat[gene == gene_nm & (model == mod & gwas_phenotype == phen),]
2772. pop_charts2 <- c()
2773. for( pop in pops){
2774. tab <- tab_all[population == pop,]
2775. head(pop)
2776. setnames(tab, c('V1', 'V5', 'V9', 'V10'), c('SNP1', 'SNP2', 'R2', 'Dprime'))
2777. init <- fread(paste0(lddir, assoc, '.txt'), header = T)
2779. g2s_rsids <- init[!is.na(weight.G2S),]$rsid
2780. gtex_rsids <- init[!is.na(weight.GTEx),]$rsid
2783. comp <- tab[, c('SNP1', 'SNP2', 'R2', 'Dprime')]
2785. comp2 <- comp[(SNP1 %in% g2s_rsids & SNP2 %in% gtex_rsids),]
2786. comp3 <- comp[(SNP2 %in% g2s_rsids & SNP1 %in% gtex_rsids) , ]
2787. comp3 <- setnames(comp3, c('SNP2' , 'SNP1'), c('SNP1' , 'SNP2'))
2788. comp4 <- setDT(rbind(comp2, comp3))
2789. setnames(comp4, c('SNP1' , 'SNP2'), c('G2S', 'GTEx'))
2792. g2s_max1 <- comp4[, .(maxr2 = max(R2)), by = 'G2S']
2793. gtex_max1 <- comp4[, .(maxr2 = max(R2)), by = 'GTEx']
2795. g2s_max <- setDT(merge(g2s_max1, sdat[, c('rsid', 'weight.G2S')], by.x = 'G2S', by.y = 'rsid'))
2796. g2s_max$group <- 'G2S'
2797. gtex_max <- setDT(merge(gtex_max1, sdat[, c('rsid', 'weight.GTEx')], by.x = 'GTEx', by.y = 'rsid'))
2798. gtex_max$group <- 'GTEx'
2800. setnames(g2s_max, c('G2S', 'weight.G2S'), c('SNP', 'weight'))
2801. setnames(gtex_max, c('GTEx', 'weight.GTEx'),c('SNP', 'weight'))
2803. lmg2s = lm(abs(g2s_max$weight)~g2s_max$maxr2)
2804. lmgtex = lm(abs(gtex_max$weight)~gtex_max$maxr2)
2806. newdat_g2s = data.table(association = assoc,
2807. gene = gene_nm,
2808. training_model = mod,
2809. gwas_phenotype = phen,
2810. population = pop,
2811. group = 'G2S',
2812. adj.r.squared = summary(lmg2s)$adj.r.squared,
2813. lm_pvalue = summary(lmg2s)$coefficients[2,4])
2815. newdat_gtex = data.table(association = assoc,
2816. gene = gene_nm,
2817. training_model = mod,
2818. gwas_phenotype = phen,
2819. population = pop,
2820. group = 'GTEx',
2821. adj.r.squared = summary(lmgtex)$adj.r.squared,
2822. lm_pvalue = summary(lmgtex)$coefficients[2,4])
2825. cordat <- rbind(cordat, newdat_g2s)
2826. cordat <- rbind(cordat, newdat_gtex)
2828. plotdata <- setDT(rbind(g2s_max, gtex_max))
2829. plotdata[, association := assoc]
2830. plotdata[, population := pop]
2831. fulldata <- setDT(rbind(fulldata, plotdata))
2833. }
2834. }
2837. write.table(cordat, file = paste0(lddir, 'LD_results2/all_results_correlation_data.txt'), col.names = T, row.names = F, quote = F, sep = '\t')
2838. write.table(fulldata, file = paste0(lddir, 'LD_results2/all_results_r2byweight_data.txt'), col.names = T, row.names = F, quote = F, sep = '\t')
2840. ld_all <- fread(paste0(lddir, 'LD_results2/all_results_r2byweight_data.txt'), header = T)
2842. lda_lm <- lm(data = ld_all, abs(weight)~maxr2 + association)
2843. lda_lm2 <- lm(data = ld_all, abs(weight)~maxr2 + association + population)
2844. lda_lm3 <- lm(data = ld_all, abs(weight)~maxr2 + population)
2846. allplot <- ggplot(ld_all, aes(x = maxr2, y = abs(weight), color = association)) +
2847. geom_smooth(method = "lm", color = "grey70") +
2848. geom_point() +
2849. theme_minimal() +
2850. facet_wrap(.~population) +
2851. guides(color = guide_legend(nrow = 3)) +
2852. theme(legend.position = 'bottom')
2856. name <- paste0('Output/revisions/ALL_LD_corplots_weights.pdf')
2858. pdf(file = name, width = 6, height = 6, pointsize = 12, bg = "white")
2859. print(allplot)
2860. dev.off()
2862. print(allplot)
2864. ```
2866. ## Correlogram and heatmap in fig 5
2868. ### Fig 5: Updated correlogram and heatmap
2869. ```{r}
2871. library(data.table)
2872. library(corrplot)
2875. dt_ <- fread("Data/all_exc_assocs_PGCancestryTWAS.txt", header = T)
2878. # exclude replication/comparison data from initial analysis.
2879. gwas <- c('PGC_AA_MDD_2023', 'PGC_AA_PTSD_2024','PGC_AA_SCZ_2022','PGC_ASIAN_SCZ_2022','PGC_EUR_PTSD_2024','PGC_HNA_PTSD_2024')
2880. mod <- c('PrediXcan_Brain_Cortex')
2882. dt <- dt_[!(gwas_phenotype %in% c('PGC_AA_MDD_2023', 'PGC_AA_PTSD_2024', 'PGC_LATINO_SCZ_2022', 'PGC_AA_SCZ_2022','PGC_ASIAN_SCZ_2022','PGC_EUR_PTSD_2024','PGC_HNA_PTSD_2024', 'PGC_EURO_SCZ_2022')) & training_model != 'PrediXcan_Brain_Cortex',]
2885. dt[training_model == 'All_Whole_Blood', `gene expression model` := 'G2S All WB']
2886. dt[training_model == 'AA_Whole_Blood', `gene expression model` := 'G2S African American WB']
2887. dt[training_model == 'MX_Whole_Blood', `gene expression model` := 'G2S Mexican American WB']
2888. dt[training_model == 'PR_Whole_Blood', `gene expression model` := 'G2S Puerto Rican WB']
2889. dt[training_model == 'PrediXcan_Whole_Blood', `gene expression model` := 'GTEx European WB']
2891. dt1 <- dt[!training_model == 'PrediXcan_Whole_Blood']
2892. dt2 <- dt[training_model == 'PrediXcan_Whole_Blood']
2896. # Clean gene names in both datasets
2897. MAtwas <- dt2
2898. MUtwas <- dt1
2901. # Get unique training models from MUtwas and add GTEX from MAtwas
2902. training_models <- unique(MUtwas$training_model)
2903. training_models <- c(training_models, 'PrediXcan_Whole_Blood')
2905. # Define a custom color palette
2906. custom_colors <- colorRampPalette(c("#1a9850","#1a9850", "#ffffbf", "#d73027"))(1000)
2908. # Open a PDF device with Letter dimensions
2909. pdf('Output/revisions/correlograms.pdf', width = 5, height = 5.5)
2910. for (model in training_models[1]) {
2911. # Filter data for the current training model
2912. if (model == 'PrediXcan_Whole_Blood') {
2913. filtered_data <- MAtwas[training_model == model, .(gwas_phenotype, gene, zscore)]
2914. } else {
2915. filtered_data <- MUtwas[training_model == model, .(gwas_phenotype, gene, zscore)]
2916. }
2918. unique_genes <- unique(filtered_data$gene)
2919. unique_phenotypes <- unique(filtered_data$gwas_phenotype)
2921. matrix_data <- matrix(NA, nrow = length(unique_genes), ncol = length(unique_phenotypes))
2922. rownames(matrix_data) <- unique_genes
2924. dz_names <- c('PGC_allADHD_2022' = 'ADHD',
2925. 'PGC_all_SCZ_2022' = 'SCZ',
2926. 'PGC_allMDD_2023' = 'MDD',
2927. 'PGC_BD1_2021' = 'BD1',
2928. 'SUD_alc_2019' = 'SUD-A',
2929. 'PGC_ALL_PTSD_2024' = 'PTSD')
2931. colnames(matrix_data) <- unique_phenotypes
2932. colnames(matrix_data) <- dz_names[colnames(matrix_data)]
2934. for (i in 1:nrow(filtered_data)) {
2935. row <- filtered_data[i, ]
2936. matrix_data[row$gene, dz_names[row$gwas_phenotype]] <- row$zscore
2937. }
2939. cor_matrix <- cor(matrix_data, use = "pairwise.complete.obs")
2940. cor_matrix[cor_matrix == 1.0] <- NA
2942. # Adjust plot margins to accommodate longer titles
2943. par(mar = c(1, 1, 5, 1))
2945. corrplot(corr = cor_matrix,
2946. method = "circle",
2947. type = "lower",
2948. #tl.pos = "tl",
2949. order = "original",
2950. tl.col = 'black',
2951. diag = FALSE,
2952. is.corr = FALSE,
2953. col.lim = range(cor_matrix, na.rm = T),
2954. col = custom_colors, # Use custom colors
2955. title = paste("Correlation plot for training model:", model),
2956. mar = c(0, 0, 2, 0)) # Adjust margins to fit the title
2957. }
2959. # Close the PDF device
2960. dev.off()
2963. ```
2966. ```{r}
2968. # and heat map too
2969. cleardups <- function(t1){
2970. otpt <- t1[2,]
2971. for(nn in 3:nrow(t1)){
2972. d1 <- t1[nn,1]
2973. d2 <- t1[nn,2]
2974. d_combo <- paste0(d1, d2)
2975. d_combo2 <- paste0(d2, d1)
2976. if( d_combo %in% paste0(otpt[[1]], otpt[[2]]) | d_combo2 %in% paste0(otpt[[1]], otpt[[2]])){
2977. #print(TRUE)
2978. } else if(d1 != d2){
2979. #print(FALSE)
2980. otpt <- setDT(rbind(otpt, t1[nn,]))
2981. }
2982. }
2983. otpt[1 != 2,]
2984. otpt[[3]] <- rank(-otpt[[3]])
2985. return(otpt)
2986. }
2989. load('Data/correlation_matrices/AA_Whole_Blood_cor_matrix.rdata')
2990. cortab <- as.data.table(cor_matrix)
2991. cortab$disease1 <- rownames(cor_matrix)
2992. cortab <- setnames(melt(cortab, id.vars = 'disease1', value.name = 'G2S_AA'), 'variable', 'disease2')
2993. cortab$disease2 <- as.character(cortab$disease2)
2994. cortab <- cleardups(cortab)
2995. AA_cor <- copy(cortab)
2997. load('Data/correlation_matrices/All_Whole_Blood_cor_matrix.rdata')
2998. cortab <- as.data.table(cor_matrix)
2999. cortab$disease1 <- rownames(cor_matrix)
3000. cortab <- setnames(melt(cortab, id.vars = 'disease1', value.name = 'G2S_All'), 'variable', 'disease2')
3001. cortab$disease2 <- as.character(cortab$disease2)
3002. cortab <- cleardups(cortab)
3003. All_cor <- copy(cortab)
3005. load('Data/correlation_matrices/MX_Whole_Blood_cor_matrix.rdata')
3006. cortab <- as.data.table(cor_matrix)
3007. cortab$disease1 <- rownames(cor_matrix)
3008. cortab <- setnames(melt(cortab, id.vars = 'disease1', value.name = 'G2S_MX'), 'variable', 'disease2')
3009. cortab$disease2 <- as.character(cortab$disease2)
3010. cortab <- cleardups(cortab)
3011. MX_cor <- copy(cortab)
3013. load('Data/correlation_matrices/PR_Whole_Blood_cor_matrix.rdata')
3014. cortab <- as.data.table(cor_matrix)
3015. cortab$disease1 <- rownames(cor_matrix)
3016. cortab <- setnames(melt(cortab, id.vars = 'disease1', value.name = 'G2S_PR'), 'variable', 'disease2')
3017. cortab$disease2 <- as.character(cortab$disease2)
3018. cortab <- cleardups(cortab)
3019. PR_cor <- copy(cortab)
3021. load('Data/correlation_matrices/PrediXcan_Whole_Blood_cor_matrix.rdata')
3022. cortab <- as.data.table(cor_matrix)
3023. cortab$disease1 <- rownames(cor_matrix)
3024. cortab <- setnames(melt(cortab, id.vars = 'disease1', value.name = 'GTEx_WB'), 'variable', 'disease2')
3025. cortab$disease2 <- as.character(cortab$disease2)
3026. cortab <- cleardups(cortab)
3027. GT_cor <- copy(cortab)
3029. multicor <- data.table()
3030. multicor <- setDT(merge(AA_cor, All_cor, by = c('disease1', 'disease2')))
3031. multicor <- setDT(merge(multicor, MX_cor, by = c('disease1', 'disease2')))
3032. multicor <- setDT(merge(multicor, PR_cor, by = c('disease1', 'disease2')))
3033. multicor <- setDT(merge(multicor, GT_cor, by = c('disease1', 'disease2')))
3035. rowmean <- rowMeans(multicor[,-c(1,2)])
3036. multicor$mean <- rowmean
3037. multicor$traits = paste0(multicor$disease1, ' - ', multicor$disease2)
3039. mcL <- setDT(melt(multicor, id.vars = c('disease1', 'disease2', 'traits')))
3041. order_vec <- mcL[variable == 'mean', c('variable', 'value', 'traits')][order(-value)]$traits
3043. mcL$traits <- factor(mcL$traits, levels = order_vec)
3049. cor_rank <- ggplot(mcL[variable != 'mean',], aes(x = variable, y = traits, fill = value)) +
3050. geom_tile(color = 'black') +
3051. geom_text(aes(label = value)) +
3052. scale_fill_distiller(palette = "RdYlGn", direction = 1) +
3053. xlab('model') +
3054. ylab('disease pair')+
3055. labs(fill = "Correlation\nRanking") +
3056. theme_minimal() +
3057. theme(axis.text.x = element_text(angle = 40, hjust = 1))
3060. name = 'Output/revisions/mod_cor_rank.pdf'
3061. pdf(file = name, width = 6.5, height = 6, pointsize = 12, bg = "white")
3062. print(cor_rank)
3063. dev.off()
3064. print(cor_rank)
3066. ```
3068. ### Figure S26: brain model vs whole blood model comparison
3069. ```{r}
3070. # Create comparison for brain vs G2S vs whole blood GTEx
3072. # Read in all significant data and subset GWAS findings of interest
3074. sig <- fread('Data/all_exc_BHsig_assocs_PGCancestryTWAS.txt', header = T)
3077. dt <- sig[!(gwas_phenotype %in% c('PGC_AA_MDD_2023', 'PGC_AA_PTSD_2024','PGC_AA_SCZ_2022','PGC_ASIAN_SCZ_2022','PGC_EUR_PTSD_2024','PGC_HNA_PTSD_2024', 'PGC_EURO_SCZ_2022')),]
3080. dt[training_model == 'All_Whole_Blood', `gene expression model` := 'G2S All WB']
3081. dt[training_model == 'AA_Whole_Blood', `gene expression model` := 'G2S African American WB']
3082. dt[training_model == 'MX_Whole_Blood', `gene expression model` := 'G2S Mexican American WB']
3083. dt[training_model == 'PR_Whole_Blood', `gene expression model` := 'G2S Puerto Rican WB']
3084. dt[training_model == 'PrediXcan_Whole_Blood', `gene expression model` := 'GTEx European WB']
3085. dt[training_model == 'PrediXcan_Brain_Cortex', `gene expression model` := 'GTEx European Brain']
3090. G2S_mods <- c('All_Whole_Blood', 'AA_Whole_Blood', 'MX_Whole_Blood', 'PR_Whole_Blood')
3092. mod_dz_tabl <- data.table()
3093. for( dz in unique(dt$gwas_phenotype)){
3094. gns <- unique(dt[gwas_phenotype == dz,]$gene_name)
3095. temp <- data.table( gwas_phenotype = dz, gene_name = gns )
3096. temp[, GTEx_WB := gene_name %in% unique(dt[gwas_phenotype == dz & training_model == 'PrediXcan_Whole_Blood',]$gene_name)]
3097. temp[, GTEx_Brain := gene_name %in% unique(dt[gwas_phenotype == dz & training_model == 'PrediXcan_Brain_Cortex',]$gene_name)]
3098. temp[, G2S := gene_name %in% unique(dt[gwas_phenotype == dz & training_model %in% G2S_mods,]$gene_name)]
3100. temp[ (GTEx_WB == F & GTEx_Brain == F) & G2S == F, type := 'none']
3101. temp[ (GTEx_WB == F & GTEx_Brain == F) & G2S == T, type := 'G2S_only']
3102. temp[ (GTEx_WB == F & GTEx_Brain == T) & G2S == F, type := 'GTEx brain only']
3103. temp[ (GTEx_WB == F & GTEx_Brain == T) & G2S == T, type := 'GTEx_brain + G2S']
3104. temp[ (GTEx_WB == T & GTEx_Brain == F) & G2S == F, type := 'GTEx_WB only']
3105. temp[ (GTEx_WB == T & GTEx_Brain == F) & G2S == T, type := 'GTEx_WB + G2S']
3106. temp[ (GTEx_WB == T & GTEx_Brain == T) & G2S == F, type := 'GTEx_WB + GTEx_brain']
3107. temp[ (GTEx_WB == T & GTEx_Brain == T) & G2S == T, type := 'all']
3109. newtemp <- temp[, .(.N), by = type]
3110. newtemp[, gwas_phenotype := dz]
3112. mod_dz_tabl <- rbind(mod_dz_tabl, newtemp)
3114. }
3117. mod_typs <- ggplot(mod_dz_tabl, aes(x = type, y = N, fill = type)) +
3118. geom_bar(stat = 'identity') +
3119. facet_wrap(gwas_phenotype~., nrow = 3) +
3120. theme(axis.text.x = element_text(angle = 60, hjust = 1))
3122. pdf(file = 'Output/revisions/mod_types_brain.pdf', height = 12, width = 10)
3123. print(mod_typs)
3124. dev.off()
3125. print(mod_typs)
3132. mod_typs2 <- ggplot(mod_dz_tabl[type %in% c('GTEx_brain + G2S', 'GTEx brain only')], aes(x = type, y = N, fill = type)) +
3133. geom_bar(stat = 'identity') +
3134. facet_wrap(gwas_phenotype~., nrow = 3, scales = 'free_y') +
3135. theme(axis.text.x = element_text(angle = 60, hjust = 1))
3137. pdf(file = 'Output/revisions/mod_types_brain2.pdf', height = 10, width = 6)
3138. print(mod_typs2)
3139. dev.off()
3140. print(mod_typs2)
3143. newmod <- mod_dz_tabl[type %in% c('GTEx_brain + G2S', 'GTEx brain only')]
3145. sums <- newmod[, .(sum = sum(N)), by = 'gwas_phenotype']
3147. newmod1 <- setDT(merge(newmod, sums, by = 'gwas_phenotype'))
3148. newmod1[, percentage := N / sum]
3152. mod_typs3 <- ggplot(newmod1, aes(x = gwas_phenotype, y = percentage, fill = type, color = type)) +
3153. geom_bar(stat = 'identity', position = 'stack') +
3154. geom_label(data = newmod1[type == 'GTEx_brain + G2S',], aes(x = gwas_phenotype, y = percentage - 0.06, label = round(percentage, 3)), fill = 'white', show.legend = F)+
3155. #facet_wrap(gwas_phenotype~., nrow = 3) +
3156. theme_minimal() +
3157. theme(axis.text.x = element_text(angle = 45, hjust = 1))
3159. pdf(file = 'Output/revisions/mod_types_brain3.pdf', height = 4, width = 7)
3160. print(mod_typs3)
3161. dev.off()
3162. print(mod_typs3)
3165. ```
3168. # Revisions Round 2
3171. ### Cell-type gene expression analysis using DEGs
3173. ```{r}
3174. library(data.table)
3175. library(ggplot2)
3178. dt_ <- fread('Data/trimmed_means.csv', header = T)
3181. dt <- dt_[!(rowSums(dt_[,-1] == 0)),]
3183. genes <- dt[[1]]
3184. expr_mat <- as.matrix(dt[, -1])
3185. rownames(expr_mat) <- genes
3186. cell_types <- colnames(expr_mat)
3188. # Row-wise z-score (across cell types, per gene)
3189. expr_scaled <- t(scale(t(expr_mat)))
3190. # scale() works column-wise, so transpose, scale, transpose back
3191. # Now each gene has mean=0, sd=1 across cell types
3193. # Build results table
3194. results <- rbindlist(lapply(cell_types, function(ct) {
3195. data.table(
3196. gene = genes,
3197. cell_type = ct,
3198. expr = expr_mat[, ct], # raw trimmed mean (log2 CPM)
3199. z_score = expr_scaled[, ct], # how extreme vs other cell types
3200. log2fc = expr_mat[, ct] - rowMeans(expr_mat[, cell_types != ct])
3201. # log2FC: already in log space so subtraction = fold change
3202. )
3203. }))
3205. # Convert z-score to p-value and correct for multiple testing
3206. results[, p_nominal := 2 * pnorm(-abs(z_score))]
3207. results[, p_adj := p.adjust(p_nominal, method = "BH")]
3209. # Filter to DEGs: significant + meaningfully upregulated
3210. degs <- results[
3211. p_adj < 0.05 & # statistically specific
3212. log2fc > 1 & # at least 2-fold above mean of other cell types
3213. expr > 1 # expressed at least log2(CPM) > 1 in this cell type
3214. ][order(cell_type, -z_score)]
3216. # How many DEGs per cell type?
3217. degs[, .N, by = 'cell_type'][order(-N)]
3218. ```
3221. ```{r}
3222. # write top findings to output file
3224. top_markers <- degs[, head(.SD, 10), by = cell_type, .SDcols = names(degs)[-2]]
3225. write.table(top_markers, 'Output/revisions2/cell/brain_degs.txt', col.names = T, row.names = F, quote = F)
3226. ```
3228. ### Cell type enrichment analysis
3229. ```{r}
3230. sdt_ <- fread('Data/all_exc_BHsig_assocs_PGCancestryTWAS.txt', header = TRUE)
3231. gwas <- c('PGC_AA_MDD_2023', 'PGC_AA_PTSD_2024','PGC_AA_SCZ_2022','PGC_ASIAN_SCZ_2022','PGC_EUR_PTSD_2024','PGC_HNA_PTSD_2024', 'PGC_EURO_SCZ_2022')
3232. mod <- c('PrediXcan_Brain_Cortex')
3234. sdt <- sdt_[!(gwas_phenotype %in% gwas | training_model %in% mod),]
3237. overlap_b <- merge(
3238. sdt[, .(`gene_name`, gwas_phenotype)],
3239. top_markers[, .(`gene`, `cell_type`)],
3240. by.x = "gene_name",
3241. by.y = "gene",
3242. allow.cartesian = TRUE
3243. )
3245. # Count overlaps for each disease × cell_type pair
3246. result_b <- overlap_b[, .(n_genes = uniqueN(gene_name), genes = paste(unique(gene_name), collapse = '|')), by = c('gwas_phenotype', 'cell_type')]
3249. write.table(result_b, 'Output/revisions2/cell/disease_cell_types_brain.txt', col.names = T, row.names = F, quote = F, sep = '\t')
3250. ```
3251. ### Cell type overlap with DEGs from brain
3253. ```{r}
3254. dt <- copy(sdt)
3256. pg <- fread('Data/PanglaoDB_markers_27_Mar_2020.tsv', header = T)
3259. overlap <- merge(
3260. dt[, .(`gene_name`, gwas_phenotype)],
3261. pg[, .(`official gene symbol`, `cell type`)],
3262. by.x = "gene_name",
3263. by.y = "official gene symbol",
3264. allow.cartesian = TRUE
3265. )
3267. # Count overlaps for each disease × cell_type pair
3268. result <- overlap[, .(n_genes = uniqueN(gene_name), genes = paste(unique(gene_name), collapse = '|')), by = c('gwas_phenotype', 'cell type')]
3270. duo <- intersect(overlap$gene_name, overlap_b$gene_name)
3272. combo <- result[sapply(strsplit(genes, "\\|"), function(g) any(g %in% duo))]
3274. write.table(combo, 'Output/revisions2/cell/disease_cell_types_brain_DuoDEGs.txt', col.names = T, row.names = F, quote = F, sep = '\t')
3276. combo
3277. ```
3280. # Genetic Ancestry and SNP predictors
3283. First I want to get all the SNPs used by GTEx and G2S and annotate them with population specific allele frequencies.
3284. ```{r}
3286. otpt <- 'Output/revisions2/allelefreq/'
3287. dir.create(otpt)
3289. # read in 1000 genomes allele frequency data and SNP annotations
3292. rsids <- fread('Data/all_rsid_locs_af_anc.annovar.hg38_multianno.txt', header = T)
3293. anno_ <- fread('Data/all_snplist_avsnp151_convert.avinput', header = F)
3295. anno <- setNames(anno_, c('Chr', 'Start', 'End', 'Ref', 'Alt', 'rsid'))
3296. freq <- setDT(merge(anno, rsids, by = c('Chr', 'Start', 'End', 'Ref', 'Alt')))
3298. # Read in all SNP data from GTEx and G2S models
3299. all_snps <- fread('Data/all_snpdata.txt', header = T)
3301. # remove multi-allelic SNPs for now.
3302. multi_allelic <- freq[, .(.N), by = 'rsid'][N > 1]$rsid
3305. #annotate GTEx and G2S SNPs and rename columns
3306. dt_ <- setDT(merge(all_snps, freq[!(rsid %in% multi_allelic),], by = 'rsid'))
3308. setnames(dt_, c('gnomad41_genome_AF', 'gnomad41_genome_AF_afr', 'gnomad41_genome_AF_ami', 'gnomad41_genome_AF_amr', 'gnomad41_genome_AF_asj', 'gnomad41_genome_AF_eas', 'gnomad41_genome_AF_fin', 'gnomad41_genome_AF_mid', 'gnomad41_genome_AF_nfe', 'gnomad41_genome_AF_sas'),
3309. c('ALL', 'AFR', 'AMISH', 'ADM_AMR', 'ASHKENAZI', 'E_ASIAN', 'FINN', 'MIDEAST', 'EUR', 'S_ASIAN'))
3311. cols <- c(names(all_snps),c('Ref', 'Alt', 'ALL', 'AFR', 'AMISH', 'ADM_AMR', 'ASHKENAZI', 'E_ASIAN', 'FINN', 'MIDEAST', 'EUR', 'S_ASIAN'))
3312. dt <- dt_[, ..cols]
3314. dt <- dt %>% mutate_at(c('ALL', 'AFR', 'AMISH', 'ADM_AMR', 'ASHKENAZI', 'E_ASIAN', 'FINN', 'MIDEAST', 'EUR', 'S_ASIAN'), as.numeric)
3317. # get magnitude of difference between African and european allele frequency.
3318. dt <- dt[, afr_eu_dif := abs(dt$AFR - dt$EUR)][!is.na(afr_eu_dif),]
3321. # first, get weights that are not shared and assess the AF difference between the two groups
3322. dt[shared == TRUE, type := 'shared']
3323. dt[shared == FALSE & !is.na(weight.G2S), type := 'G2S']
3324. dt[shared == FALSE & !is.na(weight.GTEx), type := 'GTEx']
3326. #SNP data table created
3328. ```
3331. next I want to add the weights for each SNP-gene combination.
3333. ATTN: You will have to download the following models from the kachuri et al zenodo page located at the following address:
3334. - AA.cis-eQTL.tar.gz
3335. -after unzipping the tarball, the file for AA will be ._gala.sage.AA_Whole_Blood_tw_0.5_signif.db. I renamed this to AA_Whole_Blood.db and put it in [KACHURI_MODELS_DIR]
3336. - MX.cis-eQTL.tar.gz
3337. - PR.cis-eQTL.tar.gz
3338. - All.cis-eQTL.tar.gz
3339. [https://zenodo.org/records/7735723]
3341. The gtex model can be found here
3342. - PrediXcan_Whole_Blood.db
3343. [https://zenodo.org/records/3842289]
3345. ```{r}
3346. # Extract weights from all SNPs in the models of interest
3347. aamod <- '[KACHURI_MODELS_DIR]/AA_Whole_Blood.db'
3348. allmod <- '[KACHURI_MODELS_DIR]/All_Whole_Blood.db'
3349. mamod <- '[KACHURI_MODELS_DIR]/MX_Whole_Blood.db'
3350. prmod <- '[KACHURI_MODELS_DIR]/PR_Whole_Blood.db'
3351. gtmod <- '[GTEx_MODELS_DIR]/db/PrediXcan_Whole_Blood.db'
3353. mods <- data.table(modpath = c(aamod, allmod, mamod, prmod),
3354. modnm = c('G2S_aa', 'G2S_all', 'G2S_mx', 'G2S_pr'),
3355. modlongnm = c('AA_Whole_Blood', 'All_Whole_Blood', 'MX_Whole_Blood', 'PR_Whole_Blood'))
3357. gtx <- dbConnect(RSQLite::SQLite(), gtmod)
3358. gtx_weight <- setDT(dbReadTable(gtx, 'weights'))
3359. dbDisconnect(gtx)
3361. g2s <- data.table()
3363. for(modp in mods$modpath) {
3364. # first read in the model data
3365. mod <- dbConnect(RSQLite::SQLite(), modp)
3366. mod_weight <- setDT(dbReadTable(mod, 'weights'))
3367. nm <- mods[modpath == modp,]$modlongnm
3368. mod_weight$training_model <- nm
3369. mod_weight$gene <- str_remove(mod_weight$gene, '\\..+')
3370. g2s <- rbind(g2s, mod_weight)
3371. dbDisconnect(mod)
3372. }
3375. gtx_weight$training_model <- 'GTEx'
3376. ```
3378. Now having gotten the weights and the SNPs, I want to harmonize the SNPs for the sake of merging. the ref and alt alleles may not be the same across G2S and GTEx. When this is the case, the weight is inaccurate. I need to find the matching and flipped alleles and then adjust the weight accordingly.
3381. ```{r}
3382. #check rsid matching gtex vs g2s
3384. g2s_rsid <- unique(g2s[, c('rsid', 'ref_allele','eff_allele')])
3385. gtx_rsid <- unique(gtx_weight[, c('rsid', 'ref_allele','eff_allele')])
3387. dim(g2s_rsid)
3388. length(unique(g2s_rsid$rsid))
3389. # one rsid per rsid/ref/alt grouping
3391. dim(gtx_rsid)
3392. length(unique(gtx_rsid$rsid))
3393. # same for gtex
3395. length(intersect(gtx_rsid$rsid, g2s_rsid$rsid))
3396. #52,284 shared rsids
3398. #merge across rsid, ref, and alt,
3399. rsid_both <- setDT(merge(gtx_rsid, g2s_rsid, by = c('rsid', 'ref_allele', 'eff_allele')))
3400. rsid2 <- intersect(gtx_rsid$rsid, g2s_rsid$rsid)
3402. discordant <- setDT(rbind(gtx_rsid, g2s_rsid))[rsid %in% rsid2 & !(rsid %in% rsid_both$rsid),][order(rsid)]
3404. # get all aligned SNP weights
3405. all_rsids <- setDT(rbind(gtx_rsid, g2s_rsid))[ !(rsid %in% discordant$rsid),]
3407. ## Harmonize SNPs across GTEx and G2S
3408. # Subset the used rsids by those that have consistent ref and eff alleles
3409. # then get the frequencies relying on both rsid, ref, and alt
3411. dt_same_ <- setDT(merge(dt, all_rsids, by.x = c('rsid', 'Ref', 'Alt'), by.y = c('rsid' ,'ref_allele', 'eff_allele')))
3412. dt_same_$allele_order <- 'match'
3414. dt_dif_ <- setDT(merge(dt, all_rsids, by.x = c('rsid', 'Alt', 'Ref'), by.y = c('rsid' ,'ref_allele', 'eff_allele')))
3415. dt_dif_$allele_order <- 'flip'
3417. # adjust allele frequencies (assuming biallelic SNPs)
3418. af_cols <- c('ALL', 'AFR', 'AMISH', 'ADM_AMR', 'ASHKENAZI', 'E_ASIAN', 'FINN', 'MIDEAST', 'EUR', 'S_ASIAN')
3419. dt_dif_[ , (af_cols) := lapply(.SD, function(x) 1-x), .SDcols = af_cols]
3421. # adjust weights for opposite direction of effect
3422. dt_dif_$weight.G2S <- -1 * dt_dif_$weight.G2S
3423. dt_dif_$weight.GTEx <- -1 * dt_dif_$weight.GTEx
3425. # join all
3426. dt2 <- unique(setDT(rbind(dt_same_, dt_dif_)))
3428. ## Now we can treat the allele frequencies as true values rather than just the differences
3429. dt2a <- dt2[model == 'All_Whole_Blood',]
3431. dt2a
3432. ```
3435. Now that I have harmonized and annotated the SNPs of interest with allele frequencies, I want to begin characterizing the genes according to aggregate statistics derived from the SNP predictors.
3436. ```{r}
3437. # Read in all eQTL data
3438. eqtldat <- fread('Data/eQTL_data_all.txt', header = T)
3439. eqtldat[is.na(weight.G2S), snptype := 'GTEx_SNP']
3440. eqtldat[is.na(weight.GTEx), snptype := 'G2S_SNP']
3441. eqtldat[is.na(snptype), snptype := 'shared_SNP']
3444. edat_type <- eqtldat[Significant == 'Significant', .(.N, N_genes = length(unique(gene))), by = c('gwas_phenotype', 'snptype', 'category')][order(gwas_phenotype)]
3446. # plot each gene. data is percentage of SNPs that are G2S specific
3447. eqtldat[snptype == 'GTEx_SNP', snptype_q := -1]
3448. eqtldat[snptype == 'shared_SNP', snptype_q := 0]
3449. eqtldat[snptype == 'G2S_SNP', snptype_q := 1]
3451. edat_prop <- eqtldat[ Significant == 'Significant', .(type_quant = sum(snptype_q)), by = c('gwas_phenotype', 'category', 'gene')]
3454. eqtldat[, category_collapsed := fct_collapse(category,
3455. "G2S" = c("G2S only", "G2S co-tested"),
3456. "GTEX" = c("GTEX WB co-tested"),
3457. "Shared" = c("co-significant")
3458. )]
3460. gene_features <- eqtldat[, .(
3462. # SNP counts
3463. n_snps = .N,
3465. # model-specific SNP counts
3466. n_g2s_snps = sum(!is.na(weight.G2S)),
3467. n_gtex_snps = sum(!is.na(weight.GTEx)),
3469. # proportions
3470. prop_g2s_snps = mean(!is.na(weight.G2S)),
3471. prop_gtex_snps = mean(!is.na(weight.GTEx)),
3473. # weights
3474. mean_g2s_weight_magnitude = mean(abs(weight.G2S), na.rm = TRUE),
3475. mean_gtex_weight_magnitude = mean(abs(weight.GTEx), na.rm = TRUE),
3477. # MAF
3478. mean_afr_maf = mean(AFR, na.rm = TRUE),
3479. mean_eur_maf = mean(EUR, na.rm = TRUE),
3480. mean_afr_eur_maf_diff = mean(abs(AFR -EUR), na.rm = TRUE),
3481. n_maf_diff_over_0.3 = sum(AFR-EUR > 0.3, na.rm = TRUE), #total number of SNPs with maf dif > 0.3
3482. mean_afr_maf_weighted = mean(AFR*abs(weight.G2S), na.rm = TRUE),
3483. mean_eur_maf_weighted = mean(EUR*abs(weight.GTEx), na.rm = TRUE),
3484. mean_maf_diff_weighted = mean(AFR*abs(weight.G2S), na.rm = TRUE)-mean(EUR*abs(weight.GTEx), na.rm = TRUE),
3486. # SNP type composition
3487. prop_common = mean(SNP_type == "COMMON"),
3488. prop_rare = mean(SNP_type == "RARE"),
3489. prop_afr_specific = mean(SNP_type == "AFR"),
3490. prop_eur_specific = mean(SNP_type == "EUR"),
3491. n_afr_specific = sum(SNP_type == "AFR"),
3492. n_eur_specific = sum(SNP_type == "EUR")
3494. ), by = .(gene, gwas_phenotype, category_collapsed)]
3497. dt <- gene_features[gwas_phenotype == "PGC_all_SCZ_2022",]
3498. dt$category_collapsed <- as.factor(dt$category_collapsed)
3500. dt[is.nan(mean_gtex_weight_magnitude), mean_gtex_weight_magnitude := 0]
3501. dt[is.nan(mean_g2s_weight_magnitude), mean_g2s_weight_magnitude := 0]
3502. ```
3504. Now we have data on the SNP characteristics for each gene involved in gene-trait associations. Next I want to figure out how the SNP characteristics differ for genes that were significantly associated with Schizophrenia in different models.
3505. ```{r}
3507. # Melt to long format for plotting
3508. features_of_interest <- names(dt)[-c(1:3)]
3510. gf_long <- melt(dt[category_collapsed != "None"],
3511. id.vars = "category_collapsed",
3512. measure.vars = features_of_interest)
3514. # Violin plots - distribution of each feature by category
3515. ff <- ggplot(gf_long, aes(x = category_collapsed, y = value,
3516. fill = category_collapsed)) +
3517. geom_violin(alpha = 0.7) +
3518. geom_boxplot(width = 0.1, outlier.shape = NA) +
3519. facet_wrap(~variable, scales = "free_y") +
3520. theme_minimal() +
3521. labs(title = "Feature distributions by association category",
3522. x = NULL, y = NULL) +
3523. theme(legend.position = "none",
3524. axis.text.x = element_text(angle = 45, hjust = 1))
3526. name = 'Output/revisions2/all_measures_SCZ.png'
3527. png(filename = name, width = 11, height = 6, units = "in", pointsize = 12, bg = "white", type = 'cairo', res = 600)
3528. print(ff)
3529. dev.off()
3531. print(ff)
3532. ```
3534. ### Fig S4: Odds ratios for G2S vs GTEx significance
3536. ```{r}
3537. library(coin) # for permutation tests - robust with small groups
3539. # The key scientific contrasts
3540. # G2S only vs None: what makes a gene significant in G2S?
3541. # G2S only vs Shared: what makes a gene G2S-exclusive vs both?
3542. # GTEX only vs Shared: what makes a gene GTEx-exclusive?
3544. # For each feature, test G2S vs GTEX vs Shared
3545. kruskal_results <- dt[category_collapsed != "None",
3546. lapply(.SD, function(x) kruskal.test(x ~ category_collapsed)$p.value),
3547. .SDcols = features_of_interest]
3549. # Tidy it up
3550. kruskal_long <- setDT(melt(kruskal_results, variable.name = "feature",
3551. value.name = "p_value"))
3552. kruskal_long[, p_adj := p.adjust(p_value, method = "BH")]
3553. kruskal_long[, `-log10 p_adj` := -1*log10(p_adj)]
3554. kruskal_long[order(p_adj)]
3555. kruskal_long[, Significance := fifelse(p_adj <= 0.05, "Sig", "Insig")]
3557. f1 <- ggplot(data = kruskal_long, aes(x = feature, y = `-log10 p_adj`, color = Significance, shape = Significance)) +
3558. geom_segment( aes(x=feature, xend=feature, y=0, yend=`-log10 p_adj`), color="grey") +
3559. geom_point(size=4) +
3560. theme_minimal() +
3561. coord_flip()
3563. name = 'Output/revisions2/all_measures_SCZ_kruskal.png'
3564. png(filename = name, width = 6, height = 6, units = "in", pointsize = 12, bg = "white", type = 'cairo', res = 600)
3565. print(f1)
3566. dev.off()
3568. print(f1)
3569. ```
3571. ### Fig S3: Print significantly different features with DOE
3573. ```{r}
3575. # Get significant features first
3576. sig_features <- kruskal_long[p_adj <= 0.05, as.character(feature)]
3578. dtz <- copy(dt)
3579. dtz[, (sig_features) := lapply(.SD, scale), .SDcols =sig_features]
3581. # Pairwise Wilcoxon for each significant feature
3582. # The three scientifically meaningful contrasts
3583. contrasts <- list(
3584. c("G2S", "Shared"), # G2S-exclusive vs both models
3585. c("GTEX", "Shared"), # GTEX-exclusive vs both models
3586. c("G2S", "GTEX") # the key contrast you care about most
3587. )
3589. direction_results <- rbindlist(lapply(sig_features, function(feat) {
3590. rbindlist(lapply(contrasts, function(pair) {
3592. x <- dtz[category_collapsed == pair[1], get(feat)]
3593. y <- dtz[category_collapsed == pair[2], get(feat)]
3595. wt <- wilcox.test(x, y, conf.int = TRUE)
3597. data.table(
3598. feature = feat,
3599. contrast = paste(pair[1], "vs", pair[2]),
3600. median_A = median(x, na.rm = TRUE),
3601. median_B = median(y, na.rm = TRUE),
3602. difference = median(x, na.rm = TRUE) - median(y, na.rm = TRUE),
3603. direction = fifelse(median(x, na.rm=TRUE) > median(y, na.rm=TRUE),
3604. paste("higher in", pair[1]),
3605. paste("higher in", pair[2])),
3606. p_value = wt$p.value,
3607. p_adj = p.adjust(wt$p.value, method = "BH")
3608. )
3609. }))
3610. }))
3612. # View ordered by contrast then significance
3613. direction_results[order(contrast, p_adj)]
3614. ```
3616. ### Fig S3: DOE and effect size of feature regression
3617. ```{r}
3618. library(viridis)
3620. # Dot plot: effect size on x, feature on y, faceted by contrast
3621. f3 <- ggplot(direction_results[p_adj <= 0.05],
3622. aes(color = abs(difference),
3623. y = reorder(feature, abs(difference)),
3624. shape = direction,
3625. x = -log10(p_adj))) +
3626. geom_vline(xintercept = 0, linetype = "dashed", color = "gray60") +
3627. geom_point(alpha = 1, size = 3) +
3628. facet_wrap(~contrast, ncol = 3) +
3629. scale_color_viridis_c(option = "viridis") +
3630. #scale_color_manual(
3631. # values = c("#3B8BD4", "#1D9E75"),
3632. # labels = c("Higher in B", "Higher in A")
3633. #) +
3634. theme_minimal() +
3635. labs(#title = "Direction of effect for significant SNP features",
3636. #subtitle = "Point size = -log10 adjusted p-value",
3637. color = "Median difference abs(A - B)",
3638. y = NULL,
3639. color = NULL,
3640. x = "-log10 p-adj",
3641. shape = "Model impact")
3644. name = 'Output/revisions2/all_measures_SCZ_wilcoxon.png'
3645. png(filename = name, width = 9, height = 5, units = "in", pointsize = 12, bg = "white", type = 'cairo', res = 600)
3646. print(f3)
3647. dev.off()
3649. print(f3)
3650. ```
3653. ```{r}
3654. #now look specifically for the effects when binning for coverage
3655. # Among genes where GTEx coverage is EQUAL (n_gtex_snps similar),
3656. # does prop_afr_specific still differ between G2S and GTEX?
3658. # Bin by prop_g2s-prop_gtex to control for coverage
3659. dt[category_collapsed != "None",
3660. model_coverage := cut(abs(prop_g2s_snps-prop_gtex_snps),
3661. breaks = quantile(abs(prop_g2s_snps-prop_gtex_snps),
3662. probs = seq(0,1,0.25),
3663. na.rm = TRUE),
3664. include.lowest = TRUE)]
3666. #dt[ ,':='(model_coverage=NULL)]
3668. # Now test mean_maf_diff_weighted within coverage-matched genes
3669. dt[category_collapsed %in% c("G2S", "GTEX") &
3670. model_coverage == levels(model_coverage)[1], # mindist coverage bin
3671. wilcox.test(mean_maf_diff_weighted ~ category_collapsed)]
3672. ```
3674. ```{r}
3675. dt[, g2s_binary := fifelse(
3676. category_collapsed %in% c("G2S", "Shared"), 1, 0)]
3678. dt[, g2s_binary := fifelse(
3679. category_collapsed %in% c("G2S"), 1, 0)]
3681. model_logistic <- glm(
3682. g2s_binary ~ n_snps + n_g2s_snps + n_gtex_snps +
3683. prop_g2s_snps + prop_gtex_snps + mean_g2s_weight_magnitude +
3684. mean_gtex_weight_magnitude + n_maf_diff_over_0.3 +
3685. mean_maf_diff_weighted + prop_common + prop_rare +
3686. n_eur_specific,
3687. data = dt[category_collapsed != "None"],
3688. family = binomial
3689. )
3691. summary(model_logistic)
3692. exp(coef(model_logistic)) # odds ratios
3693. exp(confint(model_logistic)) # confidence intervals
3696. summary(model_logistic)
3697. ```
3699. ## Model-wide FDR thresholding
3701. ```{r}
3702. dt_ <- fread('Data/all_exc_assocs_PGCancestryTWAS.txt', header =T)
3704. mods <- c("AA_Whole_Blood","All_Whole_Blood","MX_Whole_Blood","PrediXcan_Whole_Blood","PR_Whole_Blood")
3705. test_phens <- c("PGC_allADHD_2022", "PGC_allMDD_2023","PGC_ALL_PTSD_2024","PGC_all_SCZ_2022","PGC_BD1_2021","SUD_alc_2019")
3708. dt <- dt_[gwas_phenotype %in% test_phens & training_model %in% mods, -c('bhpval','bfpval')]
3710. dt$bfpval <- p.adjust(dt$pvalue, "bonferroni")
3711. dt$bhpval <- p.adjust(dt$pvalue, "BH")
3713. dtbh <- dt[bhpval <= 0.05,]
3714. dtbh$gg <- paste(dtbh$gene, dtbh$gwas_phenotype, sep = '_')
3716. dtbh_shared <- dtbh[training_model != "PrediXcan_Whole_Blood" & gg %in% dtbh[training_model == "PrediXcan_Whole_Blood",]$gg,]
3717. dtbh[gg %in% dtbh_shared$gg, status := 'shared']
3718. dtbh_g2s <- dtbh[training_model != "PrediXcan_Whole_Blood" & !(gg %in% dtbh[training_model == "PrediXcan_Whole_Blood",]$gg),]
3719. dtbh[gg %in% dtbh_g2s$gg, status := 'G2S only']
3720. dtbh_gtex <- dtbh[training_model == "PrediXcan_Whole_Blood" & !(gg %in% dtbh_shared$gg),]
3721. dtbh[gg %in% dtbh_gtex$gg, status := 'GTEx only']
3723. dtbhn <- dtbh[, .(.N), by = c('gwas_phenotype', 'training_model')]
3724. dtbhn[training_model == 'All_Whole_Blood', `gene expression model` := 'G2S All WB']
3725. dtbhn[training_model == 'AA_Whole_Blood', `gene expression model` := 'G2S African American WB']
3726. dtbhn[training_model == 'MX_Whole_Blood', `gene expression model` := 'G2S Mexican American WB']
3727. dtbhn[training_model == 'PR_Whole_Blood', `gene expression model` := 'G2S Puerto Rican WB']
3728. dtbhn[training_model == 'PrediXcan_Whole_Blood', `gene expression model` := 'GTEx European WB']
3730. dtbhn
3731. ```

PGC_master_analysis2.Rmd, under CC-BY-4.0 · at the source