Rup（RNA-seq 可用性评估流程）— 真核生物批量 RNA-seq 实验的质量控制

摘要

本方案可为生物湿实验研究人员（生物信息学经验有限）提供RNA测序实验的初步质量控制。

摘要

分子植物科学的现代研究方法通常需要进行大规模 RNA 测序（bulk RNA-seq）实验，例如追踪处理条件下转录组的全局变化，或鉴定调控通路中的关键组分。因此，植物科学的多个领域都依赖高质量且可重复的大规模 RNA-seq 数据来推动科学进步。然而，根据我们的经验，RNA-seq 数据集中质量控制措施的相关知识和实际应用往往不足。本文介绍了 Rup（RNA-seq 可用性评估流程），用于大规模 RNA-seq 数据的质量控制，以支持后续的基因表达分析。Rup 是一个独立运行的工具，适用于具备基础 R 语言知识的实验生物学研究人员。Rup 可帮助区分高质量、适合下游基因表达分析的测序数据与不适合进一步通用分析的数据。Rup 包含多项针对常见问题的检测功能，例如测序读段数量或比对率不足、污染序列的识别、总 RNA-seq 数据中 rRNA 组分的定量、技术重复间的相似性检验，并通过真实数据演示提供直观的可视化结果。Rup 提供了一整套工具，可在标准化转录组分析之前识别实验中的缺陷，从而提升个体研究人员及整个领域的数据质量。该工具增强了对大规模 RNA-seq 数据分析结果的信心，并为未来制定定义最低质量控制标准的指南提供了基础，从而提高已发表 RNA-seq 数据的可靠性和透明度。Rup 及测试数据可在 https://github.com/oliverrupp/rup 获取。

引言

转录组学（RNA-seq）实验可全面分析塑造表型的转录特征。该方法在植物分子遗传学中已成为不可或缺的工具，用于鉴定单个或共表达的基因以及参与特定生物学过程的通路，例如发育、植物与病原体互作或非生物胁迫抗性^1, 2 ,3。RNA-seq 技术的最新进展提高了检测特异性，并能够在单碱基分辨率下识别不同的转录本异构体和变异形式，从而可鉴定从较大的插入缺失到单核苷酸多态性（SNPs）的序列变异。RNA-seq 所获得的数据具有广泛的动态范围，有助于同时检测高丰度和低表达水平的转录本，但需要适当的实验条件以确保数据一致性。此外，RNA-seq 数据集规模庞大，分析和存储过程对计算资源要求较高，因此需要高效的数据管理和充足的计算能力。RNA-seq 在整个实验流程的每个步骤均需要高质量的起始材料，因为任一环节的不足都可能逐级放大并影响最终结果。RNA 完整性差或文库构建过程中的技术问题将导致偏差并降低结果准确性。获取高质量原始测序读段的参数可能因不同的高通量测序（NGS）平台而异。

根据我们的经验，我们建议仅使用RNA完整性数值（RIN）高于7的RNA作为RNA-seq文库构建的起始材料，因为该数值表明mRNA结构基本完整。为确保测序成功，标准文库构建方案通常需要约2 µg总RNA，浓度为50–200 ng/µL。应通过分光光度计检测确认RNA纯度，OD_260/280比值应在1.8至2.1之间，OD_260/230比值应大于1.5。测序深度不足、比对率偏低或与参考基因组错配可能会进一步扭曲基因表达和剪接的定量结果。此外，实验设计不当，例如未能区分不同组织或处理组的样本，或重复样本间变异过大，均可能引入噪声并降低可重复性。尽管许多实验室常规开展转录组分析，但对最初关键分析步骤的严格质量控制常常未在发表文献中报告，甚至完全缺失。这可能导致对转录组分析结果的过度解读，进而造成结果不可重复。

本文提供了一套用于高质量转录组分析中mRNA表达谱初始步骤的质量控制工作流程，旨在检测转录水平的变化。我们的目标是使湿实验生物学家即使在生物信息学知识有限的情况下，也能够评估其原始转录组数据。Rup（图1）适用于具备R语言基础知识的研究人员。执行本文所述的工作流程，将使研究人员能够深入理解其原始数据及其在后续分析中可能存在的局限性。据我们所知，目前尚缺乏将实用的原始转录组数据评估流程与区分高质量与低质量数据的指导原则相结合的系统方法。

图1：RNA测序 in silico 质量控制流程的工作流程图。 质量控制的输入文件来源于对植物材料进行测序所获得的RNA测序数据，以及公开可用的数据集。所提供的R语言流程通过三个主要方面评估序列质量：测序质量、比对质量和重复样本质量。计算多种质量指标，并利用常见的R语言软件包（如ggplot2和pheatmap）对统计结果进行可视化展示（例如条形图和热图）。请点击此处查看该图的放大版本。

