RNA-seq 数据中可变剪接与多聚腺苷酸化位点的鉴定

摘要

可变剪接（AS）和可变多聚腺苷酸化（APA）能够扩展转录本异构体及其产物的多样性。本文中，我们描述了用于分析整体RNA测序（bulk RNA-seq）和3'端测序实验的生物信息学流程，以检测和可视化不同实验条件下可变剪接与可变多聚腺苷酸化的差异。

摘要

除了常规的RNA-Seq分析用于检测不同实验/生物学条件下基因表达差异（DGE）外，RNA-seq数据还可用于探索外显子水平上的其他复杂调控机制。可变剪接和多聚腺苷酸化通过产生不同的转录本，在转录后水平调控基因表达，从而在基因功能多样性中发挥关键作用；若仅局限于全基因水平的分析，可能会遗漏这一重要的调控层面。本文中，我们利用Bioconductor及其他软件包和函数（包括DEXSeq、Limma包中的diffSplice以及rMATS），逐步演示如何在不同条件下识别和可视化外显子使用差异及多聚腺苷酸化位点的使用差异。

引言

多年来，RNA-seq 通常被用于估计差异基因表达和基因发现¹。此外，它还可用于估计由于基因表达不同异构体而导致的外显子水平使用变化，从而有助于更好地理解转录后水平的基因调控。大多数真核生物基因通过可变剪接（AS）产生不同的异构体，以增加 mRNA 表达的多样性。AS 事件可分为不同的模式：完全外显子跳跃（SE），即一个“盒式”外显子及其侧翼内含子完全从转录本中移除；当外显子两端存在两个或多个剪接位点时，发生可变（供体）5' 剪接位点选择（A5SS）和可变 3'（受体）剪接位点选择（A3SS）；内含子保留（RI），即内含子保留在成熟 mRNA 转录本中；以及外显子使用的相互排斥（MXE），即两个可用外显子中每次仅保留其中一个^2,3。可变多聚腺苷酸化（APA）也通过使用不同的多聚（A）位点从单个转录本生成多种 mRNA 异构体，在调控基因表达中发挥重要作用⁴。大多数多聚腺苷酸化位点（pAs）位于 3' 非翻译区（3' UTRs），产生具有不同 3' UTR 长度的 mRNA 异构体。由于 3' UTR 是识别调控元件的核心区域，不同的 3' UTR 长度可能影响 mRNA 的定位、稳定性和翻译效率⁵。目前已有一类优化用于检测 APA 的 3' 端测序检测方法，它们在实验方案细....

方案

1. 分析中使用的工具和R包的安装

1. Conda 是一款流行且灵活的软件包管理器，可在所有平台上便捷地安装软件包及其依赖项。请使用 'Anaconda' （conda 包管理器）用于安装 'conda' 可用于安装分析所需的工具/软件包。
2. Download 'Anaconda' 根据系统要求从 https://www.anaconda.com/products/individual#Downloads 下载，并按照图形安装程序中的提示进行安装。使用 ‘conda’ 在 Linux 命令行中输入以下命令。
   conda install -c daler sratoolkit
conda install -c conda-forge parallel
conda install -c bioconda star bowtie fastqc rmats rmats2sashimiplot samtools fasterq-dump cutadapt bedtools deeptools
3. 要下载本方案中使用的所有 R 软件包，请在 R 控制台（通过在 Linux 命令行中输入 `R` 启动）中键入以下代码。 'R') or Rs....

结果

完成上述逐步工作流程后，AS 和 APA 分析的输出结果及代表性数据将以表格和数据图的形式呈现，生成方式如下。

AS:
AS分析的主要输出结果（补充表1 用于差异剪接分析； 表2 对于 DEXSeq）是一份显示在不同条件下外显子使用差异的列表，以及一份其一个或多个组成外显子表现出显著整体剪接活性的基因列表，按统计显著性排序。 补充表1第2个标签页显示显著外显子，各列分别展示外显子相对于其余部分的差异倍数变化（FC）、单个外显子未经校正的p值以及经Benjamini-Hochberg校正后的p值。基于校正后p值设定阈值，可获得一组具有明确定义错误发现率（FDR）的外显子。 补充表1，表 3 显示了在总体剪接活性方面具有显著性的基因排序列表，其中一列展示了使用 Simes 方法计算的基因水平校正 p 值。类似的数据如 表2 用于 DEXSeq。 补充图1

讨论

在本研究中，我们评估了基于外显子和基于事件的方法，用于在整体RNA测序（bulk RNA-Seq）和3'端测序数据中检测可变剪接（AS）和可变多聚腺苷酸化（APA）。基于外显子的AS分析方法可同时提供差异表达外显子的列表，以及按基因水平上整体剪接活性差异的统计显著性排序的基因排名（表1-2、4-5）。对于diffSplice软件包，其通过在外显子水平拟合加权线性模型，估算某一外显子相对于同一基因内其他外显子平均对数倍数变化的差异（称为每个外显子的FC），从而确定差异使用情况。基因水平的统计显著性通过Simes方法将各个外显子水平的显著性检验结果合并为基因水平的综合检验来计算¹⁰。这种基于基因水平剪接活性差异的排序随后可用于对相关关键通路进行基因集富集分析（GSEA）¹⁰。DEXSeq采用类似策略，通过拟合广义线性模型来测量外显子的差异使用，但在过滤、标准化和离散度估计等步骤上有所不同。在比较DEXSeq和DiffSplice分析结果中排名前500位的显示AS活性的外显子和APA位点时，我们发现分别有310个外.......