搜索
查看: 8188|回复: 4

[CHIP-seq] 使用R包DiffBind进行chip-seq差异峰的分析

[复制链接]

2

主题

9

帖子

235

积分

中级会员

Rank: 3Rank: 3

积分
235
发表于 2016-12-14 01:29:36 | 显示全部楼层 |阅读模式
DiffBind是一个分析chip-seq差异的峰的R包,可以从bioconduct获得。
source("http://bioconductor.org/biocLite.R")
biocLite("DiffBind")
分析前需要准备一个CSV文件。该R包的说明书里提供了一个tamoxifen处理乳腺癌细胞,对比耐药组和敏感组的转录因子Era的CHIP-seq的结合差异。下面的文件也是以这个为例子,但是因为例子的原始数据我没有找到,所以用的chip-seq数据是来自GEO的公共数据。
tamoxifen .CSV文件的内容见附件

bamReads 这列是写明比对chip-seq数据到参考基因组后的bam文件的名称。
bamControl一列 写的是input数据比对到参考基因组后的bam文件名称,
Peaks一列是peak文件,我是通过运行macs2产生的bed文件
这些对应名称的文件要放置在你设置的工作目录下面。我的工作目录是E:/R-testdata
一、        读入数据
library(DiffBind)
setwd("E:/R-testdata")
tamoxifen <- dba(sampleSheet="tamoxifen.csv")
tamoxifen
查看tamoxifen的内容,
4 Samples, 35720 sites in matrix (55407 total):
      ID Tissue Factor  Condition  Treatment Replicate Caller Intervals
1 BT4741  BT474     ER  Resistant Full-Media         1    raw     32756
2 BT4742  BT474     ER  Resistant Full-Media         2    raw     33067
3  MCF71   MCF7     ER Responsive Full-Media         1    raw     40982
4  MCF72   MCF7     ER Responsive Full-Media         2    raw     41593
Intervals一列就是对应的peak文件的peak数目。
二、计算reads
tamoxifen <- dba.count(tamoxifen, summits=250)
summits参数的意义是chip-seq峰的宽度的一半。
返回的内容如下:
Sample: SRR2147133.sorted.bam125
Sample: SRR2147134.sorted.bam125
Sample: SRR2147125.sorted.bam125
Sample: SRR2147126.sorted.bam125
Sample: SRR2147147.sorted.bam125
Sample: SRR2147148.sorted.bam125
Re-centering peaks...
Sample: SRR2147133.sorted.bam125
Sample: SRR2147134.sorted.bam125
Sample: SRR2147125.sorted.bam125
Sample: SRR2147126.sorted.bam125
Sample: SRR2147147.sorted.bam125
Sample: SRR2147148.sorted.bam125
Warning messages:
1: In dba.multicore.init(DBA$config) :
  Parallel execution unavailable: executing serially.
2: In dba.parallel(pv) : UNKNOWN PARALLEL PACKAGE: 0
这两个Warning messages的原因我还没弄清楚。
三、分析差异
根据condition一列的Resistant和Responsive进行分析差异。
将Resistant组和Responsive组都设置一个dba.mask,分析差异的时候需要用到。
ResistantMask <- dba.mask(tamoxifen,DBA_CONDITION, "Resistant")
ResponsiveMask <- dba.mask(tamoxifen,DBA_CONDITION, "Responsive")
分析差异需要安装DESeq2,
source("http://bioconductor.org/biocLite.R")
biocLite("DESeq2")
library(DESeq2)
分析差异:
tamoxifen <- dba.contrast(tamoxifen,ResistantMask,ResponsiveMask)
tamoxifen <- dba.analyze(tamoxifen)
tamoxifen.DB <- dba.report(tamoxifen)
tamoxifen.DB
把tamoxifen.DB输出到csv文件
write.csv(tamoxifen.DB,"diffbind.csv")
查看tamoxifen.DB的内容 见附件

发现所有的peak的width一列宽度都设在了501,这个值是和tamoxifen <- dba.count(tamoxifen, summits=250)中的参数summits=250有关。
DiffBind包的dba.plotHeatmap函数画的差异峰的热图:
pdf("diffbind.pdf")

dba.plotHeatmap(tamoxifen, contrast=1, correlations=FALSE)
dev.off()


图片见附件
,输出的文件diffbind.csv里,峰的位置是以染色体,起始位置,终止位置来保存的。群主的博客里有介绍如何使用Annovar可以进行后续的注释。我还不会用,所以用了网页版的工具wANNOVAR。
小菜鸟第一次发帖,如果有错误的地方,恳请各位前辈指出。

本帖子中包含更多资源

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

x



上一篇:我bwa index 到一半了,突然被中断,
下一篇:表达芯片分析[第一讲:GEO,表达芯片与R]
回复

使用道具 举报

0

主题

1

帖子

43

积分

新手上路

Rank: 1

积分
43
发表于 2017-8-3 18:19:40 | 显示全部楼层
你好,请问您说的群主的博客进行了后续分析,群主的博客我可以去哪里找到?谢谢!
回复 支持 反对

使用道具 举报

10

主题

52

帖子

559

积分

版主

Rank: 7Rank: 7Rank: 7

积分
559
QQ
发表于 2017-8-4 16:34:55 | 显示全部楼层
闫凯robert 发表于 2017-8-3 18:19
你好,请问您说的群主的博客进行了后续分析,群主的博客我可以去哪里找到?谢谢!
...

可以点击论坛上面目录“认识管理员”,第一个“生信菜鸟团”就是。
回复 支持 反对

使用道具 举报

1

主题

2

帖子

72

积分

注册会员

Rank: 2

积分
72
发表于 2017-9-21 11:22:13 | 显示全部楼层
你好,请问能将输入的数据和文件发一份给我吗?我在试这个软件过程中遇到了些问题
回复 支持 反对

使用道具 举报

5

主题

10

帖子

166

积分

注册会员

Rank: 2

积分
166
发表于 2017-11-14 21:02:31 | 显示全部楼层
这个包好像只能处理有重复的样本,请问有没有其他好的方法处理没有重复的两组组蛋白修饰差异区域。
回复 支持 反对

使用道具 举报

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

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-10-20 07:26 , Processed in 0.031941 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.