20260910 113051 单细胞转录组细胞类型鉴定与轨迹分析全流程实现详解

import scanpy as sc

# 读取 10x 原始矩阵（未过滤），保留空液滴用于环境RNA估计
adata = sc.read_10x_h5('raw_feature_bc_matrix.h5')
adata.var_names_make_unique()

library(SoupX)
sc_ <- load10X('cellranger_outs/')   # 需要 raw + filtered 两个矩阵
sc_ <- autoEstCont(sc_)               # 估计污染比例 rho
counts_adj <- adjustCounts(sc_, roundToInt = TRUE)   # 输出非整数，取整给NB模型

import scanpy as sc
import numpy as np
from scipy.stats import median_abs_deviation

adata.var['mt'] = adata.var_names.str.startswith('MT-')    # 线粒体基因
adata.var['ribo'] = adata.var_names.str.startswith(('RPS', 'RPL'))
adata.var['hb'] = adata.var_names.str.contains(r'^HB[ABDEGMQZ]\d*(?!\w)')   # 血红蛋白亚基

sc.pp.calculate_qc_metrics(adata, qc_vars=['mt', 'ribo', 'hb'],
                           percent_top=[20], log1p=True, inplace=True)

def is_outlier(adata, metric, nmads):
    M = adata.obs[metric]
    return (M < np.median(M) - nmads * median_abs_deviation(M)) | \
           (np.median(M) + nmads * median_abs_deviation(M) < M)

adata.obs['outlier'] = (is_outlier(adata, 'log1p_total_counts', 5) |
                        is_outlier(adata, 'log1p_n_genes_by_counts', 5) |
                        is_outlier(adata, 'pct_counts_in_top_20_genes', 5))
adata.obs['mt_outlier'] = is_outlier(adata, 'pct_counts_mt', 3) | (adata.obs['pct_counts_mt'] > 8)
adata = adata[~(adata.obs['outlier'] | adata.obs['mt_outlier'])].copy()
sc.pp.filter_genes(adata, min_cells=3)   # 去掉在<3个细胞中出现的基因

# 先保存原始计数（HVG 和 doublet 检测需要）
adata.layers['counts'] = adata.X.copy()

# shifted-log 归一化（默认），sc.pp.normalize_total 将每个细胞缩放到中位数
sc.pp.normalize_total(adata)          # target_sum=None 用中位数
sc.pp.log1p(adata)                     # 自然对数(1+x)

# seurat_v3 flavor 需要原始计数，必须在 layer 中
sc.pp.highly_variable_genes(adata, n_top_genes=2000, flavor='seurat_v3', layer='counts')

import scanpy as sc

# Harmony 整合
sc.external.pp.harmony_integrate(adata, key='sample')   # 以样本为批次键
adata.obsm['X_pca_harmony'] = adata.obsm['X_pca_harmony']

sc.tl.pca(adata, n_comps=50, svd_solver='arpack')
sc.pl.pca_variance_ratio(adata, n_pcs=50, log=True)   # 肘图选择 n_pcs

sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30)
# 固定后端，保证可重复
sc.tl.leiden(adata, resolution=0.5, flavor='igraph', n_iterations=2, directed=False)
sc.tl.umap(adata, min_dist=0.3)
sc.pl.umap(adata, color=['leiden', 'CD3D', 'MS4A1', 'CD14'])

# 扫描分辨率，观察细胞流动
for res in [0.2, 0.4, 0.6, 0.8, 1.0]:
    sc.tl.leiden(adata, resolution=res, key_added=f'leiden_r{res}',
                 flavor='igraph', n_iterations=2, directed=False)
    print(res, adata.obs[f'leiden_r{res}'].nunique(), 'clusters')

# R中用clustree看跨分辨率细胞流动
# clustree(seurat_obj, prefix = 'RNA_snn_res.')

import scanpy as sc

# 注意：scanpy 默认是 t-test，必须显式用 wilcoxon
sc.tl.rank_genes_groups(adata, groupby='leiden', method='wilcoxon',
                        pts=True, corr_method='benjamini-hochberg')
markers = sc.get.rank_genes_groups_df(adata, group=None)

# 用效应量和表达比例过滤，不用 p 值
specific = markers[(markers['logfoldchanges'] > 1) &
                   (markers['pct_nz_group'] > 0.5) &
                   (markers['pct_nz_reference'] < 0.25)]
print(specific.groupby('group').head(10)[['group', 'names', 'logfoldchanges', 'pct_nz_group', 'pct_nz_reference']])

t_cell_panel = ['CD3D', 'CD3E', 'CD4', 'CD8A', 'CD8B']
sc.tl.score_genes(adata, gene_list=t_cell_panel,
                  ctrl_size=50, n_bins=25, score_name='T_cell_score')

cluster_labels = {'0': 'CD4 T', '1': 'CD14 Mono', '2': 'B', '3': 'CD8 T',
                  '4': 'NK', '5': 'FCGR3A Mono'}
adata.obs['cell_type'] = adata.obs['leiden'].map(cluster_labels).fillna('Unassigned')

# 按 sample x cell_type 聚合原始计数（不是归一化值）
cell_type = adata[adata.obs['cell_type'] == 'CD14 Mono']
pseudobulk = sc.get.aggregate(cell_type, by='sample', func='sum')
counts_df = pseudobulk.layers['sum']
# 之后喂给 DESeq2/edgeR/limma-voom，见 deseq2-basics

