使用SCALE分析單細胞ATAC-seq數(shù)據(jù)

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框架

SCALE將sc-ATAC-seq的輸入數(shù)據(jù)x(Cells-by-Peaks矩陣)建模成一個聯(lián)合分布,p(x,z,c)箕戳,c是GMM組件中對應的預定義的K個聚類某残,z是一個隱變量,是細胞在所有peak中實際可能的值陵吸,用于后續(xù)的聚類和可視化玻墅。z通過z = \mu_z + \sigma_Z \times \epsilon 計算而得,公式里面的 \mu_z \sigma_z 是編碼器網(wǎng)絡從x中學習而得壮虫,\epsilon 則是從 \mathbb{N}(0,\mathbf{I}) 抽樣而成澳厢。

從公式中我們還可以發(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))
heatmap

如果要矯正批次效應殊轴,可以通過根據(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_scorecluster_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))
聚類特異性peak

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

掃碼即刻交流
?著作權歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 序言:七十年代末遍略,一起剝皮案震驚了整個濱河市惧所,隨后出現(xiàn)的幾起案子,更是在濱河造成了極大的恐慌绪杏,老刑警劉巖纯路,帶你破解...
    沈念sama閱讀 222,681評論 6 517
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件,死亡現(xiàn)場離奇詭異寞忿,居然都是意外死亡驰唬,警方通過查閱死者的電腦和手機,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 95,205評論 3 399
  • 文/潘曉璐 我一進店門腔彰,熙熙樓的掌柜王于貴愁眉苦臉地迎上來叫编,“玉大人,你說我怎么就攤上這事霹抛〈暧猓” “怎么了?”我有些...
    開封第一講書人閱讀 169,421評論 0 362
  • 文/不壞的土叔 我叫張陵杯拐,是天一觀的道長霞篡。 經(jīng)常有香客問我,道長端逼,這世上最難降的妖魔是什么朗兵? 我笑而不...
    開封第一講書人閱讀 60,114評論 1 300
  • 正文 為了忘掉前任,我火速辦了婚禮顶滩,結(jié)果婚禮上余掖,老公的妹妹穿的比我還像新娘。我一直安慰自己礁鲁,他們只是感情好盐欺,可當我...
    茶點故事閱讀 69,116評論 6 398
  • 文/花漫 我一把揭開白布赁豆。 她就那樣靜靜地躺著,像睡著了一般冗美。 火紅的嫁衣襯著肌膚如雪魔种。 梳的紋絲不亂的頭發(fā)上,一...
    開封第一講書人閱讀 52,713評論 1 312
  • 那天粉洼,我揣著相機與錄音务嫡,去河邊找鬼。 笑死漆改,一個胖子當著我的面吹牛心铃,可吹牛的內(nèi)容都是我干的。 我是一名探鬼主播挫剑,決...
    沈念sama閱讀 41,170評論 3 422
  • 文/蒼蘭香墨 我猛地睜開眼去扣,長吁一口氣:“原來是場噩夢啊……” “哼!你這毒婦竟也來了樊破?” 一聲冷哼從身側(cè)響起愉棱,我...
    開封第一講書人閱讀 40,116評論 0 277
  • 序言:老撾萬榮一對情侶失蹤,失蹤者是張志新(化名)和其女友劉穎哲戚,沒想到半個月后奔滑,有當?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體,經(jīng)...
    沈念sama閱讀 46,651評論 1 320
  • 正文 獨居荒郊野嶺守林人離奇死亡顺少,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點故事閱讀 38,714評論 3 342
  • 正文 我和宋清朗相戀三年朋其,在試婚紗的時候發(fā)現(xiàn)自己被綠了。 大學時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片脆炎。...
    茶點故事閱讀 40,865評論 1 353
  • 序言:一個原本活蹦亂跳的男人離奇死亡梅猿,死狀恐怖,靈堂內(nèi)的尸體忽然破棺而出秒裕,到底是詐尸還是另有隱情袱蚓,我是刑警寧澤,帶...
    沈念sama閱讀 36,527評論 5 351
  • 正文 年R本政府宣布几蜻,位于F島的核電站喇潘,受9級特大地震影響,放射性物質(zhì)發(fā)生泄漏梭稚。R本人自食惡果不足惜颖低,卻給世界環(huán)境...
    茶點故事閱讀 42,211評論 3 336
  • 文/蒙蒙 一、第九天 我趴在偏房一處隱蔽的房頂上張望哨毁。 院中可真熱鬧枫甲,春花似錦、人聲如沸扼褪。這莊子的主人今日做“春日...
    開封第一講書人閱讀 32,699評論 0 25
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽话浇。三九已至脏毯,卻和暖如春,著一層夾襖步出監(jiān)牢的瞬間幔崖,已是汗流浹背食店。 一陣腳步聲響...
    開封第一講書人閱讀 33,814評論 1 274
  • 我被黑心中介騙來泰國打工, 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留赏寇,地道東北人吉嫩。 一個月前我還...
    沈念sama閱讀 49,299評論 3 379
  • 正文 我出身青樓,卻偏偏與公主長得像嗅定,于是被迫代替她去往敵國和親自娩。 傳聞我的和親對象是個殘疾皇子,可洞房花燭夜當晚...
    茶點故事閱讀 45,870評論 2 361

推薦閱讀更多精彩內(nèi)容