冠状病毒从哪里来呢？专家认为蝙蝠是病毒的来源。他们是怎样得出这个结论的呢?

传统上，大多数医学科学家将使用传统的生物信息学工具，例如BLAST序列比对。然而，在这篇文章中，我将展示如何从数据科学的角度分析冠状病毒基因组。具体来说，我们将使用数据科学方法预测COVID-19的起源，方法是将其基因组序列与其他在鸡、牛、鸭和蝙蝠等动物身上感染过的冠状病毒进行比较。根据我们的分析，我们将展示蝙蝠如何显示出与其他动物相比具有最高的基因组序列相似性。

数据集

可通过美国国家医学图书馆网站（）获得向公众提供的2019冠状病毒的基因组序列。实际上，有一个“更新”警报，您可以在访问网页时立即在顶部发现。

你会注意到病毒基因组序列有两种主要类型:核苷酸和蛋白质。在这里，我们将具体分析4种动物的核苷酸序列，即鸡、蝙蝠、牛和鸭。我们将把它们与2019年从感染者身上收集到的冠状病毒核苷酸序列进行比较。你也可以使用这里学到的方法来分析蛋白质序列。

数据采集

该分析的主要目的是确定病毒的可能来源。在实践中，我们需要用动物间传播的所有病毒扫描COVID-19冠状病毒基因组。然而，由于计算资源有限，我们将只分析传播自鸡、蝙蝠、牛和鸭的冠状病毒。

首先，点击NCBI病毒网页上的核苷酸链接，下载COVID-19冠状病毒的核苷酸序列。

您可以观察到，“Severe acute respiratory syndrome coronavirus 2”是在病毒部分的筛选结果过滤下自动选择的。有两种不同的序列类型:complete或refseq。对于本演练，请从核苷酸序列类型部分的中选择完整的序列类型。接下来，单击页面右上角的蓝色下载按钮。基因组序列将以.fasta文件格式下载到您的计算机。

下载序列数据后，让我们在文本编辑器中将其打开，以查看基因组序列的外观。

您可以观察到第一行以“>”开头，表示该病毒的单行描述。在撰写本文时，NCBI病毒数据库提供了34种COVID-19冠状病毒样本。

接下来，让我们下载四种动物（即鸡，牛，鸭和蝙蝠）的冠状病毒核苷酸序列。

要做到这一点，请从Refine结果的Virus部分删除COVID-19冠状病毒选项(即Severe acute respiratory syndrome coronavirus 2)，然后搜索“ Coronaviradae”。Coronaviradae是所有在动物中传播的冠状病毒的学名。选择此选项将为您提供过去由医学家收集的所有冠状病毒类型的核苷酸序列样品的列表。有趣的是，从搜索结果中您可以看到，首份收集的冠状病毒基因组样本是在1983年收集的！

你可以先收集鸡的核苷酸基因组序列。鸡的学名是gallus gallus。转到“ Refine Results”的“ Host”部分，并搜索Gallus gallus（鸡）。一旦选择了该选项，搜索结果列表将会更新，您可以再次点击下载按钮下载鸡的完整核苷酸序列。

对其余3种动物重复这些步骤：1) cattle (bos taurus), 2) duck (anatidae), and 3) bat (chiropter)。

数据转换

我们的第一个任务是从序列数据fasta文件中读取COVID-19和其他四种动物基因组序列，并将它们放到panda dataframe中。为了实现这一点，我们将只提取序列字符(例如，ATTAAAG)，并忽略所有以“>”开头的单行描述。下面的Python代码通过' > '来识别每个核苷酸序列，并提取两个' > '符号之间的序列字符。然后将提取的序列字符按顺序追加到list中。list包含每种动物基因组的整个字符串。

from nltk.corpus import stopwords  from sklearn.feature_extraction.text import CountVectorizerimport pandas as pddef process_file(filename,target_val):    f = open(filename) #'datasets\\corona-nucleo-chicken-complete.fasta')    lines = ""    s1 = list()    step = 0    term = 0    for line in f:        line = ''.join(line.split())        if line.startswith(">") and step==0:            line = line.split('>',1)[0].strip()            step = step + 1        if line.startswith(">") and step>=1: #and step>=1:            line = line.split('>',1)[0].strip()            s1.append(lines)            lines = ""            step = step + 1            term = 0        lines = lines + line

接下来，我们需要从包含序列字符的字符串生成特征。由于序列是由一系列字母字符表示的，所以我们需要将序列字符分割成tokens。例如，只有两个字符的token，只有三个字符的token，只有四个字符的token等等。此外，我们需要将tokens转换为数字表示格式，以便将它们输入机器学习算法进行预测。一种常见的方法是Bag-of-words数据表示，计算每个token的频率/出现次数。

要应用将n-gram和bag-of-words，我们可以使用scikit-learn机器学习库中的CountVectorizer API 。

count_vect = CountVectorizer(lowercase=False, ngram_range=(2,4),analyzer=’char’) X1 = count_vect.fit_transform(s1)

对于本演练，我们仅将n-grams设置为从2到4个tokens序列开始。例如，2-gram最多具有2个字符（例如AT），3-gram（例如ACG）和4-gram（例如ACTT）。

def generate_ngrams(s1):    count_vect = CountVectorizer(lowercase=False, ngram_range=(2,4),analyzer='char')    X1 = count_vect.fit_transform(s1)        lcount = list()    lcount = []    for i in s1:        count = len(i)        #print(count)        lcount.append(count)    count_vect_df = pd.DataFrame(X1.todense(), columns=count_vect.get_feature_names())        count_vect_df=count_vect_df.apply(lambda x: x / lcount[x.name] ,axis=1)            return count_vect_df

