使用 R 进行 miRNA-Seq 数据处理与生物信息学分析的已验证工作流程

摘要

本文介绍了一种使用 R 语言分析 miRNA-Seq 数据的实验方案。该工作流程使研究人员能够探索 miRNA 调控网络及其在多种生物学和临床问题中的重要意义。本研究旨在为 miRNA 生物信息学领域的新手和资深研究人员提供一份实用指南。

摘要

microRNA（miRNA）是关键的转录后调控因子，影响广泛的生理和病理过程。随着高通量测序技术的发展，miRNA测序（miRNA-Seq）已成为分析miRNA表达谱的有力工具。然而，对这类数据的可靠解读需要标准化且可重复的分析流程。本文介绍了一种基于R语言的经过验证的miRNA-Seq数据处理与生物信息学分析工作流程。该方案涵盖所有关键步骤，包括原始数据预处理、质量控制、序列比对、定量分析、标准化、差异表达分析、靶基因预测、功能富集分析以及调控网络构建。该工作流程设计灵活且透明，整合了广泛使用的R软件包，支持物种特异性注释和模块化定制。此外，本方案还指导用户利用精选数据库及Cytoscape等可视化工具进行下游生物学解读。该流程不仅支持稳健的统计分析，还能深入揭示miRNA-mRNA相互作用及其在疾病机制中的作用，特别适用于开展miRNA生物标志物发现、疾病建模或整合性多组学研究的初学者和资深研究人员。

引言

microRNA（miRNA）是一类短链非编码RNA分子，可通过作用于转录后阶段显著影响基因表达¹。它们通常通过与靶信使RNA（mRNA）的3'非翻译区（UTR）中的互补序列结合，导致mRNA降解或翻译抑制¹。在过去的二十年中，miRNA被 increasingly 认识为调控多种生物学过程的核心因子，包括细胞增殖、分化、凋亡、免疫反应和器官发育²。此外，miRNA表达的异常调控已被证实与多种疾病的发病机制相关，如癌症、心血管疾病、神经系统疾病和肾脏疾病³。这些发现凸显了miRNA不仅作为治疗靶点的潜力，同时也可作为临床诊断中微创性生物标志物的重要价值。

随着下一代测序（NGS）技术的出现，miRNA 的研究进入了新纪元。与仅限于已知 miRNA 的基于微阵列的方法不同，miRNA 测序（miRNA-Seq）能够对不同类型样本和条件下的已知及新 miRNA 进行全面、高通量且无偏倚的分析⁴。miRNA-Seq 具有更高的灵敏度、准确性和动态范围，使其成为研究生理与病理状态下 miRNA 表达模式及发现调控机制的首选方法⁵。然而，miRNA-Seq 数据的分析面临特定的计算挑战，包括短读长的处理、接头序列的去除、区分高度相似的 miRNA 家族成员，以及应对读数中高度冗余的问题⁶。这些特性要求建立精心设计且标准化的分析流程。

尽管已开发出多种用于miRNA-Seq数据分析的流程和软件工具，但其中许多依赖图形用户界面或固定的工作流程，限制了灵活性和可重复性⁷。相比之下，R编程环境为生物信息学分析提供了一个强大且可定制的平台⁸。R拥有丰富的软件包生态系统，支持统计建模、数据可视化以及与生物数据库的整合。这使得用户能够以透明且基于脚本的方式进行全面且可重复的分析。此外，R工作流程的模块化特性使研究人员能够根据特定实验需求，从原始数据预处理到功能解释的每一步进行个性化调整。

在本实验方案中，我们介绍了一种完全基于 R 语言实现的、经过验证且完整的 miRNA-Seq 分析流程，旨在为从事 miRNA 表达数据分析的研究人员提供一种可重复且可由用户灵活调整的解决方案。该流程从原始测序读段的质量控制和接头序列修剪开始，随后将读段比对至参考基因组或已知的 miRNA 序列。后续步骤包括读段计数的定量、标准化、差异表达分析、靶基因预测、功能富集分析以及网络可视化。该流程整合了多个广泛使用且持续维护的 R 软件包，确保了分析结果的可靠性，并具备与未来更新和扩展兼容的能力。

