计量宏基因组学数据分析的方法及进展
本文转载自“艾美达”，己获授权。
《计量宏基因组学数据分析的方法及进展》是李昂老师撰写的文章，本文共分成6个部分：
1计量宏基因组学的通用研究路线
2非冗余参考基因集的构建算法和质量控制
3基于非冗余基因菌群的丰度构建
4 宏基因组靶基因测序分析方法的进展
5机器学习算法在宏基因组临床研究中的应用。
6展望
我读了李昂老师的文章之后，感觉受益匪浅，我感觉好的东西需要与大家一起分享探讨共同进步，因此向李老师申请授权分享，获得他的准许，因此就把这篇文章分享在这里。
本综述重点讨论了非冗余参考基因集和菌群丰度构建这两部分的分析算法的进展，另外还讨论了宏基因组靶基因数据分析的进展和机器学习算法在宏基因组研究中的应用。
1计量宏基因组学的通用研究路线
宏基因组学的概念最早由Handlesman 1等人在研究土壤微生物群落时提出。最早的宏基因组集中在对环境样品的研究，如空气，水，土壤等，需要提取出高质量微生物的DNA片段，然后转染至BAC或者Cosmid载体构建DNA表达文库和筛选环境，以期待获得某一类功能基因和菌群2，3，从此，研究者们突破了环境微生物难以分离得到的限制。VenterJC 4 等人率先使用sanger测序技术对马尾藻海进行研究,从海水中发现了2000多种不同的细菌，其中148种从未被发现过，此研究不再专注于某一类功能基因，而是对整个群落的功能基因进行描述，从此奠定了宏基因组学研究的基本思路。随着二代测序技术的发展，Qin5等人在2010年完成了第一个人体肠道菌群的图谱，从124个欧洲人肠道菌群中一共找到了330万个微生物基因，研究表明人体肠道菌群可能有1000到1500种细菌，这项研究标志着“计量基因组学”时代的开始，其核心思想是将一切数据量化，包括菌群的组成，功能基因的组成，菌群之间的关系等，并在量化数据的基础上，分析菌群和人体健康之间的关系。
计量宏基因组的一般研究路线为: 1) 样品预处理和DNA提取。通常用离心法，从粪便，土壤等样品中得到菌群，破壁后提取DNA。2) 高通量测序。构建DNA文库后，用高通量测序仪进行测序，得到基因组序列。3)数据质量控制。去除DNA文库中的人工序列，低质量序列，宿主污染序列，接头污染序列等，以得到高质量测序数据用于后续分析。4) 测序序列组装。将得到的高质量测序序列进行从头组装(de-brujin de novo assemble)，得到基因组片段(genome draft)。5) 非冗余参考基因集合构建。将每个样品的基因组片段进行基因预测，然后把所有的基因进行去冗余，将测序序列和参考基因集比较后得到菌群基因组成。6) 菌群丰度组成构建。将高通量测序数据和公开的微生物基因组，包括细菌，真菌，病毒等比较后，可以得到菌群的丰度组成。7)功能基因丰度组成构建。将参考基因集用COG , eggNOG, KEGG等数据库进行注释，将属于同一个功能的基因丰度加起来，可以得到功能基因的组成。以上分析完成后得到的结果为: 菌群组成，基因组成，功能基因组成, 然后在这3类组成的基础上进行下游数据分析。研究路线为图1所示。
图1 计量宏基因组学的数据分析流程
2非冗余参考基因集的构建算法和质量控制
对于单个样品的样品，测序后会得到不同菌株的基因组，而不同菌株之间存在高相似度的直系同源基因，如果要得到“某个菌群基因在单个样品中的分布”，就必须对多个直系同源基因进行去重复，只保留一个作为代表基因。如果是多个样品，要得到“某个菌群基因在所有样品中”的丰度分布，就必须对肠道菌群基因进行定义，以保证该基因在同一个研究中的所有样品，甚至不同研究的所有样品，定义是一致的，这样才能进行不同样品的在基因组成丰度上进行比较。同时，由于测序的差异，以及不同菌株之间存在直系同源基因，又允许基因在碱基上存在差异（参数通常设置为95%相似度，90%覆盖度）。因此，去冗余得到的“非冗余参考基因集”(non redundant geneset) ，解决的是两个问题：1）不同菌种之间的冗余基因 2) 不同样品之间共有的冗余基因。
非冗余参考基因集最早由Qin 5等提出，从124个肠道菌群宏基因组数据中，一共得到了14,048,045个ORF，从中构建了330万个非冗余基因，由于个体之间的差异很大，所以大部分基因，只在很少的样品中出现，结果发现2,375,655个基因只在20%的样品中全部出现，只有294,110个基因在50%的样品中全部出现。随后对糖尿病，肥胖，肝硬化等疾病的研究6, 7, 8，都对该基因集进行了更新，并在更新的基因集的基础上对疾病和肠道菌群的关联进行了研究。在新发布的宏基因组测序数据的基础上，Li Junhua9等构建了目前最大的一个肠道菌群非冗余基因集，同时用了一种新的算法补充基因集，即从来自肠道的单个细菌的全基因组中建立基因集，和来自肠道菌群的基因集进行合并，这个算法可以补充参考基因集覆盖度不全的问题。最终研究者从1,267个肠道菌群宏基因组数据和511个原核生物基因中，构建了IGC基因集，包括9,879,896个基因。用来自中国的丹麦的样品数据比较后，找到了一些地区特异性的基因。
图2 构建非冗余参考基因集的一般流程，图片来自Li Junhua 9。
用来自不同研究团队的宏基因组数据进行一般分析后，得到3GCG基因集。同时从原核生物基因组数据库，挑选来自肠道的原核生物后，构建一般的基因集，同时进行去冗余得到SPGC，最后将3GCG和SPGC进行融合，得到IGC。
非冗余参考基因集有三个问题需要解决，1)组装基因组时，可能会将不同物种的基因组序列拼接在一起得到嵌合体基因。2)基因在去冗余时，使用的其实是uncomplete聚类算法，会将相似的基因聚成一个类，同时挑选最长的基因作为代表基因，而被代表的基因之间，可能相似度很低，也不含有代表基因的功能，但后续分析无法判定这一点。3)Nielsen等10构建的人体肠道的非冗余基因集，从4,201,877个非冗余基因中，判定330,220 (7.8%)的基因来自人，动物，或者植物。这些基因可能来自样品中的食物残渣，但没有数据能明确证实这一点。总之，非冗余参考基因集在构建后需要做质量控制。
3基于非冗余基因菌群的丰度构建
在没有参考序列的情况下，我们仍然有可能构建出未知菌群。前提是我们认为以下条件成立：来自同一个基因组(物种)的任意两个基因, 其丰度相似, 相关性分析在多个样品中也能发现, 这两个基因之间的相关系数很高, 则判定其来自同一个基因组。前提条件的基本原理如图3所示。
图3 未知菌群基因组构建的前提
由于来自同一个物种，基因丰度一定是线性关系，所以使用的是pearson相关系数，目前已经发表的研究中，基于参考基因集有两种构建未知菌群的方法:
1) 基于所有的基因丰度计算出相关系数矩阵，然后将perason相关系数高的基因聚在一起，并为这些基因来自同一个物种。由于基因通常有几百万个，无法计算所有基因的相关系数(如果基因集是4M个，那就需要计算8e12 次，对于一般服务器的计算能力是不可能实现的)，所以聚类时使用的是uncomplete linkageclusting类型的算法，最后聚类得到的菌群称为Metagenomic Species(MGS) 或者Metageomics Clusters (MGC)。Fredirk等对糖尿病的研究6, Nielsen等对未知基因组构建的方法学研究10，使用的均是这种算法，其基本流程如图4所示。
2) 仅限于“疾病－对照” 研究的研究。首先使用假设检验(常用的是Wilcoxon rank sum test)，将在疾病组或者健康组富集的基因挑选出来，然后计算相关系数矩阵。在经过挑选后，得到的差异基因通常只有几万个，对于一般的服务器，计算相关系数矩阵是可行的，然后使用complete linkage clustering 算法，以保证聚类中任意两个基因之间的相关系数都大于阈值。Qin Jujie等的糖尿病研究22，Qin Nan 8的肝硬化研究中使用的均是这种算法。这种聚类方法得到的未知物种称为metagenomic linkage group(MLG)。
图4 MGS聚类算法，基于Nielsen 10的研究整理
聚类基因集合，但定义不同，算法1得到的是所有的物种，算法2实际上得到的是在疾病或者健康组中富集的物种。从数据完整性看，算法1较好，而且对于算法2来说，如果不是“疾病-对照”的研究，无法进行类似的计算。MGS 可以是已知菌群，也可以是未知菌群，取决于聚类内的基因是否能被注释到。
MGS的聚类算法实质上是一个uncomplete mean linkage clustering的算法，1) 从NRG中随机挑选一个基因作为Seed。 2)从NRG中，挑选一个和SEED的PCC大于0.9的基因，加入seed，形成一个canopy。3)同时计算两个基因的profile中位数值，作为canopy的丰度。4)从NRG中挑选和canopy的PCC大于0.9 的基因，将基因并入canopy，同时计算新的canopy丰度。5)迭代Step2和Step3，直到所有的基因都分在不同的canopy中，得到稳定不再增长的canopy,称为 co-abundance gene groups (CAG)。同时对CAG进行质控，要求最终的CAG含有的基因个数必须大于2，同时CAG profile 90%的值,至少在3个样品中出现过(不为0)。超过700个基因的CAG被称为MGS。这个算法将已经聚为一类的基因称为canopy，将所有基因的中位数值作为整个canopy的丰度，所以只需要计算一个循环的PCC，即计算canopy和NRG中所有基因的PCC，而不需要canopy中每个基因和NGR中所有基因的PCC，因此会极大的提高计算速度，但也有可能会损失一些准确度。
无论是MGS还是MLG（以下统称为MGS），两种算法计算过程类似，都是基于基因丰度的得到的。肠道菌群一般只有40~50% 序列11的序列能被判定为来自已知的基因组，极地土壤中一般只有不到1%的序列12能找到已知基因组。MGS的最大意义是，能够最大限度的得到菌群的组成信息，尽管得到的菌群可能是未知物种，Fredirk 10 的研究发现，用MGS数据构建的模型，和仅用已知菌群构建的模型相比，对糖尿病诊断的AUC面积可以从0.71提高到0.83 。MGS的最大缺陷是只能找到基因区，而非编码基因区，非编码RNA区（tRNA, rRNA等），非编码的功能domain（如CRISPR，moble elements等）则无法找到，因此，得到的基因组不是完整的。另外，理论上讲，如果不同的菌株丰度相似（共生菌群），也有可能被聚在一个MGS中，Nielsen等10的测试发现很少存在这个问题，MGS内的基因都来自同一个物种，准确率很高。可能是肠道菌群的优势物种较少，而大部分的物种丰度较低，不同菌群之间的丰度差异较大，因此无法保证其他类型的样品也可以避免这个问题。
4宏基因组靶基因测序分析方法的进展
宏基因组靶基因测序(Amplicion sequencing, AS)指的是，用通用引物将微生物基因组中的保守基因(16S rRNA gene ,ITS gene等)扩增出来，每条序列都来自一个物种，因此可以得到菌群物种的组成。和全基因组测序相比(whole genomeshotgun sequencing, WGS)相比，AS包含的信息量少的多，只有菌群的组成，而且对菌群的分辨率较低（一般只能到属），但和WGS相比，AS测序成本极低，同时数据分析速度很快，使得大规模的测序分析成为可能,。通过上千例样品的靶基因测序数据，在婴儿肠道菌群中发现了和新生儿坏死性结肠炎相关的菌群13，以及顺产婴儿和刨腹产婴儿在全身菌群之间的差异14。在这里我们讨论近年来，针对AS测序数据的问题，数据分析中两个最重要的算法。
4.1 Usearch，序列聚类算法。
Usearch由Edgar RC 15开发并发布。AS测序在完成数据质控16，去除嵌合体17, 序列聚类18后, 每一个聚类单元包含的组成相似的序列集，被称为Operational Taxonomy Unit, OTU。AS测序的数据存在两个测序质量问题：PCR过程中组成相似的序列(来自不同物种)很可能由于PCR产物搭配错误产生嵌合体; 高通量测序仪的数据质量被高估, 混合得到的序列噪音很多，远远超出真实的序列个数21。这些问题导致OTU代表序列很可能并不代表真实的物种，且OTU的个数被严重高估, 远远超出真实的物种。
图5 Usearch 聚类算法，基于Edgar 15整理
Usearch OTU聚类的算法比较简单: 1) 将所有的序列去重复，去除丰度低的序列(singleton)，按照丰度从高到低排序。2)丰度最高的序列一定是OTU代表序列，按照0.97的阈值，依次将序列和OTU代表比较，如果相似度大于阈值，则将序列并入该OTU类中，如果和所有相似度小于0.97 ，则被判定为是一个新的OTU代表序列，若该序列被判为嵌合体则删除。最终的结果中，1)所有的OTU 代表序列两两之间的丰度都小于0.97 的相似度。2) OTU类中的序列和代表序列的相似度都大于0.97。3) 嵌合体必须去除，若序列和多个OTU代表序列相似度都大于0.97，则选择最接近的代表序列。
Usearch有一个很重要的前提假设：和低丰度序列相比，高丰度的序列更有可能是真实的物种，这更符合生物学概念，而不是cd-hit算法19，按照序列长度排序，或者是mothur 20等软件，对所有序列无差别进行聚类。因此，Usearch实际上是一种加权的uncompletedlinkage clustering 聚类算法。也有一些缺陷，由于OTU是按照序列频数排列，高丰度的序列更有可能成为代表序列，导致低丰度的序列很可能会被删掉(很可能是真实的物种)，另外，如果高丰度序列是嵌合体或者其他人工序列，会导致后面一系列的误判，所以仍然需要先去除嵌合体，而不能单纯的依靠自身的去嵌合体算法。总之，这些损失(information loss)和收益相比(information gain)比起来要小得多，Usearch仍然是目前AS数据中最好的OTU聚类算法。
4.2 Picrust，从菌群到功能组成。
AS测序数据通常都是一段保守基因，这些保守基因可以很好的区别不同的物种，因此AS测序数据可以代表菌群的物种组成，但不能直接反映出菌群的功能基因组成。Morgan G. I. Langille等开发的Picrust 23算法, 使用菌群的组成信息，以及公共数据库中已知微生物的功能基因组成，预测出菌群的功能基因组成。
图6 Picrust算法流程，图片来自Morgan 23
Picrust 算法首先构建好了如下数据库:1) 16S 在原核生物基因组中的拷贝数。2) 原核生物基因组的功能基因组成(COG,KEGG)数据库。3) 原核生物基因组的16S序列。用户首先需要将OTU代表序列和16S序列数据库对应起来，其次将OTU序列按照拷贝数进行标准化，因为通常原核生物基因组包含多个16S基因拷贝，最后将同一个功能基因的拷贝数相加起来，则可以得到最终的功能基因组成。需要注意的是Picrust使用限制较多，首先对人体微生态的数据准确度较好，环境样品(空气，水，土壤等)可能并不合适。另外，Picrsut只能用于16S靶基因，而ITS，18S等其他靶基因不能使用。
5机器学习算法在宏基因组临床研究中的应用。
宏基因组数据包括以下几种类型，基因组成，物种组成，功能基因组成。宏基因组数据有两个特点：1)高维变量，基因组成通常有上百万个变量，而物种组成和功能基因组成通常有几千个。2)变量之间并不是独立的。不同的基因来自同一个基因组，他们之间的关系一定是线性的；菌群之间可以是共生或者竞争的关系，也可以是中立的（即完全独立）24; 不同的功能基因会参与同一个代谢路径，他们的关系可以是线性，也可以是非线性，因为有些基因会参与多个功能路径。
临床研究通常有如下几个目的：诊断，预后评估，个性化治疗以及新的治疗方法开发。诊断和预后评估可以被转化为二元分类问题，即病人此时的状态是疾病还是健康，病人预后情况好转还是恶化。而个性化治疗也可以被转化为推送问题，即哪种治疗方案最适合目标疾病的治疗），如果治疗方案为两个，即是二元分类问题，如果治疗方案大于3个，则为多元分类问题。因此，我们在这里主要回顾从2010年至今，宏计量基因临床研究中的分类问题。至于新的治疗方法开发及有效性论证，通常需要做严格的临床实验和因果关系验证，在这里不做讨论。
高维数据的出现，对传统的统计学分析提出了挑战，首先是计算量很大，对算法的速度有了高要求，二是传统的统计学做多元统计分析，经常出现的前提是要求变量之间不相关。而这正是机器学习能解决的，机器学习的一般思路是: 先建立训练样品集(train set)，利用监督学习(如疾病和健康)，对每个变量(特征)进行评估，利用挑选出的特征集构建最优模型，最后用验证样品集(validation set), 评估模型的准确度。一般会将训练集进行随机抽样分为集合1和集合2，用集合1的数据构建模型，并把集合2作为独立的样品对模型进行评估优化。抽样方法通常有两种: Leave one out (LOO)和 n-fold。LOO算法是每次只抽取一个样品，剩下的所有样品都在集合1中用于构建模型，n-fold 是将样品分为n 等份，1份在集合2中，剩下的n-1 份都用于构建模型。假设有m个样品，LOO算法要构建并评估m个模型，而n-fold要进行n次独立的试验。
图7 机器学习的一般流程。
Qin 5在构建肠道菌群基因集时，利用主成分分析(PCA)算法和炎症性肠炎(IBD)和健康人群的肠道菌群组成，进行了简单的区分，结果发现在主成分1上(PC1)，IBD病人和健康人可以区别开(P=0.031)，这表明肠道菌群有作为IBD诊断手段的潜力，但样品个数太少(25个健康，14个IBD)，也缺乏独立样品的验证。PCA算法本质上是一种降维算法25，将输入变量转换为一组独立的转化变量，称为主成分(Principle Component, PC)，如果数据集发生变化，如加入新的样品，或者删除样品，得到的PC也会发生变化。如果用PCA算法构建分类模型，则需要很大的样品个数，构建的PC集较稳定，再加入新的样品时也不会发生形变 (formation)，但目前量化宏基因组研究的样品个数最多只有几百个，而变量个数远远大于样品个数，所以现在的研究水平下，PCA算法不适合用于量化宏基因组研究的模型构建。
在随后的糖尿病研究 22中，研究者在样品个数较大的基础上(345个体)，引入了机器学习中”发现集－构建模型－验证”的概念，即先用322个样品找到了52,484个基因，随后用mRMR算法26挑选出50个marker计算出T2Dindex(用在T2D富集基因的丰度减去健康人富集基因丰度)，T2D越高，表明患病概率越高）, 最终用23个独立样品进行验证。Fredirk 6等对糖尿病的研究，也尝试用肠道菌群对糖尿病进行了诊断, 从800个MGS中，使用随机森体算法(Random Forest, RF)，自适应的挑选了50个MGS构建分类模型，AUC为0.83，更重要的是，使用该模型对49个独立的糖耐受损人群(IGT)进行验证，70%的样品可以被诊断为糖尿病，而IGT被认为是糖尿病发病的前兆之一，充分说明了该模型的可靠性和稳定性。
随机森林被认为是最重要的机器学习算法之一27, RF算法本质上是随机抽取样品和特征，构建大量的分类树，然后按照多数原则进行投票判定样品分类。该算法在发现集样品中的分类效率非常高，但在验证样品中效能下降较快，因此需要的样品个数很多。随后两个独立的结肠癌研究28, 29，都使用了这一算法和肠道菌群的组成构建癌症的诊断模型，表明随机森林算法在宏基因组研究中具有很好的应用场景。
6展望
宏基因组学，包括大规模的全基因组学关联分析(GWAS)，都是随着二代测序的应用(2009年至今)而发展起来的，还是一门年轻的科学。和大基因组(人，动物，植物)相比，宏基因组非常独特，微生物的组成是随时变化的，比如一个正常人，可以在”肥胖”和”正常”两种状态之间转变。宏基因组从稳态到失衡的过程，可能和疾病发生相关的(不是因果)，而GWAS分析无法得到这一信息，这也是宏基因组最有魅力的地方。随着学科的发展，我们对宏基因组也有了更深的认识，但目前仍有很多挑战：
1)随着测序技术的发展，特别是Illumia Hiseq X10测序仪序列在2015年放开了对微生物样品的测序。测序成本的进一步降低，同一批次有可能出现上千人的新的宏基因组测序数据，这对现有的数据分析，计算能力提出了挑战。
2)宏基因组有明显的地域性特点，比如中国和瑞典的糖尿病人肠道菌群组成不同，使用瑞典人的T2D诊断模型6, 22，对中国华南地区的T2D人群诊断完全无效。因此，我们需要收集同一种疾病，在不同地区的样品。这需要多家单位，甚至国际间的科研项目协作，并建立好临床数据和测序数据的共享机制。
3)结论的特异性不好，缺乏实验验证。目前大部分的研究都是简单的“case VS control”描述性研究，结论非常单一甚至一致，比如有益菌的减少，有害菌的增多，而不同疾病的发病机制肯定不同。需要在当前的研究和结论的基础上，构建动物实验模型，进一步研究菌群和疾病进展的相关性。
我们的“肠道微生物群落研究7群“已经有38人了，如果您对肠道微生态感兴趣的话，可以加入我们，加入的方法：请扫描下面的二维码，加我为好友之后，我把你拉到群里。
猜你喜欢
热文：图表规范 DNA提取发Nature
写在后面
为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外六十位多PI，六百多名一线科研人员加入。参与讨论，获得专业指导、问题解答，欢迎分享此文至朋友圈，并扫码加创始人好友带你入群，务必备注“姓名-单位-研究方向-职务”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决推荐生信技能树-微生物组版块(http://www.biotrainee.com/forum-88-1.html) 发贴，并转发链接入群，问题及解答方便检索，造福后人。
学习16S扩增子、宏基因组科研思路和分析技术，快关注“宏基因组”