有人可能會說：單細胞分析使用Seurat，monocle等R包會更加方便。但是實際分析中,測試情況是當(dāng)細胞量大于5萬時。一般小型服務(wù)器內(nèi)存很容易不足，這時候請不要過多嘗試使用Seurat來進行分析，monocle2更是。而基于python的單細胞轉(zhuǎn)錄分析包scanpy，能很好的解決內(nèi)存不足的問題，親測整合80萬細胞量時，整個預(yù)處理流程在6小時左右能夠完成。
scanpy相關(guān)python 包安裝（安裝好python3之后，終端運行），一定要注意的是，這里的python最好不要是3.9版本往上的，否則涉及多單細胞數(shù)據(jù)整合加載bbknn包會報錯，numba將無法與當(dāng)前環(huán)境適配。

pip install -i https://pypi.doubanio.com/simple/ scanpy==1.6.1 
##注意要加上-i參數(shù)，給pip加上一個豆瓣或者清華(https://pypi.tuna.tsinghua.edu.cn/simple)的鏡像，否則下載起來你即使加了````--default-timeout=1000```也無濟于事
#我個人建議在conda中新建一個環(huán)境，執(zhí)行：
conda create -n scrna_test python=3.7.0 numba=0.55.2 pandas=1.1.5 llvmlite=0.38.1 numpy=1.21.6 bbknn=1.5.1
#然后在scanpy官網(wǎng)上下載scanpy-1.9.1-py3-none-any.whl（https://pypi.org/project/scanpy/#files）
wget https://files.pythonhosted.org/packages/51/87/a55c7992cba9b189de70eae37e9f1e2abe6fdaf3f087d30356f28698948e/scanpy-1.9.1-py3-none-any.whl
#下載好后
pip install scanpy-1.9.1-py3-none-any.whl 
#下載scanpy非常的困難，因為他對numba和 llvmlite有版本要求
#如果你在安裝scanpy沒有安裝好pandas, numpy,numba,  llvmlite,事情會變得非常復(fù)雜，報錯一個接一個。
#conda env create -n env_name pakage1=version package2=version... 非常好用,最好是通過文獻看一下單細胞的開發(fā)環(huán)境，然后把他們復(fù)制過來。
#然后也可以看一下自己的conda環(huán)境下pip list都是哪些版本，也可以進行移植。

1. 運行python3，導(dǎo)入相關(guān)包及設(shè)置一些必要路徑：

import scanpy as sc
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
sc.settings.verbosity = 3 # verbosity即冗余 。設(shè)置日志等級: errors (0), warnings (1), info (2), hints (3)，取值表示測試日志結(jié)果顯示的詳細程度，數(shù)字越大越詳細。
sc.logging.print_versions() # 輸出版本號
sc.settings.set_figure_params(dpi=80) # set_figure_params 設(shè)置圖片的分辨率/大小以及其他樣式
import os #在服務(wù)器運行，習(xí)慣性會設(shè)置一個輸出路徑，用于保存pdf圖片
os.getcwd()  #查看當(dāng)前路徑
os.chdir('./filtered_gene_bc_matrices/scanpy') #修改路徑
os.getcwd()
results_file = 'pbmc3k.h5ad' #聲明什么用于儲存分析結(jié)果

2. 讀取并查看數(shù)據(jù)：

scanpy可以讀取.h5 .h5ad 以及10X標(biāo)準(zhǔn)化(features.tsv,barcodes.tsv, matrix.mtx)格式數(shù)據(jù)

data = sc.read_10x_h5("./pbmc3K.h5",genome=None, gex_only=True)
print(data)
AnnData object with n_obs × n_vars = 2700 × 32738
    var: 'gene_ids'
