625_单细胞TCR_BCR分析¶

一句话概述： 单细胞TCR/BCR分析通过测量每个T/B细胞的抗原受体序列来追踪免疫细胞的克隆扩增和多样性，scRepertoire 2是2025年发布的最新R包，整合了Seurat/Bioconductor生态实现从VDJ数据到克隆动力学的全链路分析。

核心知识点速查表¶

概念	白话解释
TCR	T细胞受体——T细胞表面识别抗原的蛋白，由α链和β链组成
BCR	B细胞受体——B细胞表面识别抗原的蛋白（即免疫球蛋白）
CDR3	互补决定区3——TCR/BCR中变异最大的区域，决定抗原特异性
Clonotype	克隆型——拥有相同TCR/BCR序列的一组细胞
Clonal expansion	克隆扩增——识别抗原后某个克隆型的细胞大量增殖
V(D)J recombination	V(D)J重组——产生TCR/BCR多样性的基因重排机制
Repertoire	免疫组库——一个个体所有T/B细胞受体序列的集合
scRepertoire	单细胞免疫组库分析R包（2025年v2发布）
Diversity index	多样性指数——衡量免疫组库丰富程度的指标

一、安装¶

# === scRepertoire 2安装（2025年发布）===
if (!require("BiocManager")) install.packages("BiocManager")
BiocManager::install("scRepertoire")  # 从Bioconductor安装

# === 依赖包 ===
BiocManager::install("Seurat")       # 单细胞分析核心包
install.packages("ggplot2")          # 可视化
install.packages("viridis")          # 配色方案

library(scRepertoire)  # 加载scRepertoire
library(Seurat)        # 加载Seurat
library(ggplot2)       # 加载ggplot2

二、数据输入与处理¶

2.1 读取10x VDJ数据¶

# === 读取10x Genomics Cell Ranger VDJ输出 ===
# Cell Ranger输出路径中有filtered_contig_annotations.csv

# 单个样本读取
tcr_sample1 <- read.csv(
    "sample1/outs/filtered_contig_annotations.csv"  # 10x VDJ输出
)

# 多个样本批量读取
contig_list <- list(
    S1 = read.csv("sample1/outs/filtered_contig_annotations.csv"),
    S2 = read.csv("sample2/outs/filtered_contig_annotations.csv"),
    S3 = read.csv("sample3/outs/filtered_contig_annotations.csv"),
    S4 = read.csv("sample4/outs/filtered_contig_annotations.csv")
)

# === 使用combineTCR整合TCR数据 ===
combined_tcr <- combineTCR(
    contig_list,                   # 输入contig列表
    samples = c("S1", "S2", "S3", "S4"),  # 样本名
    removeNA = TRUE,               # 移除NA值
    removeMulti = TRUE             # 移除有多条链的细胞（多重比对）
)

# 查看整合后的数据结构
str(combined_tcr[[1]])  # 第一个样本的结构
head(combined_tcr[[1]])  # 前几行

2.2 读取BCR数据¶

# === 使用combineBCR整合BCR数据 ===
bcr_contigs <- list(
    S1 = read.csv("bcr_sample1/filtered_contig_annotations.csv"),
    S2 = read.csv("bcr_sample2/filtered_contig_annotations.csv")
)

combined_bcr <- combineBCR(
    bcr_contigs,                   # 输入BCR contig列表
    samples = c("S1", "S2"),       # 样本名
    removeNA = TRUE,               # 移除NA
    removeMulti = TRUE             # 移除多条链
)

# BCR特有：可以分析体细胞超突变（SHM）
# SHM反映B细胞亲和力成熟的程度

2.3 支持的输入格式¶

# scRepertoire 2支持多种输入：
# 1. 10x Genomics Cell Ranger VDJ — filtered_contig_annotations.csv
# 2. AIRR格式 — airr_rearrangement.tsv（标准化格式）
# 3. BD Rhapsody — Dominant_Contigs.csv
# 4. Omniscope（原TRUST4）— barcode_report.tsv
# 5. WAT3R — filtered_contig_annotations.csv
# 6. MIXCR — clonotypes.tsv
# 7. Dandelion（Python）— filtered_contig_dandelion.tsv

# AIRR格式示例（行业标准）：
# read.csv("airr_output.tsv", sep="\t") 然后用combineTCR

三、克隆型分析¶

3.1 克隆型频率可视化¶

# === 克隆型丰度分布 ===
clonalQuant(
    combined_tcr,                  # 输入TCR数据
    cloneCall = "aa",              # 按氨基酸序列定义克隆型
    chain = "both",                # 使用双链（alpha+beta）
    scale = TRUE,                  # 按样本大小缩放
    exportTable = FALSE            # 不导出表格
) + 
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    labs(title = "Clonotype Abundance")

# cloneCall选项：
# "gene"    — 按V-J基因组合定义克隆型
# "nt"      — 按核苷酸CDR3序列
# "aa"      — 按氨基酸CDR3序列（推荐）
# "strict"  — 按完整V-D-J基因+CDR3核苷酸

