有人可能會說：單細(xì)胞分析使用Seurat频伤，monocle等R包會更加方便。但是實際分析中,測試情況是當(dāng)細(xì)胞量大于5萬時维咸。一般小型服務(wù)器內(nèi)存很容易不足剂买，這時候請不要過多嘗試使用Seurat來進(jìn)行分析，monocle2更是癌蓖。而基于python的單細(xì)胞轉(zhuǎn)錄分析包scanpy瞬哼，能很好的解決內(nèi)存不足的問題，親測整合80萬細(xì)胞量時租副，整個預(yù)處理流程在6小時左右能夠完成坐慰。
scanpy相關(guān)python 包安裝（安裝好python3之后，終端運行）用僧，一定要注意的是结胀，這里的python最好不要是3.9版本往上的，否則涉及多單細(xì)胞數(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```也無濟(jì)于事
#我個人建議在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... 非常好用,最好是通過文獻(xiàn)看一下單細(xì)胞的開發(fā)環(huán)境颠放，然后把他們復(fù)制過來。
#然后也可以看一下自己的conda環(huán)境下pip list都是哪些版本吭敢，也可以進(jìn)行移植碰凶。
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é)果顯示的詳細(xì)程度，數(shù)字越大越詳細(xì)畜晰。
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
#細(xì)胞名和基因名可分別通過 data.obs_names 和 data.var_names 查看。 AnnData 對象可以像 dataframe 一樣進(jìn)行切片操作盏袄，例如降铸，data_subset = data[:, list_of_gene_names]
"""
AnnData各部分?jǐn)?shù)據(jù)類型
| 功能 | 數(shù)據(jù)類型 | |
|---|---|---|
| data.X | 矩陣數(shù)據(jù) | numpy,scipy sparse,matrix |
| data.obs | 觀察值數(shù)據(jù) | pandas dataframe |
| data.var | 特征和高變基因數(shù)據(jù) | pandas dataframe |
| data.uns | 非結(jié)構(gòu)化數(shù)據(jù) | dict |
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進(jìn)行細(xì)胞的過濾政敢，該函數(shù)保留至少有 min_genes 個基因（某個基因表達(dá)非0可判斷存在該基因）的細(xì)胞其徙，或者保留至多有 max_genes 個基因的細(xì)胞；
sc.pp.filter_genes進(jìn)行基因的過濾喷户，該函數(shù)用于保留在至少 min_cells 個細(xì)胞中出現(xiàn)的基因唾那，或者保留在至多 max_cells 個細(xì)胞中出現(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)：
*選擇閾值去除高表達(dá)量的細(xì)胞褪尝，閾值很大程度上取決于對自己項目的了解程度闹获，因為不同器官組織提取的單細(xì)胞，線粒體基因平均水平不一樣河哑。
## 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'
"""
# 繪制高表達(dá)的前20個基因
sc.pl.highest_expr_genes(data, n_top=20, save='_pbmc3k.png')
"""
使用violinplot度量以下質(zhì)量：
n_genes_bt_counts：每個細(xì)胞中，有表達(dá)的基因的個數(shù)羡滑；
total_counts：每個細(xì)胞的基因總計數(shù)（總表達(dá)量,umi數(shù)）菇爪；
pct_counts_mt：每個細(xì)胞中，線粒體基因表達(dá)量占該細(xì)胞所有基因表達(dá)量的百分比
pct_counts_hb:每個細(xì)胞中柒昏，血紅蛋白基因表達(dá)量占該細(xì)胞所有基因表達(dá)量的百分比
pct_counts_ribo:每個細(xì)胞中凳宙，核糖體RNA基因表達(dá)量占該細(xì)胞所有基因表達(dá)量的百分比
"""
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)
# 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 基因的表達(dá)水平被認(rèn)為主要是技術(shù)性的，因此除了計算每個細(xì)胞中線粒體基因职祷，血紅蛋白基因和核糖體基因所占的比例外氏涩，明智的做法是還要將線粒體基因，血紅蛋白基因,MALAT1基因從數(shù)據(jù)集中直接刪除有梆，然后再進(jìn)行任何進(jìn)一步的分析
#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ù)量再進(jìn)行過濾,對data進(jìn)行切片(類似于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'
"""
我們先把矩陣提取出來看一下表達(dá)量的值
mat = pd.DataFrame(data=test.X.todense(),index=test.obs_names,columns=test.var_names)
mat
注意原始值都是整數(shù)
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)化 是尖，以便細(xì)胞之間的基因表達(dá)量具有可比性
sc.pp.normalize_total(data, target_sum=1e4)
## normalize with total UMI count per cell
函數(shù)normalize_total可以對每個細(xì)胞進(jìn)行標(biāo)準(zhǔn)化，以便每個細(xì)胞在標(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é)果
沒取標(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)
5. 高變基因選冉刃凇：
高變異基因就是highly variable features（HVGs），就是在細(xì)胞與細(xì)胞間進(jìn)行比較痰催，選擇表達(dá)量差別最大的基因兜辞，Seurat使用FindVariableFeatures函數(shù)鑒定高變基因，這些基因在不同細(xì)胞之間的表達(dá)量差異很大（在一些細(xì)胞中高表達(dá)夸溶，在另一些細(xì)胞中低表達(dá)）逸吵。默認(rèn)情況下，會返回2,000個高變基因用于下游的分析缝裁。
利用FindVariableFeatures函數(shù)胁塞，會計算一個mean-variance結(jié)果，也就是給出表達(dá)量均值和方差的關(guān)系并且得到top variable features压语，這一步的目的是鑒定出細(xì)胞與細(xì)胞之間表達(dá)量相差很大的基因啸罢，用于后續(xù)鑒定細(xì)胞類型。
標(biāo)記基因 (marker gene)胎食，是一種已知功能或已知序列的基因扰才，能夠起著特異性標(biāo)記的作用。
標(biāo)記基因通常是高變基因中很小的子集厕怜；
#識別高度可變基因衩匣，并進(jìn)行過濾：
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']]
6. 歸一化(將數(shù)據(jù)縮放到單位方差)：
>>>將數(shù)據(jù)data.X縮放到單位方差和零均值蕾总，對于縮放后的數(shù)據(jù)，在值為10處截斷：
""" sc.pp.regress_out(data, keys) keys:要回歸的data.obs中的數(shù)據(jù)列 """
# 回歸 adata.obs 中的列 (columns)
# 回歸每個細(xì)胞的總計數(shù)和表達(dá)的線粒體基因的百分比的影響琅捏。
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
如果你這時候想看矩陣情況蚀浆，需要進(jìn)行轉(zhuǎn)換
s = scipy.sparse.csr_matrix(data.X)
mat = pd.DataFrame(s.todense(),index = data.obs_names,colums = data.var_names)
7. 主成分分析：
通過運行主成分分析(PCA)來降低數(shù)據(jù)的維數(shù)，該分析揭示了變化的主軸并對數(shù)據(jù)進(jìn)行去噪搜吧。
#svd_solver：指定奇異值分解SVD的方法市俊，有4個可以選擇的值：{‘a(chǎn)uto’, ‘full’, ‘a(chǎn)rpack’, ‘randomized’，}
#如果設(shè)為’arpack’滤奈，則使用scipy.sparse.linalg.svds計算SVD分解摆昧。這種方法嚴(yán)格要求 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ù)每個散點（細(xì)胞）中包含基因CST3的數(shù)量為散點標(biāo)記顏色。
計算單個pc對數(shù)據(jù)總方差的貢獻(xiàn),這可以提供給我們應(yīng)該考慮多少個 PC 以計算細(xì)胞的鄰域關(guān)系的信息(resolution選擇多少)昭躺，例如用于后續(xù)的聚類函數(shù) sc.tl.louvain()
sc.pl.pca_variance_ratio(data,log = True)
8.降維（對neighborhood graph進(jìn)行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表示來計算細(xì)胞的鄰域圖±袷猓可以在這里簡單地使用默認(rèn)值驹吮。為了可視化聚類的結(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 是指字體粗細(xì)
#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基因查找
通過文獻(xiàn)或cellmarker查找marker 基因
#先設(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)
11. 通過差異表達(dá)基因?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軸范圍。
