在完成質控(QC)后,我們已經過濾掉了低質量細胞、雙細胞和低表達基因,獲得了較為干凈的單細胞數據集單細胞分析教程 | (一)Python單細胞質控全流程。接下來,我們將進行以下關鍵步驟:
1. 數據標準化(Normalization):消除測序深度和細胞大小的影響。
2. 特征選擇(Feature Selection):挑選高變基因以減少噪聲。
3. 降維(Dimensionality Reduction):將高維數據投影到低維空間,便于可視化和分析。
4. 聚類(Clustering):將細胞分組以識別細胞類型。
5. 可視化(Visualization ):展示聚類結果。
通過這些步驟,我們可以從原始計數數據中提取生物學意義,為后續差異表達分析和功能注釋打下基礎。以下是詳細流程和代碼實現。
數據標準化(Normalization)
在進行后續內容之前,同樣需要導包等操作,如下,不再進行詳細介紹:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import scanpy as sc
import randomseed = 42
random.seed(seed)
np.random.seed(seed)adata = sc.read_h5ad("./qc_processed.h5ad")
單細胞數據的原始計數受到測序深度、細胞大小等技術因素的影響,因此需要標準化以確保細胞間表達值的可比性。常用的標準化方法是對計數進行歸一化(通常按細胞總計數歸一化到固定值,如10,000),然后進行對數化以穩定方差。
# 標準化:按總計數歸一化到10,000
sc.pp.normalize_total(adata, target_sum=1e4)# 對數化:log1p變換
sc.pp.log1p(adata)# 保存標準化后的數據到adata.raw,用于后續差異表達分析
adata.raw = adataadata
注意:
-
normalize_total?將每個細胞的計數歸一化為總和為10,000(可根據需要調整?target_sum)。
-
log1p?使用自然對數(log(x+1))變換數據,減少高表達基因的方差影響。
-
保存原始標準化數據到?adata.raw,以便后續分析(如差異表達基因計算)使用未進一步處理的數據。
在很多的代碼中可能會看到sc.pp.scale,這一步操作的意義是為了改變數據的分布模式使每個基因的表達量均值為0,方差為1。做不做這一步主要取決于normalization以及log1p之后是否還存在超高表達的基因。
此外,需要注意的時候,后面的差異表達分析等操作,最好是使用normalization以及log1p之后的結果,所以這里才需要將數據保存一下,并且scale操作可以放在選擇高變基因后再進行。
特征選擇(Feature Selection)
單細胞數據通常包含數千到數萬個基因,但并非所有基因都與生物學差異相關。選擇高變基因(Highly Variable Genes, HVG)可以減少噪聲,聚焦于生物學上重要的基因。
# 識別高變基因
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)# 可視化高變基因
sc.pl.highly_variable_genes(adata)# 過濾數據集,僅保留高變基因
adata = adata[:, adata.var.highly_variable].copy()
注意:
-
min_mean?和?max_mean?控制基因的平均表達量范圍,min_disp?控制基因的分散度(variance/mean),這些參數需要根據數據分布調整。
-
高變基因通常占總基因數的10-20%(約2000-4000個基因),具體數量取決于數據集和實驗設計。
-
可視化高變基因的散點圖可以幫助確認參數設置是否合理,理想情況下高變基因應集中在高分散度區域。
降維
高維基因表達數據難以直接分析和可視化,因此需要降維。常用的方法包括主成分分析(PCA)和非線性降維(如t-SNE或UMAP)。我們先通過PCA提取主要特征,再使用UMAP進行可視化。
# 縮放數據(零均值、單位方差)
sc.pp.scale(adata, max_value=10)# 運行PCA
sc.tl.pca(adata, svd_solver='arpack')# 可視化PCA結果
sc.pl.pca(adata, color='n_counts', show=False)
sc.pl.pca_variance_ratio(adata, log=True)# 運行UMAP
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)
sc.tl.umap(adata)# 可視化UMAP
sc.pl.umap(adata, color=['n_counts', 'percent_mt'], show=False)
![]() | ![]() |
注意:
-
sc.pp.scale?確保每個基因的表達量標準化為均值0、方差1,避免高表達基因主導PCA結果。
-
max_value=10?限制縮放后的最大值,防止極端值影響結果。
-
PCA的?n_pcs(主成分數量)通常設置為30-50,具體可通過?pca_variance_ratio?圖檢查解釋方差的比例。
-
UMAP的?n_neighbors?控制局部結構的保留程度,n_pcs?指定使用的PCA主成分數,可根據數據集大小調整。
聚類(Clustering)
通過聚類,我們可以將相似的細胞分組,初步識別可能的細胞類型。常用的方法是基于圖的聚類算法(如Louvain或Leiden)。
# 運行Leiden聚類
sc.tl.leiden(adata, resolution=0.8)# 可視化聚類結果
sc.pl.umap(adata, color='leiden', legend_loc='on data', title='Leiden Clustering')
注意:
-
resolution?參數控制聚類的細粒度,值越大,得到的簇越多(通常在0.5-1.5之間調整)。
-
Leiden算法相比Louvain更穩定,推薦使用。
-
可視化時,legend_loc='on data'?將簇編號直接標注在圖上,便于觀察。
可視化與初步注釋
為了進一步理解聚類結果,我們可以通過已知標記基因(marker genes)對簇進行初步注釋。以下以人類數據為例,假設我們關注幾種常見細胞類型(如T細胞、巨噬細胞、B細胞等)。當然前提是你已經知道了一些marker。這樣可以展示出高表達這些maker的細胞的大致分布。
# 定義標記基因(根據研究背景調整)
marker_genes = {'T_cells': ['CD3D', 'CD8A'],'Macrophages': ['CD68', 'MARCO'],'B_cells': ['CD19', 'MS4A1']
}# 可視化標記基因表達
sc.pl.umap(adata, color=marker_genes['T_cells'], title='T Cell Markers')
sc.pl.umap(adata, color=marker_genes['Macrophages'], title='Macrophage Markers')
sc.pl.umap(adata, color=marker_genes['B_cells'], title='B Cell Markers')# 點圖展示標記基因表達
sc.pl.dotplot(adata, marker_genes, groupby='leiden', dendrogram=True)
注意:
-
標記基因需要根據具體組織或實驗背景選擇,建議查閱相關文獻或數據庫(如PanglaoDB)。
-
dotplot?顯示每個簇中標記基因的平均表達量和表達細胞比例,適合初步判斷簇的細胞類型。
-
如果某些標記基因不在數據中,檢查基因名是否正確(大小寫敏感)或是否被QC過濾。
可視化(附加選擇):聚類統計
為了更好地理解聚類結果,可以統計每個簇的細胞數量并繪制柱狀圖。
import seaborn as sns
import matplotlib.pyplot as plt# 統計每個簇的細胞數量
cluster_counts = adata.obs['leiden'].value_counts().sort_index()# 繪制柱狀圖
plt.figure(figsize=(10, 6))
sns.barplot(x=cluster_counts.index, y=cluster_counts.values, color='skyblue')
plt.xlabel('Cluster')
plt.ylabel('Number of Cells')
plt.title('Cell Counts per Cluster')
plt.show()
總結與下一步
保存結果:
# 保存處理后的AnnData對象
adata.write('processed_clustered.h5ad', compression='gzip')
至此,我們完成了單細胞RNA-seq分析的核心步驟:從質控到聚類和初步注釋。后續分析可以包括:
-
差異表達分析:識別每個簇的特征基因(使用?sc.tl.rank_genes_groups)。
-
細胞類型精細注釋:結合更多標記基因或自動注釋工具(如SingleR)。
-
軌跡分析:探索細胞發育或分化路徑(使用?sc.tl.paga?或?sc.tl.diffmap)。
-
整合多組數據:處理批次效應(使用?sc.pp.combat?或?sc.external.pp.harmony)。