Chip-seq call peak和注释

我们用的MACS2来进行call peak的。

我的样本是：input，CK，treatment

================首先肯定是比对得到bam文件=========

===================下面就是call peak了=============

macs2 callpeak -t treatment.sort.bam -c Input.sort.bam -f BAMPE -g 2.1e9 -n treatment -B -q 0.05

macs2 callpeak -t CK.sort.bam -c Input.sort.bam -f BAMPE -g 2.1e9 -n CK -B -q 0.05

其中：-g：基因组大小  -f：因为我是pair-end的所以选的是BAMPE

这是运行完之后的效果。

MACS2输出结果

==================然后利用ChIPseeker做注释===============

简单的格式转化一下，我简单的提取了一下excel的信息为bed格式：

peak.xls结果

为什么要转化，不用bed文件呢，因为我发现bed文件的位置信息不太对，里面好像显示的是中心位置信息

bed文件

所以自己从excel里面提取了一个bed格式的文件

输入bed文件

library("gplots")

library("preprocessCore")

library("ChIPseeker")

library("GenomicFeatures")

txdb <- makeTxDbFromGFF("Z.may.gene.simple.gff3",format="gff3")

data <- readPeakFile("CK.figure.txt")

promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)

tagMatrix <- getTagMatrix(data, windows=promoter)

peakAnno <- annotatePeak(data, tssRegion=c(-3000, 3000),TxDb=txdb)

write.table(peakAnno,"peak.annotation.txt",sep="\t")

jpeg("chromosome.distribution.jpeg",width = 1200, height = 1200,res = 144)

covplot(data,chrs=c("1", "2","3", "4","5", "6","7", "8","9", "10"))

dev.off()

jpeg("TSS.heatmap.jpeg",width = 1200, height = 1200,res = 144)

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

dev.off()

jpeg("TSS.average.jpeg",width = 1200, height = 1200,res = 144)

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

dev.off()

jpeg("class.Tiao.jpeg",width = 1200, height = 1200,res = 144)

plotAnnoBar(peakAnno)

dev.off()

jpeg("class.pai.jpeg",width = 1200, height = 1200,res = 144)

vennpie(peakAnno)

dev.off()

基本上要注释的结果都有了。

chromosome.distribution.jpeg