专栏首页优雅R「Workshop」第四十期 常用的差异分析方法

「Workshop」第四十期 常用的差异分析方法

几种常用的差异分析方法简介

如今在生物学研究中,差异分析越来越普遍,也有许多做差异分析的方法可供选择。但是在实际应用中,大多数人不知道该使用哪种方法来处理自己的数据,所以今天我就来介绍下目前几种常用的差异分析方法及其适用场景。

1.方差分析、T检验、卡方检验、秩和检验


preview

  • 其实核心的区别在于:数据类型不一样。如果是定类和定类,此时应该使用卡方分析;如果是定类和定量,此时应该使用方差或者T检验。
  • 方差和T检验的区别在于,对于T检验的X来讲,其只能为2个类别比如男和女。如果X为3个类别比如本科以下,本科,本科以上;此时只能使用方差分析。

进一步细分

preview

T检验


t检验(student t检验)是应用t分布的特征,将t作为检验的统计量来进行统计推断方法。它对样本要求较小(例如n<30)。

主要用途:

  • 样本均数与总体均数的差异比较
  • 两样本均数的差异比较

单样本t检验

单样本t检验主要用于判断样本均数与总体均数是否存在显著差异。适用条件

已知一个总体均数 已知一个样本均数及该样本标准差 样本正态分布或近似正态总体

实际应用中,当数据量足够大时,对样本正态分布要求不再严格。只要数据分布不是严重偏态,一般来说单样本t检验都是适用的。

R语言中可以用t.test函数进行t检验

从某小学六年级抽取10名学生,其身高(单位:cm),是否认为该学校六年级平均身高130cm?

#生成数据
x <- c(130,120,130,110,130,135,129,124,130,134)
#t检验
t.test(x,mu = 130)

 One Sample t-test

data:  x
t = -1.1884, df = 9, p-value =
0.2651
alternative hypothesis: true mean is not equal to 130
95 percent confidence interval:
 121.8702 132.5298
sample estimates:
mean of x 
    127.2 
#结果显示,P=0.2651>0.05。在统计学上说明样本均数与总体均数没有差别。

独立样本t检验

独立样本t检验主要检验两个样本均数及其所代表的总体之间差异是否显著。

适用条件

  • 独立性,各观察值之间相关独立
  • 正态性,各样本均来自正态分布的总体
  • 方差齐性,各样本所在总体的方差相等

方差齐性可以用car包leveneTest函数检验

leveneTest(y=,group =)

其中,y是两组样本组成的数据,group是两组样本的分组情况。

方差齐性检验之后,才可进行独立样本t检验。

用t.test(A,B,var.equal=TRUE,paired=FALSE)

A、B为数据集,var.equal=TRUE为方差齐性。paired=FALSE非配对样本。

示例:

(虚构)有两组学生(每组10人),一组采用传统教育,一组采用素质教育。一学期后,两组学生语文成绩(满分100)如下。问两组学生成绩之间差别是否显著。

A <- c(85,84,95,73,77,65,85,93,90,91)
B <- c(87,96,77,80,79,96,93,82,84,86)
#方差齐性检验
#合并数据
y <- c(A,B)
#数据分组标签
group=as.factor(c(rep(1,10),rep(2,10)))
#载入car包
library(car)
#方差齐性检验
leveneTest(y=y,group = group)
#结果显示,P=0.5505>0.05。说明方差齐性。
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  1  0.3703 0.5505
      18  
#独立样本t检验
t.test(A,B,paired = FALSE)
#结果显示P=0.5639。说明两者没有区别。
 Welch Two Sample t-test

data:  A and B
t = -0.589, df = 16.463, p-value = 0.5639
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -10.100024   5.700024
sample estimates:
mean of x mean of y 
     83.8      86.0 

配对样本t检验

配对样本t检验同样检验两个样本均数及其所代表的总体之间差异是否显著。

独立样本t检验与配对样本t检验同属于双样本t检验,不同点在于配对样本t检验要求两个样本之间存在某些配对关系。

常见配对关系

同一样本两种不同处理方法的检验结果 同一样本前后时间点的检验结果

适用条件

正态性

示例

有20名女性分为10对,试吃两种药。经过一段时间后,药效如下。问两种药是否有区别

#生成数据
drug1 <- c(4.4,5,5.8,4.6,4.9,4.8,6,5.9,4.3,5.1)
drug2 <- c(6.2,5.2,5.5,5,4.4,5.4,5,6.4,5.8,6.2)
#配对样本t检验
t.test(drug1,drug2,paired = TRUE)
#结果显示,P=0.1575>0.05,不能说两者存在显著差别。
 Paired t-test

