本文首發於“生信補給站”,https://mp.weixin.qq.com/s/Gl6BChxSYbSHMo9oMpufPg
連鎖不平衡圖,用來可視化不同SNP之間的連鎖程度,前同事間俗稱“倒三角”圖。
本文使用自己的數據,因為安裝R包后使用內置數據集運行出結果較容易,但是自己的數據就可能會有一些不大不小的“坑”,我替你們趟了。。。
一 載入R包 數據
數據為內置CEUData保存后,進行了“細微”的處理(去掉SNP鹼基之間的“/”),這種基因型文件很常見;
library("LDheatmap")
#讀入數據
SNP <- read.csv("CEUSNP.csv",header = TRUE)
pos <- read.csv("CEUDist.csv",header= TRUE)
#查看數據
head(pos)
SNP[1:4,1:4]
二 繪制連鎖不平衡圖
2.1 直接繪制
SNPpos <- pos$x
LDheatmap(SNP, SNPpos,color = grey.colors(20))
Error in LDheatmap(SNP, SNPpos, color = grey.colors(20)) : column 1 is not a genotype object
額, 也許是因為沒有“/”的原因,加上試試?
2.2 鹼基型之間加“/“
怎么加呢? 首先想到 Tidyverse|數據列的分分合合,一分多,多合一的separate
和unite
,可是沒有分隔符。。
經高人指點 ,使用替換的方式,解決方法很多。此處使用R-do包的函數
library(do)
df <- na.omit(SNP)
#A,C,G ,T 替換為A/,C/,G/,T/
df1 = do::Replace(df,pattern = c("A:A/","C:C/","G:G/","T:T/"))
#去掉最后的/
SNPdata <- do::Trim(df1,"/")
SNPdata[1:4,1:4]
rs4615512 rs2283089 rs1894731 rs2283092
1 T/C C/C A/A T/T
2 T/C C/T A/A T/T
3 T/C C/C A/A T/T
5 T/C C/C A/A T/T
加上了,再次繪圖
LDheatmap(SNPdata, SNPpos,color = grey.colors(20))
Error in LDheatmap(SNPdata, SNPpos, color = grey.colors(20)) : column 1 is not a genotype object
額 ,還是不行,同樣的報錯。檢索報錯,嘗試轉換數據格式。
2.3 鹼基型轉為genotype
使用genetics包的函數轉化
library("genetics")
for(i in 1:ncol(SNPdata)){
SNPdata[,i]<-as.genotype(SNPdata[,i])
}
LDheatmap(SNPdata, SNPpos,color = grey.colors(20))
額 ,,,終於可以了。。。
三 圖形調整,優化
3.1 調整顏色,更改標題,標示SNP名稱
color.rgb <- colorRampPalette(rev(c("white","red")),space="rgb")
## 繪制連鎖不平衡圖
names <- c("rs1111183", "rs2237789", "rs2299531")
LDheatmap(SNPdata, SNPpos,
color=color.rgb(20),
title = "DEU Pairwise LD",
SNP.name=names,flip=TRUE)
3.2 使用grid調整SNP標記名稱字體大小、顏色
library(grid)
grid.edit(gPath("ldheatmap", "geneMap","SNPnames"),
gp = gpar(col="black",lwd = 1,cex=0.7))
所謂的”倒三角圖“完成,haploview軟件也很好看,且有block,批量也許不太友好,見仁見智了!