【让你五行代码发一篇文章】如何把DNA序列编码成深度学习的向量化输入——GP-GCN框架教程

背景

现在越来越多的工作利用深度学习来实现DNA序列的分类或者回归任务，但是如何把DNA序列编码成深度学习的向量化输入？换句话说，怎么把由A,C,G,T组成的DNA“字符串”转化成深度学习需要的向量（矩阵）格式？

目前常用的方法有one-hot encoding和$k$-mer frequency vector。One-hot encoding是把A表示成[1,0,0,0], C表示成[0,1,0,0],G表示成[0,0,1,0],T表示成[0,0,0,1],这样一条长度为$L$的DNA序列可以转成$4\ast L$的矩阵，但多数情况下参与训练的DNA序列长度都不一样，然而深度学习模型要求输入的矩阵大小一致，这时候如果长度差别不是很大，可以用cutting & padding的方法；如果长度差别很大，one-hot encoding就不太适合了。这时候可以采用$k$-mer frequency vector来表示DNA序列，也就是把这条序列中所有的$k$-mer的出现频率或者次数提取出来，组成一个向量作为输入，这种方法的很显然会丢失很多序列的信息，比如$k$-mer之间的interaction。

这里要介绍的GP-GCN framework[1]就是为了解决这些问题。[1]这篇文章主要是拿病毒的序列作为例子来验证这个框架，实际上这个框架可以用于任何DNA序列。

GP-GCN framework原理

如果对原理不感兴趣，可以跳过这部分直接看软件教程哦~

1. sequence to graph
   对于一条DNA序列$s$,我们都可以把它转化成一个graph。graph的点对应$k$-mer，边对应$k$-mer之间的interaction。点的特征是这个$k$-mer出现的频率；边的特征是这两个$k$-mer连在一起出现的频率（因为考虑不同的间隔距离，所以是一个向量）。这里考虑不同的间隔距离是为了解决在生物序列中常出现的基因变异的情况。

2. graph embedding
   接下来使用图卷积神经网络进行消息传递，将$k$-mer的信息通过序列中的信息（保存在边上）聚合到其他$k$-mer上。

3. downstream tasks
   接下来把所有$k$-mer的embeddings放在一起，接上普通的全连接层（也可以先用卷积神经网络降维），形成一个完整的神经网络用于下游分析。

软件教程

软件源码和README可以在GP-GCN Github repository找到。

安装

pip install GCNFrame

由于GP-GCN框架是用pytorch和pytorch_geometric写的，需要下载pytorch 1.7和pytorch_geometric 1.7。

数据准备

• [必要] fasta格式的序列
• [必要] 序列的classification label (需要和fasta顺序一致)
• [可选] 如果除了序列本身，还有其他的特征想要一起训练，可以提供txt文件
• 以上三种数据的example均可在Github repository里找到

代码

# This is an example to train a two-classes model.
from GCNFrame import Biodata, GCNmodel
import torch

device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')

# 1
data = Biodata(fasta_file="example_data/nature_2017.fasta", 
        label_file="example_data/lifestyle_label.txt",
        feature_file="example_data/CDD_protein_feature.txt")
# 2
dataset = data.encode(thread=20)
# 3
model = GCNmodel.model(label_num=2, other_feature_dim=206).to(device)
# 4
GCNmodel.train(dataset, model, weighted_sampling=True)
# 5
GCNmodel.test("GCN_model.pt", fasta_file="test.fasta", feature_file="test_feature.txt")

以上就是training，validation和test的全部代码~
除去前面的引用包以外一共五句，现在来一句句解释：

• #1：读入准备好的数据，包括fasta file，label file 和feature file（如果没有可以写None）
• #2：提取序列中的特征，把序列编码成graph
• #3：建立GCN 模型
• #4：开始训练，如果你的数据里每个类别sample数量不平均，可以选择weighted_sampling=True防止训练产生的bias。 这里的训练会显示training accuracy和validation accuracy，可以利用validation accuracy进行调参。训练完后会自动保存为"GCN_model.pt"。
• #5：用训练好的模型预测新序列的label。这里的fasta_file和feature_file（如果有）是你新的测试集

