Scanpy 分析 3k PBMCs：数据预处理
Scanpy 分析 3k PBMCs：数据预处理
引言
本系列讲解 使用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 = # verbosity: errors (0), warnings (1), info (2), hints (3)
sc.logging.print_header()
sc.settings.set_figure_params(dpi=, 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=)
sc.pp.filter_cells(adata, min_genes=)
sc.pp.filter_genes(adata, min_cells=)
我们来收集一些有关线粒体基因的信息，这些基因对于质量控制来说很重要。
如果细胞里线粒体基因的比例很高，那就意味着这个细胞的质量不太好。出现这种情况可能是因为细胞膜有破损，细胞质 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 < , :]
adata = adata[adata.obs.pct_counts_mt < , :].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=, 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=)
Reference
[1]
Source: https://scanpy.readthedocs.io/en/stable/tutorials/basics/clustering-2017.html
腾讯云开发者
