经典教程：全转录数据分析实战

本文介绍全转录组数据分析方法，我们将以拟南芥测序数据为例，在 UseGalaxy.cn 云平台进行数据分析实践。

概览

• 问题：
  • 哪些 miRNA 在对油菜素内酯的反应中上调？
  • 哪些基因是油菜素内酯诱导 miRNA 的潜在靶标？
• 目标：
  • 进行 miRNA 差异表达分析
  • 理解基于 quasi-mapping 的 Salmon 方法，用于使用 RNA-Seq 定量转录本的表达
  • 鉴定参与油菜素内酯介导调节网络的潜在 miRNA

简介

作为固着生物，植物在恶劣环境条件下的生存很大程度上取决于它们感知应激刺激并做出适当反应以抵消潜在的有害影响的能力。植物激素和活性氧化物的协调被认为是增强植物抗逆性的关键因素，允许对环境变化做出基因表达的精细调节(Planas-Riverola _et al._ 2019[1], Ivashuta _et al._ 2011[2])。这些分子构成了复杂的信号网络，赋予了植物对变化多端的自然环境做出反应的能力

油菜素内酯（BRs）是一类植物类固醇激素，对植物的生长发育以及控制生物和非生物胁迫具有重要作用。从结构上看，油菜素内酯是多羟基类固醇衍生物，与动物激素非常相似（图 1）。这类植物激素包括大约 60 种不同的化合物，其中油菜素酮（BL）、24-表油菜素酮（EBR）和 28-同油菜素酮（HBR）被认为是最具生物活性的。

图1：类固醇激素。

图 1：各种植物和动物类固醇激素的结构。

近期的研究表明，BR 介导的基因调控网络有可能改变农业的未来。其不仅可以缓解各种非生物胁迫条件（如干旱）的带来的对抗效果，还可以促进植物的生长并提高产量。例如，在番茄（_Solanum lycopersicum_）中，EBR 处理增强了作物对干旱的耐受性，提高了光合能力、叶片水分状态和抗氧化防御(Wang _et al._ 2018[3])。在辣椒（_Capsicum annuum_）中，BL 处理提高了作物在干旱条件下光利用效率(Hu _et al._ 2013[4])。暴露于干旱胁迫并接受 BL 处理的鹰嘴豆（_Cicer arietinum_）植株的重量显著增加(Anwar _et al._ 2018[5])。然而，BR 在增强植物对非生物胁迫的耐受性方面的作用机制仍然有待进一步研究。

微小 RNA（miRNA）是一类主要由 20-22 核苷酸组成的小 RNA（sRNA），其特征是可以调控基因在转录后水平上的表达。miRNA 与其他 sRNA 的不同之处在于，其产生自具有不完全螺旋结构的前体。与动物不同，植物 miRNA 的前处理发生在细胞核中（图 2）。miRNA 前提经过甲基化后被导出到细胞质，并与 Argonaute 1 蛋白结合形成 RNA 诱导沉默复合物 （RNA-induced silencing complex, RISC）。miRNA 本身不具有降解 mRNA 或干扰其翻译的能力，但它在定位靶 mRNA 时起到重要作用。

fig2:miRNA biosynthesis.

Figure 2: Plant miRNA biosynthesis, homeostasis and mechanisms of action (Wang _et al._ 2018[6]).

miRNA 已被发现是许多生理过程的重要调节因子，例如应激和激素反应。有条证据可以证明 miRNA 被认为是植物对周围环境响应的主要调节因子：

• 在特定环境条件下，多个 miRNA 基因受到调控
• 计算预测估计每个 miRNA 调节数百个基因
• 大多数植物 miRNA 调节编码转录因子（TFs）的基因
• 靶标不仅包括 mRNAs，还包括长非编码 RNAs（lncRNAs）

在植物中，miRNA 可以通过 RNA 降解和翻译抑制途径沉默靶标，与动物不同的是，大部分 miRNA 及其靶标碱基错配小于 4nt。这一特征已被利用于开发 miRNA 靶标预测工具，提供了一种有效的方法来阐明 miRNA 介导的调节网络，这些分析方法可以用于研究改善作物产量的方案。

在这个教程中，我们基于Park _et al._ 2020[7]的研究内容，探索油菜素内酯与 miRNA 基因沉默途径之间的相互作用，其被认为是植物在应激情况下最 广泛使用的调节机制之一。

课程安排 在本教程中，我们将涵盖以下内容：

1. 实验设计[8]
2. 数据背景[9]
3. 获取数据[10]
4. miRNA 数据分析[11]
   1. miRNA reads 的质量评估[12]
   2. miRNA 定量：MiRDeep2[13]
   3. miRNA 差异表达分析：DESeq2[14]
   4. 过滤显著差异表达的 miRNA[15]
5. mRNA 数据分析[16]
   1. mRNA reads 的质量评估[17]
   2. 基因表达定量：Salmon[18]
   3. mRNA 差异表达分析：DESeq2[19]
   4. 过滤显著差异表达的 mRNA[20]
6. miRNA 靶标的鉴定[21]
   1. 使用TargetFinder进行 miRNA 靶标预测[22]
7. 可选练习[23]
8. 结论[24]

实验设计

