讲正事儿之前，我们先按惯例闲话几句。“时近年节返故乡，过汉江，出黄冈。暂离江城，高堂迎道旁。黍米鸡粟忙遍尝，举杯酒，庆吉祥！”这是某粗通“简单的生物信息学分析”的科研工作者在春节将近，初步拟好的准备发给亲朋好友的拜年词。话说咱逢年过节给人发贺辞基本不用现成的，都是临时瞎编个打油诗，心意到了就好。豆儿妈在边上看见了，说，你难道是想告诉全世界你从武汉逃回了合肥？说得好有道理，赶紧去翻诗词大全，把刚才那个打油诗改成：“时值旧疫袭鄂乡，过汉江，出黄冈。佳节安度，一切依平常。黍米鸡粟忙遍尝，举杯酒，庆吉祥！”这样看上去是不是就不那么狼狈了？至于水平怎么样，这个，咳咳，咱是做“简单的生物信息学分析”的对吧？能糊弄个打油诗就可以了。
这样咱就讲讲这段时间的经历。往年我们家一般年三十或者二十九回合肥过年，去年一整年忙得跟打仗一样，所以早早就定了年后出去旅游的计划，因此决定今年早些回家乡。1月21日下午在实验室忙完工作回家，晚上一起看电视，看见钟南山院士讲有人传人。咱神经反射弧比较长，跟豆儿妈研究了一晚上没琢磨清楚究竟是怎么回事，因为即使人传人也还可能是有限人传人。第二天22日起床，专门看了一下新闻，没有特别的内容，只是有专家建议“武汉人最好别出去，外地人最好别来武汉”。还是没有反应过来是什么情况，所以按原计划自驾回安徽。高速上的交通还好，来来往往的车辆很多，不过我们既然是从武汉出来的，那进服务区下车都很注意戴上口罩。到了滁州住下，第二天起床发现武汉封城了，商量了半天是不是该立马开回去，后来决定既来之则安之。在滁州住了两天之后，于24日下午开回了合肥。回去了之后立马向社区报备，然后就呆家里不出门了。
老憋在家里人也着急，上蹿下跳跟猴儿似的。我妈倒是挺开心，问为啥啊？答：你从上大学开始到今天，这是头一回过年在家里呆这么长时间的。想想，好像是这么回事。这么多年晃晃悠悠，大多数都是瞎忙。不过即使是瞎忙，咱也尽全力瞎忙，争取能对得起纳税人给咱发的这份工资就好。在合肥的几天过得蛮好，该吃吃该睡睡，突然发现自己居然吃好喝好睡好就能为社会做贡献，这个大概一生也难得遇到几次。后来继续讨论，究竟是留在合肥还是回武汉，最后结论是回去。无他，在武汉工作生活了这么多年，太多的朋友、同事都在武汉，我得回来。所以30号回了武汉，一切安好，体温正常。
=======================================================
讲完废话我们讲正事儿。话说昨天咱写了篇博文《新冠病毒不是源自实验室》，写完之后被各路朋友吐槽到体无完肤，最最严厉的批评就是：说好的通俗易懂、诙谐幽默哪儿去了？这个批评的确是很严重的问题，所以我要甩锅。甩锅自古以来就是个技术活儿，遇到问责，先说自己没错，这是错误的反应。正确的该怎么说？得这么说：各位说得对啊，我上一篇写得确实不够通俗易懂，确实没有做到诙谐幽默，我有问题。我在科普性和专业性的平衡方面考虑不够，值得反省和反思。疫情当前，科学家们要主动搞好科普，例如我写博文用到的“简单的生物信息学分析”，在浙大立铭教授的博文中有较好的描述，请各位多关注。是吧，这就把锅甩出去了。当然我这段甩锅的文字纯属拍脑袋瞎编，如有雷同，概不负责。
上一篇博文里，我们用到了立铭教授吐槽过的“简单的生物信息学分析”，其中涉及到两个算法包括用于全局双序列比对的Needleman-Wunsch算法和用于局部双序列比对的Smith-Waterman算法，以及常用的序列比对工具BLAST，这些内容是几乎所有基础生信教科书里必须要讲解的、最基础的内容。所以，“简单”二字，这是名副其实，经得起历史检验的。其次，即使是教科书上的内容，那也有一定的专业性。这里，咱引用轩哥的总结做一个较为简单的科普：
=======================================================
1）老美（James Lyons-Weiler博士）发现的新冠病毒一段序列（1378个碱基），也出现在其他冠状病毒的序列中，且序列高度相似。说明这一段序列在自然界的冠状病毒中广泛存在，并非来自实验室。
2）pShuttle-SN是2005为了研究非典SARS病毒，建立的表达载体，其中包含一段来自于SARS病毒Spike基因的序列。新冠病毒的序列，实际上跟自然界其他病毒的相似度，比跟 pShuttle-SN里包含的冠状病毒序列更高。
=======================================================
科普完毕。各位如果不爱看专业分析，请忽略本文所有非加粗的文字，以及下文。下面继续我们的“简单的生物信息学分析”，因为上文有一个很重要的问题没有解答：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）进去做的比对，见下图：
见证奇迹的时候到了！右边那个就是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有可能在采集回来之后即发生泄漏。基于这个模型，我们理论上应该能够从野外找到大量的介于新冠和RaTG 13的中间过渡型，并且能够建立起明确的流行病学关系，这两个条件必须同时满足。咱上面已经讲过了，目前已经测序的新冠基因组，突变特别少，所以现有数据不支持这个猜想。发现新冠病毒的多样性和溯源，这个不是咱生信的菜，请相关专家自行解决。
最最后，话说我离开合肥回武汉的时候，家父向所在社区的疫情指挥部汇报我们行程，问：为啥回去了呢？家父说娃儿回去搞科研了，社区肃然起敬。我说，爸爸您这不是坑儿子吗？您都知道我回去就是吃好睡好玩好！昨天贴了博文，打电话给家父，说我算是搞科研了啊，以后再有人问，那就不是忽悠了对吧。
还有个最后。上篇和本篇的写作中，得到了战忽局湖北站各位同志们的大力支持、技术指导和研讨，特此致谢。
哦，还有一个最后，James的第二篇博文是2月2日贴的，轩哥说有美国的朋友打电话来，说James的猜想在美国闹翻天了。这该如何是好？简单，美国懂得“简单的生物信息学分析”的学者比国内多，自己想办法发帖子讲清楚就可以了。
