基于梅花鹿基因组SOLiD测序数据鉴定线粒体基因组的异质性点变异应用

谢留威，李春义，巴恒星※

（1.中国农业科学院特产研究所吉林省特种经济动物分子生物学省部共建国家重点实验室，吉林 长春130112；2.长春科技学院，吉林 长春130600）

通常情况下，每个哺乳动物细胞含有大约100个线粒体，每个线粒体有2~10个线粒体DNA（mtDNA）拷贝[1-3]。在没有组蛋白保护的情况下，mtDNA是高度易变的，其突变率是核DNA的6~17倍[4]。突变在细胞线粒体内累积，与野生型共存，称为mtDNA异质性。目前，mtDNA异质性已成为线粒体遗传疾病的研究热点[5]。许多研究认为，发育、衰老和世代演进的高频率mtDNA异质性是人类遗传疾病的核心[6-7]。

对于已知的mtDNA异质性突变，基于几种PCR方法可以确定样本中的突变位点及其相对频率[8-11]。然而，对于发现未知突变，大规模的高通量测序技术更适合。但最近研究报道显示，mtDNA异质性检测正从扫描已知的有限数量突变转变为全线粒体基因组筛查[12-13]。

ABI SOLiD测序平台使用短读长（35~75 bp），通过一种颜色编码框架（color space），测序不同碱基之间颜色信号。每一个测序碱基都依赖于2个连续的颜色编码值，当短读长与参考基因组比对后，如出现测序错误（即单色差异）就可以与正确的碱基变化区别分开，这为重测序研究提供了新的非常具有优势的判定方法。这种独特的编码方式具有前后碱基校对功能，使碱基测序准确率高达99.94%[14]。本研究通过对一只雄性梅花鹿血细胞全基因组SOLiD测序数据进行分析，检测mtDNA基因组中异质性突变，以期为充分利用SOLiD全基因组测序数据鉴定未知的mtDNA异质性点变异提供基本生物信息分析框架。

1 材料与方法

1.1 数据来源

利用ABI SOLiD测序平台对一只雄性梅花鹿血细胞进行全基因组双末端测序，读段长度为50 bp。另外，从GenBank中下载与全基因组测序梅花鹿属于同一亚种的mtDNA全长序列（登录号：KX689229）作为SOLiD测序数据比对的线粒体参考序列。

1.2 分析方法

低频变异检测需要高质量的测序数据，过滤掉包含未知碱基或平均质量分数低于20的读段。使用PerM v0.4软件[15]，选项设置为“-A-E-seed F4-V5-L400-e”，将高质量的读段与mtDNA参考序列进行比对。在许多真核生物中，mtDNA经常被转移到核基因组中，产生了核线粒体假基因序列（nuclear mitochondrial DNA segments，NUMTs）。不同大小NUMTs片段与mtDNA序列存在高度差异[16-17]。根据NUMTs高变异率的特点，排除包含大于5个不匹配颜色值，相当于2~3个不匹配碱基（相似度≥94%）的读段。由于NUMTs平均长度为240 bp[18]，SOLiD成对末端测序文库的平均双末端长度为831 bp（图1），这可以在一定程度上减少NUMTs的影响，同时去掉单端映射读段。最后，利用Samtools v1.2[19]中的mpileup工具和Bcftools v1.3.1软件[20]联合检测mtDNA异质性点变异，在Tablet v1.18软件中对点变异可视化显示。

图1双末端长度分布Fig.1 Insert size distribution

2 结果与分析

2.1 比对结果

过滤低质量数据后，共获得410 662 782（38 Gb）成对短读段，相当于全基因组测序深度约13（38 Gb/3 Gb）。大约0.035%的成对读段比对到mtDNA参考序列，覆盖率为99.62%，在2 265~2 329（65 bp）区域未被短读段覆盖（图2）。线粒体基因组的平均测序深度约为892，一方面表明血细胞中mtDNA的拷贝数较高，另一方面表明部分源于NUMTs的读段被比对。然而，动物中NUMTs的总量约为核基因组的0.1%。在梅花鹿基因组中，NUMTs的比例（0.035%）相对较低，表明高度可变的NUMTs读段被排除掉，这进一步提高了鉴定mtDNA异质性点突变的准确性。

图2 mtDNA上2 265~2 329（65 bp）区域未被短读段覆盖Fig.2 mtDNA region 2 265-2 329（65 bp）uncovered by short reads

2.2 mtDNA异质性点变异鉴定

本研究共检测到8个点突变（图3），包括4种转换和4种颠换（表1）。其中，5个点突变位于蛋白质编码基因（COX1、COX2、ND4和ND5），2个位于tRNASer，1个存在于12SrRNA中。为了验证NUMTs读段是否影响点突变的鉴定，通过检索GenBank中18个梅花鹿线粒体全基因组中相应位点的变化，证实这8个点突变都存在于18个mtDNA基因组中（表1），表明它们是mtDNA特有的单核苷酸多态位点，而不是由NUMTs引入的。重要的是，C6180T、T7481A和A10909T这3个位点分别导致氨基酸的变化，即：丝氨酸（Ser）脯氨酸（Phe）→亮氨酸（Leu）→脯氨酸（Phe）、亮氨酸（Leu）→谷氨酰胺（Gln）。mtDNA异质性点突变的选择性压力与组织特异性代谢率、细胞周期和生物能量需求有关[18]，这暗示梅花鹿不同类型的血细胞可能承受不同的选择压力。

图3 8个mtDNA异质性点突变可视化比对Fig.3 Visual comparison of 8 mt DNA heteroplasmic point variations

3 结论

目前，尽管SOLiD测序平台已不常用，但其已产生大量原始基因组测序数据，这些数据仍然具有较高的科学价值，特别是SOLiD平台对测序碱基进行颜色编码校对特性，在利用重测序进行点变异检测方面具有固有优势。本研究通过生物信息学分析方法充分挖掘一只雄性梅花鹿血细胞全基因组SOLiD测序中短读长数据，在mtDNA基因组中鉴定了8个异质性点突变，将为利用全基因组测序数据鉴定未知的mtDNA异质性点变异提供一个基本分析框架。