本次培训的主要目标是识别由油菜素内酯诱导表达的 miRNA 的潜在靶标。我们的初始假设是，在存在这些激素的情况下，一定会有油菜素内酯诱导的 miRNA 与其表达受到抑制的 mRNA 具有高度序列互补性（图 3）。

图3：实验设计。

图 3：实验设计

数据背景

本次培训使用的数据集可以分为三组：miRNA reads、mRNA reads 和额外数据集。

miRNA reads

miRNA 数据集包括六个 FASTQ 文件，通过使用 Illumina GAxII 测序平台获得。植物样本来自于野生型 Ws-2 苗，经过模拟处理或在收获前 90 分钟用 1 μM EBR 处理。原始数据集可在 NCBI SRA 数据库中找到，accession number 为SRP258575[25]。与之前一样，本教程将使用数据的简化版本。

mRNA reads

mRNA 数据集包括四个 FASTQ 文件，通过使用 Illumina HiSeq 2000 测序平台生成。样本来自于野生型哥伦比亚（Col-0）苗，经过模拟处理或在收获前 4 小时用 100 nM BL 处理。原始数据集可在 NCBI SRA 数据库中找到，accession number 为SRP032274[26]。为了减少分析运行时间，本教程将使用原始数据的子集。

补充数据集

除了从 NCBI 数据库获取的 RNA-Seq reads 外，我们还将使用两个来源的数据集：

• AtRTD2[27] 一个高质量的转录本参考数据集，旨在利用诸如Salmon和Kallisto等转录本定量工具的准确性来分析拟南芥的 RNA-Seq 数据。
• PmiREN[28] 一个全面的功能性植物 miRNA 数据库，包括来自多个植物物种的超过 20,000 个注释的 miRNA。

获取数据

我们分析的第一步是从 Zenodo 获取 miRNA-Seq 数据集，并将数据集组织成集合。

实践操作：检索 miRNA-Seq 和 mRNA-Seq 数据集

1. 为本教程创建一个新的历史记录

create_history

1. 从 Zenodo 导入文件：
   • 打开upload菜单
   • 点击Rule-based选项卡
   • “Upload data as”: Collection(s)
   • 复制以下表格数据，粘贴到文本框中，并按“构建数据集”按钮

SRR11611349	Control miRNA	https://zenodo.org/record/4710649/files/SRR11611349_MIRNASEQ_CTL.fastqsanger.gz	fastqsanger.gz
SRR11611350	Control miRNA	https://zenodo.org/record/4710649/files/SRR11611350_MIRNASEQ_CTL.fastqsanger.gz	fastqsanger.gz
SRR11611351	Control miRNA	https://zenodo.org/record/4710649/files/SRR11611351.MIRNASEQ_CTLfastqsanger.gz	fastqsanger.gz
SRR11611352	BR treated miRNA	https://zenodo.org/record/4710649/files/SRR11611352_MIRNASEQ_BL.fastqsanger.gz	fastqsanger.gz
SRR11611353	BR treated miRNA	https://zenodo.org/record/4710649/files/SRR11611353_MIRNASEQ_BL.fastqsanger.gz	fastqsanger.gz
SRR11611354	BR treated miRNA	https://zenodo.org/record/4710649/files/SRR11611354_MIRNASEQ_BL.fastqsanger.gz	fastqsanger.gz
SRR1019436	Control mRNA	https://zenodo.org/record/4710649/files/SRR1019436_RNASEQ_CTL.fastqsanger.gz	fastqsanger.gz
SRR1019437	Control mRNA	https://zenodo.org/record/4710649/files/SRR1019437_RNASEQ_CTL.fastqsanger.gz	fastqsanger.gz
SRR1019438	BR treated mRNA	https://zenodo.org/record/4710649/files/SRR1019438_RNASEQ_BL.fastqsanger.gz	fastqsanger.gz
SRR1019439	BR treated mRNA	https://zenodo.org/record/4710649/files/SRR1019439_RNASEQ_BL.fastqsanger.gz	fastqsanger.gz

create_history

• 从Rules 菜单选择 Add / Modify Column Definitions
• 点击 Add Definition 按钮并选择 List Identifier(s): column A

start_definition如果您选择的是上传为dataset而不是collection。关闭上传菜单，重新开始该过程，并确保您选择的是Upload data as：Collection(s)

• 点击 Add Definition 按钮并选择 Collection Name: column B
• 点击 Add Definition 按钮并选择 URL: column C
• 点击 Add Definition 按钮并选择 Type: column D
• 点击 Apply 并点击 upload 按钮

apply_definition

wait_for_data

1. 点击数据集以展开它
2. 点击Add Tags选项卡
3. 添加一个以#开头的标签
   • 以# 开头的标签将被自动识别，并对以此数据集作为输入的分析结果数据自动添加该标签。
4. 按下回车
5. 检查出现在数据集名称下的标签

接下来我们将获取剩余的数据集。

实践操作：检索额外的数据集

1. 从 Zenodo 导入文件：
   • 打开面板上的 upload菜单
   • 上传数据为：Datasets
   • 再次，复制表格数据，粘贴到文本框中，然后按“build”

