王 強(qiáng),梁 玉,范小莉,張文馨,何 歡,戴九蘭,*

1 山東省林業(yè)外資與工程項(xiàng)目管理站,濟(jì)南 250014 2 山東省林業(yè)科學(xué)研究院,濟(jì)南 250014 3 山東大學(xué)環(huán)境研究院,青島 266237

微生物群落在全球生態(tài)循環(huán)中占據(jù)重要地位,也具有重要的環(huán)境指示作用[1],從不同層次和角度對微生物群落及生態(tài)特征進(jìn)行揭示一直是學(xué)界研究熱點(diǎn)[2-4]。自Garland等1991年首次將BIOLOG方法用于微生物群落生理特征檢測[5],數(shù)十年來該方法在微生物群落和生態(tài)研究中得到大量應(yīng)用,從碳源利用功能多樣性的角度刻畫了微生物群落生理特征[6-10]。BIOLOG平板包括革蘭氏陰性板、生態(tài)板、MT板、真菌板、SF-N和SF-P板等多種類型[11],不同平板針對研究對象不同,操作方法不同,但均能產(chǎn)生大量多維數(shù)據(jù)用以表征微生物群落的功能多樣性[12-14]。

在微生物生態(tài)研究中關(guān)于數(shù)據(jù)分析方法已有學(xué)者專門討論,但重點(diǎn)多在于理論層面,且對分子生物學(xué)范疇關(guān)注較多[15-17]。相比分子方法,BIOLOG方法產(chǎn)生數(shù)據(jù)局限于生理活性,受到關(guān)注較少；但由于其數(shù)據(jù)的多維度、時間依賴、環(huán)境相關(guān)等特點(diǎn),同樣需要對分析過程進(jìn)行深入探討,如數(shù)據(jù)讀取、整理、統(tǒng)計方法選擇,統(tǒng)計結(jié)果展示,統(tǒng)計結(jié)果解釋等。然而,由BIOLOG方法的相對小眾性,缺少專為此編寫的成熟數(shù)據(jù)處理軟件,大量數(shù)據(jù)處理需要研究者手動完成,耗時費(fèi)力,容易出錯；在應(yīng)用排序、聚類、差異性分析等方法時,也沒有一定標(biāo)準(zhǔn)或參考。不僅如此,在分析中常用的canoco軟件[18-19],也容易產(chǎn)生版權(quán)問題。已有研究對BIOLOG方法統(tǒng)計分析做過探討[20],其后不斷有研究者聯(lián)系討論統(tǒng)計方法選擇及實(shí)現(xiàn)方式,可見這仍是當(dāng)前微生物生態(tài)研究領(lǐng)域的一個難點(diǎn)。

R語言是為統(tǒng)計和數(shù)據(jù)處理而誕生的軟件,具有優(yōu)良的數(shù)據(jù)存儲處理、數(shù)組運(yùn)算、統(tǒng)計分析等功能,具有自由、免費(fèi)、開源等特點(diǎn)[21-22]。其語法簡潔優(yōu)雅,功能強(qiáng)大,在生物信息學(xué)、基因數(shù)據(jù)分析、空間分析、時間序列分析等領(lǐng)域均有優(yōu)秀表現(xiàn)。當(dāng)前R官方網(wǎng)站有15000多個擴(kuò)展程序包(packages),且在不斷增加,形成了非常成熟的應(yīng)用生態(tài)[21]。

本研究采用BIOLOG ECO板研究案例,以R語言為工具,從原始文本數(shù)據(jù)開始,全過程實(shí)現(xiàn)數(shù)據(jù)讀取、整理、統(tǒng)計分析,并針對分析結(jié)果,結(jié)合研究目標(biāo)給出統(tǒng)計學(xué)解釋。研究結(jié)果對BIOLOG微生物生態(tài)研究中統(tǒng)計分析方法選擇、提高數(shù)據(jù)處理效率有直接參考作用。

1 實(shí)驗(yàn)數(shù)據(jù)與分析流程

1.1 實(shí)驗(yàn)數(shù)據(jù)

采用數(shù)據(jù)為招遠(yuǎn)礦山周邊部分土壤樣品測量數(shù)據(jù)[23],選擇6個樣品,以樣品的BIOLOG ECO板培養(yǎng)數(shù)據(jù)和重金屬有效態(tài)含量數(shù)據(jù)進(jìn)行分析。

