• 主页
  • 课程

    关于课程

    • 课程归档
    • 成为一名讲师
    • 讲师信息
    同等学历教学

    同等学历教学

    免费
    阅读更多
  • 特色
    • 展示
    • 关于我们
    • 问答
  • 事件
  • 个性化
  • 博客
  • 联系
  • 站点资源
    有任何问题吗?
    (00) 123 456 789
    weinfoadmin@weinformatics.cn
    注册登录
    恒诺新知
    • 主页
    • 课程

      关于课程

      • 课程归档
      • 成为一名讲师
      • 讲师信息
      同等学历教学

      同等学历教学

      免费
      阅读更多
    • 特色
      • 展示
      • 关于我们
      • 问答
    • 事件
    • 个性化
    • 博客
    • 联系
    • 站点资源

      未分类

      • 首页
      • 博客
      • 未分类
      • R-ChIP数据可视化

      R-ChIP数据可视化

      • 发布者 weinfoauthor
      • 分类 未分类
      • 日期 2019年10月21日
      • 评论 0评论

      2017眼看要结束,立下写《CS0: ChIPseq从入门到放弃》的flag还没完成,当时ChIPseeker是33个引用,现在已经80了,时间过得好快。

      最近放羊的Jimmy给我发来了一个截屏:

      img

      说了一篇新文章,大段在称赞ChIPseeker:

      img

      这篇文章的作者以前就跟我要ChIPseeker的文章,说是上课要给学生读,所以我要感谢他,上课和文章都给我打了广告!

      img

      最近小密圈有小伙伴问怎么对lncRNAs按基因位置进行分类,我让他去看这个ChIPseeker系列文章,不要被ChIPseeker这个名字给蒙蔽了,我当年写这个包的时候,我只是想做ChIPseq,我也没想到自己写的是个通用的包,其实应该叫genomic coordinates annotation, visualization and comparison。ChIPseeker不单用于ChIPseq,所有genomic coordination通杀。在《CS4:关于ChIPseq注释的几个问题》里我就谈到了有人做DNA breakpoints的断点位置注释。在biostars上我也展示过注释lincRNA (https://www.biostars.org/p/215069/)。另外CS4里也提到了背景注释可以用GRanges对象,这保证了各种需求都可以满足。

      如果你用ChIPseeker做注释,那么出上面的饼图,太容易,ChIPseeker只会做得更好。

      Peak可视化

      首先在分析之前,我们拿到BED文件,就可以先看看peak在整个染色体上的分布,这个在《CS2: BED文件 》一文里有展示,这里不在重复。

      另外我们可以限定在一个固定的窗口里,然后把这些peak全部align起来,看看在某个窗口上的结合谱图。比如说启动子区域,使用转录起始位点,然后指定上下游,使用getPromoters函数,就可以为我们准备好窗口。再使用getTagMatrix函数,把peak比对到这个窗口,并生成矩阵,供我们分析、可视化。

      library(ChIPseeker)
      files <- getSampleFiles()
      peak <- readPeakFile(files[[4]])
      promoter <- getPromoters(TxDb=txdb, 
          upstream=3000, downstream=3000)
      tagMatrix <- getTagMatrix(peak, windows=promoter)
      

      这时候,我们可以用tagHeatmap来画热图,颜色随意指定:

      tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red")
      

      img

      另外你可以直接从文件出发,ChIPseeker提供了peakHeatmap函数,你给文件名,就可以直接出图,两个函数都支持多种数据同时展示,方便比较,这里可以看出前三个样品结合位点不在转录起始区域,而后两个则是。

      require(TxDb.Hsapiens.UCSC.hg19.knownGene)
      txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
      peakHeatmap(files, weightCol="V5", TxDb=txdb, 
          upstream=3000, downstream=3000, 
          color=rainbow(length(files)))
      

      img

      另外我们还可以以谱图的形式来展示结合的强度:

      plotAvgProf(tagMatrix, xlim=c(-3000, 3000),
                  xlab="Genomic Region (5'->3')", 
                  ylab = "Read Count Frequency")
                  
      

      img

      上面这个图可以非常清晰地看到,这个结合位点是集中在转录起始位点的。另外ChIPseeker也提供了plotAvgProf2函数,支持以文件为输入一步做图。下面这个代码,会产生上面一样的图。

      plotAvgProf2(files[[4]], TxDb=txdb, 
      			 upstream=3000, downstream=3000,
                   xlab="Genomic Region (5'->3')", 
                   ylab = "Read Count Frequency")
      

      另外这个函数还支持使用bootstrap估计置信区间:

      plotAvgProf(tagMatrix, xlim=c(-3000, 3000), 
                  conf = 0.95, resample = 1000)
      

      img

      这个功能由微博@Puriney贡献,你们可以去微博勾搭调戏,是个颜值很高的鲜肉。

      同样这个函数也支持多个数据的比较:

      tagMatrixList <- lapply(files, getTagMatrix, 
                              windows=promoter)
      plotAvgProf(tagMatrixList, xlim=c(-3000, 3000))
      

      img

      我们可以对其进行分面展示,当然也支持估计置信区间:

      plotAvgProf(tagMatrixList, xlim=c(-3000, 3000), 
          conf=0.95,resample=500, facet="row")
      

      img

      这和前面的热图是一致的,前三个样品显然不是调控转录的。

      这里我展示的都是启动子窗口,然而ChIPseeker支持的并不止是启动子,不单单是看转录调控,比如说你是研究可变剪切的,那么蛋白结合到外显子/内含子的起始位置,显然就更感兴趣了。所以除了getPromoters函数之外,ChIPseeker还提供了getBioRegion函数,你可以指定’exon’或’intron’,它会像getPromoters一样为你准备好数据,即使这不是你感兴趣的,你也依然可以用它来准备窗口,然后可视化看一下,或许你能意外地发现你研究的蛋白竟然在外显子/内含子起始位置上有很强的结合谱,做为前期data exploration,「妹妹你大胆地往前走啊」,尽管啥都试试!

      注释信息可视化

      我在《CS3: peak注释》和《CS4:关于ChIPseq注释的几个问题》两篇文章里,重点讲了注释,注释后的信息当然也需要我们可视化展示,方便解析结果。

      peakAnno <- annotatePeak(files[[4]], tssRegion=c(-3000, 3000),
                               TxDb=txdb, annoDb="org.Hs.eg.db")
      plotAnnoPie(peakAnno)
      

      img

      首先是饼图,这也是小密圈的小伙伴问的,peak落在了什么地方,启动子?外显子?内含子?还是基因间区?饼图给出了比例,当然一个peak所在的位置可能是一个基因的外显子而同时又是另一个基因的内含子,所以annotatePeak有个参数genomicAnnotationPriority让你指定优先顺序,默认顺序如下:

      • Promoter
      • 5’ UTR
      • 3’ UTR
      • Exon
      • Intron
      • Downstream
      • Intergenic

      这里有下游但没有上游,因为Promter我们定义为转录起始位点的上下游区域,所以它包括了上游,而下游就是在基因间区里,而我又把它单独拿出来,因为这是紧接着基因的下游,和远程的间区分开,通常默认参数就可以。

      有饼图,当然有柱状图,这是躺平任调戏的柱子:

      plotAnnoBar(peakAnno)
      

      img

      正如上面说的,一个结合位点可以是一个基因的Exon和另一个基因的Intron,除了按照优先级别给出唯一的注释之外,ChIPseeker其实对所有的注释都有记录,我当时就想怎么样把这些信息给展示出来,当时在笔记本上涂鸦了这样的图:

      img

      这个图我称之为vennpie取维恩图venn plot画交集和饼图pie plot画分布的结合。

      vennpie(peakAnno)
      

      img

      这个图讲真我很喜欢,但不满意,因为无法展示全部的信息。所以当我看到有UpSetR这个包的时候,我立刻就把它应用于展示注释信息的交集:

      upsetplot(peakAnno)
      

      img

      这个可以很好的把信息全部展示出来,然而第一不够vennpie直观,第二右上角空白太多,所以我就想把vennpie强插上去,最终通过参数实现这个效果:

      upsetplot(peakAnno, vennpie=TRUE)
      

      img

      ChIPseeker还注释了最近的基因,peak离最近基因的距离分布是什么样子的?ChIPseeker提供了plotDistToTSS函数来画这个分布:

      plotDistToTSS(peakAnno,
          title="Distribution of transcription factor-binding loci\nrelative to TSS")
      

      img

      plotAnnoBar和plotDistToTSS这两个柱状图都支持多个数据同时展示,方便比较,比如:

      peakAnnoList <- lapply(files, annotatePeak, TxDb=txdb,
                      tssRegion=c(-3000, 3000), verbose=FALSE)
      plotAnnoBar(peakAnnoList)
      

      img

      plotDistToTSS(peakAnnoList)
      

      img

      ChIPseeker还提供了一个vennplot函数,比如我想看注释的最近基因在不同样本中的overlap:

      genes <- lapply(peakAnnoList, function(i) 
          as.data.frame(i)$geneId)
      vennplot(genes[2:4], by='Vennerable')
      

      img

      如果你们有画图的需求,尽管来找我合作,有偿做图也可以,请参考「生信技能树」公众号的有偿专区。

      ChIPseeker系列

      • CS0: ChIPseq从入门到放弃
      • CS1: ChIPseq简介
      • CS2: BED文件
      • CS3: peak注释
      • CS4:关于ChIPseq注释的几个问题
      • CS5: 吃着火锅,唱着歌,还把分析给做了

      请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!

      • 分享:
      weinfoauthor
      weinfoauthor

      1233

      上一篇文章

      R-ChIPseq-GEO数据挖掘
      2019年10月21日

      下一篇文章

      R-Genomic coordination的富集性分析
      2019年10月21日

      你可能也喜欢

      2-1675088548
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      30 1月, 2023
      9-1675131201
      如何快速批量修改 Git 提交记录中的用户信息
      26 1月, 2023
      8-1678501786
      肿瘤细胞通过改变CD8+ T细胞中的丙酮酸利用和琥珀酸信号来调控抗肿瘤免疫应答。
      7 12月, 2022

      留言 取消回复

      要发表评论,您必须先登录。

      搜索

      分类

      • R语言
      • TCGA数据挖掘
      • 单细胞RNA-seq测序
      • 在线会议直播预告与回放
      • 数据分析那些事儿分类
      • 未分类
      • 生信星球
      • 老俊俊的生信笔记

      投稿培训

      免费

      alphafold2培训

      免费

      群晖配置培训

      免费

      最新博文

      Nature | 单细胞技术揭示衰老细胞与肌肉再生
      301月2023
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      301月2023
      如何快速批量修改 Git 提交记录中的用户信息
      261月2023
      logo-eduma-the-best-lms-wordpress-theme

      (00) 123 456 789

      weinfoadmin@weinformatics.cn

      恒诺新知

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      链接

      • 课程
      • 事件
      • 展示
      • 问答

      支持

      • 文档
      • 论坛
      • 语言包
      • 发行状态

      推荐

      • iHub汉语代码托管
      • iLAB耗材管理
      • WooCommerce
      • 丁香园论坛

      weinformatics 即 恒诺新知。ICP备案号:粤ICP备19129767号

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      要成为一名讲师吗?

      加入数以千计的演讲者获得100%课时费!

      现在开始

      用你的站点账户登录

      忘记密码?

      还不是会员? 现在注册

      注册新帐户

      已经拥有注册账户? 现在登录

      close
      会员购买 你还没有登录,请先登录
      • ¥99 VIP-1个月
      • ¥199 VIP-半年
      • ¥299 VIP-1年
      在线支付 激活码

      立即支付
      支付宝
      微信支付
      请使用 支付宝 或 微信 扫码支付
      登录
      注册|忘记密码?