少即是多：精心构造的小数据也可以产生与大数据相当的洞察力

磐创AI

发布于 2023-11-03 15:43:09

1730

发布于 2023-11-03 15:43:09

文章被收录于专栏：磐创AI技术团队的专栏

去年的这个时候，欧密克戎变异体已经成为主要的严重急性呼吸综合征冠状病毒2型（SARS-CoV-2）变异体。

在今年的过程中，一些欧密克戎亚变异体通过突变和在某些情况下的重组交替占据主导地位。欧密克戎XBB亚变异体（在本文中缩写为Omicron XBB）是2022年8月以其高传播性而引起注意的重组亚变异体之一。

然而，写这篇文章的主要动机不是研究Omicron XBB的演化，而是要证明在生物信息学等其他数据科学子领域中，有时候精心构造的小数据可以产生与大数据相当的洞察力。

在写这篇文章时，GISAID上已经上传了约2300个Omicron XBB的全基因组序列。为了筛选我的序列，我选择了以下条件的序列：

完整的序列。
覆盖范围广的序列。
附带患者状态信息的序列。

这样，将序列缩小到了414个质量良好的可供下载的序列。

一旦有了质量良好的全基因组序列，就是时间来推断洞察力了。为了从生物信息中获取尽可能多的洞察力，我通常使用以下6个Python包。

# I'm a fan of importing all the packages I'll use in the first cell
# Useful for pinpointing which package is not executing properly

import numpy as np
import pandas as pd
import Bio     #Biopython
from collections import Counter
import matplotlib.pyplot as plt
import seaborn as sns

使用Biopython（Bio）来解析/读取fasta文件中的核苷酸序列，使用打印函数一瞥文件内部的内容。fasta序列存储在变量omicron中。

from Bio import SeqIO
# Biopython is useful for loading biological sequence

omicron = SeqIO.parse("Omicron XBB Complete Sequences.fasta", 'fasta')


for seq_record in omicron:

    print(seq_record.id)            #taking a glimpse of what's in the file
    print(repr(seq_record.seq))

在这项研究中，选择关注编码剌突蛋白的剌突基因。剌突蛋白是冠状病毒用于进入宿主细胞的关键。此外，这也展示了生物信息学科学家如何从生物序列中切割出所需的片段。

由于完整的序列的核苷酸碱基数目不同，因此使用了一个近似的剌突基因位点，使得所有剌突基因都能够被纳入，即使是具有逐渐变小的末端。

from Bio import SeqIO
# I had to recall the Biopython package, for some reason it could not work without being recalled

#THis cell is slicing the approximate spike genes from SARS-CoV-2 OmicronXBB whole genome sequences.
omicron_spike = []

for seq_record in SeqIO.parse("Omicron XBB Complete Sequences.fasta",'fasta'):
    spike_gene = seq_record[21500:25499]
    omicron_spike.append(spike_gene)


SeqIO.write(omicron_spike,'omicron_spike.fasta','fasta') 
# writes omicron_spike  contents to a fasta file with name provided.

然后，使用EBI Webserver上的Clustal对近似的剌突序列进行了比对。
使用Jalview来删除冗余/重复的序列。
使用Mega X编辑比对，使用GISAID剌突参考序列作为指南，去除剌突基因的逐渐减少的末端。
清理文件并通过Datamonkey网服务器上传到FUBAR进行选择分析。
分析后导出的.csv文件。

调查

然后使用pandas加载从FUBAR导出的.csv文件以进行进一步分析。下面简单介绍一下这些术语的含义:

site=它编码一个氨基酸，所以它相当于氨基酸的位置。

alpha = 同义替代率，表示编码的氨基酸改变后不改变编码的氨基酸的频率。

beta = 非同义替代率，表示编码的氨基酸改变后改变编码的氨基酸的频率。

其他参数与alpha和beta值相关，以提供更多有关观察到的数据集中的进化模式的细节。

selection_table = pd.read_csv('OmicronXBBdatamonkey-table.csv')
selection_table

通常情况下，负选择/纯化选择不会提供足够的信息，除非这是研究的目的，因为我们观察到的是稳定其构象的特定位点。由于这种稳定性，某些位点非常适合设计配体/药物，因为它们很少发生变化。

