单细胞分析的 Python 包 Scanpy(图文详解)

一、安装

Conda 安装使用图文详解(2021版)

scanpy 单细胞分析包图文详解 01 | 深入理解 AnnData 数据结构

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, )
image.png

过滤低质量细胞样本

过滤在少于三个细胞中表达,或一个细胞中表达少于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)
image.png

生成的三张小提琴图代表:表达基因的数量,每个细胞包含的表达量,线粒体基因表达量的百分比。

过滤

sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt')
sc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts')
image.png

过滤线粒体基因表达过多或总数过多的细胞,也就是红框标识的样本。

# 获取线粒体基因占比在 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)
image.png

获取只有特异性基因的数据集

# 获取只有特异性基因的数据集
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.pl.pca(adata, color='CST3')
image.png

检查单个 PC 对数据总方差的贡献,这可以提供给我们应该考虑多少个 PC 以计算细胞的邻域关系的信息,例如用于后续的聚类函数 sc.tl.louvain() 或 tSNE 聚类 sc.tl.tsne()。

sc.pl.pca_variance_ratio(adata, log=True)
image.png

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'])
image.png

为了绘制缩放矫正的基因表达聚类图,需要使用 use_raw=False 参数。

sc.pl.umap(adata, color=['CST3', 'NKG7', 'PPBP'], use_raw=False)
image.png

目前还没有计算出各个细胞类群,下面进行聚类

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)
image.png

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)

image.png

除了仅由 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)
image.png

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)
image.png
sc.pl.rank_genes_groups_violin(adata, groups='0', n_genes=8)
image.png

Group 0 与其余组的比较进行差异分析

adata = sc.read(results_file)
sc.pl.rank_genes_groups_violin(adata, groups='0', n_genes=8)
image.png

跨类群比较基因

sc.pl.violin(adata, ['CST3', 'NKG7', 'PPBP'], groupby='leiden')
image.png

根据已知的细胞标记,注释细胞类型

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')
image.png

可视化每个类群的标记基因

气泡图显示:

sc.pl.dotplot(adata, marker_genes, groupby='leiden');
image.png

小提琴图显示

sc.pl.stacked_violin(adata, marker_genes, groupby='leiden', rotation=90);
image.png

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、番外

我之前在处理较多数据量的时候,会有些地方不一样,具体每个数据集的处理也会有比较大的自由度,比如:

在检测线粒体基因时,这里在质控时,已经把线粒体基因直接剔除。

image.png

在做 UMAP 时,可以看到一些类群间的联系和轨迹。

image.png

做 TSNE时,可以看到类群间比较干净利索,整体比较“饱满”。

image.png
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 144,247评论 1 305
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 61,830评论 1 258
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 95,531评论 0 214
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 41,345评论 0 183
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 49,160评论 1 260
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 38,936评论 1 178
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 30,538评论 2 275
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 29,291评论 0 168
  • 想象着我的养父在大火中拼命挣扎,窒息,最后皮肤化为焦炭。我心中就已经是抑制不住地欢快,这就叫做以其人之道,还治其人...
    爱写小说的胖达阅读 29,162评论 6 237
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 32,654评论 0 214
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 29,401评论 2 217
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 30,747评论 1 232
  • 白月光回国,霸总把我这个替身辞退。还一脸阴沉的警告我。[不要出现在思思面前, 不然我有一百种方法让你生不如死。]我...
    爱写小说的胖达阅读 24,297评论 1 33
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 27,206评论 2 213
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 31,670评论 3 213
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 25,661评论 0 9
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,089评论 0 169
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 33,677评论 2 233
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 33,819评论 2 237

推荐阅读更多精彩内容