引言

本系列讲解 使用Scanpy分析单细胞(scRNA-seq)数据教程[1],持续更新,欢迎关注,转发!
数据集
本次使用的数据集包含一位健康供体的3k PBMCs,这些数据可以从10x Genomics的官方网站免费获取。如果你使用的是unix系统,只需取消以下代码的注释并运行,就可以完成数据的下载和解压操作。另外,代码的最后一行会自动创建一个目录,用于存放处理后的数据。
# !mkdir -p data # !curl https://cf.10xgenomics.com/samples/cell-exp/1.1.0/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz -o data/pbmc3k_filtered_gene_bc_matrices.tar.gz # !cd data; tar -xzf pbmc3k_filtered_gene_bc_matrices.tar.gz # !mkdir -p write数据导入
import pandas as pd import scanpy as sc sc.settings.verbosity = 3 # verbosity: 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" # the file that will store the analysis results把计数矩阵导入到一个 AnnData 对象里,这个对象能存储很多注释信息以及数据的各种不同表现形式。它还自带一种基于 HDF5 的文件格式,即 .h5ad。
adata = sc.read_10x_mtx( "data/filtered_gene_bc_matrices/hg19/", # the directory with the `.mtx` file var_names="gene_symbols", # use gene symbols for the variable names (variables-axis index) cache=True, # write a cache file for faster subsequent reading ) adata.var_names_make_unique() # this is unnecessary if using `var_names=gene_ids` in `sc.read_10x_mtx`预处理
找出在每个单细胞里计数占比最高的那些基因,范围覆盖所有细胞。
sc.pl.highest_expr_genes(adata, n_top=20)我们来收集一些有关线粒体基因的信息,这些基因对于质量控制来说很重要。
如果细胞里线粒体基因的比例很高,那就意味着这个细胞的质量不太好。出现这种情况可能是因为细胞膜有破损,细胞质 RNA 从破洞里流失了。因为线粒体比单个转录本分子大,所以相对来说更难从细胞膜的裂缝里跑出去。
通过使用 pp.calculate_qc_metrics 这个工具,我们可以非常高效地计算出很多质量控制指标。
# annotate the group of mitochondrial genes as "mt" adata.var["mt"] = adata.var_names.str.startswith("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")具体操作是通过切片 AnnData 对象来进行过滤。
adata = adata[adata.obs.n_genes_by_counts < 2500, :] adata = adata[adata.obs.pct_counts_mt < 5, :].copy()对数据矩阵进行总计数归一化(也就是库大小校正) 每个细胞的读段数调整到10,000个,这样就可以让不同细胞之间的计数具有可比性了。
sc.pp.normalize_total(adata, target_sum=1e4)对数转换
sc.pp.log1p(adata)识别highly-variable基因
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5) sc.pl.highly_variable_genes(adata)把 AnnData 对象的 .raw 属性设置为经过归一化和对数化处理的原始基因表达数据,以便后续用于差异测试和基因表达的可视化。这一步相当于把 AnnData 对象当前的状态给固定下来。
adata.raw = adata.copy()进行过滤
adata = adata[:, adata.var.highly_variable]对每个细胞的总计数和线粒体基因表达百分比的影响进行回归分析。然后将数据缩放到单位方差。
sc.pp.regress_out(adata, ["total_counts", "pct_counts_mt"])对每个基因也进行缩放,使其方差为单位方差。并且,将超过标准差10的值进行裁剪处理。
sc.pp.scale(adata, max_value=10)[1]Source: https://scanpy.readthedocs.io/en/stable/tutorials/basics/clustering-2017.html