We apply SME to hematology traits of white British individuals from
the UK Biobank. After quality control, the remaining data are 349,411
individuals and 543,813 SNPs of common variants. We select the traits
mean corpuscular hemoglobin (MCH), mean corpuscular hemoglobin
concentration, mean corpuscular volume (MCV), and hematocrit (HCT). As
external sparse data source, we leverage DNase I-hypersensitive sites
(DHSs) data measured over 12 days of ex-vivo erythroid differentiation
(Georgolopoulos et al. 2024). DHS are enriched in transcriptional
activity and used to identify regulatory DNA.
The first three traits, MCH, MCHC, and MCV are traits of the red
blood cells (RBC). Previous GWAS studies found that genes associated
with these traits are implicated in erythroid differentiation.
Therefore, we expect that genomic data that indicates regulatory regions
gathered during erythropoiesis will be informative for these traits. HCT
measures the percentage of red blood cells in the blood. The maturation
of erythroid progenitor cells is regulated by an oxygen-sensing
mechanism. We hypothesise that HCT, is not informed by functional data
of erythropoiesis.
Mask File Preparation
The external data sources used in this study represented genomic
intervals of DHS regions and LD blocks. In the following we mock this
data to illustrate how to create a mask file for sme(). See
How
To Create a Mask File for more details.
bim_data <- data.frame(
chromosome = c(1, 1, 1, 2, 2, 2, 3, 3, 3),
variant_id = c("rs1", "rs2", "rs3", "rs4", "rs5", "rs6", "rs7", "rs8", "rs9"),
cm_position = c(0, 0, 0, 0, 0, 0, 0, 0, 0),
bp_position = c(10, 20, 30, 40, 50, 60, 70, 80, 90),
allele1 = c("A", "A", "A", "G", "C", "C", "T", "T", "A"),
allele1 = c("G", "G", "G", "A", "T", "T", "A", "A", "G")
)
bim_data$index <- 1:nrow(bim_data)
# DHS intervals
hg19_dhs_regions <- data.frame(
chromosome = c(1, 2, 3),
start = c(5, 45, 85),
stop = c(15, 55, 95)
)
# LD block intervals
hg19_ld_blocks <- data.frame(
chromosome = c(1, 1, 2, 2, 3, 3, 3),
start = c(5, 25, 35, 45, 65, 75, 85),
stop = c(25, 35, 45, 65, 75, 85, 95)
)
We make use of the package GenomicRanges to efficiently
map the PLINK data to the intervals of the DHS and LD
data.
# Convert .bim to GRanges object
bim_gr <- GRanges(
seqnames = paste0("chr", bim_data$chromosome),
ranges = IRanges(start = bim_data$bp_position, end = bim_data$bp_position),
variant_id = bim_data$variant_id,
genome = "hg19"
)
# Convert DHS to GRanges object
dhs_gr <- GRanges(
seqnames = paste0("chr", hg19_dhs_regions$chromosome),
ranges = IRanges(start = hg19_dhs_regions$start, end = hg19_dhs_regions$stop),
genome = "hg19"
)
# Find overlaps of BIM variants and DHS intervals
overlaps <- findOverlaps(bim_gr, dhs_gr, maxgap = 0)
# Extract overlapping variants
dhs_data <- bim_data[queryHits(overlaps), ]
dhs_data <- dhs_data[!duplicated(dhs_data$index), ]
Of the 543,813 SNPs in our data, 4,932 are located in the DHS
regions. The DHS regions in this data are distributed along the whole
genome. To test for marginal epistasis with SME we consider the variants
in the DHS regions important.
Next we map the PLINK data to the LD blocks.
# Convert to GRanges object
ld_gr <- GRanges(
seqnames = paste0("chr", hg19_ld_blocks$chromosome),
ranges = IRanges(start = hg19_ld_blocks$start, end = hg19_ld_blocks$stop),
genome = "hg19"
)
# Find LD block of bim variants
ld_overlaps <- findOverlaps(query = bim_gr, subject = ld_gr)
With these objects, we can create the mask file. With larger data, we
recommend splitting the PLINK variants that are analyzed
into batches, create one mask file per batch, and submit one job per
batch on a High Peformance Cluster.
output_file <- tempfile()
gxg_group <- "gxg"
ld_group <- "ld"
gxg_variants <- dhs_data$index - 1 # 0-base index for C++
create_hdf5_file(output_file)
for (j in bim_data$index - 1) { # 0-base index for C++
# Write DHS mask
gxg_ds <- sprintf("%s/%d", gxg_group, j)
write_hdf5_dataset(file_name = output_file,
dataset_name = gxg_ds,
gxg_variants)
# Find LD block of focal SNP
focal_gr <- ld_gr[subjectHits(ld_overlaps[j,])]
# Find variants in LD block of focal SNP
focal_ld <- findOverlaps(query = bim_gr, subject = focal_gr)
ld_data <- bim_data[queryHits(focal_ld),]
ld_variants <- ld_data$index - 1 # 0-base index for C++
# Write LD mask
ld_ds <- sprintf("%s/%d", ld_group, j)
write_hdf5_dataset(file_name = output_file,
dataset_name = ld_ds,
ld_variants)
}
dhs_indices <- read_hdf5_dataset(file_name = output_file, dataset_name = gxg_ds)
print(sprintf("DHS indices: %s", paste(dhs_indices, collapse = ", ")))
#> [1] "DHS indices: 0, 4, 8"
With this mask file we can run SME.
sme_result <- sme(
plink_file = "/path/to/plink/data",
pheno_file = "/path/to/pheno/data",
mask_file = "/path/to/mask/file",
gxg_indices = c(1, 2, 3),
chunk_size = 250,
n_randvecs = 10,
n_blocks = 200,
n_threads = 6
)
The genome-wide association test for marginal epistasis in the red
blood cell traits MCH and MCV finds genome-wide significant statistical
epistasis (\(P
< \num{5e-8}\)) on chromosome 6 (Fig. 1). Importantly, most of
the SNPs or the genes which they map to have been previously discovered
for non-additive gene action related to erythropoiesis and RBC
traits.
