SCALE全稱是Single-Cell ATAC-seq analysis vie Latent feature Extraction, 從名字中就能知道這個軟件是通過隱特征提取的方式分析單細胞ATAC-seq數(shù)據(jù)。
在文章中夯秃,作者從開發(fā)者的角度列出了目前的scATAC-seq分析軟件座咆,chromVAR, scABC, cisTopic, scVI,發(fā)現(xiàn)每個軟件都有一定的不足之處仓洼,而從我們軟件使用者的角度介陶,其實可以考慮都試試這些工具。
SCALE結(jié)合了深度生成模型(Depp Generative Models)變分自動編碼器框架(Variational Autoencoder, VAE)與概率高斯混合模型(Gaussian Mixture Model, GMM)去學習隱特征色建,用于準確地鑒定scATAC-seq數(shù)據(jù)中的特征哺呜。
文章通過一張圖來解釋了軟件的工作機制:
SCALE將sc-ATAC-seq的輸入數(shù)據(jù)x(Cells-by-Peaks矩陣)建模成一個聯(lián)合分布,p(x,z,c)箕戳,c是GMM組件中對應的預定義的K個聚類某残,z是一個隱變量,是細胞在所有peak中實際可能的值陵吸,用于后續(xù)的聚類和可視化玻墅。z通過 計算而得,公式里面的
是編碼器網(wǎng)絡從x中學習而得壮虫,
則是從
抽樣而成澳厢。
從公式中我們還可以發(fā)現(xiàn)z其實和GMM的c有關,所以p(x,z,c)也可以寫成P(x|z)p(z|c)p(c)囚似,而p(c)是K個預定義聚類分布的離散概率分布剩拢,p(z|c)服從混合高斯分布,而p(x|z)則是服從多變量伯努利分布(multivartiable Bernoulli distribution), 通過解碼者網(wǎng)絡建模而成饶唤。
當然從一個軟件使用者的角度而言徐伐,我們不會去關心代碼,也不會關心原理募狂,我們更關心的是這個工具能做什么办素。SCALE能做以下的分析
- SCALE可以對隱特征聚類識別細胞類群
- SCALE可以降噪魏保,恢復缺失的peak
- SCALE能夠區(qū)分批次效應和生物學細胞類群之間的差異
軟件安裝
推薦使用conda的方式進行軟件安裝(我測試過了,運行沒有問題)
第一步:創(chuàng)建一個環(huán)境摸屠,名字就是SCALE
谓罗,并且啟動該環(huán)境
conda create -n SCALE python=3.6 pytorch
conda activate SCALE
第二步:從GitHub上克隆該項目
git clone git://github.com/jsxlei/SCALE.git
第三步:安裝SCALE
cd SCALE
python setup.py install
之后分析的時候,只需要通過conda activate SCALE
就能啟動分析環(huán)境季二。
考慮后續(xù)要交互的讀取數(shù)據(jù)和可視化檩咱,那么建議再安裝一個Jupyter
conda install jupyter
軟件使用
SCALE支持兩類輸入文件:
- count矩陣,行為peak胯舷,列為barcode
- 10X輸出文件: count.mtx.gz, peak.tsv, barcode.tsv
我們以官方提供的Forebrain數(shù)據(jù)集為例進行介紹刻蚯,因為這個數(shù)據(jù)相對于另外一個數(shù)據(jù)集Mouse Atlas小多了。
我們在服務器上新建一個文件夾桑嘶,用于存放從https://cloud.tsinghua.edu.cn/d/bda0332212154163a647/下載的數(shù)據(jù)
mkdir Forebrain
保證Forebrain
有下載好的數(shù)據(jù)
$ ls Forebrain
data.txt
之后運行程序
SCALE.py -d Forebrain/data.txt -k 8 --impute
軟件運行步驟為:
- 加載數(shù)據(jù): Loading data
- 模型訓練: Training Model
- 輸出結(jié)果: Saving imputed data
其中模型訓練這一步時間比較久炊汹,可以嘗試用GPU加速(我是普通CPU服務器沒有辦法)。最終會在當前文件夾看到一個output
文件夾逃顶,里面有如下內(nèi)容:
- imputed_data.txt: 每個細胞在每個特征的推斷值讨便,建議用
--binary
保存二進制格式 - model.pt: 用于重復結(jié)果的模型文件,
--pretrain
參數(shù)能夠讀取該模型 - feature.txt: 每個細胞的隱特征以政,用于聚類和可視化
- cluster_assignments.txt: 兩列霸褒,barcode和所屬類群
- tsne.txt, tsne.pdf: tSNE的坐標和PDF文件,坐標文件可以導入到R語言進行可視化
上面是命令行部分盈蛮,下面則是Python環(huán)境進行交互式操作废菱,輸入jupyter notebook
,之后在網(wǎng)頁上打開
首先是導入各種Python庫
import pandas as pd
import numpy as np
from sklearn.metrics import confusion_matrix
from matplotlib import pyplot as plt
import seaborn as sns
from scale.plot import plot_embedding, plot_heatmap
然后加載分析結(jié)果抖誉,包括聚類信息和特征信息
y = pd.read_csv('output/cluster_assignments.txt', sep='\t', index_col=0, header=None)[1].values
feature = pd.read_csv('output/feature.txt', sep='\t', index_col=0, header=None)
通過熱圖展示不同聚類細胞之間的差異圖
plot_heatmap(feature.T, y,
figsize=(8, 3), cmap='RdBu_r', vmax=8, vmin=-8, center=0,
ylabel='Feature dimension', yticklabels=np.arange(10)+1,
cax_title='Feature value', legend_font=6, ncol=1,
bbox_to_anchor=(1.1, 1.1), position=(0.92, 0.15, .08, .04))
如果要矯正批次效應殊轴,可以通過根據(jù)feature的heatmap,去掉和batch相關的feature來實現(xiàn)
我們可以展示SCALE對原始數(shù)據(jù)糾正后的值(imputed data), 該結(jié)果也能提高chromVAR鑒定motif的效果
imputed = pd.read_csv('output/imputed_data.txt', sep='\t', index_col=0)
展示聚類特異性的peak袒炉, 分析由mat_specificity_score
和cluster_specific
完成
from scale.specifity import cluster_specific, mat_specificity_score
score_mat = mat_specificity_score(imputed, y)
peak_index, peak_labels = cluster_specific(score_mat, np.unique(y), top=200)
plot_heatmap(imputed.iloc[peak_index], y=y, row_labels=peak_labels, ncol=3, cmap='Reds',
vmax=1, row_cluster=False, legend_font=6, cax_title='Peak Value',
figsize=(8, 10), bbox_to_anchor=(0.4, 1.2), position=(0.8, 0.76, 0.1, 0.015))
參數(shù)介紹
通過SCALE.py -h
可以輸出SCALE的所有可用參數(shù)
-
-d/--dataset
: 單個文件矩陣應該指定文件路徑旁理,10X輸出的多個文件則是文件目錄 -
-k
: 設定輸出結(jié)果的聚類數(shù) -
-o
: 輸出文件路徑 -
--pretrain
: 讀取之前訓練的模型 -
--lr
: 修改起始學習速率, 默認是0.002,和模型訓練有關 -
--batch_size
: 批處理大小梳杏, 默認就行韧拒,不需要修改(和批次效應處理無關) -
-g GPU
: 選擇GPU設備數(shù)目淹接,非GPU服務器用不到 -
--seed
: 初始隨機數(shù)種子十性,通常在遇到nan缺失時考慮修改 -
-encode_dim
,-decode_dim
: 編碼器和解碼器的維度,通常也不需要修改 -
-latent
隱藏層維度 -
--low
,--high
: 過濾低質(zhì)量的peak, 即出現(xiàn)比例高于或者低于某個閾值的peak塑悼,默認是0.01和0.9劲适。作者推薦保留1萬-3萬的peak用于SCALE分析。如果數(shù)據(jù)質(zhì)量很高厢蒜,且peak數(shù)不多于10萬霞势,那么可以不過濾烹植。 -
--min_peaks
: 過濾低質(zhì)量細胞,如果該細胞的peak低于閾值 -
log_transform
:log2(x+1)
的變換 -
--max_iter
: 最大迭代數(shù)愕贡,默認是30000, 可以觀察損失收斂的情況來修改草雕,也就是訓練模型這一步輸出的信息 -
-weight_decay
: 沒有說明 -
--impute
: 保存推斷數(shù)據(jù),默認開啟 -
--binary
: 推薦加上該參數(shù)固以,減少imputed data占用空間 -
--no_tsne
: 不需要保存t-SNE結(jié)果 -
--reference
: 參考細胞類型 -
-t
: 如果輸出矩陣是列為peak墩虹,行為barcode,用該參數(shù)進行轉(zhuǎn)置
對于使用者而言憨琳,我們一般只用修改-k
更改最后的聚類數(shù)诫钓,--low
, --high
, ---min_peaks
來對原始數(shù)據(jù)進行過濾,以及加上--binary
節(jié)約空間篙螟。
版權聲明:本博客所有文章除特別聲明外菌湃,均采用 知識共享署名-非商業(yè)性使用-禁止演繹 4.0 國際許可協(xié)議 (CC BY-NC-ND 4.0) 進行許可。