有一天我的師弟提了一個需求：
對於ChIP的下遊分析，我很喜歡DiffBind做差異分析然後用ChIPseeker做注釋這一套流程（因為ChIPseeker的輸入格式是GRange格式，而DiffBind的dba.report輸出也恰好是GRange格式，兩者可以無縫銜接）。我之前的套路，一直都是無腦輸出所有的DiffBind結果，然後放入ChIPSeeker裡面做注釋，完了之後轉成data.frame，然後對data.frame做subset，做後續的GO分析（）。但今天突然想對ChIPseeker的注釋結果直接做subset，然後分別對上下調的結果畫plotAnnoPie等ChipSeeker內置的畫圖函數。但發現subset無法對ChipSeeker annotatePeak函數的輸出結果做篩選。
當然，對於DiffBind的結果用subset，分別提取上下調，然後再注釋也是可以的，不過這樣感覺更麻煩。
面對需求，如果無法解決提出需求的人，那麼就只能寫代碼實現這個需求了。於是經過一個下午的努力，我給ChIPseeker加上了subset的功能。
先來介紹下如何使用。我們需要用devtools安裝最新的ChIPseeker
BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")install.packages("ggupset")install.packages("ggimage")devtools::install_github("YuLab-SMU/ChIPseeker")library(ChIPseeker)我們以測試數據集來演示
library(TxDb.Hsapiens.UCSC.hg19.knownGene)txdb <- TxDb.Hsapiens.UCSC.hg19.knownGenepeakfile <- system.file("extdata", "sample_peaks.txt", package="ChIPseeker")peakAnno <- annotatePeak(peakfile, tssRegion=c(-3000, 3000), TxDb=txdb)peakAnno此時會輸出peakAnno裡的描述統計信息
Annotated peaks generated by ChIPseeker150/150 peaks were annotatedGenomic Annotation Summary: Feature Frequency8 Promoter (<=1kb) 4.6666679 Promoter (1-2kb) 4.00000010 Promoter (2-3kb) 2.00000035' UTR 1.3333332 3' UTR 2.6666676 Other Exon 8.0000001 1st Intron 13.3333337 Other Intron 32.0000005 Downstream (<=300) 1.3333334 Distal Intergenic 30.666667peakAnno可以直接用來畫圖
plotAnnoPie(peakAnno)upsetplot(peakAnno)
Pie
upset假如你突發奇想，我能不能就看看某一條染色體的作圖結果，或者對於MACS2的結果，想根據FDR篩選下peak，然後看下peak的變化。
在之前我們無法直接對peakAnno背後的csAnno對象進行操作，所以需要先對輸入的GRanges進行過濾然後再用annotatePeak進行注釋，然後畫圖。
但是現在我們有了subset, 這一切就變得稍微簡單起來了
我們可以按照染色體編號進行篩選
peakAnno_subset <- subset(peakAnno, seqnames == "chr1")upsetplot(peakAnno_subset, vennpie = TRUE) 可以按照peak寬度進行篩選
peakAnno_subset <- subset(peakAnno, length > 500)upsetplot(peakAnno_subset, vennpie = TRUE)
filter如果篩選之後沒有任何結果，那麼就會報錯
> subset(peakAnno, FDR... < .05)Error in names(x) <- value : 'names' attribute [2] must be the same length as the vector [1]那麼問題就是，這些篩選條件的名字怎麼獲知。
只要將peakAnno轉成GRanges格式，所以可以看到的列都是你可以篩選的標準。
gr <- as.GRanges(peakAnno)head(gr)GRanges object with 2 ranges and 15 metadata columns: seqnames ranges strand | length summit tags <Rle> <IRanges> <Rle> | <integer> <integer> <integer> [1] chr10 105137981-105138593 * | 6141747 [2] chr10 42644417-42645383 * | 96866927 X.10.log10.pvalue. fold_enrichment FDR... <numeric> <numeric> <numeric> [1] 52.815.32 <NA> [2] 77.159.130.79 annotation geneChr geneStart geneEnd <character> <integer> <integer> <integer> [1] Exon (uc001kwv.3/6877, exon 3 of 11) 10105127724105148822 [2] Distal Intergenic 104282731442863493 geneLength geneStrand geneId transcriptId distanceToTSS <integer> <integer> <character> <character> <numeric> [1] 2109916877 uc010qqq.2 10257 [2] 361802441666 uc010qey.2 218110給ChIPseeker增加subset的功能，應該是我第一次在GitHub進行pull requests操作。寫這個函數，需要先去理解csAnno對象到底是什麼，以及如何保證csAnno這個對象在操作前後不會發生變化。
其實這個函數挺簡單的，就是下面幾行
##' @importFrom S4Vectors subset##' @importFrom BiocGenerics start##' @importFrom BiocGenerics end##' @method subset csAnno##' @exportsubset.csAnno <- function(x, ... ){ index <- paste(seqnames(x@anno),start(x@anno),end(x@anno), sep = "_")# subset the GRanges x@anno <- subset(x@anno, ...) index2 <- paste(seqnames(x@anno),start(x@anno),end(x@anno), sep = "_")# the tssRgion, level, hsaGenomicAnnotation keep unchanged# change the detailGenomicAnnotation x@detailGenomicAnnotation <- x@detailGenomicAnnotation[index %in% index2,]# change the annotation stat x@annoStat <- getGenomicAnnoStat(x@anno)# change peak number x@peakNum <- length(x@anno)return(x)}因為有很多功能是ChIPseeker已經有了的，比如說getGenomicAnnoStat, 還有一些是S4Vectors和BiocGenerics包提供的，比如說subset, start和end。