一谬莹、安裝
scanpy 單細胞分析包圖文詳解 01 | 深入理解 AnnData 數(shù)據(jù)結(jié)構(gòu)
pip install scanpy
conda install -y -c conda-forge leidenalg
二掷邦、使用
1钠右、準備工作
# 載入包
import numpy as np
import pandas as pd
import scanpy as sc
# 設(shè)置
sc.settings.verbosity = 3 # 設(shè)置日志等級: errors (0), warnings (1), info (2), hints (3)
sc.logging.print_header()
sc.settings.set_figure_params(dpi=80, facecolor='white')
# 用于存儲分析結(jié)果文件的路徑
results_file = 'write/pbmc3k.h5ad'
# 載入文件
adata = sc.read_10x_mtx(
'./filtered_gene_bc_matrices/hg19/', # mtx 文件目錄
var_names='gene_symbols', # 使用 gene_symbols 作為變量名
cache=True) # 寫入緩存，可以更快的讀取文件
2杀赢、預(yù)處理
顯示在所有細胞中在每個單細胞中產(chǎn)生最高計數(shù)分數(shù)的基因
sc.pl.highest_expr_genes(adata, n_top=20, )
過濾低質(zhì)量細胞樣本
過濾在少于三個細胞中表達，或一個細胞中表達少于200個基因的細胞樣本
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
過濾包含線粒體基因和表達基因過多的細胞
線粒體基因的轉(zhuǎn)錄本比單個轉(zhuǎn)錄物分子大，并且不太可能通過細胞膜逃逸骄酗。因此，檢測出高比例的線粒體基因悦冀，表明細胞質(zhì)量差（Islam et al. 2014; Ilicic et al. 2016）趋翻。
表達基因過多可能是由于一個油滴包裹多個細胞，從而檢測出比正常檢測要多的基因數(shù)盒蟆，因此要過濾這些細胞踏烙。
adata.var['mt'] = adata.var_names.str.startswith('MT-') # 將線粒體基因標記為 mt
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)
sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],
jitter=0.4, multi_panel=True)
生成的三張小提琴圖代表：表達基因的數(shù)量，每個細胞包含的表達量历等，線粒體基因表達量的百分比讨惩。
過濾
sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt')
sc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts')
過濾線粒體基因表達過多或總數(shù)過多的細胞，也就是紅框標識的樣本寒屯。
# 獲取線粒體基因占比在 5% 以下的細胞樣本
adata = adata[adata.obs.pct_counts_mt < 5, :]
# 獲取表達基因數(shù)在 2500 以下的細胞樣本
adata = adata[adata.obs.n_genes_by_counts < 2500, :]
3荐捻、檢測特異性基因
歸一化
# 歸一化，使得不同細胞樣本間可比
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
存儲數(shù)據(jù)
將 AnnData 對象的 .raw 屬性設(shè)置為歸一化和對數(shù)化的原始基因表達寡夹，以便以后用于基因表達的差異測試和可視化处面。這只是凍結(jié)了 AnnData 對象的狀態(tài)。
adata.raw = adata
識別特異性基因
# 計算
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
# 繪制特異性基因散點圖
sc.pl.highly_variable_genes(adata)
獲取只有特異性基因的數(shù)據(jù)集
# 獲取只有特異性基因的數(shù)據(jù)集
adata = adata[:, adata.var.highly_variable]
# 回歸每個細胞的總計數(shù)和表達的線粒體基因的百分比的影響菩掏。
sc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt'])
# 將每個基因縮放到單位方差魂角。閾值超過標準偏差 10。
sc.pp.scale(adata, max_value=10)
4患蹂、主成分分析（Principal component analysis）
通過運行主成分分析 (PCA) 來降低數(shù)據(jù)的維數(shù)或颊，可以對數(shù)據(jù)進行去噪并揭示不同分群的主因素砸紊。
# 繪制 PCA 圖
sc.pl.pca(adata, color='CST3')
檢查單個 PC 對數(shù)據(jù)總方差的貢獻，這可以提供給我們應(yīng)該考慮多少個 PC 以計算細胞的鄰域關(guān)系的信息囱挑，例如用于后續(xù)的聚類函數(shù) sc.tl.louvain() 或 tSNE 聚類 sc.tl.tsne()醉顽。
sc.pl.pca_variance_ratio(adata, log=True)
5、領(lǐng)域圖平挑，聚類圖（Neighborhood graph）
使用數(shù)據(jù)矩陣的 PCA 表示來計算單元格的鄰域圖游添。為了重現(xiàn) Seurat 的結(jié)果，我們采用以下值通熄。
建議使用 UMAP 唆涝，它可能比 tSNE 更忠實于流形的全局連通性，因此能更好地保留軌跡唇辨。
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)
sc.tl.umap(adata)
# 如果設(shè)置了 adata 的 .raw 屬性時廊酣，下圖顯示了“raw”（標準化、對數(shù)化但未校正）基因表達矩陣赏枚。
sc.pl.umap(adata, color=['CST3', 'NKG7', 'PPBP'])
為了繪制縮放矯正的基因表達聚類圖亡驰，需要使用 use_raw=False 參數(shù)。
sc.pl.umap(adata, color=['CST3', 'NKG7', 'PPBP'], use_raw=False)
目前還沒有計算出各個細胞類群饿幅，下面進行聚類
Leiden 圖聚類
# 計算
sc.tl.leiden(adata)
# 繪制
sc.pl.umap(adata, color=['leiden'])
6凡辱、檢索標記基因
先計算每個 leiden 分群中高度差異基因的排名，取排名前 25 的基因栗恩。
默認情況下透乾，使用 AnnData 的 .raw 屬性。
T-test
最簡單和最快的方法是 t 檢驗磕秤。
sc.tl.rank_genes_groups(adata, 'leiden', method='t-test')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)
Wilcoxon rank-sum
Wilcoxon rank-sum (Mann-Whitney-U) 檢驗的結(jié)果非常相似乳乌，還可以使用其他的差異分析包，如 MAST亲澡、limma钦扭、DESeq2 和 diffxpy。
sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)
# 保存這次的數(shù)據(jù)結(jié)果
adata.write(results_file)
邏輯回歸
sc.tl.rank_genes_groups(adata, 'leiden', method='logreg')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)
使用邏輯回歸對基因進行排名 Natranos et al. (2018)床绪，這里使用多變量方法客情，而傳統(tǒng)的差異測試是單變量 Clark et al. (2014)
除了僅由 t 檢驗發(fā)現(xiàn)的 IL7R 和由其他兩種方法發(fā)現(xiàn)的 FCER1A 之外，所有標記基因都在所有方法中都得到了重現(xiàn)癞己。
| Louvain Group | Markers | Cell Type |
|---|---|---|
| 0 | IL7R | CD4 T cells |
| 1 | CD14, LYZ | CD14+ Monocytes |
| 2 | MS4A1 | B cells |
| 3 | CD8A | CD8 T cells |
| 4 | GNLY, NKG7 | NK cells |
| 5 | FCGR3A, MS4A7 | FCGR3A+ Monocytes |
| 6 | FCER1A, CST3 | Dendritic Cells |
| 7 | PPBP | Megakaryocytes |
根據(jù)已知的標記基因膀斋，定義一個標記基因列表供以后參考：
marker_genes = ['IL7R', 'CD79A', 'MS4A1', 'CD8A', 'CD8B', 'LYZ', 'CD14',
'LGALS3', 'S100A8', 'GNLY', 'NKG7', 'KLRB1',
'FCGR3A', 'MS4A7', 'FCER1A', 'CST3', 'PPBP']
載入數(shù)據(jù)
# 使用 Wilcoxon Rank-Sum 測試結(jié)果重新加載已保存的對象
adata = sc.read(results_file)
獲取聚類分組和分數(shù)
result = adata.uns['rank_genes_groups']
groups = result['names'].dtype.names
pd.DataFrame(
{group + '_' + key[:1]: result[key][group]
for group in groups for key in ['names', 'pvals']}).head(5)
Group 1 與 Group 2 比較，進行差異分析
sc.tl.rank_genes_groups(adata, 'leiden', groups=['0'], reference='1', method='wilcoxon')
sc.pl.rank_genes_groups(adata, groups=['0'], n_genes=20)
sc.pl.rank_genes_groups_violin(adata, groups='0', n_genes=8)
Group 0 與其余組的比較進行差異分析
adata = sc.read(results_file)
sc.pl.rank_genes_groups_violin(adata, groups='0', n_genes=8)
跨類群比較基因
sc.pl.violin(adata, ['CST3', 'NKG7', 'PPBP'], groupby='leiden')
根據(jù)已知的細胞標記痹雅，注釋細胞類型
new_cluster_names = [
'CD4 T', 'CD14 Monocytes',
'B', 'CD8 T',
'NK', 'FCGR3A Monocytes',
'Dendritic', 'Megakaryocytes']
adata.rename_categories('leiden', new_cluster_names)
sc.pl.umap(adata, color='leiden', legend_loc='on data', title='', frameon=False, save='.pdf')
可視化每個類群的標記基因
氣泡圖顯示：
sc.pl.dotplot(adata, marker_genes, groupby='leiden');
小提琴圖顯示
sc.pl.stacked_violin(adata, marker_genes, groupby='leiden', rotation=90);
7却嗡、保存數(shù)據(jù)
保存壓縮文件
如果只想將其用于可視化的人共享此文件假残，減少文件大小的一種簡單方法是刪除縮放和校正的數(shù)據(jù)矩陣姥卢。
adata.write(results_file, compression='gzip')
保存為 h5ad 數(shù)據(jù)
adata.raw.to_adata().write('./write/pbmc3k_withoutX.h5ad')
讀取使用 adata = sc.read_h5ad('./write/pbmc3k_withoutX.h5ad')
導(dǎo)出數(shù)據(jù)子集
# 導(dǎo)出聚類數(shù)據(jù)
adata.obs[['n_counts', 'louvain_groups']].to_csv('./write/pbmc3k_corrected_louvain_groups.csv')
# 導(dǎo)出PCA數(shù)據(jù)
adata.obsm.to_df()[['X_pca1', 'X_pca2']].to_csv('./write/pbmc3k_corrected_X_pca.csv')
8、番外
我之前在處理較多數(shù)據(jù)量的時候赂苗，會有些地方不一樣，具體每個數(shù)據(jù)集的處理也會有比較大的自由度贮尉，比如：
在檢測線粒體基因時拌滋，這里在質(zhì)控時，已經(jīng)把線粒體基因直接剔除猜谚。
在做 UMAP 時败砂，可以看到一些類群間的聯(lián)系和軌跡。
做 TSNE時魏铅，可以看到類群間比較干凈利索昌犹，整體比較“飽滿”。
