本文將討論深度學習如何幫助細胞生物學捕捉細胞群的多樣性和複雜性。
單細胞RNA測序(scRNAseq)幾年前通過為研究細胞群異質性提供了前所未有的解決方案,徹底改變了生命科學。影響如此巨大,科學雜誌宣布scRNAseq技術成為 2018年的突破。主要進展是認識到儘管生物細胞在顯微鏡下看起來可能在形態學上相似,但它們表達的基因意義上可能非常不同,這反過來導致細胞之間的特徵差異。為了捕獲這種細胞多樣性, Human Cell Atlas 社區宣布了一項雄心勃勃的目標,即建立人體內數萬億細胞的綜合圖譜。
隨著10X Genomics單細胞基因表達平臺的發展,從幾十萬甚至數百萬個細胞中獲得全部轉錄組信息幾乎成為常規操作。因此,scRNAseq代表了目前真正的大數據,它具有卓越的統計功能,並為單細胞數據分析應用機器和深度學習開闢了新的視野。
本文將簡要概述該領域,制定主要的分析挑戰,並展示如何在單細胞RNA測序數據分析中使用Keras和TensorFlow的深度學習來解決無監督學習問題。
為什麼單細胞生物學是深度學習的理想選擇?
對一些數據進行統計分析,我們通常必須了解
特徵數量(基因,蛋白質,遺傳變異,圖像像素等)觀察數量n(樣本,細胞,序列)之間的平衡等等)
人類基因組具有大約p = 20K蛋白質編碼基因,而最近公布的10X scRNAseq數據集包括n~1.3M和n~2M個體細胞。這意味著scRNAseq的n >> p是典型的深度學習限制,即在此極限下工作,我們可以遠遠超出基於線性代數的分析,並捕獲scRNAseq數據中的高度非線性結構。
最近的研究報告了前所未有的scRNAseq樣本量,這是深度學習的理想設置。
對於其他限制,即當n << p或n~p時,貝葉斯和Frequentist統計框架更合適。
換句話說,使用scRNAseq數據我們有一個問題:雖然過擬合和缺乏普遍性是計算生物學中常見的問題,但對於scRNAseq,我們必須注意欠擬合,即如何使用大部分數據。
不僅scRNAseq目前在生命科學領域蓬勃發展,而且在單一細胞水平上提供其他類型信息(的技術變得越來越普遍。研究scATACseq的最新進展導致具有> 100K單細胞的數據集。雖然單獨的scATACseq可能無法保證發現稀有細胞群的新方法,但它提供了與scRNAseq整合的巨大潛力,從而提高了將細胞分配給特定群體的準確性。
染色質可及性區域(ATACseq)是scRNAseq分析的補充
最後,單細胞的多組學技術(CITE-SEQ,scNMTseq等),即來自同一個生物細胞多個信息源,還沒有達到很大的樣本量典型scRNAseq但對未來非常有前途的深層數據集成挑戰學習。
通過深度學習減少維度
由於scRNAseq分析的主要目標是發現新的細胞群,因此它是機器學習術語中的無監督分析。因此,用於scRNAseq的兩種最重要的分析技術是降維和聚類。
Autoencoder是一種無監督的人工神經網絡(ANN),具有有意思的的"蝴蝶"結構,通常用於減少維數。與線性技術(如主成分分析(PCA),多維尺度(MDS),因子分析(FA)等相比,自動編碼器執行非線性降維,因此可捕獲單細胞數據的高度非線性結構。
Autoencoder是一種具有"蝴蝶"結構的人工神經網絡(ANN)
在這裡,我將使用~8K臍帶血單核細胞(CBMCs))CITEseq scRNAseq數據集作為示例,證明線性和非線性自動編碼器維數減少技術之間單細胞解析度的差異。請注意,下面的代碼假設輸入文件採用表格格式,基因為列,單元格為行,文件的最後一列必須是使用你最喜歡的scRNAseq聚類技術獲得的單元格注釋。建議使用基於圖形的聚類與Louvaine社區檢測,這對於高維數據很有用,在流行的Seurat scRNAseq工作流程中實現。
import numpy as np import pandas as pd from keras.layers import Dense import matplotlib.pyplot as plt from sklearn.manifold import TSNE from keras.optimizers import Adam from sklearn.decomposition import PCA from keras.models import Sequential, Model # READ AND LOG-TRANSFORM DATA expr = pd.read_csv('MouseBrain_10X_1.3M.txt',sep=' ') X = expr.values[:,0:(expr.shape[1]-1)] Y = expr.values[:,expr.shape[1]-1] X = np.log(X + 1) # REDUCE DIMENSIONS WITH PRINCIPAL COMPONENT ANALYSIS (PCA) n_input = 50 x_train = PCA(n_components = n_input).fit_transform(X); y_train = Y plt.scatter(x_train[:, 0], x_train[:, 1], c = y_train, cmap = 'tab20', s = 10) plt.title('Principal Component Analysis (PCA)') plt.xlabel("PC1") plt.ylabel("PC2") # REDUCE DIMENSIONS WITH AUTOENCODER model = Sequential() model.add(Dense(30, activation='elu', input_shape=(n_input,))) model.add(Dense(20, activation='elu')) model.add(Dense(10, activation='elu')) model.add(Dense(2, activation='linear', name="bottleneck")) model.add(Dense(10, activation='elu')) model.add(Dense(20, activation='elu')) model.add(Dense(30, activation='elu')) model.add(Dense(n_input, activation='sigmoid')) model.compile(loss = 'mean_squared_error', optimizer = Adam()) model.fit(x_train, x_train, batch_size = 128, epochs = 500, verbose = 1) encoder = Model(model.input, model.get_layer('bottleneck').output) bottleneck_representation = encoder.predict(x_train) plt.scatter(bottleneck_representation[:,0], bottleneck_representation[:,1], c = y_train, s = 10, cmap = 'tab20') plt.title('Autoencoder: 8 Layers') plt.xlabel("Dimension 1") plt.ylabel("Dimension 2")
為了比較,還要在這裡添加t分布式隨機鄰域嵌入(tSNE)圖,這是scRNAseq區域中當前的黃金標準非線性降維技術。tSNE的一個問題是它無法處理高維數據,例如scRNAseq。因此,通常的做法是執行PCA(線性)作為預維數減少並將輸出饋送到tSNE。但是,我們可以通過使用自動編碼器以非線性方式執行預維數減少步驟來做得更好。讓我們顯示兩種策略的tSNE圖:
# TSNE ON PCA model_tsne = TSNE(learning_rate = 200, n_components = 2, random_state = 123, perplexity = 90, n_iter = 1000, verbose = 1) tsne = model_tsne.fit_transform(x_train) plt.scatter(tsne[:, 0], tsne[:, 1], c = y_train, cmap = 'tab20', s = 10) plt.title('tSNE on PCA') plt.xlabel("tSNE1") plt.ylabel("tSNE2") # TSNE ON AUTOENCODER model = Sequential() model.add(Dense(10, activation = 'elu', input_shape=(X.shape[1],))) model.add(Dense(8, activation = 'elu')) model.add(Dense(6, activation = 'elu')) model.add(Dense(4, activation = 'linear', name = "bottleneck")) model.add(Dense(6, activation = 'elu')) model.add(Dense(8, activation = 'elu')) model.add(Dense(10, activation = 'elu')) model.add(Dense(X.shape[1], activation = 'sigmoid')) model.compile(loss = 'mean_squared_error', optimizer = Adam()) model.fit(X, X, batch_size = 128, epochs = 100, shuffle = True, verbose = 1) encoder = Model(model.input, model.get_layer('bottleneck').output) bottleneck_representation = encoder.predict(X) 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, cmap = 'tab20', s = 10) plt.title('tSNE on Autoencoder: 8 Layers') plt.xlabel("tSNE1") plt.ylabel("tSNE2")
這裡一個點是一個單元格,顏色對應不同的單元格類型。你應該觀察12個細胞群,但你基本上無法從PCA圖中看到它們(細胞嚴重重疊),因為線性維數降低無法解析單個細胞結構。自動編碼器圖像看起來更好,可以清楚地檢測到不同的細胞群。tSNE通常提供更銳利的細胞團。
然而,對於這種特殊情況,自動編碼器上的tSNE似乎提供了更密集和透明的簇,特別是對於在PCA圖上的tSNE中塗抹在整個藍色簇中的紫色細胞群。因此,深度學習有望提高檢測新細胞群的解析度。
尋找可縮放的維度減少
除了難以處理高維數據外,當細胞數量達到數十萬和數百萬時,tSNE的擴展性很差。FItSNE是Barnes-Hut tSNE的一種新的有希望的修改,它似乎可以更好地擴展到大量數據。然而,在1.3M小鼠腦細胞上運行FItSNE ,在將其安裝到內存中時遇到了麻煩。特別是,為了檢查數據的全局結構,我們想獲得高難度的tSNE圖。但是,350是我能夠在計算機集群上使用256GB RAM的節點達到的最大。運行FItSNE非常簡單,類似於它們在R中執行tSNE的方式(同樣,Cluster列包含單元格注釋):
library("data.table") source("fast_tsne.R") expr <- suppressWarnings(as.data.frame(fread("10X_Mouse_Brain_1.3M.txt",sep=" "))) my_color<-as.numeric(as.character(expr$Cluster)) expr$Cluster<-NULL expr<-as.data.frame(t(expr)) N_cells<-dim(expr)[2] print(paste0("DATASET CONTAINS ",dim(expr)[1]," GENES AND ",N_cells," CELLS")) tsne_opt_perp <- fftRtsne(t(log10(expr+20)),check_duplicates=FALSE, perplexity=350,dims=2,max_iter=5000) plot(tsne_opt_perp$Y,col=my_color,xlab="tSNE1",ylab="tSNE2",cex=0.5)
均勻流形逼近和投影(UMAP)是另一種非常有趣的非線性降維技術,目前在許多方面似乎都優於tSNE。它比tSNE快,與FItSNE一樣快,但不需要那麼多的內存,它似乎捕獲了scRNAseq數據的局部和全局結構。
from umap import UMAP model = UMAP(n_neighbors = 30, min_dist = 0.3, n_components = 2) umap = model.fit_transform(X_reduced) umap_coords = pd.DataFrame({'UMAP1':umap[:, 0], 'UMAP2':umap[:, 1]}) umap_coords.to_csv('umap_coords_10X_1.3M_MouseBrain.txt', sep=' ') plt.scatter(umap[:, 0], umap[:, 1], c = Y, cmap = 'tab20', s = 1) plt.title('UMAP') plt.xlabel("UMAP1") plt.ylabel("UMAP2")
最近發布了一些基於變分自動編碼器的有趣方法。其中之一是SCVIS,它是一種神經網絡,捕獲並可視化單細胞基因表達數據中的低維結構,此外還保留了局部和全局相鄰結構。要在1.3M小鼠腦細胞上運行SCVIS,我使用以下命令行:
python scvis train --data_matrix_file 10X_1.3M.txt --out_dir out_scvis --data_label_file 10X_1.3M_MouseBrain_CellAnnotationSeurat.txt --verbose --verbose_interval 50 --show_plot
下面我提供了使用來自10X Genomics的1.3M小鼠腦細胞的4種提到的降維技術的比較,即PCA,tSNE / FItSNE,UMAP和SCVIS:
至於上述CITEseq的情況,我們可以看到,與PCA相比,非線性降維技術(SCVIS,UMAP和tSNE / FItSNE)能夠解析scRNAseq數據中的所有細胞群。在這三種技術中,UMAP是最快的,並且提供了相當好的數據低維表示。為了更精確地計算時間,SCVIS花費了大約6小時,FItSNE花了3小時和大量內存,UMAP花了大約3小時。
考慮到scRNAseq數據的增長,本文預測UMAP和Autoencoders將來會取代tSNE。
具有Keras的scRNAseq的深度自動編碼器
最後,在這裡,本文將演示如何從頭開始實現並使用Keras運行Deep Autoencoder,並不難。為了避免將整個數據集加載到內存中的困難,我選擇了前19個主要成分,我通過重新取樣發現了重要成分,即改組基因表達矩陣並檢查由置換矩陣(零假設)解釋的方差百分比。Autoencoder逐步將維度從19減少到2(自動編碼器的瓶頸),每個隱藏層都減少了一個維度。
import numpy as np import pandas as pd from keras.layers import Dense import matplotlib.pyplot as plt from keras.optimizers import Adam from sklearn.decomposition import PCA from keras.models import Sequential, Model from sklearn.preprocessing import MinMaxScaler # READ DATA AND LOG-TRANSFORM DATA expr = pd.read_csv('MouseBrain_10X_1.3M.txt', sep = ' ', header = None) X = expr.values[:,0:(expr.shape[1]-1)] Y = expr.values[:,expr.shape[1]-1] X = np.float32( np.log(X + 1) ) # REDUCE DIMENSIONS WITH PRINCIPAL COMPONENT ANALYSIS (PCA) n_input = 19 x_train = PCA(n_components = n_input).fit_transform(X) y_train = Y x_train = MinMaxScaler().fit_transform(x_train) # REDUCE DIMENSIONS WITH AUTOENCODER model = Sequential() model.add(Dense(18, activation='elu', kernel_initializer='he_uniform', input_shape=(n_input,))) model.add(Dense(17, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(16, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(15, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(14, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(13, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(12, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(11, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(10, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(9, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(8, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(7, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(6, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(5, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(4, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(3, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(2, activation='linear', kernel_initializer='he_uniform', name="bottleneck")) model.add(Dense(3, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(4, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(5, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(6, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(7, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(8, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(9, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(10, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(11, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(12, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(13, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(14, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(15, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(16, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(17, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(18, activation='elu', kernel_initializer='he_uniform')) model.add(Dense(n_input, activation='sigmoid')) model.compile(loss = 'mean_squared_error', optimizer = Adam(lr = 0.0001)) model.summary() # FIT AUTOENCODER MODEL history = model.fit(x_train, x_train, batch_size = 4096, epochs = 100, shuffle = False, verbose = 0) print("" + "Training Loss: ", history.history['loss'][-1]) plt.figure(figsize=(20, 15)) plt.plot(history.history['loss']) plt.title('Model Loss') plt.ylabel('Loss') plt.xlabel('Epoch') plt.show() encoder = Model(model.input, model.get_layer('bottleneck').output) bottleneck_representation = encoder.predict(x_train) # PLOT DIMENSIONALITY REDUCTION plt.figure(figsize=(20, 15)) plt.scatter(bottleneck_representation[:,0], bottleneck_representation[:,1], c = Y, s = 1, cmap = 'tab20') plt.title('Autoencoder: 34 Hidden Layers, 10X Genomics 1.3M Mouse Brain Cells') plt.xlabel("Dimension 1") plt.ylabel("Dimension 2") plt.show() view raw
Deep Autoencoder用於降低1.3M 10X小鼠腦細胞的維數,實現Keras
我們可以看到細胞群是完全可區分的,儘管我沒有專門尋找可能提供更好解析度的神經網絡的最佳配置。Deep Autoencoder的一個巨大優勢是它看起來非常快,因此可以擴展到大量的scRNAseq數據,在筆記本電腦只花了幾分鐘就可以使模型收斂並獲得上面的圖。
具有TensorFlow的scRNAseq的深度自動編碼器
Keras很棒,快速而且容易,但有時仍然會覺得使用TensorFlow可以更好地控制神經網絡。例如,對於Keras,人們永遠不會感覺到"手動"連接節點,這是本文對TensorFlow所理解的理解深度。後者的一個小缺點是,小批量學習這樣的有用技巧必須在TensorFlow中手動編碼,而對比則自動包含在Keras中。這裡是TensorFlow Deep Autoencoder的代碼和維數減少圖:
import numpy as np import pandas as pd import tensorflow as tf import matplotlib.pyplot as plt X = x_train # DEFINE HYPERPARAMETERS learning_rate = 0.001 training_epochs = 500 mini_batch_size = 1024 # mini_batch_size = X.shape[0]-1 display_step = 10 # how often to display loss and accuracy num_hidden_1 = 12 # 1st hidden layer num features num_hidden_2 = 8 # 2nd hidden layer num features num_hidden_3 = 4 # 3-d hidden layer num features num_bottleneck = 2 # bottleneck num features num_input = X.shape[1] # scRANAseq data input (number of genes) # TENSORFLOW GRAPH INPUT x = tf.placeholder("float") y = tf.placeholder("float") weights = { 'encoder_h1': tf.Variable(tf.random_normal([num_input, num_hidden_1])), 'encoder_h2': tf.Variable(tf.random_normal([num_hidden_1, num_hidden_2])), 'encoder_h3': tf.Variable(tf.random_normal([num_hidden_2, num_hidden_3])), 'bottleneck': tf.Variable(tf.random_normal([num_hidden_3, num_bottleneck])), 'decoder_h1': tf.Variable(tf.random_normal([num_bottleneck, num_hidden_3])), 'decoder_h2': tf.Variable(tf.random_normal([num_hidden_3, num_hidden_2])), 'decoder_h3': tf.Variable(tf.random_normal([num_hidden_2, num_hidden_1])), 'decoder_out': tf.Variable(tf.random_normal([num_hidden_1, num_input])), } biases = { 'encoder_b1': tf.Variable(tf.random_normal([num_hidden_1])), 'encoder_b2': tf.Variable(tf.random_normal([num_hidden_2])), 'encoder_b3': tf.Variable(tf.random_normal([num_hidden_3])), 'bottleneck': tf.Variable(tf.random_normal([num_bottleneck])), 'decoder_b1': tf.Variable(tf.random_normal([num_hidden_3])), 'decoder_b2': tf.Variable(tf.random_normal([num_hidden_2])), 'decoder_b3': tf.Variable(tf.random_normal([num_hidden_1])), 'decoder_out': tf.Variable(tf.random_normal([num_input])), } # CONSTRUCT AUTOENCODER MODEL print("" + "Constructing Autoencoder Model ..." + "") def encoder(x): layer_1 = tf.nn.elu(tf.add(tf.matmul(x, weights['encoder_h1']), biases['encoder_b1'])) layer_2 = tf.nn.elu(tf.add(tf.matmul(layer_1, weights['encoder_h2']), biases['encoder_b2'])) layer_3 = tf.nn.elu(tf.add(tf.matmul(layer_2, weights['encoder_h3']), biases['encoder_b3'])) bottleneck = tf.add(tf.matmul(layer_3, weights['bottleneck']), biases['bottleneck']) return bottleneck def decoder(x): layer_1 = tf.nn.elu(tf.add(tf.matmul(x, weights['decoder_h1']), biases['decoder_b1'])) layer_2 = tf.nn.elu(tf.add(tf.matmul(layer_1, weights['decoder_h2']), biases['decoder_b2'])) layer_3 = tf.nn.elu(tf.add(tf.matmul(layer_2, weights['decoder_h3']), biases['decoder_b3'])) layer_out = tf.nn.sigmoid(tf.add(tf.matmul(layer_3, weights['decoder_out']), biases['decoder_out'])) return layer_out encoder_op = encoder(x) decoder_op = decoder(encoder_op) y_pred = decoder_op y_true = x cost = tf.reduce_mean(tf.pow(y_true - y_pred, 2)) optimizer = tf.train.AdamOptimizer(learning_rate).minimize(cost) # START TRAINING AUTOENCODER print("" + "Start Training Autoencoder ..." + "") my_cost = [] tf.set_random_seed(12) with tf.Session() as sess: sess.run(tf.global_variables_initializer()) for epoch in range(training_epochs): pos = 0 idx = np.arange(X.shape[0]) my_cost_mini_batch = [] for _ in range(10000000): #print('Mini-batch {0} - {1}'.format(pos,pos + mini_batch_size)) if pos + mini_batch_size >= X.shape[0]: break batch_x = X[idx[range(pos, pos + mini_batch_size)],:] c,_ = sess.run([cost, optimizer], feed_dict={x: batch_x}) my_cost_mini_batch.append(c) pos = pos + mini_batch_size my_cost.append(np.mean(my_cost_mini_batch)) if (epoch + 1) % display_step == 0: print("Epoch:", '%04d' % (epoch + 1), "cost = ", "{:.9f}".format(np.mean(my_cost_mini_batch))) pred = sess.run(encoder(x), feed_dict={x: X}) # PLOT LOSS FUNCTION plt.figure(figsize=(20, 15)) plt.plot(range(training_epochs), my_cost) plt.xlabel("Epoch",fontsize = 20) plt.ylabel("Loss",fontsize = 20, rotation = 1) plt.show() # VISUALIZE AUTOENCODER BOTTLENECK plt.figure(figsize=(20, 15)) plt.scatter(pred[:,0], pred[:,1], c = Y, s = 1, cmap = 'tab20') plt.title('Autoencoder: 8 Hidden Layers', fontsize = 20) plt.xlabel("Dimension 1", fontsize = 20) plt.ylabel("Dimension 2", fontsize = 20) plt.show()
Deep Autoencoder用於降低1.3M 10X小鼠腦細胞的維數,TensorFlow實現
同樣,低維表示看起來很有希望。通過對配置和其他超參數的一些調整,可以獲得更好的解析度。同樣,對於Keras來說,這個TensorFlow自動編碼器實現速度非常快,與FItSNE,UMAP和SCVIS用於生成這種降維圖的時間相比只需幾分鐘。如果由於某種原因需要重新採樣程序進行分析,那麼使用FItSNE,UMAP和SCVIS就不太可行,但使用深度學習應該非常簡單。
總結
在這裡,我們了解到單細胞RNA測序(scRNAseq)正在迅速提高我們對生物細胞特徵多樣性的理解。降維可能是典型scRNAseq分析背後最重要的分析工具。由於scRNAseq技術提供的大量數據,黃金標準降維技術tSNE目前在可擴展性方面遇到困難。
通過Autoencoder和UMAP進行深度學習提供了當前最靈活和可擴展的方法來獲得scRNAseq數據的低維表示,並且很可能在未來取代tSNE。最後,我們學習了如何使用Keras和TensorFlow為巨大的10X Genomics 1.3M小鼠腦細胞scRNAseq數據集實現Deep Autoencoders。