使用 Keras 和 LSTM 进行 DNA 序列分类 (修复未知碱基问题)

本代码使用 Keras 和 LSTM 模型对 DNA 序列进行分类，并解决处理未知碱基 (N) 的问题。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.model_selection import train_test_split
from keras.models import Sequential
from keras.layers import Dense, Flatten, Conv1D, MaxPooling1D, LSTM, Reshape
from keras.utils import to_categorical

def read_fasta(file):
    seqs = []
    labels = []
    with open(file, 'r') as f:
        for line in f:
            if line.startswith('>'):
                labels.append(line.strip()[1:])
            else:
                seqs.append(line.strip())
    label_dict = {l:i for i,l in enumerate(set(labels))}
    label_indices = [label_dict[l] for l in labels]
    seqs = np.array(seqs)[np.argsort(label_indices)]
    labels = np.array(label_indices)
    return seqs, labels

def seq_to_num(seqs):
    seq_num = []
    for seq in seqs:
        seq_num.append([int(i) if i in '1234' else 0 for i in seq.replace('A', '1').replace('C', '2').replace('G', '3').replace('T', '4')])
    return np.array(seq_num)

def encode_labels(labels):
    label_dict = {l:i for i,l in enumerate(set(labels))}
    encoded_labels = [label_dict[l] for l in labels]
    return encoded_labels

def create_model():
    model = Sequential()
    model.add(Conv1D(filters=32, kernel_size=3, activation='relu', input_shape=(130, 1)))
    model.add(MaxPooling1D(pool_size=2))
    model.add(Conv1D(filters=64, kernel_size=3, activation='relu'))
    model.add(MaxPooling1D(pool_size=2))
    model.add(Flatten())
    model.add(Reshape((-1, 64)))  # 新增的Reshape层
    model.add(LSTM(64))
    model.add(Dense(64, activation='relu'))
    model.add(Dense(len(set(labels)), activation='softmax'))
    model.compile(loss='sparse_categorical_crossentropy', optimizer='adam', metrics=['accuracy'])
    return model

seqs, labels = read_fasta('CP015726.fasta')
seq_num = seq_to_num(seqs)
encoded_labels = encode_labels(labels)
X_train, X_test, y_train, y_test = train_test_split(seq_num, np.array(encoded_labels), test_size=0.2, random_state=42)
model = create_model()
history = model.fit(X_train.reshape((-1, 130, 1)), y_train, epochs=10, batch_size=32, validation_split=0.2)
score = model.evaluate(X_test.reshape((-1, 130, 1)), y_test, verbose=0)
print('Test loss:', score[0])
print('Test accuracy:', score[1])

代码解释:

seq_to_num 函数的修改:
- 原代码中，seq_to_num 函数仅考虑了 ACGT 四种碱基，如果遇到其他碱基，例如 'N'，就会报错。
- 修改后的代码中，我们将未知碱基 'N' 替换为 0，保证代码可以正常运行。
其他部分:
- read_fasta 函数用于读取 FASTA 格式的 DNA 序列文件。
- encode_labels 函数将标签转换为数字编码。
- create_model 函数创建 LSTM 模型。
- train_test_split 函数将数据分成训练集和测试集。
- model.fit 函数训练模型。
- model.evaluate 函数评估模型性能。

注意:

请根据实际情况修改 read_fasta 函数中的文件路径。
您可以调整模型参数，例如 epochs、batch_size 和 filters 等，以提高模型性能。
本代码示例仅供参考，您可以根据具体需求进行修改。