如何使用ChIPseeker进行peak注释
这篇文章主要讲解了“如何使用ChIPseeker进行peak注释”,文中的讲解内容简单清晰,易于学习与理解,下面请大家跟着小编的思路慢慢深入,一起来研究和学习“如何使用ChIPseeker进行peak注释”吧!
创新互联建站长期为上千客户提供的网站建设服务,团队从业经验10年,关注不同地域、不同群体,并针对不同对象提供差异化的产品和服务;打造开放共赢平台,与合作伙伴共同营造健康的互联网生态环境。为宜良企业提供专业的成都网站建设、网站设计,宜良网站改版等技术服务。拥有10年丰富建站经验和众多成功案例,为您定制开发。
ChIPseeker是使用的最广泛的peak注释软件之一,提供了以下多种功能
peak在染色体和TSS位点附近分布情况可视化
peak关联基因注释以及在基因组各种元件上的分布
获取GEO数据库中peak的bed文件
多个peak文件的比较和overlap分析
首先我们需要输入peak文件,支持两种格式,第一种是BED格式,最少只需要3列内容记录peak的染色体位置就可以了,示意如下
当然也可以有多余的列,只需要符合BED格式的标准即可;另外一种和MACS的peak calling输出结果类似,第一行为表头,示意如下
通过函数readPeaks
读取peak文件,用法如下
peak <- readPeakFile("peak.bed")
函数根据文件名称的后缀来判断是否为bed格式,建议BED格式的输入文件后缀统一成.bed
, 当然压缩文件也是支持的,比如.bed.gz
;如果不是BED格式的输入,文件名称则不能使用BED格式对应的后缀。
下面来详细看下几个主要功能的代码和结果展示
1. peak 在染色体上的分布
用法如下
covplot(peak, chr = c("chr1", "chr2"))
输出结果示意如下
2. peak 在TSS位点附件的分布
用法如下
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
# 定义TSS上下游的距离
promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
tagMatrix <- getTagMatrix(peak, windows=promoter)
tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red")
热图每一行代表一个基因,展示的是所有基因TSS两侧的分布,除了热图外,还可以对所有基因取均值,用折线图来展示TSS两侧分布情况,用法如下
plotAvgProf(
tagMatrix,
xlim=c(-3000, 3000),
xlab="Genomic Region (5'->3')",
ylab = "Read Count Frequency")
输出结果示意如下
3. pea关联基因注释
用法如下
peakAnno <- annotatePeak(
peak,
tssRegion = c(-3000, 3000),
TxDb = txdb,
annoDb = "org.Hs.eg.db")
write.table(
as.data.frame(peakAnno),
"peak.annotation.tsv",
sep="\t",
row.names = F,
quote = F)
注释文件内容如下
给出了关联的基因以及对应的基因组区域的类别,根据这个结果,可以提取关联基因进行下游的功能富集分析,比如提取geneid
这一列,用clusterProfiler
进行GO/KEGG等功能富集分析。
注释的结果还提供了多种可视化方式,其中饼图最为常见,用法如下
plotAnnoPie(peakAnno)
输出结果示意如下
4. 下载GEO中的peak文件
以hg19为例,首先查询对应的GEO编号信息,用法如下
> hg19 <- getGEOInfo(genome="hg19", simplify=TRUE)
> head(hg19)
series_id gsm organism
111 GSE16256 GSM521889 Homo sapiens
112 GSE16256 GSM521887 Homo sapiens
113 GSE16256 GSM521883 Homo sapiens
114 GSE16256 GSM1010966 Homo sapiens
115 GSE16256 GSM896166 Homo sapiens
116 GSE16256 GSM910577 Homo sapiens
由于列数太多,上述结果只展示了部分信息,对于每个bed文件,会列出对应的描述信息,方便筛选感兴趣的peak进行下载,可以根据GSM
编号进行下载,用法如下
downloadGSMbedFiles("GSM521889", destDir="hg19")
也可以根据下载一个基因组对应的所有peak文件,用法如下
downloadGEObedFiles(genome="hg19", destDir="hg19")
5. peak间的overlap分析
peak的overlap分析不仅可以探究生物学重复样本间的一致性,还可以进一步识别多种蛋白或者转录因子在调控网络中的作用,如果两个蛋白的chip结果overlap显著,很可能这两个蛋白构成了复合体,或者两种蛋白具有相互作用,这对于探究其调控机制有相当大的帮助。用法如下
enrichPeakOverlap(
queryPeak = peak_setA,
targetPeak = c(peak_setB, peak_setC),
TxDb = txdb,
pAdjustMethod = "BH",
nShuffle = 1000,
chainFile = NULL,
verbose = FALSE)
依次将query的peak与target中的每一个peak文件进行overlap分析,计算出一个p值代表两个peak之间overlap的程度,p值越小,overlap的程度越高。
ChIPseeker除了peak基因注释的基本功能外,整合了GEO的下载功能与peak的overlap分析,可以方便的将自己的chip_seq数据与GEO的公共数据集进行比较分析。
感谢各位的阅读,以上就是“如何使用ChIPseeker进行peak注释”的内容了,经过本文的学习后,相信大家对如何使用ChIPseeker进行peak注释这一问题有了更深刻的体会,具体使用情况还需要大家实践验证。这里是创新互联,小编将为大家推送更多相关知识点的文章,欢迎关注!
网站名称:如何使用ChIPseeker进行peak注释
URL分享:http://ybzwz.com/article/iegoco.html