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

Posted 2023-04-27

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

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

.bed文件

由上游测序数据处理而来

TxDb文件

TxDb.Hsapiens.UCSC.hg19.knownGene

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

如果实在没有TxDb呢？

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

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

makeTxDbFromUCSC：通过UCSC在线制作TxDb

makeTxDbFromBiomart: 通过ensembl在线制作TxDb

makeTxDbFromGRanges：通过GRanges对象制作TxDb

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

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

1.peakAnnoList

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

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

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

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

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

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

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

Chip-seq peak annontation

Chip-seq peak annontation

Chip-seq peak annontation

PeRl

narrowPeak/boardPeak

narrowPeak/boardPeak 是ENCODE可提供下载的两种 Chip-seq 经过参考人类基因组mapping后的关于peak的数据.
其他类型的seq数据储存个数可以参看FAQformat

narrowPeak

数据按照以下规则储存:

1. string chrom: "Reference sequence chromosome or scaffold"
2. uint   chromStart: "Start position in chromosome"
3. uint   chromEnd: "End position in chromosome"
4. string name: "Name given to a region (preferably unique). Use . if no name is assigned"
5. uint   score: "Indicates how dark the peak will be displayed in the browser (0-1000) "
6. char[1]  strand: "+ or - or . for unknown"
7. float  signalValue: "Measurement of average enrichment for the region"
8. float  pValue: "Statistical significance of signal value (-log10). Set to -1 if not used."
9. float  qValue: "Statistical significance with multiple-test correction applied (FDR -log10). Set to -1 if n-ot used."
10. int   peak: "Point-source called for this peak; 0-based offset from chromStart. Set to -1 if no point-sour-ce called."

boardPeak

数据按照以下规则储存:

1. string chrom: "Reference sequence chromosome or scaffold"
2. uint   chromStart: "Start position in chromosome"
3. uint   chromEnd: "End position in chromosome"
4. string name: "Name given to a region (preferably unique). Use . if no name is assigned"
5. uint   score: "Indicates how dark the peak will be displayed in the browser (0-1000) "
6. char[1]  strand: "+ or - or . for unknown"
7. float  signalValue: "Measurement of average enrichment for the region"
8. float  pValue: "Statistical significance of signal value (-log10). Set to -1 if<BR> not used."
9. float  qValue: "Statistical significance with multiple-test correction applied (FDR -log10). Set to -1 if n-ot used."

在接下去的peak annotation中,只演示narrowPeak.bed格式数据.

示例数据

演示数据来自ENCODE 的 H3K9me3 的Chip-seq,样本ID为ENCFF199BLM.

样本的基本信息:

• Homo sapiens liver male adult (32 years)
  • Target: H3K9me3
  • Lab: Bing Ren, UCSD
  • Project: Roadmap

narrowPeak 数据基本信息:

##   seqnames     start       end       name score strand signalValue  pValue
## 1    chr10 100134639 100134831  Peak_7974   178      .     4.41261 6.68330
## 2    chr10 100446376 100446664  Peak_4189   244      .     5.13248 8.41215
## 3    chr10 100779568 100779699 Peak_32369   101      .     3.34436 4.50936
## 4    chr10  10088147  10088346 Peak_28509   112      .     4.05830 4.84890
## 5    chr10 101149252 101149594  Peak_6146   211      .     4.89252 7.53305
## 6    chr10 101173156 101173424 Peak_32985    94      .     3.45278 4.33257
##    qValue peak
## 1 2.96490  132
## 2 4.37217  165
## 3 1.47374  105
## 4 1.69566   90
## 5 3.67439  174
## 6 1.30993  237

构建注释文件

在peak annotation中需要一个用于参考的注释文件,需要包括一下信息: 1. 染色体名 2. 转录本起始位点 3. 转录本终止位点 4. gene的名字 5. 正反链

注释文件可以直接从外部导入,也可以利用biomaRt 生成

library(ChIPpeakAnno)
library(biomaRt)
mart <- useMart("ensembl")
datasets <- listDatasets(mart)
mart <- useDataset("hsapiens_gene_ensembl",mart)
# 需要筛选的特征
props <- c("ensembl_gene_id", "external_gene_name", "transcript_biotype", "chromosome_name", "start_position", "end_position", "strand")

# 筛选染色体号
lincRNA <- subset(
  getBM(attributes=props, mart=mart, filters = "chromosome_name", values = c(1:22,"X","Y")),
  transcript_biotype == "lincRNA"
  )
# 在染色体前加"chr"保持和narrowPeak数据一致
lincRNA[,4] <- paste0("chr", lincRNA[,4])

得到的数据包含以下信息:

##   ensembl_gene_id external_gene_name transcript_biotype chromosome_name
## 1 ENSG00000276255       RP5-881P19.7            lincRNA            chr1
## 2 ENSG00000234277          LINC01641            lincRNA            chr1
## 3 ENSG00000238107      RP11-495P10.5            lincRNA            chr1
## 4 ENSG00000274020          LINC01138            lincRNA            chr1
## 5 ENSG00000225620      RP11-569A11.2            lincRNA            chr1
## 6 ENSG00000237520       RP11-443B7.2            lincRNA            chr1
##   start_position end_position strand
## 1      228073909    228076550     -1
## 2      227393591    227431035      1
## 3      148295180    148297556      1
## 4      148290889    148519604     -1
## 5      202632428    202632911      1
## 6      234957231    234959989      1

进行peak annotation

有了这个注释文件以后,我们就可以根据我们想要的筛选规则对peak进行注释,主要用到的包是 ChIPpeakAnno.用到的函数为annotatePeakInBatch.

# 将前面得到的注释文件转换为RangedData对象
library(ChIPpeakAnno)
myCustomAnno <- RangedData(
  IRanges(
    start=lincRNA[,"start_position"], 
    end=lincRNA[,"end_position"],
    names=lincRNA[,"ensembl_gene_id"]),
  space=lincRNA[,"chromosome_name"],
  strand=lincRNA[,"strand"])
# 读入需要注释的narrowPeak数据
bed_file <- read.table("ENCFF199BLM.bed", header = T, sep = "\t", stringsAsFactors = F)
#　将peak数据转换为GRanges
peaks <- toGRanges(bed_file, format="narrowPeak", colNames = colnames(bed_file))
#　根据需要进行筛选
anno <- annotatePeakInBatch(peaks, AnnotationData=myCustomAnno, 
                            output="overlapping", 
                            FeatureLocForDistance="TSS",
                            bindingRegion=c(-2000, 2000))

从anno中提取我们需要的lincRNA的ID

result_lincRNA <- [email protected]$feature
head(result_lincRNA)

## [1] "ENSG00000237579" "ENSG00000235180" "ENSG00000232259" "ENSG00000204365"
## [5] "ENSG00000226578" "ENSG00000260137"