腸道微生物測序研究中不同實驗來源差異的綜合評價以及應(yīng)用

徐珂琳， 莊 悅， 朱嗣博， 薛江莉， 蔣艷峰， 袁子宇，王久存， 索 晨,5)， 張鐵軍,5)， 呂 明,6)， 陳興棟*

(1)復(fù)旦大學(xué)公共衛(wèi)生學(xué)院生物統(tǒng)計學(xué)教研室， 上海 200032；2)復(fù)旦大學(xué)泰州健康科學(xué)研究院，江蘇， 泰州 225316；3)復(fù)旦大學(xué)生命科學(xué)學(xué)院現(xiàn)代人類學(xué)教育部重點實驗室， 上海 200433；4)復(fù)旦大學(xué)人類表型研究所， 上海 201203；5)復(fù)旦大學(xué)公共衛(wèi)生學(xué)院流行病學(xué)教研室， 上海 200032；6)山東大學(xué)齊魯醫(yī)院臨床流行病學(xué)研究中心， 濟南 250012)

人類胃腸道具有大量和高度多樣化的微生物群，是至今研究最多的生態(tài)系統(tǒng)之一。由于腸道微生物群在人體健康和內(nèi)穩(wěn)態(tài)方面的重要性，人們對其進行了詳細的研究。代謝紊亂[1-4]、炎癥性腸病[5, 6]和癌癥[7]等疾病的研究已經(jīng)揭示了腸道中特定的生物標(biāo)志物和生態(tài)系統(tǒng)。大規(guī)模的研究，例如人類微生物組計劃(Human Microbiome Project, HMP)和人類腸道宏基因組學(xué)(Metagenomics of the Human Intestinal Tract, MetaHIT)，為進一步的宏基因組學(xué)研究提供了必要的參考[8, 9]。這些研究采用16S rRNA擴增子測序(16S rRNA基因測序或16S測序)或宏基因組鳥槍(whole-metagenome shotgun, WMS)測序方法進行分類學(xué)圖譜分析。

近年來，WMS法被應(yīng)用于許多微生物組圖譜研究，在微生物組成和代謝通路的解釋方面展現(xiàn)出巨大的潛力[8, 10-13]。在WMS中，確定宏基因組分析和數(shù)據(jù)分析中的技術(shù)和生物學(xué)變異至關(guān)重要。此外，使用不同測序儀和文庫制備方法的腸道微生物測序研究中的分類學(xué)圖譜通常不一致，合并這類研究的結(jié)果仍然是一個亟待解決的問題[14, 15]。

微生物的樣本通常是收集受試者的糞便，并被隨機用于下游測序分析。研究人員發(fā)現(xiàn)，DNA提取前樣本的同質(zhì)化會影響樣品之間的變異度[15, 16]。然而，尚無數(shù)據(jù)表明樣本的不同采樣點會在多大程度上影響WMS的測序結(jié)果。

