In kevinlkx/Mapgen: Multi-function software that performs enrichment, functionally-informed genetic fine-mapping and gene mapping

knitr::opts_chunk$set(echo = TRUE, eval=TRUE, comment = "#")

Partitioning PIPs by cell types and annotation categories

Here, we show an example using AFib GWAS finemapping result and scATAC-seq data from our heart study. We assigned the likely cell type(s) through which the causal variants act in each locus using fine-mapped SNPs, and cell-type specific open chromatin regions (OCRs) obtained from scATAC-seq data.

Required input data:

• Fine-mapped summary statistics.
• A list of GRanges objects of annotation regions.

Load R packages

library(mapgen)

Load fine-mapping summary statistics from the AFib study Keep SNPs with PIP > 1e-5, rename columns, and save as a GRanges object.

finemapstats <- readRDS(system.file("extdata", "AF.finemapping.sumstats.rds", package = "mapgen"))
finemapstats <- process_finemapping_sumstats(finemapstats, 
                                             snp = 'snp', chr = 'chr', 
                                             pos = 'pos', pip = 'susie_pip', 
                                             pval = 'pval', zscore = 'zscore', 
                                             cs = 'cs', locus = 'locus',  
                                             pip.thresh = 1e-5)

We can partition PIPs into different functional annotation categories.

Load genomic annotations (hg19).

genomic.annots <- readRDS(system.file("extdata", "genomic.annots.hg19.rds", package = "mapgen"))

Load cell-type specific open chromatin regions (OCRs) (hg19).

OCRs <- readRDS(system.file("extdata", "OCRs.hg19.gr.rds", package = "mapgen"))

Create a list of the annotations, with priority in the order of OCRs, UTRS, Exons, and Introns.

annots.list <- list(OCRs = OCRs,
                    UTRs = genomic.annots$UTRs,
                    Exons = genomic.annots$exons,
                    Introns = genomic.annots$introns)

Sum PIPs within annotation categories

Unlike \code{partition_pip_regions()}(below), it is OK to have overlapping annotations here. If a SNP is in multiple annotation categories, it will be assigned to the first ordered category.

sum.pip.res <- partition_pip_annots(finemapstats, annots.list)

Sum of PIPs in each annotation category:

sum.pips <- sum.pip.res$sum.pips
head(sum.pips)

Obtain a matrix of the proportion of PIPs in each annotation category.

locus.order <- rownames(sum.pips)[with(sum.pips, order(-OCRs, UTRs, Exons, Introns, others))]
sum.pips <- sum.pips[locus.order,]
prop.pip.mat <- sum.pips/rowSums(sum.pips)

head(prop.pip.mat)

We can make a structure plot to show the proportion of PIPs in each annotation category.

categories <- c("OCRs", "UTRs", "Exons", "Introns", "others")
colors <- c(OCRs = "#E18727FF", UTRs = "#238b45", Exons =  "#bee6af", Introns = "#B09C85FF", others = "#aaaaaa")
pip_structure_plot(prop.pip.mat, categories = categories, colors = colors)

Partition PIPs into disjoint OCRs for different cell types

We can partition PIPs into disjoint OCRs for different cell types.

Load a list of GRanges objects containing disjoint OCRs for different cell types.

disjoint.OCRs <- readRDS(system.file("extdata", "disjoint.OCRs.hg19.grlist.rds", package = "mapgen"))

Sum PIPs within cell-type specific OCRs.

sum.pip.res <- partition_pip_regions(finemapstats, disjoint.OCRs)

Sum of PIPs in each cell type OCR category:

sum.pips <- sum.pip.res$sum.pips
head(sum.pips)

Select loci with total PIPs in OCR > 0.25, and compute the proportion of PIPs partitioned in each cell type category.

# reorder the loci to match the previous figure
sum.pips <- sum.pips[locus.order, ]
# filter loci with total PIPs in OCR > 0.25
sum.pips.filtered <- sum.pips[rowSums(sum.pips) > 0.25,]
prop.pip.mat <- sum.pips.filtered/rowSums(sum.pips.filtered)

head(prop.pip.mat)

We can make a structure plot to show the proportion of PIPs in each cell type category.

categories <- c("Cardiomyocyte", "Endothelial", "Fibroblast", "Lymphoid", 
                "Myeloid", "Pericyte", "Shared 2-3", "Shared 4+")
colors <- c("#b22222", "#8DD3C7", "#BEBADA", "#FB8072", 
            "#80B1D3", "#B3DE69", "royalblue", "#003C86")
pip_structure_plot(prop.pip.mat, categories = categories, colors = colors)

kevinlkx/Mapgen documentation built on April 2, 2024, 5:48 a.m.