这里是代码的output：

Encoding sequences...
Epoch 0| Loss: 0.6335| Train accuracy: 0.7480| Validation accuracy: 0.8839
Epoch 1| Loss: 0.5605| Train accuracy: 0.8165| Validation accuracy: 0.7032
Epoch 2| Loss: 0.5042| Train accuracy: 0.8469| Validation accuracy: 0.8065
Epoch 3| Loss: 0.4873| Train accuracy: 0.8344| Validation accuracy: 0.7677
Epoch 4| Loss: 0.4559| Train accuracy: 0.8703| Validation accuracy: 0.8194
Epoch 5| Loss: 0.4533| Train accuracy: 0.8763| Validation accuracy: 0.7806
Epoch 6| Loss: 0.4372| Train accuracy: 0.8931| Validation accuracy: 0.8387
Epoch 7| Loss: 0.4409| Train accuracy: 0.8842| Validation accuracy: 0.8581
Epoch 8| Loss: 0.4357| Train accuracy: 0.8858| Validation accuracy: 0.8516
Epoch 9| Loss: 0.4314| Train accuracy: 0.8987| Validation accuracy: 0.8387
Epoch 10| Loss: 0.4246| Train accuracy: 0.8992| Validation accuracy: 0.8581
Epoch 11| Loss: 0.4085| Train accuracy: 0.9180| Validation accuracy: 0.8839
Epoch 12| Loss: 0.4071| Train accuracy: 0.9290| Validation accuracy: 0.8903
Epoch 13| Loss: 0.4095| Train accuracy: 0.9170| Validation accuracy: 0.8839
Epoch 14| Loss: 0.4019| Train accuracy: 0.9241| Validation accuracy: 0.8839
Epoch 15| Loss: 0.3960| Train accuracy: 0.9342| Validation accuracy: 0.9161

参数

如果对模型比较有兴趣，想通过调参进一步提高模型准确率，可以参考下面的常用参数介绍~

Biodata的参数：

• K：$k$-mer长度（默认3，用2信息量太小，3以上容易内存过大，建议用3）
• d：gapped范围（默认3，表示pattern之间可以容忍的距离，越大包含信息量越多，但容易过拟合）
• thread：线程数（默认10）

GCNmodel.model的参数：

• K: 同上，请与Biodata保持一致
• d: 同上，请与Biodata保持一致
• gcn_layer_num：SAGEConv层数（默认4）

GCNmodel.train的参数：

• learning_rate: 训练速度（默认1e-4）
• batch_size：每批处理的样本的个数（默认64）
• epoch_n: 训练迭代次数（默认20,如果训练准确率不高，可以提高此参数）
• val_split: 选取多少比例的数据作为validation（默认10%）
• weighted_sampling：是否按照类别权重训练（默认是，可以处理不均衡数据集）
• model_name：保存的模型名称（默认"GCN_model.pt"）

GCNmodel.test的参数：

• K: 同上，请与Biodata保持一致
• d: 同上，请与Biodata保持一致

更多的参数可以在GP-GCN Github repository里找到，可以根据validation accuracy对参数进行调整。

GP-GCN框架的介绍就到这里，欢迎大家使用GP-GCN框架来把深度学习用到自己的序列和任务上！如果有任何软件方面的问题，欢迎留言；如果有对软件的建议或者需要更多的功能，也请告诉我，我们非常乐意在后面的版本增加更多的功能！

Reference

[1] Ruohan WANG, Yen Kaow NG, Xianglilan ZHANG, Jianping WANG & Shuai Cheng Li (2022). A graph representation of gapped patterns in phage sequences for graph convolutional network. bioRxiv.