关于程序员:ChIPseq-分析原始数据质控2

1. ChIPseq 简介

染色质免疫沉淀，而后进行深度测序 (ChIPseq) 是一种成熟的技术，能够在 全基因组范畴内辨认转录因子联合位点和表观遗传标记。

1.1. 试验解决

• 交联和蛋白质联合的 DNA。
• 通过抗体富集特定蛋白质或 DNA。
• 增加 末端修复、A 尾和 Illumina adapters。
• 从任一端 / 两端测序。

2. 数据格式

原始 ChIPseq 测序数据将采纳 FASTQ 格局。

在此 ChIPseq 研究中，咱们将钻研小鼠 MEL 和 Ch12 细胞系中转录因子 Myc 的全基因组联合模式。

咱们能够从 Encode 网站检索原始测序数据。在这里，咱们应用小鼠 MEL 细胞系、样品 ENCSR000EUA（反复 1）下载 Myc ChIPseq 的测序数据。

3. 数据处理

3.1. 解决筹备

一旦咱们下载了原始 FASTQ 数据，咱们就能够应用 ShortRead 包来查看咱们的序列数据品质。

首先咱们加载 ShortRead 库。

library(ShortRead)

咱们将应用 ShortRead 包中的函数查看原始测序读数。这相似于咱们为 RNAseq 执行的 QC。

不须要查看文件中的所有 reads 即可理解数据品质。咱们能够简略地查看 reads 的子样本并节俭一些工夫和内存。

请留神，当咱们进行子采样时，咱们会从整个 FASTQ 文件中检索随机 reads。这很重要，因为 FASTQ 文件通常按其在测序仪上的地位排序。

3.2. 数据读取

咱们能够应用 ShortRead 包中的函数从 FASTQ 文件中进行子采样。

在这里，应用 FastqSampler 和 yield 函数从 FASTQ 文件中随机抽取定义数量的 reads。在这里，咱们对 100 万次 reads 进行了子采样。这应该足以理解数据的品质。

fqSample <- FastqSampler("~/Downloads/ENCFF001NQP.fastq.gz", n = 10^6)
fastq <- yield(fqSample)

生成的对象是一个 ShortReadQ 对象，显示无关循环数、reads 中的碱基对和内存中的 reads 数的信息。

fastq

3.3. 数据质控

如果违心，咱们能够应用咱们相熟的拜访器函数来评估 FASTQ 文件中的信息。

• sread() – 检索 reads 序列。
• quality() – 检索 reads 品质作为 ASCII 分数。
• id() – 检索 reads 的 ID。

readSequences <- sread(fastq)
readQuality <- quality(fastq)
readIDs <- id(fastq)
readSequences

3.4. 质量检查

咱们能够为咱们的子采样 FASTQ 数据查看一些简略的质量指标。首先，咱们能够查看整体读取的品质分数。

咱们将 alphabetScore() 函数与咱们的读取品质一起应用，以检索子样本中每个读取的总和品质。

readQuality <- quality(fastq)
readQualities <- alphabetScore(readQuality)
readQualities[1:10]

而后咱们能够生成品质分数的直方图，以更好地理解分数的散布。

library(ggplot2)
toPlot <- data.frame(ReadQ = readQualities)
ggplot(toPlot, aes(x = ReadQ)) + geom_histogram() + theme_minimal()

3.5. 碱集频率

咱们能够别离应用 alphabetFrequency() 和 alphabetByCycle() 函数查看 reads 中 DNA 碱基的呈现以及整个测序周期中 DNA 碱基的呈现。在这里，咱们查看序列读取中 A、G、C、T 和 N（未知碱基）的总体频率。

readSequences <- sread(fastq)
readSequences_AlpFreq <- alphabetFrequency(readSequences)
readSequences_AlpFreq[1:3,]

一旦咱们在序列读取中取得了 DNA 碱基的频率，咱们就能够检索所有读取的总和。

summed__AlpFreq <- colSums(readSequences_AlpFreq)
summed__AlpFreq[c("A", "C", "G", "T", "N")]

3.6. 数据评估

咱们能够应用 alphabetByCycle() 函数按循环查看 DNA 碱基呈现状况。

readSequences_AlpbyCycle <- alphabetByCycle(readSequences)
readSequences_AlpbyCycle[1:4, 1:10]

咱们常常绘制此图以可视化循环中的碱基产生状况，以察看任何偏差。首先咱们将基频排列成一个数据框。

AFreq <- readSequences_AlpbyCycle["A",]
CFreq <- readSequences_AlpbyCycle["C",]
GFreq <- readSequences_AlpbyCycle["G",]
TFreq <- readSequences_AlpbyCycle["T",]
toPlot <- data.frame(Count = c(AFreq, CFreq, GFreq, TFreq), Cycle = rep(1:36, 4),
    Base = rep(c("A", "C", "G", "T"), each = 36))

当初咱们能够应用 ggplot2 绘制频率

ggplot(toPlot, aes(y = Count, x = Cycle, colour = Base)) + geom_line() + theme_bw()

咱们还能够评估周期内的均匀读取品质。这将使咱们可能确定是否存在品质随工夫降落的问题。

为此，咱们首先应用 as(read_quality,“matrix”) 函数将咱们的 ASCII 品质分数转换为数字品质分数。

qualAsMatrix <- as(readQuality, "matrix")
qualAsMatrix[1:2,]

咱们当初能够应用箱线图可视化跨周期的品质。

boxplot(qualAsMatrix[1:1000,])

在这种状况下，reads 品质分数和 read 品质随工夫的散布看起来还不错。咱们常常心愿一起拜访 FASTQ 样本，以查看是否有任何样本合乎这些指标。在这里，咱们察看到第二批低质量分数，因而将删除一些品质分数低和未知碱基高的读数。

4. 数据过滤

咱们将心愿节俭内存使用量，以容许咱们解决加载大文件。这里咱们设置了一个 FastqStreamer 对象来一次读入 100000 次读取。

fqStreamer <- FastqStreamer("~/Downloads/ENCFF001NQP.fastq.gz", n = 1e+05)

当初咱们遍历文件，过滤读取并写出过滤读取的 FASTQ。咱们正在过滤低质量的读数和具备许多非特异性 (N) 碱基调用的读数。

TotalReads <- 0
TotalReadsFilt <- 0
while (length(fq <- yield(fqStreamer)) > 0) {TotalReads <- TotalReads + length(fq)
    filt1 <- fq[alphabetScore(fq) > 300]
    filt2 <- filt1[alphabetFrequency(sread(filt1))[, "N"] < 10]
    TotalReadsFilt <- TotalReadsFilt + length(filt2)
    writeFastq(filt2, "filtered_ENCFF001NQP.fastq.gz", mode = "a")
}

本文由 mdnice 多平台公布