"""
#cell數(shù)700 基因數(shù)32738的矩陣
#也可以是data = sc.read_10x_matrix('path',var_names = 'gene_symbols', cache = True)
#使用gene_symbols作為AnnData的特征名稱
#cache=True , cache為True表示寫入緩存(cache) 便于快速讀寫
#如果是h5ad格式，可以直接read('/.h5ad')
#data.X 存儲count matrix，數(shù)據(jù)類型為稀疏矩陣scipy.sparse.csr.csr_matrix 
#data.obs 存儲關(guān)于 obervations(cells) 的 metadata，數(shù)據(jù)類型為 pandas dataframe 
#data.var 存儲關(guān)于 variables(genes) 的 metadata，數(shù)據(jù)類型為  pandas dataframe 
#AnnData.uns 存儲后續(xù)附加的其他非結(jié)構(gòu)化信息 
#data.obs_names 和 data.var_names是 index 
#細胞名和基因名可分別通過 data.obs_names 和 data.var_names 查看。 AnnData 對象可以像 dataframe 一樣進行切片操作，例如，data_subset = data[:, list_of_gene_names]
"""

AnnData數(shù)據(jù)結(jié)構(gòu)

AnnData各部分數(shù)據(jù)類型

3. 索引去重：

data.var_names_make_unique()
# # 如果在 sc.read_10x_mtx 中使用了var_names='gene_ids'，這一步就是不必要的

4. 質(zhì)量控制：

質(zhì)量控制3個指標(biāo)：測到的轉(zhuǎn)錄本總數(shù)(total_counts)、測到的基因總數(shù)(total_cells)、來源于線粒體基因的轉(zhuǎn)錄本所占比例。

4.1基本過濾

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

sc.pp.filter_cells(data,min_genes=200)
sc.pp.filter_genes(data,min_cells=3)
data
"""
AnnData object with n_obs × n_vars = 2700 × 13714
    obs: 'n_genes'
    var: 'gene_ids', 'n_cells'
"""

4.2 計算質(zhì)量控制指標(biāo)：

*選擇閾值去除高表達量的細胞，閾值很大程度上取決于對自己項目的了解程度，因為不同器官組織提取的單細胞，線粒體基因平均水平不一樣。

## startswith()方法用于檢查字符串是否是以指定子字符串開頭, 如果是則返回 True, 否則返回 False
# mitochondrial genes
data.var['mt'] = data.var_names.str.startswitch('MT-')  
# hemoglobin genes.  血紅蛋白基因
data.var['hb'] = data.var_names.str.contains('^HB[^P]')
# ribosomal genes
data.var['ribo'] = data.var_names.str.startswitch('RPS','RPL')
"""
AL627309.1       False
                 ...  
SRSF10-1         False
Name: mt, Length: 13714, dtype: bool
"""

sc.pp.calculate_qc_metrics(data, qc_vars=['mt',‘hb’,'ribo'], percent_top=None, log1p=False, inplace=True)
print(data)
#data：Anndata；
#qc_vars：用于標(biāo)識要控制的特征（基因），布爾型元素，用于作為mask使用；
#percent_top：計算與常出現(xiàn)基因的比，percent_top=[50] 計算與第 50 個最常出現(xiàn)基因的比例，None則不計算；
#inplace：決定是否將計算指標(biāo)添加到var和obs中；
#log1p：設(shè)置為False可以跳過轉(zhuǎn)換到log1p空間的過程；log1p即log(1+number)，用于壓縮數(shù)據(jù)并確保結(jié)果是一個正數(shù)；
"""
注釋中多了很多QC計算得到的信息
AnnData object with n_obs × n_vars = 2700 × 13714
    obs: 'n_genes', 'nFeature_RNA', nCount_RNA', 'total_counts_mt', 'pct_counts_mt'
    var: 'gene_ids', 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts'
"""

# 繪制高表達的前20個基因
sc.pl.highest_expr_genes(data, n_top=20, save='_pbmc3k.png')

n_top=20