SRR11611349	Control miRNA	https://zenodo.org/record/4710649/files/SRR11611349_MIRNASEQ_CTL.fastqsanger.gz	fastqsanger.gz
SRR11611350	Control miRNA	https://zenodo.org/record/4710649/files/SRR11611350_MIRNASEQ_CTL.fastqsanger.gz	fastqsanger.gz
SRR11611351	Control miRNA	https://zenodo.org/record/4710649/files/SRR11611351.MIRNASEQ_CTLfastqsanger.gz	fastqsanger.gz
SRR11611352	BR treated miRNA	https://zenodo.org/record/4710649/files/SRR11611352_MIRNASEQ_BL.fastqsanger.gz	fastqsanger.gz
SRR11611353	BR treated miRNA	https://zenodo.org/record/4710649/files/SRR11611353_MIRNASEQ_BL.fastqsanger.gz	fastqsanger.gz
SRR11611354	BR treated miRNA	https://zenodo.org/record/4710649/files/SRR11611354_MIRNASEQ_BL.fastqsanger.gz	fastqsanger.gz
SRR1019436	Control mRNA	https://zenodo.org/record/4710649/files/SRR1019436_RNASEQ_CTL.fastqsanger.gz	fastqsanger.gz
SRR1019437	Control mRNA	https://zenodo.org/record/4710649/files/SRR1019437_RNASEQ_CTL.fastqsanger.gz	fastqsanger.gz
SRR1019438	BR treated mRNA	https://zenodo.org/record/4710649/files/SRR1019438_RNASEQ_BL.fastqsanger.gz	fastqsanger.gz
SRR1019439	BR treated mRNA	https://zenodo.org/record/4710649/files/SRR1019439_RNASEQ_BL.fastqsanger.gz	fastqsanger.gz

• 从Rules菜单中选择Add / Modify Column Definitions
  • 点击Add Definition按钮，选择Name：列A
  • 点击Add Definition按钮，选择URL：列B
• 点击Apply，然后按Upload开始上传

data_done

如上所示，为了加快分析所需的时间，本教程中样本已经过下采样以降低深度（这里仅介绍下采样方式，不需要再次对数据集进行下采样）。具体做法如下： 实践操作：数据集下采样

1. 选择使用Sub-sample sequences files（Galaxy 版本 0.2.5）工具进行下采样，参数如下：
   • 选择 “Multiple datasets”以同时处理 Collection 中的每个 fastq 文件
   • Subsampling approach：选择Take every N-th sequence (or pair e.g. every fifth sequence)
   • “N”：选择100

sub_sample_start通过这种方式，我们将随机抽取原样本中 1%的 Reads

sub_sample_done

miRNA 数据分析

当我们完成了数据导入，就可以开始研究植物甾醇类激素曝露如何改变基因表达模式了。

miRNA 读数的质量评估

由于技术限制，测序被认为是一个容易出错的过程。在 Illumina 测序平台上，替换类型的错判是主要的错误来源，并可能导致分析结果不一致。另一个可能干扰我们分析的因素是接头序列污染，这可能导致未比对的 read 数增加，因为接头是合成出的序列，其不会出现在基因组序列中。

因此，序列质量控制是所有分析中必不可少的第一步。我们将使用两种流行的工具来评估我们原始读数的质量：FastQC以及MultiQC。

注解 为了在MultiQC工具中同时可视化来自两个集合的数据，需要合并FastQC生成的结果。有关质量控制的更多信息，请参阅我们的质量控制详细教程[29]。

实践操作：初步质量检查

1. 使用以下参数运行FastQC ( Galaxy version 0.72+galaxy1) ：
   • “Dataset collection”：在项目中选择Control miRNA

fastqc

1. 将输出结果重命名为 FastQC unprocessed control miRNA: RawData 以及 FastQC unprocessed control miRNA: Webpage

fastqc_rename

1. 在BR treated miRNA Collection 上重复前面的步骤，并用类似的方式进行重命名。
2. 使用以下参数选择Merge Collections工具执行结果合并：
   • “1.Input Collections”: 选择FastQC unprocessed control miRNA: RawData
   • “2.Input Collections”: 选择FastQC unprocessed BR treated mRiNA: RawData
   • 在 “Input collections” 中：

fastqc_merge

1. 使用以下参数运行MultiQC（MultiQC ( Galaxy version 1.8+galaxy1) 工具：
   • “Which tool was used generate logs?”: 选择FastQC
   • “Dataset collection”: 选择上一步中生成的结果.
   • 在 “Results” 中:
   • “Report title”: 填入miRNA initial quality check

multiqc

1. 点击眼睛图标，在线预览生成的 HTML 文件

multiqc_download

问题 根据MultiQC提供的信息，是否需要对读数进行修剪/过滤？ MultiQC生成的报告显示三个质量参数的值超出了推荐限制：每个序列的 G/C 含量、重复序列和接头比例。

图4：GC含量图 4：miRNA 样本的每个序列 GC 含量

图5：重复序列图 5：miRNA 样本中重复序列

图6：接头含量图 6：miRNA 样本的接头含量 特别值得注意的是接头的含量，某些样本中达到了 80%。考虑到接头的丰度，如果我们消除这种污染源，预计其他参数将显示明显的改善。

为了去除接头序列污染，我们将使用Trim Galore工具，这是一个围绕**Cutadapt**[30]和FastQC的包装脚本，能基于碱基质量和接头序列对测序 Reads 进行自动化的裁剪。

实践操作：去除接头序列

