190-学习一遍ChIPseeker使用

刘小泽写于2020.5.23 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

只注释基因的上游或下游

提供了ignoreDownstreamignoreUpstream,默认是FALSE

关于TxDb的知识

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

同一物种的不同版本TxDb

例如TxDb.Hsapiens.UCSC.hg19.knownGeneTxDb.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)

Yunze Liu
Yunze Liu
Bioinformatics Sharer

Co-founder of Bioinfoplanet(生信星球)

Next
Previous

Related