385_单细胞TCR_BCR免疫组库分析¶

import scirpy as ir       # pip install scirpy
import scanpy as sc
import pandas as pd
import matplotlib.pyplot as plt

# ============================================================
# 步骤1：加载 scRNA-seq 和 VDJ 数据
# ============================================================
# 加载已预处理的 scRNA-seq 数据
adata = sc.read_h5ad("preprocessed_T_cells.h5ad")

# 加载 CellRanger VDJ 输出（filtered_contig_annotations.csv）
# CellRanger VDJ 输出目录结构：
# vdj_t/
#   outs/
#     filtered_contig_annotations.csv   <- 关键文件
#     all_contig_annotations.csv

mdata = ir.io.read_10x_vdj("vdj_t/outs/filtered_contig_annotations.csv")

# ============================================================
# 步骤2：将 VDJ 数据整合到 scRNA-seq AnnData
# ============================================================
ir.pp.merge_with_ir(
    adata,
    mdata   # 将免疫组库信息添加到 adata.obs 中
)

# 检查合并结果
print(adata.obs.columns.tolist())   # 应该包含 IR_VJ_1_v_call, IR_VDJ_1_cdr3 等列
print(f"有TCR数据的细胞: {(~adata.obs['IR_VDJ_1_cdr3'].isna()).sum()}")

# ============================================================
# 步骤3：定义克隆型
# ============================================================
# 按CDR3序列定义克隆型（相同α+β链CDR3 = 同一克隆型）
ir.tl.chain_qc(adata)   # 质控：检查每个细胞有几条链

ir.pp.ir_dist(
    adata,
    metric="hamming",     # Hamming距离（允许N个碱基错误匹配）
    sequence="aa",        # 按氨基酸序列比较
    cutoff=0              # cutoff=0表示完全匹配才算同一克隆型
)

ir.tl.define_clonotypes(
    adata,
    receptor_arms="VDJ",    # 只用β链定义克隆型（更严格时用"all"）
    dual_ir="any"           # 有双链时的处理策略
)

# 查看克隆型分布
clonotype_counts = adata.obs["clone_id"].value_counts()
print(f"总克隆型数: {clonotype_counts.nunique()}")
print(f"Top5 最大克隆型:\n{clonotype_counts.head()}")

# ============================================================
# 步骤4：克隆扩增分析
# ============================================================
# 计算克隆扩增（每个克隆型有多少细胞）
ir.tl.clonal_expansion(
    adata,
    target_col="clone_id"    # 克隆型ID列
)

# 在 UMAP 上展示克隆扩增
ir.pl.clonal_expansion(
    adata,
    groupby="cell_type",     # 按细胞类型分组展示
    show_nonexpanded=True,   # 也展示未扩增的克隆
    save="clonal_expansion.pdf"
)

# 克隆扩增 vs 细胞类型（堆叠柱状图）
ir.pl.group_abundance(
    adata,
    groupby="clone_id_size",  # 按克隆大小分组
    target_col="cell_type",    # 展示各细胞类型
    max_cols=5
)

# ============================================================
# 步骤5：克隆型多样性分析
# ============================================================
# 计算各样本的TCR多样性（Shannon熵）
ir.tl.alpha_diversity(
    adata,
    groupby="sample_id",      # 按样本计算多样性
    target_col="clone_id",    # 克隆型ID
    metric="normalized_shannon"  # 归一化Shannon熵（0-1）
)

ir.pl.alpha_diversity(
    adata,
    groupby="condition",    # 按疾病状态比较多样性
    target_col="clone_id"
)

library(immunarch)   # 免疫组库分析（R）

# 加载 CellRanger VDJ 输出（也支持多种格式：MiXCR、IMGT等）
immdata <- repLoad("vdj_t/outs/")    # 自动识别格式

# 查看数据结构
names(immdata)   # $data（克隆型列表）+ $meta（元数据）

# ============================================================
# 克隆型丰度统计
# ============================================================
# 查看每个样本中Top克隆型占比
top_clones <- repTop(immdata$data, .head = 10)
vis(top_clones)    # 可视化Top克隆型

# ============================================================
# 多样性分析
# ============================================================
div_results <- repDiversity(
  .data = immdata$data,
  .method = "chao1"    # 或"shannon"、"inv.simp"
)
vis(div_results,
  .by = "condition",   # 按条件分组（在meta数据中定义）
  .meta = immdata$meta
)

# ============================================================
# V基因使用频率
# ============================================================
# 分析 V、D、J 基因的使用偏好
vdj_usage <- geneUsage(immdata$data, "hs.trbv", .norm = TRUE)  # TRBV基因（人类）
vis(vdj_usage, .by = "condition", .meta = immdata$meta)

# ============================================================
# 样本间共享克隆型分析
# ============================================================
# 分析不同样本之间有多少共同克隆型（公共免疫应答）
shared <- repOverlap(immdata$data, .method = "overlap")  # Overlap系数
vis(shared, .text.size = 1.5)   # 热图展示样本间共享程度