非典是哪一年发生的 非典是哪一年结束了

公众号：医学大数据挖掘分析

0 背景

2019年末出现的新型冠状病毒SARS-COV-2对全球人类健康构成了无时无刻不在的威胁。在2020年3月初，世界卫生组织认为COVID-19（SARS-COV-2引起的疾病）可被定为大流行病。截止到2021年3月2日，全球范围内的确诊病例已经超过1.1亿，死亡病例超过250万。在往期推文中（：利用Nextstrain全面了解新冠病毒）我们已经介绍过一个研究SARS-CoV-2的传播和演变的开源项目。

通过生信分析对新冠病毒的基因序列进行挖掘，了解新冠病毒的传播和演变，对采取有效的公共卫生措施和监测非常重要。今天，小编就跟大家继续拓展生信分析系列，分享一下如何对新冠病毒的基因序列进行分析（下文所举例子对2020年4月之前采集的病毒数据进行分析）。

图1 数据统计来源于Google

1 数据准备

数据来源

目前大部分新冠病毒的基因序列存放在GISAID（Global Initiative on Sharing All Influenza Data，https://www.gisaid.org，全球共享流感数据倡议组织）数据库上。

GISAID建立于2008年5月，是由全世界一些较权威的医学科学家组建，致力于改善流感数据的共享。目前GISAID上已经存放了超过58万条来自世界各地实验室的新冠病毒序列数据。除了序列信息外，GISAID还提供序列的meta数据，包括病毒采集的地方和日期，序列提交的时间，序列的质量，提交的实验室等。

在获取GISAID上的数据之前，需要先注册账号，注册完成后还需要经过审核，审核通过后就能获得下载权限了。

数据清洗

原始数据中有一些测序质量较低的序列需要去除，防止这些数据对下游分析造成干扰。从GISAID下载的序列带有测序质量的信息，可以根据这些信息对序列进行清洗。

以CNCB的新冠病毒数据库为例，如果序列的未知碱基大于15个，或者简并碱基的数量大于50个，则被标记为低质量序列。在分析时，我们可以把这些序列剔除掉。

展开全文

数据处理

获得干净的序列数据后，接下来就可以做比对分析了。比对分析是指把各条序列与参考序列相同位置上的碱基作比较，以确定突变信息。

比对工具有很多种，小编用的比对工具是MAFFTv7.450，使用方法可以参考官方文档（https://mafft.cbrc.jp/alignment/software/）。目前公认的参考序列的GISAID编号是EPI_ISL_402125。同时，GISAID上也有比对好的结果文件供下载。

2 分析方向

经过以上的预处理后，就可以开始进行我们的数据挖掘工作了。分析的方向有很多，这里简要地介绍一下突变分析和进化分析。

突变分析

通过比对我们知道了各个突变所处的位置，要进一步统计这些突变发生在哪个区域，需要对参考序列进行注释。注释是指确定病毒序列上的编码区，例如最重要的编码病毒S蛋白在序列中的位置。

序列的注释信息可以从NCBI中获得（https://www.ncbi.nlm.nih.gov/sars-cov-2/）。通过统计突变出现的次数，我们可以知道哪些是高频突变以及它们出现的编码区，并作进一步的深入分析，如图2是小编分析的各种高频突变在新冠基因组上的分布。

图2

比如，之前受到很大关注的突变D614G（对应图2中的A23403G突变）出现在编码新冠病毒S蛋白的区域中，而S蛋白是与细胞接触进而感染细胞的关键。为了解这个突变与病毒的传染能力是否有关，我们可以进一步分析这个突变是否改变了S蛋白的结构导致病毒与细胞的亲和力更高等。另外，对共同突变的分析可以揭示新冠病毒的多样性。

例如，小编发现根据共同突变能够把病毒分为18个病毒亚群，它们的分布在各大洲有明显的区别（图3）。

图3

进化分析

通过进化分析，我们可以追踪病毒的演变情况，推断病毒的传播路径以及估计病毒出现的时间等。

BEAST是对病毒做进化分析强有力的工具，被应用在很多包括新冠在内的病毒研究中。BEAST实际上是一组工具套件，其中包括BEAUti，BEAST和TreeAnnotator。BEAUti的作用是生成一个参数文件，BEAST读取这个参数文件后进行分析得到结果，最后TreeAnnotator对BEAST产生的结果作进一步的处理。接下来小编跟大家分享一下如何使用这款工具。

首先准备输入文件。BEAUti接受Nexus和Fasta格式的序列文件，如果输入文件不是这两种格式，则需要进行转换，这里小编使用Aliview（https://ormbunkar.se/aliview/）进行格式转换。接下来要设置一系列的参数，这里小编给出自己在分析过程中使用的参数（表1），对参数更加详细的解释请参考官方文档（https://beast.community/workshop_rates_and_dates）。

得到的参数文件输入到BEAST中然后执行程序，完成后会产生两个文件，一个是log文件（以.log结尾），另一个是trees文件（以.trees结尾）。Trees文件需要用TreeAnnotator作进一步处理，这里需要设置BurnIn参数，如果选择Burnin(as states)，则一般设置为马尔科夫链步数的10%；如果选择Burnin(as trees)，则一般设置为树总数的10%（步数与树总数的数量关系为1000:1）。产生的树文件小编推荐使用iTOL（https://itol.embl.de/）进行可视化。图4是小编分析的结果，并且把上文提到的18个病毒群使用iTOL进行了颜色标注。

图4

Log文件需要用Tracer（https://beast.community/tracer）进行读取，这里burnin会默认设置为马尔科夫链步数的10%，可以手动调整。在Traces栏目中，会显示一系列的统计结果，当统计结果的ESS大于200时，这个结果可信度更高，如果ESS<200，可以设置更加长的马尔科夫链，或者调整Burnin参数。

点击某个统计结果，会在右边栏显示详细信息。在这里age(root)统计量表示新冠病毒预估的出现时间，点击这个统计量，我们可以知道病毒出现的时间。

用小编的数据分析，发现病毒大概是在2019年11月3日（95%置信区间：[2019.10.12，2019.11.25]）出现的，比最早采集到的样本早一个月左右。要注意的是，这里小编只是举个例子，为了节省分析时间，设置了较短的步长，导致ESS小于200，所以结果不一定可靠。

图5

以上就是小编今天跟大家分享的全部内容啦，为了方便大家的理解，有一些技术细节没有详细展开，有疑问的同学欢迎后台联系小编一起讨论交流噢！

更多内容，请关注“医学大数据挖掘分析”公众号，欢迎留言联系～

邮箱：medicalda@tp-data.com

地址：广州市天河区珠江东路高德置地秋广场F座