Case study: C1S xQTL and AD GWAS

This notebook documents the analysis of xQTL case study on a targeted gene, C1S.

• Section 0: Sanity check
• Section 1: Fine-mapping for xQTL and GWAS
• Section 2: Multi-context colocalization with Bellenguez 2022
• Section 3: Refinement of colocalized loci with other AD GWAS
• Section 4: Assessment of multi-context xQTL effect sizes
• Section 5: Multi-context causal TWAS (including conventional TWAS and MR)
• Section 6: Context specific multi-gene fine-mapping
• Section 7: Epigenomic QTL and their target regions
• Section 8: Context focused validation in other xQTL data
• Section 9: Non-linear effects of xQTL
• Section 10: in silico functional studies in iPSC model
• Section 11: Functional annotations of selected loci
• Section 12: Candidate loci as trans-xQTL

Overview

FunGen-xQTL resource contains 67 xQTL contexts as well as 9 AD GWAS fine-mapped genome-wide. The overarching goal for case studies is to use these resource to raise questions and learn more about gene targets of interest.

Overall, a case study consists of the following aspects:

• Check the basic information of the gene
• Check the existing xQTL and integrative analysis results, roughly including
  • Summary table for univariate fine-mapping
  • Marginal association results
  • Multi-gene and multi-context fine-mapping
  • Multi-context colocalization with AD GWAS
  • TWAS, MR and causal TWAS
  • Integration with epigenetic QTL
  • Quantile QTL
  • Interaction QTL
  • Validation:
    • Additional xQTL data in FunGen-xQTL
    • Additional AD GWAS data-set already generated by FunGen-xQTL
  • In silico functional studies
    • Additional iPSC data-sets
  • Functional annotations of variants, particularly in relevant cellular contexts
• Creative thinking: generate hypothesis, search in literature, raise questions to discuss

Computing environment setup

Interactive analysis will be done on AWS cloud powered by MemVerge. Please contact Gao Wang to setup accounts for you to start the analysis.

Please follow instructions on https://wanggroup.org/productivity_tips/mmcloud-interactive to configure your computing environment. Here are some additional packages you need to install after the initial configuration, in order to perform these analysis: in terminal with bash:

micromamba install -n r_libs r-pecotmr

in R:

install.packages('BEDMatrix')

How to Use This Notebook

0. Before you start: Load functions from cb_plot.R and utilis.R, located at /data/interactive_analysis/rf2872/codes/. These functions and resources are packaged to streamline the analysis and ensure everything is as clean as possible. And also the codes for ColocBoost under path /data/colocalization/colocboost/R.
1. Inside of this notebook, use sed -i or Ctrl+F to replace the gene C1S with the gene you want to analyze.
2. For detailed analysis in some sections, please refer to the corresponding analysis notebooks as indicated. These companion notebooks are available under this same folder. The rest of the tasks can be completed with a few lines of code, as demonstrated in this notebook.
3. Similarly for the companion notebooks you should also use the sed -i or Ctrl+F replacing gene_name (C1S in this case) with the gene you want to investigate.

While using this notebook, you may need to generate three intermediate files from Sections 1 and 2, which will be useful for downstream analysis:

• a. Section 1:
  • Fine-mapping contexts that indicate shared signals with AD, C1S_finemapping_contexts.rds. This can be used as input for Section 8 the multi-cohort validation step
  • A subset of the xQTL-AD table, Fungen_xQTL_allQTL.overlapped.gwas.export.C1S.rds. This can be used as input for Section 12.
• b. Section 2: A variant list showing colocalization in cohorts we analyzed with ColocBoost, C1S_colocboost_res.rds. this can be used as input for Sections 7, 9, 10, and 12.

Section 0: Sanity check

Check the basic information of the gene

• To gain a preliminary understanding of this gene’s expression—specifically, whether it is cell-specific—can help us quickly determine if our results are consistent with expectations.

Useful websites:

1. check gene function, (immune) cell type specificity, tissue specifity, protein location: https://www.proteinatlas.org
2. check gene position and structure: https://www.ncbi.nlm.nih.gov/gene/
3. other collective information: https://www.genecards.org

Check the existing results which are inputs to this analysis

source('../../codes/cb_plot.R')
source('../../codes/utilis.R')
for(file in list.files("/data/colocalization/colocboost/R", pattern = ".R", full.names = T)){
          source(file)
        }
gene_name = 'C1S'

dir.create(paste0('plots/', gene_name), recursive = T)

get basic target gene information

target_gene_info <- get_gene_info(gene_name = gene_name)
target_gene_info

$gene_info

A data.table: 1 x 14

region_id	#chr	start	end	TSS	LD_matrix_id	LD_sumstats_id	LD_sumstats_id_old	TADB_index	TADB_id	gene_start	gene_end	sliding_windows	gene_name
<chr>	<chr>	<dbl>	<dbl>	<int>	<chr>	<chr>	<chr>	<chr>	<chr>	<int>	<int>	<chr>	<chr>
ENSG00000182326	chr12	4280000	8440000	6988258	chr12:3219245-4882038,chr12:4882038-6411328,chr12:6411328-8138114,chr12:8138114-10024521	12_3219245-4882038,12_4882038-6411328,12_6411328-8138114,12_8138114-10024521	12_3219245-4882038,12_4882038_6411328,12_6411328_8138114,12_8138114_10024521	TADB_947,TADB_948,TADB_949	chr12_3269770_8218574,chr12_5556863_9438425,chr12_6812512_12521686	6988259	7071032	chr12:0-5851927,chr12:529890-8218574,chr12:970785-9438425,chr12:3269770-12521686,chr12:5556863-13471233,chr12:6812512-13829975	C1S