BIOLOG ECO板培養(yǎng)自12 h開始測量,每隔12 h測量一次,至132 h共測量11次。重金屬有效態(tài)數(shù)據(jù)作為環(huán)境指標(biāo),每個樣本測定3次重復(fù)。

BIOLOG ECO板讀數(shù)為txt文本數(shù)據(jù),重金屬有效態(tài)數(shù)據(jù)以csv格式存儲。

1.2 分析方法

將BIOLOG數(shù)據(jù)導(dǎo)入并整理,進(jìn)行一元特征指標(biāo)計算、組間差異檢驗(yàn)、排序、聚類等分析。應(yīng)用常見分析方法如AWCD(Average Well Color Density, 平均吸光度)、ANOVA(Analysis of Variance, 方差分析)、MANOVA(Multivariate analysis of variance, 多元方差分析)、NP-KW test(Kruskal-Wallis test, 克魯斯卡爾-沃利斯檢驗(yàn))、PCA(Principal Component Analysis, 主成分分析)、RDA(Redundancy Analysis, 冗余分析)、CA(Canonical Analysis, 典范分析)、CCA(Canonical Correspondence Analysis, 典范對應(yīng)分析)、PRC(Principal Response Curve, 主成分響應(yīng)曲線)、PCoA(Principal Coordinate Analysis, 主坐標(biāo)分析)、Db-RDA(distanc-based Redundancy Analysis, 基于距離的冗余分析)、NMDS(Nonmetric Multidimensional Scaling, 非度量多維尺度分析)、Mantel test(蒙特爾檢驗(yàn))等；各方法所需數(shù)據(jù)形式(原始數(shù)據(jù)或距離矩陣)及相互關(guān)系如圖1所示。其中排序分析包含限制性排序和非限制性排序；非限制性排序直接探索BIOLOG數(shù)據(jù)內(nèi)部規(guī)律,限制性排序則納入環(huán)境變量(重金屬含量數(shù)據(jù)),探討其作為自變量對因變量BIOLOG數(shù)據(jù)的影響[24]。本文主要對BIOLOG數(shù)據(jù)分析進(jìn)行探討,側(cè)重對群落結(jié)構(gòu)的解釋(排序),對環(huán)境因子自身的特征探索分析及模型擬合(Y=f(X))分析方法等內(nèi)容不做涉及。

2 數(shù)據(jù)導(dǎo)入與整理

2.1 數(shù)據(jù)導(dǎo)入

