16SrRNA基因可变区与全长序列进化关系相似性分析.docx

? ? 16SrRNA基因可变区与全长序列进化关系相似性分析 ? ? 刘爽爽， 帖 云， 齐 林， 刘峰辉， 王 磊 (1. 郑州大学 信息工程学院 河南 郑州 450001； 2. 郑州大学第一附属医院 河南 郑州 450052； 3. 河南省人民医院 口腔医学中心 河南 郑州 450003) 0 引言 16S rRNA基因是微生物生态学分析中最常使用的一类分子标志物，存在于所有细菌的基因组中，具有高度的特异性和保守性，序列长度约为1 500 bp[1]。除保守区外，16S rRNA基因序列还存在V1～V9共9个可变区[2]，不同可变区的长度范围为100～300 bp，新一代测序技术可以使用短配对碱基轻易覆盖，使得 16S rRNA 序列可变区的测量更加便捷。 可变区的特异性能够反映出不同微生物的特征核苷酸序列，用于分析复杂生物环境中微生物的物种多样性[3]、相对丰度[4]、物种鉴定及进化距离等[5]。文献[6]对16S rRNA 基因两个可变区进行了比较研究，结果显示，在肠道菌群物种多样性分析及物种鉴定能力上，选择V1～V3可变区片段进行测序，得到了与全长序列更为接近的结果。文献[7]利用16S rRNA全长序列和部分基因作为热测序的靶点，在不同水平上分析了16S rRNA基因在基因组内因异质性引起的高估问题，结果表明，对于细菌使用针对V4和V5可变区的引物可以将这种高估最小化。文献[8]采用Illumina Miseq测序技术，测定了苏尼特和阿拉善双峰驼的自然发酵驼乳中微生物16S rRNA的V3、V4可变区序列，并对群落结构和物种多样性进行了比较分析。文献[9]对HIV-1包膜蛋白gp120进行分析时，找到了可变区V1可能作为传播选择靶位点的证据。文献[10]分析了大西洋鲑鱼细菌16s rRNA基因全长序列及不同可变区对微生物群落结构的影响，发现不同可变区对微生物分布和系统发育有着不同的影响。目前可变区在物种进化关系中表现如何的研究较少，但其对物种进化来源分析具有重要的指导意义。本文以核糖体数据库项目(RDP)所提供的细菌16S rRNA基因数据为基础，构建不同可变区及全长序列进化树，使用层次距离矩阵算法分析了V2、V3、V4可变区与全长序列所构建的进化树之间的距离差异值，并对可变区与全长序列进化关系的相似性进行了分析。 1 数据获取与预处理 1.1 可变区截取与筛选 原始数据采用RDP中细菌16S rRNA的全部序列。压缩文件大小为3 GB,解压后大小为76 GB，共包含约320万条16S rRNA序列，数据格式为fasta。由于V9区在实际研究中应用较少，故只选用V1～V8可变区进行相关研究。使用MEGA6软件在该序列中分别寻找各可变区两侧保守序列，保守序列及可变区位置如表1所示。 表1 可变区两侧保守序列及位置 确定了可变区位置后，使用biopython函数库(https:∥/)中的序列切片方法进行可变区片段的截取。初步截取后的序列中仍存在一些含有实际碱基数目较少且信息量较低的序列，因此需要分析可变区片段实际碱基长度并筛选出含有一定信息量的序列。使用biopython库中的seq.parse函数读取初步处理后的序列，统计每个序列中的实际碱基数目，并使用matplotlib库(https:∥/)绘制可变区各序列的实际碱基数目分析图，以便对序列进行初步筛选。对序列的初步处理操作由python脚本完成，所使用的核心函数库为biopython库。 8个可变区使用相同的方法进行实际碱基数目统计，以V2可变区为例，实际碱基数目结果如图1所示。可以看出，可变区片段中实际碱基数目出现了明显的拐点，有约70万个V2可变区片段中碱基缺失较为严重。将拐点处放大，可以观察到部分V2可变区片段实际碱基数目在80以下，表明在这些序列的测序过程中，V2可变区的测序出现了遗漏或者并未对V2可变区进行测序。因此，需要按照拐点处实际碱基数目对约300万个可变区片段进行筛选，以保留含有一定信息量的可变区片段。 图1 V2可变区实际碱基数目 筛选操作仍由python脚本完成，在完成了8个可变区片段的初步提取和筛选后，就得到了可使用命令行工具进行处理的数据。对提取出的可变区片段进行去冗余与去除嵌合体操作，以V2可变区为例，两端对齐可以发现，序列数目为2 487 500，片段长度为1 264 bp，实际碱基数目为80～120，分别使用unique.seqs与chimera.uchime函数去除冗余部分序列和包含嵌合体较多的序列，再次对序列进行总结分析，此时序列数目缩减为533 136，片段长度缩减为771 bp。在完成了所有可变区的筛选、过滤操作后，绘制8个可变区预处理前后序列碱基数目对比图，结果如图2所示。通过可变区截取、筛选等数据预处理后，从序