Using loopQTL

Introduction

loopQTL identifies SNPs associated with chromatin interaction signals across samples. It implements a band-wise joint normalization approach (referenced from Brant et al[^Brant]) that systematically tests the association between variants and their anchored loops, together with downstream functions for annotation and visualization for the user’s convenience.

This vignette walks through the basic loopQTL workflow.

[^Brant] Fletez-Brant K, Qiu Y, Gorkin DU, Hu M, Hansen KD. Removing unwanted variation between samples in Hi-C experiments. Brief Bioinform. 2024 Mar 27;25(3):bbae217. doi: 10.1093/bib/bbae217. PMID: 38711367; PMCID: PMC11074651.

Acquiring and installing loopQTL and associated packages

System requirements

R 4.6 or a newer version is required.

FIXME: update for 26/04 - recheck

On Ubuntu 24.04, successful installation required several additional packages. They can be installed by running the following command in a terminal:

sudo apt-get install r-base-dev libcurl4-openssl-dev libfreetype-dev libfribidi-dev libfontconfig-dev libharfbuzz-dev libtiff-dev libxml2-dev libssl-dev

No dependencies other than R are required on macOS or Windows.

Installation

Two versions of the package are available.

To install the stable version from Bioconductor, start R and run:

## Install BiocManager, which is required to install packages from Bioconductor
if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}

BiocManager::install("loopQTL")

The development version containing the most recent updates is available from our GitHub repository (https://github.com/rhielab/loopQTL).

To install the development version from GitHub, start R and run:

## Install prerequisite packages to install the development version from GitHub
if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}
if (!requireNamespace("remotes", quietly = TRUE)) {
    install.packages("remotes")
}

BiocManager::install(version = "devel")
BiocManager::install("rhielab/loopQTL")

Overview

  1. loopQTLsetup - Validate the multi-sample VCF header and per-sample .hic metadata, read per-sample loop BEDPEs, and build the exact-coordinate consensus loop set kept in loops(obj).
  2. loopQTLnormalize - Per-chromosome whole-genome straw pull followed by band-wise quantile normalization across samples The output is saved in phenotype(obj). Optionally caches per-SNP contact-window matrices to a file for later visualization with loopQTLplotSnpWindow() or loopQTLnormalize().
  3. loopQTLchooseK - Runs PCA on the normalized phenotype matrix and picks K via PCAForQTL’s method. This function returns the chosen K and the corresponding sample PCs.
  4. runQTL - Tabix-reads anchor-window SNPs from the VCF, runs MatrixEQTL to perform the loopQTL analysis. The output is saved in results(obj).
  5. loopQTLpairData - Return the sample-level phenotype and genotype for a single (loopID, snpID) pair - the row-level view behind any per-pair downstream plot or custom test.
  6. loopQTLannotate - Adds gene, peak (optional), and Type columns to results(obj) using a GTF for TSS mapping and an optional peak BED file for regulatory-element overlap.
  7. loopQTLplotPair - Generates boxplot of the loop’s contact-strength phenotype grouped by the SNP’s dosage class (0/1/2), for a specific SNP-loop pair.
  8. loopQTLplotLocus - Locus-level ggplot track combining consensus loops, the tested SNP, and (optionally) gene annotations, over a user-specified genomic region.
  9. loopQTLplotSnpWindow - Hi-C-style beta heatmap around a chosen SNP: for every bin-pair in the SNP’s neighborhood, plots the per-sample regression coefficient of phenotype on the SNP’s dosage.

To begin, load loopQTL and any downstream plotting deps.

library(loopQTL)
library(ggplot2)

loopQTLsetup: Validate inputs and build loopQTL object

loopQTLsetup() is the entry point of the pipeline. It:

  • validates the input files header, including .vcf and .hic files;
  • reads each per-sample BEDPE and builds a consensus loop set by exact-coordinate matching. Loops found fewer than minSamples samples will be dropped;
  • Automatically computes per-sample cis-contact depth from the .hic files when the samplesheet lacks a cis column.

The samplesheet is a tab-separated file with columns Sample, hic, loops, and optionally cis.

extdataDir <- system.file("extdata", package = "loopQTL")
vcfPath <- file.path(extdataDir, "example.vcf.gz")
samplesheetPath <- file.path(extdataDir, "example_samplesheet.tsv")
samples <- read.table(samplesheetPath, header = TRUE)