$target_LD_ids

A matrix: 1 x 4 of type chr

chr12:3219245-4882038

chr12:4882038-6411328

chr12:6411328-8138114

chr12:8138114-10024521

$target_sums_ids

A matrix: 1 x 4 of type chr

12_3219245-4882038

12_4882038-6411328

12_6411328-8138114

12_8138114-10024521

$gene_region

'chr12:4280000-8440000'

$target_TAD_ids

A matrix: 1 x 3 of type chr

chr12_3269770_8218574

chr12_5556863_9438425

chr12_6812512_12521686

gene_id = target_gene_info$gene_info$region_id
chrom = target_gene_info$gene_info$`#chr`

Take a quick look for the expression of target gene in ROSMAP bulk data, we don't want them to be too low

source('../../codes/utilis.R')
expression_in_rosmap_bulk(target_gene_info)

Section 1: Fine-mapping for xQTL and GWAS

see notebook

region_p

Bellenguez et al GWAS signals has many overlap with CS from other xQTL sources. This motivates us to look further. The figure above shows the ranges of CS to give us a loci level idea. Below, we show the variants in those CS, color-coding the variants that are shared between them in the same color. In particular, AD GWAS signals are also captured by a few xQTL data, although at this point we don't have formal statistical (colocalization) evidences for these observations yet.

pip_p

summary ADxQTL table

This summary table have been produced by conserving the loci overlap between all-source AD GWAS credible sets (CS95%mincorr,CS95% and CS70%mincorr ), and the xQTL credible sets (CS95%mincorr, CS70%mincorr, and SNPs with pip>0.2)

options(repr.matrix.max.cols=150, repr.matrix.max.rows=20)
gname=gene_name
gwas_anno<-fread('../../Landscape_analysis/res_gwas_qtl_anno_ADGWAS_xQTL_noMSBBsQTL_filtered.csv.gz')
gwas_anno[gene_name==gname]

A data.table: 1 x 57

gwas_locus_id	gwas_source	gwas_variant	locuscontext_id	gene_id	gene_name	context	n.variant.hit	n.variant.overlapping	pct.variant.overlap	top.snp.dist	gwas_region.hit	gwas_region.hit.size	gwas_n.variant.hit	gwas_pos	gwas_maf	gwas_pip	overlapping_gwas	overlapping_gwas_loci	study	region	variant_id	pip	method	betahat	sebetahat	chr	pos	maf	annotated_susie_cs	qtl_type	SingleNuc	BulkCellType	Metab	assay_type	cohort	cohort_type	cell_type	qtl_category	context_hi	junction_region	tested_molecule_event	credibleset	qtl_hit	start.hit	end.hit	region.hit	region.hit.size	overlap.with.xqtl	gwas_z	n.gwas	n.context	qtl_z	qtl_pvalue	n.hit	gwas_region.hit.union	AD_context_score
<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<int>	<int>	<dbl>	<int>	<chr>	<int>	<int>	<int>	<lgl>	<dbl>	<chr>	<chr>	<chr>	<chr>	<chr>	<dbl>	<chr>	<dbl>	<dbl>	<chr>	<int>	<dbl>	<int>	<chr>	<lgl>	<lgl>	<lgl>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<lgl>	<int>	<int>	<chr>	<int>	<lgl>	<dbl>	<int>	<int>	<dbl>	<dbl>	<int>	<chr>	<dbl>
chr12_6411328_8138114_AD_Bellenguez_2022_cs95_1	AD_Bellenguez_2022	chr12:7055497:C:T	ENSG00000182326_DLPFC_Bennett_pQTL_cs95corr_1	ENSG00000182326	C1S	DLPFC_Bennett_pQTL	18	14	0.7777778	8537	chr12-7048747-7073858	25112	21	7055497	NA	0.1856597			DLPFC_Bennett_pQTL	ENSG00000182326	chr12:7064034:T:A	0.2851886	top_loci	-0.5049025	0.0524393	chr12	7064034	0.1358173	NA	pQTL	FALSE	FALSE	FALSE	BulkBrain	Bennett	Bennett_BulkBrain	bulk_brain	others QTLs	bulk_brain_pQTL		ENSG00000182326	cs95corr_1	TRUE	7063032	7074075	chr12-7063032-7074075	11044	TRUE	4.909091	1	6	-9.628323	0	1	chr12-7048747-7073859	4.541893

here we can see than C1S is associated only with 1 AD locus : in AD bellenguez study, and 1 locuscontext_id, DLPFC_Bennett_pQTL.

We can check now how this AD locus is 'pure', i.e. non associated with others context than C1S

ad_loci<-gwas_anno[gene_name==gname]$gwas_locus_id
gwas_anno[gwas_locus_id%in%ad_loci]

