群落多样性之Alpha多样性（二）

2019-04-13 本文已影响19人布莱特杨

相信不少人看过下面这个类型的故事，我把它稍作了个改编。
超逸和艾斯无意中发现有两个比自己来的晚的同事，叫山农和辛普森，都已升职加薪，而自己却一直原地不动。
终于有一天，超逸忍不住了，冒着被解雇的危险，找老板理论：“老板，我有过迟到、早退或乱章违纪的现象吗？”
“没有！”
那是公司对我有偏见吗？”
“没有！”
“为什么比我资历浅的人，工资却比我高？”
老板说：”咱先不说这个，眼下有个事儿，你去调查下今天市场上在卖哪几种蔬菜？“
A很快去市场转了一圈，很快回来说：“报告老板，今天市场上的蔬菜主要有白菜、萝卜、番茄、土豆……！”
“价格分别是多少？”
“这个我没问！”
“都是哪些公司的产品？”
“这个您也没叫我问啊！”
“你先在这等会。”
老板打电话叫来了山农，并把同样的任务用同样的表达方式交代给了山农。
山农也去市场转了很大一圈，手中拿着一张表格，向老板汇报：“报告老板，今天市场上主要有……等公司的蔬菜，蔬菜的种类分别有……，我做了张表，蔬菜产地、价格等信息都在上面，并做了分析，推测明天XXX蔬菜价格有可能会上涨，请您过目！”
老板看完表格，满意的点点头，向超逸看去。
超逸刚接触到老板的眼神，就连忙说：“老板，谢谢您，我知道该怎么做了！”

故事讲完了。
首先，我们看超逸和艾斯这类普通员工，虽然干活很麻利，但考虑问题较为片面，让调查蔬菜种类，就只知道考查种类，而山农和辛普森这类优秀员工则比较擅长全面思考，同样调查蔬菜，他们会综合调查，考虑种类和价格。老板眼光是雪亮的，所以工资孰高孰低的原因，不言自明。

相信看过前一篇文章《群落多样性之Alpha多样性（一）》的诸位大哥们，对故事中的人名似乎有些许耳熟吧。
不必回想了，相信我们赋值一下您肯定就明白了，以下是赋值代码，‘’#‘’后为代码的注释：
超逸=Chao1
艾斯=ACE
#Chao1和ACE是两个Alpha多样性指标，仅衡量样本中物种种类数量（Richness）。
山农=Shannon
#很多文献Shannon翻译为香农，个人认为如果翻译成山农的话，跟英文匹配度很高，更接近汉语拼音的发音规则，而且感觉有味道，较为接地气。
辛普森=Simpson
#Shannon 和 Simpson也是两个Alpha多样性指标，是把物种种类数量和各个物种的丰度全部考虑在内用以兼顾衡量样本中物种种类数量（Richness）和均匀度（Evenness）的指标

Chao1和ACE前面已经具体说过是怎么回事，今天的重点是从宏基因组微生态学的角度解释山农(Shannon)和辛普森（Simpson）为啥工资高。
在讲解之前，我们需要强调几点与生物多样性有关的几条纪律，且后面会反复用到。
----------------------------------------------几条纪律----------------------------------------------------
1. 先简单回顾下前一篇文章《群落多样性之Alpha多样性（一）》提到过的OTU和标记基因序列。
OTU，即可操作分类单元，这里要求很低，只需要知道“1个OTU对应一个物种，一个物种对应一个OTU”即可。这相当于出生于某地区的人（物种）对应的身份证号前六位（OTU），比如我的身份证号前六位220524，对应的就是我家乡，“物华天宝，人杰地灵”的吉林省通化市柳河县。
标记基因序列，以下简称序列，即测序得到的能够标记细菌个体的DNA序列。这里继续要求很低，只需要知道“一个序列对应一个细菌个体，一个细菌个体对应一个序列”即可，根据序列可知其OTU归属。这相当于一个人对应的一个身份证号，比如你知道我的身份证号220524***********5，你上网一查“220524”，百度显示“吉林省通化市柳河县”。
2. $S_{obs}$ ：观察到OTU的数量，即观察到的物种数。
3. $p_i=\frac{n_i}{N}$ ：第i个OTU在样本细菌总个体数中的占比，即物种相对丰度，也可以理解为在样本中随便揪出一个细菌个体，这个个体属于第个OTU或物种的概率。其中，第个OTU的序列数，即某物种的个体数，也是个观察值；：样本中的细菌个体总数。
------------------------------------------------------------------------------------------------------------------------------------------

