ChIPseeker对ChIP-seq数据进行注释与可视化
ChIPseeker对ChIP-seq数据进行注释与可视化
这里我们主要介绍ChIPseeker包用于ChIP-seq数据的注释与可视化，主要分为以下几个部分。
01
数据准备
在用ChIPseeker包进行注释前，需要准备两个文件：
1
注释peaks的文件
该文件需要满足BED格式。BED格式文件至少得有chrom（染色体名字），chromStart（染色体起始位点）和chromEnd（染色体终止位点），其它信息如name，score，strand等可有可无。一般情况而言，可直接用做peak calling的MACS输出文件（以_peaks.bed结尾文件）。
2
有注释信息的TxDb对象
Bioconductor包提供了30个TxDb包，包含了很多物种，如人，老鼠等。当所研究的物种没有已有的TxDb时，可利用GenomicFeatures包进行制作：
makeTxDbFromUCSC：通过UCSC在线制作TxDb
makeTxDbFromBiomart：通过Ensembl在线制作TxDb
makeTxDbFromGRanges：通过GRanges对象制作TxDb
makeTxDbFromGFF：通过解析GFF文件制作TxDb如有对应物种的gff文件，用makeTxDbFromGFF来制作TxDb:
require(GenomicFeatures)
tomato <- makeTxDbFromGFF("Slycopersicum_gene.gff3")02
数据可视化
查看peak在全基因组的分布情况，用covplot可视化BED格式文件：
##安装并加载ChIPseeker包
source ("https://bioconductor.org/biocLite.R")
biocLite("ChIPseeker")
library(ChIPseeker)
library(ggplot2)
files <- getSampleFiles()> files
$`ARmo_0M`
[1] "C:/Users/huangjiao/Documents/R/win-library/3.5/ChIPseeker/extdata/GEO_sample_data/GSM1174480_ARmo_0M_peaks.bed.gz"
$ARmo_1nM
[1] "C:/Users/huangjiao/Documents/R/win-library/3.5/ChIPseeker/extdata/GEO_sample_data/GSM1174481_ARmo_1nM_peaks.bed.gz"
$ARmo_100nM
[1] "C:/Users/huangjiao/Documents/R/win-library/3.5/ChIPseeker/extdata/GEO_sample_data/GSM1174482_ARmo_100nM_peaks.bed.gz"
$CBX6_BF
[1] "C:/Users/huangjiao/Documents/R/win-library/3.5/ChIPseeker/extdata/GEO_sample_data/GSM1295076_CBX6_BF_ChipSeq_mergedReps_peaks.bed.gz"
$CBX7_BF
[1] "C:/Users/huangjiao/Documents/R/win-library/3.5/ChIPseeker/extdata/GEO_sample_data/GSM1295077_CBX7_BF_ChipSeq_mergedReps_peaks.bed.gz"files包括5个BED格式文件，对应不同样本数据。
##可视化第4个BED格式文件
covplot(files[[4]])
同时可视化多个BED格式文件：
require(GenomicFeatures)
peak=GenomicRanges::GRangesList(CBX6=readPeakFile(files[[4]]),CBX7=readPeakFile(files[[5]]))
covplot(peak, weightCol="V5") + facet_grid(chr ~ .id)
只关注其中的几条染色体：
covplot(peak, weightCol="V5", chrs=c("chr17", "chr18"), xlim=c(4e7, 5e7)) + facet_grid(chr ~ .id)
03
Peak可视化
tagHeatmap函数查看TSS（转录起始位点）上下游3kb区域peak的分布情况：
##加载对应的TxDb包
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
peak <-readPeakFile(files[[4]])##getPromoters函数准备启动子窗口
promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)##getTagMatrix函数把peaks比对到启动子窗口上
tagMatrix <- getTagMatrix(peak, windows=promoter)
tagHeatmap(tagMatrix, xlim =c(-3000, 3000), color="red")
用peakheatmap函数查看多个样本TSS上下游3kb区域内peak的分布情况：
peakHeatmap(files, TxDb=txdb, upstream=3000, downstream=3000, color=rainbow(length(files)))
这里我们可以看出ARmo_0M，ARmo_1nM和ARmo_100nM三个样本的结合位点不在TSS附近，而CBX6_BF和CBX7_BF两个样本结合位点在TSS附近。
用plotAvgProf函数查看结合的强度：
plotAvgProf(tagMatrix, xlim=c(-3000, 3000),xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency")
plotAvgProf函数同时查看多个样本结合的强度:
##lapply函数进行批量处理
tagMatrixList <- lapply(files, getTagMatrix, windows=promoter)
##conf定义置信区间，facet决定从上到下还是从左往右
plotAvgProf(tagMatrixList, xlim =c(-3000, 3000), conf=0.95, resample=500, facet="row")
04
peak注释
这里用annotatePeak函数来对peak进行注释（准备好前面那两个文件）。值得一提的是，在注释上，ChIPseeker没有物种限制，当然物种得有那两个文件。
查看peaks在基因组上的分布：
##指定tssRegion（启动子区域），一般TSS上下游区域作为启动子区域
f = getSampleFiles()[[4]]
peakAnno = annotatePeak(f, tssRegion=c(-1000, 1000), TxDb=txdb)> peakAnno
Annotated peaks generated by ChIPseeker
1331/1331 peaks were annotated
Genomic Annotation Summary:
Feature Frequency
9 Promoter 48.1592787
4 5' UTR 0.7513148
3 3' UTR 4.2073629
1 1st Exon 0.7513148
7 Other Exon 3.9068370
2 1st Intron 6.5364388
8 Other Intron 4.8835462
6 Downstream (<=300) 1.1269722
5 Distal Intergenic 29.6769346查看peaks左右5kb区域都有哪些基因：
##addFlankGeneInfo看peak附近基因，flankDistance看peak左右5kb距离
peakAnno = annotatePeak(f[[4]], tssRegion=c(-1000, 1000), TxDb=txdb,addFlankGeneInfo=TRUE,flankDistance=5000)
> as.GRanges(x) %>% head(3)
GRanges object with 3 ranges and 14 metadata columns:
seqnames ranges strand | V4 V5
<Rle> <IRanges> <Rle> | <factor> <numeric>
[1] chr1 815093-817883 * | MACS_peak_1 295.76
[2] chr1 1243288-1244338 * | MACS_peak_2 63.19
[3] chr1 2979977-2981228 * | MACS_peak_3 100.16
annotation geneChr geneStart geneEnd geneLength
<character> <integer> <integer> <integer> <integer>
[1] Distal Intergenic 1 803451 812182 8732
[2] Promoter 1 1243994 1247057 3064
[3] Promoter 1 2976181 2980350 4170
geneStrand geneId transcriptId distanceToTSS
<integer> <character> <character> <numeric>
[1] 2 284593 uc001abt.4 -2911
[2] 1 126789 uc001aed.3 0
[3] 2 440556 uc001aka.3 0
flank_txIds
<character>
[1] uc001abt.4
[2] uc001aea.2;uc001aeb.2;uc001aec.1;uc001aed.3;uc010nyi.2;uc009vjx.3;uc009vjy.1;uc001aee.2;uc001aef.2;uc001aeg.2;uc001aeh.2;uc001aei.2;uc001aek.2;uc009vjz.2;uc010nyj.2
[3] uc001aka.3;uc010nzg.1;uc001akc.3;uc001ake.3;uc001akf.3;uc009vlh.3
flank_geneIds
<character>
[1] 284593
[2] 116983;116983;116983;126789;126789;126789;54973;54973;54973;54973;54973;54973;54973;54973;54973
[3] 440556;440556;63976;63976;63976;63976
flank_gene_distances
<character>
[1] -2911
[2] -1979;-19;0;0;0;-128;8492;15729;15729;15729;15729;15729;15729;15729;15729
[3] 0;0;-4514;-4514;-4514;-4514
-------
seqinfo: 24 sequences from hg19 genome在输出结果中可以看到多出三列，分别是：flank_txIds, flank_geneIds和flank_gene_distances，可以看到在peak附近5kb区域的所有基因。
虽然找到了peak附近基因，但是基因ID我们不认识，这里用annoDb参数进行转换，利用的是TxDb对象，TxDb对象中的基因ID是哪种，annoDb参数就会将基因ID转换为对应的ID。
peakAnno = annotatePeak(f, tssRegion=c(-1000, 1000), TxDb=txdb, annoDb = 'org.Hs.eg.db')
library(GenomicRanges)
library(dplyr)
as.GRanges(peakAnno) %>% head(3)
这里，TxDb的基因ID是Entrez，annoDb参数就会将基因ID转成ENSEMBL ID，并且加上对应的注释信息。
05
注释信息可视化
用plotAnnoPie函数可视化peak分布（ceas也可以看）：
peakAnno <- annotatePeak(files[[4]], tssRegion =c(-3000, 3000),TxDb=txdb, annoDb="org.Hs.eg.db")
plotAnnoPie(peakAnno)
用plotAnnoBar函数同时查看多个样本的peak分布：
peakAnnoList <- lapply(files, annotatePeak,
TxDb =txdb,tssRegion=c(-3000, 3000))
plotAnnoBar(peakAnnoList)
用plotDistToTSS函数可视化peak离最近基因距离的分布：
plotDistToTSS(peakAnno,
title="Distribution of transcription factor-binding loci\nrelative to TSS")
用plotDistToTSS函数同时看多个样本peak最近基因距离的分布：
plotDistToTSS(peakAnnoList)
vennplot函数可视化展示注释最近的基因在不同样本中的overlap情况：
genes <- lapply(peakAnnoList, function(i)
as.data.frame(i)$geneId)
vennplot(genes[2:4])
//
今天的分享就到这里，有问题请留言吧！