3.2 克隆型比例分布¶

# === 克隆型大小分布 ===
clonalAbundance(
    combined_tcr,
    cloneCall = "aa",              # 按氨基酸序列
    chain = "both",                # 双链
    scale = FALSE                  # 不缩放（显示原始数量）
)

# === 克隆型比例（Homeostasis）===
clonalHomeostasis(
    combined_tcr,
    cloneCall = "aa",
    chain = "both",
    cloneSize = c(Rare = 1e-4,          # 稀有克隆
                  Small = 0.001,         # 小克隆
                  Medium = 0.01,         # 中等克隆
                  Large = 0.1,           # 大克隆
                  Hyperexpanded = 1)     # 超扩增克隆
)

3.3 克隆型重叠分析¶

# === 样本间克隆型重叠 ===
clonalOverlap(
    combined_tcr,
    cloneCall = "aa",              # 按氨基酸序列
    chain = "both",                # 双链
    method = "morisita"            # 重叠指数方法
) +
    labs(title = "Clonotype Overlap Between Samples")

# method选项：
# "overlap"   — 简单重叠系数
# "morisita"  — Morisita指数（推荐，考虑丰度）
# "jaccard"   — Jaccard指数（只看有无）
# "cosine"    — 余弦相似度
# "raw"       — 原始共享克隆数

四、多样性分析¶

4.1 克隆多样性指数¶

# === 计算多样性指数 ===
clonalDiversity(
    combined_tcr,
    cloneCall = "aa",              # 按氨基酸序列
    chain = "both",                # 双链
    metrics = c("shannon",         # Shannon多样性
                "simpson",         # Simpson多样性
                "chao1",           # Chao1丰富度估计
                "inv.simpson"),    # 逆Simpson
    group.by = "sample",           # 按样本分组
    n.boots = 100                  # 100次自举法
)

# 解读：
# Shannon越高 → 克隆型越多样（没有单个克隆型占主导）
# Simpson越高 → 越均匀
# 感染/肿瘤后通常Shannon降低（某些克隆大量扩增）

4.2 Rarefaction曲线¶

# === 克隆型稀释曲线 ===
clonalRarefaction(
    combined_tcr,
    cloneCall = "aa",              # 克隆定义
    chain = "both",                # 双链
    plot.type = "curve",           # 曲线类型
    n.boots = 3                    # 自举次数
) +
    labs(title = "Clonotype Rarefaction Curve",
         x = "Number of Cells",
         y = "Number of Unique Clonotypes")

# 曲线趋平 → 测序深度足够
# 曲线仍在上升 → 需要更多细胞

五、与Seurat整合分析¶

5.1 将TCR信息添加到Seurat对象¶

# === 已有Seurat对象（来自GEX数据）===
# seurat_obj <- readRDS("seurat_processed.rds")

# === 将TCR信息合并到Seurat对象 ===
seurat_obj <- combineExpression(
    combined_tcr,                  # TCR数据
    seurat_obj,                    # Seurat对象
    cloneCall = "aa",              # 克隆型定义
    chain = "both",                # 双链
    proportion = TRUE,             # 计算比例
    cloneSize = c(Single = 1,           # 单个细胞的克隆
                  Small = 5,            # 2-5个细胞的克隆
                  Medium = 20,          # 6-20个细胞的克隆
                  Large = 100,          # 21-100个细胞
                  Hyperexpanded = 500)  # >100个细胞
)

# 查看添加的元数据
head(seurat_obj@meta.data[, c("CTaa", "cloneSize", "Frequency")])

5.2 在UMAP上可视化克隆信息¶

# === UMAP上着色显示克隆大小 ===
DimPlot(seurat_obj, 
        group.by = "cloneSize",       # 按克隆大小分组
        order = c("Hyperexpanded", "Large", "Medium", "Small", "Single"),
        pt.size = 0.5) +
    scale_color_manual(values = c(
        "Hyperexpanded" = "#FF0000",   # 红色
        "Large" = "#FF6600",           # 橙色
        "Medium" = "#FFCC00",          # 黄色
        "Small" = "#66CC00",           # 绿色
        "Single" = "#CCCCCC"           # 灰色
    )) +
    labs(title = "Clonal Expansion on UMAP")

ggsave("clonal_expansion_umap.pdf", width = 10, height = 8)

# === 按细胞类型看克隆扩增 ===
occupiedscRepertoire(
    seurat_obj,
    x.axis = "cell_type",             # x轴：细胞类型
    label = FALSE                      # 不显示标签
) +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

5.3 克隆型与基因表达关联¶

# === 找扩增克隆特有的差异基因 ===
# 比较大克隆vs小克隆的基因表达差异
Idents(seurat_obj) <- "cloneSize"  # 设置分组标识

expanded_markers <- FindMarkers(
    seurat_obj,
    ident.1 = c("Hyperexpanded", "Large"),  # 扩增克隆
    ident.2 = "Single",                      # 单细胞克隆
    min.pct = 0.1,                           # 至少10%细胞表达
    logfc.threshold = 0.25                   # log2FC阈值
)

