宏基因組大數(shù)據(jù)分析的質(zhì)量控制流程規(guī)范

鄭廣勇，楊楨，曹瑞芳，劉婉，李亦學(xué),，張國(guó)慶,

1.中國(guó)科學(xué)院上海生命科學(xué)研究院生物醫(yī)學(xué)大數(shù)據(jù)中心, 上海 200031；

2.上海生物信息技術(shù)研究中心, 上海 201203

1 引言

近年來(lái)，隨著二代測(cè)序技術(shù)及各種高通量組學(xué)技術(shù)的快速發(fā)展，生物醫(yī)學(xué)研究進(jìn)入大數(shù)據(jù)時(shí)代。實(shí)驗(yàn)技術(shù)和信息技術(shù)的發(fā)展，使生物醫(yī)學(xué)數(shù)據(jù)開(kāi)始從萬(wàn)億字節(jié)的TB級(jí)躍升到千萬(wàn)億字節(jié)的EB級(jí)，引發(fā)了后基因組時(shí)代的生物醫(yī)學(xué)研究的深刻變革。生物醫(yī)學(xué)研究已經(jīng)從第一范式的實(shí)驗(yàn)科學(xué)，或快或慢地進(jìn)入第二范式理論科學(xué)、第三范式計(jì)算科學(xué)，直到最新的第四范式數(shù)據(jù)密集型科學(xué)。數(shù)據(jù)質(zhì)量對(duì)大數(shù)據(jù)分析與挖掘的價(jià)值日益凸顯。目前，測(cè)序能力的提升極大地推動(dòng)了包括微生物組學(xué)在內(nèi)的各種組學(xué)的快速發(fā)展，進(jìn)而催生了大量以微生物為主要研究對(duì)象的國(guó)際合作項(xiàng)目，使得微生物組和微生物系統(tǒng)組成為生物醫(yī)學(xué)大數(shù)據(jù)的重點(diǎn)研究方向之一。微生物組是指存在于特定環(huán)境或生態(tài)系統(tǒng)中的微生物以及它們攜帶的遺傳信息和生物學(xué)功能的總和[1]。微生物組與人類健康有極為重要的關(guān)系，人體微生物組由數(shù)百至數(shù)千種不同的微生物組成，其細(xì)胞總數(shù)可達(dá)數(shù)萬(wàn)億之多，數(shù)量遠(yuǎn)超人體自身細(xì)胞數(shù)量，其中獨(dú)特的微生物基因數(shù)量多達(dá)2 000萬(wàn)個(gè)，數(shù)目遠(yuǎn)遠(yuǎn)超過(guò)人的基因數(shù)目（大約 2.5 萬(wàn)個(gè)）。通過(guò)對(duì)人體微生物組進(jìn)行研究，可以解析人類健康、營(yíng)養(yǎng)、代謝等方面的科學(xué)問(wèn)題。

2 國(guó)內(nèi)外相關(guān)研究

