学习linux常用命令，并且实战ChIP-seq已经一个月了。现在需要用R包ChIPseeker进行注释，发现R已经忘了大半。补课R语言。。对peak的注释分为两个部分——结构注释和功能注释
结构注释：会将peak所落在基因组上的区域结构注释出来，比如说启动子区域，UTR区域，内含子区域等等。同时，也会将与peak最临近的基因注释出来，非常好用。
功能注释：其实并不是对peak进行注释的，而是对peak最临近的基因注释的，因此这部分就和传统的GO/KEGG等分析如出一辙啦。
一、安装R包ChIPseeker，各种报错及解决
按照教程在console输入：BiocManager::install("ChIPseeker")
报错：Error in loadNamespace(x) : there is no package called ‘BiocManager’
于是输入：install.packages("BiocManager") 先安装BiocManager。
再输入BiocManager::install("ChIPseeker")
运行了一会儿出现如下报错：
Warning messages:
1: In install.packages(...) :
installation of package ‘DO.db’ had non-zero exit status
2: In install.packages(...) :
installation of package ‘GO.db’ had non-zero exit status
3: In install.packages(...) :
installation of package ‘igraph’ had non-zero exit status
4: In install.packages(...) :
installation of package ‘gtools’ had non-zero exit status
还需要安装一个GO和DO包,和其他包。
BiocManager::install("GO.db")
BiocManager::install("DO.db")
BiocManager::install("topGO")
BiocManager::install("GSEABase")
BiocManager::install("Rgraphviz")
install.packages("igraph")
install.packages("gtools")
在安装igraph和gtools时候，出现提示：Do you want to install from sources the packages which need compilation?
【解决办法】
选择“NO”，而不是“YES”，便可正确安装igraph和gtools。
二、注释库
BiocManager::install("org.Hs.eg.db")
BiocManager::install("org.Mm.eg.db")
BiocManager::install("TxDb.Mmusculus.UCSC.mm10.knownGene")
BiocManager::install("TxDb.Hsapiens.UCSC.hg38.knownGene")
这些是基因组注释库，最常用的人和小鼠的。
当然也可以用gff/gtf构建一个TxDb（注释库），略。
三、bed数据
把数据从服务器下载到Mac电脑上。一定要先连上VPN，在本电脑的Terminal中进行操作。
scp -r zds209@222.28.163.113:~/ChIP-seqtest/bed/ /Users/meraner/Desktop
-r是整个文件夹里的文件都下载到bed里面。速度能到3.7M/s。
四、使用ChIPseeker
参考教程
http://events.jianshu.io/p/9aa719faa4b5
//www.greatytc.com/p/c76e83e6fa57
https://www.modb.pro/db/401713
http://www.bioconductor.org/packages/devel/bioc/vignettes/ChIPseeker/inst/doc/ChIPseeker.html #官网
https://blog.csdn.net/qq_39306047/article/details/106601095
1. 先看看整体分布吧。
library(ChIPseeker) #加载ChIPseeker包
library(ggplot2) #加载ggplot2包
library(org.Mm.eg.db) #加载参考基因组数据
library(TxDb.Mmusculus.UCSC.mm10.knownGene) #加载参考基因组注释数据包
library(clusterProfiler) #加载clusterProfiler包
library(GenomicFeatures) #加载GenomicFeatures包
eithtWG16_1<-readPeakFile("8WG16-1.sorted.bam_peaks.narrowPeak") #读取narrowPeak数据
covplot(eithtWG16_1,weightCol=5) #可视化一下数据全景
这样的图片，指示的是在每条染色体上的富集程度
还能进行多样本，及制定区域的展示。
peak=GenomicRanges::GRangesList(WG1=readPeakFile("8WG16-1.sorted.bam_peaks.narrowPeak"),WG2=readPeakFile("8WG16-2.sorted.bam_peaks.narrowPeak"))#（先指定的在右侧，这个比较奇怪，所以以后还是要改成2号在前) 使用ChIPseeker包中的readPeakFile函数将bed文件读入到R，并存储为GRanges对象
covplot(peak, weightCol="V5") + facet_grid(chr ~ .id)#画图多样本展示每条染色体的富集程度，也可以理解为coverage plot 展示峰在染色体上分布的图。
covplot(peak, weightCol="V5", chrs=c("chr17", "chr18"), xlim=c(4e7, 5e7)) + facet_grid(chr ~ .id)#截取某些片段进行展示。
2.在启动子区(TSS)附近的 ChIP peaks 的可视化
（注：这个也可以用deeptools在linux服务器里面去跑。day26都在学deeptools。）
peak <- readPeakFile("8WG16-1.sorted.bam_summits.bed")#读取bed文件
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene #指定txdb
promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000) #getPromoters函数定义启动子区域的范围
tagMatrix <- getTagMatrix(peak, windows=promoter) # getTagMatrix函数将Peaks匹配到启动子区域，然后计算出tagMatrix用于画图
tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red") #热图
plotAvgProf(tagMatrix, xlim=c(-3000, 3000),
xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency") #折线图
plotAvgProf(tagMatrix, xlim=c(-3000, 3000),conf = 0.95, resample = 1000,
xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency") #折线图上加置信区间
多个样本，在启动子区(TSS)附近的 ChIP peaks 的可视化
peak1 <- "/Users/meraner/Desktop/bed/8WG16-1.sorted.bam_summits.bed"
peak2 <- "/Users/meraner/Desktop/bed/8WG16-1.sorted.bam_summits.bed"
files <-list(peak1=peak1,peak2=peak2)#把两个bed文件指定给files, 就是一次读入多个summits文件，使用list存储，然后使用lapply注释。
tagMatrixList <- lapply(files, getTagMatrix, windows=promoter)
plotAvgProf(tagMatrixList, xlim=c(-3000, 3000),
conf=0.95,resample=500, facet="row") # 添加置信区间并分面
3. 结构注释：
library(ChIPseeker)
library(ggplot2)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
library(org.Mm.eg.db)
peak <- readPeakFile("8WG16-1.sorted.bam_summits.bed")#读取bed文件
peakAnno <- annotatePeak(peak,TxDb = TxDb.Mmusculus.UCSC.mm10.knownGene, annoDb ="org.Mm.eg.db" ,tssRegion=c(-3000, 3000)) # 进行结构注释
注释完，进行可视化，多种图可供选择
plotAnnoBar(peakAnno) #柱状图
plotDistToTSS(peakAnno, title="Distribution of transcription factor-binding loci \n relative to TSS") #展示转录因子(TF)与基因启动子区的相对位置vennpie(peakAnno)
plotAnnoPie(peakAnno) #饼图
上面这几个图都是能看出peak的分布结构的
还有另外软件可以对图进行优化。安装下面两个包。
BiocManager::install("ggupset")
install.packages("ggimage")
library(ggupset)
library(ggimage)
然后可以运行
upsetplot(peakAnno)
upsetplot(peakAnno, vennpie=TRUE)
4.多个样本结构注释
library(ChIPseeker) #加载ChIPseeker包
library(ggplot2) #加载ggplot2包
library(org.Mm.eg.db) #加载参考基因组数据
library(TxDb.Mmusculus.UCSC.mm10.knownGene) #加载参考基因组注释数据包
library(clusterProfiler) #加载clusterProfiler包
library(GenomicFeatures) #加载GenomicFeatures包
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
setwd("/Users/meraner/Desktop/bed") #制定当前工作目录，为数据存储目录
peak1 <- "/Users/meraner/Desktop/bed/8WG16-1.sorted.bam_summits.bed"
peak2 <- "/Users/meraner/Desktop/bed/8WG16-1.sorted.bam_summits.bed"
files <-list(peak1=peak1,peak2=peak2)#把两个bed文件指定给files, 就是一次读入多个summits文件，使用list存储，然后使用lapply注释。
peakAnnoList <- lapply(files, annotatePeak, TxDb=TxDb.Mmusculus.UCSC.mm10.knownGene,
tssRegion=c(-3000, 3000),
verbose=FALSE) #注释
plotAnnoBar(peakAnnoList) #画图bar
plotDistToTSS(peakAnnoList) #画图距离
#不能做多样本的饼图和韦恩图。
在看完数据全景之后，就会想知道这些peaks和什么类型的基因有关。
5.功能注释：
基因注释+富集分析
利用ChIPseeker的seq2gene 将peak的位置与所有的基因关联起来【包括 host gene (exon/intron), promoter region and flanking gene from intergenomic region】，然后用clusterProfiler拿这些基因跑ORA，做富集
下面是根据写好的R脚本，执行顺利。
6. 单样本ChIPseeker 脚本
#工作目录，样本种属来源需要根据情况调整
setwd("/Users/meraner/Desktop/bed") #制定当前工作目录，为数据存储目录
library(ChIPseeker) #加载ChIPseeker包
library(ggplot2) #加载ggplot2包
library(org.Mm.eg.db) #加载参考基因组数据
library(TxDb.Mmusculus.UCSC.mm10.knownGene) #加载参考基因组注释数据包
library(clusterProfiler) #加载clusterProfiler包
library(GenomicFeatures) #加载GenomicFeatures包
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
#peak的整体分布
peak<-readPeakFile("8WG16-1.sorted.bam_peaks.narrowPeak") #读取narrowPeak数据
covplot(peak, weightCol=5) #可视化一下数据全景
covplot(peak, weightCol=5, chrs=c("chr17", "chr18"), xlim=c(4e7, 5e7)) #截取某些染色体，某些片段进行展示
#单一样本，在启动子区(TSS)附近的 ChIP peaks 的可视化
peak <- readPeakFile("8WG16-1.sorted.bam_summits.bed") #读取bed文件
promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000) #getPromoters函数定义启动子区域的范围
tagMatrix <- getTagMatrix(peak, windows=promoter) # getTagMatrix函数将Peaks匹配到启动子区域，然后计算出tagMatrix用于画图
tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red") #一个样本在TSS附近结合区的热图
plotAvgProf(tagMatrix, xlim=c(-3000, 3000),
xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency") #折线图，
plotAvgProf(tagMatrix, xlim=c(-3000, 3000),conf = 0.95, resample = 1000,
xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency") #前面的折线图上加置信区间
#注释，结构可视化
peakAnno <- annotatePeak(peak,TxDb = txdb, annoDb ="org.Mm.eg.db" ,tssRegion=c(-3000, 3000)) # 进行结构注释
plotAnnoBar(peakAnno) #柱状图
plotDistToTSS(peakAnno, title="Distribution of transcription factor-binding loci \n relative to TSS") #展示转录因子(TF)与基因启动子区的相对位置vennpie(peakAnno)
plotAnnoPie(peakAnno) #饼图
library(ggupset)
library(ggimage)
upsetplot(peakAnno)
upsetplot(peakAnno, vennpie=TRUE)
#功能可视化-基因注释 + 富集分析
library(DOSE)
library(GO.db)
library(topGO)
library(GSEABase)
library(clusterProfiler)
require(clusterProfiler)
gene <- seq2gene(peak, tssRegion = c(-1000, 1000), flankDistance = 3000, TxDb=txdb)
enk <- enrichKEGG(gene=gene,
organism = 'mmu',
pvalueCutoff = 0.05,
pAdjustMethod = "BH",
qvalueCutoff = 0.05) #KEGG通路富集
dotplot(enk)#可视化KEGG富集
browseKEGG(enk, 'mmu04110') # 会打开浏览器调到KEGG数据库的这条通路上，并且富集的基因会以不同的颜色标出
engo_MF <- enrichGO(gene = gene,
OrgDb = org.Mm.eg.db,
ont = "MF",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05,
readable = TRUE) #GO富集分析MF
dotplot(engo_MF) #可视化GO富集
barplot(engo_MF, showCategory=15) #可视化GO富集-bar图
plotGOgraph(engo_MF) #对Go通路树状图
engo_BP <- enrichGO(gene = gene,
OrgDb = org.Mm.eg.db,
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05,
readable = TRUE) #GO富集分析BP
dotplot(engo_BP) #可视化GO富集-点图
barplot(engo_BP, showCategory=15) #可视化GO富集-bar图
enrichMap(engo_BP) #这个命令有问题,报错could not find function "enrichMap"？没搞懂，也许是R版本问题。
plotGOgraph(engo_BP) #对Go通路树状图
#输出注释信息
peakAnnotable<-as.data.frame(peakAnno)
write.table(peakAnnotable,file="result_genes",sep="\t",quote=F) #把peak相关基因导出表格
enktable<-as.data.frame(enk)
write.table(enktable,file="result_kegg.txt",sep="\t",quote=F)#把KEGG富集结果导出表格
MFtable<-as.data.frame(engo_MF)
write.table(MFtable,file="result_GO_MF.txt",sep="\t",quote=F)#把GO_MF富集结果导出表格
BPtable<-as.data.frame(engo_BP)
write.table(BPtable,file="result_GO_MF.txt",sep="\t",quote=F)#把GO_BP富集结果导出表格
7.多样本ChIPseeker分析脚本
#工作目录，样本种属来源需要根据情况调整
setwd("/Users/meraner/Desktop/bed") #制定当前工作目录，为数据存储目录
library(ChIPseeker) #加载ChIPseeker包
library(ggplot2) #加载ggplot2包
library(org.Mm.eg.db) #加载参考基因组数据
library(TxDb.Mmusculus.UCSC.mm10.knownGene) #加载参考基因组注释数据包
library(clusterProfiler) #加载clusterProfiler包
library(GenomicFeatures) #加载GenomicFeatures包
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
#peak的整体分布
peak=GenomicRanges::GRangesList("8WG16_1"=readPeakFile("8WG16-1.sorted.bam_peaks.narrowPeak"),"8WG16_2"=readPeakFile("8WG16-2.sorted.bam_peaks.narrowPeak"))
covplot(peak, weightCol="V5") + facet_grid(chr ~ .id)#coverageplot
covplot(peak, weightCol="V5", chrs=c("chr17", "chr18"), xlim=c(4e7, 5e7)) + facet_grid(chr ~ .id)#截取某些染色体，某些片段进行展示
#多个样本，在启动子区(TSS)附近的 ChIP peaks 的可视化
peak1 <- "/Users/meraner/Desktop/bed/8WG16-1.sorted.bam_summits.bed"
peak2 <- "/Users/meraner/Desktop/bed/8WG16-1.sorted.bam_summits.bed"
files <-list(peak1=peak1,peak2=peak2)#把两个bed文件指定给files, 就是一次读入多个summits文件，使用list存储，然后使用lapply注释。
promoter <- getPromoters (TxDb=txdb, upstream=3000, downstream=3000) #getPromoters函数定义启动子区域的范围
files <-list(peak1=peak1,peak2=peak2)#把两个bed文件指定给files, 就是一次读入多个summits文件，使用list存储，然后使用lapply注释。
tagMatrixList <- lapply(files, getTagMatrix, windows=promoter)
peakHeatmap(files, TxDb=txdb,
upstream=3000, downstream=3000,
color=rainbow(length(files))) #多个样本在TSS附近结合区的热图
plotAvgProf(tagMatrixList, xlim=c(-3000, 3000),
conf=0.95,resample=500, facet="row") # 添加置信区间并分面
#注释，结构可视化
peakAnnoList <- lapply(files, annotatePeak, TxDb=txdb,
tssRegion=c(-3000, 3000),
verbose=FALSE) #注释
plotAnnoBar(peakAnnoList) #画图bar
plotDistToTSS(peakAnnoList) #画图距离
#单样本可以做饼图和韦恩图，多样本的不行。如果需要，只能用单样本流程。
#功能可视化-基因注释 + 富集分析
require(clusterProfiler)
seq=lapply(files, readPeakFile)
genes <-lapply(seq, function(i)
seq2gene(i, c(-1000, 3000), 3000, TxDb=txdb))
cc <- compareCluster(geneClusters = genes, organism="mmu",
fun="enrichKEGG")
dotplot(cc, showCategory=10)
#输出注释信息
peakAnnotable<-as.data.frame(peakAnnoList)
write.table(peakAnnotable,file="result_multisample_gene",sep="\t",quote=F) #导出基因信息
cctable<-as.data.frame(cc)
write.table(cctable,file="result_multisample_kegg.txt",sep="\t",quote=F)#把KEGG富集结果导出表格
