这里是佳奥！
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。
我们下一篇再见！
