专栏首页生信菜鸟团Python 做 Nature 级的单细胞分析(图文详解)

Python 做 Nature 级的单细胞分析(图文详解)

安装

如果 conda 不熟悉的小伙伴,可以参考:https://blog.csdn.net/u011262253/article/details/88828229

数据下载:http://cf.10xgenomics.com/samples/cell-exp/1.1.0/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz

 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

主成分分析

通过运行主成分分析 (PCA) 来降低数据的维数,可以对数据进行去噪并揭示不同分群的主因素。

 # 绘制 PCA 图
 sc.pl.pca(adata, color='CST3')

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

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

5

细胞分群

使用数据矩阵的 PCA 表示来计算细胞的邻域图。

建议使用 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时,可以看到类群间比较干净,如果查看类群区别可以用该方法。需要注意与 PCA 不同的是,TSNE 的每个类群间的距离并不能代表其实际相关性。

如果想了解更多关于 AnnoData 的知识,可以参考:

用 Python 做单细胞分析 01 | 详解 AnnData 数据结构


本文分享自微信公众号 - 生信菜鸟团(bio_123456789),作者:白墨

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2021-07-21

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

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

    线粒体基因的转录本比单个转录物分子大,并且不太可能通过细胞膜逃逸。因此,检测出高比例的线粒体基因,表明细胞质量差(Islam et al. 2014; Ilic...

    白墨石
  • 细胞通讯分析结果的解读

    不过,虽然细胞通讯分析越来越普通,但它的难度并不会降低,在试图学习这个分析方法之前,大家需要自己提前了解一下:细胞通讯分析的背景知识,而且呢,还得看看细胞通讯分...

    生信技能树
  • Nature | ​心脏发生的单细胞分析揭示了器官水平发育缺陷的基础

    对单细胞技术感兴趣?点击浅蓝色字 — 中科院的算法开发博士带你真正玩转这项平均每个月都有多篇高IF文章的技术

    生信宝典
  • 单细胞转录组高级分析三:细胞通讯分析

    细胞通讯研究领域涵盖的内容很广,如上图所示包括通讯方式、功能、信号分子以及各种途径的机制。细胞之间通讯的介质有很多,例如钙离子、脂质、多肽、蛋白、外泌体以及电信...

    生信技能树jimmy
  • NBT|45种单细胞轨迹推断方法比较,110个实际数据集和229个合成数据集

    轨迹推断(Trajectory Inference,TI),是分析从千上万单细胞的组学数据中推断细胞发育轨迹的重要方法,也被称为伪时序分析 (pseudotim...

    生信宝典
  • 想分析单细胞RNA的动态变化?

    当你的才华还撑不起你的野心时,请潜下心来,脚踏实地,跟着我们慢慢进步。不知不觉在单细胞转录组领域做知识分析也快两年了,通过文献速递这个栏目很幸运聚集了一些小伙伴...

    生信技能树jimmy
  • 振聋发聩 | 新科诺奖得主William Kaelin发文批评当下论文数据“华而不实”

    William G. Kaelin,哈佛大学医学院教授、HHMI研究员,与牛津大学的Peter J.Ratcliffe、约翰霍普金斯大学的Gregg L.Sem...

    生信宝典
  • 单细胞RNA数据标准化与聚类分析

    单细胞转录组测序产生的数据是成百上千个基因在上万个细胞中的表达情况,属于高维数据,我们需要对数据进行严格的质控与过滤,将合格的数据降维到低维子空间,使数据可视化...

    生信交流平台
  • Cell 单细胞文章 | 白凡课题组与合作团队揭示儿童结肠炎及炎症性肠病的致病机制及治疗方法

    2019年11月14日,北京大学白凡研究员团队与广州医科大学附属广州市妇女儿童医疗中心儿科研究所张玉霞研究员,国家临床重点专科儿科消化团队(杨敏、耿岚岚及龚四堂...

    生信宝典
  • 单细胞高分文献精读系列(1)--近两年文献纵横

    回顾近两年的生物技术热点,单细胞分析技术无疑是各个研究领域大佬们力捧的宠儿。从去年到现在,影响因子20分以上、以单细胞分析为主体的文章已发表近100篇(文末识别...

    用户6317549
  • 僵小鱼的故事

    许多年之后,面对同一个作图需求,僵小鱼将会回想起,在微信群里提出相同问题的那个遥远的上午

    生信技能树jimmy
  • 翻车实录之Nature Medicine新冠单细胞文献|附全代码

    NGS系列文章包括NGS基础、转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这)、ChIP-seq分析 (ChIP-seq基本分析流程)、...

    生信宝典
  • 两个小鼠器官单细胞图谱

    第一个是2018年2月 发表在Cell(《细胞》)杂志的,浙江大学郭国骥团队的研究成果。其课题组自主开发了一套完全国产化的Microwell-seq高通量单细胞...

    生信技能树jimmy
  • 菜鸟团一周文献推荐(No.56)

    题目:Animal domestication in the era of ancient genomics

    生信菜鸟团
  • velocyto||sc-RNA速率:一种细胞轨迹推断方法

    18年nature发了一篇单细胞方法类文章,讲得就是如何利用RNA velocity来做细胞发育路径的推断。

    生信编程日常
  • 单细胞转录因子分析之SCENIC流程

    SCENIC (Single-Cell rEgulatory Network Inference and Clustering) is a computatio...

    生信技能树
  • 【重磅】灵长类动物脸部识别算法被破译,大脑黑箱或根本不存在

    【新智元导读】发表在 Cell 的一项研究揭示了人脸识别的具体神经元活动过程。对猕猴的实验表明,对脸部的识别是由大脑中 200 多个不同神经元共同编码完成的,每...

    新智元
  • 单细胞转录组高级分析二:转录调控网络分析

    组织内细胞异质性的基础是细胞转录状态的差异,转录状态的特异性又是由转录因子主导的基因调控网络(GRNs)决定并维持稳定的。因此分析单细胞的GRNs有助于深入挖掘...

    生信技能树jimmy
  • 单细胞转录组探索小鼠肝脏发育

    肝脏是一种多倍体器官,由具有一个或两个细胞核的肝细胞组成,每个细胞核含有2,4,8或更多单倍体染色体组。

    生信技能树jimmy

扫码关注云+社区

领取腾讯云代金券