这里是佳奥！

ChIP-seq分析进入尾声。

在文件目录复制R Project文件，直接打卡R即可定位到该文件夹，不用翻文件了。

QQ截图20220813191310.png
我们现在需要的是peaks目录下的.bed文件。（软件安装部分详见https://www.jianshu.com/p/be9f7b52d042）

结果的注释用的是Y叔的Chipseeker包。

ChIPseeker的功能分为三类:

● 注释：提取peak附近最近的基因， 注释peak所在区域

● 比较：估计ChIP peak数据集中重叠部分的显著性；整合GEO数据集，以便于将当前结果和已知结果比较

● 可视化： peak的覆盖情况；TSS区域结合的peak的平均表达谱和热图；基因组注释；TSS距离；peak和基因的重叠。

1 初步查看

##查看peaks分布位置
 covplot(peak,weightCol = "V5")

QQ截图20220813205346.png

##查看peaks分布在染色体1和2的位置
covplot(peak, chr = c("chr1", "chr2"))

QQ截图20220814171717.png

##ChIP peaks结合TSS 区域的情况
promoter <- getPromoters(TxDb = txdb, upstream = 3000, downstream = 3000)##计算ChIP peaks结合在TSS区域的情况
tagMatrix <- getTagMatrix(peak, windows = promoter)##比对map到这些区域的peak，并生成tagMatrix
tagHeatmap(tagMatrix,xlim = c(-3000,3000), color = "red")

##查看这些peaks所在基因的启动子附近的分布情况，信号强度曲线图
plotAvgProf(tagMatrix,xlim = c(-3000,3000),
            xlab = "Genomic Region (5'->3')", ylab = "Read Count Frequency")

##peaks的注释
peakAnno <- annotatePeak(peak, tssRegion = c(-3000,3000),
                         TxDb = txdb, annoDb = "org.Hs.eg.db")
peakAnno_df <- as.data.frame(peakAnno)

sampleName=strsplit(bedPeaksFile,'\\.')[[1]][1]
print(samplename)

write.csv(peakAnno_df, paste0(sampleName,'_peakAnno_df.csv'))
DT::datatable(peakAnno_df,
             extensions=list(
             #dom='t'
             scrollx=TRUE
             fixedColumes=TRUE
             ))

##查看peaks长度分布，只统计长度1000bp以下peaks
peaksLength=abs(peakAnno_df$end-peakAnno_df$strat)
peaksLength=peaksLength[peaksLengt<1000]
hist(peakLength, breaks=50, col="lightblue", xlim=c(0,1000), xlab="peak length", main="Histogram of peak length")

2 完整注释代码

bedPeaksFile = 'H3K36me3_summits.bed'; 
bedPeaksFile

##这一段连续代码得到peakAnno_df数据框，即可进行各种分析
require(ChIPseeker)
require(TxDb.Mmusculus.UCSC.mm10.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
require(clusterProfiler) 
peak <- readPeakFile( bedPeaksFile )  
keepChr= !grepl('_',seqlevels(peak))##去除带_的染色体
seqlevels(peak, pruning.mode="coarse") <- seqlevels(peak)[keepChr]
peakAnno <- annotatePeak(peak, tssRegion=c(-3000, 3000), 
                         TxDb=txdb, annoDb="org.Mm.eg.db") 
peakAnno_df <- as.data.frame(peakAnno)

QQ截图20220813212102.png

获得了名为peakAnno_df矩阵，就可以进行各种分析，如GEO篇的通路分析https://www.jianshu.com/p/361b21a1ff10。

3 可视化

##Pie and Bar plot
plotAnnoBar(peakAnno)
vennpie(peakAnno)

QQ截图20220813213403.png QQ截图20220813213420.png

##可视化TSS区域的TF binding loci
plotDistToTSS(peakAnno,title="Distribution of transcription factor-binding loci\nrelative to TSS")

QQ截图20220813213611.png

install.packages('ggupset')
upsetplot(peakAnno)

QQ截图20220814165137.png

4 多个peak的比较

##多个peak set注释时，先构建list,然后用lapply.list(name1=bed_file1,name2=bed_file2)

peaks <- list(H2Aub1='H2Aub1_summits.bed',H3K36me3='H3K36me3_summits.bed',Ring1B='Ring1B_summits.bed')
promoter <- getPromoters(TxDb=txdb, upstream=2000, downstream=2000)
tagMatrixList <- lapply(peaks, getTagMatrix, windows=promoter)
plotAvgProf(tagMatrixList, xlim=c(-2000, 2000))
plotAvgProf(tagMatrixList, xlim=c(-2000, 2000), conf=0.95,resample=500, facet="row")
tagHeatmap(tagMatrixList, xlim=c(-2000, 2000), color=NULL)

QQ截图20220814165730.png QQ截图20220814165827.png

这个热图有点怪，，，

QQ截图20220814165902.png

更多参考代码：
https://github.com/jmzeng1314/NGS-pipeline/tree/master/CHIPseq
http://www.bio-info-trainee.com/2773.html

当然，方法不止一种，使用Chipseeker包的过程中，有些关联包的更新会导致不可用，这时候要更新到GitHub最新的版本才可。

下一步便是继续寻找motif。

我们下一篇再见！