"""
使用violinplot度量以下質(zhì)量：
n_genes_bt_counts：每個細胞中，有表達的基因的個數(shù)；
total_counts：每個細胞的基因總計數(shù)（總表達量,umi數(shù)）；
pct_counts_mt：每個細胞中，線粒體基因表達量占該細胞所有基因表達量的百分比
pct_counts_hb:每個細胞中，血紅蛋白基因表達量占該細胞所有基因表達量的百分比
pct_counts_ribo:每個細胞中，核糖體RNA基因表達量占該細胞所有基因表達量的百分比
"""
sc.pl.violin(data, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt','pct_counts_ribo', 'pct_counts_hb'],
             jitter=0.4,  multi_panel = True)

1664981487608.png

# filter for percent mito
data = data[data.obs['pct_counts_mt'] < 20, :]
# filter for percent ribo > 0.05
data = data[data.obs['pct_count_ribo'] > 5,:]

由于線粒體和 MALAT1 基因的表達水平被認為主要是技術(shù)性的，因此除了計算每個細胞中線粒體基因，血紅蛋白基因和核糖體基因所占的比例外，明智的做法是還要將線粒體基因，血紅蛋白基因,MALAT1基因從數(shù)據(jù)集中直接刪除，然后再進行任何進一步的分析
#delete var_names.str.startswitch['MT-'] 
#var_names.str.startswitch('MALAT1') var_names.str.contains('^HB[^(P)]')
mito_genes = data.var_names.str.startswith('MT-')
malat1 = data.var_names.str.startswith('MALAT1')
hb_genes = data.var_names.str.contains('^HB[^(P)]')
remove = np.add(mito_genes, malat1)
remove = np.add(remove, hb_genes)
keep = np.invert(remove)
data = data[:,keep]

根據(jù)基因數(shù)量再進行過濾,對data進行切片(類似于Seurat中的subset)

data = data[data.obs['n_genes_by_counts'] < 2500, :]
data = data[data.obs['n_genes_by_counts'] > 500,:]
"""
View of AnnData object with n_obs × n_vars = 2638 × 13714
    obs: 'n_genes', 'nFeature_RNA', nCount_RNA', 'total_counts_mt', 'pct_counts_mt'
    var: 'gene_ids', 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts'
"""

我們先把矩陣提取出來看一下表達量的值

mat = pd.DataFrame(data=test.X.todense(),index=test.obs_names,columns=test.var_names)
mat

注意原始值都是整數(shù)

image.png

5. 數(shù)據(jù)標(biāo)準(zhǔn)化：

Normalize each cell by total counts over all genes, so that every cell has the same total count after normalization.
總計數(shù)標(biāo)準(zhǔn)化 ，以便細胞之間的基因表達量具有可比性

sc.pp.normalize_total(data, target_sum=1e4)
## normalize with total UMI count per cell

函數(shù)normalize_total可以對每個細胞進行標(biāo)準(zhǔn)化，以便每個細胞在標(biāo)準(zhǔn)化后沿著基因方向求和具有相同的總數(shù) ````
target_sum：

data.X
array([[ 3.,  3.,  3.,  6.,  6.],
       [ 1.,  1.,  1.,  2.,  2.],
       [ 1., 22.,  1.,  2.,  2.]], dtype=float32)
# 設(shè)置 target_sum=1 標(biāo)準(zhǔn)化后
X_norm
array([[0.14, 0.14, 0.14, 0.29, 0.29],
       [0.14, 0.14, 0.14, 0.29, 0.29],
       [0.04, 0.79, 0.04, 0.07, 0.07]], dtype=float32)

看一下做完標(biāo)準(zhǔn)化的結(jié)果

image.png

沒取標(biāo)準(zhǔn)化之前：

View of AnnData object with n_obs × n_vars = 15072 × 21655
    obs: 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'total_counts_hb', 'pct_counts_hb', 'total_counts_ribo', 'pct_counts_ribo', 'doublet_score', 'predicted_doublet'
    var: 'n_cells', 'mt', 'hb', 'ribo', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts'
    uns: 'scrublet'

