文末有赠书福利,快来参与吧!
导语
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分析的下游的绘图需要。
功能示例
01
R包安装和加载
devtools::install_github("THERMOSTATS/RVA_prod")
#GitHub安装RVA_prod没成功,可能是版本失效了,可以尝试
devtools::install_github("THERMOSTATS/RVA")
library(RVA)
02
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")
)
03
QQ图绘制
#多组QQ图绘制
plot_qq(data = d1,comp.names = c("A", "B"))
04
绘制火山图
#定义上调基因
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)
05
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]
06
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 #类型
07
基因表达绘图
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")
小编总结
这个R包非常实用,用起来是比较简单方便,对于使用者来说还是很友好的,图也算漂亮,所以还等什么呢?小伙伴们快快用起来吧?