群落多样性之Alpha多样性(二)

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

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

相信看过前一篇文章《群落多样性之Alpha多样性(一)》的诸位大哥们,对故事中的人名似乎有些许耳熟吧。
不必回想了,相信我们赋值一下您肯定就明白了,以下是赋值代码,‘’#‘’后为代码的注释:

阿尔法=Alpha
超逸=Chao1
艾斯=ACE (Abundance-based Coverage Estimator)
#Chao1和ACE是两个Alpha多样性指标,仅衡量样本中物种种类数量(Richness),。
山农=Shannon
#很多文献Shannon翻译为香农(克劳德·香农),个人认为如果翻译成山农的话,跟英文匹配度很高,更接近汉语拼音的发音规则,而且感觉有味道,较为接地气。
辛普森=Simpson
#Shannon 和 Simpson也是两个Alpha多样性指标,是把物种种类数量和各个物种的丰度全部考虑在内用以兼顾衡量样本中物种种类数量(Richness)和均匀度(Evenness)的指标。
所以,这里要祭出一个公式:richness+eveness=diversity

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的序列数,即某物种的个体数,也是个观察值;N:样本中的细菌个体总数。
---------------------------------------------------------------------------------------------------------------

步入正题,从宏基因组微生态学的角度,具体剖析一下,为什么山农(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,即等于p_i值的加权几何平均数,即
\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_1a_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=20H_{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_iN两个绝对丰度无关。
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足够;
想知道多样性(diversity)呢?那就是时候祭出Shannon和Simpson了!
不过呢?
有可能某位大哥会说,
“我就想考察均匀度(Evenness)怎么办?”
看来阿尔法老板要继续招聘新员工了,
欲知后事如何,请看下集《群落多样性之Alpha多样性(三)》。

备注:此文于2019年3月29日发于e源微生态

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

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 199,340评论 5 467
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 83,762评论 2 376
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 146,329评论 0 329
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 53,678评论 1 270
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 62,583评论 5 359
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 47,995评论 1 275
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,493评论 3 390
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,145评论 0 254
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,293评论 1 294
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,250评论 2 317
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,267评论 1 328
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 32,973评论 3 316
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,556评论 3 303
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,648评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,873评论 1 255
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,257评论 2 345
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 41,809评论 2 339

推荐阅读更多精彩内容