优秀的R包快速解决RNA-Seq下游分析

大家好，我是阿琛。关于RNA测序，差异分析之后，下游的可视化分析往往是一个十分关键而又有一定难度的内容。今天来给大家介绍一个集RNA-Seq的下游可视化和通路分析于一体的新生代R包，即RVA包。

RVA包，RNAseq Visualization Automation，最初发表于GitHub上，于2020年12月被CRAN收录发表；它可以快速地对差异分析结果进行汇总和可视化，并能在一定程度上评估相关基因集或通路的富集情况。

1.R包安装与数据准备

1.1 R包的安装

#install.packages("RVA")library(RVA)

由于该R包已经被CRAN收录，所以直接使用install.packages()代码即可完成。同时，对于安装过程中缺少的相关依赖包，大家可以根据提示，自己手动进行安装。

1.2 读取示例数据

df <-RVA::Sample_summary_statistics_tabledf1 <-RVA::Sample_summary_statistics_table1df[1:3,]df1[1:3,]

随后，读取RVA自带的示例数据集。可以看到，两个用于输入分析的数据集均是经过R的DEseq2或limma包差异分析后的统计结果表格，包含了差异倍数（logFC），P值（P.Value），校正后的P值（adj.P.Val）等等。

d1 <- list(df, df1)

随后，将两组分析结果汇总成一个列表（list），作为后续研究分析的输入数据。在该R包分析过程中，除了输入单个分析结果的数据框（data frame）外，还可以将多个分析结果组合成列表后输入。

2.Cutoff Plot

Cutoff Plot可以分析不同截断值组合情况下差异表达基因的数量，即评估具有不同FDR和倍数变化截止值时差异表达基因的数量。

2.1 以数据框形式输入

cutoff.result <-plot_cutoff(data = df,gen.plot = TRUE,gen.3d.plot = TRUE)

在输出的结果中，一共包含了3个不同板块的结果内容：

1. 表格总结了不同阈值组合下差异表达基因数目

head(cutoff.result[[1]])

可以看到，在该表格中展示了在不同的pvalue和FC组合下所得到的差异表达基因数量。

2. 通过3D图形展示分析结果

cutoff.result[[2]]

结果显示：

3. 使用2D图形展示分析结果

cutoff.result[[3]]

结果显示：

2.2 以列表形式输入

cutoff.result.list <- plot_cutoff(data = d1, comp.names = c('a', 'b'))

当然，我们也可以将两组分析结果以列表的形式进行输入，并将两组分别命名为a和b；同样，我们可以来查看一下最终的分析结果：

1. 以表格的形式汇总了整体的分析结果

head(cutoff.result.list[[1]])

可以看到，在结果中，根据分组名称，分别对不同阈值下的差异表达基因数量进行了汇总分析；

2. 对结果进行可视化展示：与之前结果不同的是，该结果无法以3D的形式进行展示

cutoff.result.list

结果显示：

3.火山图

3.1 以数据框形式输入

plot_volcano(data = df)+ ggplot2::theme_bw()

对于火山图的绘制，直接使用plot_volcano()函数即可快速进行绘制；除了数据框形式输入外，我们也可以尝试将列表d1进行输入绘制。而且，由于整个绘图是基于ggplot2包而构建的，因此我们可以根据需要对其结果进一步进行修改。

3.2 对感兴趣的基因进行标记

除了使用平时标记的显著上下调基因赋予颜色外，在这个R包中还提供了另外一种方法。在RVA包中，存在一个疾病相关基因集合，我们可以对这个集合和差异分析结果进行取交集，对共同的致病基因进行标记，查看他们的上调或下调变化水平。

#疾病相关基因集合dgs <- RVA::Sample_disease_gene_set

随后，将疾病基因集赋值给参数geneset，即可进行绘制。

plot_volcano(data = df,geneset = dgs,upcolor = "#FF0000",downcolor = "#0000FF",xlim = c(-3,3),ylim = c(0,14))

可以看到，在火山图中，两者的共同基因被标记了出来，而且在疾病中上调基因为红色，下调基为蓝色，结合火山图中该点的位置，可以快速的进行上下调关系的比对。

4.热图

在RNA-Seq下游的可视化分析过程中，除了火山图，通过热图展示整体的结果内容同样是一个十分常见的内容。在此，可以使用plot_heatmap.expr()函数来进行热图的绘制。而且，plot_heatmap.expr()函数主要基于ComplexHeatmap包的内容来进行构建，得到的图形还是比较美观的。

count <-RVA::count_tableannot <-RVA::sample_annotationcount[1:3,1:3]str(annot)

首先，读取测序样本的Count数据以及样品的相关处理信息。可以看到，在表达数据中，行名为基因的Ensembl id号，列名为样本名；而在处理信息中，包括了样本信息（Sample id），组织类型（Tissue），时间（Day），以及处理方法（Treatment）等信息内容。

plot_heatmap.expr(data = count, annot = annot, sample.id = "sample_id", annot.flags = c("tissue","day", "Treatment"), ct.table.id.type = "ENSEMBL", gene.id.type = "SYMBOL", gene.names = NULL, gene.count = 20, title = "RVA Heatmap", fill = "CPM", baseline.flag = "day", baseline.val = "0", plot.save.to = NULL, input.type = "count")

随后，将表达数据count赋值给data参数，将样品信息annot赋值给annot参数，设定样本名为“sample_id”，同时设定输入基因名为ENSEMBL，在热图中输出为SYMBOL形式，并展示差异最显著的前20个基因信息。这样，一张精美的热图就绘制完成了。