A data.table: 6 x 57

gwas_locus_id	gwas_source	gwas_variant	locuscontext_id	gene_id	gene_name	context	n.variant.hit	n.variant.overlapping	pct.variant.overlap	top.snp.dist	gwas_region.hit	gwas_region.hit.size	gwas_n.variant.hit	gwas_pos	gwas_maf	gwas_pip	overlapping_gwas	overlapping_gwas_loci	study	region	variant_id	pip	method	betahat	sebetahat	chr	pos	maf	annotated_susie_cs	qtl_type	SingleNuc	BulkCellType	Metab	assay_type	cohort	cohort_type	cell_type	qtl_category	context_hi	junction_region	tested_molecule_event	credibleset	qtl_hit	start.hit	end.hit	region.hit	region.hit.size	overlap.with.xqtl	gwas_z	n.gwas	n.context	qtl_z	qtl_pvalue	n.hit	gwas_region.hit.union	AD_context_score
<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<int>	<int>	<dbl>	<int>	<chr>	<int>	<int>	<int>	<lgl>	<dbl>	<chr>	<chr>	<chr>	<chr>	<chr>	<dbl>	<chr>	<dbl>	<dbl>	<chr>	<int>	<dbl>	<int>	<chr>	<lgl>	<lgl>	<lgl>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<lgl>	<int>	<int>	<chr>	<int>	<lgl>	<dbl>	<int>	<int>	<dbl>	<dbl>	<int>	<chr>	<dbl>
chr12_6411328_8138114_AD_Bellenguez_2022_cs95_1	AD_Bellenguez_2022	chr12:7055497:C:T	ENSG00000111254_12:4638196:4648745:clu_14532_ROSMAP_AC_sQTL_cs70corr_3	ENSG00000111254	AKAP3	ROSMAP_AC_sQTL	22	16	0.7272727	7535	chr12-7048747-7073858	25112	21	7055497	NA	0.1856597			ROSMAP_AC_sQTL_intron_chr12:4638196:4648745:clu_14532_-:ENSG00000111254	ENSG00000111254	chr12:7063032:G:A	0.0587871	top_loci	-0.3212800	0.07467939	chr12	7063032	0.1188870	NA	sQTL(leafcutter)	FALSE	FALSE	FALSE	BulkBrain	ROSMAP	ROSMAP_BulkBrain	bulk_brain	sQTL	bulk_brain_sQTL(leafcutter)	12:4638196:4648745:clu_14532	ENSG00000111254_12:4638196:4648745:clu_14532	cs70corr_3	TRUE	7059466	7074644	chr12-7059466-7074644	15179	TRUE	4.909091	1	6	-4.302124	1.691683e-05	2	chr12-7048747-7073859	4.541893
chr12_6411328_8138114_AD_Bellenguez_2022_cs95_1	AD_Bellenguez_2022	chr12:7055497:C:T	ENSG00000111254_12:4638196:4648816:clu_14532_ROSMAP_AC_sQTL_cs70corr_2	ENSG00000111254	AKAP3	ROSMAP_AC_sQTL	19	14	0.7368421	6860	chr12-7048747-7073858	25112	21	7055497	NA	0.1856597			ROSMAP_AC_sQTL_intron_chr12:4638196:4648816:clu_14532_-:ENSG00000111254	ENSG00000111254	chr12:7062357:C:T	0.1703815	top_loci	-0.3051157	0.06109810	chr12	7062357	0.1340641	NA	sQTL(leafcutter)	FALSE	FALSE	FALSE	BulkBrain	ROSMAP	ROSMAP_BulkBrain	bulk_brain	sQTL	bulk_brain_sQTL(leafcutter)	12:4638196:4648816:clu_14532	ENSG00000111254_12:4638196:4648816:clu_14532	cs70corr_2	TRUE	7059466	7074644	chr12-7059466-7074644	15179	TRUE	4.909091	1	6	-4.993866	5.918236e-07	2	chr12-7048747-7073859	4.541893
chr12_6411328_8138114_AD_Bellenguez_2022_cs95_1	AD_Bellenguez_2022	chr12:7055497:C:T	ENSG00000126746_MiGA_SVZ_eQTL_cs95corr_4	ENSG00000126746	ZNF384	MiGA_SVZ_eQTL	18	15	0.8333333	3969	chr12-7048747-7073858	25112	21	7055497	NA	0.1856597			MiGA_SVZ_eQTL	ENSG00000126746	chr12:7059466:A:G	0.0555556	top_loci	-0.1327249	0.10974247	chr12	7059466	0.1808511	NA	eQTL	FALSE	TRUE	FALSE	BulkCellType	MiGA	MiGA_BulkCellType	bulk_microglia	bulk-eQTL	bulk_microglia_eQTL		ENSG00000126746	cs95corr_4	TRUE	7059466	7074644	chr12-7059466-7074644	15179	TRUE	4.909091	1	6	-1.209422	2.265009e-01	2	chr12-7048747-7073859	4.541893
chr12_6411328_8138114_AD_Bellenguez_2022_cs95_1	AD_Bellenguez_2022	chr12:7055497:C:T	ENSG00000126749_12:6970966:6974339:clu_15629_ROSMAP_PCC_sQTL_cs95corr_1	ENSG00000126749	EMG1	ROSMAP_PCC_sQTL	20	16	0.8000000	3969	chr12-7048747-7073858	25112	21	7055497	NA	0.1856597			ROSMAP_PCC_sQTL_intron_chr12:6970966:6974339:clu_15629_+:ENSG00000126749	ENSG00000126749	chr12:7059466:A:G	0.3928376	top_loci	0.4883281	0.06985865	chr12	7059466	0.1383220	NA	sQTL(leafcutter)	FALSE	FALSE	FALSE	BulkBrain	ROSMAP	ROSMAP_BulkBrain	bulk_brain	sQTL	bulk_brain_sQTL(leafcutter)	12:6970966:6974339:clu_15629	ENSG00000126749_12:6970966:6974339:clu_15629	cs95corr_1	TRUE	7059466	7074644	chr12-7059466-7074644	15179	TRUE	4.909091	1	6	6.990231	2.744249e-12	2	chr12-7048747-7073859	4.541893
chr12_6411328_8138114_AD_Bellenguez_2022_cs95_1	AD_Bellenguez_2022	chr12:7055497:C:T	ENSG00000139197_MiGA_THA_eQTL_cs95corr_1	ENSG00000139197	PEX5	MiGA_THA_eQTL	6	4	0.6666667	0	chr12-7048747-7073858	25112	21	7055497	NA	0.1856597			MiGA_THA_eQTL	ENSG00000139197	chr12:7055497:C:T	0.2613628	top_loci	-0.5520035	0.15287697	chr12	7055497	0.2410714	NA	eQTL	FALSE	TRUE	FALSE	BulkCellType	MiGA	MiGA_BulkCellType	bulk_microglia	bulk-eQTL	bulk_microglia_eQTL		ENSG00000139197	cs95corr_1	TRUE	7048232	7057810	chr12-7048232-7057810	9579	TRUE	4.909091	1	6	-3.610770	3.052895e-04	2	chr12-7048747-7073859	4.541893
chr12_6411328_8138114_AD_Bellenguez_2022_cs95_1	AD_Bellenguez_2022	chr12:7055497:C:T	ENSG00000182326_DLPFC_Bennett_pQTL_cs95corr_1	ENSG00000182326	C1S	DLPFC_Bennett_pQTL	18	14	0.7777778	8537	chr12-7048747-7073858	25112	21	7055497	NA	0.1856597			DLPFC_Bennett_pQTL	ENSG00000182326	chr12:7064034:T:A	0.2851886	top_loci	-0.5049025	0.05243930	chr12	7064034	0.1358173	NA	pQTL	FALSE	FALSE	FALSE	BulkBrain	Bennett	Bennett_BulkBrain	bulk_brain	others QTLs	bulk_brain_pQTL		ENSG00000182326	cs95corr_1	TRUE	7063032	7074075	chr12-7063032-7074075	11044	TRUE	4.909091	1	6	-9.628323	0.000000e+00	1	chr12-7048747-7073859	4.541893