本方案的核心优势之一在于能够突破差异表达分析的结果，提供有意义的生物学解释。通过整合经过验证和预测的miRNA-mRNA相互作用 curated 数据库，该工作流程可帮助用户识别具有生物学相关性的靶基因。这些靶基因可进一步进行基因本体（Gene Ontology）和通路富集分析，以揭示受影响的生物学过程和分子通路。最后一步，可利用Cytoscape等外部工具⁹可视化miRNA-mRNA相互作用网络，从而深入了解调控格局，并识别具有潜在功能重要性的关键枢纽miRNA。

该方法已成功应用于临床研究领域，包括在肾脏疾病研究中，循环miRNA可作为诊断和预后评估的潜在生物标志物¹⁰。然而，由于该工作流程具有模块化和灵活的设计，因此适用于多种应用场景，包括疾病建模、药物反应研究、发育生物学以及比较基因组学。研究人员可轻松调整该工作流程，以适应物种特异性的注释、实验条件或额外的组学数据层次。

通过提供一种基于脚本的开源解决方案，这一以 R 语言为核心的分析流程解决了现有 miRNA-Seq 工具存在的多项主要局限性，包括定制化能力有限、依赖不透明的图形化界面、缺乏对非模式生物的支持、由于缺少版本控制而导致的可重复性差，以及难以与下游的统计和功能分析框架整合等问题。该流程使研究人员能够全面控制数据处理参数，通过版本控制的代码促进研究的可重复性，并提升生物信息学研究的透明度。随着 miRNA 在系统生物学和转化医学领域的重要性持续增加，拥有一个可靠且可灵活调整的分析框架变得愈发关键。

方案

注意：带有软件链接的材料列在材料表中。

1. 准备RNA样品并构建测序文库

注意：RNA 提取和测序应在本计算流程之外进行。miRNA 测序数据有多种分析方法。本节提供其中一种实用方法的背景说明。

1. 提取总RNA：使用专为小RNA分离优化的试剂盒（例如miRNA分离试剂盒）从生物样本中提取总RNA。严格按照制造商提供的操作说明进行操作。确保使用无RNase的耗材，并将样本置于冰上以最大限度减少降解。
2. 评估RNA完整性和浓度：取1–2 µL提取的RNA在Bioanalyzer或同等设备上进行检测。检查RNA完整性数值（RIN），确保其≥7.0，以保证测序结果可靠。使用分光光度计或荧光计记录RNA浓度。
3. 构建小RNA文库：使用商业化的small RNA-seq文库构建试剂盒，从1 µg总RNA起始制备测序文库。按照试剂盒说明书进行接头连接、逆转录和cDNA扩增。通过片段大小选择（例如18–30 nt插入片段）纯化PCR产物，以富集miRNA片段。
4. 文库测序：将文库加载至高通量测序平台。设置单端测序模式，读长约为50 bp。确保每个样本产生约1000万条原始reads，以达到足够的测序深度。
5. 导出测序数据：测序完成后，使用仪器配套的数据输出软件将原始数据导出为FASTQ文件。确认输出目录中包含序列读段及其对应的质量分数文件。将FASTQ文件存储在结构化的目录中，以便后续分析。

2. 预处理原始测序读段并进行质量控制

1. 修剪接头序列
   1. 安装并配置 Cutadapt 或 fastp。
   2. 对每个 FASTQ 文件运行接头修剪，使用以下命令：
      cutadapt -a XXXX -o trimmed_reads.fastq raw_reads.fastq
      注意：紧跟在 '-a' 后的接头序列应根据具体的 sRNA 文库制备试剂盒确定。'-o' 指定输出文件名，后接输入文件名。