# 查看top差异基因
head(expanded_markers[order(expanded_markers$p_val_adj), ], 20)

# === 可视化特定克隆型的基因表达 ===
# 找到最大的克隆型
top_clones <- names(sort(table(seurat_obj$CTaa), decreasing = TRUE))[1:5]
cat("Top 5 克隆型:\n")
print(top_clones)

六、高级分析¶

6.1 基因使用分析¶

# === V基因使用频率 ===
vizGenes(
    combined_tcr,
    x = "TRBV",                    # x轴：TRB的V基因
    y = NULL,                      # 不分组
    plot = "bar",                  # 柱状图
    chain = "TRB",                 # TRB链
    scale = TRUE                   # 缩放
) +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 6))

# === V-J基因配对使用 ===
vizGenes(
    combined_tcr,
    x = "TRBV",                    # V基因
    y = "TRBJ",                    # J基因
    plot = "heatmap",              # 热图
    chain = "TRB",                 # TRB链
    scale = TRUE                   # 缩放
)

6.2 CDR3序列特征分析¶

# === CDR3长度分布 ===
clonalLength(
    combined_tcr,
    cloneCall = "aa",              # 氨基酸序列
    chain = "TRB"                  # TRB链
) +
    labs(title = "CDR3 Length Distribution")

# CDR3长度解读：
# TCR CDR3通常12-16个氨基酸
# 过短或过长可能是异常重排
# 不同疾病可能偏好不同CDR3长度

七、常见报错与解决¶

报错信息	原因	解决方案
`No clonotypes detected`	VDJ数据质量差或路径错误	检查Cell Ranger输出是否正常
`Barcode mismatch`	GEX和VDJ的barcode不匹配	确认来自同一样本的GEM well
`Multiple chains per cell`	一个细胞有多条同类型链	用`removeMulti=TRUE`过滤
`Low VDJ recovery`	VDJ捕获效率低	检查文库质量和测序深度
`NA in clonotype`	部分细胞没有完整的V(D)J注释	用`removeNA=TRUE`过滤
`combineExpression barcode error`	Seurat和TCR的barcode格式不一致	统一barcode后缀（-1等）

八、面试高频题¶

Q1：什么是克隆型（clonotype），怎么定义？¶

答：克隆型是一组拥有相同TCR/BCR序列的细胞——它们都来自同一个前体细胞。定义克隆型有多种粒度：(1) 按CDR3氨基酸序列——最常用，因为功能相同的TCR可能有不同的核苷酸序列；(2) 按CDR3核苷酸序列——更严格，能区分趋同进化；(3) 按完整V-D-J基因+CDR3——最严格；(4) 按V-J基因组合——最宽松。一般推荐用CDR3氨基酸序列（双链）定义，兼顾灵敏度和特异性。

Q2：克隆扩增在免疫应答中意味着什么？¶

答：克隆扩增意味着某些T/B细胞识别了特定抗原后大量增殖。在肿瘤中，大量扩增的CD8+ T细胞克隆可能是肿瘤特异性T细胞（能杀肿瘤的"好兵"）；在自身免疫病中，扩增的克隆可能是攻击自身组织的"叛军"；在感染中，扩增的克隆是对抗病原体的主力军。通过分析克隆扩增，我们可以：(1) 找到肿瘤新抗原特异性T细胞；(2) 评估免疫治疗的效果（治疗后是否有新克隆扩增）；(3) 追踪感染过程中免疫应答的动态变化。

Q3：单细胞TCR/BCR分析相比传统Bulk免疫组库有什么优势？¶

答： Bulk免疫组库（Bulk TCR-seq）只给出TCR序列列表，不知道每条TCR属于什么类型的细胞。单细胞TCR/BCR分析同时测量每个细胞的TCR/BCR序列和基因表达（通过GEX），可以：(1) 知道扩增的克隆是什么细胞类型（CD4+? CD8+? Treg? Teff?）；(2) 分析克隆扩增与细胞功能状态的关系（是否耗竭？是否活化？）；(3) 追踪同一克隆在不同细胞状态间的转换。缺点是成本高、细胞通量相对低（通常几千到几万个细胞）。

Q4：scRepertoire分析的基本流程是什么？¶

答：四步流程：(1) 用combineTCR/combineBCR读取10x Cell Ranger VDJ输出，整合多样本数据；(2) 克隆型分析——统计克隆型频率分布、大小分布、样本间重叠度；(3) 多样性分析——计算Shannon、Simpson等多样性指数，画稀释曲线；(4) 与Seurat整合——用combineExpression将TCR信息添加到GEX的Seurat对象中，在UMAP上可视化克隆扩增，分析扩增克隆的基因表达特征。关键是第4步的整合——把"这个T细胞的受体是什么"和"这个T细胞表达了什么基因"关联起来。

参考资料：scRepertoire 2: Borcherding et al., F1000Research 2025 | 10x Genomics VDJ: 10xgenomics.com | AIRR标准: Rubelt et al., Nat Immunol 2017 | Seurat: Hao et al., Cell 2024