1. 使用以下参数运行Trim Galore! ( Galaxy version 0.4.3.1) 工具：
   • 在 “Reads in FASTQ format” 项目中国: 选择 Control miRNA 数据集
   • “Adapter sequence to be trimmed”: 选择 Illumina small RNA adapters
   • “Is this library paired or single-end?”: 选择Single-end
2. 将输出的 collection 重命名为 Control miRNA trimmed
3. 在 BR treated miRNA 数据集上再次运行上面的步骤.

trim_miRNA

接下来，我们将重新评估序列的质量，以检查是否已经去除了接头序列。

实践操作：处理后质量检查

1. 使用以下参数运行FastQC ( Galaxy version 0.72+galaxy1) ： “Dataset collection”: 选择Control miRNA trimmed
2. 将输出重命名为 FastQC processed control miRNA: RawData 以及 FastQC processed control miRNA: Webpage
3. 使用BR treated miRNA trimmed集合重复前述步骤。
4. 使用以下参数进行Merge Collections：
   • “1.Input Collections”: 选择FastQC processed control miRNA: RawData
   • “2.Input Collections”:选择 FastQC processed BR treated miRNA: RawData
   • 在“Input collections”中:
5. 使用以下参数运行MultiQC( Galaxy version 1.8+galaxy1) ：
   • “Which tool was used generate logs?”: FastQC
   • “Dataset collection”：选择前一步生成的输出。
   • 在 “Results” 中:
   • 在 “Report title” 中：miRNA trimmed quality check
6. 点击眼睛图标，在线查看生成的 HTML 文件

经过Trim Galore处理后，通过MultiQC生成的报告评估表明，GC 含量已成功校正，并且已消除接头污染。然而，样本仍然显示出较高的重复序列（图 7）。

图7：重复序列

图 7：miRNA Reads 中重复序列的报告

问题 你认为重复序列数量过多的原因是什么？ 导致重复序列比例的两个因素是高表达的 miRNA（Seco-Cervera 等人，2018 年），miRNA 中存在高度保守的序列 motif(Glazov _et al._ 2008[31])，以及外源序列的污染。

miRNA 定量：MiRDeep2

miRNA 的定量需要使用两种不同的工具：

• MiRDeep2 Mapper工具用于对 reads 进行预处理。
• MiRDeep2 Quantifier工具用于将深度测序 reads 比对到参考序列中的 miRNA 前体上，并确定相应 miRNA 的表达。它分为两个步骤：首先，成熟 miRNA 序列比对到前体序列（可选项，还可以将卫星序列比对到前体）。第二，将测序 Reads 比对到前体序列。

实践操作：miRNA 定量

1. MiRDeep2 Mapper ( Galaxy version 2.0.0) ：
   • “Deep sequencing reads”: 选择Control miRNA trimmed
   • “Remove reads with non-standard nucleotides”: 选择Yes
   • “Clip 3’ Adapter Sequence”: 选择Don't Clip
   • “Collapse reads and/or Map”: 选择Collapse

miRdeep

1. 将输出重命名为 Collapsed control miRNA
2. 通过提供BR treated miRNA trimmed作为输入，重复上述步骤，并将其重命名为Collapsed BR treated miRNA。
3. 使用以下参数运行MiRDeep2 Quantifier ( Galaxy version 2.0.0) ：
   • “Collapsed deep sequencing reads”: 选择Collapsed control miRNA
   • “Precursor sequences”: 选择miRNA_stem-loop_seq.fasta
   • “Mature miRNA sequences”: 选择mature_miRNA_AT.fasta
   • “Star sequences”: 选择star_miRNA_seq.fasta

miRdeep_qc

1. 将输出重命名为 MiRDeep2 control miRNA 和 MiRDeep2 control miRNA (html)。
2. 通过提供MiRDeep2 BR treated miRNA作为输入，重复第四步，并将输出重命名为 MiRDeep2 BR treated miRNA 和 MiRDeep2 BR treated miRNA（html）。

为了在差异表达分析中使用MiRDeep2 Quantifier生成的输出，需要修改数据集。

实践操作：编辑 MiRDeep2 Quantifier 的输出