2. 评估读段质量
   1. 使用 FastQC 生成质量控制报告：
      fastqc trimmed_reads.fastq
   2. 查看每个碱基的质量分数、读段长度分布以及接头污染情况：在网页浏览器中打开为每个 FASTQ 文件生成的 FastQC HTML 报告，逐步检查以下模块：
      1. 每个碱基的序列质量：确保大多数碱基位于绿色区域（Phred 分数 ≥30）。注意 3′ 端是否存在质量下降，这可能提示测序错误。
      2. 读段长度分布：确认分布符合预期的插入片段大小（例如，miRNA 为 18–30 nt）。检查是否存在异常峰。
      3. 接头含量：验证接头序列是否已有效去除。确认修剪后接头污染百分比接近零。
   3. 保存 FastQC 汇总报告，并标记任何质量指标较差的样本，以进行重新修剪或从后续分析中排除。

3. 比对序列并生成计数矩阵

1. 将测序读段比对至参考序列
   1. 下载参考基因组或成熟 miRNA 序列的 FASTA 文件（例如，从 miRBase 下载）¹¹。示例：
      wget ftp://mirbase.org/pub/mirbase/CURRENT/mature.fa
   2. 使用 Bowtie 对参考基因组建立索引。
      1. 打开终端并运行以下命令以构建索引：
         bowtie-build reference.fa reference_index
      2. 将 reference.fa 替换为实际的 FASTA 文件名。
      3. 将 reference_index 替换为索引文件所需的前缀名称。
      4. 确保 Bowtie 生成多个索引文件（例如 .ebwt）。请验证这些文件存在于工作目录中，因为比对过程需要这些文件。
   3. 使用适合短读段的参数，通过 Bowtie 进行读段比对。示例：
      bowtie -v 0 -a --best --strata reference_index trimmed_reads.fastq > aligned_reads.sam
      注意：输入文件为 trimmed_reads.fastq，输出文件为 aligned_reads.sam。'-v 0' 表示在整个读段中不允许有任何错配。'-a -best -strata' 表示丢弃比最优比对错配更多的所有比对结果。
2. 定量 miRNA 表达
   1. 使用 SAMtools 将 SAM 文件转换为 BAM 格式。
      samtools view -S -b aligned_reads.sam > aligned_reads.bam
      注意：输入文件 aligned_reads.sam 是上一步命令的结果。输出文件 aligned_reads.bam 用于后续分析。
      1. 使用 SAMtools 压缩并排序比对文件：
         samtools sort aligned_reads.bam -o aligned_reads_sorted.bam
samtools index aligned_reads_sorted.bam
      2. 确保第一条命令将 SAM 文件转换为 BAM 格式。
      3. 确保第二条命令按基因组坐标对 BAM 文件进行排序。
      4. 确保第三条命令生成索引文件（.bai），该文件是下游分析所必需的。
      5. 在进入定量分析之前，确认已成功生成排序后的 BAM 文件及其索引文件。
   2. 使用 featureCounts 或 HTSeq-count，基于 miRNA 注释 GTF 文件生成计数矩阵：
      featureCounts -a miRNA.gtf -o counts.txt aligned_reads.bam
      ​注意：featureCounts 根据 miRNA.gtf 对输入文件 aligned_reads.bam 中的读段进行定量，并输出 counts.txt 文件。

4. 在 R 中进行差异表达分析

1. 加载计数数据
   1. 将计数矩阵和样本元数据导入 R：
      library(DESeq2)
countData <- read.csv("counts.csv", row.names=1)
colData <- read.csv("metadata.csv", row.names=1)
dds <- DESeqDataSetFromMatrix(countData = countData, colData = colData, design = ~ condition)
      注：需要向 DESeq2 提供计数数据（counts.csv）和样本分组信息（metadata.csv）。此处的“condition”用于指明所提供的样本分组。具体要求请参考 DESeq2 用户手册¹²。
2. 数据标准化与转换
   1. 使用 DESeq2 的默认方法对计数数据进行标准化：
      dds <- DESeq(dds)
   2. 执行方差稳定化转换：
      vsd <- vst(dds, blind=FALSE)
   3. 使用主成分分析（PCA）可视化样本聚类：
      plotPCA(vsd, intgroup="condition")
3. 鉴定差异表达的 miRNA
   1. 提取并排序差异表达结果：
      res <- results(dds)
