查看: 511|回复: 6

[R语言] LDheatmap|SNP连锁不平衡图(LD)可视化,自己数据可实现版!

[复制链接]
  • TA的每日心情

    2016.7.20 14:55
  • 签到天数: 1 天

    连续签到: 1 天

    [LV.1]初来乍到

    版主

    Rank: 10Rank: 10Rank: 10

    主题
    29
    奥币
    876
    积分
    390
    注册时间
    2016.3.21
    在线时间
    89 小时

    活跃会员突出贡献


    发表于 2020.5.9 10:13:01 | 显示全部楼层 |阅读模式
    本帖最后由 MJPS 于 2020.5.9 10:45 编辑
    本文首发于“生信补给站”https://mp.weixin.qq.com/s/Gl6BChxSYbSHMo9oMpufPg

             连锁不平衡图,用来可视化不同SNP之间的连锁程度,前同事间俗称“倒三角”图,。本文使用自己的数据,因为安装R包后使用内置数据集运行出结果较容易,但是自己的数据就可能会有一些不大不小的“坑”,我替你们趟了。。。
    一 载入R包 数据数据为内置CEUData保存后,进行了“细微”的处理(去掉SNP碱基之间的“/”),因为这种基因型形式文件很常见;
    [AppleScript] 纯文本查看 复制代码
    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 直接绘制
    [AppleScript] 纯文本查看 复制代码
    SNPpos <- pos$x
    LDheatmap(SNP, SNPpos,color = grey.colors(20))
    [AppleScript] 纯文本查看 复制代码
    Error in LDheatmap(SNP, SNPpos, color = grey.colors(20)) : 
      column 1 is not a genotype object

    额, 也许是因为没有“/”的原因,加上试试?2.2 碱基型之间加“/“怎么加呢?
    首先想到 Tidyverse|数据列的分分合合,一分多,多合一的separate和unite,可是没有分隔符。。
    经高人指点 ,使用替换的方式,解决方法很多。
    此处使用R-do包的函数
    [AppleScript] 纯文本查看 复制代码
    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
    加上了,再次绘图
    [AppleScript] 纯文本查看 复制代码
    LDheatmap(SNPdata, SNPpos,color = grey.colors(20))

    [AppleScript] 纯文本查看 复制代码
    Error in LDheatmap(SNPdata, SNPpos, color = grey.colors(20)) : 
      column 1 is not a genotype object

    额 ,还是不行,同样的报错。检索报错,尝试转换数据格式。
    2.3 碱基型转为genotype object
    使用genetics包的函数转化
    [AppleScript] 纯文本查看 复制代码
    library("genetics")
    for(i in 1:ncol(SNPdata)){
      SNPdata[,i]<-as.genotype(SNPdata[,i])
    }
    LDheatmap(SNPdata, SNPpos,color = grey.colors(20))


    额  ,,,终于可以了。。。

    三 图形调整,优化
    3.1 调整颜色,更改标题,标示SNP名称
    [AppleScript] 纯文本查看 复制代码
    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标记名称字体大小、颜色
    [AppleScript] 纯文本查看 复制代码
    library(grid)
    grid.edit(gPath("ldheatmap", "geneMap","SNPnames"), 
        gp = gpar(col="black",lwd = 1,cex=0.7))


    所谓的”倒三角图“完成,haploview软件也很好看,且有block,批量也许不太友好,见仁见智了!

                           



    本帖子中包含更多资源

    您需要 登录 才可以下载或查看,没有帐号?立即注册

    x

    评分

    参与人数 1奥币 +10 收起 理由
    小圆 + 10 楼主V5!

    查看全部评分

    回复

    使用道具 举报

    该用户从未签到

    草履虫

    Rank: 2

    主题
    0
    奥币
    22
    积分
    2
    注册时间
    2020.5.10
    在线时间
    0 小时

    发表于 2020.5.10 04:07:15 | 显示全部楼层
    12345678910
    回复 支持 反对

    使用道具 举报

  • TA的每日心情

    2019.12.9 11:32
  • 签到天数: 3 天

    连续签到: 1 天

    [LV.2]偶尔看看I

    钵水母

    Rank: 3Rank: 3

    主题
    0
    奥币
    24
    积分
    36
    注册时间
    2019.5.12
    在线时间
    7 小时

    发表于 2020.5.11 09:07:13 | 显示全部楼层
    1234567890
    回复 支持 反对

    使用道具 举报

  • TA的每日心情

    2020.5.19 17:03
  • 签到天数: 11 天

    连续签到: 4 天

    [LV.3]偶尔看看II

    钵水母

    Rank: 3Rank: 3

    主题
    0
    奥币
    25
    积分
    47
    注册时间
    2018.10.20
    在线时间
    15 小时

    最佳新人


    发表于 2020.5.17 10:50:59 | 显示全部楼层
    zhichi yixa
    尽快发士大夫
    回复 支持 反对

    使用道具 举报

  • TA的每日心情

    2016.7.20 14:55
  • 签到天数: 1 天

    连续签到: 1 天

    [LV.1]初来乍到

    版主

    Rank: 10Rank: 10Rank: 10

    主题
    29
    奥币
    876
    积分
    390
    注册时间
    2016.3.21
    在线时间
    89 小时

    活跃会员突出贡献


     楼主| 发表于 2020.5.18 16:39:40 | 显示全部楼层
    回复 支持 反对

    使用道具 举报

    该用户从未签到

    钵水母

    Rank: 3Rank: 3

    主题
    0
    奥币
    2
    积分
    11
    注册时间
    2019.2.13
    在线时间
    1 小时

    发表于 7 天前 | 显示全部楼层
    nice
    回复

    使用道具 举报

  • TA的每日心情
    yes!
    前天 08:46
  • 签到天数: 412 天

    连续签到: 1 天

    [LV.9]以坛为家II

    帝王蝶

    Rank: 4

    主题
    1
    奥币
    1705
    积分
    420
    注册时间
    2016.9.4
    在线时间
    106 小时

    发表于 5 天前 | 显示全部楼层
    学习了!。。。。
    回复

    使用道具 举报

    您需要登录后才可以回帖 登录 | 立即注册

    本版积分规则

    快速回复 返回顶部 返回列表