標(biāo)準(zhǔn)化之后：

AnnData object with n_obs × n_vars = 15072 × 21655 
obs: 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'total_counts_hb', 'pct_counts_hb', 'total_counts_ribo', 'pct_counts_ribo', 'doublet_score', 'predicted_doublet' 
var: 'n_cells', 'mt', 'hb', 'ribo', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts' 
uns: 'scrublet'

標(biāo)準(zhǔn)化前后只是改變了data.X 沒有增加obs或者var的內(nèi)容。

#將數(shù)據(jù)壓縮到log1p的空間
sc.pp.log1p(data)

log1p 也就是log(X+1) 也就是In(X+1)

image.png

image.png

5. 高變基因選取：

高變異基因就是highly variable features（HVGs），就是在細胞與細胞間進行比較，選擇表達量差別最大的基因，Seurat使用FindVariableFeatures函數(shù)鑒定高變基因，這些基因在不同細胞之間的表達量差異很大（在一些細胞中高表達，在另一些細胞中低表達）。默認情況下，會返回2,000個高變基因用于下游的分析。
利用FindVariableFeatures函數(shù)，會計算一個mean-variance結(jié)果，也就是給出表達量均值和方差的關(guān)系并且得到top variable features，這一步的目的是鑒定出細胞與細胞之間表達量相差很大的基因，用于后續(xù)鑒定細胞類型。
標(biāo)記基因 (marker gene)，是一種已知功能或已知序列的基因，能夠起著特異性標(biāo)記的作用。
標(biāo)記基因通常是高變基因中很小的子集；

#識別高度可變基因，并進行過濾：
sc.pp.highly_variable_genes(data,min_mean = 0.0125,max_mean=3,min_disp=0.5)
sc.pl.highly_variable_genes(data)
#保存原始數(shù)據(jù)后再把data變成只有高變基因
data.raw = data
#過濾
##注意切片data[obs:var]
data = data[:,data.var['highly_variable']]

image.png

6. 歸一化(將數(shù)據(jù)縮放到單位方差)：

>>>將數(shù)據(jù)data.X縮放到單位方差和零均值，對于縮放后的數(shù)據(jù)，在值為10處截斷：

""" sc.pp.regress_out(data, keys) keys:要回歸的data.obs中的數(shù)據(jù)列 """
# 回歸 adata.obs 中的列 (columns) 
# 回歸每個細胞的總計數(shù)和表達的線粒體基因的百分比的影響。
sc.pp.regress_out(data, ['total_counts', 'pct_counts_mt'])
# 將每個基因縮放到單位方差。
sc.pp.scale(data,max_value = 10)
#保存數(shù)據(jù)
data.write(results_file)

需要注意的是，在做完歸一化后，data.X的數(shù)據(jù)格式從scipy.sparse.csr_matrix轉(zhuǎn)換為numpy.ndarray

如果你這時候想看矩陣情況，需要進行轉(zhuǎn)換

s = scipy.sparse.csr_matrix(data.X)
mat = pd.DataFrame(s.todense(),index = data.obs_names,colums = data.var_names)

image.png

7. 主成分分析：

通過運行主成分分析(PCA)來降低數(shù)據(jù)的維數(shù)，該分析揭示了變化的主軸并對數(shù)據(jù)進行去噪。

#svd_solver：指定奇異值分解SVD的方法，有4個可以選擇的值：{‘a(chǎn)uto’, ‘full’, ‘a(chǎn)rpack’, ‘randomized’，}
#如果設(shè)為’arpack’，則使用scipy.sparse.linalg.svds計算SVD分解。這種方法嚴格要求 0 < n_components < min(樣本數(shù)，特征數(shù))。
sc.tl.pca(data,svd_solver = 'arpack')
sc.pl.pca(data, color = 'CST3')

