ChIP-seq 数据质控与预处理
- ChIP-seq 技术概述
ChIP-seq（染色质免疫沉淀测序）是一种广泛应用于识别全基因组范围内转录因子结合位点和表观遗传标记的技术。
实验步骤包括：
- 通过交联固定蛋白质-DNA 复合物。
- 使用特异性抗体富集目标蛋白或 DNA 片段。
- 对富集的 DNA 片段进行末端修复、加 A 尾以及连接 Illumina 测序接头。
- 最后对 DNA 片段进行单端或双端测序。
- 数据格式
原始 ChIP-seq 测序数据通常以 FASTQ 格式存储。本文将分析小鼠 MEL 和 Ch12 细胞系中转录因子 Myc 的全基因组结合模式。我们从 ENCODE 网站下载了小鼠 MEL 细胞系样本 ENCSR000EUA（重复 1）的 Myc ChIP-seq 原始测序数据。
- 数据处理与质量控制
3.1 准备工作
使用 R 包 ShortRead 加载并检查序列数据质量。
library(ShortRead)
3.2 数据采样
为了节省时间和内存，我们可以从 FASTQ 文件中随机抽取一部分 reads 进行质量评估。以下代码从文件中随机抽取 1 百万条 reads：
fqSample <- FastqSampler("~/Downloads/ENCFF001NQP.fastq.gz", n = 1e6)
fastq <- yield(fqSample)
fastq 是一个 ShortReadQ 对象，包含有关测序循环数、reads 中的碱基对数量以及内存中 reads 数量的信息。
3.3 质量评估
使用访问器函数提取 FASTQ 文件中的信息：
sread()：获取 reads 序列。quality()：获取 reads 的质量分数（ASCII 格式）。id()：获取 reads 的 ID。
readSequences <- sread(fastq)
readQuality <- quality(fastq)
readIDs <- id(fastq)
3.4 质量分布分析
计算每个 read 的总质量分数，并生成直方图以观察其分布：
totalQualities <- alphabetScore(readQuality)
histogramData <- data.frame(ReadQuality = totalQualities)
library(ggplot2)
ggplot(histogramData, aes(x = ReadQuality)) +
geom_histogram(binwidth = 10) +
theme_minimal()
3.5 碱基频率统计
统计每种碱基（A、C、G、T 和 N）在整个序列中的频率：
baseFrequencies <- alphabetFrequency(sread(fastq))
summaryFrequencies <- colSums(baseFrequencies)
summaryFrequencies[c("A", "C", "G", "T", "N")]
3.6 按测序循环分析碱基分布
利用 alphabetByCycle() 函数按测序循环统计碱基出现频率，并绘制折线图：
cycleBaseFrequencies <- alphabetByCycle(sread(fastq))
plotData <- data.frame(
Count = c(cycleBaseFrequencies["A", ], cycleBaseFrequencies["C", ],
cycleBaseFrequencies["G", ], cycleBaseFrequencies["T", ]),
Cycle = rep(seq_len(ncol(cycleBaseFrequencies)), 4),
Base = rep(c("A", "C", "G", "T"), each = ncol(cycleBaseFrequencies))
)
ggplot(plotData, aes(x = Cycle, y = Count, color = Base)) +
geom_line() +
theme_bw()
3.7 平均质量分数随测序循环的变化
将 ASCII 质量分数转换为数字矩阵，并绘制箱线图以观察质量分数随测序循环的变化：
qualityMatrix <- as(quality(fastq), "matrix")
boxplot(t(qualityMatrix[1:1000, ]))
- 数据过滤
为减少内存占用，设置一个 FastqStreamer 对象以分批读取文件。以下代码逐批读取文件，过滤掉低质量 reads 和含有过多非特异性碱基 (N) 的 reads，并将过滤后的结果写入新的 FASTQ 文件：
streamer <- FastqStreamer("~/Downloads/ENCFF001NQP.fastq.gz", n = 1e5)
totalReads <- 0
filteredReads <- 0
while (length(batch <- yield(streamer)) > 0) {
totalReads <- totalReads + length(batch)
# 过滤质量分数低于阈值的 reads
highQuality <- batch[alphabetScore(batch) > 300]
# 过滤含有过多非特异性碱基的 reads
cleanReads <- highQuality[rowSums(alphabetFrequency(sread(highQuality))[, "N"] < 10)]
filteredReads <- filteredReads + length(cleanReads)
writeFastq(cleanReads, "filtered_ENCFF001NQP.fastq.gz", mode = "a")
}