隨著測序方法的發(fā)展，越來越多高通量、短周期和低價格的產(chǎn)品可供研究人員選擇。Illumina公司生產(chǎn)的HiSeq 2500、MiSeq、HiSeq X10和最新發(fā)布的NovaSeq系列等測序儀在WMS的測序研究中具有壓倒性價格和產(chǎn)量的優(yōu)勢。這些測序儀之間是否存在顯著差異，以及哪種是當(dāng)前宏基因組測序?qū)嵺`的最佳選擇，仍有待進一步研究討論。

近來，針對微生物分類學(xué)圖譜的幾項大型研究的結(jié)果在一定程度上闡明了這個問題。例如，微生物組質(zhì)量控制計劃(Microbiome Quality Control Project, MBQC)利用16S rRNA測序方法比較了不同的糞便樣本保存方法、DNA提取試劑盒、序列長度和生物信息學(xué)工具，并指導(dǎo)了腸道微生物組研究的實驗設(shè)計[17]。另外，Costea等[15]比較了21種DNA提取試劑盒，并推薦了用于人類糞便樣本的標(biāo)準化方法。Walker等[18]的工作主要集中在16S rRNA PCR引物的選擇，以及測序結(jié)果與熒光原位雜交法(fluorescent in situ hybridization, FISH)間的比較。這些研究主要關(guān)注基于16S擴增的測序方法的比較，但基于WMS的測序方法仍被忽視。

SRA數(shù)據(jù)庫中已上傳超過179萬條“16S”和220萬條“宏基因組”序列(截至2021年8月17日)。然而，生成WMS數(shù)據(jù)庫的規(guī)模趨勢正在以驚人的速度增長。由于16S rRNA基因測序和宏基因組鳥槍測序基于各自獨特的方法，以及獲得不同方面的信息，這可能會進一步導(dǎo)致對結(jié)果的錯誤解釋[18, 19]。16S方法可識別和量化存在于所有細菌和古細菌中的標(biāo)志基因rRNA，并使用現(xiàn)有的大型公共數(shù)據(jù)集進行比較；而WMS法則測量樣本內(nèi)所有生物體的整個基因組。過去，人們對16S和WMS方法進行了一些研究，例如文庫試劑盒[20]之間的描述性比較，或使用樣本的α和β多樣性來評估群體內(nèi)受試者之間的差異[8]。然而，目前利用配對樣本直接比較16S和WMS之間的分類圖譜研究尚為缺乏，仍然存在兩者數(shù)據(jù)集能否在同一分析中進行比較的問題。

受微陣列質(zhì)量控制(microarray quality control, MAQC)、測序質(zhì)量控制(sequencing quality control, SEQC)和新發(fā)布的MBQC項目的啟發(fā)，本文利用來自4個健康供體糞便樣本的64個測序數(shù)據(jù)樣本，研究樣品同質(zhì)化、文庫制備方法和測序儀等因素的可比性問題[17, 21-23]。從一個大樣本的不同采樣點收集生物學(xué)重復(fù)，并使用WMS方法評估3種廣泛使用的測序儀(HiSeq 2500, HiSeq X10和NovaSeq 6000)。本文還應(yīng)用HiSeq 2500這一測序儀比較了16S擴增子和WMS文庫制備方法。最后，本文構(gòu)建了一種算法來提高16S和WMS方法之間的可比性。

1 材料與方法

1.1 樣品收集和DNA提取

本研究從復(fù)旦大學(xué)泰州健康科學(xué)研究院的泰州隊列中采集供體A、B、C、D 共4份糞便樣本。簡單來說，從4個樣本中每個收集2個獨立的取樣點作為生物學(xué)重復(fù)。同時，每個生物學(xué)重復(fù)經(jīng)過樣品同質(zhì)化產(chǎn)生2個相同的等分作為技術(shù)重復(fù)。16個樣本收集后立即轉(zhuǎn)移到-80 ℃冰箱中。本研究獲得了所有4名捐贈者的知情同意。采用宏基因組DNA提取試劑盒(Tiangen, China)和溶菌酶(Sigma-Aldrich, Canada)對DNA標(biāo)本進行純化，最終產(chǎn)生了16個WMS文庫和16個16S文庫。

1.2 16S rRNA基因文庫制備

本研究采用了V3-V4引物集，該引物集在16S測序研究中受到了廣泛應(yīng)用[23, 24]。根據(jù)Illumina的說明書設(shè)計并合成16S rRNA基因V3-V4擴增引物(插入片段469 bp, V3-V4擴增片段536 bp，編碼文庫613 bp)。正向引物337F: 5′-TCGTCGG CAGCGTCAGATGTGTATAAGAGACAGCCTACGGGN GGCWGCAG-3′；反向引物805R: 5′-GTCTCGTG GGCTCGGAGATGTGTATAAGAGACAGGACTACHVG GGTATCTAATCC-3′。引物稀釋后，在文庫制備前保存于-80 ℃條件。將12.5 ng微生物DNA、擴增子引物和HiFi HotStart ReadyMix (KAPA, USA)混合，再進行V3-V4 PCR試驗。熱循環(huán)條件包括95 ℃ 3 min；循環(huán)25次的95 ℃ 30 s、 55 ℃ 30 s和72 ℃ 30 s；72 ℃ 5 min和4 ℃保溫。使用V3-V4擴增子反應(yīng)的產(chǎn)物對用于Illumina測序的編碼庫PCR進行測序。利用N7和N5指數(shù)、95 ℃ 3 min、循環(huán)8次的95 ℃ 30 s、55 ℃ 30 s和72 ℃ 30 s；以及72 ℃放置5 min和保溫4 ℃的PCR條件下生成獨立文庫。用1.1 × AMPure XP beads(Beckmann Coulter, USA)純化文庫。所有文庫均通過Qubit 3.0和電泳定量或定性質(zhì)量控制，并進一步以摩爾濃度1∶1、最終濃度2 nmol/L匯集。

1.3 全宏基因組鳥槍法文庫制備

首先使用1 ng的微生物DNA片段進行分裂(Nextera XT試劑盒，Illumina，USA)。再使用NPM聚合酶(Illumina)和Illumina N5和N7指數(shù)特異性引物進行12個周期的PCR擴增DNA片段。文庫純化使用0.8 × AMPure XP beads(Beckmann Coulter, USA)。經(jīng)2100QC和Qubit定量的文庫圖譜分析，文庫以摩爾濃度為1∶1、最終濃度為2 nmol/L匯集。

1.4 二代測序技術(shù)

16S文庫使用Illumina HiSeq 2500測序儀(Illumina, San Diego, USA)上的250 bp雙端讀流單元進行測序。WMS文庫使用HiSeq 2500測序儀的250 bp雙端讀流單元，HiSeq X10和NovaSeq 6000測序儀的150 bp雙端讀流單元進行測序。供體B的1個樣本在HiSeq X10編碼中失敗，未納入分析。

1.5 生物信息學(xué)

所有64個原始FASTQ文件首先用FASTQC進行分析，以修剪過濾低質(zhì)量的堿基。共有的讀數(shù)使用PANDAseq和zcat命令拼接。16S文件用MALT分析生成RMA6格式文件，同時WMS文件用DIAMOND軟件計算生成DAA文件。使用MEGAN 6包對指定讀數(shù)的門、綱、目、科、屬和種水平進行測定，生成不同種系水平每個樣本的計數(shù)表。使用PANDAseq對16S配對端序列進行組裝，確定共有數(shù)量，通過對共有區(qū)域的錯誤進行修正，重構(gòu)并輸出整個序列[25]。將組裝好的長讀數(shù)進一步使用MALT進行處理，將讀數(shù)與SILVA rRNA數(shù)據(jù)庫(Release v128a)進行比對以用于分類分析。利用zcat命令行對WMS的對端序列進行組合以獲得長讀數(shù)。當(dāng)長讀數(shù)生成后，使用DIAMOND工具將讀數(shù)映射到NCBI參考宏基因組數(shù)據(jù)庫[26]。相關(guān)、聚類和LEfSe分析中使用的分類圖譜(WMS生成)是通過對非原核物種進行過濾處理得到。

1.6 經(jīng)驗貝葉斯方法

批次效應(yīng)、測序儀和文庫制備法的不同是微陣列研究領(lǐng)域中遇到的常見問題，特別是當(dāng)組合來自不同實驗的多批數(shù)據(jù)或?qū)嶒灢荒芡瑫r進行時。查閱回顧既往文獻，本文發(fā)現(xiàn)現(xiàn)有的調(diào)整方法，例如奇異值分解(single value decomposition，SVD)和距離加權(quán)判別法(distance weighted discriminant，DWD)不適用于處理批次大小很小的樣本。由此，本文開發(fā)了一種經(jīng)驗貝葉斯方法，(empirical Bayes, EB)用于提高16S和WMS方法[27]之間的群體可比性。EB方法已被廣泛應(yīng)用于微陣列數(shù)據(jù)分析，特別當(dāng)樣本規(guī)模較小(<25)時，它能夠穩(wěn)健地處理高維數(shù)據(jù)[27]。EB法主要是在計算中“借用信息”，利用標(biāo)準化數(shù)據(jù)和(非)參數(shù)先驗分布估計批次效應(yīng)參數(shù)，以期得到更好的估計或更穩(wěn)定的推斷。本文將EB法推廣到調(diào)整微陣列數(shù)據(jù)的批次效應(yīng)和提高文庫制備的可比性問題中，對16S和WMS檢測到的共有分類單元進行方法效應(yīng)校正，以避免僅由一種方法檢測到的分類單元所導(dǎo)致的校正算法的偏倚。用Yijt表示方法i中樣本j的第 t種分類單元的表達式值，假設(shè):

Yijt=αt+Xβt+γit+δitεijt

(1)

(2)

(3)

1.7 統(tǒng)計學(xué)分析

為了確定不同的提取方法之間哪些物種豐度有顯著差異，對2個樣本至少2個方案中豐度非零的物種應(yīng)用了Bray-Curtis檢驗。考慮到多重檢驗問題，對結(jié)果的P值進行Bonferroni校正，校正后P值低于0.05被認為有統(tǒng)計學(xué)意義。Shannon Weaver指數(shù)、Simpson倒數(shù)指數(shù)、PCoA和Bray-Curtis矩陣使用MEGAN6(版本6.10.3)計算。使用R (v3.5.0)軟件繪制Pearson相關(guān)性、Whiskers、相關(guān)熱圖和聚類圖。采用LEfSe進行線性判別分析(LDA)。

2 結(jié)果

2.1 實驗設(shè)計

糞便樣本由來自泰州隊列[28]的4個健康捐贈者(A、B、C和D)提供，提取標(biāo)本中的DNA 用于文庫制備。最終產(chǎn)生了用于WMS庫和16S庫的64個測序數(shù)據(jù)樣本(Fig.1A)。16S樣本進一步通過HiSeq 2500測序儀進行測序，而WMS樣本采用HiSeq2500、HiSeq X10和NovaSeq 6000三個測序儀進行測序。16S樣品采用PANDAseq[25]和MALT流程處理，而WMS則采用DIAMOND流程[26]。1個基于HiSeq X10的測序結(jié)果由于條形碼技術(shù)故障而丟失，最終63個樣品通過數(shù)據(jù)質(zhì)量控制。

Fig.1 Experimental design and data description (A) Flow chart of library preparation from 64 samples. (B) The rarefaction curves of sequenced 16S (n=16) and WMS (n=48) samples. (C) Number of average detected genera by 16S (74.16±1.69) was significantly larger than that of WMS protocol (52.46±2.19). *** P < 0.001. (D) Shannon weaver index of genus level between 16S and WMS protocol had no statistically significant differences (P>0.05). (E) Bray-curtis distance differences between 16s and WMS based on PCoA

2.2 16S和WMS數(shù)據(jù)概述

計算了每個樣本的稀疏度，并繪制了校準讀數(shù)和檢測到的屬級數(shù)量之間的關(guān)系。結(jié)果表明，樣品的測序結(jié)果到達了平臺期，這保證了進一步分析數(shù)據(jù)的可靠性(Fig.1B)。16S方法檢測到的平均屬級數(shù)量(n=74.16±1.69)顯著多于WMS方法(n=52.46±2.19)(Fig.1C,P<0.001)。在Shannon Weaver指數(shù)(α多樣性)方面，來自16S和WMS方法樣本間沒有發(fā)現(xiàn)統(tǒng)計學(xué)差異(Fig.1D)。然而，基于Bray-Curtis距離(β多樣性)的PCoA圖顯示了方法不同引起偏倚，這導(dǎo)致16S方法的聚類占優(yōu)勢，特別是在供體A、C和D中(Fig.1E)。然而，當(dāng)從UPGMA、PCoA和N-J樹中考慮WMS或16S結(jié)果時，大多數(shù)樣本顯示出供體來源的聚類。在WMS方法中，HiSeq 2000、NovaSeq 6000和HiSeq X10三個測序儀的結(jié)果緊密地聚類于技術(shù)重復(fù)，也就是說有些樣品是根據(jù)它們的技術(shù)重復(fù)而不是生物學(xué)重復(fù)聚類的。

2.3 使用WMS進行測序儀的可重復(fù)性和樣本的異質(zhì)性分析

本文采用Shannon Weaver指數(shù)對3個使用WMS的測序儀在門、綱、目、科、屬5個水平上進行了比較。測序儀在各等級之間未發(fā)現(xiàn)顯著性差異。等級越低，α多樣性越高(Fig.2A)。

為了確定糞便樣本的異質(zhì)性，使用Shannon指數(shù)對配對技術(shù)重復(fù)和生物學(xué)重復(fù)進行了比較。與生物學(xué)重復(fù)(Fig.2C, Pearsonr=0.69)相比，WMS方法顯示技術(shù)重復(fù)具有更高的相關(guān)性(Fig.2B, Pearsonr=0.94)。

Fig.2 Sequencer reproducibility and Intra-specimen heterogeneity analysis using WMS (A) Shannon Weaver index of three WMS sequencers in five hierarchical ranks were calculated and compared. (B,C) Technical replicates and biological replicates in WMS protocol were compared using Shannon index and Pearson correlation coefficient. ***P<0.001. (D) Bray-Curtis index of the technical replicates and biological replicates in all five hierarchical ranks. ***P<0.001

采用Bray-Curtis 指數(shù)分析不同取樣點(生物學(xué)重復(fù))與同一取樣點(技術(shù)重復(fù))之間的差異。配對Bray-Curtis距離的結(jié)果顯示，在所有5個等級中，3個測序儀的生物學(xué)重復(fù)差異大于技術(shù)重復(fù)(Fig.2D,P<0.001)。僅使用HiSeq 2500測序儀的16S rRNA基因測序數(shù)據(jù)中也顯示類似的結(jié)果(P<0.05)。本文進一步使用基于LDA的方法計算每個測序儀的豐度差異物種。結(jié)果顯示，無特殊的物種或測序儀，表明不同測序平臺測序結(jié)果具有一致性和可重復(fù)性。

2.4 文庫制備方法不同產(chǎn)生的差異

本文使用α多樣性進一步比較測序方法、測序儀在樣本間的差異。由WMS方法間的Pearson相關(guān)性，結(jié)果顯示，測序儀之間具有高度的一致性(Fig.3A)。然而，使用HiSeq測序儀的16S和使用3種測序儀的WMS方法之間的Shannon指數(shù)存在巨大的差異，尤其是在樣本A、B和D之間(Fig.3A,B)。

基于Bray-Curtis距離，本文發(fā)現(xiàn)，在屬水平上16S與WMS方法之間的差異遠遠大于技術(shù)重復(fù)、生物學(xué)重復(fù)或測序儀間的距離(Fig.3C,P<0.001)。令人驚訝的是，測序方法間的Bray-Curtis差異指數(shù)與獨立樣本間一樣大(0.59±0.05vs0.64±0.01,P=0.22)。

Fig.3 Library preparation induced dissimilarity (A) Comparison of paired samples using alpha diversity across 2 protocols and 3 sequencers. There was a higher consistence among sequencers than in protocols. (B) Shannon Weaver index of dissimilarities between 16S and WMS. The difference was large especially in sample A, B and D. (C) Compared with the dissimilarity between technical replicates, sequencers, biological replicates, protocols or samples. ns, not significant. ***P < 0.001

2.5 16S和WMS方法中的特定分類學(xué)圖譜

本文首先比較了所有5個等級中16S和WMS間的配對分類圖譜。與高級別分類相比，屬水平的Pearson相關(guān)性較低(Fig.4A)。為了進一步揭示在屬水平中檢測到的偏倚特征，本文繪制了維恩圖，以顯示每個測序儀和方法檢測到的共有部分以及獨特的分類單元(Fig.4B)。

Fig.4 Specific taxonomic profiles in 16S and WMS protocols (A) Pairwise taxonomic assignments between 16S and WMS protocol were compared with Pearson’s correlation. (B) Venn diagram presented overlapping and unique taxa by each sequencer and protocol. (C) Stacked bar plot exhibited all samples’ taxonomic assignments in the genus level

樣品的條形圖顯示，16S方法中Faecalibacterium(糞桿菌)和Megamonas(巨單胞菌)占優(yōu)勢菌，而WMS樣品則傾向于以Prevotella(普雷沃氏菌)為優(yōu)勢菌(Fig.4C)。16S方法顯示，83個屬水平的獨特類群(計數(shù)>1)，包括Kosakonia(科薩克氏菌)、CandidatusSoleaferrea(瘤胃菌科)和Peptococcus(消化球菌屬), 而WMS鑒定了70個特定屬，包括Dialister(小桿菌屬)、Olsenella(歐陸森氏菌屬)和Akkermensia(阿克曼氏菌屬)。除了獨特的檢出類群外，有偏特征顯示出它們對2種制備方法之一的偏好，這也可能導(dǎo)致了不可重復(fù)性。本文進一步利用HiSeq 2500數(shù)據(jù)，應(yīng)用線性判別分析檢測每個方法的人工生物標(biāo)志物[29]。WMS對Clostridia(梭菌)和Chlamydia(衣原體)的鑒定效果較好，而16S方法對Cyanobacteria(藍藻菌)和Rhodobacteria(紅藻菌)的鑒定效果較好。

2.6 提高群體可比性的算法

為了降低16S和WMS間的方法效應(yīng)，本文提出了一種經(jīng)驗貝葉斯算法來提高群體可比性。結(jié)果表明，經(jīng)驗貝葉斯算法顯著增強了16S和WMS數(shù)據(jù)集中共有屬之間的方法學(xué)相關(guān)性，從r=0.45提高到r=0.70(Fig.5A)。在應(yīng)用貝葉斯過程后，文庫制備方法的差異特征被消除，微生物群特征趨于相似(Fig.5B)。此外，該算法降低了由16S方法導(dǎo)致的PCoA圖的聚類偏倚，最終樣本依據(jù)供體來源聚集 (Fig.5C)。最后，本文采用同樣的方法對泰州隊列另外2名健康捐贈者的糞便樣本進行了測序，以驗證所提出的貝葉斯方法的有效性。由通過質(zhì)量控制的30個測序數(shù)據(jù)的結(jié)果，經(jīng)過貝葉斯算法校正，微生物群特征趨于相似(從r=0.37提高到r=0.59)，且PCoA聚類偏倚有所改善。據(jù)此，此貝葉斯算法能有效提高不同測序來源群體的可比性。

Fig.5 Empirical Bayesian algorithm to improve population-wide comparability (A) Pearson correlation coefficients of samples with 16S and WMS before and after using empirical Bayes algorithm. The algorithm enhanced the correlation in overlap genus level from r=0.45 to r=0.70. (B) Microbiota relative abundance and patterns in untreated, overlap and Bayesian stage. Differences in library preparation methods were eliminated. (C) The PCoA clustering plot based on 16S and WMS protocol from different donors. Almost all the samples are clustered by their donors after the Bayesian process

2.7 基于經(jīng)驗貝葉斯算法的外部數(shù)據(jù)集驗證

本文應(yīng)用經(jīng)驗貝葉斯算法融合分析了Dubin等[30]收集的10個WMS和16S配對腸道微生物樣本(PRJNA302832)的數(shù)據(jù)集。本文的算法顯著增強了樣本間的相關(guān)性(Fig.6A)。基于Bray-Curtis的PCoA也顯示，在消除了測序方法的偏倚后，同一供體的樣本成功聚類(Fig.6B,C)。

Fig.6 Validation of the Empirical Bayesian based algorithm We applied our algorithm to a published gut microbiota dataset with 10 WMS and 16S pairwise samples (PRJNA302832). (A) A significant enhancement of sample-to-sample correlations was obtained compared with untreated samples. (B,C) The samples were successfully clustered as their donor origin in the PCoA after protocol induced effect was removed

2.8 測序儀和文庫制備的一些事實

最后，本文總結(jié)了腸道微生物組研究生成測序文庫的最小成本。NovaSeq 6000向用戶展現(xiàn)了最具成本-效益的性能和時間價值，而HiSeq 2500和HiSeq X10則需要花費更多的時間和成本。考慮到良好的重現(xiàn)性和Simpson多樣性的要求，使用WMS文庫制備法和NovaSeq測序儀是最優(yōu)的選擇。本文在一臺配備128 Intel E7 4870 Quadcores, 1TB DDR4-2400 MHz, 10K rpm SAS 12Gb的服務(wù)器上測試并估計了其性能。

3 討論

本文對樣本的異質(zhì)性、測序儀和文庫制備方法等問題進行了腸道微生物測序的質(zhì)量控制研究。結(jié)果證實，每個樣本中的微生物組成在不同的取樣點上是不同的。3種測序儀在多樣性和分類豐度方面均得到了相似的測序結(jié)果。在16S和WMS方法生成的數(shù)據(jù)集中均能觀察到獨特的分類圖譜。最終，本文設(shè)計并提供了一個有效的模型來增強兩種測序數(shù)據(jù)之間的可比性。

同質(zhì)糞便樣品已被證明可為下游處理提供相同的試樣，從而產(chǎn)生可重復(fù)的數(shù)據(jù)[15,16]。然而，以往的研究使用同質(zhì)樣本來分析多組學(xué)數(shù)據(jù)或橫向比較多個變量，而不是顯示樣本本身異質(zhì)性的影響。Codling等發(fā)現(xiàn)，糞便與粘膜樣本之間存在差異，而Friswell等在胃腸道中發(fā)現(xiàn)了特異性位點的微生物群[31,32]。Hsieh等利用16S V4測序評估了微生物樣本的取樣點和處理方法[17,33]。與本文的結(jié)果類似，樣本同質(zhì)化顯示出減少數(shù)據(jù)集中個體內(nèi)變異的巨大潛力。同質(zhì)化或多次取樣對于減少下游測序的變異是必要的。

本文使用WMS方法比較了3種廣泛采用的測序儀HiSeq 2500、HiSeq X10和NovaSeq 6000。它們的測序結(jié)果在生物學(xué)多樣性和分類圖譜方面基本一致，具有較好的可重復(fù)性。MBQC使用16S rRNA方法比較了HiSeq和MiSeq，表明測序儀引起的差異小于樣本間差異，獲得了與本文相似的結(jié)論[17]。NovaSeq 6000在運行數(shù)據(jù)量、每樣本價格和測序速度方面，都比HiSeq 2500和HiSeq X10表現(xiàn)出壓倒性的競爭力和更高的性價比。從NovaSeq 6000輸出的數(shù)據(jù)高達1.5TB，這意味著通過多路傳輸可以同時測序500～800個樣本。因此，當(dāng)測序批次對測序結(jié)果有影響時，這將有利于控制批次效應(yīng)。

MBQC等研究使用了大規(guī)模的盲法樣本集，解決了人類微生物組測序方法、樣本處理方法和微生物數(shù)據(jù)處理計算流程等問題[15,17,18]。正如Costea等所提出的，在過去5年中，超過3 000項研究(從環(huán)境到生物醫(yī)學(xué)研究)對微生物群落進行了調(diào)查，產(chǎn)生了超過160 000份發(fā)表的16S和WMS宏基因組數(shù)據(jù)樣本[15,34]。這些結(jié)果強調(diào)了生物信息學(xué)工具和DNA提取試劑盒的特異性，而本文的研究重點是測序平臺、采樣點之間的差異和文庫制備方法的特異性。

本文的結(jié)果進一步表明，在所有量化的因素中，文庫制備方法的差異對觀察到的微生物組成的影響最大；這種影響甚至與樣本間差異導(dǎo)致的同樣大。16S rRNA測序是有偏差的，因為16S rRNA基因擴增不均勻，而WMS法可能因測序不夠深導(dǎo)致無法檢測到稀有物種[35]。Steven和Poretsky等的研究通過低維聚類或描述性比較揭示了16S基因測序與鳥槍測序的差異。然而，這些研究并未顯示配對方法間的距離和相關(guān)性比較[36,37]。在更大規(guī)模的宏基因組質(zhì)量控制研究中，人們已經(jīng)揭示了16S和WMS方法之間的差異。在Sinha等人的研究中，16S與WMS之間的Spearman相關(guān)系數(shù)為0.57～0.74，表現(xiàn)出中度正相關(guān)[17]。本文的研究結(jié)果表明，兩種測序方法的Pearson相關(guān)系數(shù)在0.2～0.8之間，主要受樣本來源的影響。Jovel等[19]就16S和WMS方法間的多樣性討論了類似問題，結(jié)果表明，模擬細菌在使用16S和WMS時有檢測偏倚。有趣的是，在組成更簡單的細菌群落中，16S與WMS結(jié)果的相對豐度模式相似。當(dāng)比較Simpson倒數(shù)指數(shù)在16S與WMS之間的相關(guān)性時，本文的數(shù)據(jù)支持這一結(jié)論。Hillmann等[38]在研究中，與16S方法相比，WMS的數(shù)據(jù)可以觀察到許多特定的物種。與本文的結(jié)論類似，16S和WMS在屬水平上的一致性高于種水平。

WMS在分類圖譜中顯示出更高的預(yù)測能力，即可以獲得物種甚至菌株水平的信息，而16S在屬水平上只能達到約80%的準確率[39,40]。16S和WMS可以相互補充，發(fā)揮各自的優(yōu)勢。16S方法適用于組織等含有或被高比例宿主DNA污染的樣本[41]，而WMS法不需要給定微生物群落，在糞便或唾液樣本中均有效。因此，為了更好地提高16S和WMS方法之間的可比性，需要用算法來減少測序方法不同引起的偏倚。此前，在這一領(lǐng)域的生物信息學(xué)研究主要是嘗試提高16S rRNA基因測序分類圖譜的預(yù)測能力和系統(tǒng)進化率，例如COMPASS和EMIRGE[42-44]。本文設(shè)計的經(jīng)驗貝葉斯法在減少方法帶來的差異方面以及在內(nèi)外部數(shù)據(jù)集中都有較好的表現(xiàn)，特別是在處理使用相同測序儀的配對樣本中，原因是這種方法將數(shù)據(jù)的不能消減效應(yīng)視為加法和乘法效應(yīng)，實際是一種均值中心化算法和尺度算法的混合[27]。

本文的研究在真陽性率和qPCR驗證方面仍存在一些局限。然而，F(xiàn)ACS或qPCR試驗的分類豐度推斷，也可能是由于雜交錯誤或多余PCR產(chǎn)物的過度擴增而導(dǎo)致的系統(tǒng)誤差[18]。單分子和長讀測序(例如MinION或PacBio)改進的宏基因組測序，有希望成為得到更精確分類圖譜的替代方案[45]。

本文的研究結(jié)果表明，同質(zhì)化是樣品DNA提取前的必要步驟。測序儀對分類學(xué)變異的貢獻小于文庫制備方法。經(jīng)驗貝葉斯方法提高了16S和WMS之間的群體可比性，這意味著它在進一步融合分析已發(fā)表的16S和微生物數(shù)據(jù)集方面具有強大的潛力。