resOrdered <- res[order(res$pvalue), ]
      注：我们根据 pvalue 的数值对结果文件“res”进行重新排序。
      summary(res)
   2. 筛选显著差异表达的 miRNA（p 值 < 0.05，|log2FC| > 1）：
      sig_miRNA <- subset(res, pvalue < 0.05 & abs(log2FC) > 1)
      注：筛选显著变化的 miRNA 有多种阈值设定方式。“p 值 < 0.05，|log2FC| > 1”是广泛采用的标准，可根据具体数据调整阈值。
4. 可视化表达变化
   1. 安装并加载 EnhancedVolcano 软件包。
   2. 绘制火山图：
      library(EnhancedVolcano)
EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'pvalue',
title = '差异表达的 miRNA')

5. 预测miRNA的靶基因

1. 查询数据库
   1. 使用 TargetScan、miRDB 和 miRTarBase 等在线资源搜索特定 microRNA 并获取其靶基因。
   2. 重点关注经过实验验证的靶点，以提高可信度。
2. 在 R 中自动化预测
   1. 加载 multiMiR 软件包并查询已验证的靶点：
      library(multiMiR)
target_results <- get_multimir(mirna = c("hsa-miR-21-5p"), table = "validated")
      注：此处以 "hsa-miR-21-5p" 为例，获取其已验证的靶基因。
   2. 提取唯一的靶基因符号用于富集分析：
      genes <- unique(target_results@data$target_symbol)

6. 进行功能富集分析

1. 进行 GO 富集分析
   1. 加载富集分析工具：
      library(clusterProfiler)
library(org.Hs.eg.db)
      注：此处加载包含人类基因组注释的数据库，可用于转换常见的基因标识符。
   2. 对生物过程（Biological Processes）进行 GO 富集分析：
      ego <- enrichGO(gene = genes,
OrgDb = org.Hs.eg.db,
keyType = "SYMBOL",
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05)
dotplot(ego)
      注：在使用 enrichGO 时，我们提供一个基因列表，并明确此处基因标识类型为 'SYMBOL'。我们进行的是生物过程富集分析，对应参数 'ont ="BP"'。为进行多重检验校正，设定 pAdjustMethod = "BH"。显著性阈值选择 pvalueCutoff = 0.05。dotplot 可展示可视化结果。更多个性化选项，请参考 clusterProfiler¹³ 用户手册。
2. 进行 KEGG 通路富集分析
   1. 运行 KEGG 富集分析：
      ekegg <- enrichKEGG(gene = genes, organism = 'hsa')
dotplot(ekegg)
      注：在使用 enrichKEGG 时，我们提供一个基因列表，并指定物种为人类（"hsa"）。dotplot 可展示可视化结果。更多个性化选项，请参考 clusterProfiler¹³ 用户手册。

7. 构建并可视化miRNA-mRNA相互作用网络

1. 导出数据以进行网络可视化
2. 基于 TargetScan、miRDB 或 miRTarBase 生成的靶基因，创建 miRNA-靶基因配对的数据框
3. 将网络表格写入 CSV 文件：
   write.csv(miRNA_target_pairs, "network.csv")
4. 导入 Cytoscape
   1. 打开 Cytoscape 并导入网络表格
   2. 使用力导向布局或环形布局对网络进行可视化
   3. 分析拓扑特性（例如度中心性），以识别枢纽 miRNA

结果

我们从GSE133530下载了microRNA表达矩阵，并直接进行了差异表达分析。我们在补充文件1中提供了该数据集的示例分析R脚本。该数据集对来自四个PKD1多囊肾的16个不同大小的肾囊肿（微小囊肿：小于1-5 mL，n = 10；中等囊肿：10-25 mL之间，n = 4；大囊肿：大于50 mL，n = 4）以及微囊性组织（MCT，n = 7，包含1个重复）进行了全基因组miRNA表达谱分析。此外，从三例经诊断为孤立性肾细胞癌的肾切除标本中获取了无恶性病变的肾皮质组织，作为正常对照（n = 4）。为了鉴定变化最显著的miRNA，我们将小、中、大囊肿样本与正常对照组织进行了比较。如图1A所示，ADPKD样本的miRNA表达模式明显不同于对照样本。随后，我们利用生物分析工具R软件包DESeq2，鉴定出两组之间显著差异表达的miRNA。另一个工具Limma似乎也能提供合理的结果。图1B展示了所有显著变化miRNA的分布情况。...

