蛋白质负责我们组织、器官和身体的许多活动，它们在细胞的结构和功能中也起着核心作用。蛋白质是由 20 种称为氨基酸的构建单元组成的大分子。人体会产生数以万计的不同蛋白质，每种蛋白质由数十或数百个依次连接的氨基酸组成。该氨基酸序列决定了蛋白质的 3D 结构和构象动力学，进而决定了其生物学功能。由于正在进行的基因组测序项目，我们被来自数千个物种的大量基因组序列数据所淹没，这些数据告诉我们这些基因编码的蛋白质的氨基酸序列数据。准确分配蛋白质的生物功能是在分子水平上理解生命的关键。然而，由于许多蛋白质具有多种功能，以及它们与多个伙伴相互作用的能力，因此为任何特定蛋白质分配功能可能会变得困难。对分配给蛋白质的功能有更多的了解（可能得到数据科学的帮助）可能会导致治愈疾病并改善人类和动物的健康，包括医学和农业等不同领域的健康状况。
研究小组已经开发出许多方法来确定蛋白质的功能，包括许多基于将未解决的序列与功能已知的蛋白质数据库进行比较的方法。其他努力旨在挖掘与其中一些蛋白质相关的科学文献，而更多的方法将复杂的机器学习算法与对生物过程的理解相结合，以破译这些蛋白质的作用。然而，该领域仍然存在许多挑战，这些挑战是由模糊性、复杂性和数据集成驱动的。
背景
基因本体论 （GO） 是一个概念层次结构，它描述了基因和基因产物在不同抽象层次上的生物学功能（Ashburner et al.， 2000）。这是描述蛋白质功能的多面性质的好模型。
GO 是一个有向无环图。此图中的节点是通过它们之间的关系关系（is_a、part_of 等）连接的函数描述符（术语或类）。例如，术语“蛋白质结合活性”和“结合活性”通过is_a关系相关;然而，图中的边缘通常颠倒过来，从结合指向蛋白质结合。该图包含三个子图（子本体）：分子功能 （MF）、生物过程 （BP） 和细胞成分 （CC），由它们的根节点定义。从生物学上讲，每个子图代表蛋白质功能的不同方面：它在分子水平 （MF） 上的作用，它参与哪些生物过程 （BP） 以及它在细胞中的位置 （CC）
因此，蛋白质的功能由一个或多个子本体的子集表示。
这些注释由证据代码支持，证据代码大致可分为实验性（例如，如生物学家研究团队发表的论文中记录的那样）和非实验性。非实验项通常由计算手段推断。我们建议您阅读有关不同类型的 GO 证据代码的更多信息。
我们将使用实验确定的术语蛋白质分配作为每种蛋白质的类别标签。也就是说，如果一个蛋白质用一个术语标记，则意味着该蛋白质具有这种功能，并得到了实验证据的验证。通过处理这些带注释的术语，我们可以为每个术语生成蛋白质及其真值标签的数据集。没有术语注释并不一定意味着蛋白质没有此功能，只是这个注释在 GO 中（还）不存在。一个蛋白质可以由来自同一子本体的一个或多个术语以及来自多个子本体的术语来注释。
训练集
对于训练集，我们包括所有带有注释术语的蛋白质，这些蛋白质已通过实验或高通量证据、可追溯的作者声明（证据代码 TAS）或由策展人 （IC） 推断。有关证据代码的更多信息，请点击此处。我们使用 2022-11-17 版 UniProtKB 中的注释。参与者不需要使用这些数据，也欢迎使用他们可用的任何其他数据。
测试超集
测试超集是一组蛋白质序列，参与者被要求根据这些序列预测 GO 术语。
测试集
测试集在比赛开始时是未知的。它将包含来自测试超集的蛋白质序列（及其功能），这些超集在提交截止日期和评估时间之间获得了实验注释。
文件描述
基因本体： 本体数据在文件中。此结构是 GO 图的 2023-01-01 版本。此文件为 OBO 格式，存在许多解析库。例如，obonet 包可用于 Python。此图中的节点按术语名称编制索引，例如，三个 onotlogies 的根是：go-basic.obo
subontology_roots = {'BPO':'GO:0008150',
'CCO':'GO:0005575',
'MFO':'GO:0003674'}
Training sequences：包含训练数据集的蛋白质序列。
此文件为 FASTA 格式，这是一种用于描述蛋白质序列的标准格式。这些蛋白质都是从欧洲生物信息学研究所整理的 UniProt 数据集中检索的。
该标题包含蛋白质的 UniProt 登录 ID 和有关蛋白质的其他信息。大多数蛋白质序列是从 Swiss-Prot 数据库中提取的，但 Swiss-Prot 中未代表的蛋白质子集是从 TrEMBL 数据库中提取的。在这两种情况下，序列都来自 2022 年 12 月 14 日的 2022_05 版本。更多信息可以在这里找到。train_sequences.fasta
该文件将指示序列源自哪个数据库。例如，在 FASTA 标题中表示具有 UniProt ID 和基因名称的蛋白质取自 Swiss-Prot （）。从 TrEMBL 获取的任何序列都将位于标头中，而不是 .Swiss-Prot 和 TrEMBL 都是 UniProtKB 的一部分。train_sequences.fasta``sp|P9WHI7|RECN_MYCT``P9WHI7``RECN_MYCT``sp``tr``sp
此文件仅包含数据集中带有注释的蛋白质序列（标记的蛋白质）。要获得未标记蛋白质的全套蛋白质序列，可在此处找到 Swiss-Prot 和 TrEMBL 数据库。
Labels：包含 中蛋白质的注释项（基本实况）列表。第一列表示蛋白质的 UniProt 登录 ID，第二列是 GO 术语 ID，第三列表示术语出现在哪个本体中。train_terms.tsv``train_sequences.fasta
Taxonomy：包含蛋白质及其所属物种的列表，由“分类标识符”（taxon ID）编号表示。第一列是蛋白质 UniProt 登录 ID，第二列是分类单元 ID。他可以在此处找到有关分类法的更多信息。train_taxonomy.tsv
信息增加：包含每个 GO 术语的信息增加（权重）。这些权重用于计算加权精度和召回率，如 Evaluation 部分所述。此文件的值是使用以下代码存储库计算的。IA.txt
Test sequences：包含要求参与者提交预测的蛋白质序列。中每个序列的标题都包含蛋白质的 UniProt 登录 ID 和该蛋白质所属物种的分类群 ID。这些序列中只有一小部分会积累功能注释并构成测试集。该文件是测试超集中蛋白质的一组分类 ID。testsuperset.fasta``testsuperset.fasta``testsuperset-taxon-list.tsv
文件
- train_sequences.fasta - 训练集中蛋白质的氨基酸序列
- train_terms.tsv - 蛋白质的训练集和相应的带注释的 GO 术语
- train_taxonomy.tsv - 训练集中蛋白质的分类单元 ID
- go-basic.obo - 本体图结构
- testsuperset.fasta - 应进行预测的蛋白质的氨基酸序列
- testsuperset-taxon-list.tsv - 测试超集中蛋白质的分类单元 ID（注意：您可能需要使用
encoding="ISO-8859-1"``pandas) - IA.txt - 每个术语的信息增加。这用于加权精度和召回率（请参阅 评估）
- sample_submission.csv - 正确格式的样本提交文件
基于TensorFlow的CAFA 5蛋白功能预测
本笔记本将指导您如何在本次比赛提供的CAFA 5蛋白质功能预测数据集上使用TensorFlow训练DNN模型。
该模型的目标是根据一组蛋白质的氨基酸序列和其他数据预测其功能（又称GO术语ID）。
注意：这个笔记本不需要任何GPU。这是因为启用gpu会在虚拟机上留下更少的RAM内存，而提交步骤需要大量内存。这将影响的一点是在训练模型时。在CPU上大约需要2分钟，而在GPU上大约需要30秒。
关于数据
蛋白质序列
每一种蛋白质都由数十或数百个按顺序连接的氨基酸组成。序列中的每个氨基酸可以由一个或三个字母的代码表示。因此，蛋白质的序列通常被标记为一串字母。
train_sequences。Fasta为本次比赛提供，包含带注释的蛋白质序列（标记蛋白质）。
基因本体论
我们可以使用基因本体（GO）定义蛋白质的功能属性。基因本体（GO）从三个方面描述了我们对生物领域的理解：
- 分子功能（MF）
- 生物过程（BP）
- 蜂窝成分（CC）
在这里阅读更多关于基因本体论的内容。
文件train_terms。TSV包含train_sequences.fasta中蛋白质的注释术语列表（基本事实）。在train_terms。tsv第一列表示该蛋白质的UniProt访问ID（唯一蛋白质ID），第二列是GO术语ID，第三列表示该术语出现在哪个本体中。
数据集的标签
我们模型的目标是预测蛋白质序列的项（功能）。一个蛋白质序列可以有多种功能，因此可以分为任意数量的术语。每个词项都由一个GO词项ID唯一标识。因此，我们的模型必须预测蛋白质序列的所有GO术语id。这意味着手头的任务是一个多标签分类问题。
用于训练和测试数据的蛋白质嵌入
为了训练机器学习模型，我们不能使用ain_sequences中的字母蛋白质序列。fasta直接。它们必须被转换成矢量格式。在这个笔记本中，我们将使用蛋白质序列的嵌入来训练模型。你可以认为蛋白质嵌入类似于用于训练NLP模型的词嵌入。蛋白质嵌入是一种机器友好的方法，主要通过其序列来捕获蛋白质的结构和功能特征。一种方法是训练一个自定义的ML模型，以学习在这个笔记本中使用的数据集中蛋白质序列的蛋白质嵌入。由于该数据集使用氨基酸序列表示蛋白质，这是一种标准方法，因此我们可以使用任何公开可用的预训练蛋白质嵌入模型来生成嵌入。
蛋白质嵌入模型有很多种。为了使数据准备更容易，我们在这个笔记本中使用了Sergei Fironov使用Rost实验室的T5蛋白质语言模型创建的预计算的蛋白质嵌入。可以在这里找到预计算的蛋白质嵌入。我们已经将此数据集与比赛可用的数据集一起添加到笔记本中。
要将其添加到您的环境中，在右侧面板上，单击add Data并搜索t5embeds（确保它是正确的），然后单击它旁边的+。
导入所需的库
import tensorflow as tf
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
# Required for progressbar widget
import progressbar
加载数据集
首先加载文件train_terms。TSV，包含蛋白质的注释术语（功能）列表。我们将提取标签，也就是GO术语ID，并为蛋白质嵌入创建一个标签数据框。
train_terms = pd.read_csv("/kaggle/input/cafa-5-protein-function-prediction/Train/train_terms.tsv",sep="\t")
print(train_terms.shape)
Train_terms数据框由3列和5363863个条目组成。我们可以通过下面的代码打印出数据集的前5个元素来查看数据集的所有3个维度：
train_terms.head()
看一下train_terms的第一个条目。tsv，我们可以看到它包含蛋白质id(A0A009IHW8)， GO术语（GO:0008152）和它的aspect（BPO）。
加载蛋白质嵌入
现在，我们将使用Rost实验室的T5蛋白质语言模型加载Sergei Fironov创建的预计算的蛋白质嵌入。
如果tfembeds还没有出现在notebook的输入数据上，你可以点击add data并搜索t5embeds（确保它是正确的），然后点击它旁边的+。
用于训练的蛋白质嵌入被记录在train_embedding中。Npy及其对应的蛋白质id可以在train_ids.npy中找到。
首先，我们将加载train_ids中包含的训练数据集中的蛋白质嵌入的蛋白质id。Npy转换成numpy数组。
train_protein_ids = np.load('/kaggle/input/t5embeds/train_ids.npy')
print(train_protein_ids.shape)
train_protein_ids数组包含142246个protein_ids。使用下面的代码打印出前5个条目：
train_protein_ids[:5]
在将文件加载为numpy数组后，我们将把它们转换为Pandas dataframe。
每个蛋白质嵌入是一个长度为1024的向量。我们创建结果dataframe，使其有1024列来表示向量中1024个位置中的每个值。
train_embeddings = np.load('/kaggle/input/t5embeds/train_embeds.npy')
# Now lets convert embeddings numpy array(train_embeddings) into pandas dataframe.
column_num = train_embeddings.shape[1]
train_df = pd.DataFrame(train_embeddings, columns = ["Column_" + str(i) for i in range(1, column_num+1)])
print(train_df.shape)
(142246, 1024)
train_df数据帧包含嵌入，由1024列和142246项组成。我们可以用下面的代码打印出数据集的前5个元素，得到所有1024维的数据（由于列的长度太长，结果将被截断）：
train_df.head()
准备数据
首先，我们将从train_terms中提取所有需要的标签（GO术语ID）。tsv文件。有4万多个标签。为了简化我们的模型，我们将选择最频繁的1500个GO词条id作为标签。
让我们绘制train_terms.tsv中出现频率最高的100个GO词条id。
# 选择前1500个值进行绘图
plot_df = train_terms['term'].value_counts().iloc[:100]
figure, axis = plt.subplots(1, 1, figsize=(12, 6))
bp = sns.barplot(ax=axis, x=np.array(plot_df.index), y=plot_df.values)
bp.set_xticklabels(bp.get_xticklabels(), rotation=90, size = 6)
axis.set_title('Top 100 frequent GO term IDs')
bp.set_xlabel("GO term IDs", fontsize = 12)
bp.set_ylabel("Count", fontsize = 12)
plt.show()
现在我们将前1500个最频繁的GO词条id保存到一个列表中。
# 设定标签的限制
num_of_labels = 1500
# 以降序排列的值计数并获取前1500个`GO term ID`作为标签
labels = train_terms['term'].value_counts().index[:num_of_labels].tolist()
接下来，我们将使用选定的GO词条id过滤训练词条，从而创建一个新的dataframe。
# 只获取相关标签的train_terms数据
train_terms_updated = train_terms.loc[train_terms['term'].isin(labels)]
让我们使用饼图来绘制新的train_terms_updated数据框中的纵横比值。
pie_df = train_terms_updated['aspect'].value_counts()
palette_color = sns.color_palette('bright')
plt.pie(pie_df.values, labels=np.array(pie_df.index), colors=palette_color, autopct='%.0f%%')
plt.show()
如您所见，大多数GO术语id都以BPO（生物过程本体）作为它们的方面。
由于这是一个多标签分类问题，在labels数组中，我们将使用1或0表示蛋白质Id的每个Go术语Id的存在或不存在。首先，我们将创建一个numpy数组train_labels，其大小为标签所需的大小。为了用适当的值更新train_labels数组，我们将遍历label列表。
# 设置进度条设置。
# 这完全是为了美观。
bar = progressbar.ProgressBar(maxval=num_of_labels, \
widgets=[progressbar.Bar('=', '[', ']'), ' ', progressbar.Percentage()])
# 创建一个所需大小的空数据框用于存储标签，即train_size x num_of_labels （142246 x 1500）
train_size = train_protein_ids.shape[0] # len(X)
train_labels = np.zeros((train_size ,num_of_labels))
# 将numpy转换为pandas series以便更好地处理
series_train_protein_ids = pd.Series(train_protein_ids)
# 循环遍历每个标签
for i in range(num_of_labels):
# 对于每个标签，获取相应的train_terms数据
n_train_terms = train_terms_updated[train_terms_updated['term'] == labels[i]]
# 获取所有与当前标签相关的唯一EntryId （GO术语ID）
label_related_proteins = n_train_terms['EntryID'].unique()
# 在series_train_protein_ids pandas序列中，如果一个蛋白质与当前标签相关，则将其标记为1，否则为0。
# 用pandas序列替换train_Y的第i列。
train_labels[:,i] = series_train_protein_ids.isin(label_related_proteins).astype(float)
# 进度条百分比增加
bar.update(i+1)
# 通知进度条结束
bar.finish()
# 将train_Y numpy转换为pandas dataframe
labels_df = pd.DataFrame(data = train_labels, columns = labels)
print(labels_df.shape)
(142246, 1500)
最终的labels数据框（label_df）由1500列和142246个条目组成。我们可以用下面的代码打印出数据集的前5个元素，从而得到所有1500维的数据集（由于列的数量很大，结果将被截断）：
labels_df.head()
训练
接下来，我们将使用Tensorflow来训练带有蛋白质嵌入的深度神经网络。
INPUT_SHAPE = [train_df.shape[1]]
BATCH_SIZE = 5120
model = tf.keras.Sequential([
tf.keras.layers.BatchNormalization(input_shape=INPUT_SHAPE),
tf.keras.layers.Dense(units=512, activation='relu'),
tf.keras.layers.Dense(units=512, activation='relu'),
tf.keras.layers.Dense(units=512, activation='relu'),
tf.keras.layers.Dense(units=num_of_labels,activation='sigmoid')
])
# Compile model
model.compile(
optimizer=tf.keras.optimizers.Adam(learning_rate=0.001),
loss='binary_crossentropy',
metrics=['binary_accuracy', tf.keras.metrics.AUC()],
)
history = model.fit(
train_df, labels_df,
batch_size=BATCH_SIZE,
epochs=5
)
绘制每个epoch的模型损失和精度
history_df = pd.DataFrame(history.history)
history_df.loc[:, ['loss']].plot(title="Cross-entropy")
history_df.loc[:, ['binary_accuracy']].plot(title="Accuracy")
提交结果
在提交过程中，我们将使用Sergei Fironov使用Rost实验室的T5蛋白质语言模型创建的测试数据的蛋白质嵌入。
test_embeddings = np.load('/kaggle/input/t5embeds/test_embeds.npy')
# Convert test_embeddings to dataframe
column_num = test_embeddings.shape[1]
test_df = pd.DataFrame(test_embeddings, columns = ["Column_" + str(i) for i in range(1, column_num+1)])
print(test_df.shape)
(141865, 1024)
test_df由1024列和141865项组成。我们可以用下面的代码打印出数据集的前5个元素，得到所有1024维的数据（由于列的长度太长，结果将被截断）：
test_df.head()
现在，我们将使用该模型对测试嵌入进行预测。
predictions = model.predict(test_df)
4434/4434 [==============================] - 24s 5ms/step
df_submission = pd.DataFrame(columns = ['Protein Id', 'GO Term Id','Prediction'])
test_protein_ids = np.load('/kaggle/input/t5embeds/test_ids.npy')
l = []
for k in list(test_protein_ids):
l += [ k] * predictions.shape[1]
df_submission['Protein Id'] = l
df_submission['GO Term Id'] = labels * predictions.shape[0]
df_submission['Prediction'] = predictions.ravel()
df_submission.to_csv("submission.tsv",header=False, index=False, sep="\t")