步入正题，从宏基因组微生态学的角度，具体剖析一下，为什么山农（Shannon）的工资那么高？
Shannon的计算方式如下：
$H_{shannon} = - \sum_{i=1}^{S_{obs}} p_i ln p_i$
这个公式到底什么意思，需要把这个公式做个变换：
$H_{shannon} = - \sum_{i=1}^{S_{obs}} p_i ln p_i=- \sum_{i=1}^{S_{obs}} ln {p_i}^{p_i}=-(ln {p_1}^{p_1}+ln {p_2}^{p_2}+ln {p_3}^{p_3}+...+ln {p_{S_{obs}}}^{p_{S_{obs}}})$
$H_{shannon} =-(ln {p_1}^{p_1}{p_2}^{p_2}{p_3}^{p_3}...{p_{S_{obs}}}^{p_{S_{obs}}})=ln(\frac{1}{{p_1}^{p_1}{p_2}^{p_2}{p_3}^{p_3}...{p_{S_{obs}}}^{p_{S_{obs}}}} )=ln(\frac{1}{\prod\nolimits_{i=1}^{S_{obs}}{p_i}^{p_i}} )$
$ln \frac{n_i}{N}$ 是负数（ $n_i<N$ ），为符合人们的习惯，公式里加个负号将之修为正数。
根据上述公式，由于所有 $p_i$ 值的和等于1，所以即等于pi值的加权几何平均数，即
$\sqrt[\sum\nolimits_{i=1}^{S_{obs}}p_i]{\prod\nolimits_{i=1}^{S_{obs}}{p_i}^{p_i}}=\prod\nolimits_{i=1}^{S_{obs}}{p_i}^{p_i}$ ， $p_i$ 值本身用作几何权重(方程中的指数)。
因此括号内的项等于真正的多样性 $\frac{1}{\prod\nolimits_{i=1}^{S_{obs}}{p_i}^{p_i}}$ , $H_{shannon}$ 等于 $ln(\frac{1}{\prod\nolimits_{i=1}^{S_{obs}}{p_i}^{p_i}} )$ 。
为方便理解，这里介绍下加权几何平均数的意义，对这部分理解者可跳过此处。
----------------------------------------------------几何平均数的意义-----------------------------------------------
啥也不说，先上宝图。

引自知乎：https://www.zhihu.com/question/36176004/answer/139623544

假设a和b这两个数是两种细菌的个体数，它们构成一个菌群样本。他们的几何平均数是：
$G_2 =\sqrt{ab}$
结合上述宝图和中学数学知识可知，AE为a和b的几何平均数，AE这条垂线段越靠近B，a和b差距越大，即越不均匀。
极度均匀的情况是AE和OD重合，a=b，样本最均匀，样本的几何平均数AE最大。
如果菌群中存在3种菌，那么几何平均数为
$G_3=\sqrt[3]{abc}$
此时需要画个三维宝图解析一下，感兴趣不嫌麻烦的大哥可自绘，空间想象力好的大哥可直接脑补。
如果菌群中是n种菌，那么几何平均数为
$G_i=\sqrt[i]{abc...i}$ ,
由此可看出几何平均数可以反映数据的均匀度。
加权的意义只不过是把相同数据的频数组合放在一起而已，仅为计算方便，具体理解可见下式：
$G_3=\sqrt[2+1]{a_1^2a_2^{1}}$
$a_1$ 和 $a_2$ 的指数便是权，G就是加权几何平均数，这个式子也可画个3D的宝图解析。
如果是i维呢？
$G_{i}=\sqrt[i]{a_1a_2...a_i}$
如果是 $p_1+p_2+...+p_i$ 维呢？
$G_{p_1+p_2+...+p_i}=\sqrt[p_1+p_2+...+p_i]{a_1^{p_1}a_2^{p_2}...a_i^{p_i}}$
----------------------------------------------------------------------------------------------------------------------------
由此，我们知道了加权几何平均数可以反映样本的均一性，shannon指数最核心部分就是它。
为了更直观感受shannon指数，这里再介绍一种便于理解和感知的数学公式的方法，我称之为极限感知大法，也就是将一个极端数据带入公式去感知公式的意义。
首先，假设样本中所有物种的相对丰度都极端一致就是相等，那么所有 $p_i$ 值都等于 $\frac{1}{S_{obs}}$ ，因此Shannon取值为 $ln(S_{obs})$ 。
当类型丰度越不均匀，pi值的加权几何平均数越大，对应的Shannon越小。
然后，假设某群体中所有的个体都属于一个物种， $p_i$ 值等于1，代入公式，Shannon取值为0。

