软件16 —— ChIPseeker

一、 基本介绍

ChIPseeker可以用来对ChIP-seq数据进行注释与可视化。

二、 使用方法

(1) 安装并载入R包

ls() #列出当前工作空间中的对象
rm(list=ls()) #移除全部对象

options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")) 
options(BioC_mirror="http://mirrors.tuna.tsinghua.edu.cn/bioconductor")

if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

#BiocManager::install("ChIPpeakAnno")
BiocManager::install("ChIPseeker")
BiocManager::install("TxDb.Dmelanogaster.UCSC.dm6.ensGene")
BiocManager::install("org.Dm.eg.db")

#library(ChIPpeakAnno)
library(ChIPseeker)
library(TxDb.Dmelanogaster.UCSC.dm6.ensGene)
library(org.Dm.eg.db)
library(ggplot2)
library(clusterProfiler)

(2) bed文件处理

$ ls 11_intersect/*.bed | while read id ; do (cat $id | awk '{if($1=="2L"||$1=="2R"||$1=="3L"||$1=="3R"||$1=="4"||$1=="X"||$1=="Y") {print "chr"$1"\t"$2"\t"$3"\t"$4}}' > 12_ChIPseeker/bed/$(basename $id)) ; done  &

(3) 峰在整个染色体的分布

peak1 <- readPeakFile("12_ChIPseeker/bed/WTF_intersect.bed")  #GRanges object
peak2 <- readPeakFile("12_ChIPseeker/bed/WTM_intersect.bed")

peaks <- list(WTF=peak1, WTM=peak2)

covplot(peak1)
covplot(peaks) + facet_grid(chr ~ .id)
covplot(peaks, chrs=c("chrX"), xlim=c(1e7,1.1e7)) + facet_grid(chr ~ .id)

(4) 峰注释

txdb = TxDb.Dmelanogaster.UCSC.dm6.ensGene

peakAnno <- annotatePeak(peak1,
                         tssRegion = c(-3000, 3000),  #启动子区域
                         genomicAnnotationPriority = c("Promoter", "5UTR", "3UTR", "Exon", "Intron","Downstream", "Intergenic"),
                         flankDistance = 5000,
                         #sameStrand = T,
                         TxDb = txdb,
                         annoDb = "org.Dm.eg.db") 

peakAnno
peakAnno_df <- as.data.frame(peakAnno)
head(peakAnno_df)

write.csv(peakAnno_df,file="12_ChIPseeker/WTF_peakAnno.csv",row.names=F)

注释结果：

seqnames  start  end  width  strand  V4  annotation  geneChr  geneStart  geneEnd geneLength  geneStrand  geneId  transcriptId  distanceToTSS  ENTREZID  SYMBOL  GENENAME
chrX  138244  138423  180  *  wt_msl2_peak_1  Exon (FBtr0340053/FBgn0025836, exon 9 of 14)  6  140722  142197  1476  2  FBgn0267029  FBtr0345995  3774  19834782  lncRNA:CR45473  long non-coding RNA:CR45473

这里注释有两种，一种是genomic annotation (也就是annotation这一列），还有就是nearest gene annotation（也就是多出来的其它列）。genomic annotation注释的是peak的位置，它落在什么地方了，可以是UTR，可以是内含子或者外显子。而最近基因是peak相对于转录起始位点的距离，不管这个peak是落在内含子或者别的什么位置上，即使它落在基因间区上，我都能够找到一个离它最近的基因（即使它可能非常远）。

针对不同的问题，这两种注释的策略是不一样的。第一种策略peak所在位置，可能就是调控的根本，比如你要做可变剪切的，最近基因的注释显然不是你关注的点。而做基因表达调控的，当然promoter区域是重点，离结合位点最近的基因更有可能被调控。最近基因的注释信息虽然是以基因为单位给出，但我们针对的是转录起始位点来计算距离，针对于不同的转录本，一个基因可能有多个转录起始位点，所以注释是在转录本的水平上进行的，我们可以看到输出有一列是transcriptId。

两种注释有时候还不够，我想看peak上下游某个范围内（比如说-5k到5k的距离）都有什么基因，annotatePeak也可以做到。你只要传个参数说你要这个信息，还有什么区间内，就可以了。输出中多三列：flank_txIds、flank_geneIds和flank_gene_distances，在指定范围内所有的基因都被列出。

最近基因给出来了，但都是各种人类不友好的ID，只需要给annotatePeak传入annoDb参数就行了。如果你的TxDb的基因ID是Entrez，它会转成ENSEMBL，反之亦然，当然不管是那一种，都会给出SYMBOL，还有描述性的GENENAME。

The annotation column annotates genomic features of the peak, that is whether peak is located in Promoter, Exon, UTR, Intron or Intergenic. If the peak is annotated by Exon or Intron, more detail information will be provided. For instance, Exon (38885 exon 3 of 11) indicates that the peak is located at the 3th exon of gene 38885 (entrezgene ID) which contain 11 exons in total.

一个peak所在的位置，可能是一个基因的外显子，同时又是另一个基因的内含子。使用参数genomicAnnotationPriority指定优先顺序，默认顺序是：Promoter => 5’ UTR => 3’ UTR => Exon => Intron => Downstream => Distal Intergenic。

Downstream is defined as the downstream of gene end. Promoter定义为转录起始位点（TSS）的上下游区域。另外下游是基因间区的一部分，更确切是指紧接着基因的下游；这里的上游和下游其实都是基因间区，单独拿出来是因为和基因直接连接，是很近的区域，即近端基因间区。当然，基因间区还包含更远的间区（Distal intergenic），即远端基因间区。

(5) 比对到TSS附近的峰分布

# plot the heatmap of peaks align to flank sequences of TSS.
#和deeptools的reference-point mode差不多，但是画的图有点丑
peakHeatmap(peak1,
            TxDb=txdb, 
            upstream=3000, downstream=3000,  #指定转录起始位点上下游
            color=rainbow(length(peak1)))

peakHeatmap(peaks,  #list
            TxDb=txdb, 
            upstream=3000, downstream=3000, 
            color=rainbow(length(peaks)))

# plot the profile of peaks that align to flank sequences of TSS.
plotAvgProf2(peak1, 
             TxDb=txdb, 
             upstream=3000, downstream=3000,
             xlab="Genomic Region (5'->3')", 
             ylab = "Read Count Frequency",
             conf = 0.95, resample = 1000)

plotAvgProf2(peaks, 
             TxDb=txdb, 
             upstream=3000, downstream=3000,
             xlab="Genomic Region (5'->3')", 
             ylab = "Read Count Frequency",
             conf = 0.95, resample = 1000) + facet_grid(. ~ .id)

(6) 峰在基因组的位置

peakAnno <- annotatePeak(peak1, 
                         tssRegion=c(-3000, 3000),
                         TxDb=txdb,
                         annoDb="org.Dm.eg.db")
plotAnnoPie(peakAnno)

peakAnnoList <- lapply(peaks, annotatePeak, 
                       TxDb=txdb,
                       tssRegion=c(-3000, 3000))
plotAnnoBar(peakAnnoList)

(7) 相对于转录起始点的距离

相对于转录起始位点的距离找最近的基因。不管peak落在内含子、基因间区还是其他位置，按照peak相对于转录起始位点的距离，都能找到一个离它最近的基因。

plotDistToTSS(peakAnno,
              title="Distribution of transcription factor-binding loci\nrelative to TSS")
plotDistToTSS(peakAnnoList,
              title="Distribution of transcription factor-binding loci\nrelative to TSS")

(8) 距离最近的基因的交集和功能

genes <- lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)
names(genes)