Running univariate and bivariate LAVA
Aaron Wagen & Regina H. Reynolds
UCLAim: run univariate and bivariate tests
1 Background
Neurodegenerative diseases are often defined by the predominant neuropathology (e.g. synucleinopathies by \(\alpha\)-synuclein deposition), but more often than not, they also present with co-pathologies, suggesting that they might share common pathogenic pathways (PMID: 29878075, 30258235). This notion is supported by genetic studies, which have (i) identified shared risk loci across neurodegenerative disorders, such as APOE and BIN1 in Alzheimer’s disease (AD) and Lewy body dementia (LBD), or GBA, SNCA, TMEM175 in Parkinson’s disease (PD) and LBD and (ii) demonstrated that genetic risk scores derived from one neurodegenerative disorder can predict risk of another, as with AD and PD predicting risk of LBD (PMID: 31701892, 30953760, 33589841).
From a clinical perspective, neurodegenerative diseases are often also defined in terms of their predominant cognitive and/or motor symptoms (e.g. Alzheimer's disease by cognitive and memory impairment or Parkinson's disease by parkinsonism), but in reality, present as highly heterogenous diseases, with symptoms spanning multiple domains, including neurospychiatric symptoms (PMID: 27188934, 28332488. Indeed, a higher prevalence of depression has been observed in individuals with dementia compared to those without dementia (PMID: 30536144). Notably, a similar (albeit reversed) phenomenon has been observed in some neuropsychiatric disorders, with a higher risk of dementia diagnoses observed in individuals with schizophrenia (SCZ) versus individuals without a history of serious mental illness (PMID: 33688938, 26444987).
These clinical and/or pathological overlaps are not, however, always reflected in terms of global genetic correlations, with a recent study of global genetic correlation between neurological phenotypes demonstrating limited overlap between individual neurodegenerative disorders as well as between neurodegenerative and neuropsychiatric disorders (PMID: 29930110). One explanation for this is that global studies only consider an average across the entire genome, and thus, may fail to detect correlations confined to particular genomic regions, or where changes are balanced out over the genome. This limitation has recently been addressed with the development of Local analysis of [co]variant annotation (LAVA), which is able to evaluate local heritability over multiple traits of interest (using GWAS summary statistics) and detect local regions of shared genetic association (https://www.biorxiv.org/content/10.1101/2020.12.31.424652v2). Here, we apply LAVA to genome-wide association studies (GWASs) derived from three neurodenerative disorders (AD, PMID: 30617256;LBD, PMID: 33589841; and PD, PMID: 31701892), and three neuropsychiatric disorders (bipolar disorder, BIP, PMID: 34002096; major depressive disorder, MDD, PMID: 30718901; and SCZ, PMID: 29483656).
2 Supplementary code
Following section includes any intermediary code used in this .Rmd.
2.1 Run univariate and bivariate tests
This can either be sourced directly:
source(here::here("scripts", "02_run_univar_bivar_test.R"))Alternatively, for many phenotypes/loci it can be worth running this script directly on the cluster, using nohup, as performed below:
# Have to navigate to root project folder for script to work (as it uses here package)
cd /home/rreynolds/misc_projects/neurodegen-psych-local-corr
nohup Rscript \
/home/rreynolds/misc_projects/neurodegen-psych-local-corr/scripts/02_run_univar_bivar_test.R \
&>/home/rreynolds/misc_projects/neurodegen-psych-local-corr/logs/02_run_univar_bivar_test.log&2.2 Tidying and loading results
Once run, results can be loaded using the following code chunk:
source(here::here("scripts", "02b_tidy_univar_bivar.R"))results <- readRDS(here::here("results", "02_univar_bivar_test", "results_summary.rds"))global <-
read_delim(
file = here::here("results", "01_input_prep", "ldsc_corr", "ldsc_correlations.txt"),
delim = "\t"
) 2.3 Filtering for pleiotropic genes
ref <- import("/data/references/ensembl/gtf_gff3/v87/Homo_sapiens.GRCh37.87.gtf")
ref <- ref %>% keepSeqlevels(c(1:22), pruning.mode = "coarse")
ref <- ref[ref$type == "gene"]
# genes to check
genes <-
c("APOE", "TMEM175", "SNCA", "BIN1", "GRN", "KIF5A", "GBA")
# construct granges object with LD blocks for overlap with reference gtf
loci_gr <-
loci %>%
GenomicRanges::makeGRangesFromDataFrame(
.,
keep.extra.columns = TRUE,
ignore.strand = TRUE,
seqinfo = NULL,
seqnames.field = "CHR",
start.field = "START",
end.field = "STOP"
)
overlap <-
GenomicRanges::findOverlaps(loci_gr, ref) %>%
tibble::as_tibble()
loci_genes <-
tibble::tibble(
locus = loci_gr[overlap$queryHits]$LOC,
chr = loci_gr[overlap$queryHits] %>% GenomeInfoDb::seqnames() %>% as.character(),
locus_start = loci_gr[overlap$queryHits] %>% BiocGenerics::start(),
locus_end = loci_gr[overlap$queryHits] %>% BiocGenerics::end(),
gene_id = ref[overlap$subjectHits]$gene_id,
gene_name =
ref[overlap$subjectHits]$gene_name %>%
stringr::str_replace_all("-", "_") %>%
stringr::str_replace_all("\\.", "_"),
gene_biotype = ref[overlap$subjectHits]$gene_biotype %>% as.factor(),
gene_strand = ref[overlap$subjectHits] %>% BiocGenerics::strand() %>% as.character(),
gene_start = ref[overlap$subjectHits] %>% BiocGenerics::start(),
gene_end = ref[overlap$subjectHits] %>% BiocGenerics::end()
)
# filter df with ld block & gene info for genes of interest
loci_genes_filtered <-
loci_genes %>%
dplyr::filter(
gene_name %in% genes
)2.4 Comparison to previous data
source(here::here("scripts", "02c_run_comparison_gerring.R"))
source(here::here("scripts", "02c_run_comparison_smeland.R"))overlap_gerring <-
readRDS(here::here("results", "02_univar_bivar_test", "results_overlap_with_gerring.rds"))
overlap_smeland <-
readRDS(here::here("results", "02_univar_bivar_test", "results_overlap_with_smeland.rds"))2.5 Loading plotting functions
source(here::here("R", "plot_global_corr.R"))
source(here::here("R", "plot_chord_diagram.R"))
source(here::here("R", "plot_locus.R"))
source(here::here("R", "plot_edge_diagram.R"))3 Methods
3.1 Running univariate and bivariate tests
The detection of valid and interpretable local genetic corrlation (\(r_{g}\)) requires the presence of sufficient local genetic signal. For this reason, a univariate test was performed as a filtering step for bivariate local \(r_{g}\) analyses. Bivariate local correlations were only performed for pairs of traits which both exhibited a significant univariate local genetic signal (p < 0.05/300, where the denominator represents the total number of tested loci, i.e. LD blocks). This resulted in a total of 1603 bivariate tests spanning 275 distinct LD blocks.
3.2 Comparison to previous data
3.2.1 Gerring et al. 2022
We intersected the genomic coordinates for significant and non-significant local genetic associations (including only neuropsychiatric phenotypes, i.e. BIP, MDD and SCZ) from LAVA with existing results from a recent analysis of 10 psychiatric disorders and 10 substance abuse phenotypes, which also used LAVA (Gerring et al.; PMID: 35525699).
3.2.2 Smeland et al. 2021
We intersected the genomic coordinates for significant and non-significant local genetic associations between PD and SCZ with lead SNPs that were found to be jointly associated with PD and SCZ in a recent study using a conditional/conjunctional false discovery rate (FDR) approach to evaluate the overlap between the two traits. Lead SNPs from the study by Smeland et al. were converted from rsIDs to their GRCh37 genomic coordinates using the SNPlocs.Hsapiens.dbSNP144.GRCh37 package and overlapped with LD blocks using findOverlaps() from the GenomicRanges package.
4 Results
4.1 Tabular results
4.1.1 Global correlations
print("Global correlation column descriptions:")## [1] "Global correlation column descriptions:"tibble(
column = colnames(global),
description =
c(
"Phenotype 1",
"Phenotype 2",
"The estimated genetic correlation",
"The bootstrap standard error of the genetic correlation estimate",
"The bootstrap standard error of the genetic correlation estimate",
"P-value for genetic correlation",
"Estimated snp-heritability of the second phenotype ",
"Standard error of h2 for phenotype 2",
"Single-trait LD score regression intercept for phenotype 2",
"Standard error for single-trait LD score regression intercept for phenotype 2",
"Estimated genetic covariance between p1 and p2",
"Bootstrap standard error of cross-trait LD score regression intercept"
)
) %>%
DT::datatable(rownames = FALSE,
options = list(scrollX = TRUE),
class = 'white-space: nowrap')print("Global correlation results (p < 0.05/n_tests):")## [1] "Global correlation results (p < 0.05/n_tests):"combn <-
global$p1 %>%
unique() %>%
combn(m = 2) %>%
t() %>%
as_tibble() %>%
dplyr::rename(
p1 = V1,
p2 = V2
)
global_signif <-
global %>%
dplyr::filter(
!p1==p2
) %>%
dplyr::inner_join(combn) %>%
dplyr::filter(
p < 0.05/nrow(combn)
)
global_signif %>%
DT::datatable(rownames = FALSE,
options = list(scrollX = TRUE),
class = 'white-space: nowrap')4.1.2 Univariate results
print("Univariate column descriptions:")## [1] "Univariate column descriptions:"tibble(
column = colnames(results$univ),
description =
c(
"LD block ID",
"LD block chromosome",
"LD block start base pair",
"LD block end base pair",
"The number of SNPs within the LD block",
"The number of PCs retained within the LD block",
"Analysed phenotype",
"Observed local heritability",
"P-value from the univariate test (F-test for continuous, Chi-sq for binary)"
)
) %>%
DT::datatable(rownames = FALSE,
options = list(scrollX = TRUE),
class = 'white-space: nowrap')print("Univariate results (p < 0.05/n_ld_blocks):")## [1] "Univariate results (p < 0.05/n_ld_blocks):"results$univ %>%
dplyr::filter(p < 0.05/nrow(loci)) %>%
DT::datatable(rownames = FALSE,
options = list(scrollX = TRUE),
class = 'white-space: nowrap')4.1.3 Bivariate results
print("Bivariate column descriptions:")## [1] "Bivariate column descriptions:"tibble(
column = colnames(results$bivar),
description =
c(
"LD block ID",
"LD block chromosome",
"LD block start base pair",
"LD block end base pair",
"The number of SNPs within the LD block",
"The number of PCs retained within the LD block",
"Phenotype 1",
"Phenotype 2",
"Standardised coefficient for the local genetic correlation",
"Lower 95% confidence estimate for rho",
"Upper 95% confidence estimate for rho",
"Equivalent of taking the square of rho. Denotes the proportion of variance in genetic signal for phen1 explained by phen2 (and vice versa)",
"Lower 95% confidence estimate for r2",
"Upper 95% confidence estimate for r2",
"Simulation p-values for the local genetic correlation",
"FDR-corrected p-value"
)
) %>%
DT::datatable(rownames = FALSE,
options = list(scrollX = TRUE),
class = 'white-space: nowrap')print("Bivariate results (p < 0.05/n_bivar_tests):")## [1] "Bivariate results (p < 0.05/n_bivar_tests):"bivar_thres <- 0.05/nrow(results$bivar)
bivar_bonf <-
results$bivar %>%
dplyr::filter(p < bivar_thres)
bivar_bonf %>%
DT::datatable(rownames = FALSE,
options = list(scrollX = TRUE),
class = 'white-space: nowrap')4.2 Global correlations versus bivariate tests
4.2.1 Text
- Global cross-trait correlations revealed 4 significant global \(r_{g}\)'s (Figure 4.2). These global \(r_{g}\)'s are very much in line with what we expect from previous literature e.g. significant positive correlations between the three neuropsychiatric disorders (PMID: 29930110). Between neurodegenerative disorders, we only observe a significant postive \(r_{g}\) between LBD and PD.
- With a Bonferroni cut-off of 0.000031, we detect a total of 77 significant bivariate local \(r_{g}\)'s across 59 distinct LD blocks, of which 11 were associated with more than one trait pair. Notably, we detected local genetic correlations between pairs of phenotypes in the absence of a significant global genetic correlation (4.3).
- The phenotype with the highest number of significant bivariate local \(r_{g}\)'s was SCZ followed by BIP. Notably, the number of significant bivariate local \(r_{g}\)'s did not appear to correlate with GWAS sample size (Figure 4.4).
- As seen from the table below and Figure 4.5, the phenotype pairs with the most local \(r_{g}\)'s were (i) BIP and SCZ and (ii) AD and PD.
- For 30 local \(r_{g}\)'s, the 95% confidence intervals (CI’s) for the explained variance (r2) included 1 (Figure 4.6), which is consistent with a scenario where the genetic signal of those pairs of phenotypes in these LD blocks is completely shared.
- Within neuropsychiatric disorders, all significant bivariate local \(r_{g}\)'s were positive (matching the overall positive global \(r_{g}\)), while within neurodegenerative disorders there was a mix of positive and negative local \(r_{g}\)'s (Figure 4.5).
- All three neurodegenerative disorders shared one significant local \(r_{g}\) with SCZ, and also PD shared one negative local \(r_{g}\) with BIP.
4.2.2 Figures
data_to_plot <-
read_delim(
file = file.path(here::here("results", "01_input_prep"),
"input.info.txt"),
delim = "\t"
) %>%
dplyr::select(-filename) %>%
tidyr::pivot_longer(
cols = cases:controls,
names_to = "group",
values_to = "total"
) %>%
dplyr::mutate(
phenotype_group =
case_when(
phenotype %in% c("AD2019", "LBD2020", "PD2019.meta5.ex23andMe") ~ "Neurodegenerative",
TRUE ~ "Neuropsychiatric"
),
phenotype =
case_when(
phenotype == "AD2019" ~ "Alzheimer's disease\n(AD)",
phenotype == "LBD2020" ~ "Lewy body dementia\n(LBD)",
phenotype == "PD2019.meta5.ex23andMe" ~ "Parkinson's disease\n(PD)",
phenotype == "BIP2021" ~ "Bipolar disorder (BIP)",
phenotype == "MDD2019" ~ "Major depressive disorder\n(MDD)",
phenotype == "SCZ2018" ~ "Schizophrenia\n(SCZ)"
)
)
plot <-
data_to_plot %>%
dplyr::mutate(
n_text_col =
case_when(
total > 200009 ~ "high",
TRUE ~ "low"
),
group = stringr::str_to_title(group)
) %>%
ggplot(
aes(
x = group,
y =
fct_relevel(
phenotype,
c(data_to_plot %>% dplyr::arrange(phenotype_group, phenotype) %>% dplyr::pull(phenotype) %>% unique() %>% rev())
),
fill = total
)
) +
geom_tile(colour = "black") +
geom_text(
aes(
label = total %>% scales::comma(),
colour = n_text_col
),
size = 3
) +
coord_equal() +
scale_x_discrete(
position = "top"
) +
labs(x = "", y = "") +
scale_fill_distiller(palette = "Greys", direction = 1) +
scale_colour_manual(values = c("white", "black")) +
guides(colour = F, fill = F) +
theme_minimal() +
theme(
panel.grid = element_blank(),
)
plot
Figure 4.1: Heatmap of samples numbers in each GWAS.
plot_global_corr(
global_corr = global,
n_phenotypes = 6
)
Figure 4.2: Global genetic correlations between pairs of phenotypes. Significant negative and positive correlations are indicated by blue and red fill, respectively. Non-significant correlations (p >= 0.05/n_tests) have a grey fill.
plot_global_bivar(
global_corr =
global %>%
dplyr::mutate(
p1 = p1 %>%
stringr::str_replace_all("[:digit:]", "") %>%
stringr::str_remove("\\..*"),
p2 = p2 %>%
stringr::str_replace_all("[:digit:]", "") %>%
stringr::str_remove("\\..*")
),
bivar_corr =
results$bivar %>%
dplyr::mutate(
phen1 = phen1 %>%
stringr::str_replace_all("[:digit:]", "") %>%
stringr::str_remove("\\..*"),
phen2 = phen2 %>%
stringr::str_replace_all("[:digit:]", "") %>%
stringr::str_remove("\\..*")
),
n_phenotypes = 6
)
Figure 4.3: Comparison between the global genetic correlation estimated by LDSC (bottom) and the mean local genetic correlation from LAVA (top) across all tested LD blocks. Significant global correlations (p < 0.05/n_tests) are indicated with an asterisk ('*'). The number of significant local genetic correlations (p < 0.05/n_tests) is indicated by a number in each tile.
bivar_bonf %>%
dplyr::count(phen1) %>%
dplyr::bind_rows(
bivar_bonf %>%
dplyr::count(phen2) %>%
dplyr::rename(phen1 = phen2)
) %>%
dplyr::group_by(phen1) %>%
dplyr::summarise(
n_bivar = sum(n)
) %>%
dplyr::inner_join(
read_delim(
file = file.path(here::here("results", "01_input_prep"),
"input.info.txt"),
delim = "\t"
) %>%
dplyr::mutate(
phenotype = phenotype %>%
str_replace_all("[:digit:]", "") %>%
str_remove("\\..*"),
n_individ = cases + controls
),
by = c("phen1" = "phenotype")
) %>%
ggplot(
aes(
x = n_bivar,
y = n_individ,
colour = fct_relevel(
phen1,
fct_disease
)
)
) +
geom_point() +
scale_colour_brewer(
palette = "BrBG",
type = "div",
name = "GWAS"
) +
labs(
x = "Number of significant local genetic correlations",
y = "Sample size"
) +
theme_rhr
Figure 4.4: Plot of total number of significant bivariate local genetic correlations for each phenotype across all LD blocks versus GWAS sample size.
plot_bivar_chord_diagram(
bivar_corr = bivar_bonf,
fct_phen = fct_disease,
palette = RColorBrewer::brewer.pal(n = 6, name = "BrBG")
)
Figure 4.5: Chord diagram showing the number of significant bivariate local genetic correlations (p < 0.05/n_tests) between each of the phenotypes across all LD blocks. Positive and negative correlations are coloured red and blue, respectively.
n_upper_r2_ci <-
bivar_bonf %>%
dplyr::filter(
r2.upper == 1
) %>%
dplyr::count(phen1, phen2, name = "n_upper_r2_ci")
phen_count <-
bivar_bonf %>%
dplyr::count(phen1, phen2, name = "total_n") %>%
dplyr::left_join(n_upper_r2_ci) %>%
dplyr::mutate(
label = str_c(phen1, " vs ", phen2),
n_upper_r2_ci =
case_when(
is.na(n_upper_r2_ci) ~ 0,
TRUE ~ as.numeric(n_upper_r2_ci)
),
n_upper_r2_ci_not_one =
total_n - n_upper_r2_ci
) %>%
tidyr::pivot_longer(
cols = contains("n_upper_r2")
) %>%
dplyr::mutate(
name =
case_when(
name == "n_upper_r2_ci" ~ "TRUE",
TRUE ~ "FALSE"
)
)
phen_count %>%
ggplot(
aes(
x =
fct_reorder(
.f = label,
.x = total_n,
.fun = max,
.desc = T
),
y = value,
fill = fct_rev(name)
)
) +
geom_col(
colour = "black"
) +
scale_fill_grey() +
labs(
x = "",
y = "Number of significant\nlocal genetic correlations",
fill = "Does upper limit of r2 95% confidence\n interval include 1?"
) +
theme_rhr
Figure 4.6: Number of significant local genetic correlations between trait pairs. The fill of the bars indicates the number of significant LD blocks for which the 95% confidence interval included 1.
4.3 Signficant phenotype correlations at LD block level
4.3.1 Text
- The LD block with the highest number of significant bivariate local \(r_{g}\)'s, LD block 758 and 1719, were positioned on chromosome 4 and 11, respectively (Figure 4.7, 4.8).
- The locus on chromosome 4 contains three protein-coding genes (GPM6A, WDR17 and SPATA4), while the LD block on chromosome 11 contains ten protein-coding genes, one of which includes DRD2, which encodes dopamine receptor D2. As perhaps expected, all three neuropsychiatric disorders positively correlate with one another at LD block 1719. Somewhat unexpected is the significant negative correlation between LBD and SCZ at LD block 1719.
print("Top LD block when Bonferroni-correcting (p < 0.05/n_bivar_tests):")## [1] "Top LD block when Bonferroni-correcting (p < 0.05/n_bivar_tests):"bivar_bonf %>%
dplyr::count(locus, chr, start, stop, n_snps) %>%
dplyr::arrange(-n) %>%
dplyr::inner_join(loci %>%
janitor::clean_names() %>%
dplyr::rename(
locus = loc
)) %>%
dplyr::slice_max(order_by = n, n = 1)4.3.2 Figures
results$bivar %>%
dplyr::filter(locus %in% unique(bivar_bonf$locus)) %>%
dplyr::mutate(
rho_fill =
case_when(
p < bivar_thres ~ round(rho, 2)
)
) %>%
ggplot(
aes(
x = phen1,
y = phen2,
fill = rho_fill,
label = round(rho, 2)
)
) +
geom_tile(colour = "black") +
geom_text(
size = 2
) +
facet_wrap(vars(locus), ncol = 5) +
scale_fill_distiller(palette = "RdBu", direction = -1, na.value = "#cccccc", limits = c(-1, 1)) +
theme_rhr
Figure 4.7: Heatmaps showing the standardised coefficient for genetic correlation (rho) for each association within the LD block with significant bivariate local genetic correlations (p < 0.05/n_tests). Significant negative and positive correlations are indicated by blue and red fill, respectively. Non-significant correlations have a grey fill.
plot_edge_diagram(
bivar_corr =
bivar_bonf,
phen = fct_disease,
ncol = 3,
)
Figure 4.8: Edge diagrams for each LD block showing the standardised coefficient for genetic correlation (rho) for each significant bivariate local genetic correlations (p < 0.05/n_tests). Significant negative and positive correlations are indicated by blue and red colour, respectively.
loci_gr <-
bivar_bonf %>%
dplyr::filter(
phen1 %in% c("AD", "LBD", "PD"),
phen2 %in% c("AD", "LBD", "PD")
) %>%
dplyr::count(locus, chr, start, stop, n_snps) %>%
dplyr::arrange(locus) %>%
GenomicRanges::makeGRangesFromDataFrame(
.,
keep.extra.columns = TRUE,
ignore.strand = TRUE,
seqinfo = NULL,
seqnames.field = "chr",
start.field = "start",
end.field = "stop"
)
fig_list = vector(mode = "list", length = length(loci_gr))
for(i in 1:length(loci_gr)){
fig_list[[i]] <-
plot_locus(
locus_gr = loci_gr[i],
ref = ref
)
names(fig_list)[i] <- str_c("locus_", loci_gr[i]$locus)
}
fig_list## $locus_10
##
## $locus_335
##
## $locus_375
##
## $locus_472
##
## $locus_681
##
## $locus_711
##
## $locus_758
##
## $locus_1055
##
## $locus_1186
##
## $locus_1190
##
## $locus_1273
##
## $locus_1367
##
## $locus_1726
##
## $locus_1823
##
## $locus_1840
##
## $locus_2054
##
## $locus_2113
##
## $locus_2128
##
## $locus_2192
##
## $locus_2281
##
## $locus_2351
4.4 LD blocks with known neurodegenerative pleiotropic genes
4.4.1 Text
- Pleiotropic genes include:
- APOE (implicated in AD and LBD)
- TMEM175 (implicated in PD, LBD, and RBD)
- SNCA (implicated in PD, LBD and RBD)
- BIN1 (implicated in AD and LBD)
- GRN (implicated in AD, PD, FTD, and LBD)
- KIF5A (implicated in PD and ALS/FTD)
- GBA (implicated in PD and LBD)
- AD and PD had sufficient local genetic signal (i.e. univariate heritability) at all tested LD blocks. LBD had insufficient local heritability in the GRN-, KIF5A and TMEM175-containing LD blocks, and therefore was not carried forward in these loci for bivariate testing.
- Of the LD blocks that contained implicated genes, only the SNCA- and APOE-containing LD blocks contained local \(r_{g}\)s that passed the Bonferroni significance cut-off (p < 0.05/n_tests; Figure 4.9a). Using a more lenient FDR correction, local \(r_{g}\)s were observed across all LD blocks, except the GRN- and TMEM175-containing LD blocks (Figure 4.9b).
- BIN1-containing LD block 319
- Bonferroni: no local \(r_{g}\)s.
- FDR: positive local \(r_{g}\) between AD and LBD.
- GBA-containing LD block 113
- Bonferroni: no local \(r_{g}\)s.
- FDR: negative local \(r_{g}\) between MDD and PD and positive local \(r_{g}\) between MDD and BIP
- Bonferroni: no local \(r_{g}\)s.
- SNCA-containing LD block 681
- Bonferroni & FDR: negative local \(r_{g}\) between AD and PD.
- KIF5A-containing LD block 1794
- Bonferroni: no local \(r_{g}\)s.
- FDR: positive local \(r_{g}\)s between BIP and SCZ.
- Bonferroni: no local \(r_{g}\)s.
- APOE-containing LD block 2351
- Bonferroni: positive local \(r_{g}\) observed between AD and LBD and negative local \(r_{g}\) between PD and LBD.
- FDR: as above and a negative local \(r_{g}\) betyween AD and PD.
- BIN1-containing LD block 319
- In terms of expectations, APOE and BIN1 are what we would expect i.e. positive local \(r_{g}\) between AD and LBD, but some unexpected local \(r_{g}\)s too e.g. negative local \(r_{g}\) between AD and PD in APOE- and SNCA-containing LD blocks.
4.4.2 Tables
print("LD blocks containing pleiotropic genes:")## [1] "LD blocks containing pleiotropic genes:"loci_genes_filtered %>%
DT::datatable(rownames = FALSE,
options = list(scrollX = TRUE),
class = 'white-space: nowrap')print("Univariate results filtered for neurodegenerative phenotypes AND LD blocks containing pleiotropic genes. Significant when p < 0.05/n_ld_blocks.")## [1] "Univariate results filtered for neurodegenerative phenotypes AND LD blocks containing pleiotropic genes. Significant when p < 0.05/n_ld_blocks."results$univ %>%
dplyr::inner_join(
loci_genes_filtered %>%
dplyr::select(locus, gene_name),
by = c("locus")
) %>%
dplyr::filter(
phen %in% c("AD", "LBD", "PD")
) %>%
dplyr::mutate(
significant =
case_when(
p < 0.05/nrow(loci) ~ TRUE,
TRUE ~ FALSE
)
) %>%
dplyr::select(
locus, gene_name, significant, everything()
) %>%
DT::datatable(rownames = FALSE,
options = list(scrollX = TRUE),
class = 'white-space: nowrap')print("Bivariate results filtered for neurodegenerative phenotypes AND LD blocks containing pleiotropic genes.")## [1] "Bivariate results filtered for neurodegenerative phenotypes AND LD blocks containing pleiotropic genes."results$bivar %>%
dplyr::inner_join(
loci_genes_filtered %>%
dplyr::select(locus, gene_name),
by = c("locus")
) %>%
dplyr::filter(
phen1 %in% c("AD", "LBD", "PD") & phen2 %in% c("AD", "LBD", "PD")
) %>%
dplyr::mutate(
bonferroni =
case_when(
p < bivar_thres ~ TRUE,
TRUE ~ FALSE
),
fdr =
case_when(
fdr < 0.05 ~ TRUE,
TRUE ~ FALSE
)
) %>%
dplyr::select(
locus, gene_name, bonferroni, fdr, everything()
) %>%
DT::datatable(rownames = FALSE,
options = list(scrollX = TRUE),
class = 'white-space: nowrap')4.4.3 Figures
a <-
plot_edge_diagram(
bivar_corr =
results$bivar %>%
dplyr::inner_join(
loci_genes_filtered %>%
dplyr::select(locus, gene_name),
by = c("locus")
) %>%
dplyr::mutate(
locus = str_c(locus, ": contains ", gene_name)
),
p_threshold = bivar_thres,
multiple_corr = FALSE,
phen = fct_disease,
ncol = 3
)
b <-
plot_edge_diagram(
bivar_corr =
results$bivar %>%
dplyr::inner_join(
loci_genes_filtered %>%
dplyr::select(locus, gene_name),
by = c("locus")
) %>%
dplyr::mutate(
locus = str_c(locus, ": contains ", gene_name)
) %>%
dplyr::filter(fdr < 0.05),
multiple_corr = FALSE,
phen = fct_disease,
ncol = 3
) +
theme(legend.position = "none")
cowplot::plot_grid(
a,b,
ncol = 1,
rel_heights = c(1,2),
align = "v",
axis = "lr",
labels = letters[1:2]
)
Figure 4.9: Edge diagrams for each LD block showing local genetic correlations when corrected using (a) the stringent Bonferroni correction or (b) the lenient FDR correction. The standardised coefficient for genetic correlation (rho) is indicated for significant local correlations. Negative and positive correlations are indicated by blue and red colour, respectively.
4.5 Overlap with previous data
4.5.1 Text
4.5.1.1 Gerring et al. 2022
- We intersected significant and non-significant bivariate local \(r_{g}\)s from this study that included neuropsychiatric disease traits (n, local \(r_{g}\)s = 261; n, distinct LD blocks: 157; n, signficiant local \(r_{g}\)s = 51; n, distinct LD blocks with significant local \(r_{g}\)s = 40) with previously reported local \(r_{g}\)s from a study of 10 psychiatric disorders and 10 substance abuse phenotypes (PMID: 35525699).
- 8 LD blocks were found to contain significant local \(r_{g}\)s from both studies, 4 of which contained the same phenotype pair (LD block 457, chr3:36840137-38729767, BIP and SCZ; LD block 951, chr6:26396201-27261035, BIP and SCZ; LD block 951, chr6:26396201-27261035, MDD and SCZ; LD block 952, chr6:27261036-28666364, MDD and SCZ; LD block 2483, chr22:38718590-40378783, BIP and SCZ) (4.10).
- Note: a total of 16 local \(r_{g}\)s (across 14 LD blocks) that included only BIP, MDD or SCZ were detected in Gerring et al. -- of these, we were able to perform bivariate tests using the same phenotype pairs in 6 LD blocks, and 4 LD blocks contained 5 local \(r_{g}\)s.
- It is worth noting that there were some differences in the GWAS versions used: (i) we used a more recent BIP GWAS and (ii) while both studies use the same MDD GWAS, our study did not include 23andMe data.
4.5.1.2 Smeland et al. 2021
- We intersected LD blocks containing significant and non-significant bivariate local \(r_{g}\)s between PD and SCZ with SNPs found to be jointly associated with PD and SCZ in a recent study using a conditional/conjunctional false discovery rate (FDR) approach.
- 6 SNPs intersected with tested LD blocks, but 0 local \(r_{g}\)s within these blocks passed the threshold for significance (p < 0.05/1603). Notably, in those LD blocks containing SNPs from Smeland et al. where allelic effect directions were consistent, a positive rho was observed in local \(r_{g}\)s, and likewise, in those LD blocks containing SNPs from Smeland et al. where allelic effect directions were inconsistent, a negative rho was observed in local \(r_{g}\)s.
4.5.2 Table
print("Overlaps with Gerring et al. 2022:")## [1] "Overlaps with Gerring et al. 2022:"overlap_gerring %>%
DT::datatable(rownames = FALSE,
options = list(scrollX = TRUE),
class = 'white-space: nowrap')print("Overlaps with Smeland et al. 2021:")## [1] "Overlaps with Smeland et al. 2021:"overlap_smeland %>%
DT::datatable(rownames = FALSE,
options = list(scrollX = TRUE),
class = 'white-space: nowrap')4.5.3 Figures
plot_data <-
overlap_gerring %>%
dplyr::filter(overlapping_phen == 2) %>%
dplyr::select(
locus.rhr,
chr, start, stop,
phen1.rhr, phen2.rhr,
contains("rho"),
contains("r2"),
p.rhr,
p.gerring
) %>%
dplyr::mutate(
significant =
case_when(
p.rhr < bivar_thres & p.gerring < 0.05/24054 ~ "Both",
p.rhr >= bivar_thres & p.gerring < 0.05/24054 ~ "Gerring et. al"
)
) %>%
tidyr::pivot_longer(cols = rho.gerring:p.gerring, names_to = "feature") %>%
dplyr::mutate(
dataset =
case_when(
stringr::str_detect(feature, "rhr") ~ "rhr",
stringr::str_detect(feature, "gerring") ~ "gerring"
),
feature =
feature %>%
stringr::str_remove("\\.rhr|\\.gerring"),
feature =
case_when(
feature == "p" ~ "-log10(p)",
TRUE ~ feature
),
label =
str_c(locus.rhr, "\n", phen1.rhr, "+", phen2.rhr)
) %>%
tidyr::pivot_wider(
names_from = dataset,
values_from = value
) %>%
dplyr::mutate(
across(
.cols = rhr:gerring,
~ case_when(
feature == "-log10(p)" ~ -log10(.x),
TRUE ~ .x
)
)
)
a <-
plot_data %>%
dplyr::filter(
feature %in%
c("rho")
) %>%
ggplot(
aes(
x = rhr,
y = gerring
)
) +
geom_point(
size = 1.5,
alpha = 0.8,
shape = 21,
colour = "black",
aes(fill = significant)
) +
geom_abline(
intercept = 0,
linetype = "dashed",
colour = "grey"
) +
ggrepel::geom_label_repel(
aes(
label = label
),
size = 2,
alpha = 0.7,
min.segment.length = unit(0, 'lines')
) +
facet_wrap(vars(feature)) +
expand_limits(x = 0, y = 0) +
coord_equal() +
labs(
x = "This study",
y = "Gerring et al. 2022"
) +
theme_rhr +
theme(
axis.text.x = element_text(angle = 0, hjust = 0.5),
legend.position = "top"
)
b <-
plot_data %>%
dplyr::filter(
feature %in%
c("-log10(p)")
) %>%
ggplot(
aes(
x = rhr,
y = gerring
)
) +
geom_point(
size = 1.5,
alpha = 0.8,
shape = 21,
colour = "black",
aes(fill = significant)
) +
geom_hline(
yintercept = -log10(0.05/24054),
linetype = "dashed",
colour = "grey"
) +
geom_vline(
xintercept = -log10(bivar_thres),
linetype = "dashed",
colour = "grey"
) +
ggrepel::geom_label_repel(
aes(
label = label
),
size = 2,
alpha = 0.7,
box.padding = 0.3,
min.segment.length = unit(0, 'lines')
) +
facet_wrap(vars(feature)) +
expand_limits(x = 0, y = 0) +
coord_equal() +
labs(
x = "This study",
y = "Gerring et al. 2022"
) +
theme_rhr +
theme(
axis.text.x = element_text(angle = 0, hjust = 0.5),
legend.position = "top"
)
plot <-
ggpubr::ggarrange(
a, b,
labels = c("a", "b"),
common.legend = TRUE, align = "hv"
)
plot
Figure 4.10: Scatter plot of (a) rho and (b) -log10(p-value) for each pair of phenotypes that overlapped between our study and the study by Gerring et al. In (a) the dashed lines indicate significance thresholds in each study (Gerring et al. = 0.05/24,054; this study = 0.05/1603) and in (b) the dashed line represents the line x = y. Points in both plots are coloured by whether they passed significance thresholds in both studies or only the Gerring et al. study.
5 Session info
Show/hide
## ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────
## setting value
## version R version 4.0.5 (2021-03-31)
## os Ubuntu 16.04.7 LTS
## system x86_64, linux-gnu
## ui X11
## language (EN)
## collate en_GB.UTF-8
## ctype en_GB.UTF-8
## tz Europe/London
## date 2022-10-01
##
## ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
## ! package * version date lib source
## P abind 1.4-5 2016-07-21 [?] CRAN (R 4.0.5)
## P assertthat 0.2.1 2019-03-21 [?] CRAN (R 4.0.5)
## P backports 1.3.0 2021-10-27 [?] CRAN (R 4.0.5)
## P Biobase 2.50.0 2020-10-27 [?] Bioconductor
## P BiocGenerics * 0.36.1 2021-04-16 [?] Bioconductor
## P BiocParallel 1.24.1 2020-11-06 [?] Bioconductor
## P Biostrings 2.58.0 2020-10-27 [?] Bioconductor
## P bit 4.0.4 2020-08-04 [?] CRAN (R 4.0.5)
## P bit64 4.0.5 2020-08-30 [?] CRAN (R 4.0.5)
## P bitops 1.0-7 2021-04-24 [?] CRAN (R 4.0.5)
## P bookdown 0.22 2021-04-22 [?] CRAN (R 4.0.5)
## P broom 0.7.10 2021-10-31 [?] CRAN (R 4.0.5)
## P bslib 0.3.1 2021-10-06 [?] CRAN (R 4.0.5)
## P car 3.0-11 2021-06-27 [?] CRAN (R 4.0.5)
## P carData 3.0-4 2020-05-22 [?] CRAN (R 4.0.5)
## P cellranger 1.1.0 2016-07-27 [?] CRAN (R 4.0.5)
## P circlize * 0.4.13 2021-06-09 [?] CRAN (R 4.0.5)
## P cli 3.1.0 2021-10-27 [?] CRAN (R 4.0.5)
## P colorspace 2.0-2 2021-06-24 [?] CRAN (R 4.0.5)
## P cowplot 1.1.1 2020-12-30 [?] CRAN (R 4.0.5)
## P crayon 1.4.2 2021-10-29 [?] CRAN (R 4.0.5)
## P crosstalk 1.1.1 2021-01-12 [?] CRAN (R 4.0.5)
## P curl 4.3.2 2021-06-23 [?] CRAN (R 4.0.5)
## P data.table 1.14.2 2021-09-27 [?] CRAN (R 4.0.5)
## P DBI 1.1.1 2021-01-15 [?] CRAN (R 4.0.5)
## P dbplyr 2.1.1 2021-04-06 [?] CRAN (R 4.0.5)
## P DelayedArray 0.16.3 2021-03-24 [?] Bioconductor
## P digest 0.6.29 2021-12-01 [?] CRAN (R 4.0.5)
## P dplyr * 1.0.7 2021-06-18 [?] CRAN (R 4.0.5)
## P DT 0.19 2021-09-02 [?] CRAN (R 4.0.5)
## P ellipsis 0.3.2 2021-04-29 [?] CRAN (R 4.0.5)
## P evaluate 0.14 2019-05-28 [?] CRAN (R 4.0.5)
## P fansi 0.5.0 2021-05-25 [?] CRAN (R 4.0.5)
## P farver 2.1.0 2021-02-28 [?] CRAN (R 4.0.5)
## P fastmap 1.1.0 2021-01-25 [?] CRAN (R 4.0.5)
## P forcats * 0.5.1 2021-01-27 [?] CRAN (R 4.0.5)
## P foreign 0.8-81 2020-12-22 [?] CRAN (R 4.0.5)
## P fs 1.5.1 2021-11-30 [?] CRAN (R 4.0.5)
## P generics 0.1.1 2021-10-25 [?] CRAN (R 4.0.5)
## P GenomeInfoDb * 1.26.7 2021-04-08 [?] Bioconductor
## P GenomeInfoDbData 1.2.4 2021-07-03 [?] Bioconductor
## P GenomicAlignments 1.26.0 2020-10-27 [?] Bioconductor
## P GenomicRanges * 1.42.0 2020-10-27 [?] Bioconductor
## P ggforce 0.3.3 2021-03-05 [?] CRAN (R 4.0.5)
## P ggplot2 * 3.3.5 2021-06-25 [?] CRAN (R 4.0.5)
## P ggpubr 0.4.0 2020-06-27 [?] CRAN (R 4.0.5)
## P ggraph * 2.0.5 2021-02-23 [?] CRAN (R 4.0.5)
## P ggrepel 0.9.1 2021-01-15 [?] CRAN (R 4.0.5)
## P ggsignif 0.6.3 2021-09-09 [?] CRAN (R 4.0.5)
## P GlobalOptions 0.1.2 2020-06-10 [?] CRAN (R 4.0.5)
## P glue 1.5.1 2021-11-30 [?] CRAN (R 4.0.5)
## P graphlayouts 0.7.1 2020-10-26 [?] CRAN (R 4.0.5)
## P gridExtra 2.3 2017-09-09 [?] CRAN (R 4.0.5)
## P gtable 0.3.0 2019-03-25 [?] CRAN (R 4.0.5)
## P haven 2.4.3 2021-08-04 [?] CRAN (R 4.0.5)
## P here 1.0.1 2020-12-13 [?] CRAN (R 4.0.5)
## P highr 0.9 2021-04-16 [?] CRAN (R 4.0.5)
## P hms 1.1.1 2021-09-26 [?] CRAN (R 4.0.5)
## P htmltools 0.5.2 2021-08-25 [?] CRAN (R 4.0.5)
## P htmlwidgets 1.5.4 2021-09-08 [?] CRAN (R 4.0.5)
## P httr 1.4.2 2020-07-20 [?] CRAN (R 4.0.5)
## P igraph 1.2.7 2021-10-15 [?] CRAN (R 4.0.5)
## P IRanges * 2.24.1 2020-12-12 [?] Bioconductor
## P janitor * 2.1.0 2021-01-05 [?] CRAN (R 4.0.5)
## P jquerylib 0.1.4 2021-04-26 [?] CRAN (R 4.0.5)
## P jsonlite 1.7.2 2020-12-09 [?] CRAN (R 4.0.5)
## P knitr 1.36 2021-09-29 [?] CRAN (R 4.0.5)
## P labeling 0.4.2 2020-10-20 [?] CRAN (R 4.0.5)
## P lattice 0.20-44 2021-05-02 [?] CRAN (R 4.0.5)
## P LAVA 0.0.6 2021-09-01 [?] xgit (@7be3421)
## P lifecycle 1.0.1 2021-09-24 [?] CRAN (R 4.0.5)
## P lubridate 1.8.0 2021-10-07 [?] CRAN (R 4.0.5)
## P magrittr 2.0.1 2020-11-17 [?] CRAN (R 4.0.5)
## P MASS 7.3-54 2021-05-03 [?] CRAN (R 4.0.5)
## P Matrix 1.3-4 2021-06-01 [?] CRAN (R 4.0.5)
## P MatrixGenerics 1.2.1 2021-01-30 [?] Bioconductor
## P matrixStats 0.61.0 2021-09-17 [?] CRAN (R 4.0.5)
## P modelr 0.1.8 2020-05-19 [?] CRAN (R 4.0.5)
## P munsell 0.5.0 2018-06-12 [?] CRAN (R 4.0.5)
## P openxlsx 4.2.4 2021-06-16 [?] CRAN (R 4.0.5)
## P pillar 1.6.4 2021-10-18 [?] CRAN (R 4.0.5)
## P pkgconfig 2.0.3 2019-09-22 [?] CRAN (R 4.0.5)
## P polyclip 1.10-0 2019-03-14 [?] CRAN (R 4.0.5)
## P purrr * 0.3.4 2020-04-17 [?] CRAN (R 4.0.5)
## P R6 2.5.1 2021-08-19 [?] CRAN (R 4.0.5)
## P RColorBrewer 1.1-2 2014-12-07 [?] CRAN (R 4.0.5)
## P Rcpp 1.0.7 2021-07-07 [?] CRAN (R 4.0.5)
## P RCurl 1.98-1.5 2021-09-17 [?] CRAN (R 4.0.5)
## P readr * 2.0.2 2021-09-27 [?] CRAN (R 4.0.5)
## P readxl 1.3.1 2019-03-13 [?] CRAN (R 4.0.5)
## P reprex 2.0.0 2021-04-02 [?] CRAN (R 4.0.5)
## P rio 0.5.27 2021-06-21 [?] CRAN (R 4.0.5)
## P rlang 0.4.12 2021-10-18 [?] CRAN (R 4.0.5)
## P rmarkdown 2.11 2021-09-14 [?] CRAN (R 4.0.5)
## P rprojroot 2.0.2 2020-11-15 [?] CRAN (R 4.0.5)
## P Rsamtools 2.6.0 2020-10-27 [?] Bioconductor
## P rstatix 0.7.0 2021-02-13 [?] CRAN (R 4.0.5)
## P rstudioapi 0.13 2020-11-12 [?] CRAN (R 4.0.5)
## P rtracklayer * 1.50.0 2020-10-27 [?] Bioconductor
## P rvest 1.0.1 2021-07-26 [?] CRAN (R 4.0.5)
## P S4Vectors * 0.28.1 2020-12-09 [?] Bioconductor
## P sass 0.4.0 2021-05-12 [?] CRAN (R 4.0.5)
## P scales 1.1.1 2020-05-11 [?] CRAN (R 4.0.5)
## P sessioninfo * 1.1.1 2018-11-05 [?] CRAN (R 4.0.5)
## P shape 1.4.6 2021-05-19 [?] CRAN (R 4.0.5)
## P snakecase 0.11.0 2019-05-25 [?] CRAN (R 4.0.5)
## P stringi 1.7.6 2021-11-29 [?] CRAN (R 4.0.5)
## P stringr * 1.4.0 2019-02-10 [?] CRAN (R 4.0.5)
## P SummarizedExperiment 1.20.0 2020-10-27 [?] Bioconductor
## P tibble * 3.1.6 2021-11-07 [?] CRAN (R 4.0.5)
## P tidygraph 1.2.0 2020-05-12 [?] CRAN (R 4.0.5)
## P tidyr * 1.1.4 2021-09-27 [?] CRAN (R 4.0.5)
## P tidyselect 1.1.1 2021-04-30 [?] CRAN (R 4.0.5)
## P tidyverse * 1.3.1 2021-04-15 [?] CRAN (R 4.0.5)
## P tweenr 1.0.2 2021-03-23 [?] CRAN (R 4.0.5)
## P tzdb 0.2.0 2021-10-27 [?] CRAN (R 4.0.5)
## P utf8 1.2.2 2021-07-24 [?] CRAN (R 4.0.5)
## P vctrs 0.3.8 2021-04-29 [?] CRAN (R 4.0.5)
## P viridis 0.6.2 2021-10-13 [?] CRAN (R 4.0.5)
## P viridisLite 0.4.0 2021-04-13 [?] CRAN (R 4.0.5)
## P vroom 1.5.5 2021-09-14 [?] CRAN (R 4.0.5)
## P withr 2.4.3 2021-11-30 [?] CRAN (R 4.0.5)
## P xfun 0.27 2021-10-18 [?] CRAN (R 4.0.5)
## P XML 3.99-0.8 2021-09-17 [?] CRAN (R 4.0.5)
## P xml2 1.3.2 2020-04-23 [?] CRAN (R 4.0.5)
## P XVector 0.30.0 2020-10-27 [?] Bioconductor
## P yaml 2.2.1 2020-02-01 [?] CRAN (R 4.0.5)
## P zip 2.2.0 2021-05-31 [?] CRAN (R 4.0.5)
## P zlibbioc 1.36.0 2020-10-27 [?] Bioconductor
##
## [1] /home/rreynolds/misc_projects/neurodegen-psych-local-corr/renv/library/R-4.0/x86_64-pc-linux-gnu
## [2] /opt/R/4.0.5/lib/R/library
##
## P ── Loaded and on-disk path mismatch.