2014 JPR: 逐次累进法修正非模式菌株基因组
Iterative Genome Correction Largely Improves Proteomic Analysis of Nonmodel Organisms. |
细菌是简单而又复杂的生物，虽然基因组小，但变异相当多。例如，革兰氏阳性菌的同种菌不同菌株之间，基因组差异可达5~20%，其变异部分与其致病性、抗生素耐药性十分相关。这种较大的差异使得基因组、转录组、蛋白质组的研究十分困难。
现在的大规模测序技术可以轻易将细菌基因组打成碎片，测序得到数百万条短核苷酸序列(reads)。利用这些信息，传统上人们采用两大类方法来解决这一问题：基于mapping的基因组修正策略（Mapping-based genome correction）和 基因组从头组装策略（De novo genome assembly）。基于mapping的基因组修正策略是将reads比对到同种标准菌株已测定的基因组上，搜索单核苷酸突变(SNV)，将基因组修正过来，这种策略的好处是比较精确而且可以利用已有的精确的基因功能注释信息；而从头组装策略则是直接根据测得的reads进行组装，然后进行基因判别与自动注释。但两种方法都有非常明显的缺点。
基于mapping的基因组修正传统算法，修正效率对不同数据集波动极大，且可低至0.03%（基本修正不了）。传统较为精确的算法消耗极大的计算资源，修正一个大肠杆菌的基因组也要运算10个小时以上，消耗10G以上的内存空间。 即便如此，也最多只能修正1.6%的基因组差异，无法应对5~20%的差异。
而基因组从头组装则费时费力费钱，而且组装效果往往达不到预期，主要表现为组装错误率高。有例子显示，Alcanivorax borkumensis SK2 基因组组装结果，有高达5%错误率，即每20个碱基就有一个错。因此，需大量验证和Sanger测序校正工作。即便组装100%正确，接下来的难题是基因功能注释非常困难。用目前最好的基因自动判别和注释算法对 Alcanivorax borkumensis SK2 基因组自动注释，ORF注释假阴性率49%，假阳性率52.8%，几乎没法用。
为了解决这一问题，我们利用FANSe系列算法超高的准确性和容错性，研制了逐次累进法基因组修正算法。其原理如下图：
 |
图1： 将测序的到的reads向标准菌株基因组上mapping，寻找SNV，并修正基因组序列。这一过程反复多次，每次都能修正多一些的突变，直到将所有突变基本上都修正完毕，生成真实菌株的基因组。然后利用这一真实基因组来生成参考蛋白质组序列，用以指导蛋白质组质谱鉴定。 |
我们从环境中分离了一株菌株，16S rDNA序列测序证明其与 Bacillus pumilus SAFR-032 相同。qPCR法测定其基因组大小，3.41 Mb +/- 0.13 Mb，比Bacillus pumilus SAFR-032模式菌株的3.7Mb小，估计无大片段插入。对其进行全基因组大规模测序，使用 Illumina HiSeq-2000 pair-end 2x100nt 模式测序，获得24M reads，目标平均测序深度400~500x。
我们的方法进行7轮修正，mapping率迅速提高，稳定在73.2%，共修正182620个SNV，基因组差异大约5%。而经典的传统方法(Bowtie2+Samtools)只修正了169个SNV便无法继续修正了。差别竟达到了千倍以上！
用Sanger测序法随机选取4个基因区段，共验证测序1994个位点，结果发现我们的方法无一假阳性、无一假阴性，而作为对比的Bowtie2和Stampy都有许多假阳性和假阴性。我们还发现，突变率高达19.2%的区段，我们的方法也能完美应对。
图2：BPUM_1139基因的一个区段，对比Bowtie2, Stampy, FANSe三种算法mapping之后的修正结果，并与毛细管Sanger测序法测得的真实基因组进行比对。灰色的是与真实基因组不符的碱基。全部1994个位点的测序验证图请看原文的supplementary materials.
图3：总共用Sanger测序法测定4个基因区段，共1994个位点，比较三种方法修正的假阳性假阴性。
使用修正过的基因组产生参考蛋白质组序列，使全菌蛋白质质谱鉴定在蛋白质水平提高9%，肽段水平提高13%，且鉴定质量提高很多。用传统方法修正的参考蛋白质组序列甚至损害了质谱鉴定效率！少鉴定了8个肽段！
为了搞清楚为什么传统方法竟然会损害了质谱鉴定效率，我们仔细检查了8个肽段。以rpsC的一个肽段为例。传统方法(Bowtie2+Samtools) 认为这里有一个碱基插入，造成移框突变和截短，而 我们的方法认为这里没有突变。实际Sanger测序证明这里没有突变，质谱甚至直接鉴定到了覆盖这一位的肽段，也证明这里没有突变，传统方法鉴定错误。从道理上来说，rpsC是核糖体蛋白，截短将造成严重的核糖体功能缺陷，而我们所分离出的这株菌生长旺盛，不像是核糖体功能缺陷。对覆盖此位点的787个reads进行人工检视，确认无突变，仅有少量reads存在测序误差，确认传统算法修正错误。类似的情况还发生在rpoS, fbaA等基因上。
 |
图4： 基因rpsC的修正结果与质谱鉴定结果。绿色为蛋白质质谱鉴定到的肽段，绿色数字为Mascot质谱鉴定软件给出的肽段分数，>20为可靠鉴定。BS为传统方法修正 (Bowtie2+Samtools)，FF为我们的方法修正。 传统方法认为113942位有一个碱基的插入，导致移框突变，因而认为rpsC蛋白出现移框突变，将生成一个87氨基酸的蛋白质。而实际上质谱鉴定到了跨越113942位的肽段。 |
 |
图5： 对覆盖此位点的787个reads进行人工检视。蓝色为与原参考基因组匹配的碱基，红色为与原参考基因组不匹配的碱基。可见绝大多数碱基都能匹配上，此位无突变。 |
由此，我们可以看到，逐次累进修正基因组方法为非模式物种奠定研究基础。
- 逐次累进修正基因组方法可应对高达20%的差异
- 修正精度高，效果好，符合实验验证
- 计算快速，普通台式机甚至笔记本就能很快搞定，一轮修正耗时仅1-2小时
- 为非模式物种的基因组、转录组、蛋白质组研究提供高精度的参考序列
目前，这一算法已可在承启生物基因测序云平台上提供公众服务：
