一种宏基因组16S rRNA的高通量测序数据处理及分析流程控制方法
【技术领域】
[0001] 本发明设及药物基因组学和计算生物学领域，具体设及一种宏基因组16SrRNA的 高通量测序数据处理及分析流程控制方法。
【背景技术】
[0002] 宏基因组学在微生物研究中占据了非常重要的地位，宏基因组是W环境中微生物 的基因组的总和为研究对象。16SrRNA(smallsubunitribosomalRNA)基因是对原核微 生物进行系统化分类研究时最常用的分子标志物，广泛用于微生物生态学研究中。近年来 随着高通量测序技术及数据分析方法等不断进步，大量基于16SrRNA基因的研究促进了微 生物生态学的快速发展，例如：气候变化、水处理工程系统、大气污染、极端环境、人体肠道、 石油污染修复和生物冶金，甚至和人体健康也密切关联。然而使用16SrRNA作为分子标志 物时也存在诸多问题，例如水平基因转移、多拷贝的异质性、基因扩增效率的差异、数据分 析方法的选择等，运些问题影响了微生物群落组成和多样性分析时的准确性，尤其是与高 通量测序技术相关的大数据处理及分析流程控制，给相关科研工作者带来了挑战和困难， 成为该领域目前急需解决的问题。
【发明内容】
[0003] 本发明要解决的技术问题是克服现有技术宏基因组16SrRNA高通量测序数据处 理中不准确性、W及分析流程中步骤繁琐、费时费力等缺陷，提供一种16SrRNA的高通量测 序数据处理及分析流程控制方法。
[0004] 为解决上述技术方案，本发明提供一种宏基因组16SrRNA的高通量测序数据处理 及分析流程控制方法，其特征在于，其包括如下步骤： 阳0化](1)自定义参数配置文件的生成步骤；导入宏基因组16SrRNA高通量测序原始序 列数据，经过筛选和拼接得到理论上有效的16SrRNA可变区全长序列，在此基础上进行生 物信息学参数分析；
[0006] (2)输入步骤：用户根据需要，输入设定的各参数配置文件；
[0007] (3)分析步骤：根据参数配置文件，宏基因组高通量数据处理流程模块生成对应 的自动化分析流程；
[0008] (4)执行及输出步骤：执行所描述的自动化分析流程，获得并输出宏基因组16S rRNA分析结果报告。
[0009] 本发明的优选技术方案中，所述的步骤（1)中，具体包括如下步骤：
[0010] (A)导入宏基因组16SrRNA高通量测序原始序列文件， W11] 做对所述的宏基因组16SrRNA高通量测序原始序列文件进行质量控制与统计， 并剔除低质量序列数据，获得经过筛选的序列数据；
[0012] (C)将所述的经过筛选的数据进行拼接，组装成全长的16SrRNA可变区序列；
[0013] (D)将拼接结果进行质量控制，并去除嵌合体，得到理论上有效16SrRNA的全长 序列。
[0014] 本发明的优选技术方案中，所述的步骤（C)中，使用PANDseq拼接软件，对重叠区 域进行比对打分，比对打分值低于0. 6时将被去除，重叠区域小于5bp或者重叠区域大于2 个mismatch也就去除，根据拼接结果选择有效序列在400~480bp之间的序列用于下一步 分析
[0015] 本发明的优选技术方案中，所述的步骤值）中，先UCHIME软件在de-novo模式下 去除嵌合体序列，然后USEARCH软件在有参模式进一步去除嵌合体序列，最终得到理论上 有效的16S rRNA可变区全长序列。
[0016] 本发明的优选技术方案中，所述的步骤（1)中，生物信息学参数分析包括对于获 得的16SrRNA可变区全长序列进行聚类；包括输入指令采用使用UCLUST方法进行0TU聚 类，0TU中序列相似性设为97%，得到0TU列表及0TU代表性序列。
[0017] 本发明的优选技术方案中，所述的步骤（1)中，包括进一步对0TU代表性序列进行 物种分类分析。所述的物种分类分析包括，物种进化分析，物种丰富度分析，物种鉴定分析 和α多样性指数分析。
[0018] 本发明的优选技术方案中，系统将多样品0TU代表性序列进行聚类与差异性分 析，包括β多样性分析和多样品聚类分析。
[0019] 对每个0TU选择一条代表性序列，使用畑Ρclassifier对代表性序列进行物种分 类注释，从而得到每个样本的群落组成。
[0020] 在本发明的一个实施方案中，使用畑P classifier贝叶斯算法对97%相似水平 的0TU代表序列进行分类学分析，并在各个水平统计每个样本的群落组成，比对数据库为 Silva_11116S rRNA database化ttp://www. arb-silva. de/)。
[0021] 本发明的方法还可W对多个样品进行样品聚类分析，如采用Qiime平台，使用 UPGMA(Unweightedpairgroupmethodwitharithmeticmean)聚类方法，基于weighted uni化ac和unwei曲teduni化ac距离矩阵，将样品进行聚类。
[0022] β多样性值为两个样本间的相异系数，反映不同样本间的多样性的差异，利用各 样品序列间的进化和丰度信息计算样品间的距离，反映样品间是否有显著地微生物群落差 异。在本发明的一个实施方案中，采用Qiime平台，首先利用来自不同环境样品的0TU代表 序列构建一个进化树，化i化ac度量标准根据构建的进化树枝的长度计量两个不同环境样 品之间的差异。化i化ac分析分为wei曲teduni化ac和unwei曲teduni化ac两种度量方 法，两者之间差异在于是否计入不同环境样品的序列相对丰度。wei曲teduni化ac算法在 计算树枝长度时将序列的丰度信息进行加权计算，因此unwei曲teduni化ac可W检测样品 间变化的存在，而wei曲teduni化ac可W更进一步定量的检测样品间不同谱系上发生的变 异。
[0023] 在本发明的方法中，使用Qiime平台，采用对序列进行随机抽样的方法，W抽到的 有效序列数进行0TU的分析，并分别分别使用ACE算法、化ao算法、Shannon算法、Simpson 算法、Good's Coverage计算各α多样性指数。
[0024] Ace:用来估计群落中含有0TU数目的指数，由化ao提出，是生态学中估计物种总 数常用指数之一。（ht1:p: //www.mothur.org/wdki/Ace) 阳ο巧]
[0026] Πι:表示含有i条序列的OTU数目；
[0027] 油unf:设定的一个0TU丰度阔值；
[0028] Srare:低于或等于该丰度阔值的0TU数目；
[0029] Sgbu"d:高于该丰度阔值的0TU数目；
[0030] 化ao:是用化aol算法估计样品中所含0TU数目的指数，化ao在生态学中常用来 评估物种总数。（ht1:p://www.mothur.org/wdki/Qiao)
[0031]
阳03引 Schaol:最终评估的0TU数目； 阳03引 S"bs:实际测出的0TU数目； W34]Πι:表示含有1条序列的0TU数目； W35] Π2:表示含有2条序列的0TU数目；
[0036] 化annon:常用于反映α多样性指数，用来估算样品中微生物多样性。化annon值 越大，说明群落多样性越高。（ht1:p://www.mothur.org/wdki/Siannon)
[0037]
阳03引 S"bs:实际测出的0TU数目； W39]Πι:表示含有i条序列的0TU数目； W40]N:所有测得序列数。
[0041] Simpson:辛普森多样性指数，由EdwardHu曲Simpson(1949)提出，在生态学 中常用来定量的描述一个区域的生物多样性。Simpson指数越大，说明群落多样性越低。 (http://www.mothur.org/wiki/Simpson)
[0042]
阳〇创 S"bs:实际测出的0TU数目；
[0044]Πι:表示含有i条序列的0TU数目； W45]N:所有测得序列数。
[0046] Good'sCoverage:是指各样本文库的覆盖率，其数值越高，则样本中序列没有被 测出的概率越低。（ht1:p: //www.mothur.org/wdki/Coverage)
[0047]
W48]Πι:表示含有1条序列的OTU数目； W例 N:所有测得序列数。
[005