ChIP-seq試驗(yàn)普遍用于探究蛋白質(zhì)與核酸的結(jié)合作用。例如最常見的，期望獲得某個(gè)轉(zhuǎn)錄因子的功能，查看它結(jié)合到哪些基因啟動(dòng)子區(qū)域并調(diào)控它們的轉(zhuǎn)錄，以及更高級(jí)的分析識(shí)別該轉(zhuǎn)錄因子的motif位點(diǎn)等，就可通過ChIP-seq來實(shí)現(xiàn)。

作了ChIP-seq試驗(yàn)，經(jīng)過初步下機(jī)處理并比對(duì)參考基因組之后，獲得了目標(biāo)蛋白在基因組區(qū)域的結(jié)合峰位置。接下來，就需要對(duì)峰位置進(jìn)行注釋，獲得這些峰結(jié)合在了哪些基因上，位于基因的上下游還是內(nèi)部，以便了解目標(biāo)蛋白的功能。

那么，怎樣對(duì)ChIP-seq峰的結(jié)合位置進(jìn)行注釋呢？本篇簡(jiǎn)介一個(gè)非常好用的R包，ChIPseeker，能夠快速高效地獲得預(yù)期結(jié)果。

輸入數(shù)據(jù)準(zhǔn)備

輸入數(shù)據(jù)就是記錄ChIP-seq峰位置的bed文件，這是ChIP-seq的常見數(shù)據(jù)格式，測(cè)序公司都會(huì)給提供的。如下示例，是在人類細(xì)胞中進(jìn)行的ChIP-seq試驗(yàn)，經(jīng)下機(jī)處理并和人參考基因組hg38比對(duì)后獲得的bed文件。

bed文件中包含5列信息，以tab鍵分隔。第1列為峰值所處的染色體id，第2-3列為峰值在該染色體中的堿基位置，第4列為峰的id，第5列為信號(hào)強(qiáng)度。用Excel打開bed文件就是長(zhǎng)這樣。

bed文件內(nèi)容格式

ChIPseeker安裝

ChIPseeker包可以使用Bioconductor進(jìn)行安裝。

#安裝
BiocManager::install("ChIPseeker")

#加載
library(ChIPseeker)

準(zhǔn)備或構(gòu)建基因組注釋庫

為了對(duì)ChIP-seq的峰位置進(jìn)行注釋，查看蛋白主要結(jié)合到基因組的哪些區(qū)域，首先需要準(zhǔn)備對(duì)應(yīng)物種的基因組注釋庫。

例如，上述給定的示例bed文件來源物種是人，因此就需要準(zhǔn)備人類參考基因組注釋庫。有些R包，例如org.Hs.eg.db中提供了已經(jīng)構(gòu)建好的人類參考基因組注釋庫，可以直接調(diào)用。此外，也可以通過GenomicFeatures包，讀取基因組注釋文件gff3本地構(gòu)建。

推薦通過GenomicFeatures包本地構(gòu)建注釋庫，因?yàn)榻^大多數(shù)物種都是沒有現(xiàn)成的庫可用的，方法靈活。而且，像人類參考基因組也是有很多版本的，可能已知的庫和自己使用的基因組版本不匹配，用錯(cuò)了就尷尬了，為了避免萬一還是推薦從頭構(gòu)建一下。

#可以使用現(xiàn)有的庫，如org.Hs.eg.db包的“org.Hs.eg.db”
library(org.Hs.eg.db)

#或者自己構(gòu)建一個(gè)，指定gff3注釋文件即可
library(GenomicFeatures)
spompe <- makeTxDbFromGFF('Homo_sapiens.GRCh38.98.gff3')

注釋ChIP-seq的峰位置（單個(gè)bed文件的注釋）

好了，接下來就展示ChIPseeker注釋ChIP-seq峰的全過程。

