Scanpy（3）单细胞数据分析常规流程

单细胞数据分析常规流程

面对高效快速的要求上，使用R分析数据越来越困难，转战Python分析，我们通过scanpy官网去学习如何分析单细胞下游常规分析。

数据3k PBMC来自健康的志愿者，可从10x Genomics免费获得。在linux系统上，可以取消注释并运行以下操作来下载和解压缩数据。最后一行创建一个用于保存已处理数据的目录write，后面直接使用保存的数据，能快速加载数据。

下载数据：

$mkdir data
$cd data
$wget http://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
$tar -xzf pbmc3k_filtered_gene_bc_matrices.tar.gz
# 获得数据

1. 数据加载

import numpy as np
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')

# 声明h5ad用于存储分析结果
results_file = 'data/write/pbmc3k.h5ad'

adata = sc.read_10x_mtx(
    'data/filtered_gene_bc_matrices/hg19/',  # `.mtx`文件所在的目录
    var_names='gene_symbols',                # 用 gene 作为var
    cache=True) # 开启缓存读写

"""
注意cache=Trure
... writing an h5ad cache file to speedup reading next time
下次读取就不会从count matrix读, 会直接从cache目录下的h5ad文件读(更快)
"""

在函数 sc.read_10x_mtx 中，参数 var_names 用于指定在加载数据时使用哪个变量来作为基因的名称。在这里，如果你将 var_names='gene_ids'，它将使用基因的唯一标识符作为变量名，而如果你将 var_names='gene_symbols'，它将使用基因的符号名称作为变量名。

这两者之间的区别在于：

1. gene_ids：使用基因的唯一标识符作为变量名。这通常是一种更确切和唯一的标识，不同基因之间不存在重复。使用基因的唯一标识符作为变量名可以确保在分析中每个基因都有唯一的标识符，并且不会出现混淆或重复。

2. gene_symbols：使用基因的符号名称作为变量名。基因的符号名称通常更容易理解和记忆，因为它们通常是基于基因的功能或特征而命名的。然而，基因的符号名称不一定是唯一的，可能存在多个基因具有相同的符号名称，这可能会导致一些混淆或不一致。

因此，你可以根据具体的需求和分析的目的来选择使用哪种类型的变量名。如果需要确保每个基因都具有唯一的标识符，并且不会出现混淆或重复，那么可以使用 gene_ids。如果更关注基因的功能或特征，并且不太担心可能存在的重复符号名称，那么可以使用 gene_symbols。

---

注意，如果在函数sc.read_10x_mtx中指定参数var_names='gene_ids'时，下一个操作将是不必要的：

# 消除重复的列
adata.var_names_make_unique()

print(adata)

AnnData object with n_obs × n_vars = 2700 × 32738
    var: 'gene_ids'

adata包含2700个细胞、32738个基因的对象

2. top基因箱型图

下图计算每一个基因在所有细胞中的平均表达量，并绘制了平均表达量前30的基因箱型图。

sc.pl.highest_expr_genes(adata, n_top=30)

计算每一个基因在所有细胞中的平均表达量。所有细胞中平均分数最高n_top的基因被绘制为箱形图。

3. 质量控制

然后进行基本的过滤（质量控制），使用两个工具：

• sc.pp.filter_cells进行细胞的过滤，该函数保留至少有 min_genes 个基因（某个基因表达非0可判断存在该基因）的细胞，或者保留至多有 max_genes 个基因的细胞；
• sc.pp.filter_genes进行基因的过滤，该函数用于保留在至少 min_cells 个细胞中出现的基因，或者保留在至多 max_cells 个细胞中出现的基因；

# 基因表达低于200的细胞将要删除
sc.pp.filter_cells(adata, min_genes=200)
# 至少 3 个细胞中检测到表达的基因才会被保留下来
sc.pp.filter_genes(adata, min_cells=3)

print(adata)

AnnData object with n_obs × n_vars = 2700 × 13714
    obs: 'n_genes'
    var: 'gene_ids', 'n_cells'

# 稀疏矩阵通常用于表示高维数据，例如基因表达数据，其中大多数值都是零
print(adata.X)
# 结果如下：
(0, 29)		1.0
(0, 73)		1.0
(0, 80)		2.0
(0, 148)	1.0
(0, 163)	1.0
(0, 184)	1.0

print(adata.var)
# 结果如下：
                      gene_ids  n_cells
AL627309.1     ENSG00000237683        9
AP006222.2     ENSG00000228463        3
RP11-206L10.2  ENSG00000228327        5
RP11-206L10.9  ENSG00000237491        3
LINC00115      ENSG00000225880       18

稀疏矩阵中，每个元素由三个值组成：(i, j, value)。其中，i 表示行索引，j 表示列索引，而 value 表示在索引为 (i, j) 的位置上的值。在这个例子中，adata.X 返回的稀疏矩阵包含了多个非零元素。每一行代表一个样本或数据点，每一列代表一个特征或基因。

adata.var 是一个 DataFrame，它包含两列：gene_ids 和 n_cells。

• gene_ids 列包含基因的标识符或 ID，每行对应于一个基因。
• n_cells 列包含每个基因在数据集中出现的细胞数目，即在多少个细胞中检测到了该基因。

通过查看 adata.var，你可以获得关于数据集中基因的一些信息，比如它们的标识符以及它们在样本中的表达情况。

---

3.1 质控选做

下一步是过滤线粒体核糖体基因（质量控制的选做步骤）：这是一个很难把握的工作，需要结合自己项目的情况来做。不过通常有以下策略：

• 粗暴去除所有线粒体核糖体基因，直接去除包含”MT-”开头的基因。
• 选择阈值去除高表达量的细胞，阈值很大程度上取决于对自己项目的了解程度，因为不同器官组织提取的单细胞，线粒体基因平均水平不一样。

使用pp.calculate_qc_metrics，我们可以高效计算很多度量指标：

# 将 adata.var_names 列中以 "MT-" 开头的元素赋值为 True，并将其保存在 adata.var  Dataframe 的 mt 列中。
adata.var['mt'] = adata.var_names.str.startswith('MT-')
adata.var['mt']
"""
AL627309.1       False
                 ...  
SRSF10-1         False
Name: mt, Length: 13714, dtype: bool
"""

# 计算