作者:洲更 (上海生科院)博客:http://xuzhougeng.top/archives/Analysis-sc-ATAC-seq-with-SCALE文章链接:https://www.nature.com/articles/s41467-019-12630-7
SCALE全称是Single-Cell ATAC-seq analysis vie Latent feature Extraction, 从名字中就能知道这个软件是通过隐特征提取的方式分析单细胞ATAC-seq数据。
在文章中,作者从开发者的角度列出了目前的scATAC-seq分析软件,chromVAR, scABC, cisTopic, scVI,发现每个软件都有一定的不足之处,而从我们软件使用者的角度,其实可以考虑都试试这些工具。
SCALE结合了深度生成模型(Depp Generative Models
)变分自动编码器框架(Variational Autoencoder, VAE
)与概率高斯混合模型(Gaussian Mixture Model, GMM
)去学习隐特征,用于准确地鉴定scATAC-seq数据中的特征。
文章通过一张图来解释了软件的工作机制:
SCALE框架
SCALE将sc-ATAC-seq的输入数据x(Cells-by-Peaks矩阵)建模成一个联合分布,
p(x,z,c)
,c是GMM组件中对应的预定义的K个聚类,z是一个隐变量,是细胞在所有peak中实际可能的值,用于后续的聚类和可视化。z通过
z=μz+σZ×ϵz=μz+σZ×ϵ 计算而得,公式里面的 μzμz σzσz 是编码器网络从x中学习而得,ϵϵ 则是从 $\mathbb(0,\mathbf)$ 抽样而成。
从公式中我们还可以发现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), 通过解码者网络建模而成。
当然从一个软件使用者的角度而言,我们不会去关心代码,也不会关心原理,我们更关心的是这个工具能做什么。SCALE能做以下的分析
推荐使用conda的方式进行软件安装(我测试过了,运行没有问题)
第一步:创建一个环境,名字就是SCALE
,并且启动该环境
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
就能启动分析环境。
考虑后续要交互的读取数据和可视化,那么建议再安装一个Jupyter
conda install jupyter
SCALE支持两类输入文件:
我们以官方提供的Forebrain数据集为例进行介绍,因为这个数据相对于另外一个数据集Mouse Atlas小多了。
我们在服务器上新建一个文件夹,用于存放从https://cloud.tsinghua.edu.cn/d/bda0332212154163a647/下载的数据
mkdir -p Forebrain
保证Forebrain
有下载好的数据
$ ls Forebrain
data.txt
之后运行程序
SCALE.py -d Forebrain/data.txt -k 8 --impute
软件运行步骤为:
其中模型训练这一步时间比较久,可以尝试用GPU加速(我是普通CPU服务器没有办法)。最终会在当前文件夹看到一个output
文件夹,里面有如下内容:
--binary
保存二进制格式--pretrain
参数能够读取该模型上面是命令行部分,下面则是Python环境进行交互式操作,输入jupyter notebook
,之后在网页上打开
首先是导入各种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
然后加载分析结果,包括聚类信息和特征信息
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))
如果要矫正批次效应,可以通过根据feature的heatmap,去掉和batch相关的feature来实现
我们可以展示SCALE对原始数据纠正后的值(imputed data), 该结果也能提高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))
聚类特异性peak
通过SCALE.py -h
可以输出SCALE的所有可用参数
-d/--dataset
: 单个文件矩阵应该指定文件路径,10X输出的多个文件则是文件目录-k
: 设定输出结果的聚类数-o
: 输出文件路径--pretrain
: 读取之前训练的模型--lr
: 修改起始学习速率, 默认是0.002,和模型训练有关--batch_size
: 批处理大小, 默认就行,不需要修改(和批次效应处理无关)-g GPU
: 选择GPU设备数目,非GPU服务器用不到--seed
: 初始随机数种子,通常在遇到nan缺失时考虑修改-encode_dim
, -decode_dim
: 编码器和解码器的维度,通常也不需要修改-latent
隐藏层维度--low
, --high
: 过滤低质量的peak, 即出现比例高于或者低于某个阈值的peak,默认是0.01和0.9。作者推荐保留1万-3万的peak用于SCALE分析。如果数据质量很高,且peak数不多于10万,那么可以不过滤。--min_peaks
: 过滤低质量细胞,如果该细胞的peak低于阈值log_transform
: log2(x+1)
的变换--max_iter
: 最大迭代数,默认是30000, 可以观察损失收敛的情况来修改,也就是训练模型这一步输出的信息-weight_decay
: 没有说明--impute
: 保存推断数据,默认开启--binary
: 推荐加上该参数,减少imputed data占用空间--no_tsne
: 不需要保存t-SNE结果--reference
: 参考细胞类型-t
: 如果输出矩阵是列为peak,行为barcode,用该参数进行转置对于使用者而言,我们一般只用修改-k
更改最后的聚类数,--low
, --high
, ---min_peaks
来对原始数据进行过滤,以及加上--binary
节约空间。
假如在训练模型阶段,发现输出信息为loss=nan recon_loss=nan kl_loss=nan
,十有八九最终会报错退出, 可以如下的参数调整
--seed
-x 0.04
或 -x 0.06
--lr 0.0002
ChIP-seq基本分析流程适用于ATAC-seq的前期分析。