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

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

數(shù)據(jù)3k PBMC來自健康的志愿者，可從10x Genomics免費(fèi)獲得。在linux系統(tǒng)上，可以取消注釋并運(yùn)行以下操作來下載和解壓縮數(shù)據(jù)。最后一行創(chuàng)建一個(gè)用于保存已處理數(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用于存儲(chǔ)分析結(jié)果
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
下次讀取就不會(huì)從count matrix讀, 會(huì)直接從cache目錄下的h5ad文件讀(更快)
"""

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

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

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

2. gene_symbols：使用基因的符號(hào)名稱作為變量名?；虻姆?hào)名稱通常更容易理解和記憶，因?yàn)樗鼈兺ǔＪ腔诨虻墓δ芑蛱卣鞫?。然而，基因的符?hào)名稱不一定是唯一的，可能存在多個(gè)基因具有相同的符號(hào)名稱，這可能會(huì)導(dǎo)致一些混淆或不一致。

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

---

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

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

adata包含2700個(gè)細(xì)胞、32738個(gè)基因的對(duì)象

2. top基因箱型圖

下圖計(jì)算每一個(gè)基因在所有細(xì)胞中的平均表達(dá)量，并繪制了平均表達(dá)量前30的基因箱型圖。

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

計(jì)算每一個(gè)基因在所有細(xì)胞中的平均表達(dá)量。所有細(xì)胞中平均分?jǐn)?shù)最高n_top的基因被繪制為箱形圖。

3. 質(zhì)量控制

然后進(jìn)行基本的過濾（質(zhì)量控制），使用兩個(gè)工具：

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

# 基因表達(dá)低于200的細(xì)胞將要?jiǎng)h除
sc.pp.filter_cells(adata, min_genes=200)
# 至少 3 個(gè)細(xì)胞中檢測(cè)到表達(dá)的基因才會(huì)被保留下來
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ù)，例如基因表達(dá)數(shù)據(jù)，其中大多數(shù)值都是零
print(adata.X)
# 結(jié)果如下：
(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)
# 結(jié)果如下：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

稀疏矩陣中，每個(gè)元素由三個(gè)值組成：(i, j, value)。其中，i 表示行索引，j 表示列索引，而 value 表示在索引為 (i, j) 的位置上的值。在這個(gè)例子中，adata.X 返回的稀疏矩陣包含了多個(gè)非零元素。每一行代表一個(gè)樣本或數(shù)據(jù)點(diǎn)，每一列代表一個(gè)特征或基因。

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

• gene_ids 列包含基因的標(biāo)識(shí)符或 ID，每行對(duì)應(yīng)于一個(gè)基因。
• n_cells 列包含每個(gè)基因在數(shù)據(jù)集中出現(xiàn)的細(xì)胞數(shù)目，即在多少個(gè)細(xì)胞中檢測(cè)到了該基因。

通過查看 adata.var，你可以獲得關(guān)于數(shù)據(jù)集中基因的一些信息，比如它們的標(biāo)識(shí)符以及它們?cè)跇颖局械谋磉_(dá)情況。

---

3.1 質(zhì)控選做

下一步是過濾線粒體核糖體基因（質(zhì)量控制的選做步驟）：這是一個(gè)很難把握的工作，需要結(jié)合自己項(xiàng)目的情況來做。不過通常有以下策略：

• 粗暴去除所有線粒體核糖體基因，直接去除包含”MT-”開頭的基因。
• 選擇閾值去除高表達(dá)量的細(xì)胞，閾值很大程度上取決于對(duì)自己項(xiàng)目的了解程度，因?yàn)椴煌鞴俳M織提取的單細(xì)胞，線粒體基因平均水平不一樣。

使用pp.calculate_qc_metrics，我們可以高效計(jì)算很多度量指標(biāo)：

# 將 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
"""# 計(jì)算指標(biāo)