BIOLOG方法數(shù)據(jù)文件應(yīng)予以合理命名,便于批量讀取及整理。數(shù)據(jù)讀取后得到文本格式對象,要提取其中有效信息。本例中有效吸光度在文本的51—58行,構(gòu)建應(yīng)用bio.get()函數(shù),讀取得到含有66項(xiàng)BIOLOG數(shù)據(jù)的列表(list)；環(huán)境指標(biāo)數(shù)據(jù)通過read.csv()函數(shù)讀入(代碼見 https://yunpan.#/surl_yPE4EhVDxL3，下同)。

2.2 數(shù)據(jù)整理

列表是R語言中數(shù)據(jù)嵌套關(guān)系的一種形式,由于讀取過程而產(chǎn)生,內(nèi)容選擇可以采用$或[[]]兩種形式[25]。將列表中各樣本重新命名后轉(zhuǎn)化為數(shù)據(jù)框(data frame)形式進(jìn)行計算,31個碳源孔吸光度減去對照孔,負(fù)值歸零,得到用于分析的數(shù)據(jù)biolog.31；把樣本信息、重復(fù)信息、讀取時間信息以新列的形式加入數(shù)據(jù)集中,生成biolog.31.aid,為后續(xù)分析及作圖做準(zhǔn)備。

3 一元特征指數(shù)計算和組間差異檢驗(yàn)

3.1 一元特征指標(biāo)

本文一元特征指標(biāo)指用單一數(shù)字表征微生物群落功能特征的指標(biāo),如平均吸光度(Average Well Color Density, AWCD)、Shannon多樣性指數(shù)等。

AWCD為31個碳源吸光度的平均值,是一個有關(guān)研究幾乎必用的指標(biāo),衡量微生物利用碳源的能力[11,20]。R語言中可以使用apply()和mean()函數(shù),針對各行一次性完成。

除了以平均值表征碳源利用綜合能力之外,可以從對不同碳源利用能力多樣性的角度分析,計算多樣性指數(shù)[26-27]。大部分多樣性指數(shù)發(fā)源于植物群落研究,可以使用vegan包的vegdist()函數(shù)求取[28]。

一元特征指標(biāo)由單一測量時間信息得出,可以綜合多個時間研究其變化趨勢。圖2為培養(yǎng)期間AWCD值變化情況,同一個時間3個重復(fù)均以散點(diǎn)圖畫在圖上,之后進(jìn)行線性擬合?？梢园l(fā)現(xiàn),樣本2的AWCD值增加最快,表明群落碳源利用功能最強(qiáng)；其次為樣本3,樣本6、4、1的AWCD值依次降低但區(qū)分度不高,樣本5培養(yǎng)前期變化較小,后期超過了樣本1。

實(shí)際研究中,基于研究目標(biāo)不同,同樣也可以對多樣性指數(shù)等作圖,探索變化規(guī)律[29]。

由圖2可見,84小時附近各樣本間讀數(shù)區(qū)分度較好,可以對該時間點(diǎn)讀數(shù)信息進(jìn)行深入分析,如差異性檢驗(yàn)、排序、聚類等。

3.2 組間差異檢驗(yàn)

3.2.1一元指標(biāo)組間差異檢驗(yàn)

為驗(yàn)證84小時各樣本AWCD是否有差異,進(jìn)行方差分析,使用oneway.test()函數(shù)。結(jié)果表明,P值為0.1698,即84小時6個樣本AWCD值間不存在顯著差異。進(jìn)一步檢驗(yàn)數(shù)據(jù)是否符合正態(tài)分布,使用shapiro.test()函數(shù),P值為0.04154,即數(shù)據(jù)不滿足正態(tài)分布,因此須采用非參數(shù)檢驗(yàn)考察樣本間差異性是否顯著。

圖2 培養(yǎng)過程中AWCD變化趨勢Fig.2 Developing trend of AWCD in the culturing period

Kruskal-Wallis檢驗(yàn)是方差分析的一個非參數(shù)版本。非參數(shù),意味著不假設(shè)數(shù)據(jù)的正態(tài)分布形態(tài)背景,但它假設(shè)各組有相同形狀的分布,在BIOLOG數(shù)據(jù)差異性檢驗(yàn)中應(yīng)用較為合適[30]。使用kruskal.test()函數(shù)實(shí)現(xiàn),結(jié)果表明,P值為0.05547,接近0.05的閾值,體現(xiàn)一定差異顯著性。

3.2.2多元指標(biāo)組間差異檢驗(yàn)

圖3 使用QQ圖檢驗(yàn)BIOLOG數(shù)據(jù)多元正態(tài)性 Fig.3 QQ plot assessing multivariate normality of BIOLOG data

對多元數(shù)據(jù)的差異性檢驗(yàn)常應(yīng)用多元方差分析(multivariate analysis of variance,MANOVA)。MANOVA有兩個前提,一是多元正態(tài)性,二是方差齊性。第一個前提說明變量要滿足正態(tài)分布,可以用QQ圖來檢驗(yàn),如果符合正態(tài)性,則會呈現(xiàn)一條直線；第二個前提要求各組樣本協(xié)方差矩陣齊性,可以用M盒檢驗(yàn)完成[25]。

應(yīng)用QQ圖檢驗(yàn),發(fā)現(xiàn)數(shù)據(jù)不符合正態(tài)分布(圖3)。但MANOVA方法有一定分布耐受性,可以計算出結(jié)果,顯示數(shù)據(jù)樣本間、讀數(shù)間,均顯著差異。

對于不滿足正態(tài)分布的數(shù)據(jù)應(yīng)用非參數(shù)檢驗(yàn),采用rrcov包的Wilks.test()函數(shù),方法選擇rank；也可以用vegan包的adonis()函數(shù)。結(jié)果顯示,樣本間、讀數(shù)間均有顯著差異,表明各樣本間BIOLOG吸光度變化具有顯著差異。

4 多元分析

4.1 非限制性排序

4.1.1主成分分析

主成分分析法(Principal Component Analysis, PCA)不僅常見于BIOLOG數(shù)據(jù)分析中,在微生物生態(tài)其他類型數(shù)據(jù),如基因型、磷脂脂肪酸等分析中也經(jīng)常用到[26,31]。主成分分析有基于variance-covariance(scale1)和correlation(scale2)兩種方式[32]。第一種主要用于相同屬性的數(shù)據(jù)集,用于分析樣本間的差異；第二種主要用于不同屬性的數(shù)據(jù)集,如環(huán)境因子數(shù)據(jù)集,單位、數(shù)量級都不一樣,挖掘變量間的相互關(guān)系。

對84小時BIOLOG數(shù)據(jù)進(jìn)行主成分分析結(jié)果表明,前兩個主成分解釋63.062%的數(shù)據(jù)信息,雙標(biāo)圖(biplot)將樣本和指標(biāo)放在同一圖上,能夠清晰表現(xiàn)出各樣本之間的差異,本例中將同一樣本的3個重復(fù)相連表現(xiàn),紅色剪頭表征不同碳源吸光度的變化梯度,可以表現(xiàn)出不同樣本在碳源吸光度(A2—H4)變化梯度上的分布(圖4)。

圖4 對84小時吸光度進(jìn)行主成分分析(5個樣本重復(fù)以-1、-2、-3表示)Fig.4 Principe Component Analysis of BIOLOG data at 84 hour (-1, -2, -3 represent for replications of 5 samples)

要深入分析不同碳源對主成分的貢獻(xiàn)率,找到特征碳源,可以通過分析各變量在主成分的載荷(loading)完成[31]。可以使用vegan包中score()函數(shù)實(shí)現(xiàn)。

PCA是基于線性模型進(jìn)行的分析,需要數(shù)據(jù)基本服從正態(tài)分布,在分析前應(yīng)該對BIOLOG吸光度數(shù)據(jù)是否服從正態(tài)分布進(jìn)行檢驗(yàn)。前文對本例BIOLOG數(shù)據(jù)分析結(jié)果表明不服從多元正態(tài)分布,因此需要用更加合理的方式進(jìn)行排序分析。

圖5 BIOLOG數(shù)據(jù)84小時吸光度典范分析結(jié)果 Fig.5 Canonical Analysis of BIOLOG data at 84 hour A2—H4為31個碳源吸光度

4.1.2典范分析

典范分析(Canonical Analysis, CA)和PCA均為通過計算樣本間距離進(jìn)行的分析,不同在于PCA基于線性模型計算歐氏距離,而CA則基于單峰模型計算χ2距離(卡方距離,Chi-Square distance)。結(jié)果表明,CA應(yīng)用于BIOLOG數(shù)據(jù)效果不好,前兩個軸只能解釋39.536%的信息,且對樣本區(qū)分程度較差,如圖5所示,不同碳源指向沒有一致性。CA是基于單峰模型的分析[24],在微生物生態(tài)研究中應(yīng)用較為少見。

4.1.3主坐標(biāo)分析

主坐標(biāo)分析(Principal Coordinate Analysis, PCoA)與PCA、CA的區(qū)別在于距離計算方式。PCA基于線性歐式距離、CA基于χ2距離,而PCoA通過不同的距離矩陣(伴隨矩陣,association matrix)來進(jìn)行分析,從而增加距離計算靈活性,得到更加合理的結(jié)果。該方法在分子生物學(xué)應(yīng)用較多,基于基因序列計算距離,分析群落結(jié)構(gòu)在不同處理間的差異[33-35]。

圖6 基于比例差異和海靈格距離矩陣的主坐標(biāo)分析Fig.6 PCoA by Percentage difference and Hellinger distance dissimilarity matrices

分別基于比例差異矩陣(Percentage difference (aka Bray-Curtis) dissimilarity matrices)和海靈格距離矩陣(Hellinger distance)進(jìn)行PCoA分析(圖6)。結(jié)果表明,PCoA分析所得樣本間差異區(qū)分度較高,能發(fā)現(xiàn)明確的群落特征；其中比例差異矩陣更能體現(xiàn)出吸光度變化特點(diǎn),對于樣本2、5有更好的區(qū)分。海靈格距離矩陣在比例差異矩陣的基礎(chǔ)上進(jìn)一步進(jìn)行了開方轉(zhuǎn)換,降低了吸光度數(shù)據(jù)大小的權(quán)重,相對不利于發(fā)現(xiàn)BIOLOG數(shù)據(jù)特征?？傊?PCoA不依賴數(shù)據(jù)的分布及線性關(guān)系,受分布影響較小,結(jié)果更合理穩(wěn)健,相比PCA更適合用于BIOLOG數(shù)據(jù)分析。

4.1.4非度量多維尺度分析

PCoA在計算伴隨矩陣后,還需要計算確切的樣本間距離。但如果只需要考察樣本間相對位置,對距離要求不高,則可以采用非度量多維尺度分析/多元尺度分析(Nonmetric Multidimensional Scaling, NMDS)。該方法不僅能基于任意距離矩陣進(jìn)行計算,還能處理缺失值,應(yīng)用范圍廣泛。可以通過vegan包的metaMDA()函數(shù)實(shí)現(xiàn),若含有缺失值,可以使用isoMDS()函數(shù)[28]。

分析代碼見附件,結(jié)果如圖7所示,對表現(xiàn)樣本間特征差異同樣有較好的展現(xiàn)效果。該方法和PCoA方法常用于基于核酸序列的排序分析中[36],在基于BIOLOG板的研究中也有應(yīng)用[37-38]。本方法相比PCoA有更廣的應(yīng)用范圍,可以將BIOLOG數(shù)據(jù)之外的微生物群落信息一并納入分析,如生物量碳、基礎(chǔ)呼吸,測序信息等,理論上都可以一并納入分析。

圖7 非度量多維尺度(NMDS)分析 (NMDS1、NMDS2)Fig.7 Nonmetric Multidimensional Scaling (NMDS)

4.2 限制性排序

4.2.1環(huán)境指標(biāo)簡介

非限制性排序是基于BIOLOG數(shù)據(jù)特征進(jìn)行的排序,限制性排序則考慮環(huán)境指標(biāo)等外在因素影響。

本文采用環(huán)境指標(biāo)為土壤重金屬有效態(tài)含量,分析礦區(qū)周邊重金屬對微生物群落活性的影響,研究不同種類重金屬濃度梯度上微生物群落功能多樣性的變化。要說明的是,實(shí)際研究中,我們不僅考慮了重金屬有效態(tài),還測定了重金屬全量、若干土壤理化性質(zhì),并進(jìn)行了微生物培養(yǎng),綜合考察微生物群落模式和各方關(guān)系。本文為說明分析方法,僅以重金屬有效態(tài)作為環(huán)境因子代表進(jìn)行示例。

各樣本重金屬有效態(tài)含量見圖8,不同重金屬在不同樣本間含量不同,變化趨勢也不一樣。為發(fā)現(xiàn)重金屬含量與微生物之間關(guān)系,采用限制性排序進(jìn)一步分析。

圖8 各樣本不同重金屬有效態(tài)含量Fig.8 Available contents of heavy metal in soil samples

4.2.2冗余分析

冗余分析(Redundancy Analysis, RDA)是基于PCA的擴(kuò)展,通過自變量環(huán)境指標(biāo)對因變量生物指標(biāo)的多元回歸,揭示二者之間的聯(lián)系和關(guān)系[32]。一直有研究將本方法用于環(huán)境或其他微生物指標(biāo)與BIOLOG數(shù)據(jù)的關(guān)系分析[38-39]。

采用vegan程序包的rda()函數(shù),重金屬有效態(tài)作為自變量、84小時BIOLOG吸光度數(shù)據(jù)作為因變量進(jìn)行分析,可知重金屬有效態(tài)含量可解釋BIOLOG吸光度數(shù)據(jù)中51.89%的方差信息,作圖可直觀發(fā)現(xiàn)各樣本在重金屬含量梯度上的分布(圖9)?？梢姼鳂颖驹贑r和Fe的梯度上分散度最高,可以進(jìn)一步考慮這兩種環(huán)境因素在塑造形成目前微生物群落特征過程中發(fā)揮的影響。

RDA本質(zhì)上是自變量環(huán)境指標(biāo)對因變量生物指標(biāo)的多元回歸,自變量也不局限在一組數(shù)據(jù),可以同時考察兩組數(shù)據(jù)對BIOLOG的影響,如,將有效態(tài)重金屬和土壤理化性質(zhì)同時納入分析,即偏冗余分析方法(Partial RDA),和偏相關(guān)分析類似,可用于控制特定變量的影響。

4.2.3典范對應(yīng)分析

典范對應(yīng)分析(Canonical Correspondence Analysis, CCA)是在CA的基礎(chǔ)上,基于單峰模型計算χ2距離后,結(jié)合多元回歸分析因變量與自變量之間的關(guān)系。鑒于CA分析BIOLOG數(shù)據(jù)效果不佳,本文不再深入討論CCA方法。

實(shí)際研究中,有應(yīng)用CCA的案例。如在方差分解之后,進(jìn)一步分析不同環(huán)境因子與BIOLOG數(shù)據(jù)間的相關(guān)關(guān)系[40]；在基于分子方法的研究中也有應(yīng)用,考察化學(xué)因素與測序所得操作分類單元(Operational taxonomic units, OTUs)之間的關(guān)系[41]。

圖9 冗余分析 Fig.9 Redundancy Analysis

圖10 主響應(yīng)曲線分析Fig.10 Principle Response Curve Analysis

4.2.4主成分響應(yīng)曲線

主成分響應(yīng)曲線(Principal Response Curve, PRC)最初設(shè)計用來分析實(shí)驗(yàn)中具有重復(fù)多元數(shù)據(jù),可以考察不同設(shè)計因素對數(shù)據(jù)的影響。在BIOLOG數(shù)據(jù)的分析中,多孔吸光度為需要分析的多元數(shù)據(jù),時間和樣本均作為重復(fù)納入分析,該方法在BIOLOG研究中的應(yīng)用類似于AWCD培養(yǎng)曲線變化,但對數(shù)據(jù)的壓縮方式更加精細(xì),且能分辨不同碳源的貢獻(xiàn)[42]。本方法雖然歸入限制性排序,但不涉及環(huán)境因子,限制性是對時間及樣本而言。

PRC分析將所有BIOLOG測量數(shù)據(jù)均納入分析,按時間將不同樣本進(jìn)行區(qū)分。每個樣本的3個重復(fù)信息壓縮到一條線體現(xiàn)出彼此差異,右側(cè)不同碳源表征對差異貢獻(xiàn)的大小,分析結(jié)果見圖10?？梢园l(fā)現(xiàn)本圖與AWCD隨時間變化有一定相似性,其優(yōu)勢在于處理較為復(fù)雜的數(shù)據(jù),明確變量的貢獻(xiàn)。

4.2.5基于距離的冗余分析

RDA直接將BIOLOG數(shù)據(jù)和重金屬有效態(tài)數(shù)據(jù)納入分析,分別作為因變量和自變量,并在排序圖上以歐式距離反映樣本間差異。但有時歐氏距離并不適應(yīng)所研究生態(tài)假設(shè)或數(shù)據(jù)類型,如0—1數(shù)據(jù)無法計算歐氏距離,不能直接進(jìn)行RDA分析。針對這種情況,借鑒PCoA方法思想,先計算基于不同距離計算公式的伴隨矩陣,如Jaccard系數(shù)或S?rensen系數(shù)矩陣,然后通過PCoA提取主坐標(biāo),并將獲得的主坐標(biāo)矩陣作為響應(yīng)變量納入分析,結(jié)合環(huán)境變量,實(shí)現(xiàn)基于距離的冗余分析(distanc-based Redundancy Analysis, db-RDA)[32]。

BIOLOG數(shù)據(jù)不存在無法計算歐氏距離的情況,但仍可以應(yīng)用這種適用范圍較廣的方法,分析效果較好(圖11)。

圖11 基于距離的冗余分析Fig.11 Distance-Based Redundant Analysis (CAP, Canonical Analysis of Principal coordinates)

4.2.6環(huán)境向量擬合

環(huán)境向量擬合能得到與限制性排序表現(xiàn)形式相似的圖形結(jié)果,但內(nèi)在運(yùn)算過程不包括多元回歸和方差解釋的部分,而僅是在非限制性排序的基礎(chǔ)上,以多次置換(permutation)求最優(yōu)解的方式,將環(huán)境變量按照梯度進(jìn)行擬合[32]。

在環(huán)境數(shù)據(jù)能夠?qū)ι飻?shù)據(jù)很好的解釋時,應(yīng)優(yōu)先選用限制性排序,在環(huán)境變量能夠解釋的范疇內(nèi)解釋梯度下群落特征。但如果限制性排序能夠解釋的信息過低,則限制性排序失去意義,可以用環(huán)境向量擬合方法,對梯度下群落特征進(jìn)行探索(圖12)。

