专栏首页生信宝典R语言一键批量完成差异统计和可视化

R语言一键批量完成差异统计和可视化

撰文:文涛 南京农大

责编:刘永鑫 中科院遗传发育所

R语言一键完成差异检测从数据到展示

单因素差异分析的完整方案

关键词:正态性检验;方差齐性;非参数检验;秩和检验;多重比较;带显著性字母柱状图或箱线图

由于作者水平有限,大家可以添加我的个人微信讨论细节、bug和可改进的地方(微信号:nanjingxuezi)

方案优点

这份方案有一下几个特征,在展示之前给大家mark一下

  • 完整的差异分析思路及其R语言实现
  • 提供两种可视化方案:柱状图箱线图;差异结果使用两种表示方法:字母进行排序表示,或者两组之间连线。
  • 可视化图形标注方差齐性结果,正态性检验结果等
  • 自动识别数据指标,自动保存分析结果,文件名称标记分析指标及其分析操作内容和方式。

引子

记得从2016年入学以来,老板叫做的第一个分析便是单因素方差分析,对R 来讲也就是一个aov函数。单因素方差分析确实在很多情况下确实是大多数人的需求,基于R语言的实现也很简单。但是做完之后我便是被吐槽,没有正态性检验,没有方差齐性检验,于是之后的一天,就做了一个简单流程,当时发布在我的个人公众号:微生信生物:《R语言绘制带有显著性字母标记的柱状图》;内容是:首先QQ图,方差齐性检验,后又做了aov和多重比较。并写了个简单的循环。大家可以看到明显不够完整。今天我来的目的就是完善单因素方差分析并且在不适合方差分析的情况下的非参数检验也加入方案中,其次可视化也做了一个完善。(大家有想法的留言,我将加在本框架内,完善差异检测方案)

那么为什么选在这样一个日子来完善这样一份代码呢?刘老师NBT上线不久,其中fig6让我很是怀念之前使用R语言出带显著性标记的柱状图。于是才有了今天可视化方案中的planB。当然刘老师亲自做的NBT的分析值得我逐行运行,尤其是其中的物种分类的GraPhlan,算是到目前为止,我见过的最为漂亮的物种分类树。相信刘老师在之后的文章解读部分会为大家详细解读。

我结合之前的工作同时结合刘老师的工作,设计了一个这样的思路:试图方便自己的工作,同时希望帮到大家。

单因素差异检测完整方案实现思路

试验中经常测定的指标共同的特征是:不同的处理,有重复,需要做差异检测。因此这里我首先对数据进行正态性检验和方差齐性检验,判断符合后进行多重比较并选择喜欢的可视化方案(这里我提供了两种可视化方法,分别是:柱状图,箱线图),判断不符合后,进行非参数检验,首先进行kruskal.test检验(同时对多组进行差异检验),如果有差异,我将继续进行Wilcoxon秩和检验,之后便选择两种可视化方案中的一种进行可视化。

我写了一个骨架:

主要函数解读

两种差异表示方案及其代码

下面是进行LSD多重比较及其添加表征差异的字母柱状图代码:

#值得注意的是LSD多重比较输出的就是字母形式的结果,如果我们选择其他多重比较方法,注意提取差异显著字母格式的结果
out <- LSD.test(model,"group", p.adj="none")#进行多重比较
aa = out$group#结果显示:标记字母法
aa$group = row.names(aa)
wen1 = as.data.frame(tapply(dd,data_wt$group,mean,na.rm=TRUE))
wen2 = as.data.frame(tapply(dd,data_wt$group,sd,na.rm=TRUE))
went = cbind(wen1,wen2)
wentao = merge(aa,went, by="row.names",all=F)
colnames(wentao) = c(colnames(wentao[1:4]),"mean" ,"SD")
aa = mutate(wentao, ymin = mean - SD, ymax =  mean + SD)
a = max(aa$mean)*1.2
# ss <- round(wtx3$`Pr(>F)`[1],3)
p = ggplot(aa , aes(x = group, y = mean,colour= group)) +
    geom_bar(aes(colour= group,fill = group),stat = "identity", width = 0.4,position = "dodge") +
    geom_text(aes(label = groups,y=ymax, x = group,vjust = -0.3,size = 6))+
    geom_errorbar(aes(ymin=ymin,
                    ymax=ymax),
                    colour="black",width=0.1,size = 1)+
    scale_y_continuous(expand = c(0,0),limits = c(0,a))+
    labs(x=paste(name_i,"of all group", sep = "_"),
     y="group",
    title = paste("Normality test",p1,"Homogeneity of variance",p2,"kruskal.test",sumkrusk[3],sep = ":"))
p

字母标记箱线图代码

这部分代码来源自刘老师NBT中的fig6,也正是这份代码,让我产生了重新升级差异分析的念头。说句题外话,刘老师的代码书写思路我总有一种亲切的感觉。

out <- LSD.test(model,"group", p.adj="none")
aa = out$group
aa$group = row.names(aa)
a = max(aa$dd)*1.2

data_box = data_wt[c(1,2,i)]
    colnames(data_box) = c("ID" , "group","dd" )
