使用SCALE分析单细胞ATAC-seq数据

2024-06-23 20:18

本文主要是介绍使用SCALE分析单细胞ATAC-seq数据,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

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 = muz sigmaZ times epsilon$ 计算而得,公式里面的 $muz$ $sigmaz$ 是编码器网络从x中学习而得,$epsilon$ 则是从 $mathbb{N}(0,mathbf{I})$ 抽样而成。

从公式中我们还可以发现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能做以下的分析

  • SCALE可以对隐特征聚类识别细胞类群
  • SCALE可以降噪,恢复缺失的peak
  • 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支持两类输入文件:

  • count矩阵,行为peak,列为barcode
  • 10X输出文件: count.mtx.gz, peak.tsv, barcode.tsv

我们以官方提供的Forebrain数据集为例进行介绍,因为这个数据相对于另外一个数据集Mouse Atlas小多了。

我们在服务器上新建一个文件夹,用于存放从 下载的数据

mkdir Forebrain

保证Forebrain有下载好的数据

$ ls Forebrain 
data.txt

之后运行程序

SCALE.py -d Forebrain/data.txt -k 8 --impute

软件运行步骤为:

  • 加载数据: Loading data
  • 模型训练: Training Model
  • 输出结果: Saving imputed data

其中模型训练这一步时间比较久,可以尝试用GPU加速(我是普通CPU服务器没有办法)。最终会在当前文件夹看到一个output文件夹,里面有如下内容:

  • imputed_data.txt: 每个细胞在每个特征的推断值,建议用`--binary`保存二进制格式
  • model.pt: 用于重复结果的模型文件,--pretrain参数能够读取该模型
  • feature.txt: 每个细胞的隐特征,用于聚类和可视化
  • cluster_assignments.txt: 两列,barcode和所属类群
  • tsne.txt, tsne.pdf: tSNE的坐标和PDF文件,坐标文件可以导入到R语言进行可视化

上面是命令行部分,下面则是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))

heatmap

如果要矫正批次效应,可以通过根据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_scorecluster_specific完成

from scale.specifity import cluster_specific, mat_specificity_scorescore_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节约空间。

版权声明:本博客所有文章除特别声明外,均采用 知识共享署名-非商业性使用-禁止演绎 4.0 国际许可协议 (CC BY-NC-ND 4.0) 进行许可。

扫码即刻交流

这篇关于使用SCALE分析单细胞ATAC-seq数据的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/1088169

相关文章

使用Python从PPT文档中提取图片和图片信息(如坐标、宽度和高度等)

《使用Python从PPT文档中提取图片和图片信息(如坐标、宽度和高度等)》PPT是一种高效的信息展示工具,广泛应用于教育、商务和设计等多个领域,PPT文档中常常包含丰富的图片内容,这些图片不仅提升了... 目录一、引言二、环境与工具三、python 提取PPT背景图片3.1 提取幻灯片背景图片3.2 提取

C++如何通过Qt反射机制实现数据类序列化

《C++如何通过Qt反射机制实现数据类序列化》在C++工程中经常需要使用数据类,并对数据类进行存储、打印、调试等操作,所以本文就来聊聊C++如何通过Qt反射机制实现数据类序列化吧... 目录设计预期设计思路代码实现使用方法在 C++ 工程中经常需要使用数据类,并对数据类进行存储、打印、调试等操作。由于数据类

使用Python实现图像LBP特征提取的操作方法

《使用Python实现图像LBP特征提取的操作方法》LBP特征叫做局部二值模式,常用于纹理特征提取,并在纹理分类中具有较强的区分能力,本文给大家介绍了如何使用Python实现图像LBP特征提取的操作方... 目录一、LBP特征介绍二、LBP特征描述三、一些改进版本的LBP1.圆形LBP算子2.旋转不变的LB

Maven的使用和配置国内源的保姆级教程

《Maven的使用和配置国内源的保姆级教程》Maven是⼀个项目管理工具,基于POM(ProjectObjectModel,项目对象模型)的概念,Maven可以通过一小段描述信息来管理项目的构建,报告... 目录1. 什么是Maven?2.创建⼀个Maven项目3.Maven 核心功能4.使用Maven H

Python中__init__方法使用的深度解析

《Python中__init__方法使用的深度解析》在Python的面向对象编程(OOP)体系中,__init__方法如同建造房屋时的奠基仪式——它定义了对象诞生时的初始状态,下面我们就来深入了解下_... 目录一、__init__的基因图谱二、初始化过程的魔法时刻继承链中的初始化顺序self参数的奥秘默认

SpringBoot使用GZIP压缩反回数据问题

《SpringBoot使用GZIP压缩反回数据问题》:本文主要介绍SpringBoot使用GZIP压缩反回数据问题,具有很好的参考价值,希望对大家有所帮助,如有错误或未考虑完全的地方,望不吝赐教... 目录SpringBoot使用GZIP压缩反回数据1、初识gzip2、gzip是什么,可以干什么?3、Spr

Java程序进程起来了但是不打印日志的原因分析

《Java程序进程起来了但是不打印日志的原因分析》:本文主要介绍Java程序进程起来了但是不打印日志的原因分析,具有很好的参考价值,希望对大家有所帮助,如有错误或未考虑完全的地方,望不吝赐教... 目录Java程序进程起来了但是不打印日志的原因1、日志配置问题2、日志文件权限问题3、日志文件路径问题4、程序

Spring Boot 集成 Quartz并使用Cron 表达式实现定时任务

《SpringBoot集成Quartz并使用Cron表达式实现定时任务》本篇文章介绍了如何在SpringBoot中集成Quartz进行定时任务调度,并通过Cron表达式控制任务... 目录前言1. 添加 Quartz 依赖2. 创建 Quartz 任务3. 配置 Quartz 任务调度4. 启动 Sprin

Linux下如何使用C++获取硬件信息

《Linux下如何使用C++获取硬件信息》这篇文章主要为大家详细介绍了如何使用C++实现获取CPU,主板,磁盘,BIOS信息等硬件信息,文中的示例代码讲解详细,感兴趣的小伙伴可以了解下... 目录方法获取CPU信息:读取"/proc/cpuinfo"文件获取磁盘信息:读取"/proc/diskstats"文

Java使用SLF4J记录不同级别日志的示例详解

《Java使用SLF4J记录不同级别日志的示例详解》SLF4J是一个简单的日志门面,它允许在运行时选择不同的日志实现,这篇文章主要为大家详细介绍了如何使用SLF4J记录不同级别日志,感兴趣的可以了解下... 目录一、SLF4J简介二、添加依赖三、配置Logback四、记录不同级别的日志五、总结一、SLF4J