目前，宏基因組（metagenome）技術(shù)是微生物組研究的重要手段之一，該技術(shù)利用基因組學(xué)策略研究特定環(huán)境樣品中包含的全部微生物的遺傳組成及其功能模式[2]。宏基因組技術(shù)直接從環(huán)境樣品中提取DNA樣本，避開(kāi)了傳統(tǒng)的微生物分離培養(yǎng)方法，為研究和利用占微生物種類99%以上的不可培養(yǎng)的微生物提供了一種新的途徑和良好的策略。近年來(lái)，高通量測(cè)序技術(shù)的快速發(fā)展為微生物表型及其遺傳機(jī)制的探索提供了新的技術(shù)方案[3]。宏基因組技術(shù)在農(nóng)業(yè)、林業(yè)、環(huán)保、醫(yī)藥等領(lǐng)域都有著廣泛的應(yīng)用，提供了豐富的人體、動(dòng)植物、各類環(huán)境的科學(xué)研究數(shù)據(jù)[4-8]。2005年以來(lái)，以國(guó)際宏基因組聯(lián)盟為代表的微生物組研究計(jì)劃帶動(dòng)了很多國(guó)家的相關(guān)研究，例如美國(guó)的人類微生物組計(jì)劃（Human Microbiome Project，HMP）[9]和歐盟的人類腸道微生物聯(lián)盟（Metagenomics of the Human Intestinal Tract Consortium，MetaHIT）[10]。這些項(xiàng)目的實(shí)施，推動(dòng)了數(shù)據(jù)庫(kù)和數(shù)據(jù)挖掘等基礎(chǔ)研究工作的開(kāi)展，并對(duì)相關(guān)數(shù)據(jù)標(biāo)準(zhǔn)與質(zhì)量控制等提出了新的要求。相對(duì)傳統(tǒng)組學(xué)數(shù)據(jù)而言，宏基因組測(cè)序數(shù)據(jù)量大、涉及的微生物基因組復(fù)雜度高，而且可能存在大量未知物種，傳統(tǒng)的基于單一組學(xué)或單一物種的分析方法無(wú)法直接應(yīng)用于宏基因組研究，針對(duì)宏基因組數(shù)據(jù)的采集、整理、存儲(chǔ)以及后續(xù)研究分析等仍缺乏統(tǒng)一的規(guī)范。因此，宏基因組研究急需建立一套涵蓋樣本信息以及測(cè)序數(shù)據(jù)采集、整理、存儲(chǔ)、交換、分析的數(shù)據(jù)規(guī)范與標(biāo)準(zhǔn)。從數(shù)據(jù)流程看，宏基因組的分析過(guò)程包括元數(shù)據(jù)和測(cè)序數(shù)據(jù)質(zhì)量控制、微生物群落組成分析、群落功能生態(tài)分析、菌群差異功能分析等步驟（如圖1所示）。其中，元數(shù)據(jù)和測(cè)序數(shù)據(jù)質(zhì)量控制是后續(xù)一系列分析的基礎(chǔ)，直接影響整個(gè)分析的完整性和正確性，因而在整個(gè)數(shù)據(jù)分析過(guò)程中具有十分重要的意義。本文對(duì)宏基因組研究中的元數(shù)據(jù)及測(cè)序數(shù)據(jù)的預(yù)處理過(guò)程中的相關(guān)質(zhì)量控制標(biāo)準(zhǔn)進(jìn)行了簡(jiǎn)要探討，從而為宏基因組數(shù)據(jù)分析提供支撐。

3 宏基因組數(shù)據(jù)質(zhì)量控制

宏基因組的測(cè)序數(shù)據(jù)具有數(shù)據(jù)量大、復(fù)雜度高的特點(diǎn)，其數(shù)據(jù)分析有一些特定的要求，具體而言主要有以下幾個(gè)方面。

● 宏基因組研究中樣本收集、存儲(chǔ)、運(yùn)輸?shù)刃畔⒍紝?duì)微生物菌群差異分析有重要影響，因此需要在元數(shù)據(jù)中記錄這些信息，并在質(zhì)量控制流程中對(duì)元數(shù)據(jù)和測(cè)序數(shù)據(jù)進(jìn)行一致性檢查，以保證后續(xù)分析的正確性。

● 宏基因組研究通常依賴二代或三代測(cè)序技術(shù)，其測(cè)序速度比一代測(cè)序技術(shù)有顯著提升，但其測(cè)序長(zhǎng)度及準(zhǔn)確度卻有所下降，因此在宏基因組數(shù)據(jù)質(zhì)量控制流程中必須對(duì)數(shù)據(jù)進(jìn)行去噪、剔除接頭、切除無(wú)效序列等操作，以保證結(jié)果的可靠性。

● 與宿主相關(guān)的微生物樣品測(cè)序結(jié)果中通常包含一定數(shù)量的宿主基因片段，這些片段會(huì)對(duì)后續(xù)的微生物菌群分析產(chǎn)生干擾，因此在宏基因組數(shù)據(jù)質(zhì)量控制流程中必須去除這些宿主基因片段，以保證后續(xù)分析的有效性。

● 需要對(duì)樣本的測(cè)序深度進(jìn)行一定的評(píng)估，從而保證不同分組樣本差異分析的需求。

筆者在長(zhǎng)期的微生物組大數(shù)據(jù)分析過(guò)程中，根據(jù)宏基因組的數(shù)據(jù)特點(diǎn)和數(shù)據(jù)分析的要求，提出了一套完整的數(shù)據(jù)質(zhì)量控制流程（如圖2所示），具體包括：元數(shù)據(jù)和測(cè)序數(shù)據(jù)的一致性檢查、測(cè)序數(shù)據(jù)的質(zhì)量檢測(cè)、低質(zhì)量測(cè)序序列的過(guò)濾及切除、接頭序列及無(wú)關(guān)序列的剔除、宿主及污染序列的過(guò)濾、混合樣本的數(shù)據(jù)分割、樣本菌群的組成分析、不同分組樣本的主坐標(biāo)分析（principal coordinateanalysis，PCoA）這8個(gè)步驟，下面將對(duì)這些步驟進(jìn)行詳細(xì)論述。