既往报道的评估流程需要预处理数据（例如，以 .bam 文件形式的测序读段比对结果），未能涵盖所有评估指标，或已不再维护^4,5,6。本文所展示工作流程的优势在于其全面性，能够将最常见的质量控制问题整合到单一分析流程中。此外，我们提供了适用于所有下游分析应用的高质量数据示例，同时也提供了低质量数据的示例，并讨论了这些数据在进一步分析中的具体局限性。已有若干先前发表的研究阐述了RNA-seq分析工具的目的，并对其性能进行了比较评估^7,8。然而，RNA-seq的质量控制和方法学报告标准尚未统一，这加剧了转录组学实验在可重复性以及生物学意义解读方面的困难。

Rup 整合了标准的高质量工具、质量控制分析以及结果可视化功能。Rup 需要以原始测序读段和注释好的基因组作为输入，可在 Mac OS、Linux 以及 Windows 系统上的“Windows 子系统 for Linux”（WSL）上运行。该工具集成了用于评估测序质量的程序，可量化测序读段在基因组上的比对情况，包含对样本中 rRNA 含量的测定，并提供样本相关性分析以评估重复样本之间的一致性。测试数据是通过激光显微切割技术获取植物物种Eschscholzia californica两个不同发育阶段的中央分生组织组织，采用超低起始量文库构建方案，并在 Novaseq 6000 平台上测序获得的。

方案

Rup 可以作为一个单一脚本运行，所需输入文件最少。如本示例所示，该分析流程需要 Illumina 双端 RNA 测序数据。针对单端测序数据但包含相同分析步骤的调整版流程也已存放在 GitHub 仓库中。仅需提供基因组序列的 fasta 文件、基因模型和 rRNA 注释的 gtf 文件，以及原始测序读段的 fastq.gz 文件。请确保将基因组和注释文件分别命名为 genome.fa、annotation.gtf 和 rRNA.gtf，并放置在 reference_folder 变量指定的文件夹中。当所有测序文件（.fq.gz 格式）均保存在 read_file_folder 变量指定的目录下时，可按以下方式运行 Rup：

1. 准备工作：

1. 使用 Bioconductor 安装所有必要的 R 包。

   BiocManager::install(c("getopt", "ggplot2", "reshape2", "pheatmap", "fastqcr", "Rfastp", "Rsubread", "Rsamtools"))

2. 设置所有必要的文件。
   1. 创建一个源文件夹，用于存放所有必要文件。将参考基因组序列以名为“reference/genome.fa”的 fasta 文件形式添加到源文件夹中。提供基因模型注释的 gtf 文件：“reference/annotation.gtf”。可选地，提供 rRNA 基因注释的 gtf 文件：“reference/rRNA.gtf”。
   2. 将所有测序读段作为 fastq 文件添加到“reads”文件夹中。确保所有 fastq 文件名遵循相同的命名模式：_1.fastq.gz 和 _2.fastq.gz，分别对应正向和反向读段。

      # 定义源文件夹以及所有用于输入和输出文件的子文件夹
source_folder <- "data"
reference_folder <- file.path(source_folder, "reference")
read_file_folder <- file.path(source_folder, "reads")
results_folder <- file.path(source_folder, "results")

# 定义参考输入文件
genome_fasta_file <- file.path(reference_folder, "genome.fa")
annotation_file <- file.path(reference_folder, "annotation.gtf")
rrna_file <- file.path(reference_folder, "rRNA.gtf")

# 定义结果子文件夹
fastqc_folder <- file.path(results_folder, "fastqc")
trimmed_fastqc_folder <- file.path(results_folder, "trimmed_fastqc")
trimmed_read_folder <- file.path(results_folder, "trimmed")
bam_folder <- file.path(results_folder, "bam")
sorted_bam_folder <- file.path(results_folder, "sorted_bam")
counts_folder <- file.path(results_folder, "counts")

# 创建结果文件夹
for(folder in c(results_folder, fastqc_folder, trimmed_fastqc_folder, trimmed_read_folder, bam_folder, sorted_bam_folder, counts_folder)) {
 if(!dir.exists(folder)) {
  dir.create(folder)
 }
}

# 从输入 fastq 文件获取样本前缀
fastq_files <- list.files(read_file_folder, pattern = "*_1.f(ast)?q.gz")
sample_prefixes <- gsub("_1.f(ast)?q.gz", "", fastq_files)

3. 根据计算资源设置参数。

   n_threads <- 8 # 可用 CPU 核心数
bamSortMemory <- "1024" # 用于排序 bam 文件的最大内存