讨论

由于miRNA-Seq数据读段较短且存在冗余性，其分析面临独特挑战，因此严格的质控和预处理至关重要。工作流程中最重要的步骤之一是接头序列的切除。由于miRNA长度约为22个核苷酸，若未正确去除接头序列，其可能会在读段中占据主导地位。若切除不准确，可能导致比对错误并增加假阳性读段的数量。同样，在比对前应实施质量过滤，以剔除可能影响比对准确性的低质量碱基。比对与定量是另一个关键阶段。与通常跨越外显子的标准mRNA-Seq数据不同，miRNA读段较短，必须与已知miRNA位点完全匹配。使用Bowtie并设置严格参数可确保精确比对，尤其是在将读段比对至miRBase等数据库中的成熟miRNA序列时。选择将读段比对至基因组还是miRBase索引，应根据实验设计而定——若目标为表达定量，使用miRBase比对更为合适；若需发现新miRNA，则基因组比对更优。

该实验流程可适用于不同的实验条件。例如，当测序深度低于预期时，可放宽Bowtie中的比对参数以允许更多的错配，但这样会增加假阳性结果的风险。如果在去除接头序列后仍存在污染，用户可使用更严格的参数重新运行C...

材料

本文使用的材料清单

姓名	公司	目录编号	评论
Agilent-021827 人 miRNA 微阵列	Agilent	/	一种用于分析人源样本 microRNA 的商业化芯片
Bowtie	约翰斯·霍普金斯大学	http://bowtie-bio.sourceforge.net/index.shtml	一种用于将测序读段比对至长参考序列的软件工具
clusterProfiler（R 软件包）	Bioconductor	https://bioconductor.org/packages/clusterProfiler/	一种用于高通量生物数据分析的功能富集分析与可视化的 R 软件包
Cutadapt	开源软件	https://cutadapt.readthedocs.io	一种命令行工具，用于从高通量测序读段中去除接头序列、引物、poly-A 尾巴及其他不需要的片段
Cytoscape	Cytoscape 联盟	https://cytoscape.org/	一种用于可视化和分析复杂生物网络的开源软件平台
DESeq2（R 软件包）	Bioconductor	https://bioconductor.org/packages/DESeq2/	一种用于计数数据差异基因表达分析的 R 软件包
EnhancedVolcano（R 软件包）	Bioconductor	https://bioconductor.org/packages/EnhancedVolcano/	一种用于生成出版级火山图的 R 软件包
FastQC	Babraham 生物信息学中心	https://www.bioinformatics.babraham.ac.uk/projects/fastqc/	一种用于高通量测序数据的质量控制开源工具
featureCounts	Subread / SourceForge	http://subread.sourceforge.net/	一种用于统计比对到基因组特征的读段数量的程序
HTSeq-count	Python 软件包	https://htseq.readthedocs.io	一种命令行工具，用于统计有多少比对上的高通量测序读段与基因或外显子等基因组特征重叠。I
Illumina Human v2 MicroRNA 表达微珠芯片	Illumina	/	一种用于分析人源样本 microRNA 的商业化芯片
multiMiR（R 软件包）	Bioconductor	https://bioconductor.org/packages/multiMiR/	一种 R 软件包，提供了最大规模的整合数据库，包含预测的和经实验验证的 microRNA–靶标相互作用信息，以及其与疾病和药物的关联数据
org.Hs.eg.db（R 软件包）	Bioconductor	https://bioconductor.org/packages/org.Hs.eg.db/	一种专为人类（智人）基因组学研究设计的注释软件包
R 软件	R 项目	https://www.r-project.org/	一个用于统计计算的开源项目
Rstudio	Posit PBC	/	一种集成开发环境，有助于提高使用 R 和 Python 的工作效率
SAMtools	开源软件	http://www.htslib.org/	一种用于处理新一代测序（NGS）数据的软件包