The locus is not pure, it is associated with other genes : ZNF384 and PEX5 in MIGA eQTl, and 3 others genes in sQTL data. Interestingly only PEX5 MiGA_THA_eQTL have the same top SNP (based on pip ) than AD. this data is however to take with caution (MIGA is a little dataset and known to be probably inflated for false positive), so look at the coloc results will help.

Section 2: Multi-context colocalization with Bellenguez 2022

This is done using ColocBoost. The most updated version of ColocBoost results are under path s3://statfungen/ftp_fgc_xqtl/analysis_result/ColocBoost/2024_9/

 cb_res <- readRDS(paste0("/data/analysis_result/ColocBoost/2024_9/",gene_id,"_res.rds") )

#str(cb_res)

#save colocboost results
cb_res_table <- get_cb_summary(cb_res) 

saveRDS(cb_res_table, paste0(gene_name, "_colocboost_res.rds"))

cb <- plot_cb(cb_res = cb_res, cex.pheno = 1.5, x.phen = -0.2)

pdf('plots/C1S/sec2.colocboost_res.pdf', width = 10, height = 5)
replayPlot(cb$p)
dev.off()

# colocalized variants
cb_res_table

$max_abs_z_variant =

# effect sign for each coloc sets
get_effect_sign_csets(cb_res)

# LD between coloc sets
get_between_purity_simple(cb_res, gene.name = gene_id, path = '/data/colocalization/QTL_data/eQTL/')

Here, different colors refer to different 95% Colocalization Sets (CoS, a metric developed in ColocBoost indicating that there is 95% probabilty that this CoS captures a colocalization event). We only included ROSMAP data for this particular ColocBoost analysis. In this case, we observe cell specific eQTL, bulk sQTL colocalization on ROSMAP data with AD as two separate CoS, suggesting two putative causal signals. We did not detect colocalization with pQTL of statistical significance although from Section 1 there are some overlap with pQTL signals in fine-mapping CS, the overlapped variants in CS have small PIP.

Section 3: Refinement of colocalized loci with other AD GWAS

Here we refine the colocalization with other AD GWAS to iron out any heterogeniety between studies (heterogeniety can come from many sources), to get additional candidate loci from these more heterogenous sources as candidates to study.

