使用 R 进行 MiRNA-Seq 数据处理和生物信息学分析的经过验证的工作流程

摘要

在这里，我们提出了一种使用 R 分析 miRNA-Seq 数据的协议。该工作流程使研究人员能够探索 miRNA 调控网络及其在各种生物学和临床问题中的重要性。这项工作旨在为 miRNA 生物信息学领域的新手和经验丰富的研究人员提供实用指南。

摘要

MicroRNA （miRNA） 是影响广泛生理和病理过程的关键转录后调节因子。随着高通量测序技术的进步，miRNA-Seq 已成为分析 miRNA 表达模式的强大工具。然而，对此类数据的可靠解释需要标准化且可重复的分析管道。在这里，我们提出了一个经过验证的工作流程，用于使用 R 进行 miRNA-Seq 数据处理和生物信息学分析。该协议包含所有基本步骤，包括原始数据预处理、质量控制、比对、定量、归一化、差异表达分析、靶点预测、功能富集和调控网络构建。该工作流程旨在实现灵活性和透明度，集成了广泛采用的 R 包，并支持特定物种的注释和模块化定制。此外，还指导用户利用精选数据库和可视化工具（如 Cytoscape）进行下游生物学解释。该协议不仅支持稳健的统计分析，而且还可以对 miRNA-mRNA 相互作用及其在疾病机制中的作用有意义的见解。它特别适合进行 miRNA 生物标志物发现、疾病建模或综合多组学研究的新手和经验丰富的研究人员。

引言

MicroRNA （miRNA） 是短的非编码 RNA 分子，通过在转录后阶段¹ 起作用显着影响基因表达。它们通常通过与靶信使 RNA （mRNA） 的 3' 非翻译区 （UTR） 中的互补序列结合来发挥作用，导致 mRNA 降解或翻译抑制¹。在过去的二十年里，miRNA越来越被认为是各种生物过程的中心调节因子，包括细胞增殖、分化、细胞灭亡、免疫反应和^{器官发育2}。此外，miRNA表达失调与癌症、心血管疾病、神经系统疾病和^肾脏疾病等多种疾病的发病机制有关3。这些发现凸显了 miRNA 不仅作为治疗靶点，而且作为临床诊断中微创生物标志物的潜力。

随着下一代测序（NGS）技术的出现，miRNA的研究进入了一个新的时代。与仅限于已知 miRNA 的基于微阵列的方法不同，miRNA 测序 （miRNA-Seq） 能够对不同样品类型和条件下的已知和新型 miRNA 进行全面、高通量和无偏见的分析⁴。miRNA-Seq 提供卓越的灵敏度、准确性和动态范围，使其成为研究生理和病理环境中 miRNA 表达模式和发现调控机制的首选方法⁵。然而，miRNA-Seq数据的分析带来了特定的计算挑战，包括处理短读长，去除接头序列，区分密切相关的miRNA家族成员，以及管理读取计数的高冗余^6。这些特性需要精心设计和标准化的分析工作流程。

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

在该协议中，我们提出了一个完全在 R 中实现的经过验证且完整的 miRNA-Seq 分析工作流程，目的是为处理 miRNA 表达数据的研究人员提供可重复且用户适应性的解决方案。该工作流程从原始测序读数的质量控制和接头修剪开始，然后与参考基因组或已知 miRNA 序列进行比对。后续步骤包括读取计数的定量、归一化、差异表达分析、靶基因预测、功能富集和网络可视化。该工作流程包含几个广泛使用且维护良好的 R 包，确保可靠性以及与未来更新和扩展的兼容性。

该协议的核心优势之一在于它能够超越差异表达结果并提供有意义的生物学解释。通过整合经过验证和预测的 miRNA-mRNA 相互作用的精选数据库，该工作流程允许用户识别生物学相关的靶基因。然后可以对这些靶标进行基因本体和通路富集分析，以发现受影响的生物过程和分子通路。在最后一步中，可以使用 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。检查 RNA 完整性数 （RIN） 并确保其≥ 7.0 以实现可靠的测序。使用分光光度计或荧光计记录浓度。
3. 构建小 RNA 文库：使用商用小 RNA-seq 文库制备试剂盒从 1 μg 总 RNA 制备测序文库。按照试剂盒方案连接接头、逆转录和扩增 cDNA。通过大小选择（例如，18-30 nt插入片段）纯化PCR产物，以富集miRNA片段。
4. 序列文库：将文库加载到高通量测序平台上。设置读取长度为 ~50 bp 的单端测序的运行配置。确保每个样本生成大约 1000 万个原始读数，以达到足够的深度。
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. 查看每个碱基的质量分数、读取长度分布和接头污染：在 Web 浏览器中打开为每个 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 -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)
      注意：计数数据 （counts.csv） 和样品组 （metadata.csv） 需要提供给 DESeq2。这里的“条件”阐明了提供的样本组。对于个别要求，请参阅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 = 'Differentially Expressed miRNAs')

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")
      注意：这里我们以“has-miR-21-5p”为例，获得经过验证的靶标。
   2. 提取唯一的靶基因符号进行富集分析：
      genes <- unique(target_results@data$target_symbol)

6.进行功能富集分析

1. 执行 GO 富集
   1. 加载丰富工具：
      library(clusterProfiler)
