CHIP-Seq(7):对peak进行注释和可视化

作者: Z_bioinfo | 来源:发表于2022-04-12 17:24 被阅读0次

CHIP-Seq(7):对peak进行注释和可视化
GREAT网页版学习笔记（预测顺式调控元件）
Txdb创建和读取
Chip-seq call peak和注释
使用UPORA对peak进行注释
day27 ChIP-seq 对peak进行注释R语言ChIPs
植物chip-seq数据的可视化（Y叔的神器chip-seeke
使用GREAT对peak进行功能注释
对peak进行功能注释的GREAT
m6A图文复现08-Peak结果可视化metaPlotR

上一步骤用IDR对重复样本peaks的一致性进行了评估，同时得到了merge后的一致性的.narrowPeak文件，接下来就是对peaks的注释。

这篇主要用Y叔的R包ChIPseeker对peaks的位置（如peaks位置落在启动子、UTR、内含子等）以及peaks临近基因的注释。

1.输入文件

.bed文件

由上游测序数据处理而来

TxDb文件

TxDb.Hsapiens.UCSC.hg19.knownGene

BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

同名的包中含有相应的文件，直接引用即可，同样hg38也有同名的包，Bioconductor提供了30个TxDb包，可以供我们使用。要与测序数据比对一致

如果实在没有TxDb呢？

使用GenomicFeatures包函数来制作TxDb对象（这也是为什么说ChIPseeker支持所有物种）

制作TxDb方法示例

用人的参考基因信息来做注释，从UCSC生成TxDb:

#设置工作路径
setwd("E:/ZYH/处理流程/chip-seq/lianxi2")
install.packages("RMySQL", repos = "https://mirrors.ustc.edu.cn/CRAN/")
#查看是否安装成功：
any(grepl("RMySQL",installed.packages()))
library(GenomicFeatures)
hg19.refseq.db <- makeTxDbFromUCSC(genome="hg19", table="refGene")

makeTxDbFromUCSC：通过UCSC在线制作TxDb

makeTxDbFromBiomart: 通过ensembl在线制作TxDb

makeTxDbFromGRanges：通过GRanges对象制作TxDb

makeTxDbFromGFF：通过解析GFF文件制作TxDb

2.下载安装ChIPseeker注释相关的包

source ("https://bioconductor.org/biocLite.R")
biocLite("ChIPseeker")
# 下载人的基因和lincRNA的TxDb对象
biocLite("org.Hs.eg.db")
biocLite("TxDb.Hsapiens.UCSC.hg19.knownGene")
biocLite("TxDb.Hsapiens.UCSC.hg19.lincRNAsTranscripts")
biocLite("clusterProfiler")
#载入各种包
library("ChIPseeker")
library(clusterProfiler)
library("org.Hs.eg.db")
library(TxDb.Hsapiens.UCSC.hg19.knownGene)

3.读入peaks文件

函数readPeakFile读入peak文件
summit.bed
bedfile <- readPeakFile("E:/ZYH/处理流程/chip-seq/lianxi2/idr_peak.narrowPeak")

4.注释peaks

这里需要注意的是，启动子区域是没有明确的定义的，所以你可能需要指定tssRegion，把基因起始转录位点的上下游区域来做为启动子区域。
有了这两个输入（BED文件和TxDb对象），你就可以跑注释了，然后就可以出结果了。

#注释主体,
peakAnnoList <- annotatePeak(bedfile, tssRegion=c(-5000, 5000), TxDb=txdb, annoDb="org.Hs.eg.db")
annotatePeak传入annoDb参数,即可进行基因ID转换，将Entrez ID转化为ENSEMBL，SYMBOL，GENENAME，peakAnnoList的结果如下：
#annotatePeak传入annoDb参数,可进行基因ID转换（Entrez，ENSEMBL，SYMBOL，GENENAME）
#注释peak上下游基因：annoDb='org.Hs.eg.db'参数，注释后会多出三列信息：ENSEMBL，SYMBOL，GENENAME

5.ChIPseeker包输出文件

1.peakAnnoList

>peakAnnoList
Annotated peaks generated by ChIPseeker
32/37  peaks were annotated
Genomic Annotation Summary:
            Feature Frequency
