OS论坛专栏
叮咚。“OS论坛专栏”上线了,我们将定期挑选OS论坛(网址:www.omicshare.com/forum)的优质内容分享给大家,如代码贴、实用工具分享、软件教程、文献解读等等,给论坛优质内容更多曝光机会,欢迎大家前往OS论坛投稿~,今天推荐论坛入驻网友“生信补给站”的一篇帖子(原帖链接:https://www.omicshare.com/forum/thread-6269-1-1.html),跟大家分享连锁不平衡图的绘制。
连锁不平衡图,用来可视化不同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 object
使用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,批量也许不太友好,见仁见智了!
转自:基迪奥生物
- 本文固定链接: https://maimengkong.com/image/1079.html
- 转载请注明: : 萌小白 2022年6月30日 于 卖萌控的博客 发表
- 百度已收录