AD_cohorts <- c('AD_Jansen_2021', 'AD_Bellenguez_EADB_2022', 'AD_Bellenguez_EADI_2022',
             'AD_Kunkle_Stage1_2019', 'AD_Wightman_Excluding23andMe_2021',
             'AD_Wightman_ExcludingUKBand23andME_2021', 'AD_Wightman_Full_2021')
cb_ad <- plot_cb(cb_res = cb_res, cex.pheno = 1.5, x.phen = -0.2, add_gwas = TRUE, gene_id = gene_id, cohorts = AD_cohorts)

No pvalue cutoff. Extract all variants names.

pdf('plots/C1S/sec3.colocboost_res_allad.pdf', width = 10, height = 5)
replayPlot(cb_ad$p)
dev.off()

png: 2

Section 4: Assessment of multi-context xQTL effect sizes

Option 1: ColocBoost + MASH

Use colocboost variants and check for mash posterior contrast to see if the effect size are shared or specific or even opposite. Advantage is that colocboost result is AD GWAS informed; issue is that marginal posterior effects is not always the joint

mash_p <- mash_plot(gene_name = 'C1S')
for (plot in mash_p) {
    print(plot)
}

Option 2: mvSuSiE

Use mvSuSiE multicontext fine-mapping results --- the bubble plot to check posterior effects. Issue is that we don't have this results yet, and this is limited to one cohort at a time, without information from AD.

We should go for option 1 by default and if we want to make claim about opposite effect size we double-check with mvSuSiE multicontext analysis.

Section 5: Multi-context causal TWAS (including conventional TWAS and MR)

The most updated version of cTWAS analysis are under path s3://statfungen/ftp_fgc_xqtl/analysis_result/cTWAS/

TWAS results

We report TWAS from all contexts and methods from the pipeline. Here we will filter it down to the best performing methods and only keep contexts that are significant.

plot_TWAS_res(gene_id = gene_id, gene_name = gene_name)

MR results

This is only available for genes that are deemed significant in TWAS and have summary statistics available for effect size and standard errors in GWAS, in addition to z-scores --- current version does not support z-scores although we will soon also support z-scores in MR using MAF from reference panel.

cTWAS results

To be updated soon.

Section 6: Context specific multi-gene fine-mapping

A quick analysis: using the xQTL-AD summary table (flatten table)

We extract from xQTL-AD summary table the variants to get other genes also have CS with the variants shared by target gene and AD.

multigene_flat <- get_multigene_multicontext_flatten('Fungen_xQTL_allQTL.overlapped.gwas.export.C1S.rds', sQTL = 'no_MSBB')
multigene_flat

A data.frame: 73 x 6