圖12 對3種非限制性排序的環(huán)境向量擬合Fig.12 Environmental fitting of 3 unconstrained ordination

4.3 聚類分析

聚類分析中,首先要通過變型后的數(shù)據(jù)建立伴隨矩陣,計算距離。伴隨矩陣建立以后,主要通過三種方法進(jìn)行聚類分析：層次聚類(hierarchical clustering)、k平均值聚類(k-means partitioning)和雙向聯(lián)結(jié)聚類(two-way joining),聚類結(jié)果常以樹的形式表現(xiàn)出來。

本例中選擇層次聚類中的全連接聚類方法(Complete Linkage Agglomerative Clustering)對84小時吸光度數(shù)據(jù)進(jìn)行分析示例,結(jié)果如圖13所示。需要明確,聚類分析本身是基于數(shù)據(jù)本身對觀測進(jìn)行分組的方法,并不是嚴(yán)格的統(tǒng)計檢驗(yàn)；聚類本身不能代表群落特征,必須結(jié)合環(huán)境指標(biāo)進(jìn)行分析；在表現(xiàn)數(shù)據(jù)信息方面,聚類反映的信息含量一般比排序要少。

圖13 BIOLOG 84小時吸光度聚類分析Fig.13 Cluster analysis of 84 h BIOLOG data

