一、 基本介紹

ChIPseeker可以用來(lái)對(duì)ChIP-seq數(shù)據(jù)進(jìn)行注釋與可視化。

二、 使用方法

(1) 安裝并載入R包

ls() #列出當(dāng)前工作空間中的對(duì)象
rm(list=ls()) #移除全部對(duì)象

options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")) 
options(BioC_mirror="http://mirrors.tuna.tsinghua.edu.cn/bioconductor")

if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

#BiocManager::install("ChIPpeakAnno")
BiocManager::install("ChIPseeker")
BiocManager::install("TxDb.Dmelanogaster.UCSC.dm6.ensGene")
BiocManager::install("org.Dm.eg.db")

#library(ChIPpeakAnno)
library(ChIPseeker)
library(TxDb.Dmelanogaster.UCSC.dm6.ensGene)
library(org.Dm.eg.db)
library(ggplot2)
library(clusterProfiler)

(2) bed文件處理

$ ls 11_intersect/*.bed | while read id ; do (cat $id | awk '{if($1=="2L"||$1=="2R"||$1=="3L"||$1=="3R"||$1=="4"||$1=="X"||$1=="Y") {print "chr"$1"\t"$2"\t"$3"\t"$4}}' > 12_ChIPseeker/bed/$(basename $id)) ; done  &

(3) 峰在整個(gè)染色體的分布

peak1 <- readPeakFile("12_ChIPseeker/bed/WTF_intersect.bed")  #GRanges object
peak2 <- readPeakFile("12_ChIPseeker/bed/WTM_intersect.bed")

peaks <- list(WTF=peak1, WTM=peak2)

covplot(peak1)
covplot(peaks) + facet_grid(chr ~ .id)
covplot(peaks, chrs=c("chrX"), xlim=c(1e7,1.1e7)) + facet_grid(chr ~ .id)

(4) 峰注釋

txdb = TxDb.Dmelanogaster.UCSC.dm6.ensGene

peakAnno <- annotatePeak(peak1,
                         tssRegion = c(-3000, 3000),  #啟動(dòng)子區(qū)域
                         genomicAnnotationPriority = c("Promoter", "5UTR", "3UTR", "Exon", "Intron","Downstream", "Intergenic"),
                         flankDistance = 5000,
                         #sameStrand = T,
                         TxDb = txdb,
                         annoDb = "org.Dm.eg.db") 

peakAnno
peakAnno_df <- as.data.frame(peakAnno)
head(peakAnno_df)

write.csv(peakAnno_df,file="12_ChIPseeker/WTF_peakAnno.csv",row.names=F)

注釋結(jié)果：

seqnames  start  end  width  strand  V4  annotation  geneChr  geneStart  geneEnd geneLength  geneStrand  geneId  transcriptId  distanceToTSS  ENTREZID  SYMBOL  GENENAME
chrX  138244  138423  180  *  wt_msl2_peak_1  Exon (FBtr0340053/FBgn0025836, exon 9 of 14)  6  140722  142197  1476  2  FBgn0267029  FBtr0345995  3774  19834782  lncRNA:CR45473  long non-coding RNA:CR45473

這里注釋有兩種，一種是genomic annotation (也就是annotation這一列），還有就是nearest gene annotation（也就是多出來(lái)的其它列）。genomic annotation注釋的是peak的位置，它落在什么地方了，可以是UTR，可以是內(nèi)含子或者外顯子。而最近基因是peak相對(duì)于轉(zhuǎn)錄起始位點(diǎn)的距離，不管這個(gè)peak是落在內(nèi)含子或者別的什么位置上，即使它落在基因間區(qū)上，我都能夠找到一個(gè)離它最近的基因（即使它可能非常遠(yuǎn)）。

針對(duì)不同的問(wèn)題，這兩種注釋的策略是不一樣的。第一種策略peak所在位置，可能就是調(diào)控的根本，比如你要做可變剪切的，最近基因的注釋顯然不是你關(guān)注的點(diǎn)。而做基因表達(dá)調(diào)控的，當(dāng)然promoter區(qū)域是重點(diǎn)，離結(jié)合位點(diǎn)最近的基因更有可能被調(diào)控。最近基因的注釋信息雖然是以基因?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。

兩種注釋有時(shí)候還不夠，我想看peak上下游某個(gè)范圍內(nèi)（比如說(shuō)-5k到5k的距離）都有什么基因，annotatePeak也可以做到。你只要傳個(gè)參數(shù)說(shuō)你要這個(gè)信息，還有什么區(qū)間內(nèi)，就可以了。輸出中多三列：flank_txIds、flank_geneIds和flank_gene_distances，在指定范圍內(nèi)所有的基因都被列出。

最近基因給出來(lái)了，但都是各種人類(lèi)不友好的ID，只需要給annotatePeak傳入annoDb參數(shù)就行了。如果你的TxDb的基因ID是Entrez，它會(huì)轉(zhuǎn)成ENSEMBL，反之亦然，當(dāng)然不管是那一種，都會(huì)給出SYMBOL，還有描述性的GENENAME。

The annotation column annotates genomic features of the peak, that is whether peak is located in Promoter, Exon, UTR, Intron or Intergenic. If the peak is annotated by Exon or Intron, more detail information will be provided. For instance, Exon (38885 exon 3 of 11) indicates that the peak is located at the 3th exon of gene 38885 (entrezgene ID) which contain 11 exons in total.

一個(gè)peak所在的位置，可能是一個(gè)基因的外顯子，同時(shí)又是另一個(gè)基因的內(nèi)含子。使用參數(shù)genomicAnnotationPriority指定優(yōu)先順序，默認(rèn)順序是：Promoter => 5’ UTR => 3’ UTR => Exon => Intron => Downstream => Distal Intergenic。

Downstream is defined as the downstream of gene end. Promoter定義為轉(zhuǎn)錄起始位點(diǎn)（TSS）的上下游區(qū)域。另外下游是基因間區(qū)的一部分，更確切是指緊接著基因的下游；這里的上游和下游其實(shí)都是基因間區(qū)，單獨(dú)拿出來(lái)是因?yàn)楹突蛑苯舆B接，是很近的區(qū)域，即近端基因間區(qū)。當(dāng)然，基因間區(qū)還包含更遠(yuǎn)的間區(qū)（Distal intergenic），即遠(yuǎn)端基因間區(qū)。

(5) 比對(duì)到TSS附近的峰分布

# plot the heatmap of peaks align to flank sequences of TSS.
#和deeptools的reference-point mode差不多，但是畫(huà)的圖有點(diǎn)丑
peakHeatmap(peak1,
            TxDb=txdb, 
            upstream=3000, downstream=3000,  #指定轉(zhuǎn)錄起始位點(diǎn)上下游
            color=rainbow(length(peak1)))

peakHeatmap(peaks,  #list
            TxDb=txdb, 
            upstream=3000, downstream=3000, 
            color=rainbow(length(peaks)))

# plot the profile of peaks that align to flank sequences of TSS.
plotAvgProf2(peak1, 
             TxDb=txdb, 
             upstream=3000, downstream=3000,
             xlab="Genomic Region (5'->3')", 
             ylab = "Read Count Frequency",
             conf = 0.95, resample = 1000)

plotAvgProf2(peaks, 
             TxDb=txdb, 
             upstream=3000, downstream=3000,
             xlab="Genomic Region (5'->3')", 
             ylab = "Read Count Frequency",
             conf = 0.95, resample = 1000) + facet_grid(. ~ .id)

(6) 峰在基因組的位置

peakAnno <- annotatePeak(peak1, 
                         tssRegion=c(-3000, 3000),
                         TxDb=txdb,
                         annoDb="org.Dm.eg.db")
plotAnnoPie(peakAnno)

peakAnnoList <- lapply(peaks, annotatePeak, 
                       TxDb=txdb,
                       tssRegion=c(-3000, 3000))
plotAnnoBar(peakAnnoList)

(7) 相對(duì)于轉(zhuǎn)錄起始點(diǎn)的距離

相對(duì)于轉(zhuǎn)錄起始位點(diǎn)的距離找最近的基因。不管peak落在內(nèi)含子、基因間區(qū)還是其他位置，按照peak相對(duì)于轉(zhuǎn)錄起始位點(diǎn)的距離，都能找到一個(gè)離它最近的基因。

plotDistToTSS(peakAnno,
              title="Distribution of transcription factor-binding loci\nrelative to TSS")
plotDistToTSS(peakAnnoList,
              title="Distribution of transcription factor-binding loci\nrelative to TSS")

(8) 距離最近的基因的交集和功能

genes <- lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)
names(genes)