4. 根据测序方法设置参数。
    

   MinReadLength <- 25 # 最小读段长度，不应小于 25
minFragLength <- 0 # 最小片段长度分布
maxFragLength <- 300 # 最大片段长度分布
orientation <- "fr" # 用于读段比对的读段方向（"fr"、"rf"、"ff"）
stranded <- 0 # 链特异性测序
         # (0（非链特异性）、1（链特异性）和 2（反向链特异性）)

   
   
   注意：最小读段长度不应小于 25 bp，因为更短的读段可能导致比对失败。较大的值会增加唯一比对的读段数量，但也会在修剪过程中导致更多读段被丢弃，具体取决于测序质量。有关这些参数的更多细节见补充材料。

2. 测序质量评估

注意：此步骤将生成柱状图，显示每个样本在修剪前后的读段数量。

1. 运行 fastqc⁹，该工具用于高通量测序数据的常规质量控制，可提供有关测序读段总数、读段长度、平均 GC 含量、测序接头污染以及过表达序列（例如 rRNA 读段）的信息。根据每个碱基和每条序列的质量报告分析测序质量。使用 R 软件包 fastqcr¹⁰ 对原始测序文件运行 fastqc，并将所有输入文件置于同一目录（read_file_folder）中。使用函数 qc_aggregate 将 fastqc_folder 中的输出文件汇总为一个 R 对象。
   注意：结果将保存在 fastqc_folder 目录中。输出文件夹中由 fastqc 生成的 .html 文件可用于查看额外的统计信息。
    

   # 加载 fastqc 库
library(fastqcr)
# 对所有原始 fastq 文件运行 fastqc
fastqc(fq.dir=read_file_folder, qc.dir=fastqc_folder, threads=n_threads)

# 汇总 fastqc 统计信息
qc <- qc_aggregate(fastqc_folder, progress=F)
qc$tot.seq <- as.numeric(qc$tot.seq)
# 去除样本名称的后缀
qc$sample = gsub(".f(ast)?q.gz", "" , qc$sample)

2. 将所有样本前缀存储在向量 sample_prefixes 中，遍历这些前缀并对所有样本运行质量剪切工具 fastp。剪切操作用于去除接头序列、引物、低质量碱基、poly-A/T 序列，由 Fastp¹¹ 执行，该工具可通过 R 软件包 Rfastp¹² 获得。剪切后的 fastq 文件将存储在 trimmed_read_folder 文件夹中。
    

   # 加载 rfastp 库
library(Rfastp)

# 遍历样本前缀
for(prefix in sample_prefixes) {
 # 创建输出前缀
 # !!! Rfastp 会自动在前缀后添加 _R1.fastq.gz 和 _R2.fastq.gz
outputPrefix <- file.path(trimmed_read_folder, prefix)

# 根据样本前缀获取输入读段
read1 = file.path(read_file_folder, paste(prefix, "_1.fastq.gz", sep=""))
read2 = file.path(read_file_folder, paste(prefix, "_2.fastq.gz", sep=""))
if(!file.exists(read1)) { read1 = file.path(read_file_folder, paste(prefix, "_1.fq.gz", sep="")) }
if(!file.exists(read2)) { read2 = file.path(read_file_folder, paste(prefix, "_2.fq.gz", sep="")) }

# 使用最小读段长度运行 fastp 剪切
fastp_stats <- rfastp(read1 = read1,
     read2 = read2,
     minReadLength = minReadLength,
     outputFastq = outputPrefix,
     thread = n_threads)
}

3. 对剪切后的读段再次运行 fastqc。
    

   # 对剪切后的读段运行 fastqc
fastqc(fq.dir=trimmed_read_folder, qc.dir=trimmed_fastqc_folder, threads=n_threads)

# 汇总 fastqc 统计信息
qc_trimmed <- qc_aggregate(trimmed_fastqc_folder, progress=F)
qc_trimmed$tot.seq <- as.numeric(qc_trimmed$tot.seq)
# 修正样本名称（将 fastp 的 "R1"、"R2" 后缀改为 "1"、"2"）
qc_trimmed$sample = gsub("_R([12])$", "_\\1" , qc_trimmed$sample)

4. 通过收集所有样本在剪切前后的读段数量来评估剪切效果。使用 ggplot2¹³ 创建读段数量的柱状图。
    

   # 加载 ggplot2 库用于绘图
library(ggplot2)

# 将剪切状态添加到 fastqc 结果中
qc$Trimming = "raw"
qc_trimmed$Trimming = "trimmed"

# 将结果合并为一个向量
qc_all = rbind(qc, qc_trimmed)