gene_id	#chr	start	end	gene_name	contexts
<chr>	<chr>	<int>	<int>	<chr>	<chr>
ENSG00000008323	chr12	6310435	6310436	PLEKHG6	MiGA_THA_eQTL
ENSG00000010219	chr12	4562203	4562204	DYRK4	ROSMAP_PCC_sQTL
ENSG00000010278	chr12	6199714	6199715	CD9	ROSMAP_DLPFC_sQTL
ENSG00000010292	chr12	6493355	6493356	NCAPD2	ROSMAP_AC_sQTL,ROSMAP_PCC_sQTL,ROSMAP_DLPFC_sQTL
ENSG00000010610	chr12	6786857	6786858	CD4	ROSMAP_AC_sQTL,STARNET_eQTL
ENSG00000010626	chr12	6873568	6873569	LRRC23	MiGA_SVZ_eQTL
ENSG00000011105	chr12	3077354	3077355	TSPAN9	ROSMAP_AC_sQTL,ROSMAP_DLPFC_sQTL
ENSG00000047617	chr12	5946231	5946232	ANO2	MiGA_GFM_eQTL,Exc_Kellis_eQTL,Exc_mega_eQTL,ROSMAP_PCC_sQTL,ROSMAP_DLPFC_sQTL
ENSG00000065970	chr12	8032715	8032716	FOXJ2	STARNET_eQTL
ENSG00000067182	chr12	6342113	6342114	TNFRSF1A	ROSMAP_AC_sQTL
ENSG00000069493	chr12	9664968	9664969	CLEC2D	MiGA_SVZ_eQTL
ENSG00000078237	chr12	4307762	4307763	TIGAR	MiGA_GTS_eQTL,MiGA_SVZ_eQTL
ENSG00000089692	chr12	6772518	6772519	LAG3	MiGA_THA_eQTL
ENSG00000089693	chr12	6767474	6767475	MLF2	ROSMAP_DLPFC_sQTL
ENSG00000089818	chr12	8076938	8076939	NECAP1	STARNET_eQTL
ENSG00000110799	chr12	6124769	6124770	VWF	DLPFC_DeJager_eQTL,Exc_Kellis_eQTL,ROSMAP_PCC_sQTL,ROSMAP_DLPFC_sQTL
ENSG00000110811	chr12	6828406	6828407	P3H3	DLPFC_DeJager_eQTL,AC_DeJager_eQTL,Exc_mega_eQTL,ROSMAP_DLPFC_sQTL
ENSG00000110852	chr12	9869353	9869354	CLEC2B	BM_22_MSBB_eQTL
ENSG00000111247	chr12	4538797	4538798	RAD51AP1	MiGA_SVZ_eQTL
ENSG00000111254	chr12	4649050	4649051	AKAP3	ROSMAP_AC_sQTL
ENSG00000111262	chr12	4909904	4909905	KCNA1	AC_DeJager_eQTL
ENSG00000111319	chr12	6377729	6377730	SCNN1A	monocyte_ROSMAP_eQTL
ENSG00000111321	chr12	6375044	6375045	LTBR	ROSMAP_AC_sQTL,ROSMAP_PCC_sQTL,ROSMAP_DLPFC_sQTL
ENSG00000111639	chr12	6493840	6493841	MRPL51	ROSMAP_DLPFC_sQTL
ENSG00000111640	chr12	6534516	6534517	GAPDH	MSBB_BM36_pQTL,DLPFC_DeJager_eQTL,Exc_Kellis_eQTL,Ast_mega_eQTL,OPC_mega_eQTL,ROSMAP_AC_sQTL,ROSMAP_DLPFC_sQTL
ENSG00000111641	chr12	6568690	6568691	NOP2	MiGA_SVZ_eQTL,ROSMAP_DLPFC_sQTL
ENSG00000111642	chr12	6614523	6614524	CHD4	MiGA_THA_eQTL,ROSMAP_DLPFC_sQTL,ROSMAP_PCC_sQTL
ENSG00000111653	chr12	6663141	6663142	ING4	MSBB_BM36_pQTL
ENSG00000111664	chr12	6840210	6840211	GNB3	MiGA_THA_eQTL,BM_36_MSBB_eQTL,BM_44_MSBB_eQTL,Exc_Kellis_eQTL,Inh_mega_eQTL,ROSMAP_AC_sQTL,ROSMAP_DLPFC_sQTL,ROSMAP_PCC_sQTL
ENSG00000111665	chr12	6851291	6851292	CDCA3	MiGA_GFM_eQTL,MiGA_GTS_eQTL,MiGA_THA_eQTL
...	...	...	...	...	...
ENSG00000130035	chr12	4720399	4720400	GALNT8	Exc_DeJager_eQTL
ENSG00000130037	chr12	5043878	5043879	KCNA5	DLPFC_DeJager_eQTL
ENSG00000130038	chr12	3764818	3764819	CRACR2A	ROSMAP_AC_sQTL
ENSG00000139178	chr12	7109237	7109238	C1RL	ROSMAP_AC_sQTL,ROSMAP_DLPFC_sQTL
ENSG00000139182	chr12	7129697	7129698	CLSTN3	MSBB_BM36_pQTL,Exc_mega_eQTL,ROSMAP_AC_sQTL,ROSMAP_PCC_sQTL,ROSMAP_DLPFC_sQTL
ENSG00000139190	chr12	6470676	6470677	VAMP1	monocyte_ROSMAP_eQTL
ENSG00000139192	chr12	6451689	6451690	TAPBPL	Knight_eQTL,MiGA_THA_eQTL,ROSMAP_PCC_sQTL,ROSMAP_DLPFC_sQTL
ENSG00000139197	chr12	7188684	7188685	PEX5	Knight_eQTL,BM_22_MSBB_eQTL,OPC_DeJager_eQTL,DLPFC_DeJager_eQTL,PCC_DeJager_eQTL,AC_DeJager_eQTL,Exc_Kellis_eQTL,OPC_mega_eQTL,ROSMAP_AC_sQTL,ROSMAP_DLPFC_sQTL,ROSMAP_PCC_sQTL
ENSG00000139200	chr12	6700814	6700815	PIANP	MiGA_THA_eQTL,ROSMAP_AC_sQTL,ROSMAP_PCC_sQTL,ROSMAP_DLPFC_sQTL
ENSG00000151067	chr12	1970785	1970786	CACNA1C	monocyte_ROSMAP_eQTL,ROSMAP_AC_sQTL
ENSG00000151079	chr12	4809175	4809176	KCNA6	MiGA_GTS_eQTL,MSBB_BM36_pQTL
ENSG00000159335	chr12	6765515	6765516	PTMS	ROSMAP_AC_sQTL,ROSMAP_PCC_sQTL,ROSMAP_DLPFC_sQTL
ENSG00000166523	chr12	8540904	8540905	CLEC4E	ROSMAP_DLPFC_sQTL
ENSG00000166535	chr12	8822620	8822621	A2ML1	ROSMAP_DLPFC_sQTL
ENSG00000171847	chr12	8227617	8227618	FAM90A1	ROSMAP_PCC_sQTL
ENSG00000171860	chr12	8066358	8066359	C3AR1	MiGA_SVZ_eQTL,STARNET_eQTL
ENSG00000173262	chr12	7891147	7891148	SLC2A14	Knight_eQTL,MiGA_GFM_eQTL,MiGA_THA_eQTL,STARNET_eQTL
ENSG00000173391	chr12	10172137	10172138	OLR1	Mic_13_Kellis_eQTL
ENSG00000175899	chr12	9116228	9116229	A2M	ROSMAP_DLPFC_sQTL
ENSG00000177575	chr12	7503892	7503893	CD163	MiGA_SVZ_eQTL,ROSMAP_PCC_sQTL,ROSMAP_DLPFC_sQTL
ENSG00000177675	chr12	7479896	7479897	CD163L1	Mic_Kellis_eQTL,ROSMAP_PCC_sQTL
ENSG00000180574	chr12	10505601	10505602	EIF2S3B	MiGA_THA_eQTL
ENSG00000184344	chr12	7695774	7695775	GDF3	MiGA_GFM_eQTL
ENSG00000212124	chr12	11022619	11022620	TAS2R19	BM_22_MSBB_eQTL
ENSG00000215021	chr12	6970779	6970780	PHB2	ROSMAP_AC_sQTL,ROSMAP_PCC_sQTL,ROSMAP_DLPFC_sQTL
ENSG00000250510	chr12	6821703	6821704	GPR162	MiGA_SVZ_eQTL,ROSMAP_DLPFC_sQTL
ENSG00000255639	chr12	4649100	4649101	AC005833.1	MiGA_THA_eQTL,ROSMAP_DLPFC_sQTL
ENSG00000256537	chr12	11171193	11171194	SMIM10L1	ROSMAP_PCC_sQTL
ENSG00000285238	chr12	6607366	6607367	AC006064.6	ROSMAP_AC_sQTL,ROSMAP_DLPFC_sQTL,ROSMAP_PCC_sQTL
ENSG00000285901	chr12	4273761	4273762	AC008012.1	ROSMAP_PCC_sQTL

