首先,看一张图,对于组间差异分析有一个整体的了解:
那么问题来了,什么是组间差异检验?就是组间的差异分析以及显著性检验,应用统计学上的假设检验方法,检验组间是否有差异及其差异程度。
坦率地讲,所有的差异检验都基于一个假设:组间没有差异,变量之间没有关系
(即原假设,
H
0
H_0
H0)。也说方差分析其实研究的就是不同水平下是否有差异化的假设检验问题。而假设检验就是先对总体参数提出某种假设,然后利用样本信息判断假设是否成立的过程。
首先回顾关于假设检验的基本概念:
总体
(population):包含所研究的全部个体(数据)的集合。样本
(sample):从总体中抽取的一部分元素的集合。参数
(parameter):用来描述总体特征的概括性数字度量。统计量
(statistic):用来描述样本的概括性数字度量。对于参数来说,因为这里介绍的是组间差异检验,因此在这个水平上分为:参数检验
和非参数检验
。
那么==什么叫参数检验和非参数检验,它们之间的区别是什么呢?==要理解前面的问题,首先需要明白统计推断
的概念。
统计推断是研究如何利用样本数据来推断总体特征的统计学方法,包括参数估计
和假设检验
两大类。总体的参数一般是未知的,通常可以用样本统计量来对总体的参数进行估计,例如可以用样本均值对总体均值进行点估计,利用样本均值的分布对总体均值进行区间估计,这些都称为参数估计。
对未知参数的假设进行检验称为参数统计,所用的检验叫做参数检验
(Parameter test)。不依赖总体分布的具体形式,也不对参数进行估计或检验的统计方法,叫做非参数统计
,其检验方法就是非参数检验
(Non-parametric test)
参数检验和非参数检验的区别:
那么==什么时候用参数检验,什么时候用非参数检验呢?==非参数检验一般不直接用样本观察值作分析,统计量的计算基于原始数据在整个样本中的秩次,丢弃了观察值的具体数值,因此凡适合参数检验的资料,应首选参数检验。但是不清楚是否合适参数检验的资料,则应采用非参数检验。
Tips:这里,我们重新复习一下假设检验的4 个步骤(提出假设;构造检验统计量;根据显著水平,确定临界值和拒绝域;做出检验决策)
了解研究对象整体处于什么状态,是一件非常重要的事情。三大抽样分布( t − t- t−分布、 χ 2 \chi^2 χ2分布、 F − F- F−分布)和正态分布共同构成了现代数理统计学的基础,其中,正态分布和 t − t- t−分布是关于均值的分布; χ 2 \chi^2 χ2分布、 F − F- F−分布是关于方差的分布。
离开分布,假设检验无从谈起;离开假设检验,差异分析毫无根基。同样地,这里复习一下这几个抽样分布:
χ
2
\chi^2
χ2分布
设
X
1
,
X
2
,
.
.
.
,
X
n
X_1,X_2,...,X_n
X1,X2,...,Xn相互独立,都服从标准正态分布
N
(
0
,
1
)
N(0, 1)
N(0,1),则称随机变量
χ
2
=
X
1
2
+
X
2
2
+
.
.
.
+
X
n
2
\chi^2=X_1^2+X_2^2+...+X_n^2
χ2=X12+X22+...+Xn2所服从的分布为自由度为n的
χ
2
\chi^2
χ2分布
设随机变量 X 是自由度为 n 的
χ
2
\chi^2
χ2随机变量, 则其概率密度函数为:
g
n
(
x
)
=
{
1
2
n
2
Γ
(
n
2
)
x
n
2
−
1
e
−
n
2
x>0
0
x <= 0
g_n(x)=\begin{cases} \frac{1}{2^{\frac{n}{2}}\Gamma(\frac{n}{2})}x^{\frac{n}{2}-1}e^{-\frac{n}{2}} & \text{ x>0 }\\ 0 & \text{ x <= 0 } \end{cases}
gn(x)={22nΓ(2n)1x2n−1e−2n0 x>0 x <= 0
Γ
(
⋅
)
\Gamma(\cdot )
Γ(⋅)表示的是一个gamma函数,它是整数k的封闭形式。
χ
n
2
\chi _{n}^{2}
χn2的密度函数
g
n
(
x
)
g_{n}(x)
gn(x)形状如下图:
t
t
t分布
设
X
1
X_1
X1服从标准正态分布N(0,1),
X
2
X_2
X2服从自由度为n的
χ
2
\chi^2
χ2分布,且
X
1
X_1
X1、
X
2
X_2
X2相互独立,则称变量
t
=
X
1
/
(
X
2
/
n
)
1
/
2
t=X_1/(X_2/n)^{1/2}
t=X1/(X2/n)1/2所服从的分布为自由度为n的
t
−
分
布
t-分布
t−分布。
设随机变量
T
∼
t
n
T ∼ t_{n}
T∼tn, 则其密度函数为
t
n
(
x
)
=
Γ
(
n
+
1
2
)
Γ
(
n
2
)
n
π
(
1
+
x
2
n
)
−
n
+
1
2
,
−
∞
<
x
<
∞
t_n(x) = \frac{\Gamma (\frac{n+1}{2} )}{\Gamma(\frac{n}{2})\sqrt{n\pi }}(1+\frac{x^2}{n})^{-\frac{n+1}{2}}, -\infty <x< \infty
tn(x)=Γ(2n)nπ
Γ(2n+1)(1+nx2)−2n+1,−∞<x<∞
该密度函数的图形如下:
F分布
设
X
1
X_1
X1服从自由度为m的
χ
2
\chi^2
χ2分布,
X
2
X_2
X2服从自由度为n的
χ
2
\chi^2
χ2分布,且
X
1
X_1
X1、
X
2
X_2
X2相互独立,则称变量
F
=
(
X
1
/
m
)
/
(
X
2
/
n
)
F=(X_1/m)/(X_2/n)
F=(X1/m)/(X2/n)所服从的分布为F分布,其中第一自由度为m,第二自由度为n。一般地,这里F就是均方之比。
若随机变量
Z
∼
F
m
,
n
Z ∼F_{m,n}
Z∼Fm,n, 则其密度函数为
f
m
,
n
(
x
)
=
{
Γ
(
m
+
n
2
)
Γ
(
n
2
)
Γ
(
m
2
)
m
m
2
n
n
2
x
m
2
−
1
(
n
+
m
x
)
−
m
+
n
2
x>0
0
others
f_{m,n}(x) = \begin{cases} \frac{\Gamma(\frac{m+n}{2})}{\Gamma(\frac{n}{2})\Gamma(\frac{m}{2})}m^{\frac{m}{2}}n^{\frac{n}{2}}x^{\frac{m}{2}-1}(n+mx)^{-{\frac{m+n}{2}}} & \text{ x>0 }\\ 0 & \text{ others } \end{cases}
fm,n(x)={Γ(2n)Γ(2m)Γ(2m+n)m2mn2nx2m−1(n+mx)−2m+n0 x>0 others
自由度为 m, n 的 F 分布的密度函数如下图:
注意 F 分布的自由度 m 和 n 是有顺序的, 当
m
≠
n
m\neq n
m=n时, 若将自由度 m 和 n 的顺序颠倒一下, 得到的是两个不同的 F 分布。
补充知识点:
数据在使用前要注意采用有效的方法收集数据,如设计好抽样方案,安排好试验等等。只有有效的收集了数据,才能有效地使用数据,开展统计推断工作。获得数据后,根据问题的特点和抽样方式确定抽样分布,即统计模型.。基于统计模型,统计推断问题可以按照如下的步骤进行:
中心极限定理
或其它极限定理
找出统计量的极限分布
。其中第二步是最重要,但也是最困难的一步。统计三大分布及正态总体下样本均值和样本方差的分布,在寻求与正态变量有关的统计量精确分布时,起着十分重要作用。尤其在求区间估计和假设检验问题时可以看得十分清楚。
不管是参数检验还是非参数检验,都要基于特定的分布来做假设检验。当总体分布已知时,例如总体服从正态分布,我们可以根据给定的显著性水平(通常为0.01 或0.05)查表获得临界值。当总体分布未知时,可以先用Permutation test
构造经验分布,再根据显著性水平获得临界值。
传统的统计量检验的方法是在检验之前确定显著性水平α,也就意味着事先确定了临界值和拒绝域。这样,不论检验统计量的值是大还是小,只要它的值落入拒绝域就拒绝原假设,否则就不拒绝原假设。这种给定显著性水平的方法,无法给出观测数据与原假设之间不一致程度的精确度量。
要测量出样本观测数据与原假设中假设值的偏离程度,则需要计算pvalue值。pvalue 值,也称为观测到的显著性水平,它表示为如果原假设 H 0 H_0 H0正确时得到实际观测样本结果的概率。pvalue 值越小,说明实际观测到的数据与 H 0 H_0 H0之间的不一致的程度就越大,检验的结果就越显著。
变量较多,判断组间差异时需要多重检验的情况在宏基因组扩增子差异分析中十分常见。这种情况下,基于单次比较的检验标准将变得过于宽松,使得阳性结果中的错误率(FDR
值FalseDiscovery Rate)非常大。怎么办呢?最好的办法就提高判断的标准(p value),单次判断的犯错概率就会下降,总体犯错的概率也将下降。
在多重检验中提高判断标准的方法,我们就称之为多重检验校正。从1979 年以来,统计学家提出了多种多重检验校正的方法。相应地,对p值校正之后的叫法也不一样,比如,FDR
、Q value
、Adjusted p-value
,因此知道在多重检验时需要校正就行了。
这里我们讨论的是统计推断。换句话说,就是找差异。
在数据科学家的工具箱里,这是一款经久不衰、常用常新的瑞士军刀。几乎只要想到差异分析,就会想到箱线图。也开发出类箱线图的工具比如小提琴图(小提琴图Violin plot),例如:
可用的R包有:geom_boxplot() {ggplot2}
散点图也是一款百搭的工具,可以和箱线图结合着用,当然多元分析大多也得借助这个的散点图。比如,回归分析、排序(PCA。CA、CCA、RDA,NMDS,PCoA)、聚类(均值聚类 、划分)用散点图来反映都是比较直观地。也开发有新的散点图比如叫火山图
。示例如下:
可用的R包有:geom_point(){ggplot2}
热图可以简单地聚合大量数据,并使用一种渐进的色带来优雅地表现出来,可以很直观地展现数据的相对大小。在生物医学研究中,常用来展现基因表达或丰度数据,当然用它表达相关系数大小也是允许的。当然也有开发的热图,比如地理热图等。例如:
可用的R包有:heatmap;pheatmap
一般有进化树
和层次聚类树
,如果想表达对象之间的距离差异,最直观的的也许就是树状图了。为了用图表示亲缘关系,把分类单位摆在图上树枝顶部,根据分枝可以表示其相互关系,具有二次元和三次元。在数量分类学上用于表型分类的树状图,称为表型树状图
(phenogram),掺入系统的推论的称为系统树状图
(cladogram)以资区别。
可用的R包有:ggtree;cluster
那么了解了有哪些图可以用来展示差异,接下来就是如何寻找差异。
这里说的基于类型标签言下之意是通过统计分析,可以有针对性的找出分组间丰度变化差异显著的类别(在微生物组学分析中,即微生物的物种),并得到差异类型在不同分组间的富集情况,同时,可以比较组内差异和组间差异的大小,判断不同分组间的群落结构差异是否具有显著意义。也就是说,可以找出区别组间的一个biomarker。
这类检验一般只输出p值,它的目的很简单,就是检验比较组之间的相似性距离是否有差异。常用的分析方法有卡方检验
、Student t检验
、Wilcoxon秩和检验
等等。
如果只有两个样本比较,适合用卡方检验,不过检验出来的结果没什么可靠性,因为现阶段16s研究不做重复实在“难以服众”了。价格便宜,做重复压根没有难度,而且从生物学、统计学角度考虑,也需要做重复。
如果是两组样本(至少3重复),可以试一下Student t,Welch‘st以及Wilcoxon秩和检验。Student t检验需要样本符合正态分布,而且方差对齐。当组间样本数不同,方差也不对齐的时候,Welch’s t检验是很好的选择。
Student’s t-test: this test assumes that both groups of data are sampled from populations that follow a normal distribution and that both populations have the same variance. Welch’s t-test: this test assumes that both groups of data are sampled from populations that follow a normal distribution, but it does not assume that those two populations have the same variance. reference:https://www.statology.org/welchs-t-test/
Wilcoxon秩和检验又叫Mann-Whitney U 检验,是基于变量排名的一种统计方法,不需要样本符合正态分布,也不需要样本方差对齐,是更为广泛的检验方法,但同时也由于检验太宽松,容易带来很多假阳性。
如果是多组样本比较,可以选择one way ANOVA、TURKEY以及Kruskal-Wallis H检验等方法。one way ANOVA和TURKEY其实都是基于方差分析,只不过后者带有后验,可以知道两个分组对整体差异的贡献度。
Kruskal-Wallis H检验本质也是一种秩和检验,与前两者的区别在于,它不需要样本数和方差的对齐,应用更为广泛。Kruskal-Wallis检验又被称之为单因素非参数方差分析。
毫不客气地讲,一般秩和检验或置换检验属于非参数检验。在这类差异检验中,有两种集成方法特别值得我们注意:LEfSe 、metastats。
1. 首先在多组样本中采用的非参数检验Kruskal-Wallis秩和检验检测不同分组间丰度差异显著的特征; 2. 然后在上一步中获得的显著差异特征,用成组的Wilcoxon秩和检验进行组间差异分析(若没有亚组,该步跳过); 3. 最后用线性判别分析(LDA)对数据进行分类和评估差异显著的物种的影响力(即LDA score)。 reference:https://blog.csdn.net/weixin_42072765/article/details/108356184
得到结果展示如下,差异体现在柱形图和树状图上。
LDA值分布柱状图中展示了LDA Score大于设定值(默认设置为4)的物种,即组间具有统计学差异的Biomarker。展示了不同组中丰度差异显著的物种,柱状图的长度代表差异物种的影响大小(即为 LDA Score)。
在进化分支图中,由内至外辐射的圆圈代表了由门至属(或种)的分类级别。在不同分类级别上的每一个小圆圈代表该水平下的一个分类,小圆圈直径大小与相对丰度大小呈正比。着色原则:无显著差异的物种统一着色为黄色,差异物种Biomarker跟随组进行着色,红色节点表示在红色组别中起到重要作用的微生物类群,绿色节点表示在绿色组别中起到重要作用的微生物类群,若图中某一组缺失,则表明此组中并无差异显著的物种,故此组缺失。图中英文字母表示的物种名称在右侧图例中进行展示。
biomaker在不同组各样本中的丰度比较图:将biomaker丰度最高的样本的丰度设定为1,其他样品中该 biomarker 的丰度为相对于丰度最高样品的相对值。
+ 将丰度数据归一化成为相对丰度 + 组间T-test计算 + 显著性检验 + Permutation test 置换检验 + 重复数 ≥8 与重复数 <8的p值计算规则不同 + 重复数 ≥8:只开展单物种的置换检验 + 重复数 <8:将混合整个样本进行置换检验 + 组内某个物种的数目少于样本重复数的时候,会利用Fisher精确检验进行p值计算 + 多重检验
Metastats实际上是非参数多重检验和p值校正的整合,而LEfSe则是Metastats和LDA判别的整合。当然,由于Metastats采用的非参数t检验,只能分析两个分组;而LEfSe则因为使用的Kruskal-Wallis秩和检验可以分析两个以上的分组。
所谓基于距离也就是检验的是群落差异而不是某个物种。上面所提及的检验方法,其实都只能告诉我们那些分组是否有显著差异(可以简单理解为有无)。那如果想同时知道这些差异的程度(可以简单理解为多少),需要Anosim,Adonis以及MRPP等检验方法。这些方法不但可以输出检验显著性结果(p值),还有程度结果(R值),R值可以用来判断分组贡献度大小。Anosim、Adonis这些可用于多元统计检验的模型就非常适合了。R语言vegan包含有多种非参数检验方法,包括Anosim、Adonis、MRPP等,不同方法在统计量的选择、零模型等方面存在差异。
值得注意的是,Anosim本质是基于排名的算法,其实与NMDS的配合效果最好。如果是PCoA分析,建议配合使用Adonis检验结果。
#读取抽平后的OTU_table和环境因子信息 data=read.csv("otu_table.csv", header=TRUE, row.names=1) envir=read.table("environment.txt", header=TRUE) rownames(envir)=envir[,1] env=envir[,-1] #筛选高丰度物种并将物种数据标准化 means=apply(data, 1, mean) otu=data[names(means[means>10]),] otu=t(otu) #根据地理距离聚类 kms=kmeans(env, centers=3, nstart=22) Position=factor(kms$cluster) #进行Anosim分析 library(vegan) anosim=anosim(otu, Position, permutations=999) summary(anosim)
具体说来,Anosim分析的原理是先计算样品两两之间的距离,将样品两两之间的距离按照从大到小进行排序并计算排名(秩,r),并根据距离的归类(属于组间距离还是组内距离)来计算组间距离秩的均值
r
b
r_b
rb与组内距离秩的均值
r
w
r_w
rw之差作为统计量:
R
=
r
b
−
r
w
0.25
[
n
(
n
−
1
)
]
R=\frac{r_b-r_w}{0.25[n(n-1)]}
R=0.25[n(n−1)]rb−rw
其中:
R的范围为[-1,1]
R>0 说明组间差异大于组内差异,R<0 组间差异小于组内差异。
R 只是组间是否有差异的数值表示,并不提供显著性说明。
P值 则说明不同组间差异是否显著,该P值通过**置换检验(Permutation Test)**获得。置换检验大致原理:(假设原始分组为实验组和对照组)
Adonis
ADONIS又称置换多因素方差分析(permutational MANOVA)或非参数多因素方差分析(nonparametric MANOVA),是一种基于Bray-Curtis距离的非参数多元方差分析方法。它与Anosim的用途其实差不多,也能够给出不同分组因素对样品差异的解释度(R值)与分组显著性(P值)。不同点是应用的检验模型不同,ADONIS本质是基于F统计量的方差分析,所以很多细节与上述方差分析类似。该方法可分析不同分组因素对样本差异的解释度,并使用置换检验对分组的统计学意义进行显著性分析。ADONIS分析使用R vegan包adonis函数进行分析。
在微生物的分析中我们通常把Adonis和PCA分析结合在一起。进行完PCA分析后,我们想要检验不同的分组之间究竟是否有差异,差异是否显著,这时候我们就可以用Adonis检验。
Mantel test
尽管Mantel test通常用于确定两个距离矩阵的相关性,但也可用于检验假设或模型。通过在模型矩阵中比较组间距离与组内距离的差异程度,用以确定分组是否显著。它的原假设是两个矩阵间没有相关关系。
检验过程如下:两个矩阵都对应展开,变量两列,计算相关系数(理论上什么相关系数都可以计算,但常用pearson相关系数),然后其中一列或两列同时置换,再计算一个值,permutation 成千上万次,看实际的r值在所得r值分布中的位置,如果跟随机置换得到的结果站队较近,则不大相关,如果远远比随机由此得到显著性。
此时Mantel test和ANOSIM的工作方式相似,但其特殊形式在于,为模型矩阵选择的特定值是根据距离数值本身而非根据排位确定的。如下概括了Mantel test确定分组差异的方法:
MRPP
与Anosim类似,但是MRPP(Multi Response Permutation Procedure)是基于Bray-Curtis的参数检验,利用组内和组间差异的置换检验,确定两组或两组以上数据集有无差异的非参数过程(Mielke, 1976)。通常配合PCA、PCoA、NMDS等降维图使用,MRPP分析使用R vegan包mrpp函数。
mrpp算法首先计算整个数据集中的所有成对距离,然后计算各组内对象间的平均距离
d
i
ˉ
\bar{d_i}
diˉ,之后计算δ(组内距离的加权平均)。
δ
=
∑
i
=
1
g
C
i
d
i
ˉ
\delta =\sum_{i=1}^{g}C_i\bar{d_i}
δ=i=1∑gCidiˉ
C
i
=
n
i
N
C_i=\frac{n_i}{N}
Ci=Nni
式中,g为总分组的数量;ni表示第i个分组的对象数量,N为数据集中的总对象数量。
然后它对样本及其相关的成对距离进行置换,并根据置换数据重新计算 δ。它重复置换步骤多次。δ小于初始δ的概率即为p值。p值即代表了检验的显著性信息,p越低表明越容易接受观测δ值,MRPP结果越可信。
MRPP结果中通常会提供两种δ值:observed δ,即直接由公式计算的观测δ,值越小表明组内差异越小;以及expect δ,由置换过程得到的平均δ,值越大暗示了组间差异越大。同时会结合observed δ和expect δ再计算一个简称为A值(chance-corrected within group agreement)的统计量
A
=
1
−
o
b
s
e
r
v
e
d
δ
e
x
p
e
c
t
δ
A=1-\frac{observed \delta}{expect \delta}
A=1−expectδobservedδ,小于0表明组内差异大于组间差异,大于0表明组间差异大于组内差异。
首先,对上面差异检验的方法做一个汇总:
方法 | R值 | p值 |
---|---|---|
Anosim | R-value介于(-1,1)之间,R-value大于0,说明组间差异显著 | P< 0.05 表示统计具有显著性 |
Adonis | R2 表示不同分组对样本差异的解释度 | Pr表示P 值,小于0.05 说明本次检验的可信度高 |
Amova | — | p-value表示P 值,小于0.05 说明组间差异显著 |
Mantel test | r为相关系数,r值越大两矩阵相关性越大 | P<0.05表示统计具有显著性 |
MRPP | A值大于0说明组间差异大于组内差异 | Significance值小于0.05说明差异显著 |
然后,再盗用一张比较全面的图,对上述内容做一个总结: