hmagmaR links non-coding GWAS SNPs to their putative
target genes by combining chromatin interaction data (e.g. PLAC-seq or
Hi-C) with regulatory region annotations (e.g. promoters, enhancers,
H3K27ac peaks). The output is a MAGMA-compatible gene annotation file,
which can then be used for gene-level association testing.
The typical workflow has two stages:
AnnotationFileHmagma() (or its
down-sampled variant, SampledDownAnnotation()) builds the
gene annotation file by overlapping chromatin interactions with
regulatory regions and SNPs.GeneLevelAnalysis_hmagma() runs MAGMA
itself on that annotation file against a GWAS summary statistics file,
using an external MAGMA binary and 1000 Genomes reference panel.This vignette demonstrates stage 1 using a small example dataset
bundled with the package. Stage 2 requires external software and
reference files not included with hmagmaR, so it is shown
for reference only and not executed here.
hmagmaR ships with a small example dataset,
test_dataset.microglia.hg19, subsampled from microglia
PLAC-seq interactions (Nott et al., 2019), microglia promoter/enhancer
regions, 1000 Genomes SNPs, and a GWAS summary statistics table, so the
workflow below runs in seconds rather than requiring genome-wide
reference data.
data("test_dataset.microglia.hg19")
hic[1:3, ]
#> # A tibble: 3 × 14
#> chr1 start1 end1 chr2 start2 end2 count expected fdr ClusterLabel
#> <chr> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
#> 1 chr19 58105000 5.81e7 chr19 5.87e7 5.87e7 24 4.05 2.11e-9 chr19_484
#> 2 chr8 80985000 8.10e7 chr8 8.11e7 8.11e7 24 8.95 7.37e-4 chr8_348
#> 3 chr10 64995000 6.5 e7 chr10 6.53e7 6.53e7 13 3.53 1.40e-3 chr10_720
#> # ℹ 4 more variables: ClusterSize <dbl>, ClusterType <chr>,
#> # ClusterNegLog10P <dbl>, ClusterSummit <dbl>
promoterRegions[1:3, ]
#> chr start end
#> 1 chr1 713696 715319
#> 2 chr1 761972 763164
#> 3 chr1 761972 763164
enhancerRegions[1:3, ]
#> chr start end
#> 1 chr1 751631 751770
#> 2 chr1 752191 752824
#> 3 chr1 753253 753771
snps[1:3, ]
#> chr Position SNP
#> 1 chr8 4047042 rs6558853
#> 2 chr5 72907530 rs539414285
#> 3 chr8 84164952 rs148344334
annotated_genes[1:3, ]
#> chr start end ensg
#> 1 6 75962638 75994632 ENSG00000112697
#> 2 2 101869345 101886778 ENSG00000158435
#> 3 11 1686829 1689086 ENSG00000291437
snpgeneexon[1:3, ]
#> rsid ensg
#> 1 rs569385334 ENSG00000290976
#> 2 rs3774854 ENSG00000173040
#> 3 rs572555182 ENSG00000109654hic: chromatin interaction data with columns
chrom1, start1, end1,
chrom2, start2, end2 (plus
additional interaction-calling metadata).promoterRegions / enhancerRegions:
regulatory regions in BED-like format (chr,
start, end).snps: reference SNPs (e.g. from a .bim
file), with columns chr, Position,
SNP.annotated_genes: reference gene coordinates, with
columns chr, start, end,
ensg.snpgeneexon: SNPs falling within exons/promoters of
each gene, with columns rsid, ensg.AnnotationFileHmagma() accepts either a single
regulatoryRegions table, or separate
promoterRegions and enhancerRegions tables (in
which case promoters anchor the chromatin interactions, and enhancers
filter which SNPs are retained). Here we use the latter:
annotation_file <- tempfile("hmagmaR_example")
AnnotationFileHmagma(
hic = hic,
promoterRegions = promoterRegions,
enhancerRegions = enhancerRegions,
snps = snps,
annotated_genes = annotated_genes,
snpgeneexon = snpgeneexon,
AnnotationFile = annotation_file
)
#> Using promoterRegions and enhancerRegions
#> >> preparing features information... 2026-10-12 01:44:27 PM
#> >> Using Genome: hg19 ...
#> >> identifying nearest features... 2026-10-12 01:44:29 PM
#> >> calculating distance from peak to TSS... 2026-10-12 01:44:30 PM
#> >> assigning genomic annotation... 2026-10-12 01:44:30 PM
#> >> Using Genome: hg19 ...
#> >> Using Genome: hg19 ...
#> >> adding gene annotation... 2026-10-12 01:45:03 PM
#> Loading required package: org.Hs.eg.db
#> Loading required package: AnnotationDbi
#> Loading required package: stats4
#> Loading required package: BiocGenerics
#> Loading required package: generics
#>
#> Attaching package: 'generics'
#> The following objects are masked from 'package:base':
#>
#> as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
#> setequal, union
#>
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:stats':
#>
#> IQR, mad, sd, var, xtabs
#> The following object is masked from 'package:utils':
#>
#> data
#> The following objects are masked from 'package:base':
#>
#> anyDuplicated, aperm, append, as.data.frame, basename, cbind,
#> colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
#> get, grep, grepl, is.unsorted, lapply, Map, mapply, match, mget,
#> order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
#> rbind, Reduce, rownames, sapply, saveRDS, scale, sequence, table,
#> tapply, transform, unique, unsplit, which.max, which.min
#> Loading required package: Biobase
#> Welcome to Bioconductor
#>
#> Vignettes contain introductory material; view with
#> 'browseVignettes()'. To cite Bioconductor, see
#> 'citation("Biobase")', and for packages 'citation("pkgname")'.
#> Loading required package: IRanges
#> Loading required package: S4Vectors
#>
#> Attaching package: 'S4Vectors'
#> The following object is masked from 'package:utils':
#>
#> findMatches
#> The following objects are masked from 'package:base':
#>
#> expand.grid, I, unname
#> 'select()' returned 1:many mapping between keys and columns
#> >> assigning chromosome lengths 2026-10-12 01:45:04 PM
#> >> done... 2026-10-12 01:45:04 PM
#> Warning in .merge_two_Seqinfo_objects(x, y): Each of the 2 combined objects has sequence levels not in the other:
#> - in 'x': chr23
#> - in 'y': chrX, chrY
#> Make sure to always combine/compare objects based on the same reference
#> genome (use suppressWarnings() to suppress this warning).This writes a MAGMA-style annotation file
(<AnnotationFile>.transcript.annot), with one row per
gene: the Ensembl gene ID, its genomic coordinates, and the
tab-separated rsIDs of all SNPs linked to it.
result <- readLines(paste0(annotation_file, ".transcript.annot"))
length(result)
#> [1] 53
head(result, 5)
#> [1] "ENSG00000070831\t1:22379120:22419436\trs2501295\trs575745784\trs16860621"
#> [2] "ENSG00000074319\t11:18501858:18548503\trs75027358\trs1395320"
#> [3] "ENSG00000090104\t1:192544857:192549159\trs578121256\trs3806353"
#> [4] "ENSG00000102445\t13:46916137:47012325\trs188399416"
#> [5] "ENSG00000103187\t16:84599204:84651669\trs148207631"With this small example dataset, 53 genes were annotated with at least one interacting SNP. Because the bundled data is a random subsample of the full reference tables rather than a curated locus, this number is illustrative only - a real analysis using full genome-wide interaction, regulatory region, and SNP reference data will typically annotate many more genes.
Once an annotation file has been generated,
GeneLevelAnalysis_hmagma() calls MAGMA to perform
gene-level association testing against a GWAS summary statistics file.
This step requires:
.bed/.bim/.fam) for the
appropriate ancestry,SNP, P, and
N columns.GeneLevelAnalysis_hmagma(
magma = "/path/to/magma",
g1000 = "/path/to/g1000_eur",
gwas = "/path/to/gwas_summary_stats.txt",
AnnotationFile = paste0(annotation_file, ".transcript.annot"),
output = "hmagmaR_example"
)This produces MAGMA’s standard gene-level output files
(.genes.out, .genes.raw, .log),
which can be used for downstream gene set or pathway analyses.
SampledDownAnnotation())When comparing annotation results across samples with different
numbers of chromatin interactions, SampledDownAnnotation()
repeatedly subsamples the hic table down to a fixed number
of interactions (loopNumber), so that differences in
results reflect biology rather than differing input sizes. It otherwise
follows the same pipeline as AnnotationFileHmagma() with a
single combined regulatoryRegions table, writing one
annotation file per resample (suffixed .transcript.annot.1,
.transcript.annot.2, etc.).
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] org.Hs.eg.db_3.23.1 AnnotationDbi_1.75.2 IRanges_2.47.5
#> [4] S4Vectors_0.51.10 Biobase_2.73.2 BiocGenerics_0.59.12
#> [7] generics_0.1.4 hmagmaR_0.99.3 BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] RColorBrewer_1.1-3
#> [2] sys_3.4.3
#> [3] jsonlite_2.0.0
#> [4] tidydr_0.0.7
#> [5] magrittr_2.0.5
#> [6] ggtangle_0.1.4
#> [7] GenomicFeatures_1.65.0
#> [8] farver_2.1.2
#> [9] rmarkdown_2.32
#> [10] fs_2.1.0
#> [11] BiocIO_1.23.3
#> [12] vctrs_0.7.3
#> [13] memoise_2.0.1
#> [14] Rsamtools_2.29.0
#> [15] RCurl_1.98-1.20
#> [16] ggtree_4.3.0
#> [17] htmltools_0.5.9
#> [18] TxDb.Hsapiens.UCSC.hg19.knownGene_3.22.1
#> [19] BiocBaseUtils_1.15.1
#> [20] plotrix_3.8-14
#> [21] curl_8.0.0
#> [22] gridGraphics_0.5-1
#> [23] sass_0.4.10
#> [24] KernSmooth_2.23-27
#> [25] bslib_0.12.0
#> [26] htmlwidgets_1.6.4
#> [27] plyr_1.8.9
#> [28] cachem_1.1.0
#> [29] buildtools_1.0.0
#> [30] GenomicAlignments_1.49.2
#> [31] igraph_2.3.4
#> [32] lifecycle_1.0.5
#> [33] pkgconfig_2.0.3
#> [34] Matrix_1.7-6
#> [35] R6_2.6.1
#> [36] fastmap_1.2.0
#> [37] digest_0.6.39
#> [38] aplot_0.3.2
#> [39] enrichplot_1.99.8
#> [40] ggnewscale_0.5.2
#> [41] patchwork_1.3.2
#> [42] GenomicRanges_1.65.4
#> [43] RSQLite_3.53.3
#> [44] httr_1.4.9
#> [45] polyclip_1.10-7
#> [46] compiler_4.6.1
#> [47] bit64_4.8.6
#> [48] fontquiver_0.2.1
#> [49] withr_3.0.3
#> [50] S7_0.2.2
#> [51] BiocParallel_1.47.0
#> [52] DBI_1.3.0
#> [53] ggforce_0.5.0
#> [54] gplots_3.3.0
#> [55] MASS_7.3-66
#> [56] ChIPseeker_1.49.3
#> [57] rappdirs_0.3.4
#> [58] rjson_0.2.23
#> [59] caTools_1.18.4
#> [60] gtools_3.9.5
#> [61] tools_4.6.1
#> [62] otel_0.2.0
#> [63] ape_5.8-1
#> [64] scatterpie_0.2.6
#> [65] glue_1.8.1
#> [66] restfulr_0.0.17
#> [67] nlme_3.1-171
#> [68] GOSemSim_2.39.4
#> [69] grid_4.6.1
#> [70] cluster_2.1.8.3
#> [71] reshape2_1.4.5
#> [72] gtable_0.3.6
#> [73] tidyr_1.3.2
#> [74] utf8_1.2.6
#> [75] XVector_0.53.0
#> [76] ggrepel_0.9.8
#> [77] pillar_1.11.1
#> [78] stringr_1.6.0
#> [79] yulab.utils_0.2.5
#> [80] dplyr_1.2.1
#> [81] tweenr_2.0.3
#> [82] treeio_1.37.1
#> [83] lattice_0.23-1
#> [84] rtracklayer_1.73.0
#> [85] bit_4.6.0
#> [86] tidyselect_1.2.1
#> [87] fontLiberation_0.1.0
#> [88] GO.db_3.23.1
#> [89] maketools_1.3.2
#> [90] Biostrings_2.81.9
#> [91] knitr_1.52
#> [92] fontBitstreamVera_0.1.1
#> [93] Seqinfo_1.3.2
#> [94] xfun_0.61
#> [95] stringi_1.8.9
#> [96] UCSC.utils_1.9.0
#> [97] lazyeval_0.2.3
#> [98] ggfun_0.2.2
#> [99] yaml_2.3.12
#> [100] boot_1.3-32
#> [101] evaluate_1.0.5
#> [102] codetools_0.2-20
#> [103] cigarillo_1.3.1
#> [104] gdtools_0.5.1
#> [105] tibble_3.3.1
#> [106] BiocManager_1.30.27
#> [107] ggplotify_0.1.3
#> [108] cli_3.6.6
#> [109] systemfonts_1.3.2
#> [110] jquerylib_0.1.4
#> [111] Rcpp_1.1.2
#> [112] GenomeInfoDb_1.49.1
#> [113] png_0.1-9
#> [114] XML_3.99-0.25
#> [115] parallel_4.6.1
#> [116] ggplot2_4.0.3
#> [117] blob_1.3.0
#> [118] DOSE_4.7.3
#> [119] bitops_1.1-0
#> [120] tidytree_0.4.8
#> [121] ggiraph_0.9.6
#> [122] enrichit_0.2.6
#> [123] scales_1.4.0
#> [124] purrr_1.2.2
#> [125] crayon_1.5.3
#> [126] rlang_1.3.0
#> [127] KEGGREST_1.53.6