Other genes implicated are PROC and HS6ST1 in MiGA cohort which may share causal eQTL with C1S. Further look into the data-set --- using these genes as targets and repeating what we did above for C1S --- might be needed to establish a more certain conclusion.

Alternatively, we may be able to apply a multi-gene statistical fine-mapping test on C1S region to find these genes, as you will see in the section below.

A statistically solid approach: mvSuSiE multi-gene analysis

This multi-gene fine-mapping analysis was done for each xQTL context separately. We will need to check the results for all contexts where this gene has an xQTL, in order to identify if there are other genes also sharing the same xQTL with this target gene. We included other genes in the same TAD window along with this gene, and extended it into a sliding window approach to include multiple TADs just in case. You need to check the sliding windows belongs to that gene (TSS) on analysis repo.

sliding_windows <- target_gene_info$gene_info$sliding_windows %>% strsplit(., ',') %>% unlist %>% as.character
sliding_windows

1. 'chr12:0-5851927'
2. 'chr12:529890-8218574'
3. 'chr12:970785-9438425'
4. 'chr12:3269770-12521686'
5. 'chr12:5556863-13471233'
6. 'chr12:6812512-13829975'

The most updated version of mvSuSiE multi-gene results are under path s3://statfungen/ftp_fgc_xqtl/analysis_result/mvsusie_multi_gene_test/multi_gene/ Currently it is still WIP. You can revisit this later when we prompt you to. Here is an example for C1S:

mnm_gene <- list()
for (window in sliding_windows) {
    mnm_gene_tmp <- NULL
    mnm_gene_tmp <- tryCatch(
        readRDS(paste0('/data/analysis_result/mvsusie_multi_gene/multi_gene/ROSMAP_multi_gene.', window, '.mnm.rds')),
        error = function(e) NULL
    )
    
    if (!is.null(mnm_gene_tmp)) {
        if(target_gene_info$gene_info$region_id %in% mnm_gene_tmp$mvsusie_fitted$condition_names){
        tryCatch({
            p <- mvsusieR::mvsusie_plot(mnm_gene_tmp$mvsusie_fitted, sentinel_only = F, add_cs = T)
            print(p)  # This ensures the plot is displayed in JupyterLab
        }, error = function(e) NULL)
        } else {
            message('There is mnm result for sliding window ',window,', but not include target gene ', gene_name, ' in CS')
        }
        mnm_gene <- append(mnm_gene, list(mnm_gene_tmp))
    }
}

In this case, there is no statistical evidence for C1S sharing any of its xQTL with other genes in ROSMAP Microglia data we looked into; although we have not analyzed MiGA this way yet (which showed some potential signals from the quick analysis above).

Section 7: Epigenomic QTL and their target regions

fsusie, see notebook

Generate a crude plot to determined whether the story is interesting

This is a crude version of the case study plot which shows the fsusie Effect (colored line), the gene body (black arrow), the epi-QTL (large dots with the same color as the effects) and ADGWAS cs position (small red dots).

Only produce the refine plot if we can see either:

1. There are sharing snp between epi-QTL and AD CS
2. There are the AD CS located within one of the effect range
3. The crude plot suggest something interesting

