Ask the question the whole series has been building toward: in which immune cells does the inherited risk for critical COVID-19 sit? Fine-map the COVID-19 Host Genetics Initiative’s risk loci, test whether the likely causal variants fall in each cell type’s open chromatin, and score every cell and every donor with gchromVAR.
Every analysis since Part 5 has compared severe and healthy donors. Part 5 asked which peaks change, Part 6 which transcription factors leave deeper footprints, and Part 7 which motifs gain accessibility genome-wide. Those are questions about the consequences of disease: what severe COVID-19 does to chromatin in blood.
This tutorial asks a question about causes. People differ in their risk of becoming critically ill with COVID-19, and part of that difference is inherited. Genome-wide association studies (GWAS) have mapped dozens of risk loci, but most of the associated variants lie outside genes, in DNA whose job is regulation. A risk variant in an enhancer can only act in the cells where that enhancer is open. So if we know where the risk variants are, and we know which regions are open in which cell types, we can ask: which cell types’ regulatory DNA carries the inherited risk?
That question needs no RNA, no matched samples and no disease-state comparison. Inherited variants are the same in every cell of every donor, so the healthy donors’ cells answer it just as well as the severe donors’. It is the last question in this series, and the one that connects a single-cell chromatin dataset to human genetics.
This tutorial does four things:
- Fine-map the risk loci. Turn the COVID-19 Host Genetics Initiative’s summary statistics into loci and, within each locus, into posterior probabilities that each variant is the causal one.
- Overlap them with cell-type peaks. Count how much of that probability falls inside the peaks each cell type called in Part 8.
- Test the overlap against chance. Use permutations that respect the way GWAS variants cluster together.
- Score every cell and every donor with gchromVAR, the GWAS generalization of the chromVAR method you used in Part 7, and test cell types with donors as replicates.
Prerequisites: This tutorial assumes you have completed Part 1: From FASTQ to Peaks through Part 4: Cell Type Identification, and Part 8: Cell-Type-Specific Peak Recalling with MACS3 and Co-Accessibility with Cicero. You need two objects saved in Part 8,
unified_macs3_peaks.rdsandatac_macs3_peaks.rds, the four fragments files they point to (Part 1), and the same Pixi workspace used through Part 8. Parts 5 to 7 are not required, but the motif analysis in Step 16 ties back to the C/EBP findings of Part 6 and Part 7. Apart from installing three packages and downloading one file, everything runs in R.
🧭 Introduction: Why Bring GWAS Into a scATAC-seq Analysis?
From Risk Loci to Regulatory Elements
A GWAS tests millions of common genetic variants, one at a time, for association with a trait. For each variant it reports an effect size and a p-value. A risk locus is a region where variants pass the genome-wide significance threshold of p < 5 x 10^-8.
Two facts make GWAS results hard to interpret on their own:
- Most associated variants are non-coding. They do not change a protein. If they act at all, they usually act by changing how much of a gene is made, through enhancers and promoters — the regions an ATAC-seq experiment measures.
- Associated variants come in crowds. Variants that sit close together are inherited together, a correlation called linkage disequilibrium (LD). A single causal variant drags dozens or hundreds of correlated neighbors to significance with it. The variant with the smallest p-value is often, but not always, the causal one.
Chromatin accessibility helps with both problems. If a candidate variant sits in a peak that is open only in monocytes, it is a better candidate than its neighbors in closed chromatin, and monocytes become a candidate cell type for the disease mechanism. Done across all loci at once, that reasoning becomes a statistical test of which cell types are enriched for risk.
The COVID-19 Host Genetics Initiative, Release 7
The COVID-19 Host Genetics Initiative (HGI) is an international collaboration that meta-analyzes genome-wide association studies of COVID-19. Its data release 7 analyzed three phenotypes; the most severe, critical illness (respiratory support or death), included 21,194 cases and identified 30 genome-wide significant loci (COVID-19 Host Genetics Initiative, 2023).
We use the release 7 file for phenotype A2 — in the HGI’s words, very severe respiratory confirmed COVID-19 versus the population — across all ancestries, in the version that leaves out the 23andMe cohort, because that is the version released publicly with full per-variant results. The release 7 results page lists 18,152 cases and 1,145,546 controls from 34 contributing studies for this analysis, combined by fixed-effects inverse-variance meta-analysis.
The HGI distributes this analysis in several forms. We use the GRCh38 filtered file, which keeps only the variants with p < 1 x 10^-5. It is small (about 1.3 MB) and in the same genome build as our data.
Reading the Summary Statistics File
The column names follow the HGI’s meta-analysis software. The descriptions below are the HGI’s own, from the release 7 results page:
| Column | Meaning |
|---|---|
#CHR, POS | Chromosome and position (GRCh38). Chromosome X is coded as 23 |
REF, ALT | Reference allele, and the alternative allele, which is the effect allele: the effect size is for ALT |
SNP | #CHR:POS:REF:ALT, a unique variant identifier |
all_meta_N | Number of studies that had the variant after quality filtering — not a sample size |
all_inv_var_meta_beta | Effect size on the log odds ratio scale |
all_inv_var_meta_sebeta | Standard error of the effect size |
all_inv_var_meta_p | p-value |
all_inv_var_meta_cases, _controls, _effective | Total cases, total controls and effective sample size for this variant |
all_inv_var_het_p | p-value of Cochran’s Q test for heterogeneity between studies |
lmso_inv_var_beta, _se, _pval | The same meta-analysis with the most significant study left out |
all_meta_AF | Allele frequency in the meta-analysis |
rsid | dbSNP identifier, where one exists |
Beginner’s Tip: a positive
betameans theALTallele raises the odds of critical illness, a negativebetathat it lowers them. The odds ratio isexp(beta): abetaof 0.73 is an odds ratio of about 2.1.
What Is Fine-Mapping?
A GWAS finds regions, not variants. Variants that sit close together on a chromosome tend to be inherited together, a correlation called linkage disequilibrium (LD). So when one variant truly affects disease risk, dozens or hundreds of its neighbors show almost the same association, and the GWAS cannot tell them apart on p-values alone.
Statistical fine-mapping is the general name for methods that take an associated region and estimate, for each variant in it, how likely it is to be the causal one (Schaid et al., 2018). The HGI summary statistics are not fine-mapped. We fine-map them ourselves, in Step 9.
Five terms are all you need to follow the rest of this tutorial:
| Term | Meaning |
|---|---|
| Locus | A region containing one association signal. We define loci in Step 8 |
| Lead variant | The variant with the smallest p-value in a locus |
| LD | Correlation between nearby variants, which makes them carry nearly the same association |
| Posterior inclusion probability (PIP) | The estimated probability that a variant is the causal one. PIPs add up to 1 within a locus |
| 95% credible set | The smallest group of variants whose PIPs add up to at least 0.95. One variant means the locus is resolved; hundreds mean LD has left it unresolved |
The method used here is the simplest one that works from summary statistics alone: Wakefield’s approximate Bayes factor (ABF). It assumes one causal variant per locus, and under that assumption needs only each variant’s effect size and standard error — no LD information. Where a locus holds two independent signals, it gives almost all the probability to the stronger one. FINEMAP and SuSiE allow several causal variants, but need an LD matrix from the same population as the GWAS, which a 34-study, multi-ancestry meta-analysis does not provide. Treat the PIPs in this tutorial as a ranking within each locus, not as exact probabilities.
What Kind of GWAS File Does This Workflow Need?
Nothing about the HGI file is special. It is a standard table of GWAS summary statistics: one row per variant, with its position, its two alleles, an effect size, a standard error and a p-value. The HGI applied two kinds of processing before release — quality filters within each study, and the meta-analysis that combined the 34 studies — and the file we download keeps only variants with p < 1 x 10^-5. No fine-mapping, annotation or cell-type information has been added.
So the workflow applies to any GWAS for which you have these seven pieces of information per variant:
| Needed | In the HGI file | In a GWAS Catalog file |
|---|---|---|
| Chromosome | #CHR | chromosome |
| Position (same genome build as your peaks) | POS | base_pair_location |
| Effect allele | ALT | effect_allele |
| Other allele | REF | other_allele |
| Effect size (beta, or an odds ratio to convert) | all_inv_var_meta_beta | beta or odds_ratio |
| Standard error of the effect | all_inv_var_meta_sebeta | standard_error (or ci_lower / ci_upper) |
| p-value | all_inv_var_meta_p | p_value |
Everything else in the HGI file — the number of studies, the heterogeneity test, the leave-most-significant-study-out results — is used only for optional quality checks. The section Running This Workflow on Your Own GWAS Results, after the analysis, shows how to read a file in the GWAS Catalog’s standard format, how to handle genome-wide and genome-build differences, and how to start from fine-mapping results you already have.
Three Ways to Ask “Which Cell Type?”
- Overlap. Sum the PIPs of the credible-set variants that fall inside each cell type’s peaks.
- Permutation enrichment. Ask whether that sum is larger than expected by chance, by moving the variants to random positions many times with regioneR (Gel et al., 2016) and recomputing it.
- gchromVAR. Score how accessible the risk-weighted peaks are in each sample — each donor pseudobulk, or each single cell — relative to matched background peaks (Ulirsch et al., 2019).
The first two treat a peak as open or closed in a cell type. The third uses how strongly each sample’s chromatin is open at the risk regions, and it is the one that gives a score per donor, which makes a replicate-level test possible.
How gchromVAR Differs From Part 7’s chromVAR
gchromVAR is a generalization of chromVAR (Ulirsch et al., 2019), so most of what Part 7 explained still applies: each sample’s accessibility at a set of peaks is compared with its expected accessibility, and the difference is scaled against GC- and accessibility-matched background peaks to give a bias-corrected z-score. What changes is the input and the question.
| Part 7: chromVAR | Part 9: gchromVAR | |
|---|---|---|
| Peak annotation | Binary: does the peak contain the motif? | Continuous: summed PIPs of the fine-mapped variants inside the peak |
| Annotations tested | 633 JASPAR2020 motifs | One trait: COVID-19 critical illness |
| Samples scored | Four donor pseudobulks of CD14 monocytes | Donor pseudobulks of every cell-type label, then every single cell |
| Question | Which TFs’ targets change accessibility in severe disease? | Which cells’ open chromatin holds the inherited risk variants? |
| Downstream test | limma, borrowing variance across 633 motifs | One trait leaves nothing to borrow across: a t-test on donor-level scores |
Three consequences are worth spelling out:
- The weights are the point. chromVAR’s
computeDeviations()accepts only yes/no annotations. A peak holding one variant with PIP 0.9 should count far more than a peak holding one of 40 variants in tight LD, each at PIP 0.02. gchromVAR’scomputeWeightedDeviations()uses the PIPs as weights; a binary version would treat those two peaks alike. - Per-cell scores are noisy. One trait scored on sparse single cells gives noisy z-scores. That is why we test on donor pseudobulks, and why SCAVENGE (Yu et al., 2022) adds network propagation on top of gchromVAR for single-cell work.
- Two Part 7 settings need a fresh decision. Part 7 resized peaks to 300 bp; here we keep the Part 8 peaks at their native width, so that a variant near the edge of a peak is not dropped. And Part 7 used 2,000 background sets on the advice of a benchmark of TF activity methods. We keep 2,000 for consistency, but that benchmark did not test GWAS weights, so the number is a choice, not a validated setting.
Why Not SCAVENGE?
SCAVENGE (Yu et al., 2022) was built for exactly this kind of data, so it is worth saying why this tutorial does not use it. SCAVENGE is not a different scoring method: its first stage is gchromVAR in every cell, which is our Step 15. Its second stage takes the top-scoring cells as seeds and spreads their signal across a nearest-neighbor graph of all cells, so that relevant cells missed because of sparse data still receive a high trait relevance score.
We stop after the first stage for two reasons:
- This tutorial asks which cell types carry risk, with donors as replicates. Pseudobulk gchromVAR answers that directly (Step 14). SCAVENGE’s scores are per cell, and propagation deliberately shares signal between neighboring cells, so they are not independent observations that a donor-level test can rest on.
- It adds choices that need their own justification: the number of neighbors, the embedding and batch correction behind the graph, and the fraction of seed cells all change the scores.
SCAVENGE is most appropriate when the question is about cell states rather than cell types: which cells within a population carry the trait relevance, or along which part of a differentiation trajectory. That is how Yu et al. separated risk-enriched from risk-depleted CD14 monocytes, and it would be the natural way to sharpen the noisy per-cell UMAP in Step 20. Use it for exploration and display, and keep any statistical test at the level of donors.
What This Approach Has Revealed: Published Examples
Combining GWAS with single-cell chromatin data has answered real biological questions:
- Blood cell traits. Ulirsch et al. (2019) introduced g-chromVAR while fine-mapping 16 blood-cell traits in about 115,000 UK Biobank participants. Enrichment was strongest in lineage-restricted populations but also appeared in progenitors, and in single-cell data, progenitor populations split into trait-enriched subpopulations that differed in master-regulator activity.
- COVID-19 severity. Yu et al. (2022) developed SCAVENGE and applied it to PBMC scATAC-seq from COVID-19 and healthy donors. COVID-19 severity variants were enriched in monocytes and dendritic cells, and within CD14 monocytes the method separated risk-enriched from risk-depleted cell states. This is the closest precedent to this tutorial: same disease, same tissue, same assay.
- Type 1 diabetes. Chiou et al. (2021) profiled 131,554 nuclei from pancreas and blood. Risk was most enriched in T cell regulatory elements, and also in classical monocytes and plasmacytoid dendritic cells. Unexpectedly, several signals mapped to exocrine pancreas: a risk variant in a ductal-specific CFTR enhancer reduced enhancer activity and CFTR expression in ductal cells.
- Alzheimer’s and Parkinson’s diseases. Corces et al. (2020) profiled about 70,600 brain nuclei. Alzheimer’s heritability was enriched in microglial peaks, and machine-learning variant scoring nominated specific candidate causal variants, such as a microglia-specific variant in a SPI1 motif.
Most of these studies had genome-wide summary statistics or their own fine-mapping. Ours is a lighter analysis by design, and the conclusions should be correspondingly modest.
What We Expect to See in This Dataset
Two published results set expectations before we run anything.
- SCAVENGE found COVID-19 severity variants enriched in monocytes and dendritic cells of PBMC scATAC-seq. If our data agree, CD14 monocytes and dendritic cells should rank high. Dendritic cells are a caution here: Part 8 called their peaks from 87 nuclei, and those peaks were by far the weakest of the ten cell types.
- The strongest association in the file is at 3p21.31, and its lead variant, rs17713054, is the variant that Downes et al. (2021) identified as the probable causal variant there. Its risk A allele creates a second CEBPB motif in an enhancer that upregulates LZTFL1. The enhancer is accessible in lung epithelial and endothelial cells; Downes et al. found none of the 28 associated variants at 3p21.31 in open chromatin across 24 immune cell populations, and reported the enhancer as not accessible in primary CD14+ monocytes and CD4+ T cells. They concluded that pulmonary epithelial cells, rather than immune cells, are likely responsible for the 3p21.31 risk.
The second point has two consequences for us. The strongest risk locus may well fall outside every PBMC peak — which would be an expected result, not a failure. (It does not, and Steps 10, 12 and 22 show what it does instead: a useful reminder that expectations are hypotheses too.) And the transcription factor whose motif the risk allele creates, CEBPB, is the same factor that showed the largest footprint change in severe-disease monocytes in Part 6, and whose C/EBP family ranked near the top of Part 7’s genome-wide test (without reaching significance). Step 16 checks the motif gain directly, and Step 22 looks at whether the enhancer is open in any of our cell types.
The Part 9 Workflow
COVID-19 HGI release 7, A2 (p < 1e-5) unified_macs3_peaks.rds (Part 8)
| atac_macs3_peaks.rds (Part 8)
v |
[1] Read, check build against hg38, QC Steps 5-7
| |
v |
[2] Loci: greedy clumping, lmso check Step 8 |
| |
v |
[3] Wakefield ABF -> PIPs, 95% credible sets Step 9|
| |
+--------------------+---------------------+
|
v
[4] PIP of credible-set variants in each cell type's peaks Step 10
[5] Permutation tests: circular vs independent nulls Step 11
[6] Specificity (private peaks), per-locus contributions, Step 12
local z-score
|
v
[7] gchromVAR on donor-by-cell-type pseudobulks Step 13
[8] Cell types tested with donors as replicates Step 14
[9] gchromVAR in every cell Step 15
|
v
[10] REF vs ALT motif scores for candidate variants Step 16
|
v
[11] Visualization: Manhattan, credible sets, enrichment, Steps 17-22
UMAP, donors, coverage at two loci
⚙️ Setting Up the Analysis Environment
What Is New in Part 9
Most of what this tutorial needs is already in the Pixi workspace: Signac and Seurat throughout, chromVAR, motifmatchr, TFBSTools, JASPAR2020 and BSgenome.Hsapiens.UCSC.hg38 from Parts 4 and 7, and ggrepel from Part 7. Three packages are new:
| Package | Source | Used for |
|---|---|---|
bioconductor-regioner | Bioconda | Permutation tests on genomic regions (Steps 11 and 12) |
r-remotes | conda-forge | Installing gchromVAR from GitHub |
gchromVAR | GitHub | PIP-weighted deviations (Steps 13 and 15) |
#-----------------------------------------------
# STEP 1: Add regioneR and gchromVAR to the Pixi workspace
#-----------------------------------------------
# Run from a LOGIN node (package downloads need internet access)
cd /projects/mylab/shared/scatac-analysis
pixi add bioconductor-regioner r-remotes
# gchromVAR from GitHub, into the workspace's R library.
# upgrade = "never" and dependencies = FALSE leave every Pixi-managed package untouched.
pixi run Rscript -e 'remotes::install_github("caleblareau/gchromVAR", upgrade = "never", dependencies = FALSE)'
# Confirm every package this tutorial loads is visible to R
pixi run Rscript -e 'sapply(c("Signac", "regioneR", "gchromVAR", "chromVAR", "motifmatchr", "JASPAR2020", "TFBSTools", "BSgenome.Hsapiens.UCSC.hg38"), requireNamespace, quietly = TRUE)'
Expected output: all eight package names followed by TRUE.
Downloading the GWAS Summary Statistics
The filtered A2 file is a plain tab-separated text file on the HGI’s public Google Cloud bucket. We keep it next to the other reference files from Part 2, so every later tutorial or project can reuse it.
#-----------------------------------------------
# STEP 2: Download the COVID-19 HGI release 7 A2 summary statistics (p < 1e-5)
#-----------------------------------------------
# Run from a LOGIN node (compute nodes often have no internet access)
mkdir -p /projects/mylab/shared/reference/gwas
cd /projects/mylab/shared/reference/gwas
curl -O https://storage.googleapis.com/covid19-hg-public/freeze_7/results/20220403/main/sumstats/COVID19_HGI_A2_ALL_leave_23andme_20220403_1e-5.tsv
# One header line plus one line per variant
wc -l COVID19_HGI_A2_ALL_leave_23andme_20220403_1e-5.tsv
Output:
8307 COVID19_HGI_A2_ALL_leave_23andme_20220403_1e-5.tsv
Reading this output. One header line plus 8,306 variants, all with p < 1 x 10^-5. That is all the shell work; everything from here on runs in one R session.
Loading Libraries and Creating the Part 9 Directory Structure
Launch R from the workspace on a compute node with several CPUs. Steps 13 and 15 recount nothing from the fragments files, but they hold the Part 8 count matrix and a 2,000-column background-peak matrix in memory, and the permutations in Steps 11 and 12 run on parallel workers that each carry a copy of the session. We suggest at least 4 CPUs and 32 GB of memory.
cd /projects/mylab/shared/scatac-analysis
pixi run R
#-----------------------------------------------
# STEP 3: Load packages and set options
#-----------------------------------------------
library(Signac)
library(Seurat)
library(GenomicRanges)
library(GenomeInfoDb) # seqlengths()
library(BSgenome.Hsapiens.UCSC.hg38) # reference sequence: build check, GC content, motifs
library(regioneR) # permutation tests on genomic regions
library(chromVAR) # GC bias and background peaks
library(gchromVAR) # PIP-weighted deviations
library(motifmatchr) # motif scores for each allele
library(JASPAR2020)
library(TFBSTools)
library(SummarizedExperiment)
library(Matrix)
library(BiocParallel)
library(ggplot2)
library(ggrepel)
library(patchwork)
# L'Ecuyer-CMRG gives every parallel worker its own reproducible random stream,
# so the permutations in Steps 11 and 12 come out the same on every run
RNGkind("L'Ecuyer-CMRG")
set.seed(1234)
register(MulticoreParam(workers = 4)) # computeWeightedDeviations() in Steps 13 and 15
options(mc.cores = 4) # permTest() in Steps 11 and 12
#-----------------------------------------------
# STEP 4: Define and create the Part 9 directories
#-----------------------------------------------
base_dir <- "/projects/mylab/shared/scatac_tutorial"
recall_dir <- file.path(base_dir, "peak_recalling") # Part 8 output
gwas_dir <- file.path(base_dir, "gwas_enrichment") # this tutorial
obj_dir <- file.path(gwas_dir, "objects")
plot_dir <- file.path(gwas_dir, "plots")
table_dir <- file.path(gwas_dir, "tables")
invisible(lapply(c(obj_dir, plot_dir, table_dir),
dir.create, recursive = TRUE, showWarnings = FALSE))
# Downloaded in Step 2
gwas_file <- "/projects/mylab/shared/reference/gwas/COVID19_HGI_A2_ALL_leave_23andme_20220403_1e-5.tsv"
Your project tree now looks like this:
/projects/mylab/shared/
├── reference/
│ ├── blacklist/ # ENCODE blacklist (Part 2)
│ └── gwas/ # HGI A2 summary statistics (this tutorial)
└── scatac_tutorial/
├── results/ # Cell Ranger ATAC output (Part 1) <- fragments files
├── downstream_qc/ # Part 2
├── integrated/ # Part 3
├── annotation/ # Part 4
├── differential/ # Part 5
├── footprinting/ # Part 6
├── tf_activity/ # Part 7
├── peak_recalling/ # Part 8 <- input
└── gwas_enrichment/ # this tutorial
├── objects/
├── plots/
└── tables/
📥 Loading and Validating the Summary Statistics
Step 5: Reading the File
#-----------------------------------------------
# STEP 5: Read the summary statistics
#-----------------------------------------------
# check.names = FALSE keeps "#CHR" as it is instead of turning it into "X.CHR"
raw <- read.delim(gwas_file, check.names = FALSE)
dim(raw)
colnames(raw)
range(raw$all_inv_var_meta_p)
table(raw$all_meta_N >= 30)
colSums(is.na(raw[, c("rsid", "all_inv_var_het_p", "lmso_inv_var_pval", "all_meta_AF")]))
Output:
[1] 8306 18
[1] "#CHR" "POS"
[3] "REF" "ALT"
[5] "SNP" "all_meta_N"
[7] "all_inv_var_meta_beta" "all_inv_var_meta_sebeta"
[9] "all_inv_var_meta_p" "all_inv_var_meta_cases"
[11] "all_inv_var_meta_controls" "all_inv_var_meta_effective"
[13] "all_inv_var_het_p" "lmso_inv_var_beta"
[15] "lmso_inv_var_se" "lmso_inv_var_pval"
[17] "all_meta_AF" "rsid"
[1] 6.1225e-220 9.9963e-06
FALSE TRUE
3011 5295
rsid all_inv_var_het_p lmso_inv_var_pval all_meta_AF
112 6 4073 17
Reading this output.
- 8,306 variants, with p-values from 6 x 10^-220 up to the 1 x 10^-5 filter.
- About a third of the variants are missing from some studies. A variant typed in few studies rests on fewer cases, which matters in Step 8.
- The leave-most-significant-study-out columns are filled only for genome-wide significant variants, the ones Step 8 checks. Some variants have no rsid, so we never join on rsid.
Step 6: Checking the Genome Build Against the Reference Sequence
Our peaks are in hg38. If the GWAS coordinates were in hg19 (GRCh37), every overlap in this tutorial would be wrong while looking perfectly normal. The file name and the download page both say GRCh38, but the strongest check uses the data itself: if the coordinates are hg38, the hg38 sequence at each position must spell the REF allele.
#-----------------------------------------------
# STEP 6: Check the genome build against the hg38 reference sequence
#-----------------------------------------------
# Short working names for the columns we use; chromosome 23 is chrX
gwas <- data.frame(
chr = paste0("chr", ifelse(raw[["#CHR"]] == 23, "X", raw[["#CHR"]])),
pos = raw$POS,
ref = raw$REF,
alt = raw$ALT, # effect allele: beta is for ALT
rsid = raw$rsid,
n_studies = raw$all_meta_N,
beta = raw$all_inv_var_meta_beta, # log odds ratio
se = raw$all_inv_var_meta_sebeta,
p = raw$all_inv_var_meta_p,
het_p = raw$all_inv_var_het_p,
lmso_p = raw$lmso_inv_var_pval, # p-value without the most significant study
af = raw$all_meta_AF
)
gwas$id <- paste(gwas$chr, gwas$pos, gwas$ref, gwas$alt, sep = ":")
gwas$type <- ifelse(nchar(gwas$ref) == 1 & nchar(gwas$alt) == 1, "SNV", "indel")
# Every variant as a genomic range covering its REF allele
var_gr <- GRanges(gwas$chr, IRanges(start = gwas$pos, width = nchar(gwas$ref)))
# Does the hg38 sequence at POS match the REF allele?
gwas$ref_match <- as.character(getSeq(BSgenome.Hsapiens.UCSC.hg38, var_gr)) == gwas$ref
table(type = gwas$type, ref_match = gwas$ref_match)
gwas[which(gwas$rsid == "rs17713054"), c("id", "rsid", "beta", "p")]
Output:
ref_match
type TRUE
indel 437
SNV 7869
id rsid beta p
596 chr3:45818159:G:A rs17713054 0.72868 6.1225e-220
Reading this output.
- Every
REFallele matches hg38, so the coordinates are GRCh38, the same build as our peaks.- The strongest association is rs17713054, the causal variant Downes et al. (2021) identified at 3p21.31. Each copy of its risk allele roughly doubles the odds of critical illness.
Step 7: Keeping Variants That Pass the Build Check
#-----------------------------------------------
# STEP 7: Keep variants whose REF allele matches hg38
#-----------------------------------------------
gwas <- gwas[gwas$ref_match, ]
var_gr <- GRanges(gwas$chr, IRanges(start = gwas$pos, width = nchar(gwas$ref)))
nrow(gwas)
sum(gwas$p < 5e-8)
Output:
[1] 8306
[1] 4237
Reading this output. Nothing is removed from this file, but keep the step for others. About half the variants are genome-wide significant, but they are not independent findings: LD makes neighboring variants significant together, which the next step deals with.
🧬 From Variants to Loci to Credible Sets
Step 8: Defining Loci Around the Lead Variants
A locus groups the variants that most likely share one underlying signal. Without an LD reference we define loci by distance, with greedy clumping: take the most significant variant as the first lead; every genome-wide significant variant within 500 kb of it belongs to it; the most significant remaining variant becomes the next lead, and so on. Every variant in the file, significant or not, then joins the strongest lead within 500 kb, if there is one.
We add one quality check that the HGI file makes possible. The lmso columns repeat the meta-analysis without its most significant study. A lead whose association disappears without one study rests on that study alone. The HGI used a leave-most-significant-biobank-out analysis to flag likely false positives in release 7 (COVID-19 Host Genetics Initiative, 2023); we keep a locus only if its lead still has lmso_p < 0.05. That threshold is our choice, and a lenient one: it asks only that the signal survive the removal of one study.
#-----------------------------------------------
# STEP 8: Define loci around genome-wide significant lead variants
#-----------------------------------------------
window <- 5e5 # 500 kb on each side of a lead variant
# Greedy clumping: walk down the genome-wide significant variants from the smallest
# p-value; a variant becomes a new lead unless it lies within 500 kb of an earlier lead
gw_order <- order(gwas$p)
gw_order <- gw_order[gwas$p[gw_order] < 5e-8]
lead_idx <- Reduce(function(leads, i) {
near <- gwas$chr[leads] == gwas$chr[i] & abs(gwas$pos[leads] - gwas$pos[i]) <= window
if (any(near)) leads else c(leads, i)
}, gw_order, integer(0))
leads <- gwas[lead_idx, ]
leads$locus <- sprintf("L%02d", seq_len(nrow(leads))) # L01 = strongest lead
leads$label <- paste0(leads$chr, ":", leads$pos, " (", ifelse(is.na(leads$rsid), "no rsid", leads$rsid), ")")
# Assign every variant to the strongest lead within 500 kb -- the same rule the
# clumping used, so a variant always joins the lead that "claimed" it
lead_win <- resize(var_gr[lead_idx], width = 2 * window + 1, fix = "center")
hits <- findOverlaps(var_gr, lead_win)
strongest <- tapply(subjectHits(hits), queryHits(hits), min) # leads are ordered strongest first
gwas$locus <- NA_character_
gwas$locus[as.integer(names(strongest))] <- leads$locus[strongest]
# Locus QC: does the lead stay significant without its most significant study?
leads$n_variants <- as.vector(table(factor(gwas$locus, levels = leads$locus)))
leads$lmso_pass <- leads$lmso_p < 0.05
nrow(leads)
table(lmso_pass = leads$lmso_pass)
leads[!leads$lmso_pass, c("locus", "label", "type", "n_studies", "p", "lmso_p")]
# Keep loci whose lead passes; variants outside them are not fine-mapped
keep_loci <- leads$locus[leads$lmso_pass]
fm <- gwas[!is.na(gwas$locus) & gwas$locus %in% keep_loci, ]
leads <- leads[leads$lmso_pass, ]
c(loci = length(keep_loci), variants = nrow(fm), not_in_a_locus = sum(is.na(gwas$locus)))
Output:
[1] 48
lmso_pass
FALSE TRUE
5 43
locus label type n_studies p lmso_p
428 L22 chr2:223837219 (no rsid) indel 2 1.7735e-12 0.5905000
4126 L28 chr12:26114393 (rs137859282) indel 4 2.1102e-11 0.9943757
329 L41 chr2:57089368 (rs12614007) SNV 33 1.2333e-08 0.5283491
2381 L42 chr6:14714491 (rs57696775) indel 4 1.4740e-08 0.1896000
2300 L45 chr5:90862356 (no rsid) indel 3 2.6475e-08 0.6074215
loci variants not_in_a_locus
43 7653 646
Reading this output.
- 48 lead variants become 43 loci after QC. The HGI reported 30 critical-illness loci; the difference comes from their inclusion of 23andMe, their own locus definitions, and our 500 kb windows splitting multi-signal regions such as 3p21.31.
- The five loci removed each rest on one study: without it, their association disappears.
Why 500 kb, and not 1 Mb? With 1 Mb windows, the two strongest 3p21.31 signals fall into one locus, and the single-causal-variant model of Step 9 would give the weaker one almost no probability. Neither window is “correct”: state the one you use, and assign variants to leads with the same rule the clumping used.
Step 9: Fine-Mapping Each Locus With Approximate Bayes Factors
#-----------------------------------------------
# STEP 9: Fine-map each locus with Wakefield approximate Bayes factors
#-----------------------------------------------
# Wakefield (2009) approximate Bayes factor, assuming one causal variant per locus.
# prior_sd = 0.2 on the log-odds scale is coloc's default for case-control traits.
abf_pip <- function(beta, se, prior_sd = 0.2) {
z <- beta / se
V <- se^2
r <- prior_sd^2 / (prior_sd^2 + V)
labf <- 0.5 * (log(1 - r) + r * z^2)
exp(labf - max(labf)) / sum(exp(labf - max(labf))) # subtract the max to avoid overflow
}
fm$pip <- ave(seq_len(nrow(fm)), fm$locus,
FUN = function(i) abf_pip(fm$beta[i], fm$se[i]))
# 95% credible set: the fewest variants whose PIPs add up to at least 0.95
fm <- fm[order(fm$locus, -fm$pip), ]
fm$cum_pip <- ave(fm$pip, fm$locus, FUN = cumsum)
fm$in_cs <- (fm$cum_pip - fm$pip) < 0.95
cs_summary <- do.call(rbind, lapply(split(fm, fm$locus), function(d) data.frame(
locus = d$locus[1],
lead = leads$label[leads$locus == d$locus[1]],
n_variants = nrow(d),
cs_size = sum(d$in_cs),
top_id = d$id[1],
top_rsid = d$rsid[1],
top_pip = round(d$pip[1], 3)
)))
summary(cs_summary$cs_size)
table(cut(cs_summary$cs_size, c(0, 1, 5, 20, 50, Inf), labels = c("1", "2-5", "6-20", "21-50", ">50")))
head(cs_summary[order(-cs_summary$top_pip), ], 10)
cs_summary[cs_summary$locus == "L01", ]
saveRDS(fm, file.path(obj_dir, "finemapped_variants.rds"))
A few lines deserve a closer look:
exp(labf - max(labf)). A log Bayes factor of several hundred, as at 3p21.31, overflowsexp()toInf. Subtracting the largest value before exponentiating leaves the ratios, and therefore the PIPs, unchanged.ave(seq_len(nrow(fm)), fm$locus, FUN = ...)runsabf_pip()once per locus on the row numbers of that locus, so PIPs sum to 1 within each locus.prior_sd = 0.2is how large we expect a true effect to be, as a standard deviation on the log odds scale. It is the default that thecolocpackage uses for case-control traits (Giambartolomei et al., 2014); a quantitative trait needs a different value (see Running This Workflow on Your Own GWAS Results).(fm$cum_pip - fm$pip) < 0.95keeps a variant if the variants ranked above it have not yet reached 0.95. It includes the variant that crosses the threshold, which is the usual definition of a 95% credible set.
Output:
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.00 2.00 6.00 38.33 19.00 1124.00
1 2-5 6-20 21-50 >50
9 12 13 7 2
locus lead n_variants cs_size top_id
L02 L02 chr3:46347209 (rs146375688) 358 1 chr3:46347209:C:T
L05 L05 chr19:10352442 (rs34536443) 118 1 chr19:10352442:G:C
L19 L19 chr11:1219991 (rs35705950) 5 1 chr11:1219991:G:T
L31 L31 chr3:46871689 (rs79152783) 1 1 chr3:46871689:A:G
L03 L03 chr19:4717660 (rs12610495) 97 1 chr19:4717660:A:G
L06 L06 chr6:41515007 (rs2496644) 32 1 chr6:41515007:A:C
L09 L09 chr17:49863260 (rs3848456) 25 1 chr17:49863260:C:A
L13 L13 chr2:60480453 (rs1123573) 29 1 chr2:60480453:A:G
L36 L36 chr4:105897896 (rs34712979) 163 1 chr4:105897896:G:A
L11 L11 chr3:45218795 (rs17278585) 3 2 chr3:45218795:C:T
top_rsid top_pip
L02 rs146375688 1.000
L05 rs34536443 1.000
L19 rs35705950 1.000
L31 rs79152783 1.000
L03 rs12610495 0.998
L06 rs2496644 0.995
L09 rs3848456 0.994
L13 rs1123573 0.985
L36 rs34712979 0.970
L11 rs17278585 0.948
locus lead n_variants cs_size top_id
L01 L01 chr3:45818159 (rs17713054) 1121 7 chr3:45818159:G:A
top_rsid top_pip
L01 rs17713054 0.542
Reading this output.
- Resolution varies widely. Nine loci resolve to a single variant; a few keep dozens or hundreds of variants, where LD leaves the causal variant unresolved.
- One chromosome 17 locus holds two-thirds of all credible-set variants, each with a tiny PIP. Its total probability is the same as any other locus, which is why later steps weight variants by PIP instead of counting them.
- Read each PIP next to
n_variants. A PIP of 1 among hundreds of candidates is a resolved signal; a PIP of 1 with almost no competitors is not.- At 3p21.31, rs17713054 leads a 7-variant credible set. Statistics alone cannot fully separate it from its neighbors.
🔬 Do Risk Variants Fall in Cell-Type Peaks?
We now bring in the chromatin. Part 8 called peaks separately in ten cell types and merged them into 155,853 unified peaks, recording in called_in which cell types called each one. That gives every cell type its own peak set, and the peaks that only one cell type called — its private peaks, 68,021 of the 155,853 — describe what is specific to it.
Cell type names. As in Part 8, peak sets use the cleaned names Signac gave the MACS3 runs:
CD14_monocyte,NK_T_cell_unresolved, and so on. Steps 13 onwards use the object’s own labels (CD14 monocyte,NK/T cell unresolved).
Step 10: Overlapping Credible-Set Variants With the Part 8 Peaks
#-----------------------------------------------
# STEP 10: Overlap credible-set variants with the Part 8 cell-type peaks
#-----------------------------------------------
unified <- readRDS(file.path(recall_dir, "objects", "unified_macs3_peaks.rds"))
# One peak set per cell type: every unified peak that cell type called (Part 8, Step 6)
cell_types <- sort(unique(unlist(unified$called_in)))
peaks_by_ct <- setNames(lapply(cell_types, function(ct) unified[any(unified$called_in == ct)]),
cell_types)
private_by_ct <- lapply(peaks_by_ct, function(pk) pk[pk$n_cell_types == 1])
# Credible-set variants as genomic ranges, carrying their PIP and locus
cs <- fm[fm$in_cs, ]
cs_gr <- GRanges(cs$chr, IRanges(start = cs$pos, width = nchar(cs$ref)),
pip = cs$pip, locus = cs$locus, id = cs$id, rsid = cs$rsid)
overlap_tab <- data.frame(
cell_type = cell_types,
n_peaks = vapply(peaks_by_ct, length, integer(1)),
genome_Mb = vapply(peaks_by_ct, function(pk) sum(width(pk)) / 1e6, numeric(1)),
n_variants = vapply(peaks_by_ct, function(pk) sum(overlapsAny(cs_gr, pk)), integer(1)),
n_loci = vapply(peaks_by_ct, function(pk) length(unique(cs_gr$locus[overlapsAny(cs_gr, pk)])), integer(1)),
pip_in_peaks = vapply(peaks_by_ct, function(pk) sum(cs_gr$pip[overlapsAny(cs_gr, pk)]), numeric(1)),
pip_private = vapply(private_by_ct, function(pk) sum(cs_gr$pip[overlapsAny(cs_gr, pk)]), numeric(1)),
row.names = NULL
)
overlap_tab$genome_Mb <- round(overlap_tab$genome_Mb, 1)
overlap_tab$pip_in_peaks <- round(overlap_tab$pip_in_peaks, 2)
overlap_tab$pip_private <- round(overlap_tab$pip_private, 2)
c(cs_variants = length(cs_gr), in_any_peak = sum(overlapsAny(cs_gr, unified)),
pip_total = round(sum(cs_gr$pip), 2), pip_in_any_peak = round(sum(cs_gr$pip[overlapsAny(cs_gr, unified)]), 2))
overlap_tab[order(-overlap_tab$pip_in_peaks), ]
# Is the 3p21.31 lead, rs17713054, inside any peak?
cs_gr[cs_gr$rsid %in% "rs17713054"]
overlapsAny(cs_gr[cs_gr$rsid %in% "rs17713054"], unified)
Output:
cs_variants in_any_peak pip_total pip_in_any_peak
1648.00 69.00 41.79 3.04
cell_type n_peaks genome_Mb n_variants n_loci pip_in_peaks
2 CD14_monocyte 98141 54.2 54 17 2.33
6 CD8_T_effector 63767 41.0 38 14 2.31
1 B_cell 66410 41.9 40 15 2.06
4 CD4_T_memory 56361 38.0 38 14 1.43
9 NK_cell 52397 36.3 32 13 1.43
8 Dendritic_cell 30332 21.2 19 9 1.15
7 CD8_T_naive 52217 35.9 37 13 0.89
5 CD4_T_naive 56733 37.6 40 13 0.88
10 NK_T_cell_unresolved 46503 34.5 31 11 0.79
3 CD16_monocyte 24525 21.7 16 6 0.28
pip_private
2 0.53
6 0.00
1 0.00
4 0.00
9 0.00
8 0.01
7 0.00
5 0.00
10 0.00
3 0.00
GRanges object with 1 range and 4 metadata columns:
seqnames ranges strand | pip locus id
<Rle> <IRanges> <Rle> | <numeric> <character> <character>
[1] chr3 45818159 * | 0.542234 L01 chr3:45818159:G:A
rsid
<character>
[1] rs17713054
-------
seqinfo: 18 sequences from an unspecified genome; no seqlengths
[1] TRUE
The overlap table has one row per cell type:
| Column | Meaning |
|---|---|
n_peaks | Unified peaks this cell type called |
genome_Mb | Total width of those peaks, in megabases |
n_variants | Credible-set variants inside them |
n_loci | Loci with at least one credible-set variant inside them |
pip_in_peaks | Summed PIP of those variants: how much of the fine-mapped probability this cell type’s open chromatin holds |
pip_private | The same, counting only the cell type’s private peaks |
Reading this output.
- Most of the fine-mapped probability lies outside PBMC open chromatin: only about 7 percent falls in any peak. Where the rest acts — another tissue, a cell state absent from blood, or no regulatory element — this analysis cannot say.
- Cell types with larger peak sets capture more probability by size alone, so this ranking is not yet an enrichment. Step 11 corrects for size.
- Probability in cell-type-private peaks belongs almost entirely to CD14 monocytes.
- rs17713054 is inside a peak, contrary to the expectation in the introduction. Steps 12 and 22 show where.
Step 11: Permutation Tests — Is the Overlap More Than Chance?
The question is simple: if the credible-set variants were somewhere else in the genome, how much probability would fall in each cell type’s peaks? regioneR’s permTest() answers it by moving the variants many times, recomputing the statistic each time, and comparing the observed value with that null distribution. It returns a p-value (the fraction of permutations at least as extreme as the observed value) and a z-score (how many null standard deviations the observed value sits above the null mean).
How the variants are moved matters as much as the test itself:
- Circular randomization (
circularRandomizeRegions()) rotates all variants on a chromosome by the same random distance, wrapping around the end. The variants keep their spacing, so a credible set of 20 variants within 40 kb stays a cluster of 20 variants within 40 kb. - Independent placement moves each variant on its own. A credible set becomes 20 scattered variants, and the null forgets that the real variants come in clusters.
When both the variants and the peaks cluster — and both do, since variants come in LD blocks and peaks in regulatory hubs — whole clusters hit or miss peaks together. That can make the real null more variable than independent placement implies, in which case the independent null gives larger z-scores than the data deserve. We run both and compare.
#-----------------------------------------------
# STEP 11: Permutation test -- is the overlap more than chance?
#-----------------------------------------------
# The genome to randomize over: autosomes and chrX, the chromosomes in the GWAS file
std_chr <- paste0("chr", c(1:22, "X"))
genome_gr <- GRanges(std_chr, IRanges(1, seqlengths(BSgenome.Hsapiens.UCSC.hg38)[std_chr]))
# One evaluation function per cell type: PIP-weighted overlap with that cell type's peaks
make_evaluators <- function(peak_list) {
lapply(peak_list, function(pk) {
force(pk)
function(A, ...) sum(A$pip[overlapsAny(A, pk)])
})
}
# Null 2: every variant spun independently along its own chromosome.
# regioneR's randomizeRegions() drops metadata columns (the PIPs), so we write our own.
independentSpin <- function(A, genome, ...) {
chr_end <- end(genome)[match(as.character(seqnames(A)), as.character(seqnames(genome)))]
new_start <- (start(A) + floor(runif(length(A)) * chr_end)) %% chr_end + 1
ranges(A) <- IRanges(new_start, width = width(A))
A
}
run_perm <- function(A, peak_list, randomize, ntimes = 1000) {
set.seed(1234)
pt <- permTest(A = A, ntimes = ntimes,
randomize.function = randomize,
evaluate.function = make_evaluators(peak_list),
genome = genome_gr,
alternative = "greater", verbose = FALSE)
data.frame(
cell_type = names(pt),
observed = vapply(pt, function(x) x$observed, numeric(1)),
perm_mean = vapply(pt, function(x) mean(x$permuted), numeric(1)),
z = vapply(pt, function(x) x$zscore, numeric(1)),
p = vapply(pt, function(x) x$pval, numeric(1)),
row.names = NULL
)
}
# Null 1: circular randomization rotates each chromosome's variants together,
# so their spacing (and the LD clustering behind it) is preserved
perm_circ <- run_perm(cs_gr, peaks_by_ct, circularRandomizeRegions)
# Null 2: the same variants, each moved on its own
perm_unif <- run_perm(cs_gr, peaks_by_ct, independentSpin)
perm_circ$p_adj <- p.adjust(perm_circ$p, method = "BH")
perm_tab <- merge(perm_circ, perm_unif[, c("cell_type", "z")], by = "cell_type",
suffixes = c("_circular", "_uniform"))
perm_tab$fold <- perm_tab$observed / perm_tab$perm_mean
perm_tab[order(-perm_tab$z_circular),
c("cell_type", "observed", "perm_mean", "fold", "z_circular", "z_uniform", "p", "p_adj")]
Output:
cell_type observed perm_mean fold z_circular z_uniform
6 CD8_T_effector 2.3106770 0.7147326 3.2329253 2.9370 2.7836
1 B_cell 2.0550584 0.7398860 2.7775340 2.3147 2.2152
2 CD14_monocyte 2.3338217 0.9061136 2.5756392 2.2272 2.2574
8 Dendritic_cell 1.1509758 0.3939806 2.9214020 1.7723 1.8063
9 NK_cell 1.4261931 0.6439172 2.2148705 1.5179 1.3726
4 CD4_T_memory 1.4263476 0.6712369 2.1249542 1.3967 1.3191
7 CD8_T_naive 0.8945864 0.6443260 1.3884065 0.4772 0.4128
5 CD4_T_naive 0.8773882 0.6746530 1.3005029 0.3743 0.2922
10 NK_T_cell_unresolved 0.7880967 0.6219535 1.2671313 0.3264 0.2430
3 CD16_monocyte 0.2772593 0.4136878 0.6702139 -0.3198 -0.3580
p p_adj
6 0.01198801 0.1132201
1 0.02897103 0.1132201
2 0.03396603 0.1132201
8 0.07692308 0.1618382
9 0.08091908 0.1618382
4 0.10289710 0.1714952
7 0.27172827 0.3330003
5 0.29370629 0.3330003
10 0.29970030 0.3330003
3 0.47252747 0.4725275
A few details of the code:
evaluate.functionis a list.permTest()accepts several evaluation functions and applies all of them to the same random placement. All ten cell types are therefore tested against identical permutations, which makes their z-scores directly comparable.force(pk)makes each function keep its own peak set.independentSpin()is our own randomizer. regioneR’srandomizeRegions()places regions independently, but returns them without their metadata columns, so thepipcolumn — and with it every evaluation — would be lost. Our function applies a separate random rotation to each variant, within its own chromosome. The two nulls therefore differ in exactly one respect: variants move together, or each on its own.- No mask.
circularRandomizeRegions()uses a mask only to reject rotations that put more thanmax.mask.overlapof the regions in masked sequence, and we set no such limit. Rotated variants can therefore land in assembly gaps, where no peak can exist. That lowers the null mean slightly and favors enrichment slightly, but equally for every cell type, so the comparison between cell types is unaffected.
Reading this output.
- After correcting for size, CD8 effector T cells, B cells and CD14 monocytes show the most enrichment (about 2.5- to 3-fold over chance), but none survives correction for ten tests.
- The two null models give nearly identical z-scores, so keeping the variants’ spacing did not change the conclusion here. With another GWAS or unweighted counts they can disagree, so it is worth running both.
- The cell types are not independent tests: a variant in a widely shared peak counts for every cell type that called it.
Step 12: Specificity and Robustness
Step 11 asks whether a cell type’s peaks hold more risk than random DNA. Three follow-up checks ask whether that result is specific to the cell type, whether it rests on one locus, and whether it depends on the exact variant positions.
#-----------------------------------------------
# STEP 12: Specificity and robustness
#-----------------------------------------------
# (a) Cell-type-private peaks: peaks only that one cell type called
perm_priv <- run_perm(cs_gr, private_by_ct, circularRandomizeRegions)
perm_priv$n_private <- vapply(private_by_ct, length, integer(1))[perm_priv$cell_type]
perm_priv$p_adj <- p.adjust(perm_priv$p, method = "BH")
perm_priv[order(-perm_priv$z), c("cell_type", "n_private", "observed", "perm_mean", "z", "p", "p_adj")]
# (b) Which loci carry each cell type's overlap? PIP in peaks, per locus and cell type
locus_contrib <- sapply(peaks_by_ct, function(pk) {
hit <- overlapsAny(cs_gr, pk)
tapply(cs_gr$pip * hit, factor(cs_gr$locus, levels = leads$locus), sum)
})
locus_contrib <- locus_contrib[rowSums(locus_contrib) > 0, , drop = FALSE]
round(locus_contrib, 2)
# (c) Does the signal sit exactly on the variants? Shift them and re-evaluate
top_ct <- perm_tab$cell_type[which.max(perm_tab$z_circular)]
set.seed(1234)
pt_top <- permTest(A = cs_gr, ntimes = 1000, randomize.function = circularRandomizeRegions,
evaluate.function = make_evaluators(peaks_by_ct[top_ct]),
genome = genome_gr, alternative = "greater", verbose = FALSE)
lz_top <- localZScore(A = cs_gr, pt = pt_top, window = 5000, step = 100)
Output:
cell_type n_private observed perm_mean z p
2 CD14_monocyte 37276 0.5288973107 0.184194508 1.2166 0.0989011
8 Dendritic_cell 4992 0.0110269230 0.016869679 -0.0622 0.1128871
10 NK_T_cell_unresolved 544 0.0000000000 0.002370076 -0.0663 1.0000000
4 CD4_T_memory 1614 0.0003906244 0.005845323 -0.0960 0.1918082
3 CD16_monocyte 735 0.0000000000 0.002414597 -0.1045 1.0000000
7 CD8_T_naive 1618 0.0000000000 0.005306159 -0.1070 1.0000000
5 CD4_T_naive 2505 0.0005986108 0.011561625 -0.1320 0.2257742
9 NK_cell 2363 0.0000000000 0.008700532 -0.1468 1.0000000
6 CD8_T_effector 4803 0.0000000000 0.017956828 -0.2178 1.0000000
1 B_cell 11571 0.0028574002 0.064436098 -0.3640 0.5664336
p_adj
2 0.5644356
8 0.5644356
10 1.0000000
4 0.5644356
3 1.0000000
7 1.0000000
5 0.5644356
9 1.0000000
6 1.0000000
1 1.0000000
B_cell CD14_monocyte CD16_monocyte CD4_T_memory CD4_T_naive CD8_T_effector
L01 0.54 0.00 0.00 0.54 0.00 0.54
L04 0.00 0.15 0.00 0.15 0.15 0.15
L10 0.00 0.29 0.00 0.00 0.00 0.00
L12 0.00 0.00 0.00 0.13 0.00 0.13
L16 0.01 0.01 0.00 0.01 0.01 0.01
L17 0.00 0.04 0.00 0.04 0.04 0.00
L18 0.88 0.88 0.00 0.00 0.00 0.88
L20 0.02 0.03 0.01 0.02 0.02 0.02
L21 0.07 0.07 0.00 0.00 0.07 0.00
L23 0.14 0.14 0.00 0.14 0.14 0.14
L25 0.10 0.10 0.00 0.10 0.10 0.10
L26 0.00 0.00 0.00 0.00 0.00 0.00
L27 0.04 0.04 0.04 0.04 0.04 0.04
L32 0.01 0.07 0.01 0.01 0.01 0.01
L33 0.09 0.26 0.11 0.09 0.09 0.09
L35 0.01 0.06 0.00 0.02 0.08 0.06
L38 0.07 0.09 0.09 0.07 0.07 0.07
L39 0.00 0.04 0.00 0.00 0.00 0.00
L40 0.00 0.00 0.00 0.00 0.00 0.00
L47 0.06 0.06 0.00 0.06 0.06 0.06
CD8_T_naive Dendritic_cell NK_cell NK_T_cell_unresolved
L01 0.00 0.00 0.54 0.00
L04 0.15 0.00 0.15 0.15
L10 0.00 0.00 0.00 0.00
L12 0.13 0.00 0.13 0.13
L16 0.01 0.00 0.01 0.01
L17 0.04 0.00 0.00 0.00
L18 0.00 0.88 0.00 0.00
L20 0.02 0.01 0.02 0.02
L21 0.00 0.00 0.00 0.00
L23 0.14 0.00 0.14 0.14
L25 0.10 0.10 0.10 0.00
L26 0.00 0.00 0.00 0.00
L27 0.04 0.04 0.04 0.04
L32 0.01 0.01 0.01 0.01
L33 0.09 0.00 0.09 0.09
L35 0.03 0.00 0.06 0.06
L38 0.07 0.09 0.07 0.07
L39 0.00 0.00 0.00 0.00
L40 0.00 0.01 0.00 0.00
L47 0.06 0.00 0.06 0.06
Reading this output.
- (a) Private peaks: once shared peaks are removed, almost no enrichment is left. Only CD14 monocytes keep some, and it is not significant. Most risk variants in open chromatin sit in peaks that several cell types share.
- (b) Per-locus table: only 20 of the 43 loci contribute, and two variants carry most of the lymphocyte enrichment of Step 11: rs17713054 (L01), in a peak called by lymphocytes and NK cells, and rs11614702 (L18), in a peak called by B cells, CD14 monocytes, CD8 effector T cells and dendritic cells.
- (c) Local z-score: plotted in Step 19.
🧫 Single-Cell Resolution With gchromVAR
Steps 10 to 12 treat a peak as open or closed in a cell type: a variant counts for a cell type if that cell type called the peak. gchromVAR uses the counts: how many of each sample’s fragments fall in the risk-weighted peaks, compared with background peaks of the same GC content and overall accessibility. It gives one score per sample, so it can score donor pseudobulks — and test them as replicates — or single cells.
Step 13: gchromVAR on Donor-by-Cell-Type Pseudobulks
We sum each cell type’s cells within each donor, as Part 7 did for CD14 monocytes, but now for every label in the object. Each peak’s weight is the summed PIP of all fine-mapped variants inside it — not only credible-set variants, since gchromVAR’s weighting already makes low-PIP variants count little.
#-----------------------------------------------
# STEP 13: gchromVAR on donor-by-cell-type pseudobulks
#-----------------------------------------------
# The Part 8 object: the ten cell types, with the "macs3" assay of unified peaks
atac <- readRDS(file.path(recall_dir, "objects", "atac_macs3_peaks.rds"))
DefaultAssay(atac) <- "macs3"
counts <- LayerData(atac[["macs3"]], layer = "counts")
peak_gr <- granges(atac[["macs3"]])
names(peak_gr) <- rownames(atac[["macs3"]]) # "chr1-10000-10400", as in the count matrix
# Sum the cells of each cell type in each donor -> peaks x (cell types x donors)
group <- paste(atac$cell_type, atac$sample_id, sep = " | ")
pb_counts <- sapply(split(colnames(counts), group),
function(cells) Matrix::rowSums(counts[, cells, drop = FALSE]))
pb_info <- data.frame(
pseudobulk = colnames(pb_counts),
cell_type = sub(" \\| .*", "", colnames(pb_counts)),
sample_id = sub(".* \\| ", "", colnames(pb_counts)),
n_cells = as.vector(table(group)[colnames(pb_counts)]),
row.names = colnames(pb_counts)
)
pb_info$condition <- factor(ifelse(grepl("^severe", pb_info$sample_id), "Severe", "Healthy"),
levels = c("Healthy", "Severe"))
se <- SummarizedExperiment(
assays = list(counts = pb_counts),
rowRanges = peak_gr,
colData = pb_info
)
se <- addGCBias(se, genome = BSgenome.Hsapiens.UCSC.hg38)
# chromVAR needs every peak to have fragments and a GC value (an all-N peak has none)
se <- se[rowSums(assay(se, "counts")) > 0 & !is.na(rowData(se)$bias), ]
# Weight of each peak: summed PIP of all fine-mapped variants inside it
fm_gr <- GRanges(fm$chr, IRanges(start = fm$pos, width = nchar(fm$ref)), pip = fm$pip)
hits <- findOverlaps(rowRanges(se), fm_gr)
peak_w <- numeric(nrow(se))
w_sum <- tapply(fm_gr$pip[subjectHits(hits)], queryHits(hits), sum)
peak_w[as.integer(names(w_sum))] <- w_sum
# The weight matrix needs a column name: without one, gchromVAR 0.3.2 fails
# trying to name the columns of a sparse matrix
w_mat <- Matrix(peak_w, ncol = 1, sparse = TRUE,
dimnames = list(rownames(se), "COVID19_A2_critical"))
w_se <- SummarizedExperiment(
assays = list(weights = w_mat),
rowRanges = rowRanges(se),
colData = DataFrame(trait = "COVID19_A2_critical", row.names = "COVID19_A2_critical")
)
c(peaks = nrow(se), weighted_peaks = sum(peak_w > 0), total_weight = round(sum(peak_w), 2))
set.seed(1234)
bg_peaks <- getBackgroundPeaks(se, niterations = 2000)
gdev <- computeWeightedDeviations(se, w_se, background_peaks = bg_peaks)
pb_info$z <- as.vector(assays(gdev)$z)
saveRDS(gdev, file.path(obj_dir, "gchromvar_pseudobulk.rds"))
# One row per cell type, one column per donor
round(xtabs(z ~ cell_type + sample_id, data = pb_info), 2)
Output:
peaks weighted_peaks total_weight
155853.00 315.00 3.14
sample_id
cell_type healthy1 healthy2 severe1 severe2
B cell -0.48 -0.42 -0.93 -1.08
CD14 monocyte 1.80 2.42 3.47 0.87
CD16 monocyte 0.50 1.48 0.95 1.05
CD4 T memory -1.82 -1.75 -1.56 -1.74
CD4 T naive -1.91 -1.90 -1.42 -2.09
CD8 T effector -0.65 -0.93 -0.34 -0.62
CD8 T naive -1.89 -1.97 -0.94 -0.33
Dendritic cell 0.19 4.23 5.81 2.26
Low quality -0.11 0.48 0.76 -1.95
NK cell -0.70 -0.54 -0.11 0.41
NK/T cell unresolved -1.11 -1.18 -0.68 -1.09
Unresolved 0.55 -0.65 0.00 -0.36
Reading this output.
- Only 315 of the 155,853 peaks carry risk weight, which is why one trait’s scores are noisier than a motif’s in Part 7.
- The table includes
Low qualityandUnresolvedbecause the Part 8 object holds all twelve Part 4 labels. They serve only as reference rows;Unresolvedhas nosevere1cells, so that entry is empty (printed as 0.00).- The myeloid rows are positive and the lymphoid rows negative. This is a different picture from Step 11, and Step 14 explains why.
Step 14: Testing Cell Types With Donors as Replicates
#-----------------------------------------------
# STEP 14: Test cell types with donors as replicates
#-----------------------------------------------
# Is each cell type's mean z, across its four donors, above the average pseudobulk (z = 0)?
ct_test <- do.call(rbind, lapply(split(pb_info, pb_info$cell_type), function(d) {
tt <- t.test(d$z, mu = 0, alternative = "greater")
data.frame(cell_type = d$cell_type[1], n_donors = nrow(d), min_cells = min(d$n_cells),
mean_z = mean(d$z), sd_z = sd(d$z), t = unname(tt$statistic), p = tt$p.value)
}))
ct_test$p_adj <- p.adjust(ct_test$p, method = "BH")
ct_test[order(-ct_test$mean_z), ]
# Exploratory: severe vs healthy within each cell type (two donors per group)
cond_test <- do.call(rbind, lapply(split(pb_info, pb_info$cell_type), function(d) {
tt <- t.test(z ~ condition, data = d, var.equal = TRUE)
data.frame(cell_type = d$cell_type[1],
healthy = mean(d$z[d$condition == "Healthy"]),
severe = mean(d$z[d$condition == "Severe"]),
t = unname(tt$statistic), p = tt$p.value)
}))
cond_test$difference <- cond_test$severe - cond_test$healthy
cond_test$p_adj <- p.adjust(cond_test$p, method = "BH")
cond_test[order(cond_test$p), c("cell_type", "healthy", "severe", "difference", "t", "p", "p_adj")]
Output:
cell_type n_donors min_cells mean_z
Dendritic cell Dendritic cell 4 8 3.1226349
CD14 monocyte CD14 monocyte 4 508 2.1375852
CD16 monocyte CD16 monocyte 4 27 0.9958403
Unresolved Unresolved 3 39 -0.1523052
Low quality Low quality 4 72 -0.2059996
NK cell NK cell 4 99 -0.2352262
CD8 T effector CD8 T effector 4 306 -0.6368786
B cell B cell 4 269 -0.7268802
NK/T cell unresolved NK/T cell unresolved 4 677 -1.0130899
CD8 T naive CD8 T naive 4 27 -1.2830641
CD4 T memory CD4 T memory 4 170 -1.7187528
CD4 T naive CD4 T naive 4 73 -1.8319313
sd_z t p p_adj
Dendritic cell 2.4340101 2.5658356 0.041396222 0.16558489
CD14 monocyte 1.0941906 3.9071533 0.014889348 0.08933609
CD16 monocyte 0.4017213 4.9578665 0.007876704 0.08933609
Unresolved 0.6269585 -0.4207622 0.642584873 0.99995998
Low quality 1.2166512 -0.3386338 0.621399274 0.99995998
NK cell 0.4958715 -0.9487384 0.793626769 0.99995998
CD8 T effector 0.2425602 -5.2513024 0.993275476 0.99995998
B cell 0.3235143 -4.4936510 0.989715845 0.99995998
NK/T cell unresolved 0.2253131 -8.9927299 0.998548680 0.99995998
CD8 T naive 0.7912106 -3.2432936 0.976134234 0.99995998
CD4 T memory 0.1139634 -30.1632437 0.999959979 0.99995998
CD4 T naive 0.2846471 -12.8715962 0.999493958 0.99995998
cell_type healthy severe difference
B cell B cell -0.45261253 -1.0011479 -0.548535372
CD8 T naive CD8 T naive -1.93227463 -0.6338536 1.298421045
NK cell NK cell -0.62008354 0.1496312 0.769714735
CD8 T effector CD8 T effector -0.79312637 -0.4806309 0.312495506
CD4 T memory CD4 T memory -1.78802325 -1.6494823 0.138540966
NK/T cell unresolved NK/T cell unresolved -1.14311014 -0.8830697 0.260040488
Dendritic cell Dendritic cell 2.21124490 4.0340249 1.822780045
Low quality Low quality 0.18446461 -0.5964638 -0.780928446
CD4 T naive CD4 T naive -1.90706140 -1.7568012 0.150260213
Unresolved Unresolved -0.04902828 -0.3588592 -0.309830916
CD14 monocyte CD14 monocyte 2.10681239 2.1683580 0.061545579
CD16 monocyte CD16 monocyte 0.99439293 0.9972877 0.002894796
t p p_adj
B cell 6.779411291 0.02107255 0.2528706
CD8 T naive -4.189049993 0.05253577 0.3152146
NK cell -2.856630971 0.10381017 0.4152407
CD8 T effector -1.573800621 0.25618742 0.6673260
CD4 T memory -1.393467552 0.29813720 0.6673260
NK/T cell unresolved -1.263785759 0.33366302 0.6673260
Dendritic cell -0.678117561 0.56763427 0.9269692
Low quality 0.564257640 0.62941784 0.9269692
CD4 T naive -0.452544005 0.69522692 0.9269692
Unresolved 0.297689357 0.81580324 0.9789639
CD14 monocyte -0.045950197 0.96752544 0.9958396
CD16 monocyte -0.005883709 0.99583963 0.9958396
Reading this output.
- Only the three myeloid cell types are above the average pseudobulk — dendritic cells, CD14 monocytes and CD16 monocytes — but nothing passes 0.05 after adjustment. Dendritic cells have the highest mean but the most variable donors, one built from 8 cells.
- This agrees with Yu et al. (2022), who found COVID-19 severity variants enriched in monocytes and dendritic cells.
- It disagrees with Step 11 because the two methods ask different questions. The permutation test asks whether a cell type called a peak at a variant; gchromVAR asks how many of a sample’s fragments fall in risk-weighted peaks, whoever called them. A weakly open peak, like the one at rs17713054 (Step 22), counts fully in the first and adds little to the second.
- The second table finds no difference between severe and healthy donors in any cell type. The myeloid enrichment is present in both, as expected for an inherited property.
Step 15: gchromVAR in Every Cell
For display, we compute the same score in every single cell. These per-cell z-scores are what the UMAP in Step 20 shows; they are not tested.
#-----------------------------------------------
# STEP 15: gchromVAR in every cell
#-----------------------------------------------
# The same peaks as Step 13, now one column per cell
cell_se <- SummarizedExperiment(
assays = list(counts = counts[match(rownames(se), rownames(counts)), ]),
rowRanges = rowRanges(se)
)
# Background peaks depend only on GC content and total fragments per peak,
# which are the same here as in Step 13, so bg_peaks can be reused
cell_dev <- computeWeightedDeviations(cell_se, w_se, background_peaks = bg_peaks)
atac$covid_z <- as.vector(assays(cell_dev)$z)
summary(atac$covid_z)
sort(tapply(atac$covid_z, atac$cell_type, median), decreasing = TRUE)
saveRDS(atac[[]][, c("cell_type", "sample_id", "condition", "covid_z")],
file.path(obj_dir, "gchromvar_cell_scores.rds"))
Output:
Min. 1st Qu. Median Mean 3rd Qu. Max.
-1.03195 -0.45055 -0.30353 -0.02026 0.06250 15.40179
Dendritic cell CD14 monocyte CD16 monocyte
-0.06859377 -0.08330201 -0.17737879
Low quality Unresolved NK cell
-0.26684499 -0.26779108 -0.28177233
CD8 T effector NK/T cell unresolved CD4 T memory
-0.29296003 -0.30649605 -0.36816635
B cell CD4 T naive CD8 T naive
-0.37724510 -0.38146551 -0.38592680
Reading this output.
- Per-cell scores are noisy and skewed: most cells have few fragments in the weighted peaks, so they sit just below 0.
- The per-cell ranking matches Step 14, with dendritic cells, CD14 monocytes and CD16 monocytes on top.
Why
bg_peakscan be reused. chromVAR picks background peaks from each peak’s GC content and total fragment count, which are the same here as in Step 13.
🧩 Do the Risk Alleles Change Transcription Factor Motifs?
A variant inside an open peak is a candidate; a variant that also creates or destroys a transcription factor motif is a candidate with a mechanism. For each candidate SNV we score all 633 JASPAR2020 motifs on 41 bp of sequence around it, once with the REF allele and once with the ALT allele, and look at what changes.
Step 16: Scoring Both Alleles With motifmatchr
#-----------------------------------------------
# STEP 16: Do the alleles change transcription factor motifs?
#-----------------------------------------------
# The same 633 JASPAR2020 motifs as Parts 4 to 7
pfm <- getMatrixSet(JASPAR2020, opts = list(species = 9606, collection = "CORE", all_versions = FALSE))
motif_names <- setNames(TFBSTools::name(pfm), TFBSTools::ID(pfm))
# SNVs to score: rs17713054 (Downes et al. 2021) plus credible-set SNVs with PIP > 0.1 inside a peak
cs$in_peak <- overlapsAny(cs_gr, unified)
to_score <- cs[cs$type == "SNV" & (cs$rsid %in% "rs17713054" | (cs$in_peak & cs$pip > 0.1)), ]
to_score <- to_score[!duplicated(to_score$id), ]
# 41 bp around each variant, once with the REF and once with the ALT allele
flank <- 20
ref_seq <- getSeq(BSgenome.Hsapiens.UCSC.hg38, GRanges(to_score$chr, IRanges(to_score$pos - flank, to_score$pos + flank)))
alt_seq <- ref_seq
subseq(alt_seq, start = flank + 1, width = 1) <- DNAStringSet(to_score$alt)
names(ref_seq) <- names(alt_seq) <- to_score$id
score_motifs <- function(seqs, p_cutoff) {
m <- matchMotifs(pfm, seqs, genome = BSgenome.Hsapiens.UCSC.hg38, bg = "genome",
out = "scores", p.cutoff = p_cutoff)
list(score = as.matrix(motifScores(m)), match = as.matrix(motifMatches(m)))
}
# Best score of every motif anywhere in each sequence: p.cutoff = 1 lets every
# position count, so a score exists for both alleles even without a match
ref_all <- score_motifs(ref_seq, 1)
alt_all <- score_motifs(alt_seq, 1)
# Matches at two thresholds: 5e-5, the genome-wide default used since Part 4,
# and 5e-4, a lenient threshold for scoring one candidate site
ref_strict <- score_motifs(ref_seq, 5e-5)$match; alt_strict <- score_motifs(alt_seq, 5e-5)$match
ref_lenient <- score_motifs(ref_seq, 5e-4)$match; alt_lenient <- score_motifs(alt_seq, 5e-4)$match
# A motif "changes" if it matches one allele but not the other, at either threshold
ids <- colnames(ref_all$score)
allele_tab <- do.call(rbind, lapply(seq_len(nrow(to_score)), function(i) {
lenient <- ref_lenient[i, ] != alt_lenient[i, ]
strict <- ref_strict[i, ] != alt_strict[i, ]
changed <- which(lenient | strict)
if (length(changed) == 0) return(NULL)
data.frame(id = to_score$id[i], rsid = to_score$rsid[i], pip = round(to_score$pip[i], 3),
motif = motif_names[ids[changed]],
ref_score = round(ref_all$score[i, changed], 2),
alt_score = round(alt_all$score[i, changed], 2),
lenient = lenient[changed], # match status differs at p < 5e-4
strict = strict[changed], # match status differs at p < 5e-5
row.names = NULL)
}))
allele_tab$delta <- allele_tab$alt_score - allele_tab$ref_score
allele_tab$change <- ifelse(allele_tab$delta > 0, "gained", "lost") # by the ALT allele
allele_tab <- allele_tab[order(allele_tab$id, -abs(allele_tab$delta)), ]
nrow(to_score)
allele_tab[allele_tab$rsid %in% "rs17713054", ]
table(strict = allele_tab$strict)
head(allele_tab[allele_tab$strict, ], 15)
Output:
[1] 7
id rsid pip motif ref_score alt_score lenient
3 chr3:45818159:G:A rs17713054 0.542 CEBPG -2.01 7.69 TRUE
4 chr3:45818159:G:A rs17713054 0.542 CEBPG(var.2) 2.30 9.63 TRUE
5 chr3:45818159:G:A rs17713054 0.542 CEBPD 2.84 9.49 TRUE
6 chr3:45818159:G:A rs17713054 0.542 NFIL3 4.19 10.49 TRUE
1 chr3:45818159:G:A rs17713054 0.542 CEBPB 4.03 9.71 TRUE
2 chr3:45818159:G:A rs17713054 0.542 CEBPE 5.23 10.06 TRUE
strict delta change
3 FALSE 9.70 gained
4 FALSE 7.33 gained
5 FALSE 6.65 gained
6 FALSE 6.30 gained
1 FALSE 5.68 gained
2 FALSE 4.83 gained
strict
FALSE TRUE
61 3
id rsid pip motif ref_score alt_score
19 chr13:112886111:C:T rs3742238 0.285 ZNF460 3.70 -3.78
13 chr13:112886111:C:T rs3742238 0.285 PLAG1 9.66 3.84
23 chr21:33949755:C:T rs75912728 0.132 SMAD2::SMAD3::SMAD4 11.22 -0.10
lenient strict delta change
19 TRUE TRUE -7.48 lost
13 FALSE TRUE -5.82 lost
23 TRUE TRUE -11.32 lost
How the scoring works:
- 41 bp per allele. Twenty bases on each side of the variant are enough for any JASPAR motif that overlaps it; a motif elsewhere in the window scores the same in both alleles and cannot appear as a change.
p.cutoff = 1gives a score for both alleles. With a threshold, motifmatchr reports a motif’s best score only where it matches, and 0 elsewhere — so a “change” could simply be a score compared with 0. Withp.cutoff = 1every position counts, andref_score,alt_scoreanddeltaare the real best log-odds scores of the motif in each sequence.- Two thresholds.
strictuses p < 5e-5, the default that Parts 4 to 7 used to scan hundreds of thousands of peaks, where a strict threshold is what keeps false matches manageable.lenientuses p < 5e-4, a reasonable threshold when you score one candidate site rather than the genome — but ten times more permissive, so its calls need support from other evidence. A motif appears in the table if its match status differs between the alleles at either threshold.
Reading this output.
- rs17713054 is the positive control: its risk
Aallele creates C/EBP-family motifs, the “second CEBPB motif” Downes et al. (2021) reported. The gains hold only at the lenient threshold, so the site is modest by genome-wide standards; Downes et al. also supported it with allele-specific accessibility.- Only three changes hold at the strict threshold, all losses by a risk allele: at rs3742238 (in the CD14-monocyte-private peak of Step 12) and rs75912728.
- The full results are in
tables/allelic_motif_changes.csv. Its largest block belongs to rs11614702, the highest-PIP variant in any peak:
📊 Visualization
The same R session continues. Six figures go from the whole GWAS down to single loci.
The GWAS Signal: A Manhattan Plot of the 1e-5 File
#-----------------------------------------------
# STEP 17: Manhattan plot of the 1e-5 file
#-----------------------------------------------
chr_len <- seqlengths(BSgenome.Hsapiens.UCSC.hg38)[std_chr]
offset <- setNames(c(0, cumsum(as.numeric(chr_len))[-length(chr_len)]), std_chr)
man <- gwas
man$x <- man$pos + offset[man$chr]
man$chr <- factor(man$chr, levels = std_chr)
axis_pos <- offset + chr_len / 2
p_man <- ggplot(man, aes(x = x, y = -log10(p), colour = as.integer(chr) %% 2 == 0)) +
geom_point(size = 0.6) +
geom_hline(yintercept = -log10(5e-8), linetype = "dashed", colour = "grey40") +
geom_text_repel(data = leads[1:8, ], aes(x = pos + offset[chr], y = -log10(p), label = rsid),
inherit.aes = FALSE, size = 3, min.segment.length = 0) +
scale_colour_manual(values = c("#2E86AB", "#A23B72"), guide = "none") +
scale_x_continuous(breaks = axis_pos, labels = sub("chr", "", std_chr), expand = c(0.01, 0)) +
scale_y_log10() +
labs(x = "Chromosome", y = "-log10(p), log scale",
title = "COVID-19 HGI release 7, A2 critical illness: variants with p < 1e-5") +
theme_classic(base_size = 11) +
theme(axis.text.x = element_text(size = 7))
ggsave(file.path(plot_dir, "01_manhattan.png"), p_man, width = 11, height = 4.5, dpi = 300)
Reading this figure. The y-axis is log-scaled because the p-values span more than 200 orders of magnitude; the dashed line is genome-wide significance.
- Nothing sits below -log10(p) = 5: that floor is the 1e-5 filter, not biology.
- The 3p21.31 region on chromosome 3 dominates, led by rs17713054.
Fine-Mapping Resolution per Locus
#-----------------------------------------------
# STEP 18: Credible-set size and top PIP per locus
#-----------------------------------------------
cs_plot <- cs_summary
cs_plot$lead <- factor(cs_plot$lead, levels = cs_plot$lead[order(cs_plot$cs_size)])
# Points, not bars: on a log axis a bar of length 1 would be invisible
p_cs <- ggplot(cs_plot, aes(x = cs_size, y = lead)) +
geom_segment(aes(x = 0.8, xend = cs_size, yend = lead), colour = "grey75") +
geom_point(aes(colour = top_pip), size = 2.8) +
scale_x_log10(limits = c(0.8, NA)) +
scale_colour_gradient(low = "#BDBDBD", high = "#A23B72", limits = c(0, 1), name = "Top PIP") +
labs(x = "Variants in the 95% credible set (log scale)", y = NULL,
title = "Fine-mapping resolution per locus") +
theme_classic(base_size = 10)
ggsave(file.path(plot_dir, "02_credible_sets.png"), p_cs, width = 7, height = 8, dpi = 300)
Reading this figure. One row per locus: the dot shows credible-set size (log scale), colored by the top variant’s PIP. Small credible sets go with high top PIPs and large ones with low, the chromosome 17 locus at the top being the extreme.
Enrichment by Cell Type, Under Three Nulls
#-----------------------------------------------
# STEP 19: Enrichment per cell type, circular vs uniform randomization
#-----------------------------------------------
perm_long <- rbind(
data.frame(cell_type = perm_tab$cell_type, z = perm_tab$z_circular, null = "Circular (keeps spacing)"),
data.frame(cell_type = perm_tab$cell_type, z = perm_tab$z_uniform, null = "Uniform (independent)"),
data.frame(cell_type = perm_priv$cell_type, z = perm_priv$z, null = "Circular, private peaks only")
)
perm_long$cell_type <- factor(perm_long$cell_type,
levels = perm_tab$cell_type[order(perm_tab$z_circular)])
perm_long$null <- factor(perm_long$null, levels = c("Circular (keeps spacing)",
"Uniform (independent)",
"Circular, private peaks only"))
p_perm <- ggplot(perm_long, aes(x = z, y = cell_type, colour = null)) +
geom_vline(xintercept = 0, colour = "grey60") +
geom_point(size = 3) +
scale_colour_manual(values = c("#2E86AB", "#F18F01", "#A23B72"), name = "Null model") +
labs(x = "Permutation z-score (PIP-weighted overlap)", y = NULL,
title = "Credible-set variants in cell-type peaks") +
theme_classic(base_size = 11) +
theme(legend.position = "bottom")
ggsave(file.path(plot_dir, "03_enrichment_by_celltype.png"), p_perm, width = 8, height = 5, dpi = 300)
png(file.path(plot_dir, "04_local_zscore.png"), width = 7, height = 4.5, units = "in", res = 300)
plot(lz_top)
dev.off()
Reading this figure. One point per null model for each cell type; the gray line marks chance.
- The circular and independent nulls (blue and orange) nearly coincide.
- With private peaks only (purple), the enrichment disappears except for a little in CD14 monocytes: the lymphocyte enrichment comes from shared peaks.
Reading this figure. The z-score for CD8 effector T cells as all variants are shifted together by up to 5 kb. There is no sharp peak at 0, and some shifts give higher z-scores than the real positions. The enrichment rests on a few loci (Step 12), not on risk variants sitting precisely in CD8 effector T-cell peaks.
gchromVAR Scores on the UMAP
#-----------------------------------------------
# STEP 20: gchromVAR scores on the UMAP and per cell type
#-----------------------------------------------
p_umap <- FeaturePlot(atac, features = "covid_z", reduction = "umap.harmony",
min.cutoff = "q02", max.cutoff = "q98", order = TRUE) +
scale_colour_gradient2(low = "#2E86AB", mid = "grey90", high = "#A23B72", midpoint = 0,
name = "gchromVAR z") +
ggtitle("Per-cell COVID-19 critical-illness score")
ct_order <- names(sort(tapply(atac$covid_z, atac$cell_type, median)))
p_violin <- ggplot(atac[[]], aes(x = factor(cell_type, levels = ct_order), y = covid_z)) +
geom_violin(fill = "grey85", colour = NA, scale = "width") +
geom_boxplot(width = 0.15, outlier.shape = NA) +
geom_hline(yintercept = 0, linetype = "dashed", colour = "grey40") +
coord_flip(ylim = quantile(atac$covid_z, c(0.01, 0.99))) +
labs(x = NULL, y = "gchromVAR z (per cell)") +
theme_classic(base_size = 11)
ggsave(file.path(plot_dir, "05_gchromvar_umap_violin.png"), p_umap | p_violin,
width = 13, height = 5.5, dpi = 300)
Reading this figure. The highest-scoring cells concentrate in the CD14 monocyte clusters (compare Part 4’s annotated UMAP), with scattered high scores elsewhere, and the violins repeat the ranking of Step 15. Per-cell scores are noisy; SCAVENGE’s network propagation was designed to sharpen exactly this kind of picture.
Donor-Level Scores per Cell Type
#-----------------------------------------------
# STEP 21: Donor-level gchromVAR z per cell type
#-----------------------------------------------
pb_info$cell_type_f <- factor(pb_info$cell_type, levels = ct_test$cell_type[order(ct_test$mean_z)])
p_donor <- ggplot(pb_info, aes(x = z, y = cell_type_f)) +
geom_vline(xintercept = 0, colour = "grey60") +
geom_point(aes(colour = condition, shape = sample_id), size = 2.8) +
stat_summary(fun = mean, geom = "point", shape = 3, size = 4, colour = "black") +
scale_colour_manual(values = c(Healthy = "#2E86AB", Severe = "#C73E1D")) +
labs(x = "gchromVAR z (donor pseudobulk)", y = NULL, colour = NULL, shape = NULL,
title = "COVID-19 critical-illness enrichment by cell type and donor") +
theme_classic(base_size = 11)
ggsave(file.path(plot_dir, "06_gchromvar_donors.png"), p_donor, width = 8, height = 5, dpi = 300)
Reading this figure. One point per donor (shape = donor, color = condition); the large black cross is the mean, and the small colored cross is
severe2. The myeloid rows sit right of zero in every donor, and healthy and severe donors mix within every row: no cell type separates by condition.
Looking at the Variants on the Chromatin
Two loci, chosen by rule rather than by eye: the credible-set variant with the highest PIP that sits inside a peak, and rs17713054, the strongest association in the file.
#-----------------------------------------------
# STEP 22: Look at the variants on the chromatin
#-----------------------------------------------
# The credible-set variant with the highest PIP that sits in a peak
top_hit <- cs_gr[overlapsAny(cs_gr, unified)]
top_hit <- top_hit[which.max(top_hit$pip)]
top_hit
DefaultAssay(atac) <- "macs3"
# A 1 bp variant is too thin to see in a 40 kb window: draw each one 201 bp wide
cs_track <- resize(cs_gr, width = 201, fix = "center")
p_top <- CoveragePlot(
object = atac,
region = top_hit,
group.by = "cell_type",
ranges = cs_track,
ranges.title = "Credible set",
region.highlight = resize(top_hit, width = 200, fix = "center"),
extend.upstream = 20000,
extend.downstream = 20000
)
ggsave(file.path(plot_dir, "07_coverage_top_variant.png"), p_top, width = 9, height = 8, dpi = 300)
# rs17713054 at 3p21.31: the strongest association in the file
rs177 <- cs_gr[cs_gr$rsid %in% "rs17713054"]
p_rs177 <- CoveragePlot(
object = atac,
region = rs177,
group.by = "cell_type",
ranges = cs_track,
ranges.title = "Credible set",
region.highlight = resize(rs177, width = 200, fix = "center"),
extend.upstream = 20000,
extend.downstream = 20000
)
ggsave(file.path(plot_dir, "08_coverage_rs17713054.png"), p_rs177, width = 9, height = 8, dpi = 300)
Output:
GRanges object with 1 range and 4 metadata columns:
seqnames ranges strand | pip locus id
<Rle> <IRanges> <Rle> | <numeric> <character> <character>
[1] chr12 132481571 * | 0.881677 L18 chr12:132481571:G:A
rsid
<character>
[1] rs11614702
Reading this output. The highest-PIP credible-set variant inside a peak is rs11614702 (L18), one of the two variants that carry most of the enrichment in Step 12.
Reading this figure. Tracks show accessibility per cell type, genes, Part 8 peaks and credible-set variants (each drawn 201 bp wide so it is visible); the gray band marks rs11614702. Its peak is open mainly in dendritic cells and CD14 monocytes — the peak gchromVAR rewards in Step 14. The nearest gene, FBRSL1, is not necessarily its target; linking an element to a gene needs other data.
Reading this figure. Everything at 3p21.31 is weak (note the y-axis, 0 to 35). A small peak at rs17713054 is open in NK cells and, more faintly, in some lymphocytes, but not in monocytes. Downes et al. (2021) found no immune-cell accessibility here and placed the effect in lung epithelium. One dataset of four donors cannot resolve the difference; check donor support with
CoveragePlot(..., group.by = "sample_id")before reading anything into it.
💾 Saving Results for Reuse and Supplementary Material
#-----------------------------------------------
# STEP 23: Save tables and session information
#-----------------------------------------------
write.csv(leads, file.path(table_dir, "loci_lead_variants.csv"), row.names = FALSE)
write.csv(cs_summary, file.path(table_dir, "credible_set_summary.csv"), row.names = FALSE)
write.csv(fm, file.path(table_dir, "finemapped_variants.csv"), row.names = FALSE)
write.csv(overlap_tab, file.path(table_dir, "credible_set_overlap_by_celltype.csv"), row.names = FALSE)
write.csv(perm_tab, file.path(table_dir, "permutation_enrichment.csv"), row.names = FALSE)
write.csv(perm_priv, file.path(table_dir, "permutation_enrichment_private_peaks.csv"), row.names = FALSE)
write.csv(locus_contrib, file.path(table_dir, "pip_in_peaks_by_locus.csv"))
write.csv(pb_info[, c("pseudobulk", "cell_type", "sample_id", "condition", "n_cells", "z")],
file.path(table_dir, "gchromvar_pseudobulk_z.csv"), row.names = FALSE)
write.csv(ct_test, file.path(table_dir, "gchromvar_celltype_test.csv"), row.names = FALSE)
write.csv(cond_test, file.path(table_dir, "gchromvar_condition_test.csv"), row.names = FALSE)
write.csv(allele_tab, file.path(table_dir, "allelic_motif_changes.csv"), row.names = FALSE)
writeLines(capture.output(sessionInfo()), file.path(gwas_dir, "sessionInfo_part9.txt"))
Together with the objects saved along the way, these files let anyone reproduce every table and figure without rerunning the fine-mapping, the permutations or gchromVAR:
| File | Contents |
|---|---|
objects/finemapped_variants.rds | All 7,653 fine-mapped variants with locus, PIP and credible-set membership |
objects/gchromvar_pseudobulk.rds | gchromVAR deviations for the 47 donor-by-label pseudobulks |
objects/gchromvar_cell_scores.rds | Per-cell gchromVAR z-scores with cell type, donor and condition |
tables/loci_lead_variants.csv | The 43 loci that passed QC, with lead variants and lmso p-values |
tables/credible_set_summary.csv | Credible-set size and top variant per locus |
tables/finemapped_variants.csv | The fine-mapped variants as a table |
tables/credible_set_overlap_by_celltype.csv | The Step 10 overlap table |
tables/permutation_enrichment.csv | Step 11: circular and independent nulls |
tables/permutation_enrichment_private_peaks.csv | Step 12: private peaks |
tables/pip_in_peaks_by_locus.csv | Step 12: PIP in peaks by locus and cell type |
tables/gchromvar_pseudobulk_z.csv | Step 13: one z-score per pseudobulk |
tables/gchromvar_celltype_test.csv | Step 14: cell types tested with donors as replicates |
tables/gchromvar_condition_test.csv | Step 14: the exploratory severe-vs-healthy test |
tables/allelic_motif_changes.csv | Step 16: motif gains and losses per candidate variant |
sessionInfo_part9.txt | Exact package versions |
🔁 Running This Workflow on Your Own GWAS Results
Nothing in Steps 5 to 23 depends on COVID-19 or on the HGI. This section collects what changes when you bring another GWAS.
What the HGI File Is, and What It Is Not
The file used here is ordinary GWAS summary statistics. Before release, the HGI applied an allele frequency filter (0.001) and an imputation quality filter (INFO 0.6) within each study and combined the 34 studies by fixed-effects inverse-variance meta-analysis; the version we downloaded then keeps only variants with p < 1 x 10^-5. It is not fine-mapped, and it carries no annotation or cell-type information. Everything after the download — defining loci, fine-mapping, overlap, permutations, gchromVAR — is done in this tutorial, so it can be done for any GWAS.
The Seven Columns You Need
Steps 6 onwards work on one data frame, gwas. Whatever your source file looks like, the job is to fill these columns:
gwas column | Meaning | Notes |
|---|---|---|
chr | Chromosome, as chr1 … chr22, chrX | Must match the naming of your peaks |
pos | Position of the variant | Same genome build as your peaks (hg38 here) |
ref | The reference allele, as in the genome | Not always the file’s “other allele” — see below |
alt | The alternative allele | beta must be the effect of this allele |
beta | Effect size; log odds ratio for a case-control trait | If the file gives an odds ratio, beta = log(OR) |
se | Standard error of beta | If missing, recover it from the confidence interval |
p | p-value | Used for locus definition and filtering |
rsid and af are optional and only used for labels and tables. n_studies, het_p and lmso_p are specific to a meta-analysis; most files have no leave-one-study-out result, and setting lmso_p to 0 makes every locus pass the Step 8 check.
The effect allele is the one that matters. A GWAS reports the effect of one allele, whatever it calls it. The HGI’s effect allele is always
ALT, soREFmatches the genome and the sign ofbetarefers toALT. Other files report the effect of whichever allele was tested, which may be the reference allele. Step 6’s build check would then reject those variants as mismatches; worse, if you swap the alleles without flipping the sign ofbeta, every direction of effect in Step 16 becomes wrong.
✅ Best Practices for GWAS Enrichment in Single-Cell Chromatin Data
1. Check the genome build against the sequence, not the label. A file name, a README and a download page can all be wrong or misread. Comparing REF alleles with the reference genome (Step 6) takes one line and settles the question. Lifted-over files deserve the same check, because liftover can move a variant without changing its alleles.
2. Join on chromosome, position and alleles, never on rsid. 112 of the 8,306 variants here have no rsid, mostly indels, and rsids are merged and retired between dbSNP builds. Coordinates and alleles are what overlap with peaks anyway.
3. Know which allele the effect refers to. In HGI files beta is for ALT. Check the convention against a published result before interpreting directions: here rs9305744, the TMPRSS2 variant the HGI paper reports as protective (odds ratio 0.92), has beta = -0.086 for its A allele, an odds ratio of 0.92.
4. Use the meta-analysis’s own robustness columns. The leave-most-significant-study-out p-value removed five of 48 loci here (Step 8), four of them indels typed in two to four studies. A genome-wide significant p-value that rests on one study is a hypothesis, not a locus.
5. Define loci and assign variants with the same rule, and state the window. Clumping by one rule and assigning by another (nearest lead, say) can move a strong signal into a neighboring locus. With 1 Mb windows, two 3p21.31 signals 529 kb apart would have shared one locus and one single-causal-variant model.
6. Read every PIP next to the number of variants it competed with. A PIP of 1 out of 357 candidates and a PIP of 1 out of one candidate are different findings (Step 9). Under the single-causal-variant assumption, and for a meta-analysis, treat PIPs as a ranking within a locus rather than calibrated probabilities (Kanai et al., 2022).
7. Weight variants by PIP, not by count. One locus on chromosome 17 holds 1,124 of the 1,648 credible-set variants. Counting variants would let that single region dominate every overlap; PIP weighting gives each locus the same total say.
8. Randomize with a null that keeps the structure of the data. GWAS variants come in LD clusters, and peaks come in regulatory clusters. Circular randomization moves each chromosome’s variants together and keeps their spacing; moving them independently does not (Step 11). Whichever null you use, say which, because it changes the z-scores.
9. Separate “regulatory DNA” from “this cell type”. Risk variants are enriched in open chromatin in general, so almost any cell type’s peaks will look enriched against the whole genome. Cell-type-private peaks (Step 12) and comparisons between cell types are what speak to cell-type specificity.
10. Look at which loci carry an enrichment. A cell-type result driven by one locus is a statement about that locus. The per-locus table in Step 12 shows this directly, and a leave-one-locus-out rerun of Step 11 tests it.
11. Make the donor the unit of replication. Per-cell gchromVAR scores are good for a UMAP; tests belong on donor pseudobulks (Step 14), the same principle that separated 6,445 from 96 differential peaks in Part 5.
12. Be careful with small cell types. Dendritic cells contributed 87 nuclei and the weakest peak set of Part 8; in the pseudobulk, one donor has eight of them. A high or low score from such a group is a reason to look again, not a conclusion.
13. Keep genetic enrichment and disease-state differences apart. Which cell types carry risk variants (Steps 10 to 14) and whether those regions change accessibility in severe disease (the second test in Step 14) are different questions. The first is a property of the genome and the cell type; the second is a property of the patients.
14. Report the strongest locus, whatever it shows. rs17713054 is the strongest association in the file, and published evidence places its mechanism in lung epithelium. Here it turned out to sit in a weak peak called in NK cells and three lymphocyte types (Steps 10, 12 and 22) — the opposite of what the introduction predicted. Report that, with how weak the peak is, rather than leaving the top locus out because it does not fit.
15. Know which labels you are scoring. The Part 8 object holds all twelve Part 4 labels, so Steps 13 to 15 also score Low quality and Unresolved cells. Keep them as reference rows, or subset to the ten cell types before Step 13, but never interpret them as cell types; and remember that every gchromVAR z-score is relative to the average of whatever pseudobulks you include.
16. Run more than one kind of test, and expect them to disagree. A peak-overlap permutation test and gchromVAR ranked the cell types differently here (Steps 11 and 14): one counts whether a cell type called a peak at a variant, the other how many of its fragments fall in risk-weighted peaks. Report both, and explain the difference instead of choosing the one you prefer.
17. Use the right tool when you have more data. With genome-wide summary statistics, use S-LDSC for heritability enrichment (Finucane et al., 2015). With an ancestry-matched LD reference, use a multi-signal fine-mapping method. The lightweight workflow here is for exploring a filtered file, not a substitute for those.
18. Record versions, including where each package came from. gchromVAR was installed from GitHub, so record its commit: packageDescription("gchromVAR")$RemoteSha returns it after a remotes::install_github() install, and the session information saved in Step 23 lists the rest.
⚠️ Common Pitfalls and How to Avoid Them
| Pitfall | Why it happens | How to avoid it |
|---|---|---|
| Overlapping hg19 variants with hg38 peaks | Builds are not in the column names; many public files are GRCh37 | Check REF against the hg38 sequence (Step 6) |
| Zero overlaps on one chromosome | The HGI codes chromosome X as 23, peaks use chrX | Convert chromosome names before building ranges (Step 6) |
| Matching variants by rsid | rsids are missing for indels and change between dbSNP builds | Use chromosome, position, REF and ALT (Step 6) |
Reading all_meta_N as a sample size | The name suggests it | It is the number of studies; sample sizes are in the cases and controls columns |
| Trusting a locus carried by one study | Meta-analysis p-values hide where the signal comes from | Require the lead to survive lmso (Step 8) |
| Assigning variants to the nearest lead after greedy clumping | It sounds natural | Assign to the strongest lead within the window, the clumping rule (Step 8) |
| Running S-LDSC on a p < 1e-5 file | The file looks like summary statistics | S-LDSC needs all variants and an LD panel; use the genome-wide file |
| Reading a PIP of 1 as a solved locus | PIPs sum to 1 even with one candidate | Read n_variants next to every PIP (Step 9) |
| Expecting ABF to find a second signal | One causal variant per locus is built into the method | Split loci with clear separate peaks; use LD-based methods when available |
| Counting credible-set variants | Counting is the obvious overlap statistic | Weight by PIP (Steps 10 to 12) |
| Uniform randomization of clustered variants | overlapPermTest(), regioneR’s convenience function, uses randomizeRegions() | Use circularRandomizeRegions() and compare (Step 11) |
| Losing the PIP weights in permutations | randomizeRegions() returns ranges without metadata columns | Use a randomizer that keeps mcols() (Step 11) |
| Calling every enriched cell type “causal” | All regulatory DNA is enriched against the genome | Test private peaks and compare cell types (Step 12) |
| Testing per-cell scores | 30,000 cells look like 30,000 replicates | Test donor pseudobulks (Step 14) |
| Reading a severe-vs-healthy difference as genetic | Both questions use the same scores | Report the condition test separately, as exploratory (Step 14) |
| Scoring indels with a one-base substitution | The motif step replaces one base | Score SNVs only, or build each allele’s sequence explicitly |
| Reporting motif changes without motif matches | delta_score exists for every motif | Require a match in at least one allele (Step 16) |
📚 References
- COVID-19 Host Genetics Initiative. A second update on mapping the human genetic architecture of COVID-19. Nature. 2023;621(7977):E7-E26. doi:10.1038/s41586-023-06355-3
- COVID-19 Host Genetics Initiative. Results: data release 7. https://www.covid19hg.org/results/r7/
- Downes DJ, Cross AR, Hua P, et al. Identification of LZTFL1 as a candidate effector gene at a COVID-19 risk locus. Nature Genetics. 2021;53(11):1606-1615. doi:10.1038/s41588-021-00955-3
- Schaid DJ, Chen W, Larson NB. From genome-wide associations to candidate causal variants by statistical fine-mapping. Nature Reviews Genetics. 2018;19(8):491-504. doi:10.1038/s41576-018-0016-z
- Wakefield J. Bayes factors for genome-wide association studies: comparison with P-values. Genetic Epidemiology. 2009;33(1):79-86. doi:10.1002/gepi.20359
- Benner C, Spencer CC, Havulinna AS, Salomaa V, Ripatti S, Pirinen M. FINEMAP: efficient variable selection using summary data from genome-wide association studies. Bioinformatics. 2016;32(10):1493-1501. doi:10.1093/bioinformatics/btw018
- Wang G, Sarkar A, Carbonetto P, Stephens M. A simple new approach to variable selection in regression, with application to genetic fine mapping. Journal of the Royal Statistical Society Series B. 2020;82(5):1273-1300. doi:10.1111/rssb.12388
- Giambartolomei C, Vukcevic D, Schadt EE, Franke L, Hingorani AD, Wallace C, Plagnol V. Bayesian test for colocalisation between pairs of genetic association studies using summary statistics. PLoS Genetics. 2014;10(5):e1004383. doi:10.1371/journal.pgen.1004383
- Kanai M, Elzur R, Zhou W, Daly MJ, Finucane HK. Meta-analysis fine-mapping is often miscalibrated at single-variant resolution. Cell Genomics. 2022;2(12):100210. doi:10.1016/j.xgen.2022.100210
- Finucane HK, Bulik-Sullivan B, Gusev A, et al. Partitioning heritability by functional annotation using genome-wide association summary statistics. Nature Genetics. 2015;47(11):1228-1235. doi:10.1038/ng.3404
- Cerezo M, Sollis E, Ji Y, et al. The NHGRI-EBI GWAS Catalog: standards for reusability, sustainability and diversity. Nucleic Acids Research. 2025;53(D1):D998-D1005. doi:10.1093/nar/gkae1070
- NHGRI-EBI GWAS Catalog. Preparing summary statistics files (GWAS-SSF format). https://www.ebi.ac.uk/gwas/docs/summary-statistics-format
- Gel B, Diez-Villanueva A, Serra E, Buschbeck M, Peinado MA, Malinverni R. regioneR: an R/Bioconductor package for the association analysis of genomic regions based on permutation tests. Bioinformatics. 2016;32(2):289-291. doi:10.1093/bioinformatics/btv562
- Ulirsch JC, Lareau CA, Bao EL, Ludwig LS, Guo MH, Benner C, Satpathy AT, Kartha VK, Salem RM, Hirschhorn JN, Finucane HK, Aryee MJ, Buenrostro JD, Sankaran VG. Interrogation of human hematopoiesis at single-cell and single-variant resolution. Nature Genetics. 2019;51(4):683-693. doi:10.1038/s41588-019-0362-6
- gchromVAR (version 0.3.2). https://github.com/caleblareau/gchromVAR
- Schep AN, Wu B, Buenrostro JD, Greenleaf WJ. chromVAR: inferring transcription-factor-associated accessibility from single-cell epigenomic data. Nature Methods. 2017;14(10):975-978. doi:10.1038/nmeth.4401
- Yu F, Cato LD, Weng C, Liggett LA, Jeon S, Xu K, Chiang CWK, Wiemels JL, Weissman JS, de Smith AJ, Sankaran VG. Variant to function mapping at single-cell resolution through network propagation. Nature Biotechnology. 2022;40(11):1644-1653. doi:10.1038/s41587-022-01341-y
- Chiou J, Geusz RJ, Okino ML, Han JY, Miller M, Melton R, Beebe E, Benaglio P, Huang S, Korgaonkar K, Heller S, Kleger A, Preissl S, Gorkin DU, Sander M, Gaulton KJ. Interpreting type 1 diabetes risk with genetics and single-cell epigenomics. Nature. 2021;594(7863):398-402. doi:10.1038/s41586-021-03552-w
- Corces MR, Shcherbina A, Kundu S, et al. Single-cell epigenomic analyses implicate candidate causal variants at inherited risk loci for Alzheimer’s and Parkinson’s diseases. Nature Genetics. 2020;52(11):1158-1168. doi:10.1038/s41588-020-00721-x
- Stuart T, Srivastava A, Madad S, Lareau CA, Satija R. Single-cell chromatin state analysis with Signac. Nature Methods. 2021;18(11):1333-1341. doi:10.1038/s41592-021-01282-5
- Fornes O, Castro-Mondragon JA, Khan A, et al. JASPAR 2020: update of the open-access database of transcription factor binding profiles. Nucleic Acids Research. 2020;48(D1):D87-D92. doi:10.1093/nar/gkz1001
- Akimov VE, Tychinin DI, Antonova OA, et al. Remodeling of the chromatin landscape in peripheral blood cells in patients with severe Delta COVID-19. Frontiers in Immunology. 2024;15:1415317. doi:10.3389/fimmu.2024.1415317 (scATAC-seq data: NCBI GEO GSE282769)
This tutorial completes the NGS101.com single-cell ATAC-seq analysis series for beginners.
Leave a Reply