data:  drug1 and drug2
t = -1.5417, df = 9, p-value = 0.1575
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.0609306  0.2009306
sample estimates:
mean of the differences 
                  -0.43 

方差检验


方差分析(analysis of variance ,ANOVA)就是通过检验各总体的均值是否相等来判断分类型自变量对数值型因变量是否有显著影响。

示例

我们使用的是R里内置的“npk”数据集,该数据集由24行和5列数据组成,第一列代表区组(共6个),N、P和K分别代表氮、磷和钾元素的使用情况,yield代表豌豆产量,该数据集主要是用来研究不同肥料对豌豆产量的影响。

fit <- aov(yield ~ N, data=npk)
fit <- aov(yield ~ N + block, data=npk)

卡方检验


卡方检验就是统计样本的实际观测值与理论推断值之间的偏离程度,实际观测值与理论推断值之间的偏离程度就决定卡方值的大小,卡方值越大,越不符合;卡方值越小,偏差越小,越趋于符合,若两个值完全相等时,卡方值就为0,表明理论值完全符合。

适用条件

1.所有的理论数T≥5并且总样本量n≥40,用Pearson卡方进行检验. 2.如果理论数T<5但T≥1,并且n≥40,用连续性校正的卡方进行检验. 3.如果有理论数T<1或n<40,则用Fisher’s检验.

示例

判断5种品牌啤酒的爱好者有无显著差异:

img

x<-c(210,312,170,85,223)
chisq.test(x)

img

x<-matrix(c(46,18,6,8),ncol=2,nrow=2)
chisq.test(x)
chisq.test(x)$expected ###查看理论值
fisher.test(x)  ##进行fisher检验

秩和检验


秩和检验是对原假设的非参数检验,在不需要假设两个样本空间都为正态分布的情况下,测试它们的分布是否完全相同。

library(stats)
data("mtcars")
wilcox.test(mpg~am,data = mtcars)

几种常用的R包

目前常用差异分析的R包有edgeR、limma、DESeq2

img

img

三种包的区别:

          1.limma包做差异分析要求数据满足正态分布或近似正态分布,如基因芯片、TPM格式的高通量测序数据。
          2.通常认为Count数据不符合正态分布而服从泊松分布。对于count数据来说,用limma包做差异分析,误差较大
          3.DESeq2、和 EdgeR都是基于count,然后两个都是NB(negative binomial)但是在估计dispersion parameter的方法上面不一样。
          4.limma,edgeR,DESeq2三大包基本是做转录组差异分析的金标准,大多数转录组的文章都是用这三个R包进行差异分析。
          5.edgeR差异分析速度快,得到的基因数目比较多,假阳性高(实际不差异,结果差异)。DESeq2差异分析速度慢,得到的基因数目比较少,假阴性高(实际差异,结果不差异)。
          6.需要注意的是制作分组信息的因子向量是,因子水平的前后顺序,在R的很多模型中,默认将因子向量的第一个水平看作对照组

如果数据量大并且要求比较conservative的话可以所有方法都用下,然后取并集。

数据预处理


library("edgeR")
expr = read.csv("mRNA_exprSet.csv",sep = ',',header=T)  
head(expr)
expr = avereps(expr[,-1],ID = expr$X) # 对重复基因名取平均表达量,然后将基因名作为行名
expr = expr[rowMeans(expr)>1,] # 去除低表达的基因
library(stringr)
tumor <- colnames(expr)[as.integer(substr(colnames(expr),14,15)) < 10]
normal <- colnames(expr)[as.integer(substr(colnames(expr),14,15)) >= 10]

tumor_sample <- expr[,tumor]
normal_sample <- expr[,normal]

exprSet_by_group <- cbind(tumor_sample,normal_sample)
group_list <- c(rep('tumor',ncol(tumor_sample)),rep('normal',ncol(normal_sample)))

edgeR