options(repr.plot.width = 40, repr.plot.height = 40)

 ggplot() + theme_bw() + facet_grid(cs_coverage_0.95 + study + region ~ ., labeller = labeller(.rows = function(x) gsub("([_:,-])", "\n", x)), scale = "free_y") +

      theme(text = element_text(size = 20), strip.text.y = element_text(size = 25, angle = 0.5)) +
     # xlim(view_win) +
      ylab("Estimated effect") +
   #   geom_line(data = haQTL_df %>% mutate(study = "haQTL effect") %>% filter(CS == 5),
    #            aes_string(y = "fun_plot", x = "x", col = "CS"), size = 4, col = "#00AEEF") +
  geom_line(data = effect_of_interest ,
                aes_string(y = "fun_plot", x = "x", col = "cs_coverage_0.95"), size = 2) +  
    geom_point(data = effect_of_interest ,
                aes_string(y = "pip", x = "pos", col = "cs_coverage_0.95"), size = 10) +
    theme(text = element_text(size = 40), strip.text.y = element_text(size = 15, angle = 0.5), 
            axis.text.x = element_text(size = 40), axis.title.x = element_text(size = 40)) +
      xlab("Position") +
      ylab("Estimated\neffect") +
      geom_segment(arrow = arrow(length = unit(1, "cm")), aes(x = gene_start, xend = gene_end, y = 1, yend = 1), size = 6,
                  data = tar_gene_info$gene_info, alpha = 0.3) +
      geom_text(aes(x = (gene_start + gene_end) / 2, y = 1 , label = gene_name), size = 10, 
              data = tar_gene_info$gene_info)+
        geom_point(aes(x = pos, y = pip  ) ,color = "red", data = flatten_table%>%filter( str_detect(study,"AD_") , cs_coverage_0.95 != 0  )%>%mutate(AD_study = study%>%str_replace_all("_","\n" ))%>%select(-study,-region,-cs_coverage_0.95) ) 

Section 8: Context focused validation in other xQTL data

see notebook

add fake version for now, so you don't have to refer to above link

Background: our "discovery set" is ROSMAP but we have additional "validation" sets including:

• STARNET
• MiGA
• KnightADRC
• MSBB
• metaBrain
• UKB pQTL

TODO:

• Get from Carlos WashU CSF based resource (pQTL and metabolomic QTL)

This section shows verification of findings from these data-sets. In principle we should check them through sections 1-6 more formally. In practice we will start with colocalization via colocboost --- since our study is genetics (variant and loci level) focused. We can selectively follow them up for potentially intereting validations. We therefore only demonstrate validation via colocboost as a starting point.

finempping_contexts <- readRDS(paste0(gene_name, '_finemapping_contexts.rds')) # from sec1

finempping_contexts <- get_norosmap_contexts(finempping_contexts)

cb_ad <- plot_cb(cb_res = cb_res, cex.pheno = 1.5, x.phen = -0.2, add_QTL = TRUE, cohorts = finempping_contexts, gene_id = gene_id)

No pvalue cutoff. Extract all variants names.

In conclusion from what's shown above, when we check the association signals in STARNET and MiGA on colocalization established from ROSMAP and AD GWAS, we see additional evidences.

Section 9: Non-linear effects of xQTL

see notebook

APOE interaction

options(repr.plot.width=6, repr.plot.height=6)

ggplot(C1S_int_res, aes(x = variant_id, y = qvalue_interaction)) +
  geom_point(alpha = 0.7, size = 6) +
  labs(title = "qvalue for C1S csets in interaction association nalysis",
       x = "Gene Name",
       y = "qvalue_interaction",
       size = "qvalue_interaction") +
  theme_minimal(base_size = 14) +
  theme(panel.background = element_blank(),
        panel.grid.major = element_line(color = "grey80"),
        legend.position = NULL,
        axis.text.x = element_text(angle = 45, hjust = 1))  + ylim(0,1)
  # scale_color_manual(values = colorRampPalette(brewer.pal(8, "Set1"))(length(unique(flat_var$gene_name))))
ggsave('plots/C1S/sec11.interaction_association_C1S_lessPIP25.pdf', height = 5, width = 8) 

In conclusion, there is no interaction QTL with APOE identified.

Section 10: in silico functional studies in iPSC model

see notebook

vars_p

apoe_p

Section 11: Functional annotations of selected loci

see notebook

TODO

• Touch base with Ryan on the snATAC annotations
• Run this by Pavel to see if there are additional comments on how we do this

func_p

Section 12: Candidate loci as trans-xQTL

see notebook

options(repr.plot.width=12, repr.plot.height=6)
if(!is.null(flat_var)){
   p =  ggplot(flat_var, aes(x = gene_name, y = pip, size = pip)) +
      geom_point(alpha = 0.7) +
      labs(title = paste0("PIP values for trans fine mapped Genes in ", gene_name ," csets with AD"),
           x = "Gene Name",
           y = "PIP",
           size = "PIP",
           color = "CS Coverage 0.95 Min Corr") +
      theme_minimal(base_size = 14) +
      theme(panel.background = element_blank(),
            panel.grid.major = element_line(color = "grey80"),
            legend.position = NULL,
            axis.text.x = element_text(angle = 45, hjust = 1))  
      # scale_color_manual(values = colorRampPalette(brewer.pal(8, "Set1"))(length(unique(flat_var$gene_name))))
    ggsave(paste0('plots/CR1/sec12.trans_fine_mapping_',gene_name,'.pdf'),p, height = 5, width = 8)
    p
    } else{
    message('There are no detectable trans signals for ', gene_name)
}

Creative thinking: generate hypothesis, search in literature, raise questions to discuss

You can now generate some preliminary hypotheses based on the above results. Next, you should search for evidence in the literature to support or refine these hypotheses and identify additional analyses needed to confirm them.