圖1 宏基因組數(shù)據(jù)主要分析流程

圖2 宏基因組數(shù)據(jù)質(zhì)量控制流程

3.1 元數(shù)據(jù)和測(cè)序數(shù)據(jù)的一致性檢查

元數(shù)據(jù)是對(duì)數(shù)據(jù)進(jìn)行描述的數(shù)據(jù)（data about data），主要對(duì)數(shù)據(jù)的屬性進(jìn)行表述。元數(shù)據(jù)主要包括數(shù)據(jù)來(lái)源、數(shù)據(jù)收集整理模式以及數(shù)據(jù)可靠性等信息，是實(shí)驗(yàn)數(shù)據(jù)的重要組成部分。就宏基因組研究而言，元數(shù)據(jù)提供了項(xiàng)目技術(shù)設(shè)計(jì)、材料來(lái)源、實(shí)驗(yàn)方案、結(jié)果描述等至關(guān)重要的信息[11]。目前國(guó)際基因組標(biāo)準(zhǔn)聯(lián)盟（Genomic Standards Consortium，GSC）已開(kāi)發(fā)了針對(duì)基因組[12]、宏基因組[13]及擴(kuò)增子[14]研究的多種元數(shù)據(jù)標(biāo)準(zhǔn)框架，框架包括組學(xué)數(shù)據(jù)的檢查列表以及取樣環(huán)境描述包（environmental package）。其中，宏基因組數(shù)據(jù)標(biāo)準(zhǔn)框架的制定原則指出，人體樣本和環(huán)境樣本的最少信息集合應(yīng)當(dāng)以樣本為出發(fā)點(diǎn)，保證樣本的基因型和表型關(guān)聯(lián)分析、不同分組樣本的差異分析、樣本菌群組成的機(jī)理研究分析的需求，因此筆者建議宏基因組的元數(shù)據(jù)中應(yīng)當(dāng)包含如下基本信息。

● 人體樣本：個(gè)人基本生理信息、生活行為方式、膳食結(jié)構(gòu)、營(yíng)養(yǎng)狀況、既往病史等信息。

● 環(huán)境樣本：樣本獲取過(guò)程中的采集地點(diǎn)、大氣、水文、溫度、壓力、運(yùn)輸方法、存儲(chǔ)媒介等信息。

同時(shí)，筆者建議使用國(guó)際宏基因組數(shù)據(jù)標(biāo)準(zhǔn)框架制定的標(biāo)準(zhǔn)詞匯來(lái)描述元數(shù)據(jù)，從而給出精確的元數(shù)據(jù)語(yǔ)義信息，進(jìn)而建立可共享的、可被機(jī)器處理的本體支持，利于未來(lái)的數(shù)據(jù)資源整合。在元數(shù)據(jù)完整性檢查完成后，需要開(kāi)展元數(shù)據(jù)與測(cè)序數(shù)據(jù)的一致性檢查，重點(diǎn)檢查測(cè)序數(shù)據(jù)的樣本是否符合元數(shù)據(jù)描述的樣本，查看是否存在數(shù)據(jù)遺漏、差錯(cuò)匹配、錯(cuò)誤標(biāo)注等現(xiàn)象，這些都會(huì)對(duì)后續(xù)的數(shù)據(jù)分析造成嚴(yán)重的影響。這種一致性檢查是數(shù)據(jù)完整性、有效性的重要質(zhì)量控制步驟。

3.2 測(cè)序數(shù)據(jù)的質(zhì)量檢測(cè)

