Allele-specific
accessibility
using ChrAccR
using ChrAccR
Irem B. Gunduz
2026-11-01
Source:vignettes/allelespec.Rmd
allelespec.Rmd
Introduction
At a heterozygous site, both alleles share the same nucleus, the same
transcription factors and the same sequencing library. If one allele is
covered by more ATAC-seq reads than the other, the difference can
therefore be attributed to the DNA sequence itself. This is referred to
as allele-specific chromatin accessibility (ASC). ChrAccR
stores allele-specific counts in a DsASC object, a subclass
of DsAcc that holds a reference and an alternative allele
count matrix (sites x samples) next to the usual peak accessibility
matrix.
This vignette covers counting alleles at heterozygous sites, cleaning the counts, testing for allelic imbalance, comparing imbalance between groups of samples and relating imbalance to predicted transcription factor binding. Steps that require a computing cluster are shown but not run. The scripts we used for these steps are available in the AlleleSpeC repository (https://github.com/EpigenomeInformatics/allele_spec).
In addition to the usual ChrAccR dependencies, the
allele-specific functions require the following packages:
BiocManager::install(c("VariantAnnotation", "Rsamtools", "HDF5Array", "DelayedArray",
"JASPAR2020", "TFBSTools", "motifmatchr",
"BSgenome.Hsapiens.UCSC.hg38"))
install.packages("ashr")First, let us load the package and prepare the plot style used in the remainder of this vignette.
Preparing the data
Genotyping and mapping-bias correction
Two steps are performed outside of ChrAccR before any
alleles are counted. First, genotypes are called for each donor, jointly
across all of that donor’s libraries, and only biallelic heterozygous
single-nucleotide variants are retained. Second, reference mapping bias
is removed. Reads carrying the alternative allele align slightly worse
than reads matching the reference genome, which shifts every site
towards the reference allele. WASP (van de Geijn et al. 2015)
corrects for this by remapping each read with its alleles swapped and
discarding reads that do not map to the same location. We recommend
checking the corrected libraries: the pooled reference allele fraction
should be close to 0.5, and values above about 0.55 indicate remaining
bias.
The sample annotation table should contain the columns
sampleId, bamFilename and donor,
as well as any columns used for grouping, such as cellType
and stimulus:
meta <- fread("sampleAnnot.csv")
meta[, donor := as.character(donor)]
# point to the WASP-corrected BAM files
meta[, bamFilename := file.path(waspDir, paste0(sampleId, "_wasp.bam"))]
meta <- meta[file.exists(bamFilename)]Building a shared site list
Each donor is heterozygous at a different set of sites. The
buildMasterSNPs function reads the VCF files of all donors
and returns a single sorted list of sites, including a
donorMembership column that records which donors are
heterozygous at each site. Counting all samples against this shared list
keeps the sites aligned across datasets.
masterGr <- buildMasterSNPs(donors = unique(meta$donor), vcfDir = vcfDir, genome = "hg38")
saveRDS(masterGr, file.path(ascDir, "masterSNPs.rds"))Counting alleles
The DsASC.gatk function counts reads supporting the
reference and the alternative allele at each site. If a sample’s donor
is not heterozygous at a site, the count is set to NA
rather than zero, so that homozygous sites are not mistaken for balanced
ones. For large datasets, we recommend counting one cell type at a time
and setting diskDump = TRUE, which writes the counts to
disk sample by sample and keeps memory usage low.
setConfigElement("minMapq", 20L)
setConfigElement("minBaseq", 20L)
setConfigElement("maxDepth", 5000L)
ct <- meta$cellType[1]
dsAsc <- DsASC.gatk(
sampleAnnot = as.data.frame(meta[cellType == ct,
.(sampleId, bamFilename, donor, cellType, stimulus)]),
vcfDir = vcfDir,
genome = "hg38",
diskDump = TRUE,
masterGr = masterGr
)Restricting to accessible regions
Allelic imbalance is only informative within accessible chromatin.
The following code removes sex chromosomes and lowly covered sites,
retains sites within consensus peaks and adds the peak accessibility
matrix from a DsATAC dataset:
dsAtac <- loadDsAcc(atacDir)
peakGr <- getCoord(dsAtac, "cons_peaks")
peakMat <- getCounts(dsAtac, "cons_peaks")
names(peakGr) <- paste0("peak_", seq_along(peakGr))
rownames(peakMat) <- names(peakGr)
# allelic balance is not defined on the sex chromosomes
dsAsc <- filterChroms(dsAsc, remove = c("chrX", "chrY", "chrM"))
# keep sites with at least 5 reads in at least one sample
dsAsc <- filterLowCovg(dsAsc, thresh = 5, reqSamples = 1)
# keep sites within consensus peaks
dsAsc <- annotateASCSites(dsAsc, peakGr)
dsAsc <- filterASCByPeaks(dsAsc)
dsAsc <- setAccessibility(dsAsc, peakMat[, getSamples(dsAsc)], peakGr)Loading the example data
The ChrAccRex data package contains a small example
dataset for allele-specific accessibility. If you haven’t already
installed ChrAccRex use
devtools::install_github("EpigenomeInformatics/ChrAccRex")The example data used in this vignette was designed for illustrating
the functionality of ChrAccR rather than a full analysis.
It contains allele-specific counts for two cell types, each profiled
resting (U) and stimulated (S), taken from the
human immune cell atlas of (Calderon and Pritchard 2019),
together with precomputed motif scores for the covered sites. It can be
loaded using:
example <- ChrAccRex::loadExample("dsAsc_ia_example")
dsDemo <- example$ds
deltaDemo <- example$delta
annotDemo <- as.data.table(getSampleAnnot(dsDemo))
dsDemo> DsASC allele-specific chromatin accessibility dataset
> contains:
> * 34 samples
> * 211 heterozygous SNP sites
> * accessibility: not loaded
> * allele assays: ref, alt
table(annotDemo$cellType, annotDemo$stimulus)>
> S U
> Bmem 8 8
> TeffNaive 10 8
DsASC data structure
The allele-specific count matrices can be retrieved using
getRefCounts and getAltCounts.
getAllelicBalance returns the fraction of reads supporting
the alternative allele, setting entries with low coverage to
NA:
> [1] 211 34> Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
> 0.000 0.375 0.488 0.482 0.588 1.000 5384Testing for allelic imbalance
Cleaning the counts
Before testing, ChrAccR provides functions for removing
common artifacts. ascDropHomozygous masks sites that were
called heterozygous but appear homozygous once all samples of a donor
are combined; such sites would otherwise show up as strong but false
imbalance. ascFilterRegions removes sites in blacklisted
regions, where genotyping and mapping-bias correction are less reliable.
ascPoolReplicates sums the counts of technical replicates,
increasing the number of reads per test. The example counts have already
been processed with these steps.
# mask sites that look homozygous within a donor
cleaned <- ascDropHomozygous(refMat, altMat, annotDemo,
minMinor = 2L, minMinorFrac = 0.05, minTotal = 10L)
# remove sites in blacklisted regions
keepIds <- ascFilterRegions(rownames(cleaned$ref), blacklistGr)
cleaned$ref <- cleaned$ref[keepIds, ]
cleaned$alt <- cleaned$alt[keepIds, ]
# combine technical replicates
pooled <- ascPoolReplicates(cleaned$ref, cleaned$alt, annotDemo,
groupCols = c("cellType", "stimulus", "donor"))Statistical testing
Using the calcASCStatistics function, we can test each
covered site in each sample for allelic imbalance. It performs a
two-sided binomial test against an expected allele fraction of 0.5 and
corrects the p-values within each sample. minAllele
requires a minimum number of reads on both alleles.
ascStats <- calcASCStatistics(dsDemo, minCoverage = 10, minAllele = 2)> 2026-11-01 14:57:36 2.0 STATUS STARTED Calculating ASC Statistics (binomial test vs p = 0.5)
> 2026-11-01 14:57:36 2.0 INFO Running two-sided binomial test against p = 0.5
> 2026-11-01 14:57:36 2.1 INFO Tested 1747 site-sample observations; 47 significant at FDR < 0.1.
> 2026-11-01 14:57:36 2.1 STATUS COMPLETED Calculating ASC Statistics (binomial test vs p = 0.5)
ascStats[, cellType := annotDemo$cellType[match(sampleId, annotDemo$sampleId)]]
# the most significant sites
ascStats[order(fdr)][1:5, .(snpId, sampleId, cellType, ref, alt, total, log2FC, fdr)]> snpId sampleId cellType ref alt total log2FC
> <char> <char> <char> <num> <num> <num> <num>
> 1: chr20:25390923_C_T Bmem_S_1002_b Bmem 6 45 51 2.716207
> 2: chr20:844525_A_G TeffNaive_U_1001_a TeffNaive 29 3 32 -2.906891
> 3: chr20:62938103_G_C Bmem_U_1001_a Bmem 24 2 26 -3.058894
> 4: chr20:2301403_A_C Bmem_U_1001_b Bmem 6 32 38 2.237039
> 5: chr20:844525_A_G TeffNaive_S_1001_b TeffNaive 21 3 24 -2.459432
> fdr
> <num>
> 1: 1.337625e-06
> 2: 3.450620e-04
> 3: 1.468658e-03
> 4: 3.578356e-03
> 5: 7.793037e-03The results can be visualized as a volcano plot, showing effect size against significance, and as a Manhattan plot along the genome:
plotASCVolcano(dsDemo, stats = ascStats, sigCut = 0.1)
plotASCManhattan(dsDemo, stats = ascStats, sigCut = 0.1)
For datasets with several donors, filterForRecurrence
retains sites that are significant in a minimum number of donors:
recurrent <- filterForRecurrence(ascStats, minDonors = 2, fdrCutoff = 0.1)Comparing groups of samples
To identify sites that are imbalanced in one group of samples but not in another, we combine the reads of each group and test each site once per group. Here, the two groups are the two cell types, but any grouping from the sample annotation can be used in the same way.
groupA <- annotDemo[cellType == "TeffNaive", sampleId]
groupB <- annotDemo[cellType == "Bmem", sampleId]
# binomial test on the reads combined across a group of samples
groupStats <- function(sids) {
r <- rowSums(refMat[, sids, drop = FALSE])
a <- rowSums(altMat[, sids, drop = FALSE])
n <- r + a
p <- pmin(1, 2 * pmin(pbinom(r, n, 0.5), pbinom(r - 1, n, 0.5, lower.tail = FALSE)))
data.table(snpId = rownames(refMat), refFrac = r / n, total = n, fdr = p.adjust(p, "BH"))
}
groups <- merge(groupStats(groupA), groupStats(groupB), by = "snpId", suffixes = c(".A", ".B"))
# imbalanced in one group, balanced in the other despite good coverage
siteA <- groups[total.B >= 20 & fdr.B > 0.5][order(fdr.A)][1, snpId]
siteB <- groups[total.A >= 20 & fdr.A > 0.5][order(fdr.B)][1, snpId]
groups[snpId %in% c(siteA, siteB)]> Key: <snpId>
> snpId refFrac.A total.A fdr.A refFrac.B total.B
> <char> <num> <num> <num> <num> <num>
> 1: chr20:35619758_C_G 0.2658228 79 0.000813577 0.4677419 124
> 2: chr20:63521166_C_G 0.4939759 83 1.000000000 0.7681159 69
> fdr.B
> <num>
> 1: 0.8338637190
> 2: 0.0001128858To visualize such a site, we combine the reads of each cell type and condition, matching the way the sites were tested, and plot the reference allele fraction with its exact binomial confidence interval:
plotGroupBalance <- function(snp) {
d <- annotDemo[, .(ref = sum(refMat[snp, sampleId], na.rm = TRUE),
alt = sum(altMat[snp, sampleId], na.rm = TRUE)),
by = .(cellType, stimulus)]
d <- d[ref + alt > 0]
ci <- t(mapply(function(r, n) binom.test(r, n)$conf.int, d$ref, d$ref + d$alt))
d[, `:=`(refFrac = ref / (ref + alt), lo = ci[, 1], hi = ci[, 2])]
ggplot(d, aes(x = refFrac, y = cellType, colour = stimulus)) +
geom_vline(xintercept = 0.5, linetype = "dashed", colour = "grey60") +
geom_errorbar(aes(xmin = lo, xmax = hi), width = 0,
position = position_dodge(width = 0.5)) +
geom_point(size = 2, position = position_dodge(width = 0.5)) +
scale_x_continuous(limits = c(0, 1)) +
labs(x = "Proportion of reads mapping to the reference allele", y = NULL, title = snp)
}
plotGroupBalance(siteA)
plotGroupBalance(siteB)
The per-library view of the same site, with one point for each
sample, can be obtained using plotASCBalance:
plotASCBalance(dsDemo, siteB, colorBy = "stimulus", orderBy = "cellType",
sigStats = ascStats, sigCut = 0.1)Sharing across cell types
Because all cell types were counted against the same site list, their
datasets can be combined using mergeDsASCArray. The
ascSharing function takes the sites that are significant in
one group and estimates, using ashr (Stephens
2017), the fraction of these sites that also show imbalance
in another group. estimateSharedImbalance performs the same
comparison directly on two DsASC objects.
merged <- mergeDsASCArray(ascDir)
aggA <- ascAggregate(refMat, altMat, groupA)
aggB <- ascAggregate(refMat, altMat, groupB)
sharing <- ascSharing(aggA$ref, aggA$alt, aggB$ref, aggB$alt, fdrA = 0.01, minReadsB = 4)
shared <- estimateSharedImbalance(dsGroupA, dsGroupB, fdrCutoff = 0.01)Allelic imbalance and transcription factor binding
If a variant affects transcription factor binding, the allele
preferred by the factor should also be the more accessible one. The
example dataset contains the difference in motif score between the
alternative and the reference allele (delta) for a set of
transcription factor motifs. We label each site by the allele its motif
prefers if the score changes by more than 1, and treat all other sites
as background:
prefDelta <- 1
prefDt <- data.table::copy(deltaDemo)
prefDt[, preference := fifelse(delta > prefDelta, "Alternative",
fifelse(delta < -prefDelta, "Reference", "No preference"))]Next, we choose a cell type and the conditions to compare, combine the reads of that cell type’s samples for each condition and remove sites with fewer than 4 reads:
focusCellType <- "TeffNaive"
conditions <- c(U = "Resting", S = "Stimulated")
focusSamples <- annotDemo[cellType == focusCellType, sampleId]
condFrac <- function(cond) {
sids <- annotDemo[stimulus == cond & sampleId %in% focusSamples, sampleId]
r <- rowSums(refMat[, sids, drop = FALSE])
a <- rowSums(altMat[, sids, drop = FALSE])
data.table(snpId = rownames(refMat), refFrac = r / (r + a), total = r + a,
condition = conditions[[cond]])
}
fracDt <- rbindlist(lapply(names(conditions), condFrac))
fracDt <- fracDt[total >= 4]
plotDt <- merge(fracDt, prefDt, by = "snpId", allow.cartesian = TRUE)
plotDt[, preference := factor(preference, levels = c("Alternative", "No preference", "Reference"))]
plotDt[, condition := factor(condition, levels = unname(conditions))]Sites at which the motif prefers the reference allele are expected to show a reference allele fraction above 0.5, and sites at which it prefers the alternative allele a fraction below 0.5:
ggplot(plotDt, aes(x = condition, y = refFrac, fill = preference)) +
geom_hline(yintercept = 0.5, linetype = "dotted") +
geom_boxplot(outlier.size = 0.3, position = position_dodge(width = 0.8)) +
facet_wrap(~ tf, ncol = 3) +
scale_fill_manual(values = c("Alternative" = "#1B9E77",
"No preference" = "#D95F02",
"Reference" = "#7570B3"), name = "Binding preference") +
scale_y_continuous(limits = c(0, 1)) +
labs(x = NULL, y = "Proportion of reads mapping to the reference allele") +
theme(legend.position = "bottom")
The difference between the two groups can be tested for each motif and condition using a Wilcoxon rank-sum test:
plotDt[preference %in% c("Reference", "Alternative"),
.(n_ref = sum(preference == "Reference"),
n_alt = sum(preference == "Alternative"),
p = tryCatch(wilcox.test(refFrac ~ droplevels(preference))$p.value,
error = function(e) NA_real_)),
by = .(tf, condition)][order(p)]> tf condition n_ref n_alt p
> <char> <fctr> <int> <int> <num>
> 1: TCF7 Resting 19 27 0.002239061
> 2: BACH1 Stimulated 29 41 0.005364915
> 3: CTCF Resting 62 56 0.006602805
> 4: GATA3 Stimulated 24 17 0.011064639
> 5: FOSL1 Resting 22 32 0.017051864
> 6: JUNB Resting 22 29 0.032372026
> 7: BACH1 Resting 29 41 0.054122890
> 8: JUNB Stimulated 22 29 0.148286599
> 9: IRF4 Resting 50 36 0.189144576
> 10: SPI1 Stimulated 40 49 0.193889109
> 11: CTCF Stimulated 62 57 0.211255822
> 12: FOSL1 Stimulated 22 32 0.224424321
> 13: NFKB1 Stimulated 49 45 0.289262565
> 14: SPI1 Resting 40 49 0.338661670
> 15: TCF7 Stimulated 19 27 0.468244592
> 16: IRF4 Stimulated 50 36 0.671143241
> 17: GATA3 Resting 24 17 0.894726793
> 18: NFKB1 Resting 49 44 0.990790378Session Info
Sys.Date()> [1] "2026-11-01"> R version 4.4.1 (2024-06-14)
> Platform: x86_64-conda-linux-gnu
> Running under: Debian GNU/Linux 12 (bookworm)
>
> Matrix products: default
> BLAS/LAPACK: /icbb/projects/share/software/packages/miniconda3/envs/alleleSpec/lib/libopenblasp-r0.3.30.so; LAPACK version 3.12.0
>
> locale:
> [1] C
>
> time zone: Europe/Berlin
> tzcode source: system (glibc)
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> other attached packages:
> [1] ggplot2_4.0.1 data.table_1.17.8 ChrAccR_0.9.26
>
> loaded via a namespace (and not attached):
> [1] SummarizedExperiment_1.36.0 gtable_0.3.6
> [3] xfun_0.54 bslib_0.9.0
> [5] htmlwidgets_1.6.4 rhdf5_2.50.0
> [7] Biobase_2.66.0 lattice_0.22-7
> [9] generics_0.1.4 rhdf5filters_1.18.0
> [11] vctrs_0.6.5 tools_4.4.1
> [13] bitops_1.0-9 ChrAccRex_0.4
> [15] stats4_4.4.1 parallel_4.4.1
> [17] tibble_3.3.0 pkgconfig_2.0.3
> [19] Matrix_1.7-4 RColorBrewer_1.1-3
> [21] S7_0.2.1 desc_1.4.3
> [23] S4Vectors_0.44.0 lifecycle_1.0.4
> [25] GenomeInfoDbData_1.2.13 compiler_4.4.1
> [27] farver_2.1.2 Rsamtools_2.22.0
> [29] textshaping_1.0.1 Biostrings_2.74.0
> [31] codetools_0.2-20 muLogR_0.2
> [33] GenomeInfoDb_1.42.0 muRtools_0.9.6
> [35] htmltools_0.5.8.1 sass_0.4.10
> [37] yaml_2.3.11 pillar_1.11.1
> [39] pkgdown_2.2.0 crayon_1.5.3
> [41] jquerylib_0.1.4 BiocParallel_1.40.0
> [43] DelayedArray_0.32.0 cachem_1.1.0
> [45] abind_1.4-8 mime_0.13
> [47] tidyselect_1.2.1 digest_0.6.39
> [49] dplyr_1.1.4 labeling_0.4.3
> [51] fastmap_1.2.0 grid_4.4.1
> [53] cli_3.6.5 SparseArray_1.6.0
> [55] magrittr_2.0.4 S4Arrays_1.6.0
> [57] dichromat_2.0-0.1 withr_3.0.2
> [59] UCSC.utils_1.2.0 scales_1.4.0
> [61] rmarkdown_2.30 XVector_0.46.0
> [63] httr_1.4.7 matrixStats_1.5.0
> [65] ragg_1.5.0 HDF5Array_1.34.0
> [67] evaluate_1.0.5 knitr_1.50
> [69] GenomicRanges_1.58.0 IRanges_2.40.0
> [71] rlang_1.1.6 glue_1.8.0
> [73] BiocGenerics_0.52.0 jsonlite_2.0.0
> [75] R6_2.6.1 Rhdf5lib_1.28.0
> [77] MatrixGenerics_1.18.0 GenomicAlignments_1.42.0
> [79] systemfonts_1.3.1 fs_1.6.6
> [81] zlibbioc_1.52.0