# 绘制剪切前后读段数量的柱状图
read_number_plot <- ggplot(qc_all, aes(x=sample, y=tot.seq, fill=Trimming)) +
geom_bar(stat="identity", position = "dodge") +
xlab("样本") + ylab("读段数量") +
theme(axis.text.x = element_text(angle = 90, hjust = 0)) +
theme(text = element_text(size = 18)) +
ggtitle("剪切前后读段数量")
print(read_number_plot)

图 2：测序与修剪结果。修剪前（红色）和修剪后（绿色）的读段数。样本 s2_r1 在修剪前读段数已较低，而修剪去除了样本 s2_r2 大部分的读段。其余所有样本在修剪过程中读段损失均处于可接受范围。请点击此处查看此图的放大版本。

3. 比对质量

注意：此步骤通过计算柱状图，以区分唯一比对的读段（例如，来自蛋白质编码基因的转录本）、多重比对的读段（例如，rRNA 读段）以及未比对的读段（例如，污染序列）。

1. 使用 Rsubread¹⁴ 使用软件包将测序读段比对至参考基因组。首先，构建参考基因组fasta文件（genome_fasta_file）的索引。
   注意：此步骤针对每个参考基因组仅执行一次。
    

   # 加载 Rsubread 库以进行读段比对和读段计数
library(Rsubread)

# 从参考基因组序列创建 subread 索引
buildindex(file.path(参考文件夹, "subread.index"), 基因组_fasta_文件)

2. 遍历所有样本，并使用比对工具将测序读段比对至参考基因组() Rsubread 包的功能。该函数会在输出文件夹中生成 .bam 文件。
    

   # 加载 Rsamtools 库以进行 bam 文件的排序和索引
library(Rsamtools)

# 遍历样本名称
用于（前缀 在 样本前缀) {
# 创建 bam 输出文件名
output_bam = file.path(bam_folder, paste(prefix, ".bam", sep=""))

# 将测序读段比对到参考基因组索引
映射统计信息 <- 将索引对齐（index=file.path(参考文件夹,"subread.index")，
     # 根据前缀创建正向和反向读段文件名
     # “_R1.fastq.gz”和“_R2.fastq.gz”由 fastp 强制规定
     readfile1=file.path(修剪后读段文件夹, paste(前缀, "_R1.fastq.gz", sep=""))
     readfile2 = 文件路径(修剪后读段文件夹, paste(前缀, "_R2.fastq.gz", sep = "")),
     输出文件=output_bam, # 设置输出文件名
     类型=0, # 读段为RNA测序
     minFragLength=minFragLength, # 允许的最短片段长度
     maxFragLength=maxFragLength, # 允许的最小片段长度
     PE_方向=方向, # 读取方向
     nthreads=n_threads,
     # 使用GTF注释文件以支持比对映射
     使用注释=TRUE,
     annot.ext=注释文件,
     isGTF=真实,
     nBestLocations=2 # 区分唯一比对读段与多重比对读段

# 创建排序后输出文件的名称（.bam 后缀将自动添加）
output_sorted_bam = file.path(sorted_bam_folder, prefix)

# 对BAM文件进行排序和索引
sortBam(output_bam, output_sorted_bam, maxMemory=bamSortMemory, nThreads=n_threads)
indexBam(paste0(output_sorted_bam, ".bam"))
}

3. 使用 Rsubread 软件包中的 featureCounts() 函数对每个基因的测序读段进行计数。确保注释文件（annotation_file）为 GTF 格式。仅对与基因组唯一比对的读段进行计数。
    

   # 收集所有 bam 文件
bam文件 <- list.files(bam_folder, pattern = "*.bam$")

