單細胞數(shù)據(jù)分析常規(guī)流程

面對高效快速的要求上，使用R分析數(shù)據(jù)越來越困難，轉戰(zhàn)Python分析，我們通過scanpy官網(wǎng)去學習如何分析單細胞下游常規(guī)分析。

數(shù)據(jù)3k PBMC來自健康的志愿者，可從10x Genomics免費獲得。在linux系統(tǒng)上，可以取消注釋并運行以下操作來下載和解壓縮數(shù)據(jù)。最后一行創(chuàng)建一個用于保存已處理數(shù)據(jù)的目錄write，后面直接使用保存的數(shù)據(jù)，能快速加載數(shù)據(jù)。

下載數(shù)據(jù)：

$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
# 獲得數(shù)據(jù)

1. 數(shù)據(jù)加載

import numpy as np
import pandas as pd
import scanpy as scsc.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 作為varcache=True) # 開啟緩存讀寫"""
注意cache=Trure
... writing an h5ad cache file to speedup reading next time
下次讀取就不會從count matrix讀, 會直接從cache目錄下的h5ad文件讀(更快)
"""

在函數(shù) sc.read_10x_mtx 中，參數(shù) var_names 用于指定在加載數(shù)據(jù)時使用哪個變量來作為基因的名稱。在這里，如果你將 var_names='gene_ids'，它將使用基因的唯一標識符作為變量名，而如果你將 var_names='gene_symbols'，它將使用基因的符號名稱作為變量名。

這兩者之間的區(qū)別在于：

1. gene_ids：使用基因的唯一標識符作為變量名。這通常是一種更確切和唯一的標識，不同基因之間不存在重復。使用基因的唯一標識符作為變量名可以確保在分析中每個基因都有唯一的標識符，并且不會出現(xiàn)混淆或重復。

2. gene_symbols：使用基因的符號名稱作為變量名?；虻姆柮Q通常更容易理解和記憶，因為它們通常是基于基因的功能或特征而命名的。然而，基因的符號名稱不一定是唯一的，可能存在多個基因具有相同的符號名稱，這可能會導致一些混淆或不一致。

因此，你可以根據(jù)具體的需求和分析的目的來選擇使用哪種類型的變量名。如果需要確保每個基因都具有唯一的標識符，并且不會出現(xiàn)混淆或重復，那么可以使用 gene_ids。如果更關注基因的功能或特征，并且不太擔心可能存在的重復符號名稱，那么可以使用 gene_symbols。

---

注意，如果在函數(shù)sc.read_10x_mtx中指定參數(shù)var_names='gene_ids'時，下一個操作將是不必要的：

# 消除重復的列
adata.var_names_make_unique()print(adata)AnnData object with n_obs × n_vars = 2700 × 32738var: 'gene_ids'

adata包含2700個細胞、32738個基因的對象

2. top基因箱型圖

下圖計算每一個基因在所有細胞中的平均表達量，并繪制了平均表達量前30的基因箱型圖。

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

計算每一個基因在所有細胞中的平均表達量。所有細胞中平均分數(shù)最高n_top的基因被繪制為箱形圖。

3. 質量控制

然后進行基本的過濾（質量控制），使用兩個工具：

• sc.pp.filter_cells進行細胞的過濾，該函數(shù)保留至少有 min_genes 個基因（某個基因表達非0可判斷存在該基因）的細胞，或者保留至多有 max_genes 個基因的細胞；
• sc.pp.filter_genes進行基因的過濾，該函數(shù)用于保留在至少 min_cells 個細胞中出現(xiàn)的基因，或者保留在至多 max_cells 個細胞中出現(xiàn)的基因；

# 基因表達低于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 × 13714obs: 'n_genes'var: 'gene_ids', 'n_cells'

# 稀疏矩陣通常用于表示高維數(shù)據(jù)，例如基因表達數(shù)據(jù)，其中大多數(shù)值都是零
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.0print(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 返回的稀疏矩陣包含了多個非零元素。每一行代表一個樣本或數(shù)據(jù)點，每一列代表一個特征或基因。

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

• gene_ids 列包含基因的標識符或 ID，每行對應于一個基因。
• n_cells 列包含每個基因在數(shù)據(jù)集中出現(xiàn)的細胞數(shù)目，即在多少個細胞中檢測到了該基因。

通過查看 adata.var，你可以獲得關于數(shù)據(jù)集中基因的一些信息，比如它們的標識符以及它們在樣本中的表達情況。

---

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
"""# 計算指標