stat = out$groups
data_box$stat=stat[as.character(data_box$group),]$groups
max=max(data_box[,c("dd")])
min=min(data_box[,c("dd")])
x = data_box[,c("group","dd")]
y = x %>% group_by(group) %>% summarise_(Max=paste('max(',"dd",')',sep=""))
y=as.data.frame(y)
rownames(y)=y$group
data_box$y=y[as.character(data_box$group),]$Max + (max-min)*0.05
p = ggplot(data_box, aes(x=group, y=data_box[["dd"]], color=group)) +
    geom_boxplot(alpha=1, outlier.size=0, size=0.7, width=0.5, fill="transparent") +
    labs(x=paste(name_i," group", sep = "_"),
    y="group",
    title = paste("Normality test",p1,"Homogeneity of variance",p2,sep = ":"))+
    geom_text(data=data_box, aes(x=group, y=y, color=group, label= stat)) +
    geom_jitter( position=position_jitter(0.17), size=1, alpha=0.7)+theme(legend.position="none")
p
if (length(unique(data_box$group))>3){    p=p+theme(axis.text.x=element_text(angle=45,vjust=1, hjust=1))}
        FileName <- paste(name_i,"_aov_LSD_box", ".pdf", sep = "")
        ggsave(FileName, p, width = 8, height = 8)

ggpubr + 箱线图 + 连线差异标注

由于两组之间的连线需要指定两组信息,这里我又想将所有组之间的差异展示出来,所以使用combn函数得到分组信息两两匹配的结果,并使用tapply结合函数将矩阵改变为列表。完成这一工作。

wtq = levels(data_wt$group)
lis = combn(levels(data_wt$group), 2)
x <-lis
my_comparisons <- tapply(x,rep(1:ncol(x),each=nrow(x)),function(i)i)
p = ggplot(data_box, aes(x=group, y=data_box[["dd"]], color=group)) +
          geom_boxplot(alpha=1, outlier.size=0, size=0.7, width=0.5, fill="transparent") +
          labs(x=paste(name_i,"of all group", sep = "_"),
               y="group",
               title = paste("Normality test",p1,"Homogeneity of variance",p2,sep = ":"))+
          # geom_text(data=data_box, aes(x=group, y=y, color=group, label= stat)) +
          geom_jitter( position=position_jitter(0.17), size=1, alpha=0.7)+theme(legend.position="none")+
          stat_compare_means()+
          stat_compare_means(comparisons=my_comparisons,label = "p.signif",hide.ns = F) # Add pairwise

p

实战

测序数据和Rmd代码,生信宝典后台回复“anova”获取。

写在后面

值得注意的是,我们在方差分析之后可以选择t检验两两比对差异,并使用P值矫正。同样非参数检验也使用类似的方法进行两两比对,但是就两组之间的显著性结果转化为字母标记,我并没有实现,大家如果有想法,多多赐教。

本文分享自微信公众号 - 生信宝典(Bio_data)

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-11-04

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • PCA主成分分析实战和可视化 | 附R代码和测试数据

    一文看懂PCA主成分分析中介绍了PCA分析的原理和分析的意义(基本简介如下,更多见博客),今天就用数据来实际操练一下。

    生信宝典
  • 单细胞分析Seurat使用相关的10个问题答疑精选!

    作为一个刚刚开始进行单细胞转录组分析的菜鸟,R语言底子没有,有时候除了会copy外,如果你让我写个for循环,我只能cross my fingers。。。。

    生信宝典
  • NGS基础 - 参考基因组和基因注释文件

    参考基因组和基因注释文件获取 通常测序生成的reads要与参考基因组或参考转录组进行比对,或Pseudo-alignment。所以首先需要获取参考基因组和参考转...

    生信宝典
  • SAP CRM WebUI上Opportunity reason字段的后台配置

    有客户抱怨SAP CRM Web UI上reason这个类型为drop down list的字段在SAP CRM Fiori的My Opportunity应用里...

    Jerry Wang
  • answer

    2.跟我预想中不一样,我这边特地用swift和oc分别敲了一遍,oc的时候array2中只有obj2并且name为test,而swift中,array2依然是o...

    清墨
  • Focus related

    Jerry Wang
  • 基于helium自动化测试的方法进行代码仓库梳理和备份

    helium是一款基于 Selenium 实现的网页自动化工具,他的 API 比 Selenium 更简介,当然也因为他是基于 Selenium 之上构建的,因...

    机械视角
  • BZOJ4819: [Sdoi2017]新生舞会(01分数规划)

    Description 学校组织了一次新生舞会,Cathy作为经验丰富的老学姐,负责为同学们安排舞伴。有n个男生和n个女生参加舞会 买一个男生和一个女生一起跳舞...

    attack
  • YOLO V3网络结构解析

    论文翻译:https://zhuanlan.zhihu.com/p/34945787

    算法发
  • MOS文章实验:ORA-01722 from Queries with Dependent Predicates

    今天读了一篇MOS文章,《ORA-01722, ORA-01839, ORA-01841, ORA-01847 or ORA-01858 from Queri...

    bisal

扫码关注云+社区

领取腾讯云代金券