Thanks a million!!!!
• 0 views
•
link
written
by
selplat21 •
Plotting LD heat map with plink.ld result
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)
• 2,234 views
•
link
written
by
selplat21 •
1 answer
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)
Log in to answer this question.