使用ChIPpeakAnno進行peak注釋

  • 2019 年 12 月 19 日
  • 筆記

ChIPpeakAnno是一個bioconductor上的R包，針對peak calling之後的下游分析，提供了以下多種功能

  1. 查找與peak區域最相鄰的基因, 也支援自定義查找的特徵，可以是exon,miRNA等
  2. peak相鄰基因的GO富集分析
  3. 提取peak及其周圍區域的序列

在ChIPpeakAnno中，無論是peak區間資訊還是基因組的注釋資訊，都通過toGRanges方法轉化為R語言中的GRanges對象，以peak為例，bed格式的內容如下

通過如下程式碼可以導入該資訊

library(ChIPpeakAnno)  bed <- "peaks.bed"  gr <- toGRanges(bed, format="BED", header=FALSE)

除了BED格式外，該方法也支援導入GTF格式的資訊，只需要修改format參數即可。導入peak資訊和基因組注釋資訊後就可以進行後續分析了。

1. 進行peak之間的overlap分析

當導入了多個樣本的peak資訊時，可以進行venn分析，用法如下

# 導入A樣本的peak  bedA    <- "sampleA_peaks.bed"  sampleA <- toGRanges(bedA, format="BED", header=FALSE)  # 導入B樣本的peak  bedB    <- "sampleB_peaks.bed"  sampleB <- toGRanges(bedB, format="BED", header=FALSE)  # 求交集  ol <- findOverlapsOfPeaks(sampleA, sampleB)  # 繪製venn圖  makeVennDiagram(ol)

結果示意如下

在進行venn分析時，會發現venn圖上的個數加起來並不是輸入的peak區間的總數，在默認

2. 提取peak周圍的序列

用法如下

library(BSgenome.Hsapiens.UCSC.hg19)  seq <- getAllPeakSequence(sampleA, upstream=20, downstream=20, genome=Hsapiens)  write2FASTA(seq, "sampleA.peaks.fa")

3. 進行peak motif分析

提取到peak序列之後，可以進行motif分析，用法如下

# 用1號染色體的鹼基分布當做背景  freqs <- oligoFrequency(Hsapiens$chr1, MarkovOrder=3)  # oligoLength規定了motif的長度  os <- oligoSummary(seq, oligoLength=6, MarkovOrder=3,                     quickMotif=TRUE, freqs=freqs)  zscore <- sort(os$zscore)  # 繪製所有6個鹼基組合的頻率分布圖  h <- hist(zscore, breaks=100, xlim=c(-50, 50), main="Histogram of Z-score")  # 頻率最大的鹼基組合即為motif的結果  text(zscore[length(zscore)], max(h$counts)/10,       labels=names(zscore[length(zscore)]), adj=1)

結果示意如下

還可以通過motifStack這個R包繪製motif的sequence logo, 用法如下

library(motifStack)  pfms <- mapply(function(.ele, id)      new("pfm", mat=.ele, name=paste("SAMPLE motif", id)),      os$motifs, 1:length(os$motifs))  motifStack(pfms[[1]])

輸出結果示意如下

4. 進行peak注釋

首先是peak在基因組各個特徵區間的分布比例，用法如下

library(TxDb.Hsapiens.UCSC.hg19.knownGene)  aCR<-assignChromosomeRegion(sampleA, nucleotideLevel=FALSE,                             precedence=c("Promoters", "immediateDownstream",                                           "fiveUTRs", "threeUTRs",                                           "Exons", "Introns"),                             TxDb=TxDb.Hsapiens.UCSC.hg19.knownGene)  barplot(aCR$percentage, las=3)

輸出結果如下所示

然後進行peak關聯基因的注釋，用法如下

# 準備基因組注釋資訊  library(EnsDb.Hsapiens.v75)  annoData <- toGRanges(EnsDb.Hsapiens.v75, feature="gene")  # 進行  overlaps.anno <- annotatePeakInBatch(sampleA,                                       AnnotationData=annoData,                                       output="nearestLocation"  )  library(org.Hs.eg.db)  overlaps.anno <- addGeneIDs(overlaps.anno,                              "org.Hs.eg.db",                              IDs2Add = "entrez_id")  pie1(table(overlaps.anno$insideFeature))

輸出結果示意如下

在使用annotatePeakInBatch進行注釋時，默認查找距離peak最近的基因，也可以修改output的值，overlapping代表與peak區域存在overlap的基因，設置成這個值之後就會將與peak區間存在overlap的基因作為關聯基因了，此外還有多種取值，適用不同條件，具體可以參考函數的幫助文檔。

5. 進行peak關聯基因的富集分析

進行完基因注釋之，得到peak關聯的基因，就可以進行後續的功能富集分析，用法如下

over <- getEnrichedGO(overlaps.anno, orgAnn="org.Hs.eg.db",                       maxP=.05, minGOterm=10,                       multiAdjMethod="BH", condense=TRUE)

ChIPpeakAnno提供了一條完整的peak下游分析功能，包括基因注釋，富集分析，motif分析等等，是一個非常強大的工具，以上只是基本用法，更多用法和細節請參考官方文檔。

·end·

Tags:

VirMach 便宜 VPS