已知上述ChIP-seq的bed文件獲取過程中，比對(duì)的參考基因組來自hg38，因此首先通過本地準(zhǔn)備的人類參考基因組hg38的gff3注釋文件，構(gòu)建本地注釋庫，隨后讀取ChIP-seq的bed文件進(jìn)行注釋。

#推薦手動(dòng)構(gòu)建TxDb注釋文件
spompe <- makeTxDbFromGFF('Homo_sapiens.GRCh38.98.gff3')

#讀取chipseq峰的bed文件
peak <- readPeakFile('ChIP-seq.bed')

#注釋，TSS的范圍可自定義
peakAnno <- annotatePeak(peak, tssRegion = c(-3000, 3000), TxDb = spompe)

#輸出結(jié)果
write.table(peakAnno, file = 'peak.txt',sep = '\t', quote = FALSE, row.names = FALSE)

ChIPseeker輸出表格內(nèi)容展示

結(jié)果中，注釋了ChIP-seq峰結(jié)合的染色體位置、區(qū)段、基因名稱、基因位置等信息。

對(duì)于基因組區(qū)域，基因啟動(dòng)子區(qū)、外顯子或內(nèi)含子區(qū)均有。后續(xù)就可以根據(jù)這些信息作一些篩選，例如ChIP-seq試驗(yàn)使用的蛋白如果是轉(zhuǎn)錄因子，那么通常就是在啟動(dòng)子區(qū)發(fā)揮作用，只關(guān)注啟動(dòng)子區(qū)的結(jié)合位點(diǎn)即可。


同時(shí)，ChIPseeker包中自帶了一些可視化函數(shù)，可以進(jìn)行一些簡(jiǎn)單的統(tǒng)計(jì)幫助我們了解數(shù)據(jù)分布，例如ChIP-seq峰所在基因不同位置的數(shù)量占比等。

plotAnnoBar(peakAnno)
vennpie(peakAnno)
plotAnnoPie(peakAnno)
plotDistToTSS(peakAnno)

一些簡(jiǎn)單統(tǒng)計(jì)圖

注釋ChIP-seq的峰位置（多個(gè)bed文件的批量處理）

很多情況下，ChIP-seq試驗(yàn)不止作了一次，可能存在多個(gè)bed文件有待處理。如果這時(shí)，一個(gè)個(gè)分別注釋就會(huì)略顯繁瑣。好在ChIPseeker也提供了批量注釋的方法，如果有多個(gè)bed文件，可以放在一個(gè)list里面統(tǒng)一執(zhí)行，更加高效。

參考以下示例，兩個(gè)bed文件的注釋，參考基因組仍然是人類hg38，因此繼續(xù)使用上文構(gòu)建好的庫執(zhí)行。

#讀取chipseq峰的bed文件
peak1 <- readPeakFile('ChIP-seq1.bed')
peak2 <- readPeakFile('ChIP-seq2.bed')
peaks <- list(peak1 = peak1, peak2 = peak2)

#注釋
peakAnnoList <- lapply(peaks, annotatePeak, TxDb = spompe, tssRegion = c(-3000, 3000), addFlankGeneInfo = TRUE, flankDistance = 5000)

#分別輸出結(jié)果
write.table(peakAnnoList[1], file = 'peak1.txt',sep = '\t', quote = FALSE, row.names = FALSE)
write.table(peakAnnoList[2], file = 'peak2.txt',sep = '\t', quote = FALSE, row.names = FALSE)

#可視化，兩個(gè)可以放在一起比較
plotAnnoBar(peakAnnoList)
vennpie(peakAnnoList[[1]])
plotAnnoPie(peakAnnoList[[1]])
plotDistToTSS(peakAnnoList)

多組比較統(tǒng)計(jì)圖

輸出表格略，內(nèi)容和單個(gè)運(yùn)行的結(jié)果是一致的。對(duì)于輸出結(jié)果圖，則是將兩組放在一起比較的圖。