VIDHOP, viral host prediction with Deep Learning 论文阅读笔记

VIDHOP, viral host prediction with Deep Learning 论文阅读笔记

github : https://github.com/flomock/vidhop

摘要

Zoonosis即人类和动物都可相互传染致病，如寨卡病毒、埃博拉病毒和新冠病毒等，为了预防全球化带来的病毒传染加快的问题，本文提出一种基于病毒的基因组序列来推测病毒宿主的预测方法（input = 病毒基因碱基序列， Y = 宿主种类），并且定义了一种基于预测宿主数来计算的平均准确度的计算公式，本模型可以用于transfer到其他的病毒上面去，直接分类准确度起伏较大，基于作者定义的准确度比较平稳。

原理

模型总览

如上图所示，先从国家某些生物网站上下载数据——本文中从欧洲核苷酸档案 (ENA) 数据库 收集了带有宿主标签的甲型流感病毒、狂犬病狂犬病病毒和轮状病毒 A 的所有核苷酸序列，再从德国生物技术信息中心 (NCBI) 提供的分类信息来管理宿主标签——然后将其按照6：2：2来划分成训练集、验证集和测试集，在将数据进行预处理输入模型得到预测结果并且进行分析，下文将对模型每一个步骤进行深入的介绍。

数据预处理

数据预处理的要点在于每一个病毒基因序列的长短都是不一样的，当数据集一旦变大，这个现象带来的影响就被放大。
所以就要用到上图所示的方法对数据集进行修剪，修剪的长度等于基因从小到大排序的95%分位的长度，这样就可以很好地把数据统一起来。

但是有个问题是，如果最短基因序列的长度过于短，这样修剪出来的数据放到模型里面训练就会欠拟合，所以本文针对这一问题提出了几种数据扩展方案：

1. 简单重复：将原基因本体重复拼接在后面以达到规定长度
2. 简单重复并填充占位符：在上一个方法的基础上在任意位置随机插入2-8个占位符
3. 随机重复：随机抽取原基因的一部分重复拼接到原基因后面
4. 随机重复并填充占位符：同上类比
5. 填充占位符：在原基因的末尾全部填充占位符到规定长度
6. 裁剪：将全部基因的长度修剪到最短基因的长度

处理好之后的数据集在输入到模型的时候为了加速训练的过程，将其拆分成小片分批依次输入到模型里面去。最后要对多批次的输出结果进行处理得到一个最终的预测结果。对碱基进行one-hot编码，如：A = [1, 0, 0, 0, 0], T = [0, 0, 0, 1, 0], N = [0, 0, 0, 0, 1]。

本文还提出了一种叫做在线学习(Online)的方法，这种方法并不需要在输入模型训练的时候对数据进行预处理，而是在训练的过程中用修改过的数据对模型施加一定的影响，具体操作方法见实际代码，文中并没有加以详细描述。

在训练的时候，可以划分好一定数量的验证集和测试集，基于剩下的数据，每一个训练epoch都抽取一定数量（小于剩下样本的总数）不同的样本进行训练，这叫做随机重复欠采样，可以避免因为某个类别的数据过多引起的训练“注意力偏移”。

DNN模型结构

本文着重提出了两种DNN模型，一种是纯LSTM模型，另一种是CNN+LSTM模型，他们的差别在于前者是由三层LSTM接两层Dense，后者是两层全连接层+两层LSTM+两层Dense，最后一层Dense用于分类。LSTM可以很好的应对复杂的数据集和处理任务，CNN则相对于LSTM来说计算速度是其四倍。对于类别种类很多的模型来说，数据量很大，所以就更加需要之前提到的输入切片。

预测结果处理

由于输入用的是切片输入，一个长的基因被切成好几段长分别输入到模型里面，每一段都会有相应的预测结果，所以要把全部结果整合到一起去才能代表整个一个基因的预测结果，所以本文提出了几种方法：

1. 不处理：直接输出每个基因切片的预测结果
2. 票选法： 对所有子序列使用多数“投票”来确定预测结果，每个子序列的票就是其所属类别。这种方法得出的结果偏差较大比较“离散”，不具有一般连续性。
3. 均值法：将所有子序列的概率按类取平均值，最后得出平均值大的那一类作为最终结果。
4. 标准差：在上以种方法的基础上，用其标准差对每个子序列进行加权。 具有更多不同预测类别的子序列获得更高的权重。
   最后的结果要通过以上方法处理之后才能得到。

实验结果

概况

