神经网络预测RNA序列中RBP结合位点-RNA序列二分类

一、简介

本文使用源自CLIP-seq的数据集ALKBH5_Baltz2012数据集进行实验，旨在使用神经网络来预测RNA序列中的RNA-蛋白质结合位点，实现RNA序列数据的分类任务。

二、数据预处理

在这里我的原始RNA序列的数据集是直接被分成正类和负类的，所以处理起来比较方便。需要这个数据集的可以留言。

处理思路：

（1）用“T”替换序列中的“U”

（2）序列长度不足501的使用N填充到501；序列长度大于501的只取前501个核苷酸

（3）序列—矩阵的转化

下面上代码

# 读取 .fa 文件，返回序列和标签
def read_seq_graphprot(seq_file, label):
    seq_list = []
    labels = []
    seq = ''
    with open(seq_file, 'r') as fp:
        for line in fp:
            if line[0] == '>':
                name = line[1:-1]
            else:
                seq = line[:-1].upper()
                seq = seq.replace('T', 'U')
                seq_list.append(seq)
                labels.append(label)
    return seq_list, labels
 
def read_data_file(posifile, negafile):
    data = dict()
    seqs1, labels1 = read_seq_graphprot(posifile, label = 1)  # 正样本标签是1
    seqs2, labels2 = read_seq_graphprot(negafile, label = 0)  # 负样本标签是0
    seqs = seqs1 + seqs2
    labels = labels1 + labels2
    # print(labels)
    data["seq"] = seqs
    data["Y"] = np.array(labels)
    return data
 
# 若长度不足501，使用N填充到501；若长度大于501，只取前501个核苷酸
def padding_sequence(seq, max_len = 501, repkey = 'N'):
    seq_len = len(seq)
    if seq_len < max_len:
        gap_len = max_len -seq_len
        new_seq = seq + repkey * gap_len
    else:
        new_seq = seq[:max_len]
    return new_seq

# 把处理好的等长（501）RNA序列 转变成 矩阵
def get_RNA_seq_concolutional_array(seq, motif_len = 4):
    seq = seq.replace('U', 'T')
    alpha = 'ACGT'
    row = (len(seq) + 2*motif_len - 2)
    new_array = np.zeros((row, 4))
    for i in range(motif_len-1):
        new_array[i] = np.array([0.25]*4)

    for i in range(row-3, row):
        new_array[i] = np.array([0.25]*4)
        
    #pdb.set_trace()
    for i, val in enumerate(seq):
        i = i + motif_len-1
        if val not in 'ACGT':
            new_array[i] = np.array([0.25]*4)
            continue
        try:
            index = alpha.index(val)
            new_array[i][index] = 1
        except:
            pdb.set_trace()
        #data[key] = new_array
    return new_array

def get_bag_data_1_channel(data, max_len):
    bags = []
    seqs = data["seq"]
    labels = data["Y"]
    for seq in seqs:
        bag_seq = padding_sequence(seq, max_len = max_len)
        bag_subt = []
        tri_fea = get_RNA_seq_concolutional_array(bag_seq)
        # print tri_fea
        bag_subt.append(tri_fea.T)
        # print tri_fea.T
        bags.append(np.array(bag_subt))
        # print bags
    return bags, labels

# channel
def get_data(posi, nega, window_size):
    # 读取正、负样本文件，处理得到标准化的data和label
    data = read_data_file(posi, nega)
    # 得到处理好的样本矩阵和对应标签
    train_bags, label = get_bag_data_1_channel(data, max_len = window_size)
    return train_bags, label

三、网络模型

利用上面的方法对RNA序列处理成矩阵后，就可以得到RNA的序列矩阵了。这里我统一处理成了单通道，窗口长度是501，最终训练数据的维度是（2170,1,4,507）。处理成其他形式的效果也大同小异。可以看一下这位学者的处理思路。

下面言归正传，网络模型我试了比较多，但是感觉深层的网络反而得不到好的效果，甚至训练出来的网络过拟合或者loss一直0.69（概率都是0.5）。这可能跟RNA序列的特殊性有关，所以试来试去还是用了非常简单的模型，具体结构看下面这个图吧。

这个网络甚至比LeNet还要简单，但是效果确实比一些更深的网络好。下面上代码。

# 简化模型
model = nn.Sequential(
    nn.Conv2d(1, 32, kernel_size=3),
    nn.ReLU(),  # ([1, 32, 2, 505])
    nn.MaxPool2d(kernel_size=2,padding=(1,0)),  # ([1, 32, 2, 252])
    nn.Flatten(),  # ([1, 16*2*126])
    nn.Linear(32 * 2 * 252, 120),
    nn.BatchNorm1d(120),
    nn.ReLU(),
    # nn.Dropout(0.1),  # 添加Dropout层，可加可不加
    nn.Linear(120, 60),
    nn.ReLU(),
    nn.Linear(60, 2)
)

三、参数调整

使用随机搜索方法找到了准确率最高的超参数，最终使用验证集得到准确率是70.4%。用了Adagrad优化方法，学习率是0.01，epochs是10，batchsize是128。

训练过程就不展示了，最后有需要完整代码的可以去这里。