一旦数据被转换成数值作为特征，我们就需要对它们进行标记以进行预测。

df1 = process_file(‘datasets\\corona-seq-chicken.fasta’,”chicken”)df2 = process_file(‘datasets\\corona-seq-duck.fasta’,”duck”)df3 = process_file(‘datasets\\corona-seq-cattle.fasta’,”cattle”)df4 = process_file(‘datasets\\corona-seq-bat.fasta’,”bat”)

在这里，您可以看到每个fasta序列文件都分配有自己的标签（即，鸡，鸭，牛和蝙蝠）。标记是在数据转换过程的最后完成的。

def generate_ngrams(s1):.....count_vect_df = generate_ngrams(s1) count_vect_df[‘target’] = target_valreturn count_vect_df

如果您已经正确地完成了所有工作，您应该能够以panda dataframe格式的数值形式查看所有特性。

鸡核苷酸序列的样本特征

探索性分析

以dataframe格式表示所有序列后，我们可以探索冠状病毒数据集的百分比分布。我们可以看到，我们的数据集中有64％是鸡的核苷酸序列，31%来自蝙蝠，剩下的4%和1%分别来自牛和鸭。

import matplotlib.pyplot as pltplot_size = plt.rcParams[“figure.figsize”]plot_size[0] = 8plot_size[1] = 6plt.rcParams[“figure.figsize”] = plot_sizedf7[‘target’].value_counts().plot(kind=’pie’, autopct=’%1.0f%%’)

现在，我们准备通过将标记的核苷酸数据集应用于机器学习算法来构建预测模型。我们将使用XGBoost梯度提升算法来构建我们的预测模型。python代码如下：

from sklearn.model_selection import train_test_split from xgboost import XGBClassifierfrom xgboost import plot_importanceimport xgboost# create a train/test splitX_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=7, shuffle=True)model = XGBClassifier()model.fit(X_train, y_train)

使用我们建立的机器学习模型，我们可以进一步分析以确定核苷酸序列的哪些区域在确定哪些序列属于哪些动物方面具有最重要的影响因素。我们可以观察到“ CWG”，“ AAAY”和“ AAC”是区分4种动物基因组的前3个序列区域。

COVID-19冠状病毒的数据转换和文本表示

接下来，我们需要将COVID-19核苷酸序列文件加载到我们的程序中以进行预测。

与动物冠状病毒序列相似，我们需要对COVID-19冠状病毒进行数据转换和数字表示。n-gram和Bag-of-words算法都要用到。

首先，我们将从Fasta文件中加载COVID-19冠状病毒的核苷酸。

cov = process_file(‘datasets\\corona-covid-19.fasta’,“COVID-19”)

我们可以看到有409个特征的34个COVID-19样本。注意特征的数量与训练特征不匹配。

因为我们要进行预测(而不是构建模型)，所以可以删除目标标签

cov = cov.drop(‘target’, axis=1)

我们可以观察到，COVID-19冠状病毒的特征数量为409个，与动物冠状病毒1786个特征相比要少很多。由于这种不匹配，在数据集与我们的预测模型完全兼容之前，我们需要执行另外两种类型的数据转换。

1.动物冠状病毒中哪些特征是COVID-19冠状病毒所没有的?

这意味着找到动物核苷酸序列中的特征，但它们在COVID-19冠状病毒中不存在。找到此类特征后，我们将需要在预测数据集中将其填充为“unknown”

mc = X_train.columns.difference(cov.columns)

我们可以看到，“AAAM”、“AAR”和“AAAY”是我们在冠状病毒中发现的一些特征，但在COVID-19冠状病毒中并不存在。

添加这些特征，以使其与我们的训练数据集匹配。因此，我们可以将这些值填充为unknown。由于特征必须为数字形式，因此我们将使用“ -999”数字值表示“unknown”。

for newcol in mc: cov[newcol]=-999

2.COVID-19有哪些其他动物冠状病毒所没有的新特征?

最重要的是，我们还需要找到COVID-19中存在但训练数据集中不存在的特征。找到此类特征后，我们将需要删除这些特征以使其与我们的训练数据集匹配。

rf = cov.columns.difference(X_train.columns)cov = cov.drop(rf, axis=1)

我们可以看到“ AWGC”，“ CS”，“ CST”和“ CSTG”是COVID-19病毒的一些新核苷酸区域，是四种动物冠状病毒中从未出现过的。

在从预测数据集中删除所有这些特性之后，我们的预测数据集最终与我们的训练数据集匹配。

关键时刻

现在我们可以预测COVID-19冠状病毒的可能来源了！下一步是预测我们收集的34个COVID-19冠状病毒样本的可能来源。这可以通过使用XGBoost predict API轻松实现。

model.predict(cov)

根据预测值可以看出，我们的模型预测了全部36个样本为Bat。甚至没有一个样本被预测为其他动物。为了进一步分析这个问题，我们将使用XGBoost的predict_proba API来获得每个预测的概率。

print(model.classes_)similarities = model.predict_proba(cov)np.round(similarities, 3)

根据预测概率，我们可以清楚地看到bat的预测概率极高，预测值为0.957。这意味着，COVID-19病毒的起源有95.7%的可能性来自蝙蝠，而来自鸭子的可能性只有4.1%。病毒来自牛或鸡的概率只有0.1%。因此，我们可以得出结论，蝙蝠是2019年新型冠状病毒的可能来源。