劉小澤寫于2020.5.23-24
Y叔的原文在：https://mp.weixin.qq.com/s/3CMj0xejiV-FSMC-Vxd_-w
0 ChIPseeker的誕生
Y叔一開始使用ChIPpeakAnno進(jìn)行注釋，但使用UCSC genome browser檢驗結(jié)果的時候，發(fā)現(xiàn)對不上；另外之前在使用ChIPpeakAnno過程中寫了一些可視化函數(shù)。后來經(jīng)過漫長的半夜宿舍苦戰(zhàn)，寫出了ChIPseeker
1 ChIP-seq簡介
ChIP是指染色質(zhì)免疫沉淀，它通特異結(jié)合抗體將DNA結(jié)合蛋白免疫沉淀，可以用于捕獲蛋白質(zhì)（如轉(zhuǎn)錄因子，組蛋白修飾）的DNA靶點。之前結(jié)合芯片就有ChIP-on-chip，后來二代測序加持誕生了ChIP-seq。優(yōu)點是：不再需要設(shè)計探針（探針往往存在著一定的偏向性）
2007年來自三個不同的實驗室，幾乎是同時間出來（最長差不了3個月），分別發(fā)CNS，一起定義了這個ChIPseq技術(shù)
- 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和蛋白質(zhì)交聯(lián)(cross-linking)、超聲(sonication)將染色體隨機(jī)切割、利用抗原抗體的特異性識別(IP)、把目標(biāo)蛋白相結(jié)合的DNA片段沉淀下來，反交聯(lián)釋放DNA片段，最后是測序(sequencing)
分析流程示例圖1：
分析流程示例圖2：
原始數(shù)據(jù)=》質(zhì)控=》比對=》拿到DNA片段在染色體上的位置信息=》peak calling （去除背景噪音）=》拿到peaks（protein binding site）=》下游分析（可視化、找相關(guān)基因、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 - 正負(fù)鏈，對于ChIPseq數(shù)據(jù)來說，一般沒有正負(fù)鏈信息。
- 7 thickStart - 畫矩形的起點
- 8 thickEnd - 畫矩形的終點
- 9 itemRgb - RGB值
- 10 blockCount - 子元件（比如外顯子）的數(shù)目
- 11 blockSizes - 子元件的大小
- 12 blockStarts - 子元件的起始位點
一般只用前5個足矣（MACS的輸出結(jié)果也是前5個字段）
第5列score的含義是：the summit height of fragment pileup. 也即是片段堆積的峰高
3 使用covplot可視化BED數(shù)據(jù)
一般拿到數(shù)據(jù)后，會先可視化一下數(shù)據(jù)的全景
# 自帶示例數(shù)據(jù)（這也是Bioconductor包的一個特點，提交R包需要有說明書和測試數(shù)據(jù)）
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]])
還支持多個文件同時畫
只要轉(zhuǎn)為GRanges對象即可
# 比如要畫第4、5個文件（MACS生成的BED文件包含常規(guī)的5列）
peak=GenomicRanges::GRangesList(CBX6=readPeakFile(files[[4]]),CBX7=readPeakFile(files[[5]]))
畫圖
covplot(peak, weightCol="V5") + facet_grid(chr ~ .id)
取小區(qū)間，例如只取幾條染色體，還能定義染色體的區(qū)間大學(xué)
covplot(peak, weightCol="V5", chrs=c("chr17", "chr18"),
xlim=c(4e7, 5e7)) + facet_grid(chr ~ .id)
在看完數(shù)據(jù)全景之后，就會想知道這些peaks和什么類型的基因有關(guān)
4 annotatePeak進(jìn)行peaks的注釋
需要使用BED文件（作為query）+注釋文件（作為target）
重點是如何獲取注釋文件
注釋信息一般要包含基因的起始終止，基因的外顯子、內(nèi)含子及它們的起始終止、非編碼區(qū)域位置、功能元件的位置等
ChIPseeker沒有物種限制，但前提是物種本身有這些注釋信息（不能說物種連參考基因組也沒有，那就真的是巧婦難為無米之炊）
需要一個TxDb對象，例如TxDb.Hsapiens.UCSC.hg19.knownGene，然后ChIPseeker就會從中提取信息
# 三步走（提供TxDb注釋、提供bed文件、進(jìn)行注釋）
require(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb = TxDb.Hsapiens.UCSC.hg19.knownGene
f = getSampleFiles()[[4]]
x = annotatePeak(f, tssRegion=c(-1000, 1000), TxDb=txdb)
看到這里有個參數(shù)tssRegion ，它指定了啟動子區(qū)域（而啟動子區(qū)域是沒有明確定義的，需要自己指定，這里指定了上下游1kb）
看一下大體結(jié)果：
> 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
看一下詳細(xì)結(jié)果：
as.GRanges(x) %>% head(3)
可以轉(zhuǎn)為數(shù)據(jù)框，方便輸出：
tmp=as.data.frame(x)
關(guān)于注釋的類型：
注釋類型一：genomic annotation（annotation這一列）
指peak在基因組的位置：落在什么地方，例如外顯子、內(nèi)含子或是UTR
注釋類型二：nearest gene annotation（annotation后面的列）
指peak最近的基因：不管peak落在內(nèi)含子、基因間區(qū)還是其他位置，按照peak相對于轉(zhuǎn)錄起始位點的距離，都能找到一個離它最近的基因【一般做基因表達(dá)調(diào)控的，會關(guān)注promoter區(qū)域，離結(jié)合位點最近的基因更可能被調(diào)控】
這個距離是根據(jù)轉(zhuǎn)錄起始位點來計算，一個基因具有多個轉(zhuǎn)錄本，因此一個基因可能有多個轉(zhuǎn)錄起始位點。注釋的結(jié)果就會看到有一列是轉(zhuǎn)錄本ID
注釋類型三：flankDistance（三列: flank_txIds, flank_geneIds和flank_gene_distances）
指peak上下游某個范圍內(nèi)（比如-5kb《=》5kb范圍內(nèi)）都有什么基因
# 傳個參數(shù)flankDistance
x2 = annotatePeak(f, tssRegion=c(-1000, 1000), TxDb=txdb, addFlankGeneInfo=TRUE, flankDistance=5000)
讓基因名變得友好
上面得到的結(jié)果都是以geneId（Entrez ID）給出，如果想要Symbol名稱，可以再傳參數(shù)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，那么結(jié)果就會是Entrez ID、SYMBOL、GENENAME三列）
按正負(fù)鏈分開注釋
一般ChIPseq數(shù)據(jù)通常情況下是沒有正負(fù)鏈信息的（有特殊的實驗可以有）
但如果要做，可以先給peaks分別賦予正負(fù)鏈的信息，然后指定參數(shù)sameStrand=TRUE 并分別做兩次
這個參數(shù)的意思是：（logical）whether find nearest/overlap gene in the same strand
只注釋基因的上游或下游
提供了ignoreDownstream和 ignoreUpstream，默認(rèn)是FALSE
關(guān)于TxDb的知識
上面一起操作的前提是物種本身有這些注釋信息，而注釋信息主要是用TxDb
同一物種的不同版本TxDb
例如TxDb.Hsapiens.UCSC.hg19.knownGene和TxDb.Hsapiens.UCSC.hg38.knownGene 的注釋結(jié)果是不同的，不能混用。用哪個取決于上游分析比對使用的哪個版本的基因組
不同的版本中基因坐標(biāo)是不一樣的，如果硬要替換，可以使用liftOver將基因組版本坐標(biāo)進(jìn)行轉(zhuǎn)換
支持多少物種？
Bioconductor上有30個左右TxDb，也只能覆蓋一小部分物種（https://bioconductor.org/packages/3.11/data/annotation/），但UCSC和Ensemble的基因組都可以被ChIPseeker支持，因此所有物種都支持
除了基因注釋還能注釋lincRNA
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()函數(shù)依賴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怎么查看呢？
最詳細(xì)的操作在官方文檔：https://bioconductor.org/packages/release/bioc/vignettes/GenomicFeatures/inst/doc/GenomicFeatures.pdf
不管是從Bioconductor下載的還是自己制作的，都是一個GenomicFeatures對象
如果簡單對名稱操作，會返回這個注釋文件的基本信息。要把TxDb當(dāng)成一個數(shù)據(jù)庫來對待，而不是一個簡單的數(shù)據(jù)框或者矩陣。因此它的提取方法也會比較特別
如果想看其中包含的類目，可以用
columns(txdb)如果想指定提取轉(zhuǎn)錄本或外顯子信息，可以：
transcripts(txdb) 或者 exons(txdb)如果想看全部的信息，可以：
AnnotationDbi::select(glycine, columns=columns(glycine), keys=keys(glycine), keytype=c("GENEID"))
需要注意，如果使用這個select的時候，同時加載了tidyverse，那么同名的select就會發(fā)生沖突導(dǎo)致報錯，這時可以用顯式指定的形式來規(guī)范（如下圖）
可視化
peak在整個染色體的分布
見：第3部分=》 使用covplot可視化BED數(shù)據(jù)
peak在某個窗口的結(jié)合譜圖
一般有兩種方式：一是直接使用BED文件，二是一步步手動進(jìn)行
第一種：直接使用BED文件
peakHeatmap(f, weightCol="V5", TxDb=txdb,
upstream=3000, downstream=3000,
color=rainbow(length(f)))
其實看運行日志也能看出來做了什么，首先根據(jù)轉(zhuǎn)錄起始位點指定上下游（也就是熱圖的窗口區(qū)間范圍），然后把peaks比對到這個窗口，并生成矩陣以進(jìn)行可視化
稍微查看一下這個peakHeatmap函數(shù)，就會發(fā)現(xiàn)以上說的幾步：
當(dāng)然，如果是多個文件也是可以的
files=getSampleFiles()
peakHeatmap(files, TxDb=txdb,
upstream=3000, downstream=3000,
color=rainbow(length(files)))
第二種：一步步手動進(jìn)行
如果說第一種提供了一個打包好的計算過程，那么第二種就是把第一種拆分運行
promoter <- getPromoters(TxDb=txdb,
upstream=3000, downstream=3000)
tagMatrix <- getTagMatrix(f,
windows=promoter)
tagHeatmap(tagMatrix, xlim=c(-3000, 3000),
color="red")
看看結(jié)合的強度
第一種：直接使用BED文件
plotAvgProf2(files[[4]], TxDb=txdb,
upstream=3000, downstream=3000,
xlab="Genomic Region (5'->3')",
ylab = "Read Count Frequency",
conf = 0.95, resample = 1000)# 添加置信區(qū)間
第二種：手動進(jìn)行
使用上面的tagMatrix計算結(jié)果
plotAvgProf(tagMatrix, xlim=c(-3000, 3000),
xlab="Genomic Region (5'->3')",
ylab = "Read Count Frequency")
支持多個數(shù)據(jù)比較
tagMatrixList <- lapply(files, getTagMatrix,
windows=promoter)
# 添加置信區(qū)間并分面
plotAvgProf(tagMatrixList, xlim=c(-3000, 3000),
conf=0.95,resample=500, facet="row")
這個結(jié)果和上面peakHeatmap的結(jié)果一致，前3個樣本不是調(diào)控轉(zhuǎn)錄的
除了關(guān)注轉(zhuǎn)錄起始位點（研究轉(zhuǎn)錄調(diào)控），還能看蛋白與外顯子/內(nèi)含子起始位置的結(jié)合譜，使用getBioRegion函數(shù)，可以指定'gene', 'transcript', 'exon', 'intron'
注釋結(jié)果之注釋類型一：genomic annotation
指peak在基因組的位置：落在什么地方，例如外顯子、內(nèi)含子或是5’ /3‘UTR
餅圖
peakAnno <- annotatePeak(files[[4]],
tssRegion=c(-3000, 3000),
TxDb=txdb, annoDb="org.Hs.eg.db")
plotAnnoPie(peakAnno)
柱狀圖
plotAnnoBar(peakAnno)
注意第一個問題：關(guān)于上圖中的各種Features分類
看到這里的分類有下游（Downstream）但沒有上游，這是因為Promoter定義為了轉(zhuǎn)錄起始位點（TSS）的上下游區(qū)域，包含了上游；另外這個下游是是基因間區(qū)的一部分，更確切是指緊接著基因的下游；這里的上游和下游其實都是基因間區(qū)，單獨拿出來是因為和基因直接連接，是很近的區(qū)域=》近端基因間區(qū)
當(dāng)然，基因間區(qū)還包含更遠(yuǎn)的間區(qū)（Distal intergenic）=》遠(yuǎn)端基因間區(qū)
默認(rèn)下游的范圍是3kb，但是可以自己調(diào)整
# 比如調(diào)成500
options(ChIPseeker.downstreamDistance = 500)
還有一個需求就是：自定義分類
# 依然是設(shè)置options，用于總結(jié)結(jié)果
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所在的位置，可能是一個基因的外顯子，同時又是另一個基因的內(nèi)含子。為了解決這個問題，有以下幾種方案：
- 第一種：使用參數(shù)
genomicAnnotationPriority指定優(yōu)先順序
默認(rèn)順序是：Promoter => 5’ UTR => 3’ UTR => Exon => Intron => Downstream => Distal Intergenic
- 第二種：餅圖+韋恩圖
vennpie(peakAnno)
優(yōu)點是：直觀；缺點是：無法顯示全部的信息
- 第三種：UpSetR + vennpie
upsetplot(peakAnno, vennpie=TRUE)
多個文件的區(qū)域注釋
peakAnnoList <- lapply(files, annotatePeak,
TxDb=txdb,tssRegion=c(-3000, 3000))
plotAnnoBar(peakAnnoList)
注釋結(jié)果之注釋類型二：nearest gene annotation
指peak最近的基因：不管peak落在內(nèi)含子、基因間區(qū)還是其他位置，按照peak相對于轉(zhuǎn)錄起始位點的距離，都能找到一個離它最近的基因
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的位置與所有的基因關(guān)聯(lián)起來【包括 host gene (exon/intron), promoter region and flanking gene from intergenomic region】，然后用clusterProfiler拿這些基因跑ORA，做富集
require(clusterProfiler)
bedfile=getSampleFiles()
# 將bed文件讀入（readPeakFile是利用read.delim讀取，然后轉(zhuǎn)為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)
歡迎關(guān)注我們的公眾號～_～
我們是兩個農(nóng)轉(zhuǎn)生信的小碩，打造生信星球，想讓它成為一個不拽術(shù)語、通俗易懂的生信知識平臺。需要幫助或提出意見請后臺留言或發(fā)送郵件到jieandze1314@gmail.com
