去年的这个时候，欧密克戎变异体已经成为主要的严重急性呼吸综合征冠状病毒2型（SARS-CoV-2）变异体。
在今年的过程中，一些欧密克戎亚变异体通过突变和在某些情况下的重组交替占据主导地位。欧密克戎XBB亚变异体（在本文中缩写为Omicron XBB）是2022年8月以其高传播性而引起注意的重组亚变异体之一。
然而，写这篇文章的主要动机不是研究Omicron XBB的演化，而是要证明在生物信息学等其他数据科学子领域中，有时候精心构造的小数据可以产生与大数据相当的洞察力。
在写这篇文章时，GISAID上已经上传了约2300个Omicron XBB的全基因组序列。为了筛选我的序列，我选择了以下条件的序列：
- 完整的序列。
- 覆盖范围广的序列。
- 附带患者状态信息的序列。
这样，将序列缩小到了414个质量良好的可供下载的序列。
一旦有了质量良好的全基因组序列，就是时间来推断洞察力了。为了从生物信息中获取尽可能多的洞察力，我通常使用以下6个Python包。
代码语言：javascript复制# 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中。
代码语言：javascript复制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))在这项研究中，选择关注编码剌突蛋白的剌突基因。剌突蛋白是冠状病毒用于进入宿主细胞的关键。此外，这也展示了生物信息学科学家如何从生物序列中切割出所需的片段。
由于完整的序列的核苷酸碱基数目不同，因此使用了一个近似的剌突基因位点，使得所有剌突基因都能够被纳入，即使是具有逐渐变小的末端。
代码语言：javascript复制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值相关，以提供更多有关观察到的数据集中的进化模式的细节。
代码语言：javascript复制selection_table = pd.read_csv('OmicronXBBdatamonkey-table.csv')
selection_table
通常情况下，负选择/纯化选择不会提供足够的信息，除非这是研究的目的，因为我们观察到的是稳定其构象的特定位点。由于这种稳定性，某些位点非常适合设计配体/药物，因为它们很少发生变化。
主要对正选择/适应性选择感兴趣，因为它给我们提供了病毒如何进化的想法，当某个突变逐渐出现在病毒群体中时，它应该提供比没有该突变的病毒更有优势。
如果对负选择感兴趣，可以取消注释涉及负选择的代码。
代码语言：javascript复制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')
正选择位点的贝叶斯因子
喜欢使用贝叶斯因子来理解正选择，因为它可以给我们明显的突出峰值，表示强烈的正选择位点。
代码语言：javascript复制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表明根据所选择的特征，显示出强正选择或强负选择的证据，因此下面的代码块有助于我们找到这种模式。
代码语言：javascript复制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也可以进行转换，只是需要更长的时间。
代码语言：javascript复制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
用于描述患者状态的一些不同术语意味着相同的含义，这是因为这些序列是由世界各地的不同实验室生成和上传的。
因此，以手动方式完成，因为找不到更好的使用代码的方法。同样，如果处理数百万条患者记录，这将是耗时的。
代码语言：javascript复制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
所以从上面生成的饼状图可以清楚地看出，Omicron XBB感染者中有34%住院，3%因Omicron XBB死亡，31%已出院。感染Omicron XBB的COVID-19患者中，30%是步行或没有住院治疗。
医疗专业人员可以根据这些统计数据推断出Omicron XBB的致命性或对健康的负担，从而帮助他们采取健康措施来遏制病毒的传播。
从理论上讲，使用大数据或数百万个生物序列，由于可以选择的样本数量较多，洞察力的准确性会增加。
然而，随着数据规模的增加，查找个别的“未知”或错误变得冗长和耗时。通过积累这些未知或错误，数据会积累噪音，可能会干扰下游处理中使用的算法的顺利执行。
除了流畅地执行精心设计的小数据之外，生成见解也更快，因此你可以有更多的时间来分析数据，从中提取有用的信息。唯一要付出的代价就是精心制作。
☆ END ☆