目前，宏基因組研究大都采用二代測(cè)序技術(shù)。以illumina為代表的二代測(cè)序技術(shù)基本都運(yùn)用邊合成邊測(cè)序的策略，在堿基鏈合成的過(guò)程中，隨著合成鏈的增長(zhǎng)，DNA聚合酶的效率會(huì)不斷下降，特異性也逐漸變差，從而造成堿基合成錯(cuò)誤率增高。此外，測(cè)序儀在開(kāi)始進(jìn)行合成反應(yīng)時(shí)，也會(huì)由于反應(yīng)不夠穩(wěn)定帶來(lái)質(zhì)量值的波動(dòng)。測(cè)序數(shù)據(jù)的質(zhì)量好壞會(huì)影響下游的分析，不同測(cè)序平臺(tái)的測(cè)序錯(cuò)誤率存在較大差異，因此在數(shù)據(jù)分析前，需要確定原始數(shù)據(jù)是通過(guò)哪種測(cè)序平臺(tái)產(chǎn)生的、其錯(cuò)誤分布規(guī)律如何、是否存在一定測(cè)序偏向性、是否受序列中堿基含量分布（GC含量）影響等[15]。目前對(duì)原始測(cè)序數(shù)據(jù)進(jìn)行質(zhì)量檢測(cè)，通常使用FastQC和PRINSEQ[16]這兩種檢測(cè)軟件包。FastQC軟件包既有圖形操作界面，又有命令行操作模式，能同時(shí)支持多種文件格式，包括FastQ、SAM和BAM格式等。其測(cè)序質(zhì)量報(bào)告主要包括：讀段（reads）各位置的堿基質(zhì)量值分布、堿基的總體質(zhì)量值分布、reads各個(gè)位置上堿基分布比例、GC含量分布、reads各個(gè)位置的非確定堿基數(shù)目、是否含有測(cè)序接頭序列等。PRINSEQ是另一款常用的測(cè)序質(zhì)量檢測(cè)軟件，其報(bào)告內(nèi)容主要包括：reads數(shù)目以及長(zhǎng)度頻率分布、堿基質(zhì)量頻率分布、序列復(fù)雜度、GC含量、重復(fù)序列及接頭信息等。

3.3 低質(zhì)量測(cè)序序列的過(guò)濾及切除

原始測(cè)序數(shù)據(jù)通常包含測(cè)序分值較低的堿基序列，其處理方式包括過(guò)濾和切除。對(duì)于全長(zhǎng)質(zhì)量都比較低的序列，可直接進(jìn)行過(guò)濾處理；而對(duì)于只有部分片段質(zhì)量較低的序列，則可以通過(guò)片段切除的方法來(lái)處理。最基本的切除方法為設(shè)定特定的切除長(zhǎng)度或剩余長(zhǎng)度，這種方法會(huì)去除部分測(cè)序質(zhì)量較好的片段，從而造成一定的信息損失，因此更為通用的方案是從序列任意一端開(kāi)始，逐個(gè)切除低于質(zhì)量閾值的堿基。如目前較為流行的低質(zhì)量序列切除與過(guò)濾軟件Trimmomatic[17]及PRINSEQ均支持從任意一端開(kāi)始切除低質(zhì)量堿基。另一種方案為滑窗策略，將窗口內(nèi)的堿基質(zhì)量與設(shè)定的閾值進(jìn)行比較，如果滑窗內(nèi)的堿基數(shù)值低于質(zhì)量閾值，則切除整個(gè)滑窗的堿基。用戶可根據(jù)實(shí)際情況設(shè)置滑窗大小及質(zhì)量閾值。另外，值得注意的是，測(cè)序時(shí)如果無(wú)法判定一個(gè)堿基是哪種堿基，通常標(biāo)記為N字符，不同的組裝比對(duì)軟件對(duì)于N字符的處理方式完全不同，有的用4種堿基隨機(jī)替代N字符，有的則是用固定的某個(gè)堿基替代N字符，由于N字符會(huì)導(dǎo)致錯(cuò)誤的組裝和比對(duì)，因此含N字符較多的序列同樣需要在組裝和比對(duì)前進(jìn)行切除或過(guò)濾操作。

3.4 接頭序列及無(wú)關(guān)序列的剔除