开篇的故事中除了山农，辛普森(Simpson)的工资也很高，接下来我们还是从宏基因组微生态学的角度说明下原因。
Simpson指数的计算方法如下：
$D_{simpson} =\sum_{i=1}^{S_{obs}}p_i^2$
这个公式相对来说比shannon更好理解一些。 $p_i=\frac{n_i}{N_i}$ ，可理解为从当前的菌群中随机挑选1个细菌，这个细菌属于第 $i$ 个物种的概率。那么 $p_i^2$ 就是从当前的菌群中随机挑选1个细菌，然后把这个细菌放回去，再从这个菌群中随机挑选1个细菌，这2个细菌都属于第i个物种的概率。然后把所有 $p_i^2$ 加到一起的意义就是在当前的菌群中随机挑选（有放回抽样）2个细菌，这两个细菌属于同一个物种的概率。
我们继续采用极限感知大法：
一个极端就是，让群落物种丰度极低且分布极端不均匀，只包含有1种细菌，其他细菌都是0，即 $n_1=N$ ，此时
$D_{simpson} =(\frac{n_1}{N}) ^2=1$
另一个极端，让群落物种丰度极端均匀，菌群包含 $S_{obs}$ 种细菌，每种细菌的个数是1，即 $S_{obs}=N$ ，此时
$D_{simpson} =(\frac{1}{N}) ^2\times S_{obs}=\frac{1}{S_{obs} } 或\frac{1}{N}$ 。
由此可见，Simpson值在0-1之间，值越小，多样性越高，均匀度均匀。
不过这怎么看着这么别扭呢，为了解决这个问题，通常用Inverse Simpson index(计算方法为 $1-D_{simpson}$ )或者Gini–Simpson index(计算方法为 $\frac{1}{D_{simpson} }$ ）替代Simpson。
搜底斯奈，这下能看出点规律了吧。
另外，对于Simpson指数的计算，还存在另外一个版本：
$D_{simpson} =\frac {\sum_{i=1}^{S_{obs}} {n_i \left ( n_i - 1 \right )}}{N \left( N-1 \right )}$
两个版本原理基本一致，唯一的不同就是这个版本在菌群种随机挑选2个细菌的时候是无放回抽样，而前面那个版本是有放回抽样。
那到底用哪个版本呢？
最科学建议是：想用哪个就用哪个！
为什么？
最充分的理由是：看心情！
如果你实在是有选择困难症，建议抛硬币占卜一下，看天意吧。

电视剧《甄嬛传》中甄嬛曾吟过一句诗，“逆风如解意，容易莫摧残。”，阶段性地俘获了雍正的心。
这句诗的大意是“北风如果能够理解梅花的心意，就请不要再摧残她了。”
可见解意很重要。对待公式也要充分解意，不然有人提问，答不上来，就是对公式的摧残。
极限感知大法固然能对公式有个初步的意会，然而真正更直观的解意可用计算和比较的方法。
比如有这么道判断题，Shannon和Simpson指数是否与细菌的绝对丰度有关？
通过公式的推导我们可以解答这类问题，不过用具体的数字代入计算会更直观一些。
如果对公式充分理解的话，计算部分可直接跳过。
---------------------------------------------规避各个因素后的计算------------------------------------------
这里我列举出一组数据：
A组：2, 3, 6, 9
B组：20, 30, 60, 90
C组：5, 5, 5, 5
D组：5, 5, 5, 5, 5
E组：4, 4, 4, 4, 4
F组：17, 1, 1, 1
求各组数据的Shannon和Simpson。
可直接代入公式。
A组。
$N_A=2+3+6+9=20$ $H_{shannon\_A}=\frac{2}{20} ln(\frac{2}{20} )+\frac{3}{20} ln(\frac{3}{20} )+\frac{6}{20} ln(\frac{6}{20} )+\frac{9}{20} ln(\frac{9}{20} )]=1.235$
$D_{simpson\_A} =\frac{2}{20}\times \frac{2}{20}+\frac{3}{20} \times \frac{3}{20}+\frac{6}{20} \times \frac{6}{20}+\frac{9}{20} \times \frac{9}{20}=0.325$
B组。
$N_B=20+30+60+90=200$
$H_{shannon\_B}=\frac{20}{200} ln(\frac{20}{200} )+\frac{30}{200} ln(\frac{30}{200} )+\frac{60}{200} ln(\frac{60}{200} )+\frac{90}{200} ln(\frac{90}{200} )]=1.235$

