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

单细胞数据分析常规流程

面对高效快速的要求上，利用R分析数据越来越困难，转战Python分析，我们通过scanpy官网去学习如何分析单细胞卑鄙常规分析。
数据3k PBMC来自健康的志愿者，可从10x Genomics免费获得。在linux体系上，可以取消表明并运行以下操作来下载和解压缩数据。最后一行创建一个用于保存已处理惩罚数据的目录write，后面直接利用保存的数据，能快速加载数据。
下载数据：

1. $mkdir data
2. $cd data
3. $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
4. $tar -xzf pbmc3k_filtered_gene_bc_matrices.tar.gz
5. # 获得数据

复制代码

1. 数据加载

1. import numpy as np
2. import pandas as pd
3. import scanpy as sc
5. sc.settings.verbosity = 3             # verbosity: errors (0), warnings (1), info (2), hints (3)
6. sc.logging.print_header()
7. sc.settings.set_figure_params(dpi=80, facecolor='white')
9. # 声明h5ad用于存储分析结果
10. results_file = 'data/write/pbmc3k.h5ad'
12. adata = sc.read_10x_mtx(
13.     'data/filtered_gene_bc_matrices/hg19/',  # `.mtx`文件所在的目录
14.     var_names='gene_symbols',                # 用 gene 作为var
15.     cache=True) # 开启缓存读写
17. """
18. 注意cache=Trure
19. ... writing an h5ad cache file to speedup reading next time
20. 下次读取就不会从count matrix读, 会直接从cache目录下的h5ad文件读(更快)
21. """

复制代码

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

• gene_ids：利用基因的唯一标识符作为变量名。这通常是一种更确切和唯一的标识，差别基因之间不存在重复。利用基因的唯一标识符作为变量名可以确保在分析中每个基因都有唯一的标识符，并且不会出现混淆或重复。
  
  
• gene_symbols：利用基因的符号名称作为变量名。基因的符号名称通常更容易明白和记忆，因为它们通常是基于基因的功能或特征而定名的。然而，基因的符号名称不一定是唯一的，大概存在多个基因具有相同的符号名称，这大概会导致一些混淆或不一致。
  

因此，你可以根据详细的需求和分析的目的来选择利用哪种类型的变量名。如果需要确保每个基因都具有唯一的标识符，并且不会出现混淆或重复，那么可以利用 gene_ids。如果更关注基因的功能或特征，并且不太担心大概存在的重复符号名称，那么可以利用 gene_symbols。
<hr> 留意，如果在函数sc.read_10x_mtx中指定参数var_names='gene_ids'时，下一个操作将是不必要的：

1. # 消除重复的列
2. adata.var_names_make_unique()
4. print(adata)
6. AnnData object with n_obs × n_vars = 2700 × 32738
7.     var: 'gene_ids'

复制代码

adata包罗2700个细胞、32738个基因的对象
2. top基因箱型图

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

1. 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 个细胞中出现的基因；
  

1. # 基因表达低于200的细胞将要删除
2. sc.pp.filter_cells(adata, min_genes=200)
3. # 至少 3 个细胞中检测到表达的基因才会被保留下来
4. sc.pp.filter_genes(adata, min_cells=3)
6. print(adata)
8. AnnData object with n_obs × n_vars = 2700 × 13714
9.     obs: 'n_genes'
10.     var: 'gene_ids', 'n_cells'

复制代码

1. # 稀疏矩阵通常用于表示高维数据，例如基因表达数据，其中大多数值都是零
2. print(adata.X)
3. # 结果如下：
4. (0, 29)                1.0
5. (0, 73)                1.0
6. (0, 80)                2.0
7. (0, 148)        1.0
8. (0, 163)        1.0
9. (0, 184)        1.0
11. print(adata.var)
12. # 结果如下：
13.                       gene_ids  n_cells
14. AL627309.1     ENSG00000237683        9
15. AP006222.2     ENSG00000228463        3
16. RP11-206L10.2  ENSG00000228327        5
17. RP11-206L10.9  ENSG00000237491        3
18. 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，你可以获得关于数据集中基因的一些信息，比如它们的标识符以及它们在样本中的表达情况。
<hr> 3.1 质控选做

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

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

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

1. # 将 adata.var_names 列中以 "MT-" 开头的元素赋值为 True，并将其保存在 adata.var  Dataframe 的 mt 列中。
2. adata.var['mt'] = adata.var_names.str.startswith('MT-')
3. adata.var['mt']
4. """
5. AL627309.1       False
6.                  ...  
7. SRSF10-1         False
8. Name: mt, Length: 13714, dtype: bool
9. """
11. # 计算指标

复制代码

免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。