中国科学报
2020-02-08 08:27:35
作者 | 薛宇（华中科技大学生命科学与技术学院教授）1月21日下午在实验室忙完工作回家，晚上一起看电视，看见钟南山院士讲新冠肺炎有人传人。第二天22日起床，专门看了一下新闻，没有特别的内容，只是有专家建议“武汉人最好别出去，外地人最好别来武汉”。还是没有反应过来是什么情况，所以按原计划自驾回安徽。
高速上的交通还好，来来往往的车辆很多，不过我们既然是从武汉出来的，那进服务区下车都很注意戴上口罩。到了滁州住下，第二天起床发现武汉封城了，商量了半天是不是该立马开回去，后来决定既来之则安之。在滁州住了两天之后，于24日下午开回了合肥。回去了之后立马向社区报备，然后就呆家里不出门了。在合肥的几天过得蛮好，该吃吃该睡睡。后来继续讨论，究竟是留在合肥还是回武汉，最后结论是回去。
无他，在武汉工作生活了这么多年，太多的朋友、同事都在武汉，我得回来。所以30号回了武汉，一切安好，体温正常。讲完废话我们讲正事儿。前几天写了篇博文《新冠病毒不是源自实验室》，写完之后被各路朋友吐槽到体无完肤。最最严厉的批评就是：说好的通俗易懂、诙谐幽默哪儿去了？于是，有了下面的科普版本。
上一篇博文里，我们用到了“简单的生物信息学分析”，其中涉及到两个算法包括用于全局双序列比对的Needleman-Wunsch算法和用于局部双序列比对的Smith-Waterman算法，以及常用的序列比对工具BLAST，这些内容是几乎所有基础生信教科书里必须要讲解的、最基础的内容。所以，“简单”二字，这是名副其实，经得起历史检验的。其次，即使是教科书上的内容，那也有一定的专业性。
这里，咱引用轩哥的总结做一个较为简单的科普：1）老美（James Lyons-Weiler博士）发现的新冠病毒一段序列（1378个碱基），也出现在其他冠状病毒的序列中，且序列高度相似。说明这一段序列在自然界的冠状病毒中广泛存在，并非来自实验室。2）pShuttle-SN是为了研究非典SARS病毒，建立的表达载体，其中包含一段来自于SARS病毒Spike基因的序列。科普完毕。
各位如果不爱看专业分析，请忽略本文所有非加粗的文字，以及下文。下面继续我们的“简单的生物信息学分析”，因为上文有一个很重要的问题没有解答：James从哪儿整出来的那段1378个碱基的片段？从理论上来讲，咱和James用的都是BLAST，搜索的数据库都是GenBank，为啥出来的结果不一样？那究竟是大菜鸟算错了，还是James算错了？这是一道送分题，答案必须是James算错了。错在什么地方？
各位请看下图：NCBI的BLAST（https://blast.ncbi.nlm.nih.gov/Blast.cgi），在程序选择上有三个选项。第一个megablast运算速度快，但是准确性不高；第二个非连续的megablast运算速度要慢一些，但是准确度要高一些；第三个blastn是BLAST程序包里最经典的一个工具，速度慢但是准确度最高。
BLAST的默认选项是megablast，James应该是没有调参数直接贴新冠病毒的基因组序列（https://www.ncbi.nlm.nih.gov/nuccore/NC_045512.2）进去做的比对，见下图：新冠病毒的序列，实际上跟自然界其他病毒的相似度，比跟pShuttle-SN里包含的冠状病毒序列更高。见证奇迹的时候到了！右边那个就是James找到的INS1378。
但问题是左边看起来新冠还有一个插入，为啥James就不讲呢？究竟是为啥呢？答：应该是视力不佳没有看见。三种不同的BLAST工具，算法原理都是先匹配个种子序列，然后两端延伸，megablast使用的种子序列最长，blastn用的最短。这样在序列比对中，如果一段序列两端跟其他序列的相似性不那么高，就可能匹配不上种子序列，这样整段序列就被忽略掉，看起来像是个插入。
无巧不成书，新冠病毒这两段看起来像是专属的插入，正好符合这样这个情况。这两段插入，两端的序列跟其他病毒的相似性比较低，所以megablast就检测不到了。但是如果我们换成第二个非连续的megablast，如下图：可以看见，上图那两段较长的插入片段就没有了。程序调成blastn也是类似的结果，这里我们就不贴图了。
这样我们就很清楚，James是公司的技术人员，不是咱生信领域的专家，对于“简单的生物信息学分析”不熟悉，不懂得BLAST的计算原理，不清楚不同BLAST程序的细微区别，所以犯了一个小小的、微不足道的技术性错误，导致他基于INS1378做得所有推测都不成立。
James犯得第二个小小的错误，是当找到一段疑似的、可能仅在某个基因组中存在的插入片段时，理论上生信学者的第一反应，是把这段序列重新比对回数据库（也就是回帖），看看是不是的确不能匹配到其他基因组上，而不是瞎找一条载体序列来做比对。
咱在很多年前写过一篇博文《如何成为顶级生物信息学家》，按研究水平将生信学者分成五个等级，其中1级（Level 1）是给数据、能分析；2级（Level 2）是想新招、玩数据。James不是专业的科研人员，但做分析的态度严谨，推导结论也比较小心，就生信的水平而言，接近1级水准。
咱从02年开始搞生信，06年底博士毕业的时候就是妥妥的2级水平，所以搞清楚一位准1级的生信爱好者究竟犯了什么技术性错误，不是特别困难的事情。最后，我们在讲讲上贴的评论里，有朋友提到的两个主要问题：1）“关于突变速率的推断应该是基于自然条件下，如果实验室培养这个速率应该更快”。看上一篇博文的朋友请注意，利用不同的计算模型，估算出来的病毒突变速率不太相同，都会有一个范围。
就目前的测序数据来看，新冠病毒基因组的突变特别少，变异速率要比如SARS-CoV的速率低。我们写博客不能搞得过于复杂，简单粗暴用个理论最大值算一下即可。2）朋友发来的“某人看过你的分析后评论：这个分析的第二条完蛋了，薛宇说两个病毒分叉是6.6年，可是RaTG13貌似正好是2013年采集的病毒”。根据这个评论，我们可以建立的理论模型是RaTG13有可能在采集回来之后即发生泄漏。
基于这个模型，我们理论上应该能够从野外找到大量的介于新冠和RaTG13的中间过渡型，并且能够建立起明确的流行病学关系，这两个条件必须同时满足。咱上面已经讲过了，目前已经测序的新冠基因组，突变特别少，所以现有数据不支持这个猜想。发现新冠病毒的多样性和溯源，这个不是咱生信的菜，请相关专家自行解决。
哦，还有一个最后，James的第二篇博文是2月2日贴的，轩哥说有美国的朋友打电话来，说James的猜想在美国闹翻天了。这该如何是好？简单，美国懂得“简单的生物信息学分析”的学者比国内多，自己想办法发帖子讲清楚就可以了。编辑 | 宗华 赵路
排版 | 李言
请按下方二维码3秒识别
阅读原文