PCA將高維數(shù)據(jù)data.X聚類到2維空間后得到的只是一些平面下的散點，但我們可以根據(jù)每個散點（細胞）中包含基因CST3的數(shù)量為散點標(biāo)記顏色。

計算單個pc對數(shù)據(jù)總方差的貢獻,這可以提供給我們應(yīng)該考慮多少個 PC 以計算細胞的鄰域關(guān)系的信息(resolution選擇多少)，例如用于后續(xù)的聚類函數(shù) sc.tl.louvain()

sc.pl.pca_variance_ratio(data,log = True)

8.降維（對neighborhood graph進行embedding）

sc.pp.neighbors(data,n_neighbors=10,n_pcs=25)
"""
computing neighbors
    using 'X_pca' with n_pcs = 40
    finished: added to `.uns['neighbors']`
    `.obsp['distances']`, distances for each pair of neighbors
    `.obsp['connectivities']`, weighted adjacency matrix (0:00:05)
"""
sc.tl.umap(data)

9. 聚類

scanpy提供leiden和louvain兩種圖聚類算法，圖聚類在開始之前要先找到鄰居。
首先計算neighborhood graph，我們使用數(shù)據(jù)矩陣的PCA表示來計算細胞的鄰域圖?？梢栽谶@里簡單地使用默認值。為了可視化聚類的結(jié)果我們做一下umap降維，看看分群結(jié)果。

sc.tl.leiden(data,resolution = 0.3)
#用umap可視化聚類結(jié)果
sc.pl.umap(data,color = 'leiden', frameon=False, title='',
           legend_loc='right margin', legend_fontsize=12,
           legend_fontweight='normal', legend_fontoutline=6,
           palette=None, cmap=None)
#platelet是指給clutesr指定一組顏色。例如color=['red','blue'....]
#cmap是指顏色選擇
#legend_fontweight 是指字體粗細
#legend_loc='on data'
#legend_fontoutline:圖例字體輪廓的線寬，單位為pt。使用withStroke的路徑效果繪制白色輪廓。
#. 當(dāng)frameon=True的時候, 圖示會被繪制在一個patch實體上;
# 否則, 如果frameon=False, 則圖示會被直接繪制在圖片上。
#這里, 討論是否將圖示繪制在一個patch實體上的意義在于,
#當(dāng)把它繪制在一個patch實體上時, 
#我們才可以使用facecolor, edgecolor, framealpha, fancybox等參數(shù)來設(shè)置圖示的背景(不是圖片的背景)的顏色, 邊框顏色, 透明度, 以及形狀, 而當(dāng)frameon=False的時候這些參數(shù)就會失效

10. Marker基因查找

通過文獻或cellmarker查找marker 基因

PBMC common marker genes

#先設(shè)置分辨率以及長寬
sc.settings.set_figure_params(dpi=50, dpi_save=600, figsize=(5,5))
marker_names = ['IL7R','CCR7',
                'CD14','LYZ',
               'IL7R','S100A4',
               'MS4A1',
               'CD8A',
               'FCGR3A','MS4A7',
               'GNLY','NKG7',
               'FCER1A','CST3',
               'PPBP']
sc.pl.umap(data,color = marker_names, ncols=2)

1665051072562.png

11. 通過差異表達基因?qū)ふ襇arker基因

讓我們計算每個cluster中高度差異基因的排名，最簡單和最快的方法是t-test，當(dāng)然還有wilcoxon。

sc.settings.set_figure_params(dpi=80,dpi_save=600,figsize=(10,10))
sc.tl.rank_genes_groups(data,'leiden',method='t-test')
#繪圖指定每個cluster前多少個基因，每個cluster之間是否共享
sc.pl.rank_genes_groups(data,n_genes=25,sharey=False,fontsize=15)

#sharely ： 控制是否應(yīng)共享每個panel的y軸，sharey =False表示每個panel都有自己的y軸范圍。

image.png