2018, Vol. 16
地球上的各种环境中，如土壤，海洋或人体肠道等，都生存着成千上万的微生物，包含细菌、真菌、病毒和古菌等。微生物组(Microbiome)指某一特定生态环境中全部微生物及其遗传物质的总和[1-2]。人类微生物组的组成与功能具有非常丰富的多样性，并会随着环境的变化而动态变化。如在人类口腔中定殖的微生物就有超过600种，口腔微生物菌群失衡会引起龋齿、口臭等口腔健康问题[3]。随着高通量测序技术的快速发展，宏基因组方法可以系统研究不同环境中微生物群落的组成与多样性变化[4]。大量研究表明，微生物组结构变化与人类疾病，如糖尿病、癌症、消化不良等密切相关[5]。微生物组学已成为最热门的研究方向之一，目前有多项国际合作大的研究项目正在开展，如人类微生物组计划(HMP)、地球微生物组计划(EMP)和中国微生物组计划等[2，6]。
16S rRNA基因全长约1 542 bp，因在结构和功能上都具有高度的保守性，存在于所有的细菌中，因此常被用于微生物分类鉴定的标志物。16S rRNA基因有9个高变区(V1~V9)中，V4区是最准确且能识别到属水平的可变区，而V3~V4区测序长度(约465 bp)比V4区更长(约290 bp)，相较单独V4区测序更加准确[7]。新一代高通量测序平台Illumina MiSeq克服了传统高通量测序方法的读长较短与误差较大的缺陷，因而用于16S rRNA的V3~V4区测序可获得更准确的微生物物种分类信息。
高通量测序项目一般会产生海量的测序原始数据，对于原始数据的过滤筛选，加工以及统计分析，对计算机的软硬件要求较高，使没有计算机背景的研究人员难以分析数据和解释研究结果[8]。尤其是，常用的微生物组分析软件(如QIIME、Mothur等)需要在Linux系统下运行，而普通的生物研究人员完全掌握Linux系统的使用非常困难[9]。因此，如何使大多数科研人员能进行生物信息分析是生物大数据时代的一个挑战[10]。尽管现在已经开发有一些基于网页的数据分析平台，如MG-RAST[11]，SEED2[12]等，但是这些平台往往有许多限制，如数据不能太大，网络传输速度慢，分析任务等待时间长等。尤其是网络分析平台不适合那些需要对样本数据严格保密的项目。
微软在2016年推出Windows的Linux子系统－Windows Subsystem for Linux (WSL)，原生支持在windows10上运行Linux程序，包括大多数的Linux命令与工具，如bash, grep, ssh等。这样大量在Linux平台运行的生物信息软件都可以直接安装在WSL下运行。WSL使Linux软件安装与使用变得容易，只需要在Windows系统下，学习Linux系统的基本命令，不需要专门在电脑上或用虚拟机(如VirtualBox)安装Linux系统。
为方便研究人员分析微生物组数据，我们开发了一套在WSL-Ubuntu下分析样品16S rRNA基因扩增子高通量测序的简易流程。本流程基于Win10的Linux子系统，避免了安装Linux的复杂过程，并简化了QIIME[13]、VSEARCH[14]及相关软件的安装过程，只要运行一个安装脚本就自动安装并配置运行所需的环境变量，方便没有计算机背景的研究人员进行相关微生物组数据分析。
1 材料与方法 1.1 唾液样本采集与DNA测序数据本实验数据来源于一项针对在校大学生的口腔微生物研究项目。实验简要说明如下：首先在实验的前期阶段，我们进行样本信息的搜集，筛选和整理，后送实验室提取DNA和质量检验，其中，唾液样本的采集参考陈瑜等人的方法[15]，DNA的提取使用QIAGEN DNA Mini KIT试剂盒提取，经PCR扩增16S DNA的V3~V4区后，进行琼脂糖凝胶电泳和核酸蛋白定量检测来检测提取的DNA的纯度及浓度。对于检测合格的样本DNA送测序公司采用Illumina HiSeq测序平台进行16S rDNA双端(Paired-End)测序。公司返回的测序结果为FASTQ格式文件的压缩包，包含所有的样本正反两个方向的序列信息。
为方便学习本分析流程，本文只使用部分唾液样本(10个)的微生物组数据，测序数据可从本流程网站下载(http://www.ligene.cn/win16s/)。整个分析过程只需要约2 h。
1.2 电脑系统与软件安装 1.2.1 装有最新版Windows 10系统的笔记本电脑为了运行本分析流程，首先需要一台Windows 10系统的电脑，并且将该系统更新到最新版本(版本1 709及以上)。升级过程可以选择电脑自动升级，也可以在微软网站下载升级程序，进行手动升级(Windows易升软件: https://www.microsoft.com/en-us/software-download/windows10)。在系统升级完毕后，就可开启系统的Windows Subsystem for Linux(WSL)功能：在控制面板-“应用”-“程序与功能”-“启用或关闭Windows功能”-“Windows Subsystem for Linux”中开启“Windows Subsystem for Linux”, 从而开启WSL的功能。
1.2.2 WSL+Ubuntu安装完成WSL的开启之后，在Windows10自带的微软应用商店“Microsoft Store”中搜索并下载Ubuntu。当Ubuntu app下载完毕之后，点击“Launch”启动安装。然后窗口会提示你创建一个Linux的账户，以后在Linux环境下进行数据分析都是在这个用户名下。这个用户名和密码与电脑微软系统的账户和密码没有任何关系，因此不必设置为相同的。安装完成之后，即可以通过开始菜单的Ubuntu程序图标来运行Ubuntu Linux系统。
1.2.3 Xming安装目前WSL不支持图形界面，为了实现数据分析结果的可视化，需要安装一个Xming软件，通过Xming官网(http://sourceforge.net/projects/xming/)下载安装包(最新版本6.9.0.31)，进行安装即可。
1.2.4 安装win16s流程此流程是基于VSEARCH (v2.7.1)与QIIME (v1.9.1)等开源软件，因此要在WSL-Ubuntu的命令行终端中输入不同软件的安装指令，所有命令集见附件1(流程安装)。为方便初学者安装，我们提供一个bash脚本(win16s_install.sh)，下载此文件后，在ubuntu bash终端使用“bash win16s_install.sh”命令运行此脚本，即可一次安装所有流程运行依赖软件，并能自动配置合适的环境变量。
Win16s所有软件与相关脚本可以从网址下载：http://www.ligene.cn/win16s/.
1.3 win16s分析流程本分析流程的基本步骤如图 1所示。主要步骤为：样本双端测序运用VSEARCH软件进行序列拼接，OTU聚类，后通过比对RDP分类数据库注释物种，最后用QIIME相关脚本进行α/β多样性分析(如PCoA等)。
|
图 1 16S rRNA基因扩增子测序数据分析流程图 Figure 1 Flowchart of 16S rRNA gene amplicon data analysis pipeline |
本结果是基于对10个唾液样本的16S rRNA基因测序数据分析的教学过程。首先，每个学生按流程分析样本数据，可以得到样本的OTU聚类与物种分类信息；在查看样品的细菌种类后，统计每个样品的差异菌种。然后，根据问卷调查结果对样本进行分组：荤组与素组(见表 1)，通过微生物组多样性相关分析命令，统计不同组样品间的菌群结构差异。
|
表 1 样本信息表 Table 1 Basic information of the samples |
为了分析不同饮食习惯大学生的唾液微生物的种群差异性，我们首选通过调查问卷来筛选出合适的志愿者，然后再统一组织开展取样工作。取样人群的基本信息，包括样本编号，性别，年龄，BMI值，取样方法，所属组别(见表 1)。为了避免长时间的储存对样本产生的影响，对于取得的唾液样本将直接用于后续研究(具体参见材料与方法1.1部分)。
2.2 数据质控与预处理为保证数据质量和避免错误结果，16s rDNA扩增子测序数据需要质控及预处理后，再用于下游的分析。首先常用FastQC软件对测序数据进行质量评估，由于一般公司的测序报告文件会附带评估报告，质量太差会重测，此步非本实验必须步聚。然后，利用常用的NGS质控工具(如Trimmomatic[16])进行预处理，如过滤低质量reads(Q < 20)，切除read末端低质量碱基与未识别碱基N、切除测序接头(Adapters)、去除长度太短的序列等。
本实验数据为Illumina HiSeq2500平台Paired-end测序(PE250)，并已去除引物与Barcodes后的原始序列，每个样品有两个正反向测序reads数据文件(注意文件名书写格式，如SAA_R1.fastq与SAA_R2.fastq)。正向Reads(R1)和反向Reads(R2)两段序列有部分重叠可进行拼接。一般默认用fastq-join方法拼接，所得序列再过滤低质量reads与切除末端低质量碱基。由于16S rRNA的V3~V4区长度为465 bp，本实验过滤长度小于380 bp的序列片段。原始序列经过滤后，最终每个样本得到约15 000条有效reads(见表 2)用于后续分析，好的质量控制可以得到准确的分类结果。
|
|
表 2 唾液细菌测序的基本信息 Table 2 Basic sequencing information of bacteria in saliva |
OTU(Operational Taxonomic Units)是在系统发生学研究中，为了便于进行分析，人为给某一个分类单元(品系，种，属，分组等)设置的同一标志。要了解一个样品测序结果中的菌种、菌属等数目信息，就需要对序列进行聚类操作(Clustering)。通过聚类将序列按照彼此的相似性97%以上归为一个OTU。一般建议采用开放参考(Open reference)聚类方式，但是PICRUSt预测需要以默认Greengene数据库的封闭参考(Closed reference)聚类分析结果。
本实验对样品数据进行质量控制、序列拼接及去除嵌合体后，选择序列相似度值97%进行OTU聚类。细菌的16S rRNA分类一般使用97%相似度，而真核生物的18S rRNA分类一般使用96%相似度，使用不同的算法或相似度阀值可能得到完全不同的OTU数[17]。得到OTU后，采用RDP算法与Greengenes rDNA数据库比对，注释各OTU的分类单元。大多数OTU序列能注释到属水平，一部分OTU可以鉴定到种水平。
本实验共得到了212个OTU，在门和属两个水平上对同类的OTU进行分类鉴定，归属于12个门，94个属(见表 2)。通过OTU聚类及分类鉴定得到的结果，可以作为后续菌群统计分析的基础。
2.4 细菌种群的丰富度不同人体的口腔微生物种群存在着一定的差异，通过分析可以得到不同分类水平的细菌种群结构及其丰富度，可以清晰地得到不同样本在同一分类水平的差异。分析流程在进行可视化操作后，在taxa_summary文件夹中生成子文件夹taxa_summary_plots, 打开其中的文件bar_charts.html, 即可观察到各样品自门水平(L2)到属水平(L6)的分类及其组成比例。图 2是在门和属水平对菌群的差异性的结果。具体菌群分类与丰度信息可以在多样性分析输出目录(“taxa_summary/”)，查看在不同分类水平下的细菌种群的丰富度。由结果可知，部分菌种在所有样本中都存在；而少数菌种则只在部分样本中存在。
|
图 2 实验样品菌群分类与相对丰度 Figure 2 Taxonomy and relative abundance of the bacterial populations across samples |
通过对不同样本分组的种群丰富度分析，本实验共检出优势菌涉及12个门(见图 2(a))，其中Firmicutes、Proteobacteria、Bacteroidetes、Actinobacteria不存在显著性差异，而Fusobacteria存在显著性组间差异(wilcox.test, P < 0.05)。本实验共检出优势菌涉及94个属(见图 2(b))，荤食组中Streptococcus比素食组显著增加，而Fusobacterium比素食组显著减少(wilcox.test, P < 0.05)。
2.5 其它统计分析最后得到的biom文件(otus_table_tax.biom)可用于后续统计分析。通常可以计算两类多样性指数：Alpha多样性指数(α-diversity)与Beta多样性指数(β-diversity)。α多样性指数度量的是群落内部的物种多样性，β多样性指数度量的是群落间的物种多样性。按照生态学的观点，多样性高的生态群落，抵抗力和稳定性都可能更强。这里以QIIME流程中的core_diversity_analyses.py脚本分析，并通过最小测序量样本的reads(Counts)数作为-e参数(Rarefaction抽样数)。core_diversity_analyses.py分析运行一系列标准的alpha与beta多样性分析，可以查看分析输出文件夹(core_diversity/)下的文件查看相关结果。最后，可以将OTU表的biom格式文件转换成文本文件，用于其它统计分析软件的分析，如STAMP[18]或R软件包phyloseq, vena等。
2.5.1 样本菌群Alpha多样性分析本实验统计了样本α多样性分析中的不同多样性指数。图 3(a)是系统发育多样性指数(PD)，两组的多样性非常相似，虽然两条曲线(PD值)区分明显，但误差棒有重叠，即没有显著统计差异(P > 0.05)。另外，荤食组和素食组Simpson多样性指数分别为(0.94±0.04)与(0.95±0.02)，Shannon多样性指数分别为(5.54±0.63)与(5.78±0.30)，可见两组样本有着相似的多样性指数。
|
图 3 样本分组的菌群差异 Figure 3 Bacterial population difference between sample groups |
本实验利用QIIME的脚本对得到的OTU表做差异距离矩阵，后用主坐标分析(Principal coordinates analysis, PCoA)将距离矩阵可视化。PCoA结果如图 3(b)所示, 其中一个点代表一个样本，相同的颜色表示是同一组的样本，两个点之间的距离越近，说明两个样品的微生物群落差异越小。由PCoA分析结果(见图 3(b))可知两组样本均比较分散，主成分区分度不是特别明显，但在各自象限均占主要优势，说明两组样本的主成分有一定差异，但没有显著的差异。这可能与本研究的采样对象样本数量太少有关，还需要扩大取样人群范围进一步研究。
生活习惯和饮食习惯的不同，形成具有差异的口腔微生物群落结构，存在不同的优势菌和劣势菌。它们影响着口腔的环境，对人体造成不同的影响。通过分类统计分析研究不同组别的菌群组成差异，进一步研究菌群差异与样本之间存在的关系，是微生物组研究的意义和目的。
3 结论通过整合开源软件VSEARCH与QIIME等分析流程，我们开发了在Windows系统下分析16S rRNA基因扩增子测序数据的简易流程。本流程不仅适用于16S rRNA数据分析，还可以用于18S rRNA或ITS等标记基因的分析。以后随着研究的深入，我们将提供更多元化的分析脚本，不断更新现有分析流程，并同时在网上公布，供初学者学习微生物组分析方法[19]。
WSL是微软在Windows 10中推出的Linux子系统功能, 允许Linux ELF64二进制文件运行在Windows系统上。随着微软对WSL的改进，越来越多的生物信息软件将支持在WSL运行，从而方便广大生物学背景的人员在Windows系统下进行生物信息数据分析[20]。然而WSL也有一定局限，如它只支持64位程序，一些老的32位程序将不能在WSL下运行，如QIIME默认的uclust程序不能在WSL下编译成功。好在目前大多数32位程序都有相应的64位替代程序，如uclust、USEARCH[20]可以用免费的64位程序VSEARCH代替。另外，在WSL下，程序对硬盘的读写效率还比较低，程序运行的效率还是比在原生Linux系统下运行的速度慢。但WSL适合作为数据分析入门的学习环境，当研究人员习惯了Linux命令操作，可以很方便迁移到原生Linux计算机上进行实际项目数据的分析研究工作。
| [1] |
WHITESIDE S A, RAZVI H, DAVE S, et al. The microbiome of the urinary tract-a role beyond infection[J]. Journal of Urology, 2015, 12(2): 81. DOI:10.1016/j.juro.2015.09.053 (  0) |
| [2] |
刘双江, 施文元, 赵国屏. 中国微生物组计划:机遇与挑战[J]. 中国科学院院刊, 2017, 32(3): 241-250. LIU Shuangjiang, SHI Wenyuan, Zhao Guoping. China microbiome initiative:Opportunity and challenges[J]. Bulletin of Chinese Academy of Sciences, 2017, 32(3): 241-250. DOI:10.16418/j.issn.1000-3045.2017.03.004 ( 0) |
| [3] |
周学东, 徐健, 施文元. 人类口腔微生物组学研究:现状、挑战及机遇[J]. 微生物学报, 2017, 57(6): 806-821. ZHOU Xuedong, XU Jian, SHI Wenyuan. Human oral microbiome:Progress, challenge and opportunity[J]. Acta Microbiologica Sinica, 2017, 57(6): 806-821. DOI:10.13343/j.cnki.wsxb.20170063 ( 0) |
| [4] |
WOOLEY J C, GODZIK A, FRIEDBERG I. A primer on metagenomics[J]. PLoS Computational Biology, 2010, 6(2): e1000667. DOI:10.1371/journal.pcbi.1000667 ( 0) |
| [5] |
DEWEERDT S. Microbiome:Microbial mystery[J]. Nature, 2015, 521(7551): S10. DOI:10.1038/521S10a ( 0) |
| [6] |
DUBILIER N, MCFALL-NGAI M, ZHAO L P. Create a global microbiome effort[J]. Nature, 2015, 526: 631-634. DOI:10.1038/526631a ( 0) |
| [7] |
WU X, ZHANG H, CHEN J, et al. Comparison of the fecal microbiota of dholes high-throughput Illumina sequencing of the V3-V4 region of the 16S rRNA gene[J]. Applied Microbiology and Biotechnology, 2016, 100(8): 3577-3586. DOI:10.1007/s00253-015-7257-y ( 0) |
| [8] |
WANG J T, DALY J N, WILLNER D L, et al. Do you kiss your mother with that mouth? An authentic large-scale undergraduate research experience in mapping the human oral microbiome[J]. Journal of Microbiology and Biology Education, 2015, 16(1): 50-60. DOI:10.1128/jmbe.v16i1.816 ( 0) |
| [9] |
STEVENS J L, DEHORITY R, GOLLER C C. Using QⅡME to interpret environmental microbial communities in an upper level metagenomics course[J]. CourseSource, 2017(4): 1-6. DOI:10.24918/cs.2017.3 ( 0) |
| [10] |
MARX V. The big challenges of big data[J]. Nature, 2013, 498(7453): 255-260. DOI:10.1038/498255a ( 0) |
| [11] |
MEYER F, PAARMANN D, D'SOUZA M, et al. The metagenomics RAST server-a public resource for the automatic phylogenetic and functional analysis of metagenomes[J]. BMC Bioinformatics, 2008, 9(1): 386-386. DOI:10.1186/1471-2105-9-386 ( 0) |
| [12] |
VETROVSKY T, BALDRIAN P, MORAIS D. SEED 2:A user-friendly platform for amplicon high-throughput sequencing data analyses[J]. Bioinformatics, 2018. DOI:10.1093/bioinformatics/bty071 ( 0) |
| [13] |
ROGNES T, FLOURI T, NICHOLS B, et al. VSEARCH:A versatile open source tool for metagenomics[J]. PeerJ, 2016, 4(10): e2584. DOI:10.7287/peerj.preprints.2409v1 ( 0) |
| [14] |
CAPORASO J G, KUCZYNSKI J, STOMBAUGH J, et al. QⅡME allows analysis of high-throughput community sequencing data[J]. Nature Methods, 2010, 7(5): 335-336. DOI:10.1038/nmeth.f.303 ( 0) |
| [15] |
陈瑜, 刘婉薇, 李良芳, 等. 功能性消化不良患者唾液菌群的研究[J]. 重庆医学, 2017, 46(13): 1789-1796. CHEN Yu, LIU Wanwei, LI Liangfang, et al. Analysis on saliva microbiome in patients with functional dyspepsia[J]. Chongqing Medicine, 2017, 46(13): 1789-1796. DOI:10.3969/j.issn.1671-8348.2017.13.020 ( 0) |
| [16] |
BOLGER A M, LOHSE M, USADEL B. Trimmomatic:A flexible trimmer for Illumina sequence data[J]. Bioinformatics, 2014, 30(15): 2114-2120. DOI:10.1093/bioinformatics/btu170 ( 0) |
| [17] |
杨晨雪, 季吟秋, 王晓阳, 等. 基于18S rDNA的metabarcoading技术分析土壤小型动物多样性3种方法的比较[J]. 中国科学:生命科学, 2012, 42(12): 993-1001. YANG Chenxue, JI Yinqiu, WANG Xiaoyang, et al. Testing three pipelines for 18S rDNA-based metabarcoding of soil faunal diversity[J]. Scientia Sinica Vitae, 2012, 42(12): 993-1001. DOI:10.1360/zc2012-42-12-993 ( 0) |
| [18] |
PARKS D H, TYSON G W, HUGENHOLTZ P, et al. STAMP:Statistical analysis of taxonomic and functional profiles[J]. Bioinformatics, 2014, 30(21): 3123-3124. DOI:10.1007/978-1-4614-6418-1_780-1 ( 0) |
| [19] |
MORAIS D, ROESCH L, REDMILE-GORDON M, et al. BTW-bioinformatics through windows:An easy-to-install package to analyze marker gene data[J]. PeerJ, 2018, 6: e26581v1. DOI:10.7287/peerj.preprints.26581v1 ( 0) |
| [20] |
EDGAR R C. Search and clustering orders of magnitude faster than BLAST[J]. Bioinformatics, 2010, 26(19): 2460-2461. DOI:10.1093/bioinformatics/btq461. ( 0) |