import scanpy as sc
sc.pp.neighbors(adata, n_neighbors=15, use_rep='X_pca')
sc.tl.leiden(adata, resolution=1.0, flavor='igraph', n_iterations=2, directed=False)
sc.tl.paga(adata, groups='leiden')
sc.pl.paga(adata, threshold=0.03, color='leiden')   # 修剪低连通边
sc.tl.umap(adata, init_pos='paga')                  # 全局拓扑保留

sc.tl.diffmap(adata, n_comps=15)
# 根细胞用已知 marker 锚定，不在 UMAP 上用眼睛选
adata.uns['iroot'] = np.flatnonzero(adata.obs['cell_type'] == 'HSC')[0]
sc.tl.dpt(adata, n_dcs=10, n_branchings=0)   # n_branchings=0 -> 纯伪时间

import palantir
dm_res = palantir.utils.run_diffusion_maps(adata, n_components=5)
ms_data = palantir.utils.determine_multiscale_space(dm_res)
pr_res = palantir.core.run_palantir(ms_data, early_cell='HSC_cell_id',
                                    terminal_states=None, num_waypoints=1200)
# pr_res.pseudotime, pr_res.entropy, pr_res.branch_probs

import cellrank as cr
pk = cr.kernels.PseudotimeKernel(adata, time_key='dpt_pseudotime').compute_transition_matrix()
ck = cr.kernels.ConnectivityKernel(adata).compute_transition_matrix()
combined = 0.8 * pk + 0.2 * ck     # 权重是研究者选择，需扫描

g = cr.estimators.GPCCA(combined)
g.compute_macrostates(n_states=10, cluster_key='leiden')
g.predict_terminal_states(method='stability')
g.predict_initial_states(n_states=1)
g.compute_fate_probabilities()
g.compute_lineage_drivers()

import scvelo as scv
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
scv.tl.recover_dynamics(adata)
scv.tl.velocity(adata, mode='dynamical')   # 默认是 'stochastic'，显式传 'dynamical'
scv.tl.velocity_graph(adata)
scv.tl.velocity_confidence(adata)   # 先检查置信度

import liana as li
import scanpy as sc

adata = sc.read_h5ad('adata_annotated.h5ad')
li.mt.rank_aggregate(adata, groupby='cell_type', resource_name='consensus',
                     expr_prop=0.1, use_raw=False, n_perms=1000, verbose=True)

res = adata.uns['liana_res']
# 输出 magnitude_rank 和 specificity_rank（不是单一 liana_rank）
robust = res[(res['specificity_rank'] < 0.05) & (res['magnitude_rank'] < 0.05)]

from liana.method import cellphonedb
for resource in ['consensus', 'cellphonedb', 'cellchatdb']:
    cellphonedb(adata, groupby='cell_type', resource_name=resource,
                expr_prop=0.1, use_raw=False, key_added=f'cpdb_{resource}', verbose=False)
# 比较 adata.uns['cpdb_*'] 的 top 对——跨资源存活的才可靠

library(CellChat)
cellchat <- createCellChat(object = seurat_obj, group.by = 'cell_type')
cellchat@DB <- CellChatDB.human
cellchat <- subsetData(cellchat)
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)
cellchat <- computeCommunProb(cellchat, type = 'triMean')   # 保守 trim 均值
cellchat <- filterCommunication(cellchat, min.cells = 10)
cellchat <- computeCommunProbPathway(cellchat)
cellchat <- aggregateNet(cellchat)
cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = 'netP')
# netVisual_aggregate(signaling='WNT'), netVisual_bubble(), netAnalysis_signalingRole_heatmap()

library(nichenetr)
geneset_oi <- FindMarkers(seurat_obj, ident.1 = 'activated_T', ident.2 = 'naive_T') %>%
    filter(p_val_adj < 0.05, avg_log2FC > 0.5) %>% rownames()
background <- get_expressed_genes('T_cell', seurat_obj, pct = 0.10)

expressed_ligands <- intersect(unique(lr_network$from),
                               get_expressed_genes(c('Macrophage', 'Dendritic'), seurat_obj, 0.10))
expressed_receptors <- intersect(unique(lr_network$to), background)
potential_ligands <- lr_network %>%
    filter(from %in% expressed_ligands, to %in% expressed_receptors) %>%
    pull(from) %>% unique()

ligand_activities <- predict_ligand_activities(
    geneset = geneset_target, background_expressed_genes = background,
    ligand_target_matrix = ligand_target_matrix, potential_ligands = potential_ligands)

原始 10x 矩阵
  │
  ├─ Step 1: 数据读取 + 环境RNA去除（SoupX）+ 质控（MAD自适应, 逐样本）
  │
  ├─ Step 2: 归一化（shifted-log）+ 高变基因选择（seurat_v3需raw counts）
  │
  ├─ Step 3: 批次整合（Harmony/scVI，整合后不用于DE）
  │
  ├─ Step 4: 聚类（PCA elbow → kNN图 → Leiden，固定后端，分辨率扫描）
  │
  ├─ Step 5: Marker检测（Wilcoxon, 双态校正感知）+ 手动注释（marker面板）
  │
  ├─ Step 6: 轨迹推断（PAGA拓扑 → Slingshot/DPT/Palantir/CellRank2，伪时间）
  │
  └─ Step 7: 细胞通讯（LIANA共识 + CellChat通路 + NicheNet活性）