主要对正选择/适应性选择感兴趣，因为它给我们提供了病毒如何进化的想法，当某个突变逐渐出现在病毒群体中时，它应该提供比没有该突变的病毒更有优势。

如果对负选择感兴趣，可以取消注释涉及负选择的代码。

sites_positive = []
# sites_negative = []

for idx, val in enumerate(selection_table['beta-alpha']):
    if val > 1:
        sites_positive.append(idx)

    else:
        pass



number_of_positive_sites = str(len(sites_positive))
# number_of_negative_sites = str(len(sites_negative))


print(f'{number_of_positive_sites} sites show adaptive selection')
# print(f'{number_of_negative_sites} sites show purifying selection')

正选择位点的贝叶斯因子

喜欢使用贝叶斯因子来理解正选择，因为它可以给我们明显的突出峰值，表示强烈的正选择位点。

site = selection_table['Site']
bayes_factor = selection_table['BayesFactor[alpha<beta]']

plt.figure(figsize=(16,10), dpi=1200)
ax = plt.axes()
ax.plot(site, bayes_factor)
plt.show;

显示强正选择和强负选择位点

上面的图表显示了强正选择位点呈尖峰状，负选择位点呈小/平坦峰状的一般趋势。

根据FUBAR的说法，后验概率大于0.9表明根据所选择的特征，显示出强正选择或强负选择的证据，因此下面的代码块有助于我们找到这种模式。

strong_beta = selection_table['Prob[alpha<beta]']
strong_alpha = selection_table['Prob[alpha>beta]']
# The Posterior probability where alpha (synonymous rates) < beta (non-synonymous rates/var(strong_beta)) is more than 0.9 were selected as evidence of strong positive selection 
# The Posterior probability where alpha (synonymous rates) > beta (non-synonymous rates/var(strong_alpha)) is more than 0.9 were selected as evidence of strong negative selection 

strong_positive_selection = selection_table.loc[strong_beta > 0.9] 
strong_negative_selection = selection_table.loc[strong_alpha > 0.9]

# print(Counter(strong_positive_selection))

#Then we combine the two strong_selection dataframes.
strong_selection_dfs = [strong_negative_selection,strong_positive_selection]
strong_selection_sites = pd.concat(strong_selection_dfs)
strong_selection_sites

讨论上述表格和图表所展示的趋势和数据，将值得撰写一篇独立的文章，并且这将偏离本文的主旨。

Omicron XBB对患者的影响

嗯，我们经过漫长的旅程才从我们的数据中得到了数字和趋势，但到目前为止，只有少数了解我们如何处理数据的人认为生成的信息很酷。如果它不能帮助我们理解COVID-19，它将变成一种爱好，而爱好无法为我们提供资金，哈哈。

下一部分简要介绍了这些突变如何对健康产生影响。

请记住，用于选择这些序列的过滤器之一是患者数据，并且我将根据使用情况演示使用这些数据的一种方法。

首先选择了患者状态，但还有其他字段需要考虑，如性别、地点、采集日期等其他参数。

患者数据文件以.tsv（制表符分隔值）格式下载，但发现很难使用，所以使用在线工具将其转换为.csv（逗号分隔值）文件。如果离线工作，MS Excel也可以进行转换，只是需要更长的时间。

patient_data = pd.read_csv('Omicron XBB Patient Data.csv')
patient_data

patient_status  #some of the fields can be combined as they mean the same to reduce noise

用于描述患者状态的一些不同术语意味着相同的含义，这是因为这些序列是由世界各地的不同实验室生成和上传的。

因此，以手动方式完成，因为找不到更好的使用代码的方法。同样，如果处理数百万条患者记录，这将是耗时的。

generalised_patient_status = ('Hospitalised','Ambulatory','Deceased','Released','Unknown')
patient_numbers = (135,117,12,122,8) 

labels = list(generalised_patient_status)
size = list(patient_numbers)

sns.set()
plt.figure(figsize=(16,10))
plt.pie(size, labels=labels, autopct='%1.0f%%');
plt.title('Patient statuses after testing postive for Omicron XBB')
plt.figure(dpi=1200); #for image quality