本文首發(fā)于“生信補(bǔ)給站”,https://mp.weixin.qq.com/s/Gl6BChxSYbSHMo9oMpufPg
連鎖不平衡圖,用來可視化不同SNP之間的連鎖程度,前同事間俗稱“倒三角”圖。 本文使用自己的數(shù)據(jù),因為安裝R包后使用內(nèi)置數(shù)據(jù)集運(yùn)行出結(jié)果較容易,但是自己的數(shù)據(jù)就可能會有一些不大不小的“坑”,我替你們趟了。。。 一 載入R包 數(shù)據(jù)數(shù)據(jù)為內(nèi)置CEUData保存后,進(jìn)行了“細(xì)微”的處理(去掉SNP堿基之間的“/”),這種基因型文件很常見; library("LDheatmap") #讀入數(shù)據(jù) SNP <- read.csv("CEUSNP.csv",header = TRUE) pos <- read.csv("CEUDist.csv",header= TRUE) #查看數(shù)據(jù) 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|數(shù)據(jù)列的分分合合,一分多,多合一的separate和unite,可是沒有分隔符。。 經(jīng)高人指點(diǎn) ,使用替換的方式,解決方法很多。此處使用R-do包的函數(shù) 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
額 ,還是不行,同樣的報錯。檢索報錯,嘗試轉(zhuǎn)換數(shù)據(jù)格式。 2.3 堿基型轉(zhuǎn)為genotype使用genetics包的函數(shù)轉(zhuǎn)化 library("genetics") for(i in 1:ncol(SNPdata)){ SNPdata[,i]<-as.genotype(SNPdata[,i]) } LDheatmap(SNPdata, SNPpos,color = grey.colors(20))
額 ,,,終于可以了。。。 三 圖形調(diào)整,優(yōu)化3.1 調(diào)整顏色,更改標(biāo)題,標(biāo)示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調(diào)整SNP標(biāo)記名稱字體大小、顏色library(grid) grid.edit(gPath("ldheatmap", "geneMap","SNPnames"), gp = gpar(col="black",lwd = 1,cex=0.7)) 所謂的”倒三角圖“完成,haploview軟件也很好看,且有block,批量也許不太友好,見仁見智了!
|