1. 使用以下参数运行Cut columns from a table：
   • “Cut columns”: c1,c2
   • “Delimited by”: Tab
   • “From”: `MiRDeep2

cut_col

1. 将输出重命名为 control miRNA counts
2. 使用以下参数运行Cut columns from a table：
   • “Cut columns”_：c1,c2
   • “Delimited by”_：Tab
   • “From”_：MiRDeep2 BR treated miRNA
3. 将输出重命名为 BR treated miRNA counts

miRNA 差异表达分析：DESeq2

DESeq2是一种基于负二项广义线性模型的差异基因表达分析工具。DESeq2在内部校正了文库大小的差异，因此不需要对输入数据集进行预处理归一化。

注释 最好使用每种实验条件的至少三个重复样本，以确保足够的统计功效。

实践操作：使用 DESeq2 进行差异表达分析

1. 使用以下参数运行DESeq2 ( Galaxy version 2.11.40.6+galaxy1):
   • 在 “Factor” 中:
   • “Specify a factor name, e.g. effects_drug_x or cancer_markers”: effects_brassinolide
   • 在 “Factor level” 中:
   • “Specify a factor level”: control
   • “Counts file(s)”: control miRNA counts
   • “Specify a factor level”: brassinolide
   • “Counts file(s)”: BR treated miRNA counts
   • “Insert Factor level”_
   • “Insert Factor level”_
   • “Insert Factor”_
   • “How”: Select datasets per level
   • “Files have header?”: No
   • “Choice of Input data”: Count data (e.g. from HTSeq-count, featureCounts or StringTie)

DESeq2

1. 将输出重命名为 DESeq2 results miRNA 和 DESeq2 plots miRNA
2. 点击眼睛图标，查看 DESeq2 plots miRNA 文件

DESeq2_Plot

DESeq2 生成了两个输出：一个带有归一化计数的表格，以及一份图形摘要。为了评估我们样本的相似性，我们将检查主成分分析（PCA）图。PCA 允许评估数据中最高变异性的主要方向。因此，相同条件下的样本应该聚集在一起。

图8:PCA miRNA.

图 8：来自对照组和 BR 处理组的 miRNA 样本的表达数据的 PCA 图。

如图所示，前两个主成分仅解释了总体变异的 47%和 19%。这表明植物生长素对 miRNA 调控的影响是有限的（图 8）。

过滤差异表达显著的 miRNA

接下来，我们将提取那些由于 BR 处理且在统计上显著差异表达（DE）的基因，即选择那些调整后的 p 值小于或等于 0.05 的基因。0.05 的阈值表示假阳性结果的概率小于 5%。

p 值是衡量观察到的差异可能仅由随机机会引起的概率的指标。较小的 p 值表明，如果没有真实差异存在，获得当前数据的可能很小。0.05 的 p 值阈值表示结果是假阳性的概率为 5%。p-adj（也称为 q 值）是一个调整后的 p 值，考虑到假发现率（FDR）。当我们从一个小样本集中检测成千上万的变量时，应用 FDR 是必要的。

实践操作：过滤差异表达的 miRNA

1. 使用Filter（Galaxy Version 1.1.1）对结果表格进行过滤，具体参数如下：
   • “Filter”: DESeq2 results miRNA
   • “With following condition”: c7<0.05

filter_DEseq

empty_filter

问题 有多少 miRNA 显示出统计上显著的差异表达？ 不幸的是，我们没有检测到任何差异表达的 miRNA。这是下采样数据集没有足够的数据来进行差异表达检测导致的。

为了获得合理的结果，我们需要分析完整数据集。您可以按照上述教程使用完整数据集进行分析，也可以将我们从完整数据集生成的 DESeq2 分析结果导入到您的历史数据中。

实践操作：获取完整 miRNA 数据集上的 DESeq2 分析结果

从 Zenodo 导入文件：

• 打开 上传 菜单
• 单击 Paste/Fetch 按钮
• 复制 Zenodo 链接并点击 Start https://zenodo.org/record/4663389/files/miRNA_DESeq2_results_complete_dataset.tabular

upload_complete_miSeq

1. 根据样本 ID（例如 miRNA_DESeq2_results_complete_dataset）重命名每个数据集
2. 为 miRNA 数据添加相关的标签：#miRNA #BR #control

重复导入的 DESeq2 结果上的过滤步骤。

实践操作：从完整数据集中筛选和排序差异表达的 miRNA

1. 使用Filter(Galaxy Version 1.1.1)对于数据进行筛选，参数如下：
   • 参数文件 _“Filter”_：miRNA_DESeq2_results_complete_dataset
   • “With following condition”_：c7<0.05

filter_p_complete

1. 将输出重命名为 差异表达的miRNA
2. Filter 数据，Filter在任何列上进行筛选（Galaxy 版本 1.工具进行 1.，：
   • 参数文件 _Filter_：差异表达的miRNA
   • “With following condition”_：c3>0

filter_upgrade

1. 将输出重命名为 Upregulated miRNAs
2. 使用 Sort(Galaxy Version 1.1.0) 对数据进行排序，参数如下：
   • “Sort Dataset”_：Upregulated miRNAs
   • “on column”: Column: 3
   • “everything in”: Descending order

sort_order

问题

1. 有多少 miRNA 表达差异显著？
2. 有多少 miRNA 表达上调且具有统计学意义，上调最多的 miRNA 是什么？
3. 在 442 个 miRNA 中，有 39 个有显著的表达差异。
4. 有 16 个上调的 miRNA，其中Ath-miR156g上调幅度最大。

mRNA 数据分析

在完成 miRNA 的差异表达分析之后，下一阶段是分析 mRNA 表达如何响应芸苔素的变化。

mRNA 读取质量评估

与前一节相同，我们首先将评估我们序列的质量。

实践操作：初始质量检查

1. FastQC ( Galaxy version 0.72+galaxy1)，使用以下参数：
   • “Dataset collection”_：Control mRNA

fastQC_mRNA

1. 将输出重命名为 FastQC unprocessed control mRNA: RawData 和 FastQC unprocessed control mRNA: Webpage
2. 重复上一步操作，使用 BR treated mRNA collection。
3. Merge Collections，使用以下参数：
   • “1.Input Collections”: FastQC unprocessed control mRNA: RawData
   • “2.Input Collections”: FastQC unprocessed BR treated mRNA: RawData
   • 在 “Input collections” 中：
4. 使用MultiQC ( Galaxy version 1.8+galaxy1) ，参数如下:
   • “Which tool was used generate logs?”: FastQC
   • “Dataset collection”: select the output generated in the previous step.
   • 在 “Results” 中:
   • “报告标题”_：mRNA初始质量检查

multiQC_mRNA

1. 单击眼睛图标，检查生成的 HTML 文件

问题 是否有任何统计数据表明需要处理样本以改善其质量？ 所有统计数据都在可接受范围内。然而，接头含量显示我们的 reads 中存在 Illumina 通用接头，可以去除以避免可能在后期阶段的干扰（见图 10）。

图10：mRNA样本的接头含量。图 10：mRNA 样本的质量评估

注释 尽管我们的样本有接头序列，但在这种情况下我们不会对其进行裁剪。我们将在下一节中解释原因。

基因表达的定量：Salmon

在完成 Read 的质量评估之后，我们可以继续对基因表达进行定量。这一步的目标是确定每个 Read 来自哪个转录本以及与每个转录本相关联的 Read 总数。在本教程中，我们将使用 Salmon 工具(Patro _et al._ 2017[32])来对 mRNA 转录本定量。

fig11:Salmon标志。

Salmon的一个特点是它不需要执行基于碱基的比对，这是工具如STAR和HISAT2的耗时步骤。Salmon 依赖于 quasi-mapping 概念，这是一种新的比对技术，可以快速而准确地将 RNA-SeqRead 比对到目标转录组。与标准比对不同，quasi-mapping 试图找到每个 Read 的最佳比对，并通过在目标和查询位置之间找到最小集合的动态大小的、右极大的匹配上下文来实现(Srivastava _et al._ 2016[33])。

Salmon的 quasi-mapping 方法需要一个参考索引来确定准确比对之前的位置和方向信息。它允许以一种优化转录本识别和定量使用的格式提供转录组。

在确定每个 Read 的最佳比对后，Salmon在建模特定于样本的参数和偏差后生成最终的转录本表达量估计值。比对到多个转录本的 Read 将在所有比对之间分配计数，从而避免了对不同基因异构体的信息丢失。

quasi-mapping 算法利用了两个主要的数据结构，转录组 T 的广义后缀数组（SA）和一个哈希表（h），将 T 中的每个 k-mer 比对到其 SA 区间（默认 k = 31）。在 quasi-mapping 过程中，从左到右扫描 Read，直到遇到一个出现在 h 中的 k-mer（ki，从 Read 的位置 i 开始）。在哈希表中查找 ki，并检索 SA 间隔，给出包含 k-mer ki 的所有后缀。然后，通过找到与参考后缀匹配的 Read 的最长子字符串来计算最大可比对前缀（MMPi）。由于测序错误，MMPs 可能无法跨越完整的 Read。在这种情况下，基于 MMPi 的 SA 间隔的最长公共前缀确定下一个信息位置（NIP）。从 NIP 之前的 k 个碱基继续后缀数组搜索，直到最终 NIP 小于 Read 末尾的 k 个碱基之前。最后，对于每个 Read，该算法报告了它比对到的转录本、位置和链信息(Srivastava _et al._ 2016[34])。

fig12：quasi-mapping算法。

图 10：使用 k=3 进行 quasi-mapping 的演示。ki 的哈希表查找返回后缀数组区间[b, e)。读取中位置 6 处的碱基'G'是测序错误的结果。因此，MMPi 是'ATTGA'，MMPi 的 SA 区间是[b'，e')。下一个 k-mer 从 NIP 之前的 k 个碱基开始，这是区间[b'，e')的最长公共前缀之后的基。最后，在上述示例中，读取很可能比对到后缀数组[e'-1, e')。

实践操作：使用 Salmon 量化基因表达

1. Salmon quant( Galaxy version 0.14.1.2+galaxy1)，使用以下参数：
   • “Select a reference transcriptome from your history or use a built-in index?”: Use one from the history
   • 在 “Data input” 中:
   • “Validate mappings”: Yes
   • “Transcripts fasta file”: transcriptome.fasta
   • 在 “Salmon index” 中:
   • “FASTQ/FASTA file”: Control mRNA
   • “Is this library mate-paired?”: Single-end
   • “Select salmon quantification mode:”: Reads
   • “File containing a mapping of transcripts to genes”: annotation_AtRTD2.gtf

salmon_quant

1. 将输出重命名为 Salmon control mRNA (Quantification) 和 Salmon control mRNA (Gene Quantification)
2. 通过使用BR treated mRNA数据集重复上述过程

注释：quasi-mapping 序列要求 使用此方法时不需要裁切 Read，因为如果 Read 中有不在哈希表中的 k-mer，则不会计数。Read 的量化结果仅取决于参考转录组的质量。

Salmon为每个输入样本生成两个输出：

• Quantification：按转录本汇总的量化结果
• Gene quantification：按基因汇总的量化结果

每个输出包含一个具有五列的表格数据集：

问题 你能解释为什么我们之前没有裁切 Read 吗？ 原因是，由于包含接头序列的 k-mer 不出现在转录组中（从中生成哈希表），它们被忽略了。

mRNA 的差异表达分析：DESeq2

现在，让我们分析哪些基因在对油菜素内酯的反应中显示出统计学上显著的差异表达。

实践操作：使用 DEseq2 进行差异表达分析

1. DESeq2( Galaxy version 2.11.40.6+galaxy1)，使用以下参数：
   • “Program used to generate TPMs”: Salmon
   • “Gene mapping format”: GTF/GFF3
   • “GTF/GFF3 annotation file”: annotation_AtRTD2.gtf
   • 在 “Factor” 中:
   • “Specify a factor name, e.g. effects_drug_x or cancer_markers”: effects_brassinolide
   • 在 “Factor level” 中:
   • “Specify a factor level”: control
   • “Counts file(s)”: Salmon control mRNA (Gene Quantification)
   • “Specify a factor level”: brassinolide
   • “Counts file(s)”: Salmon BR treated mRNA (Gene Quantification)
   • “Insert Factor level”_
   • “Insert Factor level”_
   • “Insert Factor”_
   • “How”: Select datasets per level
   • “Choice of Input data”: TPM values (e.g. from kallisto, sailfish or salmon)

DESeq_mRNA

1. 将输出重命名为 DESeq2 results mRNA 和 DESeq2 plots mRNA

fig13：mRNA样本的PCA图。

图 11：来自控制和 BR 处理样本的差异表达数据的 PCA 图。

问题 从 PCA 图中可以得出关于样本相似性的什么结论？ 从图提供的信息可以得出结论，即属于相同实验条件的样本之间存在高度相似性，因为第一主成分（x 轴）能够解释 81%的差异，并且样本位于 x 轴的两侧。

过滤显著差异表达的 mRNA

为了完成 mRNA 差异表达的分析，我们将提取那些在对油菜素内酯的反应中显示出显著表达诱导的转录本。在继续进行进一步分析之前，类似于 miRNA 数据分析，导入从完整 mRNA 数据集生成的 DESeq2 结果。

实践操作：检索完整 mRNA 数据集上的 DESeq2 分析结果

从 Zenodo 导入文件：

• 点击 upload 菜单
• 点击 Paste/Fetch 按钮
• 复制 Zenodo 链接并按“Start” https://zenodo.org/record/4663389/files/mRNA_DESeq2_results_complete_dataset.tabular

根据样本 id 重命名每个数据集（例如 mRNA_DESeq2_results_complete_dataset）

添加所有与 mRNA 数据分析相关的标签：#mRNA #BR #control

现在我们继续进行差异表达基因分析。

实践操作：提取显著差异表达基因

1. Filter(Galaxy Version 1.1.1)工具进行过滤，参数如下：
   • “Filter”: mRNA_DESeq2_results_complete_dataset
   • “With following condition”: c7<0.05
2. 将其重命名为 Differentially expressed mRNAs

mRNA_filter_p

1. Filter(Galaxy Version 1.1.1)工具进行过滤，参数如下：
   • “Filter”_：Differentially expressed mRNAs
   • “With following condition”_：c3>1

mRNA_upgrade

1. 将其重命名为 Upregulated mRNAs
2. Filter(Galaxy Version 1.1.1)工具进行过滤，参数如下：
   • “Filter”_：Differentially expressed mRNAs
   • “With following condition”_：c3<-1

mRNA_downgrade

1. 将其重命名为 Downregulated mRNAs

问题

1. 有多少基因显示出统计学上显著差异？
2. 其中有多少基因显著上调（至少两倍变化）？下调呢？
3. 最显著的差异表达下调基因是什么，其生物功能是什么？
4. 在两种实验条件之间表达差异的基因总数为 4176。
5. 其中，有 328 个基因在 BR 处理下显著下调，778 个基因上调。
6. 最显著的差异表达基因是 AT3G30180（BR60X2）。BR60X2 编码一种细胞色素 p450 酶，该酶催化油菜素内酯的生产中的最后一步反应。它能够将 6-去氧铸杆酮转化为铸杆酮，进行 C-6 氧化，并将铸杆酮进一步转化为油菜素内酯（来源：TAIR database[35]）。

miRNA 靶标的识别

为了预测哪些 miRNA 靶向哪些 mRNA，首先我们需要它们的转录组序列，以 FASTA 格式。现在我们将获取由油菜素内酯诱导的 miRNA 序列。

实践操作：获取上调 miRNA 的序列 注释 以下工具可以在Text Manipulation和Filter and Sort部分找到。

1. Cut columns from a table，使用以下参数：
   • “Cut columns”: c1
   • “Delimited by”: Tab
   • “From”: Upregulated miRNAs

cut_upgrade_miRNA_id

1. 将输出重命名为 Upregulated miRNA ids
2. 使用Filter FASTA( Galaxy version 2.1)，参数见下：
   • “List of IDs to extract sequences for”: Upregulated miRNA ids
   • “Match IDs by”: Default: ID is expected at the beginning: >ID
   • “FASTA sequences”: mature_miRNA_AT.fasta
   • “Criteria for filtering on the headers”: List of IDs

extract_mature_miRNA

1. Filter FASTA ( Galaxy version 2.1)使用以下参数：
   • “List of IDs to extract sequences for”: Upregulated miRNA ids
   • “Match IDs by”: Default: ID is expected at the beginning: >ID
   • “FASTA sequences”: star_miRNA_seq.fasta
   • “Criteria for filtering on the headers”: List of IDs

extract_star_miRNA

1. Concatenate datasets tail-to-head (Galaxy Version 1.0.0)使用以下参数：
   • 插入数据集
   • “Select”: output of Filter FASTA tool on star_miRNA_seq.fasta
   • “Concatenate Dataset”: output of Filter FASTA tool on mature_miRNA_AT.fasta

cat_miRNA

1. 将其重命名为 Upregulated miRNA sequences
2. 点击眼睛图标，检查Upregulated miRNA sequences文件

为了识别上调 miRNA 的潜在靶标，有必要获取 FASTA 格式的所有下调 mRNA 序列。

实践操作：获取下调 mRNA 的基因序列

1. Cut columns from a table，使用以下参数：
   • “Cut columns”: c1
   • “Delimited by”: Tab
   • “From”: Downregulated mRNAs

cut_downgrade_mRNA_id

1. 将输出重命名为 Downregulated mRNA ids
2. Filter FASTA ( Galaxy version 2.1)使用以下参数：
   • *“List of IDs to extract sequences for: Downregulated mRNA ids
   • “Match IDs by”: Custom regex pattern
   • “Regex search pattern for ID”: GENE=(AT.{7})
   • “FASTA sequences”_：transcriptome.fasta（输入数据集）
   • “Criteria for filtering on the headers”: List of IDs

extract_mRNA

1. 将其重命名为 Downregulated mRNA sequences
2. 点击眼睛图标，检查Downregulated mRNA sequences文件

使用TargetFinder进行 miRNA 靶标预测

我们现在准备开始搜索 miRNA 靶标基因。为此，我们将使用TargetFinder工具，该工具用于植物中 miRNA 靶标预测(Srivastava _et al._ 2014[36], Ma _et al._ 2017[37])。

实践操作：识别潜在的 miRNA 靶标

1. TargetFinder ( Galaxy version 1.7.0+galaxy1) 使用以下参数：
   • “Input small RNA sequences file”: Upreguled miRNA sequences
   • “Target sequence database file_”: Downregulated mRNA sequences
   • “Prediction score cutoff value”: 5.0
   • “Output format”: Tab-delimited format

TargetFinder

1. 点击眼睛图标，检查TargetFinder的输出。

恭喜！我们已经识别出以下 5 个潜在的参与 brassinosteroid-miRNA 调控网络的基因。

fig15:TargetFinder结果。

图 12：TargetFinder 结果。

最后，我们可以访问 TAIR 数据库中关于所识别基因的所有信息：

• AT5G10180[38]：拟南芥硫酸盐转运蛋白 68，AST6
• AT3G09220[39]：LAC7，漆酶 7
• AT2G46850[40]：蛋白激酶超家族蛋白
• AT5G64260[41]：EXL2，EXORDIUM LIKE 2
• AT3G63200[42]：类脂质转移蛋白 9，PLA IIIB

AT4G14365 和 AT1G26890 都是特征不太明确的基因。对于 AT5G50570，实验研究表明，该基因通过 miR156 介导的信号通路参与了Medicago sativa的耐涝性(Feyissa _et al._ 2021[43])。另一方面，根据Gao _et al._ 2017[44]，SPL13 可能是植物营养生长的负调控因子。miR156 与 SPL 转录因子之间的相互作用已在多个禾本科家族成员中被提出(Yue _et al._ 2021[45])。

注释：假设检验 我们可以根据我们的结果提出一个假设：抑制 AT2G46850 基因可以使植物具有更好的抗旱性。有可能验证它吗？是的！我们可以这么做：获取AT2G46850 mutant seeds[46]和wild type seeds[47]，在两种控制条件下培养：浇水和干旱胁迫，并在 33 天后分析植物重量（图 13）。

fig16:植物培养。图 13：拟南芥生长条件(de Ollas _et al._ 2019[48])。

可选练习

作为额外内容，您可以尝试使用 NCBI GEO 数据库中存储的序列，使用访问号GSE119382来重复工作流程。在这种情况下，我们将比较在两种实验条件下过表达油菜素内酯受体 BRL3 的突变体的基因表达模式：对照和干旱胁迫。所需的数据集在数据库中可用：

实践操作：从数据库导入数据

进入Shared data（顶部面板）并点击Data Libraries

在搜索框中输入以下标识符：4710649

选择以下文件：

https://zenodo.org/record/4710649/files/SRR7779222_BRL3_mRNA_drought.fastqsanger.gz
https://zenodo.org/record/4710649/files/SRR7779223_BRL3_mRNA_drought.fastqsanger.gz
https://zenodo.org/record/4710649/files/SRR7779224_BRL3_mRNA_drought.fastqsanger.gz

点击导出到历史记录和作为一个集合

在弹出窗口中点击继续

为其命名为BRL3 mRNA drought，然后点击创建列表

使用剩余的文件重复之前的程序：

https://zenodo.org/record/4710649/files/SRR7779228_BRL3_mRNA_watered.fastqsanger.gz
https://zenodo.org/record/4710649/files/SRR7779229_BRL3_mRNA_watered.fastqsanger.gz

最后为其命名为BRL3 mRNA control并点击创建列表

我们将使用前一次分析中获得的上调 miRNA 来识别潜在的靶标。

Upregulated miRNA   https://zenodo.org/record/4710649/files/upregulated_miRNA_BR_complete_dataset.fasta

结论

在这个教程中，我们分析了 RNA 测序数据以提取有关潜在受油菜素内酯调控的基因的信息。为此，选择的方法是识别与油菜素内酯响应上调的 miRNA 互补的基因。最终结果已经鉴定出五个潜在的 miRNA 靶标。