# Inferring Linkage Disequilibrium blocks from genotypes

#### Shubham Chaturvedi, Pierre Neuvial, Nathalie Vialaneix

#### 2025-11-19

    # IMPORTANT: this vignette is not created if snpStats is not installed
    if (!require("snpStats")) {
      knitr::opts_chunk$set(eval = FALSE)
    }

    ## Loading required package: snpStats

    ## Loading required package: survival

    ## Loading required package: Matrix

    ##
    ## Attaching package: 'Matrix'

    ## The following object is masked from 'package:S4Vectors':
    ##
    ##     expand

# Introduction

In this vignette we demonstrate the use of `snpClust` function in the `adjclust` package. `snpClust` performs adjacency-constrained hierarchical clustering of single nucleotide polymorphisms (SNPs), where the similarity between SNPs is defined by linkage disequilibrium (LD).

This function implements the algorithm described in [1]. It is an extension of the algorithm described in [3,4]. Denoting by \(p\) the number of SNPs to cluster and assuming that the similarity between SNPs whose indices are more distant than \(h\), its time complexity is \(O(p (\log(p) + h))\), and its space complexity is \(O(hp)\).

    library("adjclust")

# Loading and displaying genotype data

The beginning of this vignette closely follows the “LD vignette” of the SnpStats package [2]. First, we load genotype data.

    data("ld.example", package = "snpStats")

We focus on the `ceph.1mb` data.

    geno <- ceph.1mb[, -316]  ## drop one SNP leading to one missing LD value
    p <- ncol(geno)
    nSamples <- nrow(geno)
    geno

    ## A SnpMatrix with  90 rows and  602 columns
    ## Row names:  NA06985 ... NA12892
    ## Col names:  rs5993821 ... rs5747302

These data are drawn from the International HapMap Project and concern 602 SNPs1 over a 1Mb region of chromosome 22 in sample of 90 Europeans.

We can compute and display the LD between these SNPs.

    ld.ceph <- snpStats::ld(geno, stats = "R.squared", depth = p-1)
    image(ld.ceph, lwd = 0)

# Adjacency-constrained Hierarchical Agglomerative Clustering

The `snpClust` function can handle genotype data as an input:

    fit <- snpClust(geno, stats = "R.squared")

    ## Warning in run.snpClust(x, h = h, stats = stats): Forcing the LD similarity to
    ## be smaller than or equal to 1

    ## Note: 135 merges with non increasing heights.

Note that due to numerical errors in the LD estimation, some of the estimated LD values may be slightly larger than 1. These values are rounded to 1 internally.

The above figure suggests that the LD signal is concentrated close to the diagonal. We can focus on a diagonal band with the bandwidth parameter `h`:

    fitH <- snpClust(geno, h = 100, stats = "R.squared")

    ## Warning in run.snpClust(x, h = h, stats = stats): Forcing the LD similarity to
    ## be smaller than or equal to 1

    ## Note: 133 merges with non increasing heights.

    fitH

    ##
    ## Call:
    ## snpClust(geno, h = 100, stats = "R.squared")
    ##
    ## Cluster method   : snpClust
    ## Number of objects: 602

# Output

The output of the `snpClust` is of class `chac`. In particular, it can be plotted as a dendrogram silently relying on the function `plot.dendrogram`:

    plot(fitH, type = "rectangle", leaflab = "perpendicular")

    ## Warning:
    ## Detected reversals in dendrogram: mode = 'corrected', 'within-disp' or 'total-disp' might be more relevant.

Moreover, the output contains an element named `merge` which describes the successive merges of the clustering, and an element `gains` which gives the improvement in the criterion optimized by the clustering at each successive merge.

    head(cbind(fitH$merge, fitH$gains))

    ##      [,1] [,2]
    ## [1,]   -1   -2
    ## [2,] -255 -256
    ## [3,] -488 -489
    ## [4,] -487    3
    ## [5,] -486    4
    ## [6,] -234 -235

# Other types of input

In this section we show how the `snpClust` function can also be applied directly to LD values.

    h <- 100
    ld.ceph <- snpStats::ld(geno, stats = "R.squared", depth = h, symmetric = TRUE)
    image(ld.ceph, lwd = 0)

Note that we have forced the `snpStats::ld` function to return a symmetric matrix. We can apply `snpClust` directly to this LD matrix (of class `Matrix::dsCMatrix`):

    fitL <- snpClust(ld.ceph, h)

    ## Note: forcing the diagonal of the LD similarity matrix to be 1

    ## Warning in run.snpClust(x, h = h, stats = stats): Forcing the LD similarity to
    ## be smaller than or equal to 1

    ## Note: 133 merges with non increasing heights.

