动手机器学习之序列数据的预处理
随着生物信息学的快速发展，DNA 序列数据呈爆炸式增长。这些序列蕴含着丰富的生物学信息。然而，对一些处理数值数据的机器学习模型来说，原始的 DNA 序列是一串字符，无法直接作为输入。因此，将 DNA 序列转化为机器学习模型能够理解的数值矩阵形式成为一个关键步骤。
本案例以识别 E.Coli 启动子序列的二分类器来介绍转换序列数据的一种简单方式。
数据预览
import numpy as np
import pandas as pd
url = 'https://archive.ics.uci.edu/ml/machine-learning-databases/molecular-biology/promoter-gene-sequences/promoters.data'
names = ['Class', 'id', 'Sequence']
data = pd.read_csv(url, names = names)
print(data.iloc[0])
print("\n", data.head)
输出：
Class +
id S10
Sequence \t\ttactagcaatacgcttgcgttcggtggttaagtatgtataat...
Name: 0, dtype: object
<bound method NDFrame.head of Class id Sequence
0 + S10 \t\ttactagcaatacgcttgcgttcggtggttaagtatgtataat...
1 + AMPC \t\ttgctatcctgacagttgtcacgctgattggtgtcgttacaat...
2 + AROH \t\tgtactagagaactagtgcattagcttatttttttgttatcat...
3 + DEOP2 \taattgtgatgtgtatcgaagtgtgttgcggagtagatgttagaa...
4 + LEU1_TRNA \ttcgataattaactattgacgaaaagctgaaaaccactagaatgc...
.. ... ... ...
101 - 799 \t\tcctcaatggcctctaaacgggtcttgaggggttttttgctga...
102 - 987 \t\tgtattctcaacaagattaaccgacagattcaatctcgtggat...
103 - 1226 \t\tcgcgactacgatgagatgcctgagtgcttccgttactggatt...
104 - 794 \t\tctcgtcctcaatggcctctaaacgggtcttgaggggtttttt...
105 - 1442 \t\ttaacattaataaataaggaggctctaatggcactcattagcc...
[106 rows x 3 columns]>
序列的提取与分割
DataFrame 可以看作是由多个 Series 组成的，，每个 Series 代表 DataFrame 的一列。我们要想办法关于序列的信息提取到 DataFrame 的每一列中。
classes = data.loc[:, 'Class']
print(classes[:5])
sequences = list(data.loc[:, 'Sequence'])
print(sequences[:5])
输出：
0 +
1 +
2 +
3 +
4 +
Name: Class, dtype: object
['\t\ttactagcaatacgcttgcgttcggtggttaagtatgtataatgcgcgggcttgtcgt', '\t\ttgctatcctgacagttgtcacgctgattggtgtcgttacaatctaacgcatcgccaa', '\t\tgtactagagaactagtgcattagcttatttttttgttatcatgctaaccacccggcg', '\taattgtgatgtgtatcgaagtgtgttgcggagtagatgttagaatactaacaaactc', '\ttcgataattaactattgacgaaaagctgaaaaccactagaatgcgcctccgtggtag']
字符串列表 sequences 包含 DNA 序列原始数据中完整的字符串，需要进一步分割成字符，即单个碱基。
dataset = {}
for i, seq in enumerate(sequences):
nucleotides = seq
# 删除原始数据中的制表符
nucleotides = [x for x in nucleotides if x != '\t']
# 写入序列是否为启动子的信息
nucleotides.append(classes[i])
# 合成 series
dataset[i] = nucleotides
print(dataset[0])
输出：
['t', 'a', 'c', 't', 'a', 'g', 'c', 'a', 'a', 't', 'a', 'c', 'g', 'c', 't', 't', 'g', 'c', 'g', 't', 't', 'c', 'g', 'g', 't', 'g', 'g', 't', 't', 'a', 'a', 'g', 't', 'a', 't', 'g', 't', 'a', 't', 'a', 'a', 't', 'g', 'c', 'g', 'c', 'g', 'g', 'g', 'c', 't', 't', 'g', 't', 'c', 'g', 't', '+']
行列的转置
训练数据通常以矩阵的每一行作为样本，因此，需要对由 Series 合成的数据进行转置。
# 行列转置
df = pd.DataFrame(dataset).transpose()
# 类别列的重命名
df.rename(columns = {57: 'Class'}, inplace = True)
print(df.iloc[:5])
输出：
0 1 2 3 4 5 6 7 8 9 ... 48 49 50 51 52 53 54 55 56 Class
0 t a c t a g c a a t ... g c t t g t c g t +
1 t g c t a t c c t g ... c a t c g c c a a +
2 g t a c t a g a g a ... c a c c c g g c g +
3 a a t t g t g a t g ... a a c a a a c t c +
4 t c g a t a a t t a ... c c g t g g t a g +
[5 rows x 58 columns]
数值转换：One-hot 编码
One-hot 编码的原理是，将每个碱基用一个四维向量表示。例如，四个维度分别为a、c、g、t，其中对应碱基的位置为 1，其余为 0，则 A 可以表示为 [1,0,0,0]。
# 生成整型虚拟变量，即 0 或 1
df_encoded = pd.get_dummies(df, dtype=int)
# 删除多余的类别列
df = df_encoded.drop(columns=['Class_-'])
df.rename(columns = {'Class_+': 'Class'}, inplace = True)
df.iloc[:5]
输出：
0_a 0_c 0_g 0_t 1_a 1_c 1_g 1_t 2_a 2_c ... 54_t 55_a 55_c 55_g 55_t 56_a 56_c 56_g 56_t Class
0 0 0 0 1 1 0 0 0 0 1 ... 0 0 0 1 0 0 0 0 1 1
1 0 0 0 1 0 0 1 0 0 1 ... 0 1 0 0 0 1 0 0 0 1
2 0 0 1 0 0 0 0 1 1 0 ... 0 0 1 0 0 0 0 1 0 1
3 1 0 0 0 1 0 0 0 0 0 ... 0 0 0 0 1 0 1 0 0 1
4 0 0 0 1 0 1 0 0 0 0 ... 1 1 0 0 0 0 0 1 0 1
5 rows × 229 columns
DNA 序列每一个位置上都被“扩张”为一个四维向量，DataFrame 的行数增加至原来的 4 倍左右。
df_encoded 的行数为 57 × 4 + 1 × 2 = 230，其中 57 为序列的长度，1 × 2 是分类类别的编码。
训练模型
以 AdaBoost 为例，训练简单的模型。
from sklearn.ensemble import AdaBoostClassifier
scoring = 'accuracy'
seed = 1
clf = AdaBoostClassifier()
kfold = model_selection.KFold(n_splits=5, random_state = seed, shuffle = True)
cv_results = model_selection.cross_val_score(clf, X_train, y_train, cv=kfold, scoring=scoring)
msg = "AdaBoost: %f (std %f) \n" % (cv_results.mean(), cv_results.std())
print(msg)
clf.fit(X_train, y_train)
y_pred = clf.predict(X_test)
print("Classification report: \n", classification_report(y_test, y_pred))
输出：
AdaBoost: 0.902941 (std 0.085040)
Classification report:
precision recall f1-score support
0 1.00 0.93 0.97 15
1 0.88 1.00 0.93 7
accuracy 0.95 22
macro avg 0.94 0.97 0.95 22
weighted avg 0.96 0.95 0.96 22
总结
将 DNA 序列转化为机器学习可用的矩阵数据，是生物信息学中不可或缺的预处理环节。多余字段的清洗、字符串的分割、以及行列信息的调整等，是处理序列数据的基础。除此之外，在进行数值转换时，不同方法各具优劣，需根据具体任务与数据特点加以权衡。虽然本案例中的 One-hot 编码简单直观，保留了碱基的原始信息，但在处理长序列时可能引发维度灾难。若操作对象为氨基酸序列，维度扩张会更加严重，造成计算开销的增大。不同方法的选择应结合模型的训练效果及数据的表现进行优化，以在保证精度的前提下减少计算成本。