学习一遍ChIPseeker的使用

https://www.jianshu.com/p/c76e83e6fa57
刘小泽写于2020.5.23-24
Y叔的原文在：https://mp.weixin.qq.com/s/3CMj0xejiV-FSMC-Vxd_-w

0 ChIPseeker的诞生

Y叔一开始使用ChIPpeakAnno进行注释，但使用UCSC genome browser检验结果的时候，发现对不上；另外之前在使用ChIPpeakAnno过程中写了一些可视化函数。后来经过漫长的半夜宿舍苦战，写出了ChIPseeker

1 ChIP-seq简介

ChIP是指染色质免疫沉淀，它通特异结合抗体将DNA结合蛋白免疫沉淀，可以用于捕获蛋白质（如转录因子，组蛋白修饰）的DNA靶点。之前结合芯片就有ChIP-on-chip，后来二代测序加持诞生了ChIP-seq。优点是：不再需要设计探针（探针往往存在着一定的偏向性）

2007年来自三个不同的实验室，几乎是同时间出来（最长差不了3个月），分别发CNS，一起定义了这个ChIPseq技术

Johnson DS, Mortazavi A et al. (2007) Genome-wide mapping of in vivo protein–DNA interactions. Science 316: 1497–1502
Robertson G et al.(2007) Genome-wide profiles of STAT1 DNA association using chromatin immunoprecipitation and massively parallel sequencing. Nature Methods 4: 651–657
Schmid et al. (2007) ChIP-Seq Data reveal nucleosome architecture of human promoters. Cell 131: 831–832

主要有4步：Cross-linking、Sonication、IP、Sequencing

简而言之是：DNA和蛋白质交联(cross-linking)、超声(sonication)将染色体随机切割、利用抗原抗体的特异性识别(IP)、把目标蛋白相结合的DNA片段沉淀下来，反交联释放DNA片段，最后是测序(sequencing)

分析流程示例图1：

分析流程示例图2：

原始数据=》质控=》比对=》拿到DNA片段在染色体上的位置信息=》peak calling （去除背景噪音）=》拿到peaks（protein binding site）=》下游分析（可视化、找相关基因、motif分析等等）

2 必须知晓的BED文件

全称是：Browser Extensible Data，为基因组浏览器而生

包括3个必须字段和9个可选字段：

3个必须

1 chrom - 染色体名字
2 chromStart - 染色体起始位点（起始于0，而不是1）许多软件忽略了这一点，存在一个碱基的位移（如peakAnalyzer, ChIPpeakAnno存在这个问题），Homer、ChIPseeker没有
3 chromEnd - 染色体终止位点

9个可选

4 name - 名字
5 score - 分值(0-1000), 用于genome browser展示时上色。
6 strand - 正负链，对于ChIPseq数据来说，一般没有正负链信息。
7 thickStart - 画矩形的起点
8 thickEnd - 画矩形的终点
9 itemRgb - RGB值
10 blockCount - 子元件（比如外显子）的数目
11 blockSizes - 子元件的大小
12 blockStarts - 子元件的起始位点

一般只用前5个足矣（MACS的输出结果也是前5个字段）

第5列score的含义是：the summit height of fragment pileup. 也即是片段堆积的峰高

3 使用`covplot`可视化BED数据

一般拿到数据后，会先可视化一下数据的全景

# 自带示例数据（这也是Bioconductor包的一个特点，提交R包需要有说明书和测试数据）
library(ChIPseeker)
library(ggplot2)

files <- getSampleFiles()
# 有5个文件
> basename(unlist(files))
[1] "GSM1174480_ARmo_0M_peaks.bed.gz"                   
[2] "GSM1174481_ARmo_1nM_peaks.bed.gz"                  
[3] "GSM1174482_ARmo_100nM_peaks.bed.gz"                
[4] "GSM1295076_CBX6_BF_ChipSeq_mergedReps_peaks.bed.gz"
[5] "GSM1295077_CBX7_BF_ChipSeq_mergedReps_peaks.bed.gz"