## Fix the input file paths to be absolute since we don't know what the working
## directory will be when the vignette is built. It is not necessary to do this
## during normal use of loopQTL if the working directory is set correctly.
samples$hic <- unlist(
    lapply(samples$hic, function(x) {
        file.path(extdataDir, x)
    })
)
samples$loops <- unlist(
    lapply(samples$loops, function(x) {
        file.path(extdataDir, x)
    })
)

## Display the format of the samplesheet
head(samples)
#>     Sample
#> 1  sample1
#> 2  sample3
#> 3  sample5
#> 4  sample7
#> 5 sample10
#> 6 sample12
#>                                                                          hic
#> 1  /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example_hic/sample1.hic
#> 2  /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example_hic/sample3.hic
#> 3  /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example_hic/sample5.hic
#> 4  /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example_hic/sample7.hic
#> 5 /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example_hic/sample10.hic
#> 6 /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example_hic/sample12.hic
#>                                                                                  loops
#> 1  /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example_loops/sample1_loops.bedpe
#> 2  /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example_loops/sample3_loops.bedpe
#> 3  /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example_loops/sample5_loops.bedpe
#> 4  /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example_loops/sample7_loops.bedpe
#> 5 /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example_loops/sample10_loops.bedpe
#> 6 /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example_loops/sample12_loops.bedpe
#>         cis
#> 1  62392910
#> 2 106429140
#> 3 309557957
#> 4  91097755
#> 5 294491793
#> 6  68389886
## Run loopQTLsetup
obj <- loopQTLsetup(
    samples    = samples,
    vcf        = vcfPath,
    resolution = 10000,
    minSamples = 2
)
#> Proceeding with 10 sample(s) present in both samplesheet and VCF.
#> exact-coord consensus across 10 samples ...
#> 979 unique loop tuples; 282 survive minSamples = 2
obj
#> loopQTL object:
#>  samples   : 10 ( sample1, sample3, sample5, sample7 , ... )
#>  loops     : 282 consensus intra-chromosomal loops
#>  resolution: 10000 bp
#>  phenotype : not yet computed 
#>  results   : not yet computed 
#>  meta      :
#>       resolution = 10000 
#>       minSamples = 2 
#>       n_trans_dropped = 0 
#>       n_trans_per_sample = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 
#>       cisSource = userProvided 
#>       cis_chroms_in = chr1 
#>       normalization = NONE 
#>       hicPaths = /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example_hic/sample1.hic, /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example_hic/sample3.hic, /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example_hic/sample5.hic, /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example_hic/sample7.hic, /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example_hic/sample10.hic, /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example_hic/sample12.hic, /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example_hic/sample14.hic, /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example_hic/sample17.hic, /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example_hic/sample20.hic, /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example_hic/sample23.hic

loopQTLnormalize: Cross-normalization on Hi-C signal

loopQTLnormalize() normalizes contact strength across samples so that per-loop phenotypes are directly comparable. For each chromosome it pulls the whole-genome straw matrix per sample, log-CPM adjusts, then does band-wise quantile normalization across samples per genomic-distance band. Normalized values at the consensus-loop bin-pairs are written to phenotype(obj). Parallelizes per chromosome via BiocParallel.

Set saveWindows = TRUE to additionally output a per-SNP contact-window matrix as .rds file. This file is needed by loopQTLplotSnpWindow() downstream. Please note that this option does not change phenotype(obj) or affect the runQTL() step.

## Get the path to a temporary directory in which to save the output interact
## file
tempDirectory <- tempdir()

obj <- loopQTLnormalize(
    obj,
    BPPARAM     = BiocParallel::MulticoreParam(workers = 2),
    saveWindows = TRUE,
    windowsOut  = file.path(tempDirectory, "loop_snp_windows.rds"),
    window      = 2000000,
    snpWindow   = 10000
)
#> saveWindows: tabix-reading SNPs at 35 anchor regions (snpWindow = 10000 bp) ...
#> saveWindows: 39 SNP(s), 27210 unique cells across 1 chrom(s); estimated ~2.1 MB in memory (8 bytes x cells x samples)
#> saveWindows: wrote /tmp/RtmpRg9f6I/loop_snp_windows.rds (0.4 MB)