实验用了两种模型针对三个不同的数据集分别进行，最多训练了1500个epochs，当acc在300个epochs还没上升的情况下，模型停止训练。对于轮状病毒A数据集来说，数据集总共有40000个病毒基因序列，对应6个宿主类别。因为这几个宿主类别之间关联性比较强，实验测得6个不同的宿主导致预期的随机准确度约为 16.67%，导致两个模型直接测得的预测准确率都比较高，分别是85.28%和82.88%，而在甲流病毒数据集中宿主类别有36个，36个不同的宿主导致预期的随机准确度只有2.78%，所以直接测得的准确率只有50%左右。所以可以得出一个结论：宿主数量越多类型差别越大，随机准确率越低。

实验结果如下表所示

不难看出，在诸多预处理数据扩展方案中，最有效的是简单重复并填充占位符和随机重复并填充占位符，其他的方案效果不佳。

为了通过考虑类的数量获得更好的可比性，所以本文引入了一个新的准确率计算公式：
$$\text{average accuracy}=\frac{2\cdot \text{ accuracy }+\mid \text{ classes }\mid -2}{\mid \text{ classes }\mid }$$
通过这个公式就可以修正分类准确率，如下表示
与其他模型的对比请看原文

总结

本文提出了一种使用病毒碱基序列编码成数据用于训练的模型来预测分类宿主，并用了几个不同的数据集进行实验，结果表明当分类的类别增大的时候，想要准确预测宿主类别就变得比较困难，本文中使用的模型结构比较简单，并没有用上现在最新的模型。本文的主要贡献其一也就是本文花大量篇幅来讲解的数据处理部分，在对基因序列进行修剪之后怎么样对它们进行处理来让他们变成可用的数据；其二就是为我们引出了一个新的方向，关于使用NLP的方法来处理基因序列，从而可以训练模型来预测给定基因对应的宿主类型，在现在新冠流行的大环境下是一个比较有意义的事情。

2021.07.12更新

实验流程及对应代码解读

代码结构

