学习一遍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数据
covplot
可视化BED数据一般拿到数据后,会先可视化一下数据的全景
还支持多个文件同时画
只要转为GRanges对象即可
画图
取小区间,例如只取几条染色体,还能定义染色体的区间大学
在看完数据全景之后,就会想知道这些peaks和什么类型的基因有关
4 annotatePeak进行peaks的注释
需要使用BED文件(作为query)+注释文件(作为target)
重点是如何获取注释文件
注释信息一般要包含基因的起始终止,基因的外显子、内含子及它们的起始终止、非编码区域位置、功能元件的位置等
ChIPseeker没有物种限制,但前提是物种本身有这些注释信息(不能说物种连参考基因组也没有,那就真的是巧妇难为无米之炊)
需要一个TxDb对象,例如TxDb.Hsapiens.UCSC.hg19.knownGene
,然后ChIPseeker就会从中提取信息
看到这里有个参数tssRegion
,它指定了启动子区域(而启动子区域是没有明确定义的,需要自己指定,这里指定了上下游1kb)
看一下大体结果:
看一下详细结果:
可以转为数据框,方便输出:
关于注释的类型:
注释类型一: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范围内)都有什么基因
让基因名变得友好
上面得到的结果都是以geneId(Entrez ID)给出,如果想要Symbol名称,可以再传参数annoDb
会再增加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
如何自己制作TxDb?
使用GenomicFeatures
包来制作TxDb对象
makeTxDbFromUCSC: 通过UCSC在线制作TxDb
makeTxDbFromBiomart: 通过ensembl在线制作TxDb
makeTxDbFromGRanges:通过GRanges对象制作TxDb
makeTxDbFromGFF:通过解析GFF文件制作TxDb
比如在线从UCSC生成TxDb:
然后可以对比一下:
再比如自己下载GTF然后生成TxDb
以大豆(glycine_max)为例
有了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文件
其实看运行日志也能看出来做了什么,首先根据转录起始位点指定上下游(也就是热图的窗口区间范围),然后把peaks比对到这个窗口,并生成矩阵以进行可视化
稍微查看一下这个peakHeatmap
函数,就会发现以上说的几步:
当然,如果是多个文件也是可以的
第二种:一步步手动进行
如果说第一种提供了一个打包好的计算过程,那么第二种就是把第一种拆分运行
看看结合的强度
第一种:直接使用BED文件
第二种:手动进行
使用上面的tagMatrix计算结果
支持多个数据比较
这个结果和上面peakHeatmap
的结果一致,前3个样本不是调控转录的
除了关注转录起始位点(研究转录调控),还能看蛋白与外显子/内含子起始位置的结合谱,使用getBioRegion
函数,可以指定'gene', 'transcript', 'exon', 'intron'
注释结果之注释类型一:genomic annotation
指peak在基因组的位置:落在什么地方,例如外显子、内含子或是5’ /3‘UTR
饼图
柱状图
注意第一个问题:关于上图中的各种Features分类
看到这里的分类有下游(Downstream)但没有上游,这是因为Promoter定义为了转录起始位点(TSS)的上下游区域,包含了上游;另外这个下游是是基因间区的一部分,更确切是指紧接着基因的下游;这里的上游和下游其实都是基因间区,单独拿出来是因为和基因直接连接,是很近的区域=》近端基因间区
当然,基因间区还包含更远的间区(Distal intergenic)=》远端基因间区
默认下游的范围是3kb,但是可以自己调整
还有一个需求就是:自定义分类
注意第二个问题:peak的位置可能不是唯一的
这是因为,一个peak所在的位置,可能是一个基因的外显子,同时又是另一个基因的内含子。为了解决这个问题,有以下几种方案:
第一种:使用参数
genomicAnnotationPriority
指定优先顺序
默认顺序是:Promoter => 5’ UTR => 3’ UTR => Exon => Intron => Downstream => Distal Intergenic
第二种:饼图+韦恩图
优点是:直观;缺点是:无法显示全部的信息
第三种:UpSetR + vennpie
多个文件的区域注释
注释结果之注释类型二:nearest gene annotation
指peak最近的基因:不管peak落在内含子、基因间区还是其他位置,按照peak相对于转录起始位点的距离,都能找到一个离它最近的基因
同样也支持多个文件
距离最近的基因在不同样本的交集
基因注释 + 富集分析
利用ChIPseeker的seq2gene
将peak的位置与所有的基因关联起来【包括 host gene (exon/intron), promoter region and flanking gene from intergenomic region】,然后用clusterProfiler拿这些基因跑ORA,做富集
Last updated