单细胞分析的 Python 包 Scanpy(图文详解)
文章目录
- 一、安装
- 二、使用
- 1、准备工作
- 2、预处理
- 过滤低质量细胞样本
- 3、检测特异性基因
- 4、主成分分析(Principal component analysis)
- 5、领域图,聚类图(Neighborhood graph)
- 6、检索标记基因
- 7、保存数据
- 8、番外
一、安装
如果没有conda 基础,参考: Conda 安装使用图文详解(2021版)
pip install scanpy
conda install -y -c conda-forge leidenalg
二、使用
1、准备工作
# 载入包
import numpy as np
import pandas as pd
import scanpy as sc# 设置
sc.settings.verbosity = 3 # 设置日志等级: errors (0), warnings (1), info (2), hints (3)
sc.logging.print_header()
sc.settings.set_figure_params(dpi=80, facecolor='white')# 用于存储分析结果文件的路径
results_file = 'write/pbmc3k.h5ad'# 载入文件
adata = sc.read_10x_mtx('./filtered_gene_bc_matrices/hg19/', # mtx 文件目录var_names='gene_symbols', # 使用 gene_symbols 作为变量名cache=True) # 写入缓存,可以更快的读取文件
2、预处理
显示在所有细胞中在每个单细胞中产生最高计数分数的基因
sc.pl.highest_expr_genes(adata, n_top=20, )
过滤低质量细胞样本
过滤在少于三个细胞中表达,或一个细胞中表达少于200个基因的细胞样本
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
过滤包含线粒体基因和表达基因过多的细胞
线粒体基因的转录本比单个转录物分子大,并且不太可能通过细胞膜逃逸。因此,检测出高比例的线粒体基因,表明细胞质量差(Islam et al. 2014; Ilicic et al. 2016)。
表达基因过多可能是由于一个油滴包裹多个细胞,从而检测出比正常检测要多的基因数,因此要过滤这些细胞。
adata.var['mt'] = adata.var_names.str.startswith('MT-') # 将线粒体基因标记为 mt
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)
sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],jitter=0.4, multi_panel=True)
生成的三张小提琴图代表:表达基因的数量,每个细胞包含的表达量,线粒体基因表达量的百分比。
过滤
sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt')
sc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts')
过滤线粒体基因表达过多或总数过多的细胞,也就是红框标识的样本。
# 获取线粒体基因占比在 5% 以下的细胞样本
adata = adata[adata.obs.pct_counts_mt < 5, :]
# 获取表达基因数在 2500 以下的细胞样本
adata = adata[adata.obs.n_genes_by_counts < 2500, :]
3、检测特异性基因
归一化
# 归一化,使得不同细胞样本间可比
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
存储数据
将 AnnData 对象的 .raw 属性设置为归一化和对数化的原始基因表达,以便以后用于基因表达的差异测试和可视化。这只是冻结了 AnnData 对象的状态。
adata.raw = adata
识别特异性基因
# 计算
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
# 绘制特异性基因散点图
sc.pl.highly_variable_genes(adata)
获取只有特异性基因的数据集
# 获取只有特异性基因的数据集
adata = adata[:, adata.var.highly_variable]
# 回归每个细胞的总计数和表达的线粒体基因的百分比的影响。
sc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt'])
# 将每个基因缩放到单位方差。阈值超过标准偏差 10。
sc.pp.scale(adata, max_value=10)
4、主成分分析(Principal component analysis)
通过运行主成分分析 (PCA) 来降低数据的维数,可以对数据进行去噪并揭示不同分群的主因素。
# 绘制 PCA 图
sc.tl.pca(adata, svd_solver='arpack')
sc.pl.pca(adata, color='CST3')
检查单个 PC 对数据总方差的贡献,这可以提供给我们应该考虑多少个 PC 以计算细胞的邻域关系的信息,例如用于后续的聚类函数 sc.tl.louvain() 或 tSNE 聚类 sc.tl.tsne()。
sc.pl.pca_variance_ratio(adata, log=True)
5、领域图,聚类图(Neighborhood graph)
使用数据矩阵的 PCA 表示来计算细胞的邻域图。为了重现 Seurat 的结果,我们采用以下值。
建议使用 UMAP ,它可能比 tSNE 更忠实于流形的全局连通性,因此能更好地保留轨迹。
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)
sc.tl.umap(adata)
# 如果设置了 adata 的 .raw 属性时,下图显示了“raw”(标准化、对数化但未校正)基因表达矩阵。
sc.pl.umap(adata, color=['CST3', 'NKG7', 'PPBP'])
为了绘制缩放矫正的基因表达聚类图,需要使用 use_raw=False 参数。
sc.pl.umap(adata, color=['CST3', 'NKG7', 'PPBP'], use_raw=False)
目前还没有计算出各个细胞类群,下面进行聚类
Leiden 图聚类
# 计算
sc.tl.leiden(adata)
# 绘制
sc.pl.umap(adata, color=['leiden'])
6、检索标记基因
先计算每个 leiden 分群中高度差异基因的排名,取排名前 25 的基因。
默认情况下,使用 AnnData 的 .raw 属性。
T-test
最简单和最快的方法是 t 检验。
sc.tl.rank_genes_groups(adata, 'leiden', method='t-test')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)
Wilcoxon rank-sum
Wilcoxon rank-sum (Mann-Whitney-U) 检验的结果非常相似,还可以使用其他的差异分析包,如 MAST、limma、DESeq2 和 diffxpy。
sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)
# 保存这次的数据结果
adata.write(results_file)
逻辑回归
sc.tl.rank_genes_groups(adata, 'leiden', method='logreg')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)
使用逻辑回归对基因进行排名 Natranos et al. (2018),这里使用多变量方法,而传统的差异测试是单变量 Clark et al. (2014)
除了仅由 t 检验发现的 IL7R 和由其他两种方法发现的 FCER1A 之外,所有标记基因都在所有方法中都得到了重现。
Louvain Group | Markers | Cell Type |
---|---|---|
0 | IL7R | CD4 T cells |
1 | CD14, LYZ | CD14+ Monocytes |
2 | MS4A1 | B cells |
3 | CD8A | CD8 T cells |
4 | GNLY, NKG7 | NK cells |
5 | FCGR3A, MS4A7 | FCGR3A+ Monocytes |
6 | FCER1A, CST3 | Dendritic Cells |
7 | PPBP | Megakaryocytes |
根据已知的标记基因,定义一个标记基因列表供以后参考:
marker_genes = ['IL7R', 'CD79A', 'MS4A1', 'CD8A', 'CD8B', 'LYZ', 'CD14','LGALS3', 'S100A8', 'GNLY', 'NKG7', 'KLRB1','FCGR3A', 'MS4A7', 'FCER1A', 'CST3', 'PPBP']
载入数据
# 使用 Wilcoxon Rank-Sum 测试结果重新加载已保存的对象
adata = sc.read(results_file)
获取聚类分组和分数
result = adata.uns['rank_genes_groups']
groups = result['names'].dtype.names
pd.DataFrame({group + '_' + key[:1]: result[key][group]for group in groups for key in ['names', 'pvals']}).head(5)
Group 1 与 Group 2 比较,进行差异分析
sc.tl.rank_genes_groups(adata, 'leiden', groups=['0'], reference='1', method='wilcoxon')
sc.pl.rank_genes_groups(adata, groups=['0'], n_genes=20)
sc.pl.rank_genes_groups_violin(adata, groups='0', n_genes=8)
Group 0 与其余组的比较进行差异分析
adata = sc.read(results_file)
sc.pl.rank_genes_groups_violin(adata, groups='0', n_genes=8)
跨类群比较基因
sc.pl.violin(adata, ['CST3', 'NKG7', 'PPBP'], groupby='leiden')
根据已知的细胞标记,注释细胞类型
new_cluster_names = ['CD4 T', 'CD14 Monocytes','B', 'CD8 T','NK', 'FCGR3A Monocytes','Dendritic', 'Megakaryocytes']
adata.rename_categories('leiden', new_cluster_names)
sc.pl.umap(adata, color='leiden', legend_loc='on data', title='', frameon=False, save='.pdf')
可视化每个类群的标记基因
气泡图显示:
sc.pl.dotplot(adata, marker_genes, groupby='leiden');
小提琴图显示
sc.pl.stacked_violin(adata, marker_genes, groupby='leiden', rotation=90);
7、保存数据
保存压缩文件
如果只想将其用于可视化的人共享此文件,减少文件大小的一种简单方法是删除缩放和校正的数据矩阵。
adata.write(results_file, compression='gzip')
保存为 h5ad 数据
adata.raw.to_adata().write('./write/pbmc3k_withoutX.h5ad')
读取使用 adata = sc.read_h5ad(’./write/pbmc3k_withoutX.h5ad’)
导出数据子集
# 导出聚类数据
adata.obs[['n_counts', 'louvain_groups']].to_csv('./write/pbmc3k_corrected_louvain_groups.csv')
# 导出PCA数据
adata.obsm.to_df()[['X_pca1', 'X_pca2']].to_csv('./write/pbmc3k_corrected_X_pca.csv')
8、番外
我之前在处理较多数据量的时候,会有些地方不一样,具体每个数据集的处理也会有比较大的自由度,比如:
在检测线粒体基因时,这里在质控时,已经把线粒体基因直接剔除。
在做 UMAP 时,可以看到一些类群间的联系和轨迹。
做 TSNE时,可以看到类群间比较干净利索,整体比较“饱满”。
其他 Scanpy 的使用教程:
scanpy 单细胞分析包图文详解 01 | 深入理解 AnnData 数据结构
单细胞分析的 Python 包 Scanpy(图文详解)相关推荐
- java python c++比喻图_Java/Python/PHP/C++图文详解它们之间的尿性
PHP:没有优点 Java:就是库多 Python:语法清晰 语法清晰 C:能操纵底层,最细粒度优化没有之一 C++:啥都有,啥都有,啥都有 等等等之类不在例举.直接上图吧. SQL: ps:千万别惹 ...
- php尿,Java/Python/PHP/C 图文详解它们之间的尿性
PHP:没有优点 Java:就是库多 Python:语法清晰 语法清晰 C:能操纵底层,最细粒度优化没有之一 C++:啥都有,啥都有,啥都有 等等等之类不在例举.直接上图吧. SQL: ps:千万别惹 ...
- python爬虫图片实例-【图文详解】python爬虫实战——5分钟做个图片自动下载器...
我想要(下)的,我现在就要 python爬虫实战--图片自动下载器 之前介绍了那么多基本知识[Python爬虫]入门知识(没看的赶紧去看)大家也估计手痒了.想要实际做个小东西来看看,毕竟: talk ...
- python详细安装教程环境配置-python3.6环境安装+pip环境配置教程图文详解
1.python安装可以跨平台 2.有两个版本2.7和3.6,第三方库适用2.7版,两个版本不兼容 windows安装: 第一种方法官网安装: 在官网下载安装包如图: 图下点击是默认下载32位所以我们 ...
- 全网最全的Windows下Anaconda2 / Anaconda3里Python语言实现定时发送微信消息给好友或群里(图文详解)...
不多说,直接上干货! 缘由: (1)最近看到情侣零点送祝福,感觉还是很浪漫的事情,相信有很多人熬夜为了给爱的人送上零点祝福,但是有时等着等着就睡着了或者时间并不是卡的那么准就有点强迫症了,这是也许程序 ...
- spark最新源码下载并导入到开发环境下助推高质量代码(Scala IDEA for Eclipse和IntelliJ IDEA皆适用)(以spark2.2.0源码包为例)(图文详解)...
不多说,直接上干货! 前言 其实啊,无论你是初学者还是具备了有一定spark编程经验,都需要对spark源码足够重视起来. 本人,肺腑之己见,想要成为大数据的大牛和顶尖专家,多结合源码和操练编程. ...
- python利器的使用-图文详解python开发利器之ulipad的使用实践
Ulipad是一个国人limodou编写的专业Python编辑器,它基于wxpython开发的GUI(图形化界面).下面这篇文章主要介绍了python开发利器之ulipad的使用实践,文中介绍的非常详 ...
- ftp linux包,图文详解Ubuntu搭建Ftp服务器的方法(包成功)
一.今天下午由于课程的要求不得已做了Ubuntu搭建Ftp服务器的实验,但是实验指导书还是N年前的技术,网上搜了一大把,都是模模糊糊的! 在百般困难中终于试验成功,特把经验分给大家 希望大家少走弯路! ...
- python path包的使用详解
python path包的使用详解 files函数的使用 files函数的使用 import argparse from path import Pathparser = argparse.Argum ...
最新文章
- 超图桌面版制作一幅简单专题图示例
- 成功解决ImportError: cannot import name 'pywrap_tensorflow'
- 各类神经网络知识收集
- 使用 Proto 构建了一个简单但功能强大的 lambda 库的测试程序
- 安装JDK-- Java基础
- 阿里云PolarDB发布重大更新 支持Oracle等数据库一键迁移上云
- Ken Block 漂移大叔,程序实现精准漂移算法。
- iPhone 12蓝色版疑似翻车:眼前的蓝不是蓝......
- Ubuntu 11.10 系统启动默认进入终端
- jsonView插件的安装方法
- Windows 10开启Teredo隧道连接IPV6
- Sklearn官方文档中文整理4——随机梯度下降和最近邻篇
- 微信小程序 录音之获取、保存、读取
- Java Web图书管理系统(MVC框架)-包含源码
- 无分类编址(超网)中的网络前缀
- python制作小提琴图
- 中文***测试专用Linux系统—MagicBox(魔方系统)
- matlab中离散信号模型
- 从翻译到产品本地化 Airbnb语言经理的出海建议
- 【论文阅读】结合空洞卷积的 FuseNet变体网络高分辨率遥感影像语义分割
热门文章
- java计算机毕业设计景区在线购票系统MyBatis+系统+LW文档+源码+调试部署
- Hive结合Apache Ranger进行数据脱敏
- 数据库.net3.5安装失败错误0X80070422
- Pact broker应用-overview
- iOS-Provisioning Profile#40;Certificate#41;与Code Signing详解
- 表单提交方式get和post的区别
- 计算机网络之XSS攻击
- 豆瓣电影最新API接口(亲测可用)
- plc和pc串口通讯接线_三菱FX系列PLC与电脑之间串口RS232通讯协议简易解析
- echarts gallery最新地址