高通量測(cè)序過(guò)程一般會(huì)使用測(cè)序接頭序列，接頭序列的存在可能會(huì)導(dǎo)致基因組組裝和轉(zhuǎn)錄本組裝出現(xiàn)問(wèn)題，因此需要在分析數(shù)據(jù)之前予以去除。此外，其他的測(cè)序標(biāo)簽以及引物片段也需要去除。去除接頭及標(biāo)簽序列是一項(xiàng)比較困難的任務(wù)，首先，這些序列可能存在測(cè)序錯(cuò)誤，需要考慮應(yīng)對(duì)錯(cuò)配、插入缺失片段（indels）以及不確定堿基（N字符）的情況；其次，如果測(cè)序的目標(biāo)序列較短，reads可能會(huì)延伸到3’端的接頭序列。而這種“讀穿”的情況會(huì)導(dǎo)致reads中含有部分3’端的接頭序列無(wú)法被識(shí)別。此外，某些公共來(lái)源的測(cè)序數(shù)據(jù)可能根本無(wú)法知道接頭序列信息。目前較為常用的接頭序列切除軟件包括Trimmomatic[17]、TagCleaner[18]和Cutadapt等。這些軟件均可以應(yīng)對(duì)錯(cuò)配，并允許用戶指定測(cè)序和標(biāo)簽序列的最小重疊，TagCleaner還可以應(yīng)對(duì)indels以及不確定堿基。假如接頭序列未知，可以先用TagCleaner軟件預(yù)測(cè)，再進(jìn)行后續(xù)的切除。另外，值得注意的是，在測(cè)序文庫(kù)的制備過(guò)程中，聚合酶鏈?zhǔn)椒磻?yīng)（PCR）過(guò)度擴(kuò)增可能導(dǎo)致重復(fù)序列的出現(xiàn)，在后續(xù)分析過(guò)程中，這些重復(fù)序列應(yīng)當(dāng)剔除。此外，某些低復(fù)雜度的序列由于缺乏信息量，很難可靠地比對(duì)到參考序列上，因此也應(yīng)當(dāng)予以剔除。

3.5 宿主及污染序列的過(guò)濾

宏基因組測(cè)序是對(duì)樣本中所有DNA分子進(jìn)行測(cè)序，因此獲得的數(shù)據(jù)中可能含有不屬于微生物的序列，例如與人類健康相關(guān)的宏基因組研究多使用與人體相關(guān)的組織或樣本（如口腔菌斑、唾液、皮膚及糞便樣本等），這些樣本中可能存在人類基因組序列污染。此外，研究樣本也可能存在其他生物體或載體的序列污染，從而對(duì)后續(xù)分析造成影響，因此在質(zhì)量控制過(guò)程中需要剔除宿主序列以及可能的污染序列。最直接的方法就是把序列比對(duì)至宿主基因組及可能的污染源序列上，然后剔除這些序列。FastQ Screen、BWA[19]、Bowtie[20]、SOAP等工具可以把序列比對(duì)至用戶懷疑的污染源序列上，如果存在一致序列則予以剔除。

3.6 混合樣本的數(shù)據(jù)分割

測(cè)序時(shí)，為了區(qū)別不同樣本來(lái)源的序列，需要在待測(cè)序列中加入一段具有特定序列的編碼序列（barcode）。在質(zhì)量控制過(guò)程中，為了獲得不同樣本的完整數(shù)據(jù)，需要根據(jù)碼序列對(duì)混合樣本進(jìn)行分割，將一個(gè)FastQ文件中的序列分別存儲(chǔ)到多個(gè)樣本文件中。

3.7 樣本菌群的組成分析

在獲得測(cè)序數(shù)據(jù)之后，首要工作是對(duì)相關(guān)樣本進(jìn)行菌群組成分析。在早期宏基因組研究中，通常需要對(duì)測(cè)序數(shù)據(jù)進(jìn)行拼接，獲得contig序列，之后通過(guò)一些常用的序列比對(duì)算法（如BLAST算法、BLAT算法等），將拼接后的序列比對(duì)至微生物參考基因組序列，從而獲得與序列相關(guān)的物種分類信息。然而，不同于其他高等生物基因組的拼接，由于微生物的多樣性，宏基因組研究中測(cè)序數(shù)據(jù)的拼接往往存在一定的難度，導(dǎo)致序列準(zhǔn)確度不高。此外，傳統(tǒng)的BLAST算法在對(duì)海量的拼接后數(shù)據(jù)進(jìn)行比對(duì)時(shí)往往耗時(shí)較長(zhǎng)，BLAST算法速率為目前二代測(cè)序常用的比對(duì)軟件（如BWA、bowtie和SOAP等）速率的數(shù)百分之一至數(shù)十分之一，通常某些樣本完成所有序列比對(duì)可能花費(fèi)數(shù)天甚至更久的時(shí)間。因此筆者建議在實(shí)際分析過(guò)程中，使用二代測(cè)序的比對(duì)軟件，把測(cè)序片段（不經(jīng)過(guò)拼接）直接比對(duì)到參考物種的標(biāo)識(shí)基因，從而快速獲取物種分類信息及豐度信息，為后續(xù)的樣本主坐標(biāo)分析提供足夠的信息。例如，目前國(guó)際人類微生物組計(jì)劃項(xiàng)目中采用Metaphlan軟件進(jìn)行數(shù)據(jù)處理，該軟件可以直接將測(cè)序序列比對(duì)至微生物參考數(shù)據(jù)庫(kù)中（涵蓋了美國(guó)生物技術(shù)信息中心基因組數(shù)據(jù)庫(kù)中的2 887個(gè)微生物基因組數(shù)據(jù)），從而獲得測(cè)序樣本的菌群組成信息[21]。

