ChIP-Seq 分析流程-下游(2)

生信菜鸟团

发布于 2025-02-25 19:52:34

8700

代码可运行

文章被收录于专栏：生信菜鸟团生信菜鸟团

运行总次数：0

代码可运行

文接上回：

https://mp.weixin.qq.com/s/7gADGKEthliI-1viN1FC7w

这次我们主要聊一下 peaks 基因组注释和富集分析。

peaks 注释

Peaks 注释主要是看一下我们发现的 peaks 在哪些基因上（附近），同时看一下这些peaks 主要在哪些功能元件上。

首先加载我们需要的R包

# 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 <- list.files("./", pattern= ".bed", full.names=T)
samplefiles <- as.list(samplefiles)
names(samplefiles) <- c("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 <- list.files("./", pattern= ".bed", full.names=T)
samplefiles <- as.list(samplefiles)
names(samplefiles) <- c("Nanog", "Pou5f1")

# 步骤 2: 读取文件并自动识别转换
bed_data_list <- lapply(names(samplefiles), function(sample_name) {
  file <- samplefiles[[sample_name]]
  
  # 读取 BED 文件
  bed_data <- read.table(file, header=FALSE, stringsAsFactors=FALSE)
  
  # 自动识别并转换染色体名称
  bed_data$V1 <- gsub("^NC_00(\\d{3})\\.\\d+$", "chr\\1", bed_data$V1)  # 转换为 chr1-22
  bed_data$V1 <- gsub("^NC_0000(\\d)\\.\\d+$", "chr\\1", bed_data$V1)  # 转换为 chr1-9
  bed_data$V1 <- gsub("^NC_000023\\.\\d+$", "chrX", bed_data$V1)  # 转换为 chrX
  bed_data$V1 <- gsub("^NC_000024\\.\\d+$", "chrY", bed_data$V1)  # 转换为 chrY
  
  # 输出到新文件
  output_file <- paste0(sample_name, "_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 <- TxDb.Hsapiens.UCSC.hg19.knownGene

就可以查看注释结果了

> 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)

相对于 TSS 的 TF 结合位点分布

plotDistToTSS(peakAnnoList, title="Distribution of transcription factor-binding loci \n relative to TSS")

导出数据

# 将entrez ID 转换为 symbol ID 然后导出数据
nanog_annot <- data.frame(peakAnnoList[["Nanog"]]@anno)

# Get the entrez IDs
entrez <- nanog_annot$geneId

# Return the gene symbol for the set of Entrez IDs
annotations_edb <- AnnotationDbi::select(EnsDb.Hsapiens.v75,
                                         keys = entrez,
                                         columns = c("GENENAME"),
                                         keytype = "ENTREZID")

# Change IDs to character type to merge
annotations_edb$ENTREZID <- as.character(annotations_edb$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 文件，并没有任何条目被富集到，图片是示例图，不过一般正常的样本都可以富集得到。

GO

# Run GO enrichment analysis 
ego <- enrichGO(gene = entrez, 
                    keyType = "ENTREZID", 
                    OrgDb = org.Hs.eg.db, 
                    ont = "BP", 
                    pAdjustMethod = "BH", 
                    qvalueCutoff = 0.05, 
                    readable = TRUE)

# Dotplot visualization
dotplot(ego, showCategory=50)

KEGG

ekegg <- enrichKEGG(gene = entrez,
                 organism = 'hsa',
                 pvalueCutoff = 0.05)

dotplot(ekegg)

比较样本间的富集情况

我们的数据集由两个不同的转录因子峰值调用组成，因此能够在同一图中并排比较功能富集结果会很有用。

# Create a list with genes from each sample
genes = lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)

# Run KEGG analysis
compKEGG <- compareCluster(geneCluster = genes, 
                         fun = "enrichKEGG",
                         organism = "human",
                         pvalueCutoff  = 0.05, 
                         pAdjustMethod = "BH")
dotplot(compKEGG, showCategory = 20, title = "KEGG Pathway Enrichment Analysis")