ChIP-seq分析：功能富集与基序识别技术
1. 必要软件包加载
我们使用rGREAT包中提供的Bioconductor接口进行后续分析。
# 安装并加载rGREAT包
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("rGREAT")
library(rGREAT)
2. GO功能富集分析
通过submitGreatJob函数可将我们的峰数据提交至GREAT服务器进行功能富集分析。此函数返回一个GreatJob对象，包含服务器结果的引用。
# 提交分析任务
analysisJob <- submitGreatJob(peakRegions, genome = "mm10",
version = "3.0.0", checkInterval = 1)
# 查看可用的分析类别
availableCategories(analysisJob)
使用getEnrichmentTables函数可获取特定类别的富集结果表。下面我们检索"Regulatory Motifs"类别中的结果。
# 获取功能富集结果表
enrichmentResults <- getEnrichmentTables(analysisJob, category = "Regulatory Motifs")
names(enrichmentResults)
现在我们可以查看"MSigDB Predicted Promoter Motifs"基因集中的富集情况。
# 提取启动子基序富集结果
promoterMotifs <- enrichmentResults[["MSigDB Predicted Promoter Motifs"]]
head(promoterMotifs, 4)
3. 基序分析流程
3.1 基序识别概述
转录因子ChIP-seq分析中，识别峰内富集的基序是关键步骤。虽然可以在R/Bioconductor中进行从头基序发现，但计算成本较高。本节介绍使用MEME-ChIP套件识别新基序的方法。
3.2 序列提取
首先，我们需要加载相应基因组的BSgenome对象。这里使用小鼠mm10基因组。
# 加载BSgenome包
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("BSgenome.Mmusculus.UCSC.mm10")
library(BSgenome.Mmusculus.UCSC.mm10)
我们有一个以峰最高信号点为中心的GRanges对象。
# 查看峰中心点信息
peakCenters
使用getSeq函数提取峰序列，返回DNAStringSet对象。
# 提取峰序列
regionSequences <- getSeq(BSgenome.Mmusculus.UCSC.mm10, peakCenters)
names(regionSequences) <- paste0(seqnames(peakCenters), ":", start(peakCenters),
"-", end(peakCenters))
# 查看前两条序列
regionSequences[1:2]
3.3 FASTA文件生成
使用writeXStringSet函数将序列写入FASTA格式文件，满足MEME-ChIP的输入要求。
# 写入FASTA文件
writeXStringSet(regionSequences, file = "transcriptionFactorPeaks.fa")
3.4 MEME-ChIP分析
现在"transcriptionFactorPeaks.fa"文件包含适合基序分析的序列。在实际工作中，通常会在本地安装MEME套件运行分析，但本节介绍通过门户网站提交的方法。
3.5 结果解析
从FIMO输出中可获取MEME-ChIP识别的基序位置信息。FIMO将基序位置报告为GFF3文件，可在IGV中可视化。
3.6 GFF文件解析
使用rtracklayer包可将GFF文件解析为R对象。
# 加载rtracklayer包
library(rtracklayer)
motifAnnotations <- import("~/Downloads/fimo.gff")
3.7 GFF3文件优化
为便于IGV可视化，我们给序列添加更合理的名称并导出更新后的GFF文件。
# 修改命名规则
motifAnnotations$Name <- paste0(seqnames(motifAnnotations), ":", start(motifAnnotations),
"-", end(motifAnnotations))
motifAnnotations$ID <- paste0(seqnames(motifAnnotations), ":", start(motifAnnotations),
"-", end(motifAnnotations))
# 导出优化后的GFF文件
export.gff3(motifAnnotations, con = "~/Downloads/fimoOptimized.gff")
3.8 已知基序扫描
使用Biostrings包的matchPattern函数可扫描序列中的已知基序。通常，我们可以使用数据库（如JASPAR）中的已知基序集。
# 加载JASPAR数据库
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("JASPAR2020")
library(JASPAR2020)
3.9 从JASPAR获取基序矩阵
使用TFBSTools包可获取特定转录因子的位置权重矩阵(PWM)。
# 加载TFBSTools包
library(TFBSTools)
# 获取Myc转录因子的PWM
motifMatrix <- getMatrixByName(JASPAR2020, name = "MYC")
motifMatrix
3.10 基序匹配分析
使用motifmatchr包可扫描序列中的基序并返回位置信息。
# 安装并加载motifmatchr
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("motifmatchr")
library(motifmatchr)
# 在峰序列中扫描Myc基序
matchedMotifs <- matchMotifs(motifMatrix, peakCenters,
BSgenome.Mmusculus.UCSC.mm10, out = "positions")
matchedMotifs
3.11 匹配基序导出
可将峰内的基序位置导出为BED格式，用于后续IGV可视化或元图分析。
# 导出匹配的基序位置
export.bed(matchedMotifs[[1]], con = "identifiedMotifs.bed")