3.8 不同分組樣本的主坐標(biāo)分析

主坐標(biāo)分析是一種微生物組學(xué)數(shù)據(jù)分析中常用的降維及可視化方法，主要用于研究數(shù)據(jù)相似性或差異性，在微生物群落相關(guān)研究中，不同樣本之間及不同環(huán)境之間微生物組成差異往往較大，所獲得的物種豐度表通常由稀疏矩陣構(gòu)成，因而不太適合使用其他組學(xué)數(shù)據(jù)通常采用的主成分分析方法。在主坐標(biāo)分析中，首先對(duì)物種豐度組成的距離矩陣進(jìn)行分解，獲得一系列的特征值和特征向量，然后對(duì)特征向量進(jìn)行排序，選擇前幾位的特征向量作為主要坐標(biāo)，并將樣品投影到這些向量的坐標(biāo)軸上進(jìn)行可視化展示。在對(duì)基于不同分類水平的物種豐度信息進(jìn)行PCoA時(shí)，樣本的物種組成越相似，它們?cè)赑CoA圖上的距離越小。通過(guò)不同分組樣本的主坐標(biāo)分析，可以檢測(cè)樣本的測(cè)序質(zhì)量以及深度是否提供了足夠的信息來(lái)區(qū)別不同的分組樣本，因而不同分組樣本的主坐標(biāo)分析是宏基因組數(shù)據(jù)分析質(zhì)量控制流程中非常重要的步驟。

4 宏基因組數(shù)據(jù)質(zhì)量控制應(yīng)用實(shí)例

本文以一套已公開(kāi)發(fā)表的腸道微生物宏基因組數(shù)據(jù)為例進(jìn)行質(zhì)量控制流程結(jié)果展示，該數(shù)據(jù)為正常人群與肝硬化人群腸道微生物隊(duì)列研究數(shù)據(jù)，相應(yīng)宏基因組測(cè)序數(shù)據(jù)下載自歐洲生物信息研究中心核酸數(shù)據(jù)庫(kù)（ID：ERP005860）。本文選取部分有代表性的樣本，對(duì)元數(shù)據(jù)和測(cè)序數(shù)據(jù)進(jìn)行包括元數(shù)據(jù)整理、元數(shù)據(jù)和測(cè)序數(shù)據(jù)一致性檢查、測(cè)序質(zhì)量評(píng)估、低質(zhì)量序列過(guò)濾、測(cè)序接頭片段去除、宿主基因剔除、樣本菌群組成分析、樣本主坐標(biāo)分析在內(nèi)的質(zhì)量控制檢測(cè)。其原始數(shù)據(jù)測(cè)序片段的質(zhì)量評(píng)估結(jié)果如圖3所示，個(gè)別樣本總體測(cè)序質(zhì)量偏低，后續(xù)分析需加以控制或予以剔除。

質(zhì)量控制前后樣本原始數(shù)據(jù)及干凈數(shù)據(jù)質(zhì)量統(tǒng)計(jì)見(jiàn)表1（括號(hào)內(nèi)為質(zhì)量控制后的數(shù)值）。質(zhì)量控制前后各樣本重復(fù)序列百分比、序列平均長(zhǎng)度及讀段總數(shù)等均有明顯變化，這表明了對(duì)原始測(cè)序數(shù)據(jù)進(jìn)行相關(guān)質(zhì)量控制的必要性。隨后，筆者對(duì)各樣本的菌群組成及微生物豐度進(jìn)行評(píng)估，并在此基礎(chǔ)上開(kāi)展了樣本的主坐標(biāo)分析，結(jié)果如圖4所示，正常樣本與疾病樣本可觀察到一定的分布差異，表明了質(zhì)量控制流程的有效性。

