Plotting LD heat map with plink.ld result

Question

[selplat21]

Hello,

I have an output from plink as a plink.ld file that looks like this:

              CHR_A  BP_A SNP_A             CHR_B  BP_B SNP_B       R2
1 JANIGQ010000021.1 15082     . JANIGQ010000021.1 15155     . 0.896000
2 JANIGQ010000021.1 18234     . JANIGQ010000021.1 18262     . 0.936724
3 JANIGQ010000021.1 18234     . JANIGQ010000021.1 18276     . 0.883963
4 JANIGQ010000021.1 18234     . JANIGQ010000021.1 18297     . 0.938125
5 JANIGQ010000021.1 18262     . JANIGQ010000021.1 18276     . 1.000000
6 JANIGQ010000021.1 18262     . JANIGQ010000021.1 18297     . 1.000000

There are some old threads on this topic, but the software are out of date in those threads. Does anyone have any advice on how to format the LD file output from plink to plot an LD heat map.

This is not straightforward with LDheatmap or snp.plotter as the plink.ld requires quite a bit of reformatting and often does not result in a square matrix.

An example set of code for plotting in R would be incredibly helpful.

I'm interested in viewing a particular part of this chromosome, which should also help to reduce the memory requirements for plotting:

# Read the data
ld_data <- read.table("Chr5.ld", header = TRUE)

# Filter the data for the desired range
filtered_ld_data <- subset(ld_data, BP_A >= 48540000 & BP_A <= 48560000 & BP_B >= 48540000 & BP_B <= 48560000)

Answer 1

[bk11]

Here is the how you can plot LD Heatmap of SNPs from PLINK LD result in R using LDheatmap.

rm(list=ls())
# Load necessary libraries
library(LDheatmap)

# Read the LD file
ld_data <- read.table("plink.ld", header = TRUE)
head(ld_data)
CHR_A    BP_A      SNP_A CHR_B    BP_B      SNP_B          R2
1      1  846808  rs4475691     1  854250  rs7537756 8.33739e-01
2      1  846808  rs4475691     1  861808 rs13302982 1.06904e-03
3      1  846808  rs4475691     1  998395  rs7526076 2.89505e-04
4      1  846808  rs4475691     1 1040026  rs6671356 5.98311e-03
5      1  846808  rs4475691     1 1062638  rs9442373 3.61126e-04
6      1  846808  rs4475691     1 1110019 rs11260542 6.11007e-04
7      1  846808  rs4475691     1 1113121 rs12092254 6.74151e-04
8      1  846808  rs4475691     1 1130727 rs10907175 5.31309e-04
9      1  846808  rs4475691     1 1158277  rs3813199 2.24842e-05
10     1  854250  rs7537756     1  861808 rs13302982 4.98467e-03

#Selecting a list of SNP for LD Plot
snps=c("rs4475691","rs7537756","rs13302982","rs7526076","rs6671356","rs9442373","rs11260542")
ld_data_SelectSNPs=subset(ld_data, SNP_B %in% snps)

# Create a matrix of r^2 values for the SNPs
# Assuming your SNPs are ordered in the file and you want to use all of them
snp_names <- unique(c(ld_data_SelectSNPs$SNP_A, ld_data_SelectSNPs$SNP_B))
n <- length(snp_names)
ld_matrix <- matrix(NA, nrow = n, ncol = n)
rownames(ld_matrix) <- snp_names
colnames(ld_matrix) <- snp_names

# Fill the matrix with r^2 values
for (i in 1:nrow(ld_data_SelectSNPs)) {
  snp_a <- ld_data_SelectSNPs$SNP_A[i]
  snp_b <- ld_data_SelectSNPs$SNP_B[i]
  r2 <- ld_data_SelectSNPs$R2[i]

  ld_matrix[snp_a, snp_b] <- r2
  ld_matrix[snp_b, snp_a] <- r2
}

# Fill the diagonal with 1 (r^2 with itself)
diag(ld_matrix) <- 1

# Create a physical map (the position of the SNPs)
# Assuming the SNPs are in the order of their physical positions
physical_map <- ld_data_SelectSNPs[, c("BP_A", "BP_B")]
physical_map <- as.numeric(physical_map$BP_A)
names(physical_map) <- snp_names

# Plot the LD heatmap
LDheatmap(ld_matrix, genetic.distances=NULL, distances="physical",
          LDmeasure="r", title="Pairwise LD", add.map=TRUE, add.key=TRUE,
          geneMapLocation=0.15, geneMapLabelX=NULL, geneMapLabelY=NULL,
          SNP.name=snps, color=NULL, newpage=TRUE,
          name="ldheatmap", vp.name=NULL, pop=FALSE, flip=FALSE, text=FALSE)

[selplat21]

Thanks a million!!!!