4.4 蒙特爾檢驗(yàn)

蒙特爾檢驗(yàn)(Mantel test)是對兩個矩陣相關(guān)關(guān)系的檢驗(yàn),檢驗(yàn)相關(guān)性是否顯著。這種方法多用于植物生態(tài)學(xué)研究,考察不同植物群落間相關(guān)關(guān)系。微生物生態(tài)研究中,在基于測序的方法中也有采用,如在限制性排序后,進(jìn)一步分析環(huán)境因子與細(xì)菌與真菌的關(guān)系[43]。

蒙特爾檢驗(yàn)的原假設(shè)是兩個矩陣間沒有相關(guān)關(guān)系,較小的P值證明拒絕原假設(shè),即矩陣間存在相關(guān)關(guān)系。對84小時BIOLOG吸光度數(shù)據(jù)和重金屬有效態(tài)數(shù)據(jù)進(jìn)行相關(guān)性分析,結(jié)果表明,不論采用Spearman還是Pearson方法,P值均為0.002,即二組數(shù)據(jù)均存在顯著相關(guān)關(guān)系。這也說明上述針對兩組數(shù)據(jù)做進(jìn)行的各種分析是合理有依據(jù)的。

5 討論

本文在R語言環(huán)境下進(jìn)行了AWCD值、多樣性指數(shù)計算,組間差異檢驗(yàn),排序、聚類和蒙特爾檢驗(yàn)等分析。其中AWCD值、多樣性指數(shù)描述樣本內(nèi)吸光度的特征,屬于生態(tài)研究中的樣方內(nèi)α多樣性,排序則進(jìn)行樣本間特征比較,考察β多樣性。在傳統(tǒng)生態(tài)學(xué)研究中,這些多樣性指標(biāo)可以用來考察物種更新速率、物種周轉(zhuǎn)變化等規(guī)律；微生物生態(tài)中得到這些指標(biāo),也能對群落研究有所借鑒和啟發(fā)。分析過程中快速、準(zhǔn)確的計算方法,能解放出更多精力用于思考科學(xué)問題,會對研究工作走向深入提供助力。

實(shí)際研究中還有很多可以用于BIOLOG數(shù)據(jù)分析的方法,篇幅所限沒有納入。如針對培養(yǎng)曲線進(jìn)行方程擬合并分析參數(shù)、Pathway 富集分析、相似性分析(Analysis of similarities,Anosim)、典范變異分析(Canonical variate analysis, CVA)、去趨勢對應(yīng)分析(Detrended correspondence analysis, DCA)、BIOLOG指標(biāo)作為因變量與環(huán)境指標(biāo)的建模等,均可結(jié)合研究目標(biāo)選用[44-45]。如對原生函數(shù)圖不滿意,可提取分析結(jié)果中關(guān)鍵數(shù)據(jù),使用R語言中g(shù)gplot2程序包[46],或?qū)隣rigin等作圖軟件中作圖。附件中也給出使用ggplot2包對PCA結(jié)果作圖示例,供參考。