covplot(files[[5]])

还支持多个文件同时画

只要转为GRanges对象即可

# 比如要画第4、5个文件（MACS生成的BED文件包含常规的5列）
peak=GenomicRanges::GRangesList(CBX6=readPeakFile(files[[4]]),CBX7=readPeakFile(files[[5]]))

画图

covplot(peak, weightCol="V5") + facet_grid(chr ~ .id)

取小区间，例如只取几条染色体，还能定义染色体的区间大学

covplot(peak, weightCol="V5", chrs=c("chr17", "chr18"), 
        xlim=c(4e7, 5e7)) + facet_grid(chr ~ .id)

在看完数据全景之后，就会想知道这些peaks和什么类型的基因有关

4 annotatePeak进行peaks的注释

需要使用BED文件（作为query）+注释文件（作为target）

重点是如何获取注释文件

注释信息一般要包含基因的起始终止，基因的外显子、内含子及它们的起始终止、非编码区域位置、功能元件的位置等

ChIPseeker没有物种限制，但前提是物种本身有这些注释信息（不能说物种连参考基因组也没有，那就真的是巧妇难为无米之炊）

需要一个TxDb对象，例如TxDb.Hsapiens.UCSC.hg19.knownGene，然后ChIPseeker就会从中提取信息

# 三步走（提供TxDb注释、提供bed文件、进行注释）
require(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb = TxDb.Hsapiens.UCSC.hg19.knownGene
f = getSampleFiles()[[4]]
x = annotatePeak(f, tssRegion=c(-1000, 1000), TxDb=txdb)

看到这里有个参数tssRegion ，它指定了启动子区域（而启动子区域是没有明确定义的，需要自己指定，这里指定了上下游1kb）

看一下大体结果：

> x
Annotated peaks generated by ChIPseeker
1331/1331  peaks were annotated
Genomic Annotation Summary:
             Feature  Frequency
9           Promoter 48.1592787
4             5' UTR  0.7513148
3             3' UTR  4.2073629
1           1st Exon  0.7513148
7         Other Exon  3.9068370
2         1st Intron  6.5364388
8       Other Intron  4.8835462
6 Downstream (<=300)  1.1269722
5  Distal Intergenic 29.6769346

看一下详细结果：

as.GRanges(x) %>% head(3)

可以转为数据框，方便输出：

tmp=as.data.frame(x)

关于注释的类型：

注释类型一：genomic annotation（annotation这一列）

指peak在基因组的位置：落在什么地方，例如外显子、内含子或是UTR

注释类型二：nearest gene annotation（annotation后面的列）

指peak最近的基因：不管peak落在内含子、基因间区还是其他位置，按照peak相对于转录起始位点的距离，都能找到一个离它最近的基因【一般做基因表达调控的，会关注promoter区域，离结合位点最近的基因更可能被调控】

这个距离是根据转录起始位点来计算，一个基因具有多个转录本，因此一个基因可能有多个转录起始位点。注释的结果就会看到有一列是转录本ID

注释类型三：flankDistance（三列: flank_txIds, flank_geneIds和flank_gene_distances）

指peak上下游某个范围内（比如-5kb《=》5kb范围内）都有什么基因

# 传个参数flankDistance
x2 = annotatePeak(f, tssRegion=c(-1000, 1000), TxDb=txdb, addFlankGeneInfo=TRUE, flankDistance=5000)

让基因名变得友好

上面得到的结果都是以geneId（Entrez ID）给出，如果想要Symbol名称，可以再传参数annoDb

library(org.Hs.eg.db)
x3 = annotatePeak(f, tssRegion=c(-1000, 1000), TxDb=txdb, 
                  addFlankGeneInfo=TRUE, flankDistance=5000,
                  annoDb = "org.Hs.eg.db")
tmp3=as.data.frame(x3)

会再增加3列：ENSEMBL、SYMBOL、GENENAME（如果这里使用的TxDb是Ensemble ID，那么结果就会是Entrez ID、SYMBOL、GENENAME三列）

按正负链分开注释

一般ChIPseq数据通常情况下是没有正负链信息的（有特殊的实验可以有）

但如果要做，可以先给peaks分别赋予正负链的信息，然后指定参数sameStrand=TRUE 并分别做两次

这个参数的意思是：（logical）whether find nearest/overlap gene in the same strand

只注释基因的上游或下游

提供了ignoreDownstream和 ignoreUpstream，默认是FALSE

关于TxDb的知识

上面一起操作的前提是物种本身有这些注释信息，而注释信息主要是用TxDb

同一物种的不同版本TxDb

例如TxDb.Hsapiens.UCSC.hg19.knownGene和TxDb.Hsapiens.UCSC.hg38.knownGene 的注释结果是不同的，不能混用。用哪个取决于上游分析比对使用的哪个版本的基因组

不同的版本中基因坐标是不一样的，如果硬要替换，可以使用liftOver将基因组版本坐标进行转换

支持多少物种？

Bioconductor上有30个左右TxDb，也只能覆盖一小部分物种（https://bioconductor.org/packages/3.11/data/annotation/），但UCSC和Ensemble的基因组都可以被ChIPseeker支持，因此所有物种都支持

除了基因注释还能注释lincRNA

比如就可以利用：https://bioconductor.org/packages/3.11/data/annotation/html/TxDb.Hsapiens.UCSC.hg19.lincRNAsTranscripts.html

require("TxDb.Hsapiens.UCSC.hg19.lincRNAsTranscripts")
linc_txdb=TxDb.Hsapiens.UCSC.hg19.lincRNAsTranscripts
x=annotatePeak(peak, TxDb=linc_txdb)
as.GRanges(x)

如何自己制作TxDb？

使用GenomicFeatures包来制作TxDb对象

makeTxDbFromUCSC：通过UCSC在线制作TxDb
makeTxDbFromBiomart: 通过ensembl在线制作TxDb
makeTxDbFromGRanges：通过GRanges对象制作TxDb
makeTxDbFromGFF：通过解析GFF文件制作TxDb

比如在线从UCSC生成TxDb：

require(GenomicFeatures)
# makeTxDbFromUCSC()函数依赖RMariaDB这个包
# BiocManager::install('RMariaDB')
hg19.refseq.db <- makeTxDbFromUCSC(genome="hg19", table="refGene")
# 可能会遇到一个报错：namespace ‘DBI’ 1.0.0 is already loaded, but >= 1.1.0 is required =>自己升级
# remove.packages("DBI", lib="~/Library/R/3.6/library")
# packageurl <- "https://cran.r-project.org/src/contrib/DBI_1.1.0.tar.gz"
# install.packages(packageurl, repos=NULL, type="source")

然后可以对比一下：

再比如自己下载GTF然后生成TxDb

以大豆（glycine_max）为例

# 下载
download.file('ftp://ftp.ensemblgenomes.org/pub/plants/release-47/gff3/glycine_max/Glycine_max.Glycine_max_v2.1.47.chr.gff3.gz',destfile = 'Glycine_max_v2.1.47.chr.gff3.gz')
# 解压
R.utils::gunzip('Glycine_max_v2.1.47.chr.gff3.gz')
# 制作
glycine <- makeTxDbFromGFF("Glycine_max_v2.1.47.chr.gff3")

有了TxDb怎么查看呢？

最详细的操作在官方文档：https://bioconductor.org/packages/release/bioc/vignettes/GenomicFeatures/inst/doc/GenomicFeatures.pdf

不管是从Bioconductor下载的还是自己制作的，都是一个GenomicFeatures对象

如果简单对名称操作，会返回这个注释文件的基本信息。要把TxDb当成一个数据库来对待，而不是一个简单的数据框或者矩阵。因此它的提取方法也会比较特别

如果想看其中包含的类目，可以用columns(txdb)
如果想指定提取转录本或外显子信息，可以：transcripts(txdb) 或者 exons(txdb)
如果想看全部的信息，可以：AnnotationDbi::select(glycine, columns=columns(glycine), keys=keys(glycine), keytype=c("GENEID"))

需要注意，如果使用这个select的时候，同时加载了tidyverse，那么同名的select就会发生冲突导致报错，这时可以用显式指定的形式来规范（如下图）

可视化

peak在整个染色体的分布

见：第3部分=》使用covplot可视化BED数据

peak在某个窗口的结合谱图

一般有两种方式：一是直接使用BED文件，二是一步步手动进行

第一种：直接使用BED文件

peakHeatmap(f, weightCol="V5", TxDb=txdb, 
            upstream=3000, downstream=3000, 
            color=rainbow(length(f)))

其实看运行日志也能看出来做了什么，首先根据转录起始位点指定上下游（也就是热图的窗口区间范围），然后把peaks比对到这个窗口，并生成矩阵以进行可视化

稍微查看一下这个peakHeatmap函数，就会发现以上说的几步：

当然，如果是多个文件也是可以的

files=getSampleFiles()
peakHeatmap(files, TxDb=txdb, 
            upstream=3000, downstream=3000, 
            color=rainbow(length(files)))

第二种：一步步手动进行

如果说第一种提供了一个打包好的计算过程，那么第二种就是把第一种拆分运行

promoter <- getPromoters(TxDb=txdb, 
                  upstream=3000, downstream=3000)
tagMatrix <- getTagMatrix(f, 
                          windows=promoter)
tagHeatmap(tagMatrix, xlim=c(-3000, 3000), 
           color="red")

看看结合的强度

第一种：直接使用BED文件

plotAvgProf2(files[[4]], TxDb=txdb, 
             upstream=3000, downstream=3000,
             xlab="Genomic Region (5'->3')", 
             ylab = "Read Count Frequency",
             conf = 0.95, resample = 1000)# 添加置信区间

第二种：手动进行

使用上面的tagMatrix计算结果

plotAvgProf(tagMatrix, xlim=c(-3000, 3000),
            xlab="Genomic Region (5'->3')", 
            ylab = "Read Count Frequency")

支持多个数据比较

tagMatrixList <- lapply(files, getTagMatrix, 
                        windows=promoter)
# 添加置信区间并分面
plotAvgProf(tagMatrixList, xlim=c(-3000, 3000), 
            conf=0.95,resample=500, facet="row")

这个结果和上面peakHeatmap的结果一致，前3个样本不是调控转录的

除了关注转录起始位点（研究转录调控），还能看蛋白与外显子/内含子起始位置的结合谱，使用getBioRegion函数，可以指定'gene', 'transcript', 'exon', 'intron'

注释结果之注释类型一：genomic annotation

指peak在基因组的位置：落在什么地方，例如外显子、内含子或是5’ /3‘UTR

饼图

peakAnno <- annotatePeak(files[[4]], 
                         tssRegion=c(-3000, 3000),
                         TxDb=txdb, annoDb="org.Hs.eg.db")
plotAnnoPie(peakAnno)

柱状图

plotAnnoBar(peakAnno)

注意第一个问题：关于上图中的各种Features分类

看到这里的分类有下游（Downstream）但没有上游，这是因为Promoter定义为了转录起始位点（TSS）的上下游区域，包含了上游；另外这个下游是是基因间区的一部分，更确切是指紧接着基因的下游；这里的上游和下游其实都是基因间区，单独拿出来是因为和基因直接连接，是很近的区域=》近端基因间区

当然，基因间区还包含更远的间区（Distal intergenic）=》远端基因间区

默认下游的范围是3kb，但是可以自己调整

# 比如调成500
options(ChIPseeker.downstreamDistance = 500)

还有一个需求就是：自定义分类

# 依然是设置options，用于总结结果
f2=getSampleFiles()[[5]]
options(ChIPseeker.ignore_1st_exon = T)
options(ChIPseeker.ignore_1st_intron = T)
options(ChIPseeker.ignore_downstream = T)
options(ChIPseeker.ignore_promoter_subcategory = T)
x=annotatePeak(f2)
plotAnnoPie(x)

注意第二个问题：peak的位置可能不是唯一的

这是因为，一个peak所在的位置，可能是一个基因的外显子，同时又是另一个基因的内含子。为了解决这个问题，有以下几种方案：

第一种：使用参数genomicAnnotationPriority指定优先顺序

默认顺序是：Promoter => 5’ UTR => 3’ UTR => Exon => Intron => Downstream => Distal Intergenic

第二种：饼图+韦恩图

vennpie(peakAnno)

优点是：直观；缺点是：无法显示全部的信息

第三种：UpSetR + vennpie

upsetplot(peakAnno, vennpie=TRUE)

多个文件的区域注释

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

注释结果之注释类型二：nearest gene annotation

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

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

同样也支持多个文件

plotDistToTSS(peakAnnoList)

距离最近的基因在不同样本的交集

# 先得到基因列表
genes <- lapply(peakAnnoList, function(i) 
    as.data.frame(i)$geneId)
> names(genes)
[1] "ARmo_0M"    "ARmo_1nM"   "ARmo_100nM" "CBX6_BF"    "CBX7_BF"  
                
# 然后作图(需要借助Vennerable包)
devtools::install_github("js229/Vennerable")
library(Vennerable)
vennplot(genes[2:4], by='Vennerable')

基因注释 + 富集分析

利用ChIPseeker的seq2gene 将peak的位置与所有的基因关联起来【包括 host gene (exon/intron), promoter region and flanking gene from intergenomic region】，然后用clusterProfiler拿这些基因跑ORA，做富集

require(clusterProfiler)
bedfile=getSampleFiles()
# 将bed文件读入（readPeakFile是利用read.delim读取，然后转为GRanges对象）
seq=lapply(bedfile, readPeakFile)

genes=lapply(seq, function(i) 
    seq2gene(i, c(-1000, 3000), 3000, TxDb=txdb))
cc = compareCluster(geneClusters = genes, 
                    fun="enrichKEGG", organism="hsa")
dotplot(cc, showCategory=10)

学习一遍ChIPseeker的使用

0 ChIPseeker的诞生

1 ChIP-seq简介

分析流程示例图1：

分析流程示例图2：

2 必须知晓的BED文件

3个必须

9个可选

3 使用covplot可视化BED数据

还支持多个文件同时画

4 annotatePeak进行peaks的注释

重点是如何获取注释文件

看一下大体结果：

看一下详细结果：

可以转为数据框，方便输出：

关于注释的类型：

注释类型一：genomic annotation（annotation这一列）

注释类型二：nearest gene annotation（annotation后面的列）

注释类型三：flankDistance（三列: flank_txIds, flank_geneIds和flank_gene_distances）

让基因名变得友好

按正负链分开注释

只注释基因的上游或下游

关于TxDb的知识

同一物种的不同版本TxDb

支持多少物种？

除了基因注释还能注释lincRNA

如何自己制作TxDb？

有了TxDb怎么查看呢？

可视化

peak在整个染色体的分布

peak在某个窗口的结合谱图

第一种：直接使用BED文件

第二种：一步步手动进行

看看结合的强度

第一种：直接使用BED文件

第二种：手动进行

支持多个数据比较

注释结果之注释类型一：genomic annotation

饼图

柱状图

注意第一个问题：关于上图中的各种Features分类

注意第二个问题：peak的位置可能不是唯一的

多个文件的区域注释

注释结果之注释类型二：nearest gene annotation

同样也支持多个文件

距离最近的基因在不同样本的交集

基因注释 + 富集分析

3 使用`covplot`可视化BED数据