ChIPseeker包南方醫(yī)科大學(xué)Y叔大牛寫(xiě)的許多有名的生信R包之一赋元，其最初設(shè)計(jì)用于chip-seq的macs peak calling結(jié)果分析以及可視化，后來(lái)逐漸也適用于相關(guān)的peak分析。
參考鏈接：https://www.bioconductor.org/packages/release/bioc/vignettes/ChIPseeker/inst/doc/ChIPseeker.html掖鱼；
以及Y叔自己的微信公眾號(hào)教程：https://mp.weixin.qq.com/mp/appmsgalbum?__biz=MzI5NjUyNzkxMg==&action=getalbum&album_id=1300625300497268737&scene=173&from_msgid=2247488238&from_itemidx=1&count=3#wechat_redirect
1、關(guān)于ChIP-seq
- 之前在學(xué)習(xí)macs文章時(shí)，有了解過(guò)疆前；這里再簡(jiǎn)單學(xué)習(xí)一下。
- 如下圖聘萨，DNA上的蛋白結(jié)合位點(diǎn)往往是基因表達(dá)調(diào)控的關(guān)鍵位置竹椒，ChIP技術(shù)就是針對(duì)性的挑選出這些位置。
- DNA和蛋白質(zhì)交聯(lián)(cross-linking)米辐，超聲(sonication)將染色體隨機(jī)切割胸完，利用抗原抗體的特異性識(shí)別(IP)书释，把目標(biāo)蛋白相結(jié)合的DNA片段沉淀下來(lái)，反交聯(lián)釋放DNA片段赊窥，最后是測(cè)序(sequencing)爆惧。
-
MACS軟件通過(guò)一定原理算法，對(duì)測(cè)序比對(duì)結(jié)果識(shí)別出有意義的peak锨能。ChIPseeker包就是銜接這一步之后開(kāi)始做的扯再。
ChIP
2、準(zhǔn)備工作 preparation
- 因?yàn)閙acs結(jié)果為bed輸出格式址遇，所以需要了解bed熄阻，即bedtools軟件
- 了解GRanges、TxDb這兩種常見(jiàn)的生信基礎(chǔ)數(shù)據(jù)對(duì)象
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
- 安裝R包倔约，找到示例數(shù)據(jù)
BiocManager::install("ChIPseeker")
library(ChIPseeker)
files <- getSampleFiles()
print(files)
#bed轉(zhuǎn)為Granges對(duì)象
peak <- readPeakFile(files[[4]])
peak
-
如下圖秃殉，即為ChIPseeker包分析所需的peak GRange對(duì)象。
分割線(xiàn)左邊三列分別為所在染色體信息浸剩，起止位點(diǎn)钾军，正負(fù)鏈情況；
右邊兩列分別為peak name與score（我認(rèn)為可以理解與reads數(shù)正相關(guān)）
peak
3乒省、ChIPseeker基礎(chǔ)peak可視化
3.1巧颈、概況covplot()
觀察所有peak在染色體的分布、表達(dá)情況
#依據(jù)第五列score袖扛，表明峰的高低情況
covplot(peak, weightCol="V5")
#篩選指定染色體的指定區(qū)域的分布情況
covplot(peak, weightCol="V5", chrs=c("chr17", "chr18"), xlim=c(4.5e7, 5e7))
3.2砸泛、針對(duì)某一feature的分布情況
- heatmap
常見(jiàn)的分析是觀察不同peak分布在TSS的promoter區(qū)域情況
#自己定義promoter區(qū)域，上下游3000bp
promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
#不理解這個(gè)函數(shù)也沒(méi)關(guān)系蛆封，是為下一步做熱圖提供matrix
tagMatrix <- getTagMatrix(peak, windows=promoter)
tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red")
如下圖結(jié)果唇礁，每一行代表一個(gè)promoter區(qū)域，紅線(xiàn)的即為peak分布
#一鍵繪圖惨篱，效果同上
peakHeatmap(peak, TxDb=txdb, upstream=3000, downstream=3000, color="red")
- 峰圖
上面的熱圖是描繪了所有的promoter情況盏筐，可以繪制一個(gè)峰圖描述所有分布的平均情況。
plotAvgProf(tagMatrix, xlim=c(-3000, 3000),
xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency")
#加一個(gè)置信區(qū)間
plotAvgProf(tagMatrix, xlim=c(-3000, 3000), conf = 0.95, resample = 1000)
we developed getBioRegion function to support centering all peaks to the start region of Exon/Intron. Users can also create heatmap or average profile of ChIP peaks binding to these regions.
4砸讳、ChIPseeker peak annotation
4.1 what's peak annotation
-
簡(jiǎn)單理解peak 注釋就是peak落在染色體的哪一個(gè)位置上琢融。常見(jiàn)的基因結(jié)構(gòu)組成如下圖所示。
basic structure of gene - 此外ChIPseeker的peak注釋時(shí)還提供另外一種注釋方法，具體在注釋結(jié)果時(shí)再具體了解（nearest gene annotation）。
4.2 annotatePeak()
（1）just do it
- ChIPseeker包主要用
annotatePeak()注釋peak官地。需要提供兩個(gè)文件：一是peak文件康嘉，可以是bed或者Granges；另一個(gè)是對(duì)應(yīng)物種的TxDb對(duì)象（提供原始注釋信息） - 此外promoter的區(qū)間可以自己定義，默認(rèn)設(shè)置為T(mén)SS上下游3k區(qū)域
peak
#共計(jì)1331個(gè)peak
txdb
peakAnno <- annotatePeak(files[[4]], tssRegion=c(-3000, 3000), TxDb=txdb)
peakAnno
如下圖香浩，如果在R里直接觀察結(jié)果德挣，它會(huì)告訴我們ChIPseq的位點(diǎn)落在基因組上什么樣的區(qū)域平绩，分布情況如何圈匆。（即第一種注釋方法genomic annotation）
在注釋時(shí)，有的peak可能同時(shí)落在兩個(gè)或者更多的gene feature里（例如是一個(gè)基因的外顯子而同時(shí)又是另一個(gè)基因的內(nèi)含子）捏雌，但只能注釋其中一個(gè)跃赚。默認(rèn)按照Promoter、5’ UTR腹忽、3’ UTR来累、Exon、Intron窘奏、Downstream嘹锁、Intergenic順序先后注釋。
- 一般會(huì)將上述的結(jié)果輸出為GRanges格式着裹、或者data.frame格式领猾；便于查看，同時(shí)也能了解到
annotatePeak第二種nearest gene annotation結(jié)果骇扇。
class(peakAnno)
peakAnno.df <- as.data.frame(peakAnno)
peakAnno.gr <- as.GRanges(peakAnno)
head(peakAnno.gr, 3)
如下圖摔竿，右上角為genomic annotation結(jié)果、下面為nearest gene annotation結(jié)果少孝。
- nearest gene annotation最近基因注釋?zhuān)菏莗eak相對(duì)于轉(zhuǎn)錄起始位點(diǎn)的距離继低，不管這個(gè)peak是落在內(nèi)含子或者別的什么位置上，即使它落在基因間區(qū)上稍走，我都能夠找到一個(gè)離它最近的基因（即使它可能非常遠(yuǎn)）袁翁。
- 如果peak和TSS有overlap，genomic annotation就是promoter婿脸，距離就是0粱胜，而最近基因也是同一個(gè)，所以在這種情況下狐树，兩種注釋都指向同一個(gè)基因焙压。
-
最近基因的注釋信息雖然是以基因?yàn)閱挝唤o出，但我們針對(duì)的是轉(zhuǎn)錄起始位點(diǎn)來(lái)計(jì)算距離抑钟，針對(duì)于不同的轉(zhuǎn)錄本涯曲，一個(gè)基因可能有多個(gè)轉(zhuǎn)錄起始位點(diǎn)，所以注釋是在轉(zhuǎn)錄本的水平上進(jìn)行的在塔，我們可以看到輸出有一列是transcriptId.
head(peakAnno.gr, 3)
另外一種思路：注意上述nearest gene annotation默認(rèn)找的是最近的TSS幻件，即first anno與second anno對(duì)應(yīng)的可以不是同一個(gè)基因。如果我想說(shuō)只要和基因有overlap就是最近基因心俗，那么這兩種注釋的基因應(yīng)該是一致的，只需把
overlap="TSS"(default)設(shè)置為overlap="all"
5、ChIPseeker基于注釋的peak可視化
（1）genomic annotation可視化
- 餅圖或柱狀圖可視化組成比例
plotAnnoPie(peakAnno)
plotAnnoBar(peakAnno)
- 考慮到注釋到多種feature的可能
vennpie(peakAnno)
upsetplot(peakAnno, vennpie=TRUE)
如下圖可以清楚地看到絕大多數(shù)的peak同時(shí)落到了多種feature里
（2）nearest gene annotation結(jié)果可視化
- 可視化the distance from the peak (binding site) to the TSS of the nearest gene
-
plotDistToTSScan calculate the percentage of binding sites upstream and downstream from the TSS of the nearest genes, and visualize the distribution.
plotDistToTSS(peakAnno,
title="Distribution of transcription factor-binding loci\nrelative to TSS")
ChIPseeker包暫時(shí)先學(xué)習(xí)到這里城榛，還有很多深入的功能揪利，比如富集分析等，之后有機(jī)會(huì)再學(xué)習(xí)狠持。
感覺(jué)到國(guó)人寫(xiě)的R包疟位，然后看中文的原版說(shuō)明書(shū)還是比較輕松的~~