head(phenotype(obj)[, 1:5])
#>        sample1   sample3   sample5   sample7  sample10
#> [1,] -2.567672 -2.612742 -2.307309 -3.015544 -2.194319
#> [2,] -2.575766 -2.440413 -2.277193 -2.619341 -2.462711
#> [3,] -3.627543 -2.499697 -2.661846 -2.202957 -2.365405
#> [4,] -2.325173 -2.012506 -2.262684 -2.618017 -2.188440
#> [5,] -2.172129 -2.513661 -2.382541 -2.993291 -2.266896
#> [6,] -2.897595 -2.780669 -2.342360 -2.502227 -2.290136

loopQTLchooseK: Auto-select the number of PCs for covariates

loopQTLchooseK() runs PCA on the normalized phenotype matrix and picks a number of principal components K to use as covariates in the downstream regression. This function returns a list with the chosen K and the sample × K matrix of PCs.

Kres <- loopQTLchooseK(obj, Kmethod = "elbow")
message("PCAForQTL chose K = ", Kres$K, " via ", Kres$method)
#> PCAForQTL chose K = 2 via elbow
head(Kres$PCs)
#>                PC1       PC2
#> sample1  -9.816768  5.526077
#> sample3  -8.345106  7.758669
#> sample5   9.981544  2.972938
#> sample7  -9.375970 -9.421029
#> sample10  9.985602 -1.938478
#> sample12 -6.847075 -3.910076

runQTL: Anchor-restricted association testing

runQTL() reads all anchor-window SNPs from the VCF in a single targeted tabix pass, converts them to dosages (0/1/2), and runs MatrixEQTL twice - once treating anchor A as the gene position, once anchor B. The two result tables are outer-joined on SNP-loop level. Covariates are supplied via the covariates argument - typically the phenotype PCs chosen in the previous step. RINT is applied per loop before regression by default.

obj <- runQTL(
    obj,
    window     = 10000,
    covariates = Kres$PCs
)
#> 282 of 282 genes matched
#> 39 of 39 SNPs matched
#> 100.00% done, 192 cis-eQTLs
#> 282 of 282 genes matched
#> 39 of 39 SNPs matched
#> 100.00% done, 18 cis-eQTLs
head(results(obj))
#>                            loopID     snpID chrom      pos ref alt anchor
#> 1 L0000187_chr1_46030000_46220000 rs3795318  chr1 46029352   G   A      A
#> 2 L0000179_chr1_46020000_46230000 rs3795318  chr1 46029352   G   A      A
#> 3 L0000177_chr1_46020000_46180000 rs3795318  chr1 46029352   G   A      A
#> 4 L0000176_chr1_46020000_46170000 rs3795318  chr1 46029352   G   A      A
#> 5 L0000175_chr1_46020000_46140000 rs3795318  chr1 46029352   G   A      A
#> 6 L0000184_chr1_46030000_46180000 rs3795318  chr1 46029352   G   A      A
#>         beta        se           t           p nUsed  maf
#> 1 -0.9445405 0.2490755 -3.79218596 0.009050187    10 0.25
#> 2 -0.4646524 0.4237463 -1.09653429 0.314893409    10 0.25
#> 3 -0.0124666 0.3851023 -0.03237217 0.975225264    10 0.25
#> 4  0.1563867 0.4240769  0.36876967 0.724961982    10 0.25
#> 5 -0.1409811 0.5063318 -0.27843620 0.790024483    10 0.25
#> 6  0.1397946 0.4565184  0.30621896 0.769785411    10 0.25

loopQTLpairData: Sample-level view of one loop × SNP pair

loopQTLpairData() returns the per-sample phenotype and genotype (dosage 0/1/2) for a single SNP-loop pair. Useful for inspecting a specific hit or exporting sample-level values for reporting.

pair <- loopQTLpairData(
    obj,
    loopID = "L0000187_chr1_46030000_46220000",
    snpID  = "rs11584814"
)
head(pair)
#>     sample                          loopID      snpID genotype  phenotype
#> 1  sample1 L0000187_chr1_46030000_46220000 rs11584814        1 -1.0364334
#> 2  sample3 L0000187_chr1_46030000_46220000 rs11584814        0  0.3853205
#> 3  sample5 L0000187_chr1_46030000_46220000 rs11584814        0  1.6448536
#> 4  sample7 L0000187_chr1_46030000_46220000 rs11584814        0 -1.6448536
#> 5 sample10 L0000187_chr1_46030000_46220000 rs11584814        1  0.6744898
#> 6 sample12 L0000187_chr1_46030000_46220000 rs11584814        0 -0.1256613
#>         cis
#> 1  62392910
#> 2 106429140
#> 3 309557957
#> 4  91097755
#> 5 294491793
#> 6  68389886

