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) 進行許可。