2      Other Intron     3.125
1 Distal Intergenic    96.875

2.查看具体信息，会显示出每个peak的位置，所在的区域
as.GRanges(peakAnnoList) %>% head(5)

>as.GRanges(peakAnnoList) %>% head(1)
GRanges object with 3 ranges and 29 metadata columns:
      seqnames              ranges strand |          V4        V5          V6
         <Rle>           <IRanges>  <Rle> | <character> <integer> <character>
  [1]     chr1 121483893-121485575      * |           .      1000           .
             V7        V8        V9       V10       V11       V12       V13       V14
      <integer> <numeric> <integer> <integer> <numeric> <numeric> <integer> <integer>
  [1]        -1  114.8300        -1      1257   4.12846   4.12846 121483919 121485533
            V15       V16       V17       V18       V19       V20        annotation
      <numeric> <integer> <integer> <integer> <numeric> <integer>       <character>
  [1]  114.8300      1147 121483892 121485575  233.4480      1259 Distal Intergenic
        geneChr geneStart   geneEnd geneLength geneStrand      geneId transcriptId
      <integer> <integer> <integer>  <integer>  <integer> <character>  <character>
  [1]         1 121306424 121313686       7263          1      647121   uc009wht.1
      distanceToTSS         ENSEMBL      SYMBOL               GENENAME
          <numeric>     <character> <character>            <character>
  [1]        177469 ENSG00000231752       EMBP1   embigin pseudogene 1

输出文件解读:
genomic annotation注释：annotation列，注释的是peak的位置，它落在什么地方了，可以是UTR，可以是内含子或者外显子

nearest gene annotation（最近基因）注释：多出来的其它列，注释的是peak相对于转录起始位点（TSS）的距离，不管这个peak是落在内含子或者别的什么位置上，即使它落在基因间区上，都能够注释到一个离它最近的基因（即使它可能非常远）

针对于不同的转录本，一个基因可能有多个转录起始位点，所以注释是在转录本的水平上进行的，可以看到输出有一列是transcriptId

两种注释策略面对不同的研究对象，第一种策略，peak所在的位置可能就是调控本身（比如你要做可变剪切的，最近基因的注释显然不是你关注的点）；而做基因表达调控的，当然promoter区域是重点，离结合位点最近的基因更有可能被调控。

在R里打输出的对象，它会告诉我们ChIPseq的位点落在基因组上什么样的区域，分布情况如何。

6.可视化,展示peak在基因组特征区域的分布以及转录因子在TSS区域的结合。

library(ggplot2)
#饼图
plotAnnoPie(peakAnnoList)
#柱状分布图
plotAnnoBar(peakAnnoList)

plotDistToTSS(peakAnnoList,title="Distribution of transcription factor-binding loci \n relative to TSS")
#covplot函数可以直接接受macs2-callpeak的输出bed文件进行画图。
#函数返回的图中，每个竖线表示一个peaks在染色体的位置，线的宽度表示peaks在染色体的分布范围。
covplot(bedfile)
#具体到指定染色体上的指定范围
covplot(bedfile, weightCol="V5", chrs=c("chr17", "chr18"), xlim=c(4e7, 5e7))+ facet_grid(chr ~ .id)
#covplot选项：
#chrs    指定想要展示染色体，默认为全部
#xlim    指定x轴坐标范围，即染色体的范围
#如果需要ggplot刻面展示，不同样本是用id来表示

image.png

7.功能富集分析

提取peakAnnolist中的基因，结合clusterProfiler包对peaks内的邻近基因进行富集注释。

# Run GO enrichment analysis 
ego <- enrichGO(unique(peakAnnoList@anno$SYMBOL), keyType = 'SYMBOL', OrgDb = org.Hs.eg.db, ont = "BP", pAdjustMethod = "BH", qvalueCutoff = 0.05)
# Dotplot visualization
dotplot(ego, showCategory=50)
# Multiple samples KEGG analysis
kk <- enrichKEGG(unique(peakAnnoList@anno$geneId),organism = 'hsa',  pvalueCutoff = 0.05)

dotplot(kk, showCategory = 20, title = "KEGG Pathway Enrichment Analysis")