loopQTLannotate: Annotate results with target genes and regulatory peaks

loopQTLannotate() adds three columns to results(obj):

  • gene - comma-separated gene names whose TSS (± TSSwindow) sits in the SNP’s partner anchor.
  • peak - logical: does the SNP overlap any peak from the user-supplied peak file . If the distal TSS peaks are preferred, then you can set `removeTSS = TRUE``.
  • Type - annotation of each SNP-loop pair such as WithinTSS, rs-gene, rs-<peakAnnotation>, rs-O, rs(<peakAnnotation>)-gene, or rs(<peakAnnotation>)-<peakAnnotation> summarizing the SNP–loop types.
peaksPath <- system.file("extdata", "example_annotation.bed", package = "loopQTL")
gtfPath <- system.file("extdata", "example.gtf.gz", package = "loopQTL")
obj <- loopQTLannotate(
    obj,
    peak            = peaksPath,
    peakAnnotation  = "Enhancer",
    TSSwindow       = 2000,
    removeTSS       = TRUE,
    gtf             = gtfPath
)
#> reading GTF: /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example.gtf.gz
#> 189 gene TSS regions (±2000 bp)
#> reading peaks: /tmp/RtmpSRlPA1/Rinst10c367521ed4/loopQTL/extdata/example_annotation.bed
#> dropped 60 peaks overlapping TSS regions; 3218 peaks retained for SNP overlap.
head(results(obj))
#>                            loopID     snpID chrom      pos ref alt anchor
#> 1 L0000187_chr1_46030000_46220000 rs3795318  chr1 46029352   G   A      A
#> 2 L0000179_chr1_46020000_46230000 rs3795318  chr1 46029352   G   A      A
#> 3 L0000177_chr1_46020000_46180000 rs3795318  chr1 46029352   G   A      A
#> 4 L0000176_chr1_46020000_46170000 rs3795318  chr1 46029352   G   A      A
#> 5 L0000175_chr1_46020000_46140000 rs3795318  chr1 46029352   G   A      A
#> 6 L0000184_chr1_46030000_46180000 rs3795318  chr1 46029352   G   A      A
#>         beta        se           t           p nUsed  maf
#> 1 -0.9445405 0.2490755 -3.79218596 0.009050187    10 0.25
#> 2 -0.4646524 0.4237463 -1.09653429 0.314893409    10 0.25
#> 3 -0.0124666 0.3851023 -0.03237217 0.975225264    10 0.25
#> 4  0.1563867 0.4240769  0.36876967 0.724961982    10 0.25
#> 5 -0.1409811 0.5063318 -0.27843620 0.790024483    10 0.25
#> 6  0.1397946 0.4565184  0.30621896 0.769785411    10 0.25
#>                           gene  peak    Type
#> 1                      POMGNT1 FALSE rs-gene
#> 2                              FALSE    rs-O
#> 3 TSPAN1,AL672043.1,AL358075.4 FALSE rs-gene
#> 4 TSPAN1,AL672043.1,AL358075.4 FALSE rs-gene
#> 5            PIK3R3,AL358075.2 FALSE rs-gene
#> 6 TSPAN1,AL672043.1,AL358075.4 FALSE rs-gene

loopQTLplotPair: Boxplot of loop strength by SNP dosage

loopQTLplotPair() renders the per-sample boxplot behind a chosen result SNP-loop pair. Please note that the plot is a ggplot object that you can save with ggsave().

p <- loopQTLplotPair(
    obj,
    loopID = "L0000187_chr1_46030000_46220000",
    snpID  = "rs11584814"
)
p

loopQTLplotLocus: Locus-level track view of loops, SNPs, and genes

loopQTLplotLocus() renders a locus-scale ggplot combining the consensus loops, the tested SNP, and (when a GTF is supplied) gene annotations, over a user-specified region. Please note that to sse values = "cpm" or values = "raw" to show per-loop contact strength, you need to set saveWindows = TRUE to save the per-SNP contact-window matrices from loopQTLnormalize.

gtfPath <- system.file("extdata", "example.gtf.gz", package = "loopQTL")
p <- loopQTLplotLocus(
    obj,
    chrom = "chr1",
    start = 45900000,
    end = 46300000,
    snpID = "rs11584814",
    values = "cpm",
    gtf = gtfPath
)
p

loopQTLplotSnpWindow: Hi-C-style beta heatmap around a SNP

loopQTLplotSnpWindow() renders a Hi-C-style heatmap centered on a chosen SNP: for every bin-pair in the SNP’s cached contact window (from loopQTLnormalize(saveWindows = TRUE)), the color encodes the per-sample regression coefficient (beta) of the bin-pair’s contact strength on the SNP’s dosage. Set normal = TRUE and covariates = "auto" to match the model used in runQTL().

p <- loopQTLplotSnpWindow(
    obj,
    chrom      = "chr1",
    start      = 45900000,
    end        = 46300000,
    snpID      = "rs11584814",
    normal     = TRUE,
    covariates = "auto",
    Kmax       = 2
)

## Delete the output file, since this is just an example.
## Do not use this line if running on real data, as it will delete your created
## file.
invisible(file.remove(metadata(obj)$windowsOut))

p

Session info

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] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] ggplot2_4.0.3    loopQTL_0.99.0   BiocStyle_2.41.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] KEGGREST_1.53.6             SummarizedExperiment_1.43.0
#>  [3] gtable_0.3.6                rjson_0.2.23               
#>  [5] xfun_0.61                   bslib_0.12.0               
#>  [7] Biobase_2.73.2              lattice_0.23-1             
#>  [9] vctrs_0.7.3                 tools_4.6.1                
#> [11] bitops_1.1-0                generics_0.1.4             
#> [13] curl_8.0.0                  stats4_4.6.1               
#> [15] parallel_4.6.1              AnnotationDbi_1.75.2       
#> [17] RSQLite_3.53.3              blob_1.3.0                 
#> [19] BiocBaseUtils_1.15.1        Matrix_1.7-6               
#> [21] BSgenome_1.81.1             data.table_1.18.6.1        
#> [23] RColorBrewer_1.1-3          cigarillo_1.3.1            
#> [25] S7_0.2.2                    S4Vectors_0.51.10          
#> [27] lifecycle_1.0.5             compiler_4.6.1             
#> [29] farver_2.1.2                Rsamtools_2.29.0           
#> [31] Biostrings_2.81.9           Seqinfo_1.3.2              
#> [33] codetools_0.2-20            InteractionSet_1.41.0      
#> [35] htmltools_0.5.9             sys_3.4.3                  
#> [37] buildtools_1.0.0            sass_0.4.10                
#> [39] RCurl_1.98-1.20             yaml_2.3.12                
#> [41] crayon_1.5.3                jquerylib_0.1.4            
#> [43] MatrixEQTL_2.4              BiocParallel_1.47.0        
#> [45] DelayedArray_0.39.6         cachem_1.1.0               
#> [47] abind_1.4-8                 digest_0.6.39              
#> [49] restfulr_0.0.17             VariantAnnotation_1.59.4   
#> [51] labeling_0.4.3              maketools_1.3.2            
#> [53] fastmap_1.2.0               grid_4.6.1                 
#> [55] cli_3.6.6                   SparseArray_1.13.2         
#> [57] patchwork_1.3.2             S4Arrays_1.13.0            
#> [59] GenomicFeatures_1.65.0      XML_3.99-0.24              
#> [61] withr_3.0.3                 scales_1.4.0               
#> [63] bit64_4.8.6                 rmarkdown_2.32             
#> [65] XVector_0.53.0              httr_1.4.9                 
#> [67] matrixStats_1.5.0           bit_4.6.0                  
#> [69] png_0.1-9                   memoise_2.0.1              
#> [71] evaluate_1.0.5              knitr_1.52                 
#> [73] BiocIO_1.23.3               GenomicRanges_1.65.4       
#> [75] IRanges_2.47.5              rtracklayer_1.73.0         
#> [77] rlang_1.3.0                 Rcpp_1.1.2                 
#> [79] glue_1.8.1                  DBI_1.3.0                  
#> [81] BiocManager_1.30.27         BiocGenerics_0.59.12       
#> [83] jsonlite_2.0.0              strawr_0.0.92              
#> [85] R6_2.6.1                    GenomicAlignments_1.49.2   
#> [87] MatrixGenerics_1.25.0