写在前面的话：
参考使用的文件资料是由哈佛陈生物信息学核心 (HBC) 教学团队成员开发的。另外也看了多个公众号文章和书籍。
Website: https://hbctraining.github.io/Intro-to-ChIPseq/
Github: https://github.com/hbctraining/Intro-to-ChIPseq
文接上回：
https://mp.weixin.qq.com/s/7gADGKEthliI-1viN1FC7w
这次我们主要聊一下 peaks 基因组注释和富集分析。
peaks 注释
Peaks 注释主要是看一下我们发现的 peaks 在哪些基因上（附近），同时看一下这些peaks 主要在哪些功能元件上。
# Load libraries
library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(EnsDb.Hsapiens.v75)
library(clusterProfiler)
library(AnnotationDbi)
library(org.Hs.eg.db)
library(tidyverse)
-
接着我们需要准备我们想注释的 bed 格式文件，这个bed 文件内容取决于你想研究的问题，如果你想注释一下组间差异 peaks，那你就准备这个数据的 bed 文件。如果你想每个样本都看一下，那就是我们最开始使用 macs2 查找到的peaks bed 文件。
这里我们就看一个两个转录因子 "Nanog", "Pou5f1" 差异peaks吧
# Load data
samplefiles "./", pattern= ".bed", full.names=T)
samplefiles names(samplefiles) "Nanog", "Pou5f1")
但这里我出现了一个意外，就是我使用的是 ENSEMBL 数据库下载的参考基因组和注释文件，正常来说，与参考基因组比对后的bed文件前三列应该是
染色体
起始位点
终止位点
, 那正常而言，染色体标识应该是
chr1
Chr1
或者是
1
。但是这里情况实在是不同。！！！！
NC_000001.11
这是什么染色体号表示方式，生物数据库表示方式的不同真的是令人伤心的一件事。
NC_000001.11 11874 12227 DDX11L1
NC_000001.11 12613 12721 DDX11L1
NC_000001.11 13221 14409 DDX11L1
NC_000001.11 29321 29370 WASH7P
NC_000001.11 24738 24891 WASH7P
NC_000001.11 18268 18366 WASH7P
NC_000001.11 17915 18061 WASH7P
NC_000001.11 17606 17742 WASH7P
NC_000001.11 17233 17368 WASH7P
所以下边继续走我就会报错，我不得不转换一下ID，这里祝你的数据分析不会出现这个问题。
# 加载必要的库
library(GenomicRanges)
# 步骤 1: 读取 BED 文件
samplefiles "./", pattern= ".bed", full.names=T)
samplefiles names(samplefiles) "Nanog", "Pou5f1")
# 步骤 2: 读取文件并自动识别转换
bed_data_list function(sample_name) {
file
# 读取 BED 文件
bed_data
# 自动识别并转换染色体名称
bed_data$V1 "^NC_00(\\d{3})\\.\\d+$", "chr\\1", bed_data$V1) # 转换为 chr1-22
bed_data$V1 "^NC_0000(\\d)\\.\\d+$", "chr\\1", bed_data$V1) # 转换为 chr1-9
bed_data$V1 "^NC_000023\\.\\d+$", "chrX", bed_data$V1) # 转换为 chrX
bed_data$V1 "^NC_000024\\.\\d+$", "chrY", bed_data$V1) # 转换为 chrY
# 输出到新文件
output_file "_converted.bed")
write.table(bed_data, file=output_file, quote=FALSE, sep="\t", row.names=FALSE, col.names=FALSE)
# 返回修改后的数据框
return(bed_data)
})
# 查看转换后的数据框
print(bed_data_list)
然后我就重新读进去，接着后续的分析
txdb
就可以查看注释结果了
> peakAnnoList
$Nanog
Annotated peaks generated by ChIPseeker
70/113 peaks were annotated
Genomic Annotation Summary:
Feature Frequency
2 3' UTR 2.857143
4 Other Exon 2.857143
1 1st Intron 10.000000
5 Other Intron 35.714286
3 Distal Intergenic 48.571429
$Pou5f1
Annotated peaks generated by ChIPseeker
1137/1903 peaks were annotated
Genomic Annotation Summary:
Feature Frequency
8 Promoter 1.93491645
4 5' UTR 0.35180299
3 3' UTR 1.67106420
1 1st Exon 0.08795075
6 Other Exon 2.46262093
2 1st Intron 13.80826737
7 Other Intron 28.67194371
5 Distal Intergenic 51.01143360
plotAnnoBar(peakAnnoList)
plotDistToTSS(peakAnnoList, title="Distribution of transcription factor-binding loci \n relative to TSS")
# 将entrez ID 转换为 symbol ID 然后导出数据
nanog_annot "Nanog"]]@anno)
# Get the entrez IDs
entrez $geneId
# Return the gene symbol for the set of Entrez IDs
annotations_edb keys = entrez,
columns = c("GENENAME"),
keytype = "ENTREZID")
# Change IDs to character type to merge
annotations_edb$ENTREZID $ENTREZID)
# Write to file
nanog_annot %>%
left_join(annotations_edb, by=c("geneId"="ENTREZID")) %>%
write.table(file="results/Nanog_peak_annotation.txt", sep="\t", quote=F, row.names=F)
富集分析
一旦我们获得了峰值调用的基因注释，我们就可以进行功能富集分析，以识别这些基因中主要的生物学主题，方法是结合来自生物本体知识，例如基因本体（Gene Ontology）、KEGG 和 Reactome。
和普通 RNA-Seq 富集分析一样，我们也是获得了一个基因 list, 挖掘它的机制，哈哈哈。这里我不展示自己富集结果了，因为我选择的是一个很小的bed 文件，并没有任何条目被富集到，图片是示例图，不过一般正常的样本都可以富集得到。