`snpClust` also handles inputs of class `base::matrix`:

    gmat <- as(geno, "matrix")
    fitM <- snpClust(geno, h, stats = "R.squared")

    ## Note: 133 merges with non increasing heights.

# References

[1] Ambroise C., Dehman A., Neuvial P., Rigaill G., and Vialaneix N. (2019). Adjacency-constrained hierarchical clustering of a band similarity matrix with application to genomics. Algorithms for Molecular Biology, 14, 22.

[2] Clayton D. (2015). snpStats: SnpMatrix and XSnpMatrix classes and methods. R package version 1.20.0

[3] Dehman A., Ambroise C., Neuvial P. (2015). Performance of a blockwise approach in variable selection using linkage disequilibrium information. BMC Bioinformatics, 16, 148.

[4] Randriamihamison N., Vialaneix N., and Neuvial P. (2021). Applicability and interpretability of Ward’s hierarchical agglomerative clustering with or without contiguity constraints. Journal of Classification, 38, 363–389.

# Session information

    sessionInfo()

    ## R version 4.5.1 (2025-06-13)
    ## Platform: x86_64-pc-linux-gnu
    ## Running under: Ubuntu 24.04.3 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.26.so;  LAPACK version 3.12.0
    ##
    ## locale:
    ##  [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C
    ##  [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=fr_FR.UTF-8
    ##  [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8
    ##  [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C
    ##  [9] LC_ADDRESS=C               LC_TELEPHONE=C
    ## [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C
    ##
    ## time zone: Europe/Paris
    ## tzcode source: system (glibc)
    ##
    ## attached base packages:
    ## [1] stats4    stats     graphics  grDevices utils     datasets  methods
    ## [8] base
    ##
    ## other attached packages:
    ##  [1] snpStats_1.46.0      Matrix_1.6-5         survival_3.8-3
    ##  [4] adjclust_0.6.11      HiTC_1.40.0          GenomicRanges_1.48.0
    ##  [7] GenomeInfoDb_1.32.4  IRanges_2.30.1       S4Vectors_0.34.0
    ## [10] BiocGenerics_0.42.0
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] SummarizedExperiment_1.26.1 gtable_0.3.6
    ##  [3] capushe_1.1.2               rjson_0.2.23
    ##  [5] xfun_0.51                   bslib_0.9.0
    ##  [7] ggplot2_3.5.2               Biobase_2.56.0
    ##  [9] lattice_0.22-6              vctrs_0.6.5
    ## [11] tools_4.5.1                 bitops_1.0-9
    ## [13] generics_0.1.4              parallel_4.5.1
    ## [15] tibble_3.2.1                pkgconfig_2.0.3
    ## [17] RColorBrewer_1.1-3          sparseMatrixStats_1.8.0
    ## [19] lifecycle_1.0.4             GenomeInfoDbData_1.2.8
    ## [21] compiler_4.5.1              farver_2.1.2
    ## [23] Rsamtools_2.12.0            Biostrings_2.64.1
    ## [25] codetools_0.2-20            htmltools_0.5.8.1
    ## [27] sass_0.4.9                  RCurl_1.98-1.16
    ## [29] yaml_2.3.10                 pillar_1.10.2
    ## [31] crayon_1.5.3                jquerylib_0.1.4
    ## [33] MASS_7.3-60.0.1             BiocParallel_1.30.4
    ## [35] DelayedArray_0.22.0         cachem_1.1.0
    ## [37] viridis_0.6.5               tidyselect_1.2.1
    ## [39] digest_0.6.37               dplyr_1.1.4
    ## [41] restfulr_0.0.15             splines_4.5.1
    ## [43] fastmap_1.2.0               grid_4.5.1
    ## [45] cli_3.6.5                   magrittr_2.0.3
    ## [47] XML_3.99-0.18               scales_1.4.0
    ## [49] rmarkdown_2.29              XVector_0.36.0
    ## [51] matrixStats_1.5.0           gridExtra_2.3
    ## [53] evaluate_1.0.5              knitr_1.49
    ## [55] BiocIO_1.6.0                viridisLite_0.4.2
    ## [57] rtracklayer_1.56.1          rlang_1.1.6
    ## [59] Rcpp_1.0.14                 dendextend_1.19.0
    ## [61] glue_1.8.0                  rstudioapi_0.17.1
    ## [63] jsonlite_2.0.0              R6_2.6.1
    ## [65] MatrixGenerics_1.8.1        GenomicAlignments_1.32.1
    ## [67] zlibbioc_1.42.0

1. We have dropped SNP rs2401075 because it produced a missing value due to the lack of genetic diversity in the considered sample.↩︎