$D_{simpson\_B} =\frac{20}{200}\times \frac{20}{200}+\frac{30}{200} \times \frac{30}{200}+\frac{60}{200} \times \frac{60}{200}+\frac{90}{200} \times \frac{90}{200}=0.325$
数据占比相同的情况下，AB两组的两个参数相等，原因是这两个参数只与 $p_i=\frac{n_i}{N}$ 有关，与 $n_i$ 和 $N$ 两个绝对丰度无关。
C组。
$N_C=5+5+5+5=20$
$H_{shannon\_C}=-[\frac{5}{20} ln(\frac{5}{20})\times 4]=1.386$
$D_{simpson\_C} =(\frac{5}{20}\times \frac{5}{20})\times 4=0.25$
D组。
$N_D=5+5+5+5+5=25$
$H_{shannon\_D}=-[\frac{5}{25} ln(\frac{5}{25})\times 5]=1.609$
$D_{simpson\_D} =(\frac{4}{20}\times \frac{4}{20})\times 5=0.2$
E组。
$N_E=4+4+4+4+4=20$
$H_{shannon\_E}=-[\frac{4}{20} ln(\frac{4}{20})\times 5]=1.609$
$D_{simpson\_E} =(\frac{4}{20}\times \frac{4}{20})\times 5=0.2$
F组。
$N_F=17+1+1+1=20$
$H_{shannon\_F}=-[\frac{17}{20} ln(\frac{17}{20})+\frac{1}{20} ln(\frac{1}{20})\times 3]=0.5875$
$D_{simpson\_F} =\frac{17}{20}\times \frac{17}{20}+(\frac{1}{20}\times \frac{1}{20})\times 3=0.7300$
C和D规避了均匀度和 $n_i$ 的干扰，物种数量越多，Shannon越大，Simpson越小，与 $n_i$ 无关。
C和E规避了均匀度和 $N$ 的干扰，物种数量越多，Shannon越大，Simpson越小，与 $N$ 无关。
D和E基本上与A和B的比较情况一致，故不再多言。
C和F对比，N相同的情况下，不均匀的情况下，Shannon降低，Simpson升高。
注：这部分磨叽了点，本纠结要不要把这部分放上，还是不纠结了，一起充分感受一下。
--------------------------------------------------------------------------------------------------------------------------

综上所述可见，倘若菌群中几乎所有的个体都属于一个物种，而其他物种非常罕见，即使物种类别有很多，Shannon也会趋近于0，Simpsion也会趋于1。当数据集中只有一种类型时，Shannon正好等于0，Simpsion正好等于1。

末了，我们再回头想想前面那个小故事，为什么公司的老板没炒掉超逸（Chao1）和艾斯（ACE）呢？
因为经营一家公司，山农(Shannon)和辛普森(Simpson)这样全面考虑问题的优秀员工公司必然需要，但是超逸和艾斯这样，虽说考虑问题不全面但有一定执行力的员工我们也需要，分工不同。
我们做群落Alpha多样性分析也是一样，各类指标都有需求。
当我们只需要知道这堆细菌种有多少物种，Chao1和ACE足够；
想均匀度（Evenness）呢，那就是时候祭出Shannon和Simpson了。

参考文献：
[1] https://mothur.org/wiki/Shannon
[2] https://en.wikipedia.org/wiki/Diversity_index#cite_note-Simpson1949-7
[2] Simpson, E. H. (1949). Measurement of diversity. Nature.163: 688.
[2] http://www.countrysideinfo.co.uk/simpsons.htm

群落多样性之Alpha多样性（二）

猜你喜欢

热点阅读