对于edgeR的分析流程而言,我们需要输入的数据包括:

  1. 表达矩阵(counts
  2. 分组信息(group
  3. 拟合信息(design):指明如何根据样本的分组进行建模

edgeR默认使用 trimmed mean of M-values (TMM) 计算文库的scale factor进行normalization,以最大程度地缩小样本间基因表达量的log-fold change。这是因为TMM 法认为样本间大部分的基因都没有发生差异表达,而那些真正差异表达的基因并不会受到normalization的严重影响。如此一来,便将那些由于测序引起的差异表达基因的表达量给校正了,消除了一部分的假阳性。

data = exprSet_by_group
group_list = factor(group_list)
design <- model.matrix(~0+group_list)
rownames(design) = colnames(data)
colnames(design) <- levels(group_list)
DGElist <- DGEList( counts = data, group = group_list)
## Counts per Million or Reads per Kilobase per Million
keep_gene <- rowSums( cpm(DGElist) > 1 ) >= 2 ## 自定义
table(keep_gene)
DGElist <- DGElist[ keep_gene, , keep.lib.sizes = FALSE ]

DGElist <- calcNormFactors( DGElist )
DGElist <- estimateGLMCommonDisp(DGElist, design)##为所有基因都计算同样的dispersion
DGElist <- estimateGLMTrendedDisp(DGElist, design)##根据不同基因的均值--方差之间的关系来计算dispersion,并拟合一个trended model
DGElist <- estimateGLMTagwiseDisp(DGElist, design)##计算每个基因的dispersion,并利用经验性贝叶斯方法(empirical bayes)压缩至Trend Didspersion。个人认为这一项相当于GLM中每个基因的beta值
#####################负二项式广义对数线性模型
fit <- glmFit(DGElist, design)
results <- glmLRT(fit, contrast = c(-1, 1)) 
####################类似然负二项式广义对数线性模型
fit <- glmQLFit(dge, design, robust = TRUE)        #拟合模型
lrt <- glmQLFTest(fit) 

nrDEG_edgeR <- topTags(results, n = nrow(DGElist))
nrDEG_edgeR <- as.data.frame(nrDEG_edgeR)
head(nrDEG_edgeR)

padj = 0.01 # 自定义
foldChange= 2 # 自定义
nrDEG_edgeR_signif  = nrDEG_edgeR[(nrDEG_edgeR$FDR < padj & 
                       (nrDEG_edgeR$logFC>foldChange | nrDEG_edgeR$logFC<(-foldChange))),]
nrDEG_edgeR_signif = nrDEG_edgeR_signif[order(nrDEG_edgeR_signif$logFC),]

如果没有生物学重复

img

library(edgeR)
##跟DESeq2一样,导入数据,预处理(用了cpm函数)
exprSet<- read.table(file = "test.txt", sep = "\t", header = TRUE, row.names = 1, stringsAsFactors = FALSE)
group_list <- factor(c(rep("Contral",1),rep("Treat",1)))

##设置分组信息,并做TMM标准化
exprSet <- DGEList(counts = exprSet, group = group_list)
bcv = 0.4  #设置BCV值
et <- exactTest(exprSet, dispersion=bcv^2)
write.csv(topTags(et, n = nrow(exprSet$counts)), 'result.csv', quote = FALSE)   #输出主要结果

limma,DEseq2的代码实现可以在https://blog.csdn.net/weixin_43700050/article/details/98085127找到。

在这里插入图片描述

参考:

1.https://www.jianshu.com/p/bb0bd72bc428.

2.https://zhuanlan.zhihu.com/p/57756620.

3.https://blog.csdn.net/weixin_43700050/article/details/98085127.

4.https://blog.csdn.net/u010608296/article/details/114135253?utm_medium=distribute.pc_relevant.none-task-blog-2%7Edefault%7EBlogCommendFromMachineLearnPai2%7Edefault-2.control&dist_request_id=1328761.11124.16172439205851981&depth_1-utm_source=distribute.pc_relevant.none-task-blog-2%7Edefault%7EBlogCommendFromMachineLearnPai2%7Edefault-2.control.edgeR/limma/DESeq2

5.https://www.jianshu.com/p/17ae5bee7bb0

6.https://www.jianshu.com/p/517167c50a5f

本文分享自微信公众号 - 优雅R(elegant-r),作者:宁伟

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

原始发表时间:2021-04-20

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

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 参会见闻系列:ACL 2018,在更具挑战的环境下理解数据表征及方法评价

    AI 科技评论按:本篇属于「顶会见闻系列」。每年这么多精彩的人工智能/机器学习会议,没去现场的自然可惜,在现场的也容易看花眼。那么事后看看别的研究员的见闻总结,...

    AI科技评论
  • 看过 CVPR 2018 workshop 后,发现有一个我不认识的 Lady Gaga

    当地时间 6.18-22 日,CVPR 2018 将在美国盐湖城举办。目前距离大会开幕还有四个月,随着春节期间大会主办方公布接收论文名单,引来大家对 CVPR ...

    AI研习社
  • 动态 | KDD首推Health Day,探讨如何用AI、数据挖掘改变未来医疗 | KDD 2018

    伦敦当地时间8月20日, ACM SIGKDD (知识发现与数据挖掘会议)正式在伦敦开幕,AI 科技评论也来到了现场为大家报道现场的盛况。

    AI科技评论
  • 最常用的四大数据分析方法

    描述型分析:发生了什么?这是最常见的分析方法。在业务中,这种方法向数据分析师提供了重要指标和业务的衡量方法。

    加米谷大数据
  • 学界 | 跟着大神回顾ACL 2018:大会亮点一览

    很高兴看到很多论文都在从方法上研究现有模型以及它们捕获的内容,而不是一直在引入更新的模型。进行这样的研究最常用的办法是自动创建一个侧重于泛化行为的某个特定方面的...

    机器之心
  • 图说 | 刚刚结束的AI盛会NIPS 2017上,你需要知道的所有细节

    机器之心
  • 近九千人齐聚 NeurIPS 2018,四篇最佳论文,十大研究领域,1010 篇论文被接收

    AI 科技评论消息:NeurIPS 2018 于 12 月 3 日—12 月 8 日在加拿大蒙特利尔会展中心(Palais des Congrès de Mon...

    AI科技评论
  • Kaggle 图像识别新赛来袭,还有一家中国企业提供赞助

    AI 研习社此前介绍过 CVPR 2018 workshop 上的多个比赛,详情参见看过 CVPR 2018 workshop 后,发现有一个我不认识的 Lad...

    AI研习社
  • CVPR 2019细粒度图像分类竞赛中国团队DeepBlueAI获冠军 | 技术干货分享

    近日,在Kaggle上举办的CVPR 2019 Cassava Disease Classification挑战赛公布了最终结果,国内团队 DeepBlueAI...

    磐创AI
  • 「Workshop」第三期:生存分析

    生存函数:个体存活到某个时间点t的概率,或者说到时间t为止,感兴趣的事件(T)没有发生的概率:

    王诗翔呀
  • 不忘初心,砥砺前行|VALSE之行收获

    随着VALSE2019的精彩落幕,SIGAI迎来了建号一周年庆。回顾这365天的时间,我们用130篇优质原创技术文章,与人工智能领域数万的专业人士、学生以及爱好...

    SIGAI学习与实践平台
  • CVPR 2019细粒度图像分类竞赛中国团队DeepBlueAI获冠军 | 技术干货分享

    近日,在Kaggle上举办的CVPR 2019 Cassava Disease Classification挑战赛公布了最终结果,国内团队 DeepBlueAI...

    新智元
  • 吃饺子不如撸代码!今年冬至 workshop 干货都在这了

    当然有啦!12月22日,中国专业IT社区CSDN在中关村创业大街为广大开发者带来了一场精彩绝伦的区块链实战开发Workshop,本次活动作为CSDN Block...

    区块链大本营
  • 犀牛鸟视野|现场报道 SIGMOD 2019 数据库顶级会议 (上篇)

    ACM SIGMOD/PODS 2019 数据管理国际会议 6 月30日到7 月5日在荷兰首都阿姆斯特丹召开。会场位于阿姆斯特丹市的一座建筑Beurs van...

    腾讯高校合作
  • 犀牛鸟视野|现场报道 SIGMOD 2019数据库领域顶级会议 (上篇)

    ACM SIGMOD/PODS 2019 数据管理国际会议6月30日到7月5日在荷兰首都阿姆斯特丹召开。会场位于阿姆斯特丹市的一座建筑Beurs van Be...

    腾讯高校合作
  • 「Workshop」第一期:我理解的(生信)数据分析核心基础

    我在简书和公众号上已经分享了很多之前学习的数据分析笔记和文章,覆盖了各方面的内容,数据分析方面以后不会再个人分享特别基础的东西了。接下来我会让师弟师妹们定期分享...

    王诗翔呀
  • Kaggler 看过来,CVPR 2018 Workshop 植物识别赛来袭

    朱晓霞
  • 动物森友会首届「AI 顶会」ACAI 2020开幕,4大主题、17场演讲,戴着口罩做Presentation

    筹备了三个月的动物森友会首届「AI 顶会」ACAI 2020,已于北京时间 7 月 24 日凌晨正式开幕。

    机器之心
  • 推荐 | 10个好用的Web日志安全分析工具

    首先,我们应该清楚,日志文件不但可以帮助我们溯源,找到入侵者攻击路径,而且在平常的运维中,日志也可以反应出很多的安全攻击行为。

    辞令

扫码关注云+社区

领取腾讯云代金券