作者 | Nikolay Oskolkov
来源 | Medium
编辑 | 代码医生团队
本文将讨论如何利用多种生物信息源，OMIC数据，以便通过深度学习实现更准确的生物系统建模。
生物和生物医学研究已大大过去十年中从技术进步提供DNA序列（代受益组学），基因表达（转录组学），蛋白质丰度（prote 组学）和生物信息等诸多层面通常被称为组学。尽管个别OMIC层能够回答许多重要的生物学问题，但它们的组合以及随之而来的互补性的协同效应有望对细胞，组织和生物体等生物系统的行为产生新的见解。因此，OMIC的整合代表了生物学和生物医学的当代挑战。
在本文中将使用Keras进行深度学习，并展示如何集成多个OMIC数据，以揭示在各个OMIC中不可见的隐藏模式。
单细胞产生大数据
数据集成的问题对于数据科学来说并不是全新的问题。想象一下知道一个人看某些图像，阅读某些文本并听某些音乐。图像，文本和声音是非常不同类型的数据，但是可以尝试组合这些类型的数据以便构建例如更好的推荐系统，其实现捕获人的兴趣的更高准确度。至于生物学和生物医学，数据整合的概念最近刚到这里，然而它的生物学角度正在积极发展，产生了一些有趣的方法，如mixOmics，MOFA，Similarity Network Fusion（SNF），OnPLS / JIVE / DISCO，贝叶斯网络等
综合OMIC方法
上面列出的所有综合OMIC方法面临的一个问题是维数的诅咒，即无法在具有有限数量的统计观察的高维空间中工作，这是生物数据分析的典型设置。单细胞OMIC技术非常有用，因为它们提供了数十万甚至数百万的统计观察（细胞），因此提供了真正的大数据集成理想。
单细胞多OMIC技术
令人兴奋的是，像CITEseq和scNMTseq这样的多OMIC单细胞技术分别从完全相同的细胞提供两个和三个水平的生物信息。
将CITEseq数据与深度学习集成
将进行单细胞转录（scRNAseq）和蛋白质组学的无监督集成（scProteomics）从CITEseq数据，8个617脐带血单核细胞（CBMC），采用自动编码器，其非常适合用于捕获单细胞组学的高度非线性性质数据。首先从这里下载CITEseq数据，用Pandas和log-transformationing读取它们，这相当于一个温和的规范化。通常行是单元格，列是mRNA或蛋白质特征，最后一列对应于单元格注释。
https://satijalab.org/seurat/v3.0/multimodal_vignette.html
代码语言：javascript复制import numpy as np
import pandas as pd
from keras.models import Model
import matplotlib.pyplot as plt
from sklearn.manifold import TSNE
from keras.utils import plot_model
from keras.layers import Input, Dense
from keras.layers.merge import concatenate
scRNAseq = pd.read_csv('scRNAseq.txt',sep='t')
scProteomics = pd.read_csv('scProteomics.txt',sep='t')
X_scRNAseq = scRNAseq.values[:,0:(scRNAseq.shape[1]-1)]
Y_scRNAseq = scRNAseq.values[:,scRNAseq.shape[1]-1]
X_scProteomics = scProteomics.values[:,0:(scProteomics.shape[1]-1)]
Y_scProteomics = scProteomics.values[:,scProteomics.shape[1]-1]
X_scRNAseq = np.log(X_scRNAseq 1)
X_scProteomics = np.log(X_scProteomics 1)现在将使用Keras功能API构建一个带有4个隐藏层的Autoencoder模型。Autoencoder有两个输入，一个用于每层信息，即scRNAseq和scProteomics，以及相应的两个输出，旨在重建输入。在第二隐藏层中连接两个输入层之前，它们在第一隐藏层中分别线性变换（相当于PCA降维）。最后，合并的OMIC通过Autoencoder的瓶颈进行处理，最后根据Autoencoders 典型的“蝶形”对称，将尺寸逐渐重建为初始尺寸。
无监督的CITEseq数据集成
在下面的自动编码器的代码中，重要的是要注意第一个隐藏层对scRNAseq的977到50个基因进行了严重的维数降低，同时它使得scProteomics几乎不受影响，即将维度从11减少到10.进一步的瓶颈连接后总共60个维度减少到50个潜在变量，这些变量代表mRNA和蛋白质特征的组合。
代码语言：javascript复制# Input Layer
ncol_scRNAseq = X_scRNAseq.shape[1]
input_dim_scRNAseq = Input(shape = (ncol_scRNAseq, ), name = "scRNAseq")
ncol_scProteomics = X_scProteomics.shape[1]
input_dim_scProteomics = Input(shape = (ncol_scProteomics, ), name = "scProteomics")
# Dimensions of Encoder for each OMIC
encoding_dim_scRNAseq = 50
encoding_dim_scProteomics = 10
# Encoder layer for each OMIC
encoded_scRNAseq = Dense(encoding_dim_scRNAseq, activation = 'linear',
name = "Encoder_scRNAseq")(input_dim_scRNAseq)
encoded_scProteomics = Dense(encoding_dim_scProteomics, activation = 'linear',
name = "Encoder_scProteomics")(input_dim_scProteomics)
# Merging Encoder layers from different OMICs
merge = concatenate([encoded_scRNAseq, encoded_scProteomics])
# Bottleneck compression
bottleneck = Dense(50, kernel_initializer = 'uniform', activation = 'linear',
name = "Bottleneck")(merge)
#Inverse merging
merge_inverse = Dense(encoding_dim_scRNAseq encoding_dim_scProteomics,
activation = 'elu', name = "Concatenate_Inverse")(bottleneck)
# Decoder layer for each OMIC
decoded_scRNAseq = Dense(ncol_scRNAseq, activation = 'sigmoid',
name = "Decoder_scRNAseq")(merge_inverse)
decoded_scProteomics = Dense(ncol_scProteomics, activation = 'sigmoid',
name = "Decoder_scProteomics")(merge_inverse)
# Combining Encoder and Decoder into an Autoencoder model
autoencoder = Model(input = [input_dim_scRNAseq, input_dim_scProteomics],
output = [decoded_scRNAseq, decoded_scProteomics])
# Compile Autoencoder
autoencoder.compile(optimizer = 'adam',
loss={'Decoder_scRNAseq': 'mean_squared_error',
'Decoder_scProteomics': 'mean_squared_error'})
autoencoder.summary()可以为来自不同统计分布的OMIC 分配不同的损失函数，例如组合分类和连续数据，可以分别应用分类交叉熵和均方误差。通过自动编码器进行数据集成的另一个好处是，所有OMIC都相互了解，因为每个节点/特征的权重都是通过彼此上下文中的反向传播来更新的。最后训练Autoencoder并将瓶颈输入tSNE进行可视化：
代码语言：javascript复制# Autoencoder training
estimator = autoencoder.fit([X_scRNAseq, X_scProteomics],
[X_scRNAseq, X_scProteomics],
epochs = 100, batch_size = 128,
validation_split = 0.2, shuffle = True, verbose = 1)
print("Training Loss: ",estimator.history['loss'][-1])
print("Validation Loss: ",estimator.history['val_loss'][-1])
plt.plot(estimator.history['loss'])
plt.plot(estimator.history['val_loss'])
plt.title('Model Loss')
plt.ylabel('Loss')
plt.xlabel('Epoch')
plt.legend(['Train','Validation'], loc = 'upper right')
plt.show()
# Encoder model
encoder = Model(input = [input_dim_scRNAseq, input_dim_scProteomics],
output = bottleneck)
bottleneck_representation = encoder.predict([X_scRNAseq, X_scProteomics])
# tSNE on Autoencoder bottleneck representation
model_tsne_auto = TSNE(learning_rate = 200, n_components = 2, random_state = 123,
perplexity = 90, n_iter = 1000, verbose = 1)
tsne_auto = model_tsne_auto.fit_transform(bottleneck_representation)
plt.scatter(tsne_auto[:, 0], tsne_auto[:, 1], c = Y_scRNAseq, cmap = 'tab20', s = 10)
plt.title('tSNE on Autoencoder: Data Integration, CITEseq')
plt.xlabel("tSNE1")
plt.ylabel("tSNE2")
plt.show()
CITEseq数据集成的效果：查看单个OMIC中不可见的模式
将使用单个OMIC获得的tSNE图与结合数据的Autoencoder瓶颈上的tSNE进行比较，可以立即看到集成在某种程度上平均并强化了各个OMIC。例如仅使用scRNAseq数据很难发现紫色簇，因为它与蓝色细胞群不同，但是在整合后，紫色细胞群很容易区分。这是数据集成的力量！
将scNMTseq数据与深度学习集成
虽然CITEseq包括两个单细胞水平的信息（转录组学和蛋白质组学），另一个奇妙的技术scNMTseq提供来自相同生物细胞的三个OMIC：1）转录组学（scRNAseq），2）甲基化模式（scBSseq），和3）开放染色质地区（scATACseq）。
scNMTseq与Autoencoder的数据集成
Autoencoder的体系结构类似于用于CITEseq的体系结构，只有一个特点：在输入层使用Dropout正则化。这是因为我们只有约120个单元序列，而特征空间的维数是数万个，所以我们需要应用正则化来克服维数的诅咒。请注意对于CITEseq来说这不是必需的，有~8K单元和~1K特征，所以情况正好相反。整体scNMTseq并不是一个简单的数据集成案例，这只是单细胞多OMIC时代的开始，很多细胞将很快从这项激动人心的技术中到来，所以最好做好准备。
代码语言：javascript复制import numpy as np
import pandas as pd
from umap import UMAP
from keras.models import Model
import matplotlib.pyplot as plt
from keras.layers.merge import concatenate
from keras.layers import Input, Dense, Dropout
################## READ AND TRANSFORM DATA ##################
scRNAseq = pd.read_csv('scRNAseq.txt',sep='t')
scBSseq = pd.read_csv('scBSseq.txt',sep='t')
scATACseq = pd.read_csv('scATACseq.txt',sep='t')
X_scRNAseq = scRNAseq.values[:,0:(scRNAseq.shape[1]-1)]
Y_scRNAseq = scRNAseq.values[:,scRNAseq.shape[1]-1]
X_scBSseq = scBSseq.values[:,0:(scBSseq.shape[1]-1)]
Y_scBSseq = scBSseq.values[:,scBSseq.shape[1]-1]
X_scATACseq = scATACseq.values[:,0:(scATACseq.shape[1]-1)]
Y_scATACseq = scATACseq.values[:,scATACseq.shape[1]-1]
X_scRNAseq = np.log(X_scRNAseq 1)
X_scBSseq = np.log(X_scBSseq 1)
X_scATACseq = np.log(X_scATACseq 1)
######################## AUTOENCODER ########################
# Input Layer
ncol_scRNAseq = X_scRNAseq.shape[1]
input_dim_scRNAseq = Input(shape = (ncol_scRNAseq, ), name = "scRNAseq")
ncol_scBSseq = X_scBSseq.shape[1]
input_dim_scBSseq = Input(shape = (ncol_scBSseq, ), name = "scBSseq")
ncol_scATACseq = X_scATACseq.shape[1]
input_dim_scATACseq = Input(shape = (ncol_scATACseq, ), name = "scATACseq")
encoding_dim_scRNAseq = 30
encoding_dim_scBSseq = 30
encoding_dim_scATACseq = 30
# Dropout on Input Layer
dropout_scRNAseq = Dropout(0.2, name = "Dropout_scRNAseq")(input_dim_scRNAseq)
dropout_scBSseq = Dropout(0.2, name = "Dropout_scBSseq")(input_dim_scBSseq)
dropout_scATACseq = Dropout(0.2, name = "Dropout_scATACseq")(input_dim_scATACseq)
# Encoder layer for each OMIC
encoded_scRNAseq = Dense(encoding_dim_scRNAseq, activation = 'elu',
name = "Encoder_scRNAseq")(dropout_scRNAseq)
encoded_scBSseq = Dense(encoding_dim_scBSseq, activation = 'elu',
name = "Encoder_scBSseq")(dropout_scBSseq)
encoded_scATACseq = Dense(encoding_dim_scATACseq, activation = 'elu',
name = "Encoder_scATACseq")(dropout_scATACseq)
# Merging Encoder layers from different OMICs
merge = concatenate([encoded_scRNAseq, encoded_scBSseq, encoded_scATACseq])
# Bottleneck compression
bottleneck = Dense(50, kernel_initializer = 'uniform', activation = 'linear',
name = "Bottleneck")(merge)
#Inverse merging
merge_inverse = Dense(encoding_dim_scRNAseq encoding_dim_scBSseq
encoding_dim_scATACseq,
activation = 'elu', name = "Concatenate_Inverse")(bottleneck)
# Decoder layer for each OMIC
decoded_scRNAseq = Dense(ncol_scRNAseq, activation = 'sigmoid',
name = "Decoder_scRNAseq")(merge_inverse)
decoded_scBSseq = Dense(ncol_scBSseq, activation = 'sigmoid',
name = "Decoder_scBSseq")(merge_inverse)
decoded_scATACseq = Dense(ncol_scATACseq, activation = 'sigmoid',
name = "Decoder_scATACseq")(merge_inverse)
# Combining Encoder and Decoder into an Autoencoder model
autoencoder = Model(input = [input_dim_scRNAseq, input_dim_scBSseq,
input_dim_scATACseq],
output = [decoded_scRNAseq, decoded_scBSseq, decoded_scATACseq])
# Compile Autoencoder
autoencoder.compile(optimizer = 'adam',
loss={'Decoder_scRNAseq': 'mean_squared_error',
'Decoder_scBSseq': 'binary_crossentropy',
'Decoder_scATACseq': 'binary_crossentropy'})
autoencoder.summary()
# Autoencoder training
estimator = autoencoder.fit([X_scRNAseq, X_scBSseq, X_scATACseq],
[X_scRNAseq, X_scBSseq, X_scATACseq], epochs = 130,
batch_size = 16, validation_split = 0.2,
shuffle = True, verbose = 1)
print("Training Loss: ",estimator.history['loss'][-1])
print("Validation Loss: ",estimator.history['val_loss'][-1])
plt.plot(estimator.history['loss']); plt.plot(estimator.history['val_loss'])
plt.title('Model Loss'); plt.ylabel('Loss'); plt.xlabel('Epoch')
plt.legend(['Train','Validation'], loc = 'upper right')
# Encoder model
encoder = Model(input = [input_dim_scRNAseq, input_dim_scBSseq,
input_dim_scATACseq], output = bottleneck)
bottleneck_representation = encoder.predict([X_scRNAseq, X_scBSseq, X_scATACseq])
############### UNIFORM MANIFOLD APPROXIMATION AND PROJECTION (UMAP) ###############
model_umap = UMAP(n_neighbors = 11, min_dist = 0.1, n_components = 2)
umap = model_umap.fit_transform(bottleneck_representation)
plt.scatter(umap[:, 0], umap[:, 1], c = Y_scRNAseq, cmap = 'tab10', s = 10)
plt.title('UMAP on Autoencoder: Data Integration, scNMTseq')
plt.xlabel("UMAP1"); plt.ylabel("UMAP2")
将转录组学与scNMTseq的表观遗传学信息相结合
出于好奇将自动编码器的瓶颈结合起来，将三个scNMTseq OMIC组合成均匀流形逼近和投影（UMAP）非线性降维技术，该技术在大量数据的可扩展性方面似乎优于tSNE。当scRNAseq与来自相同细胞（scBSseq和scATACseq）的表观遗传学信息组合时，可以立即看到基因表达意义上的同源蓝色簇分裂成两个簇。因此已经捕获了细胞之间的新异质性当仅查看基因表达scRNAseq数据时隐藏了这一点。这可以成为一种利用生物学的整体复杂性对细胞进行分类的新方法吗？如果是这样，那么问题就出现了：什么是细胞群或细胞类型？
结论
由于最近的技术进步，多种分子和临床信息来源在生物学和生物医学中变得越来越普遍。因此数据整合是合乎逻辑的下一步，它通过利用数据的整体复杂性提供对生物过程的更全面的理解。深度学习框架非常适合数据集成，因为当多种数据类型相互学习信息时，它通过反向传播真正“整合”更新参数。展示了数据集成可以导致数据中新模式的发现，这些模式以前没有在各个数据类型中看到过。
在github上查看这篇文章的代码。
https://github.com/NikolayOskolkov/DeepLearningDataIntegration
推荐阅读
“Keras之父发声：TF 2.0 Keras 深度学习必知的12件事”