圖3 測(cè)序片段不同位置堿基質(zhì)量分布

5 宏基因組數(shù)據(jù)質(zhì)量控制應(yīng)用效果

上述宏基因組數(shù)據(jù)質(zhì)量控制實(shí)例表明，本文提出的質(zhì)量控制流程可以對(duì)宏基因組研究中的元數(shù)據(jù)和測(cè)序數(shù)據(jù)質(zhì)量進(jìn)行有效評(píng)估，從而為后續(xù)分析提供干凈的數(shù)據(jù)。具體而言主要有以下幾個(gè)方面。

● 元數(shù)據(jù)和測(cè)序數(shù)據(jù)的一致性檢查保證了取樣樣本和測(cè)序樣本的對(duì)應(yīng)關(guān)系，并賦予測(cè)序樣本表型信息，從而為測(cè)序樣本的基因型和表型關(guān)聯(lián)分析、不同分組樣本的差異分析、樣本菌群組成的機(jī)理研究分析提供支持。

表1 質(zhì)量控制前后序列質(zhì)量統(tǒng)計(jì)

● 宏基因組數(shù)據(jù)質(zhì)量控制流程中的質(zhì)量檢測(cè)、低質(zhì)量測(cè)序序列的過(guò)濾及切除、接頭序列及無(wú)關(guān)序列的剔除等步驟可以為后續(xù)的功能分析提供干凈的分析數(shù)據(jù)，保證結(jié)果的可靠性。

● 宿主及污染序列的過(guò)濾，可以去除宿主基因片段對(duì)后續(xù)菌群組成分析的影響，保證分析的有效性。

● 不同分組樣本的菌群組成和主坐標(biāo)分析，可以在一定程度上評(píng)估樣本的測(cè)序深度是否足夠，即在當(dāng)前的測(cè)序深度和質(zhì)量條件下是否提供了足夠的信息來(lái)滿足分組樣本差異分析的需求。

圖4 主坐標(biāo)分析結(jié)果

6 結(jié)束語(yǔ)

近年來(lái)，隨著各種微生物組項(xiàng)目的開(kāi)展，已有大量的宏基因組數(shù)據(jù)發(fā)布。如何充分利用和挖掘這些數(shù)據(jù)，對(duì)其進(jìn)行更為深入的二次分析，從而獲得新的發(fā)現(xiàn)，是一個(gè)極為重要的問(wèn)題。宏基因組數(shù)據(jù)多分散在不同的數(shù)據(jù)庫(kù)中或者不同的研究者手中，其數(shù)據(jù)收集和分析標(biāo)準(zhǔn)存在較大的差異，從而給數(shù)據(jù)的整合帶來(lái)較大的障礙。如果能從生態(tài)類型（biotype）、數(shù)據(jù)類型（datatype）等角度系統(tǒng)整合這些數(shù)據(jù)，并提供統(tǒng)一的質(zhì)量控制評(píng)價(jià)標(biāo)準(zhǔn)，將為宏基因組研究提供更大的便利。因此制定規(guī)范合理的數(shù)據(jù)質(zhì)量控制標(biāo)準(zhǔn)，從而有效整合這些不同來(lái)源和類型的宏基因組數(shù)據(jù)，是未來(lái)的發(fā)展方向和趨勢(shì)。本文對(duì)宏基因組數(shù)據(jù)分析過(guò)程中的質(zhì)量控制流程進(jìn)行了探討，這將加速宏基因組學(xué)的相關(guān)研究。在此基礎(chǔ)上，整合微生物的分類、進(jìn)化、生態(tài)以及相關(guān)組學(xué)的數(shù)據(jù)，構(gòu)建統(tǒng)一的微生物組數(shù)據(jù)倉(cāng)庫(kù)，并輔以微生物云服務(wù)平臺(tái)，將解決微生物組研究各種數(shù)據(jù)分散在不同地方的局面，為微生物組大數(shù)據(jù)提供科學(xué)的管理機(jī)制和運(yùn)行范式，從而為我國(guó)微生物組學(xué)研究提供技術(shù)支撐及基礎(chǔ)大數(shù)據(jù)平臺(tái)。