library(org.Hs.eg.db)
      注意：在这里，我们加载一个包含人类基因组注释的数据库，这对于转换常见的基因标识符很有用。
   2. 对生物过程运行GO富集分析：
      ego <- enrichGO(gene = genes,
OrgDb = org.Hs.eg.db,
keyType = "SYMBOL",
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05)
dotplot(ego)
      注意：对于 enrichGO 的使用，我们提供了一份基因列表，并澄清它们在此处被提名为“符号”。我们进行了生物过程富集，对应于参数'ont =“BP”'。我们执行多次测试校正，因此我们指定 pAdjustMethod = “BH”。对于显著阈值，我们选择 pvalueCutoff = 0.05。dotplot 呈现可视化结果。有关更多单独选项，请参阅 clusterProfiler¹³ 的手册。
2. 进行 KEGG 通路富集
   1. 运行 KEGG 富集：
      ekegg <- enrichKEGG(gene = genes, organism = 'hsa')
dotplot(ekegg)
      注意：对于enrichKEGG的使用，我们提供了基因列表，并将生物体阐明为人类（“has”）。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 脚本。该数据集对16个不同大小的肾囊肿（最小囊肿：小于1-5 mL，n = 10;中型囊肿：10-25 mL之间，n = 4;大囊肿：大于50 mL，n = 4）和来自四个PKD1多囊肾的微囊性组织（MCT，n = 7，包括1个重复）进行了全局miRNA分析。此外，从三个诊断为孤立性肾细胞癌的肾切除术标本中获得了无恶性肿瘤的肾皮质组织，并作为正常对照 （n = 4）。为了找出最明显变化的miRNA，我们比较了来自小、中和大囊肿的样本与正常对照组织。如 图1A 所示，ADPKD样品的miRNA表达模式明显与对照样品不同。随后，我们利用生物分析仪工具（R 包 DESeq2）在两组之间产生了显着不同表达的 miRNA。另一个工具 Limma 似乎也提供了合理的结果。 图1B 显示了所有显着变...

讨论

由于读数尺寸小且冗余，miRNA-Seq 数据的分析面临明显的挑战，因此严格的质量控制和预处理至关重要。工作流程中最重要的步骤之一是适配器修剪。由于 miRNA 大约有 22 个核苷酸长，因此如果去除不当，接头序列很容易主导读数。未能执行准确的修剪可能会导致误对和误报读数膨胀。同样，在对齐之前应实施质量过滤，以消除可能影响测绘精度的低质量基础。映射和量化是另一个重要阶段。与标准 mRNA-Seq 数据不同，标准 mRNA-Seq 数据的读数通常跨越外显子，miRNA 读数紧凑，必须与已知的 miRNA 位点完美对齐。使用具有严格参数的 Bowtie 可确保精确比对，尤其是在映射到 miRBase 等数据库中的成熟 miRNA 序列时。在基因组比对与 miRBase 指数比对之间做出选择应取决于实验设计 - miRBase 比对适用于表达定量，而当需要发现新的 miRNA 时，基因组比对更好。

该工作流程适用于不同的实验设置。例如，当排序深度低于预期时，Bowtie 中的对齐参数可能...

材料

本文使用的材料清单

姓名	公司	目录编号	评论
Agilent-021827 人 miRNA 微阵列	安捷伦	/	用于分析人类样本 microRNA 的商业阵列
领结	约翰·霍普金斯大学	http://bowtie-bio.sourceforge.net/index.shtml	用于将测序读数与长参考序列对齐的软件工具
clusterProfiler（R 包）	生物导体	https://bioconductor.org/packages/clusterProfiler/	专为高通量生物数据的功能富集分析和可视化而设计的 R 包。
剪接	开源	https://cutadapt.readthedocs.io	一种命令行工具，可从高通量测序读数中去除接头序列、引物、poly-A 尾部和其他不需要的片段。
细胞景观	Cytoscape 联盟	https://cytoscape.org/	一个开源软件平台，专为复杂生物网络的可视化和分析而设计。
DESeq2（R 封装）	生物导体	https://bioconductor.org/packages/DESeq2/	专为计数数据的差异基因表达分析而设计的 R 包
增强型火山（R 包）	生物导体	https://bioconductor.org/packages/EnhancedVolcano/	一个 R 包，旨在创建出版质量的火山图。
快速质量控制	巴伯拉罕生物信息学	https://www.bioinformatics.babraham.ac.uk/projects/fastqc/	用于高通量测序数据的开源质量控制工具。
功能计数	子阅读 / SourceForge	http://subread.sourceforge.net/	用于计数映射到基因组特征的读数的程序
HTSeq 计数	Python 包	https://htseq.readthedocs.io	一种命令行工具，用于计算有多少比对的高通量测序读数重叠基因组特征，例如基因或外显子。我
Illumina Human v2 MicroRNA 表达磁珠芯片	Illumina	/	用于分析人类样本 microRNA 的商业阵列
multiMiR（R包）	生物导体	https://bioconductor.org/packages/multiMiR/	一个 R 包，提供最大的预测和实验验证的 microRNA 集成集合–靶向相互作用及其与疾病和药物的关联。
组织。Hs.eg.db（R 包）	生物导体	https://bioconductor.org/packages/org.Hs.eg.db/	专为人类（智人）基因组学研究而设计的注释包。
R 软件	R项目	https://www.r-project.org/	统计计算开源项目
Rstudio 工作室	假设 PBC	/	集成开发环境有助于提高 R 和 Python 的工作效率
SAM工具	开源	http://www.htslib.org/	用于处理下一代测序 （NGS） 数据的软件包。