随着深度学习技术的快速发展，利用深度学习技术来研究生物信息学问题的论文也呈现出快速增长的态势(doi: 10.1093/bib/bbw068)。

paper.png

为了继续提高自己的深度学习技术，我选择了“深度学习在生物信息学中的应用”作为这段时间刻意练习的主题。

今天练习的目标是一篇2016年发表在Bioinformatics上的论文：
DeepChrome: deep-learning for predicting gene expression from histone modifications. 虽然论文中的代码可以在Github上查看 ( https://github.com/QData/DeepChrome )，但作者用的是我不熟悉的Torch，没有太大参考价值。我决定用自己熟悉的开发工具来重现作者的工作。

我的刻意练习策略分为以下几个步骤：

1. 研究文章中的算法
2. 下载全部或者部分数据
3. 数据预处理
4. 训练神经网络模型
5. 和文章中的结果进行比较

1. 研究文章中的算法

首先，我们来看看DeepChrome究竟是做什么的？下面是文章的摘要部分，我用黄色highlight出最关键的一句话。

deepchrom_paper_1.png

简单来说，就是用组蛋白修饰（histone modification）数据作为输入，用卷积神经网络（convolutional neural network）来对基因的表达量进行分类（classify）。注意，文章中用的是
classify，而不是
predict。它并不是要预测基因的表达量（回归问题），而是将基因表达量分为高和低两个类别（二元分类问题）。

接下来，看看文章中的算法是怎样的：

deepchrom_paper_2.png

(1) 取基因的转录起始位点（TSS）上下游各5Kb
(2) 将10Kb序列分成长度为100bp的bins
(3) 统计每个bins上5种组蛋白修饰的量
(4) 得到的5×100的矩阵就是基因的表观遗传特征矩阵
(5) 将基因表达量分为高表达（+1）和低表达（-1）两类，得到一个由1和-1组成的向量。分类的标准就是看表达量是否高于基因表达量分布的中值（media）
(6) 用特征矩阵和分类向量做为input，训练卷积神经网络
(7) 神经网络的结构和hyperparameters：

deepchrom_paper_3.png

• 一个temporal convolution layer（filter的数目是Nout，长度是k）
• 激活函数是ReLU
• 用的是Maxpooling（pooling size是m）
• Pooling以后用了一个dropout layer（chosen probability是0.5）
• 用了有2个hiden layer的MLP
• 最后是一个softmax layer

2. 下载全部或者部分数据

表观基因组数据（包括基因表达量数据）来自Integrative Analysis of 111 reference human epigenomes (Nature, Feb. 2015)，下载地址是：http://egg2.wustl.edu/roadmap/web_portal/index.html

DeepChrome文章中使用了56种cell-type的数据，为了节省时间和存储空间，我只下载了6种cell-type (Liver,Lung,Ovary,Pancreas,Thymus,Spleen）的数据。文章只考虑5种组蛋白修饰，因为这5种修饰的数据是56种cell-type都有的。我并不局限于这5种组蛋白修饰，而是将6种cell-type对应的所有组蛋白修饰和DNase数据都下载下来。文章中用read count作为组蛋白修饰的量，我使用的是MACS软件输出的fold-enrichment signal（http://egg2.wustl.edu/roadmap/data/byFileType/signal/consolidated/macs2signal/foldChange/ ）

具体代码见EpiNet（ https://github.com/bioxfu/EpiNet ）：01.download.sh

3. 数据预处理

数据下载后，根据文章中的算法进行预处理：
(1) 根据gene注释找到TSS上下游5Kb，并切成100个bins。用bedtools轻松搞定，具体代码见：02.cut_bins.sh
(2) 根据表观基因组数据（bigwig格式）和bins的location（bed格式），统计每个bins上的组蛋白修饰（或DNase）数据量。用bwtool搞定，具体代码见：03.bins_sum.sh
(3) 合并每个cell-type的bins数据，生成tsv文件。每一行表示一个bin，每一列表示一种组蛋白修饰（或DNase）的量，具体代码见：04.merge_summary.sh
(4) 将gene表达量进行分类：1表示高表达，0表示低表达。具体代码见：05.gene_expr_binary_class.R
(5) 将tsv文件转换成numpy array，并保存在npy文件中。具体代码见：06.tsv_to_npy.sh

4. 训练神经网络模型

Keras（ https://keras.io/ ）是基于ThensorFlow和Theano的一个深度学习开发工具。它的最大好处就是可以快速开发，用它来构建和训练卷积神经网络十分方便。按照它的话说就是：Being able to go from idea to result with the least possible delay is key to doing good research. 我花了一天时间就基本掌握了Keras的API，从此投入Keras的怀抱（后端还是用的TensorFlow，只是用Keras进行快速开发）。

keras.png

(1) 加载所需的packages，并定义卷积神经网络，我给它起名为：EpiNet

code1.png

(2) 加载训练数据。需要注意的是，有些gene的特征矩阵值全部为0，训练的时候会出现问题，所以一开始就把它们去掉。训练前进行normalization可以提高预测的准确率，我用的方法是将特征矩阵中的所有值除以矩阵中的最大值。

code2.png

(3) 将所有数据分成训练集和测试集，并对模型进行初始化

code3.png

(4) 对模型进行训练

code4.png

(5) 对模型进行评估

code5.png

可以看到，用Ovary数据训练的模型的准确率为78%，AUC值为0.86。用同样的方法对其他cell-type的数据进行训练，得到的结果如下：

(‘AUC of Liver:’, 0.84168648346955832)

(‘AUC of Lung:’, 0.85599895550382787)

(‘AUC of Pancreas:’, 0.86048219116434088)

(‘AUC of Spleen:’, 0.84516774396835082)

(‘AUC of Thymus:’, 0.84395698331240487)

5. 和文章中的结果进行比较

这是DeepChrome文章中各个cell-type的AUC值，黄色标注的是我用到的6个cell-type：

deepchrom_auc.png

看上去我的EpiNet要比DeepChrome效果好。其实我们的网络结构基本上是一致的，可能是我的input要更好一些：我用了6到7种组蛋白数据，使用的是fold-enrichment signal而不是read count。由于这只是一个练习，不是为了发表论文，具体原因我就不再深究了。

EpiNet的代码可以在Github上下载 ( https://github.com/bioxfu/EpiNet )：

github.png

欢迎继续关注Rapp带来的【深度学习】系列文章！

qrcode.jpg