# 统计每个基因的读段数量
     gene_feature_counts = featureCounts(file.path(bam_folder, bam_files),
     annot.ext = 注释文件， # 基因模型
     isGTF注释文件 = TRUE,
     countMultiMappingReads = 错误, # 仅计数唯一序列读段
     链特异性 = 链特异性， # 用于链特异性数据
     isPairedEnd = TRUE, # 双端测序
# 以下参数可用于控制假阳性读段分配
     requireBothEndsMapped = 真实,
     checkFragLength = 真实,
     minFragLength = minFragLength,
     maxFragLength = maxFragLength,
     nthreads = n_threads

4. 统计比对到rRNA基因的reads数量，以评估样本中的rRNA含量。允许在此步骤中计入多重比对的reads。在GTF文件（rrna_file）中指定rRNA基因位点。
    

   # 计算比对到 rRNA 基因的读段数量
rrna_feature_counts = featureCounts(file.path(bam_folder, bam_files),
     annot.ext = rrna_file, # rRNA基因模型文件
     isGTF注释文件 = TRUE,
     countMultiMappingReads = TRUE, # rRNA 读段通常会比对到多个位点
     分数 = 真实,
     strandSpecific = 链特异性，
     isPairedEnd = 真实,
     nthreads = n_threads

5. 收集由 featureCounts 生成的读段分配统计信息。这些统计信息包括每个样本在不同类别中的比对数（例如，已分配、未比对、多重比对）。
    

   # 加载 reshape2 库以进行数据重构
library(reshape2)

# 获取任务统计信息
基因计数统计 <- gene_feature_counts$stat
# 将样本名称设为列名
gene_feature_counts$counts 的列名 = gsub(".bam", "", gene_feature_counts$counts 的列名)

# 将分配类型设置为行名称
gene_count_stats 的行名 <- gene_count_stats$Status
# 选择包含分配统计信息的列
基因计数统计 <- gene_count_stats[,seq(2, ncol(gene_count_stats))]
# 将样本名称设为列名
gene_count_stats 的列名 <- gsub(".bam", "", colnames(gene_count_stats))

# 为绘图转换数据框
转化后的统计量 <- melt(t(基因计数统计))

# 设置列名
transformed_stats 的列名 <样本，组别，比对结果

# 调整 groupa 和样本顺序以用于绘图
transformed_stats$Group <- factor(转换后统计$组别,
          levels = rev(levels(transformed_stats$Group)[order(levels(transformed_stats$Group))])
transformed_stats$样本 <- factor(转换后统计$样本,
          levels = rev(levels(transformed_stats$Sample)[order(as.character(transformed_stats$Sample))]))

# 移除没有任何比对结果的分类单元
转化后的统计量 <- transformed_stats[transformed_stats$Alignments > 0,]
# 统计数据涉及编码蛋白质的基因
transformed_stats$Reference = "Genes"

6. 收集rRNA比对的统计信息。
    

   # 收集 rRNA 分配统计信息
rrna_count_stats <- rrna_feature_counts$stat
# 将样本名称设为列名
rrna_count_stats 的列名 <- gsub(".bam", "", colnames(rrna_count_stats))

# 此情况下仅分配的读段数量重要
rrna_counts = as.data.frame(t(rrna_count_stats[rrna_count_stats$Status == "Assigned",2:ncol(rrna_count_stats)]))
colnames(rrna_counts) = "比对数"

# 设置绘图的名称和类别
rrna_counts$样本 = 行名(rrna_counts)
rrna_counts$Group = "rRNA"
rrna_counts$参考 = "rRNA"

7. 将读段比对统计结果绘制成柱状图。
    

   # 合并蛋白编码基因和rRNA基因的统计结果
mapping_stats = rbind(转换后的统计数据, rrna_counts)

# 绘制分配统计图
mapping_stats_plot = ggplot(data = mapping_stats, aes(x = 样本, y = 比对数)) +
 geom_col(aes(fill = 组别), width = 0.7) +
 theme_bw() + facet_wrap(~Reference) +
 ylab("比对数") + xlab("样本") +
 ggtitle("读段比对数量") +
 theme(text = element_text(size = 18)) +
 theme(axis.text.x = element_text(angle = 90, hjust = 0))
print(mapping_stats_plot)

8. 根据分配给基因的读段数量将基因分类，并将结果以柱状图形式绘制。
    

   # 获取每个基因的读段计数
gene_count_matrix = gene_feature_counts$counts
# 将样本名称设为列名
gene_count_matrix 的列名 = gsub(".bam", "", gene_count_matrix 的列名)

# 按特定读段数类别对基因进行分组与计数
read_count_classes = data.frame("无reads"=colSums(gene_count_matrix == 0),
     "至少1个读段"=colSums(基因计数矩阵 >= 1 & 基因计数矩阵 < 10),
     "至少10个读段"=colSums(基因计数矩阵 >= 10 & 基因计数矩阵 < 100),
     "至少100条读段"=colSums(基因计数矩阵 >= 100 & 基因计数矩阵 < 1000),
     "至少1000条读段"=colSums(基因计数矩阵 >= 1000))

# 为绘图设置数据框
read_count_classes$样本 = 行名(read_count_classes)
melt_rcc = melt(读取计数分类)
melt_rcc$variable = as.character(melt_rcc$variable)

# 对基因组进行排序
melt_rcc$variable = factor(melt_rcc$variable, levels=c("no_reads",
         至少1条读段
         至少10次读取
         至少100次读取
         至少1000条读段

# 将数据绘制成柱状图
基因覆盖度图 <- ggplot(melt_rcc, aes(x=样本, y=值, fill=变量)) +
geom_bar(stat="identity") +
ylab("基因数量") + xlab("样本") +
ggtitle("每个基因的读段数") +
guides(fill=guide_legend(title="已分配读段的数量")) +
主题（文本 = 元素_文本（大小 = 18））+
theme(axis.text.x = element_text(angle = 90, hjust = 0))
print(基因覆盖度图)

图3：读段比对概览与统计结果。 大多数样本显示出较高数量的已分配读段（左侧，棕色）和较低数量的rRNA读段（右侧，蓝色）。样本s2_r3具有异常高数量的多位置比对读段（左侧，绿色），对应于较高的rRNA读段数（右侧）。样本s2_r4显示出大量未比对读段，同时rRNA读段数处于预期水平，提示可能存在来自其他生物体的读段污染。 请点击此处查看该图的放大版本。

4. 重复实验的质量。

1. 确保来自同一样本或组织的重复样本之间的相关性高于来自不同样本或组织的重复样本之间的相关性。使用原始读段数或经每百万转录本（Transcripts Per Million, TPM）¹⁵标准化后的读段数，绘制重复样本相关性的聚类热图以可视化数据。
    

   # 加载 pheatmap 程序包
library(pheatmap)

# 计算 TPM 值
Length_kb = gene_feature_counts$annotation$Length / 1000 # 获取以 kb 为单位的基因长度
RPK = gene_feature_counts$counts / Length_kb # 根据基因长度对读段数进行标准化
scaling_factors = colSums(RPK) / 1e6 # 基于标准化后读段数的总和计算缩放因子
TPM = RPK / scaling_factors # 将标准化后的读段数缩放至 1e6

# 绘制样本/重复样本相关性热图
# 皮尔逊相关系数基于 log2 转换后的 TPM 值计算
pheatmap(cor(log2(TPM+1)), fontsize = 18, main="样本相关性热图")

图4：基因读数分类。 样本中的所有基因根据分配到的读数数量被分为五类之一。红色显示的是未分配到任何读数的基因数量。总分配读数较少的样本（s2_r1 至 s1_r4）具有较多分配到10至100个读数的基因，而分配到超过1000个读数的基因数量较少。这些样本中可能遗漏低表达水平的基因。 请点击此处查看该图的放大版本。

结果

该流程完全以 R 脚本实现，并已在 Linux 和 Mac OS 操作系统上进行测试。Windows 用户可使用 Windows Subsystem for Linux (WSL)。代码和测试数据作为 GitHub 仓库提供：https://github.com/oliverrupp/rup。测序数据可在 ENA EBI 项目 PRJEB96400 下获取。

通过将两个真实样本人工构建为十个样本，以展示在使用 Rup 进行分析的 bulk RNA-seq 质量控制过程中可能遇到的各种问题。样本 s2_r1 被设计为总读段数较低，样本 s2_r2 包含大量低质量读段，将在剪切过程中被剔除。样本 s2_r3 包含大量 rRNA 读段，而样本 s2_r4 包含未能比对到参考基因组的污染读段。样本 s1_r5 和 s2_r5 的名称被互换，以说明重复样本间相关性低的情况。

第2个实验方案部分能够在修剪前或修剪后识别出读数较低的样本（图2）。柱状图显示，样本s1_r1在修剪前后读数均较低，此处初始读数原本就偏低。样本s1_r2在修剪后读数偏低，提示其中含有大量接头序列、测序错误、引物序列或被测序的poly-A/T片段，在修剪过程中被去除。输入RNA的降解也可能导致高质量读段数量减少。

第3部分的实验方案指出了读段分配中存在的问题。图3显示了样本s2_r3（绿色）中多位置比对读段数量升高，以及核糖体RNA（rRNA）读段数量较多的情况。样本s2_r4的污染情况表现为大量读段未能比对到参考基因组（粉色）。这些读段并不对应于rRNA，而是来源于非目标生物体的序列。图4展示了转录组中分配读段数量过低所引发的普遍问题。在比对到参考基因组的读段比例较低的样本（s2_r1至r4）中，仅有约三分之一的基因被分配到超过100个读段；而在其他样本中，约有一半的基因属于此类。在s2_r1至r4样本中，低表达水平基因的读段可能无法被检测到，因此使用这些样本进行的比较表达分析将极不可靠，应予以避免。

第4部分方案可用于识别重复样本中的离群值。相同样本/条件/组织的重复样本之间预期具有比其他样本/条件/组织的重复样本之间更高的相关性。图5展示了两个样本（S1和S2）各自五个重复样本的相关性热图。图上方和侧边的树状图显示了两个聚类簇，每个簇包含五个重复样本。左侧簇包含样本1的四个重复样本和样本2的一个重复样本（s2_r5），右侧簇包含样本2的四个重复样本和样本1的一个重复样本（s1_r5）。在此情况下，若将样本名称s1_r5和s2_r5重新对调，则每个簇均包含单一样本的全部重复样本，这可能提示存在重复样本标记错误。重复样本未能聚在一起的其他原因可能包括样本/条件/组织之间缺乏差异性，或重复样本仅因测序质量因素而聚类。后一种情况可能发生在所有在剪切和比对后具有异常低或异常高读段数的重复样本自行形成一个簇时。

图 5：样本相关性热图。该样本相关性热图基于 log₂ 转换后的 TPM 值构建，显示了两个明显聚类，每个聚类包含五个样本。生物学/技术学重复样本之间的相关性应高于来自其他组织/处理组的样本。左侧聚类包含样本 1 的四个重复和样本 2 的一个重复（s2_r5），右侧聚类包含样本 2 的四个重复和样本 1 的一个重复（s1_r5）。应对这些单个重复样本进行检查，以确认是否存在可能的样本混淆或批次效应，从而解释其在热图中的聚类情况。请点击此处查看该图的放大版本。

表1：RNA-seq质量评估流程的比较。 详细分析请参见补充文件1。请点击此处下载该表格。

补充文件 1：参数和参考选择。 这些分析展示了参数和参考选择对整体质控结果的影响。请点击此处下载该文件。

讨论

差异基因表达分析的质量高度依赖于两个因素：每个样本及重复中测序的读段数量¹⁶，以及每个样本的重复次数^17,18。本文介绍了用户友好的Rup分析流程，用于确定每个重复中适用于基因表达定量的读段数量。通过多种评估指标，研究人员可以了解为何某些重复的比对读段数偏低，并识别样本间相关性方面存在的问题。尽管Rup最初是为植物RNA-seq样本的质量评估而开发，但它同样适用于其他真核生物，且在用于非植物样本时无需额外调整（补充文件1）。

Rup 的第一步是确定读段修剪前后总的读段数量。测序读段的数量决定了可检测基因的数量，如果读段数量过低，许多差异表达的基因将无法被检测到。在质量修剪和过滤过程中被丢弃的读段占比较大时，可能表明 RNA 存在降解，这可能导致对高度降解基因的表达水平被低估。然而，一个样本是否可用于后续分析，取决于目标生物体以及研究目的。例如，根据我们的经验，要检测大多数植物基因的表达，需要 3000 万至 5000 万条读段，而对于真菌，可能仅需 1000 万条读段就已足够。因此，Rup 不会设定排除问题样本的质量阈值，而是提供一系列指标，用于识别可能遇到的各种问题。

第二步（比对质量）评估测序读段与参考基因组比对的准确性和可靠性，以确定其是否适用于后续分析。并非所有测序得到的读段都可用于基因丰度计算；在大多数情况下，无法比对到基因组或转录组的读段，或比对到基因组多个位置的读段将被忽略。¹⁹ 且不计入总体读段计数。分析流程第二步的结果可用于理解为何某些读段未被用于丰度计算。未比对读段数量较高可能表明RNA提取过程中存在污染（例如植物病原体或植食性生物的污染）或参考基因组不完整。参考基因组中基因模型不完整可能导致大量“无特征”读段。多重比对读段占比较高可能源于文库中rRNA基因数量过多，提示在测序文库构建过程中rRNA去除不充分（若这些多重比对读段确实来自rRNA，可通过向分析流程提供rRNA注释文件以进一步分析，流程将自动计算可能的rRNA读段数量）。

第三个同样重要的质量指标是同一样本重复之间的相关性。通常情况下，同一样本重复之间的相关性应高于不同样本重复之间的相关性。重复样本之间相关性较低可能表明重复间存在较大的生物学变异、样本/条件/组织之间高度相似、其他批次效应，甚至样本混淆或标签错误。Rup 可计算所有样本间的成对相关性，并生成样本相关性的聚类热图。此外，还可对样本进行主成分分析（PCA），以识别在后续下游分析中需要关注的潜在批次效应。尽管批次效应可在后续分析中进行校正，但如果检测到的差异过大，更合适的做法可能是剔除存在问题的样本。

起始材料直接影响测序质量：组织取样、RNA提取和文库构建是减少RNA-seq质量问题的关键步骤。应尽可能保持一致的实验条件，例如使用生长培养箱。应及时采取病虫害防控措施，仅选择健康个体进行取样。建议在同一天和/或同一时间点采集样品，以降低昼夜节律引起的转录组变异。市面上有多种RNA提取试剂盒，针对目标物种选择合适的试剂盒可提高mRNA质量。文库构建方案中应包含polyA+ RNA富集步骤，以最大限度减少rRNA比例。

Rup 可作为差异基因表达分析中初步的质量控制步骤。该质量控制因多个关键原因而必需，因为它可确保输入数据的质量（识别低质量的重复样本/样品，设定测序深度和比对率的质量阈值），并发现诸如测序错误、批次效应或样本标签错误等技术问题。如果输入的 RNA 质量足够高，Rup 可通过修剪低质量的测序读段或识别标签错误的样本提供帮助。然而，Rup 无法弥补输入 RNA 质量差的问题。在 RNA 质量较低或测序数据存在错误的情况下，可能需要重新采集样本、调整 RNA 提取方案和/或重复测序。对于因污染导致大量未比对读段的样本，同样适用此处理原则。然而，根据我们的经验，并非所有 RNA 提取都能因材料不足而简单重复。在此情况下，这些样本不适用于某些下游应用：含有较少唯一比对读段的样本不应用于差异基因表达分析，但仍可用于转录本存在性分析。此时，在组织/处理/条件下转录本的缺失情况不能纳入转录本丰度的分析与定量中。

Rup 存在一些局限性。例如，它不能直接检测 RNA 降解，因为在文库构建和测序之前需要对 RNA 完整性进行直接测定。然而，本流程的代码仓库中包含一个使用 RSeQC 模块 geneBodyCoverage.py 和 tin.py 来识别 RNA 降解的脚本。此外，Rup 不检测 GC 含量和转录本长度偏差。目前参考基因组的大小限制为 4 Gb，根据我们的经验，读段比对模块在多倍体中会高估多重比对的读段数量，这一点在补充文件中通过 E. californica 单倍体与二倍体基因组版本的读段比对比较得到了证实，因此该流程更推荐使用单倍体基因组作为输入。 Rup 的输出质量在很大程度上依赖于参考基因组和注释的质量，不完整或高度碎片化的基因组序列可能导致未比对读段数量被高估。此外，基因模型注释不完整会导致未分配读段数量被高估。基因组序列和基因注释的完整性与重复性可通过 BUSCO²⁰ 等工具进行评估。

Rup 是一个独立工具，适用于 RNA-seq 实验中所有必需的初始质量控制步骤。与 RSeQC 和 RNA-SeQC 不同，Rup 无需对原始测序读段进行预处理。RNA-SeQC 无法进行测序质量分析，而 RSeQC 和 RNA-SeQC 均不计算样本相关性热图（表 1）。此外，RSeQC 的标准化输出为 FPKM 格式，且 RSeQC 和 RNA-SeQC 均未将输出可视化作为所有模块的标准功能。因此，这两种替代工具未能检测若干质量控制问题，包括读段数量过低、修剪后读段比例过高以及重复样本离群值的识别。Rup、RSeQC 和 RNA-SeQC 的比较详见 补充文件 1。此外，Rup 可与 RNA-SeQC 或 RSeQC 互补使用；例如，本流程生成的已排序 BAM 文件可作为这些工具的输入文件。

目前，Rup 已针对少量样本进行了优化，但最耗时的步骤（如质量修剪和读段比对）可以预先计算，例如在计算机集群或云基础设施上完成。对于大于 4 Gb 的基因组，这一步是必需的，因为 Rsubread 比对工具仅支持小于 4 Gb 的基因组。rRNA 注释由 barrnap 完成，该工具可识别高度保守的 rRNA 基因。根据我们的经验，rRNA 基因注释无需完全覆盖，因为即使某些异常的 rRNA 基因缺失，仅对数据集中 rRNA 存在情况的总体概览也足以用于数据集质量评估。未来版本的 Rup 可能会切换至其他比对方法，例如伪比对工具 salmon²¹，以减少运行时间。此外，还可向 Rup 中添加更多标准化和校正方法，如 GC 偏好性或长度偏好性校正。

总之，Rup 提供了确保 RNA-seq 数据分析可靠性和可重复性的关键信息。该流程全面报告了主要的 RNA-seq 质控指标，并生成可直接用于下游分析的输出文件。它被设计为一个独立工具，便于生物信息学知识有限的研究人员通过直观的可视化方式评估其原始测序数据的质量。

材料

本文使用的材料清单

姓名	公司	目录编号	评论
fastqcr	R	0.1.3
FUJITSU ESPRIMO D958 台式机	FUJITSU		该分析流程在 Ubuntu Linux 24.02 系统（32 Gb 内存）上开发和测试
getopt	R	1.20.4
ggplot2	R	3.5.2
MacBook Pro	Apple		该分析流程在 macOS 15.4.1 系统（16 Gb 内存）上测试
pheatmap	R	1.0.13
R		4.4.3
reshape2	R	1.4.4
rfastp	bioconductor	1.16.0
rsamtools	bioconductor	2.22.0
rsubread	bioconductor	2.20.0