[1]張國(guó)慶, 寧康, 職曉陽(yáng), 等. 建設(shè)微生物組大數(shù)據(jù)中心發(fā)揮長(zhǎng)期科學(xué)影響[J]. 中國(guó)科學(xué)院院刊, 2017(3): 280-289.ZHANG G Q, NING K, ZHI X Y, et al.Development of comprehensive microbiome big data warehouse/center for long-term scientific impact[J].Bulletin of Chinese Acadamy of Sciences, 2017(3): 280-289.

[2]HANDELSMAN J, RONDON M R, BRADY S F, et al. Molecular biological access to the chemistry of unknown soil microbes: a new frontier for natural products[J]. Chemistry& Biology, 1998, 5(10): 245-249.

[3]WANG J, JIA H. Metagenome-wide association studies: fine-mining the microbiome[J]. Nat Rev Microbiol, 2016,14(8): 508-522.

[4]FANG H, CAI L, YU Y, et al. Metagenomic analysis reveals the prevalence of biodegradation genes for organic pollutants in activated sludge[J]. Bioresource Technology, 2013, 129(2): 209-218.

[5]DOS SANTOS H F, CURY J C, DO CARMO F L, et al. Mangrove bacterial diversity and the impact of oil contamination revealed by pyrosequencing: bacterial proxies for oil pollution[J]. Plos One, 2011,6(3): e16943.

[6]QIN N, YANG F, LI A, et al. Alterations of the human gut microbiome in liver cirrhosis[J].Nature, 2014, 513(7516): 59-64.

[7]HE Z, PICENO Y, DENG Y, et al. The phylogenetic composition and structure of soil microbial communities shifts in response to elevated carbon dioxide[J].Isme Journal, 2012, 6(2): 259-272.

[8]CHEN Y E, TSAO H. The skin microbiome: current perspectives and future challenges[J]. Journal of the American Academy of Dermatology, 2013,69(1): 143-155.

[9]TURNBAUGH P J, LEY R E, HAMADY M,et al. The human microbiome project[J].Nature, 2007, 449(7164): 804-810.

[10]EHRLICH S D. Metagenomics of the intestinal microbiota: potential applications[J]. Gastroenterologie Clinique Et Biologique, 2010, 34(4S1):S23-S28.

[11]PAGANI I, LIOLIOS K, JANSSON J,et al. The Genomes OnLine Database(GOLD) v.4: status of genomic and metagenomic projects and their associated metadata[J]. Nucleic Acids Research, 2012(40): 571-579.

[12]FIELD D, GARRITY G, GRAY T, et al.The minimum information about a genome sequence (MIGS) specification[J]. Nature Biotechnology, 2008, 26(5): 541-547.

[13]KOTTMANN R, GRAY T, MURPHY,et al. A standard MIGS/MIMS compliant XML Schema: toward the development of the Genomic Contextual Data Markup Language (GCDML)[J]. Omics-a Journal of Integrative Biology, 2008, 12(2): 115-121.

[14]YILMAZ P, KOTTMANN R, FIELD D, et al.Minimum information about a marker gene sequence (MIMARKS) and minimum information about any (x) sequence (MIxS)specifications[J]. Nature Biotechnology,2011, 29(5): 415-420.

[15]CARAUS I, ALSUWAILEM A A, NADON R, et al. Detecting and overcoming systematic bias in high-throughput screening technologies: a comprehensive review of practical issues and methodological solutions[J]. Briefings in Bioinformatics, 2015, 16(6): 97416.

[16]SCHMIEDER R, EDWARDS R. Quality control and preprocessing of metagenomic datasets[J]. Bioinformatics, 2011, 27(6):863-864.

[17]BOLGER A M, LOHSE M, USADEL B.Trimmomatic: a flexible trimmer for Illumina sequence data[J]. Bioinformatics,2014, 30(15): 2114-2120.

[18]SCHMIEDER R, YAN W L, ROHWER F,et al. TagCleaner: identification and removal of tag sequences from genomic and metagenomic datasets[J]. Bmc Bioinformatics, 2010, 11(1): 1-14.

[19]LI H, DURBIN R. Fast and accurate longread alignment with Burrows-Wheeler transform[J]. Bioinformatics, 2010, 26(5):589-595.

[20]LANGMEAD B, TRAPNELL C, POP M,et al. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome[J]. Genome Biology, 2009,10(3): R25.

[21]SEGATA N, WALDRON L, BALLARINI A,et al. Metagenomic microbial community profiling using unique clade-specific marker genes[J]. Nature Method, 2012, 9(8):811-814.