不仅仅是火山图，你可以获得更多可视化结果！

作图丫

发布于 2022-03-29 13:59:06

5150

发布于 2022-03-29 13:59:06

文章被收录于专栏：作图丫

文末有赠书福利，快来参与吧！

导语

GUIDE ╲

今天小编给大家带来一个很好用的RNA-seq 可视化的R包-RVA( RNAseq Visualization Automation)。“RVA”是一个功能集合，可有效地可视化RNAseq差异表达的分析结果，并利用Fisher精确测试方便有效地评估基因集或通路富集。该包用于RNA-seq分析中的下游可视化和通路富集分析真的是很实用和方便了。

R包链接 https://github.com/THERMOSTATS/RVA

功能介绍

RVA包能够实现cutoff预置下保留差异基因数目的可视化，并且能够轻松绘出QQ图、火山图、基因富集分析以及基因表达热图、boxplot。整体来说满足了RNA-seq分析的下游的绘图需要。

功能示例

R包安装和加载


devtools::install_github("THERMOSTATS/RVA_prod")
#GitHub安装RVA_prod没成功，可能是版本失效了，可以尝试
devtools::install_github("THERMOSTATS/RVA")
library(RVA)

cutoff阈值的差异基因数目展示，方便观察保留的差异基因数目

#加载测试数据
df <- RVA::Sample_summary_statistics_table
df1 <- RVA::Sample_summary_statistics_table1 
d1 <- list(df, df1)

‍

每一行为一个基因，基因ID可以是ACCNUM, ALIAS, ENSEMBL, ENSEMBLPROT, ENSEMBLTRANS, ENTREZID, ENZYME, EVIDENCE, EVIDENCEALL, GENENAME, GO, GOALL, IPI, MAP, OMIM中的一种。

plot_cutoff(data = df,
  comp.names = NULL,
  FCflag = "logFC",#Foldchange所在的列名
  FDRflag = "adj.P.Val",#FDR所在的列明
  FCmin = 1.2,#Foldchange最小值
  FCmax = 2,#Foldchange最大值
  FCstep = 0.1,#展示Foldchange的步长
  p.min = 0,#FDR的最小值
  p.max = 0.2,#FDR的最大值
  p.step = 0.01,#展示Foldchange的步长
  plot.save.to = NULL,
  gen.3d.plot = TRUE,
  gen.plot = TRUE)

此Bar图可以清楚看出每个阈值下保留的差异基因的数目，方便进行决策。

同时也可多组差异基因进行绘图展示。

plot_cutoff(data = d1,
            comp.names = c("A", "B")
            ）

QQ图绘制

#多组QQ图绘制
plot_qq(data = d1,comp.names = c("A", "B"))

绘制火山图

#定义上调基因
up<-row.names(df)[which((df$logFC>=1) & (df$adj.P.Val<=0.05))]
#定义下调基因
down<-row.names(df)[which((df$logFC<=-1) & (df$adj.P.Val<=0.05))]
plot_volcano(data = df,highlight.1=up,
             highlight.2=down,             
             upcolor = "#FF0000",            
             downcolor = "#0000FF",            
             xlim = c(-2,2),            
             ylim = c(0,10),highlight.FC.cutoff=4)

这个有个可能是bug，highlight.FC.cutoff=2/4/16等数是highlight线是能对应上的，不过其他数测试时有问题，有兴趣的小伙伴可以测试下，可以反馈给我，让我也明白下！

同时也可对疾病相关基因进行标记，命令行如下：

plot_volcano(data = d1,
             geneset=dgs,
             comp.names = c('a', 'b'), 
             upcolor = "#FF0000",
             downcolor = "#0000FF", 
             xlim = c(-1.5,1.5),
             ylim = c(0,10),
             highlight.FC.cutoff=4)

Pathway富集分析

plot_pathway(
  data = df,
  comp.names = NULL,#分组定义
  gene.id.type = "ENSEMBL",#基因ID类型
  FC.cutoff = 1.3, #Foldchange cutoff阈值
  FDR.cutoff = 0.05, #FDR cutoff阈值
  FCflag = "logFC",
  FDRflag = "adj.P.Val", 
  Fisher.cutoff = 0.1, #富集cutoff阈值
  Fisher.up.cutoff = 0.1, #上调基因富集cutoff阈值
  Fisher.down.cutoff = 0.1, #下调基因富集cutoff阈值
  plot.save.to = NULL, #保存图
  pathway.db = "rWikiPathways" # 数据库选择
  )

简单的运行一下：

pathway.result <- plot_pathway(data = df, pathway.db = "KEGG", gene.id.type = "ENSEMBL")

pathway.result[[1]] #上下调基因的富集table
pathway.result[[2]] #所有差异基因的富集table

#上下调基因的富集bar图
pathway.result[[3]]

##所有基因的富集bar图
pathway.result[[4]]

#上下调/所有差异基因的富集merge bar图
pathway.result[[5]]

同时类似的，这个富集也是可以同时进行多组差异基因的富集分析。

list.pathway.result <- plot_pathway(data = list(df,df1),
                       comp.names=c("A","B"),
                       pathway.db = "KEGG", 
                       gene.id.type = "ENSEMBL")

#差异基因富集绘图
list.pathway.result[4]

#上下调基因富集绘图
list.pathway.result[3]

heatmap 绘图

count <- RVA::count_table[,1:50]

count[1:6,1:5]

annot <- RVA::sample_annotation[1:50,]
head(annot)

hm.expr <- plot_heatmap.expr(data = count, #count 文件
                             annot = annot, #注释文件
                             sample.id = "sample_id", #样品ID
                             annot.flags = c("day", "Treatment"), #按照day 和 Treatment 进行分组注释
                             ct.table.id.type = "ENSEMBL", #输入的基因ID类型
                             gene.id.type = "SYMBOL", #选择展示的基因ID类型
                             gene.names = NULL, #选择展示的基因ID
                             gene.count = 10, #展示的基因数量
                             title = "RVA Heatmap", #title 名
                             fill = "CPM", count 计数方法
                             baseline.flag = "day", #当fill = "CFB"时起作用 
                             baseline.val = "0", #当fill = "CFB"时起作用
                             plot.save.to = NULL, #保存
                             input.type = "count") count #类型

基因表达绘图

anno <- RVA::sample_annotation
head(anno)

ct <- RVA::sample_count_cpm
ct[1:6,1:5]

gene.result <- plot_gene(ct, 
               anno,
               gene.names = c("AAAS", "A2ML1", "AADACL3"),
               ct.table.id.type = "ENSEMBL",
               gene.id.type = "SYMBOL",
               treatment = "Treatment",
               sample.id = "sample_id",
               time = "day",
               log.option = TRUE,
               plot.save.to = NULL,
               input.type = "cpm")