vidhop
|-- DataParsing
|   `-- DataParsing_main.py
|-- cli.py
|-- training
|   |-- DataGenerator.py
|   |-- make_dataset_out
|   |   |-- X_test.csv
|   |   |-- X_train.csv
|   |   |-- X_val.csv
|   |   |-- Y_test.csv
|   |   |-- Y_train.csv
|   |   `-- Y_val.csv
|   |-- make_datasets.py
|   `-- train_new_model.py
|-- vidhop_main.py
`-- weights
    |-- influ_weights.best.acc.normal_repeat_spacer_run2.hdf5
    |-- rabies_weights.best.acc.random_repeat_run2_design_7.hdf5
    `-- rota_weights.best.acc.online_design_7.hdf5

整个项目的结构如上图所示，其中：

1. DataParsing 部分用于预处理数据
2. DataGenerator.py 用于提供训练模型的数据
3. make_dataset.py 用于从数据集中抽取并划分train、val和test的数据，并写入csv文件，生成用于预处理的数据
4. train_new_model.py 用于训练模型，里面包含了训练模型、处理数据等全部过程
5. vidhop_main.py 用于加载模型来进行测试
6. weights 文件夹内保存的是这个模型的预训练权重

数据预处理——DataParsing.py

class CircularList(list):       # 获取list
    def __getitem__(self, x):
        if isinstance(x, slice):
            return [self[x] for x in self._rangeify(x)]

        index = operator.index(x)
        try:
            return super().__getitem__(index % len(self))
        except ZeroDivisionError:
            raise IndexError('list index out of range')

    def _rangeify(self, slice):
        start, stop, step = slice.start, slice.stop, slice.step
        if start is None:
            start = 0
        if stop is None:
            stop = len(self)
        if step is None:
            step = 1
        return range(start, stop, step)

获取待处理的list

def encode_string(maxLen=None, x=[], y=[], y_encoder=None, repeat=True, use_spacer=False, online_Xtrain_set=False,
                  randomrepeat=False):
    """
    One hot encoding for classes
    to convert the "old" exported int data via OHE to binary matrix
    http://machinelearningmastery.com/multi-class-classification-tutorial-keras-deep-learning-library/

    for dna ony to int values
    """

    def pad_n_repeat_sequences(sequences, maxlen=None, dtype='int32',
                               padding='post', truncating='post', value=0.):
        """extended version of pad_sequences()"""
        if not hasattr(sequences, '__len__'):
            raise ValueError('`sequences` must be iterable.')
        lengths = []
        for x in sequences:
            if not hasattr(x, '__len__'):
                raise ValueError('`sequences` must be a list of iterables. '
                                 'Found non-iterable: ' + str(x))
            lengths.append(len(x))
        num_samples = len(sequences)
        if maxlen is None:
            maxlen = np.max(lengths)  # sequences是基因序列，x是每一个序列的长度，求出最大值maxLen得到填充值

        # take the sample shape from the first non empty sequence
        # checking for consistency in the main loop below.
        sample_shape = tuple()
        for s in sequences:
            if len(s) > 0:
                sample_shape = np.asarray(s).shape[1:]
                break

        # make new array and fill with input seqs
        x = (np.ones((num_samples, maxlen) + sample_shape) * value).astype(dtype)    # np.ones((2, 1) + (1, 2)) 的 shape(2, 1, 1, 2)？？   maybe 3维
        for idx, s in enumerate(sequences):
            if not len(s):
                continue  # empty list/array was found
            if truncating == 'pre':    # 加在序列的前面
                trunc = s[-maxlen:]
            elif truncating == 'post':  # 加在序列的后面
                trunc = s[:maxlen]
            else:
                raise ValueError('Truncating type "%s" not understood' % truncating)

            # check `trunc` has expected shape
            trunc = np.asarray(trunc, dtype=dtype)
            if trunc.shape[1:] != sample_shape:
                raise ValueError(
                    'Shape of sample %s of sequence at position %s is different from expected shape %s' %
                    (trunc.shape[1:], idx, sample_shape))

            if repeat:
                # repeat seq multiple times
                repeat_seq = np.array([], dtype=dtype)
                while len(repeat_seq) < maxLen:
                    if use_spacer:
                        spacer_length = random.randint(1, 50)
                        spacer = [value for i in range(spacer_length)]
                        repeat_seq = np.append(repeat_seq, spacer)
                        if randomrepeat:
                            random_start = random.randint(0, len(trunc))
                            repeat_seq = np.append(repeat_seq,
                                                   CircularList(trunc)[random_start:random_start + len(trunc)])  
                                                   ## 序列+间隔+序列
                            # 随机位置插入一段序列
                        else:
                            repeat_seq = np.append(repeat_seq, trunc)
                    else:
                        if randomrepeat:
                            random_start = random.randint(0, len(trunc))
                            repeat_seq = np.append(repeat_seq,
                                                   CircularList(trunc)[random_start:random_start + len(trunc)])
                        else:
                            repeat_seq = np.append(repeat_seq, trunc)
                x[idx, :] = repeat_seq[-maxLen:]

            else:
                if padding == 'post':
                    x[idx, :len(trunc)] = trunc
                elif padding == 'pre':
                    x[idx, -len(trunc):] = trunc
                else:
                    raise ValueError('Padding type "%s" not understood' % padding)

        return x
 #  ↑ 数据处理部分

    encoder = LabelEncoder()

    if len(x) > 0:
        a = "ATGCN-"

        encoder.fit(list(a))  # fit将a中的6个元素编码成0-5的数字
        out = []
        if type(x)==str:
            dnaSeq = re.sub(r"[^ACGTUacgtu]", 'N', x)
            encoded_X = encoder.transform(list(dnaSeq)) # transform将原始序列变成编码序列
            out.append(encoded_X)
        else:
            for i in x:
                dnaSeq = re.sub(r"[^ACGTUacgtu]", 'N', i)
                # dnaSeq = i[0]
                encoded_X = encoder.transform(list(dnaSeq))
                out.append(encoded_X)

        if online_Xtrain_set:
            X_train_categorial = []
            for seq in out:
                X_train_categorial.append(np.array(to_categorical(seq, num_classes=len(a)), dtype=np.bool))
            return X_train_categorial
        else:
            out = pad_n_repeat_sequences(out, maxlen=maxLen, dtype='int16', truncating='pre', value=0)

        return np.array(to_categorical(out, num_classes=len(a)), dtype=np.bool)
    else:
        if y_encoder != None:
            encoder.fit(y)
            if np.array(encoder.classes_ != y_encoder.classes_).all():
                warning(f"Warning not same classes in training and test set")
            useable_classes = set(encoder.classes_).intersection(y_encoder.classes_)    #  将X和Y放在一起
            try:
                assert np.array(encoder.classes_ == y_encoder.classes_).all()
            except AssertionError:
                warning(
                    f"not all test classes in training data, only