单细胞数据分析笔记4: 免疫组库分析（scRepertoire）
scRepertoire 旨在从 10x Genomics Cell Ranger 管道获取filter contig（要求barcode没有任何标签（前缀）输出，处理该数据以根据两条 TCR 或 Ig 链分配克隆型，并分析克隆型动态。
可分为：
1）仅克隆型分析功能，例如独特的克隆型或克隆空间定量；
2）使用 Seurat、SingleCellExperiment 或 Monocle 3 软件包与 mRNA 表达数据的交互
1. 
1.1 clonotypes.csv (克隆型相关信息)
TCR_BCR-1
| Column | Description |
|---|---|
| clonotype_id | 该样本中共有序列所属的克隆型 ID（相同clonotype_id在不同样本之间不代表同一克隆型） |
| frequency | 具有该克隆型的cell barcodes的数量（虽然受测序深度等方面的限制，可能会miss一些高频克隆型，但也侧面部分反映了这个BCR/TCR的丰度） |
| proportion | 具有该克隆型的cell barcodes的分数 |
| cdr3s_aa | 以分号分隔的链：序列对列表，其中链是 TRA、TRB、TRG、TRD、IGK、IGL 或 IGH，序列是该链的 CDR3 氨基酸序列 |
| cdr3s_nt | 以分号分隔的链：序列对列表，其中链是 TRA、TRB、TRG、TRD、IGK、IGL 或 IGH，序列是该链的 CDR3 核苷酸序列 |
1.2 filtered_contig_annotations.csv (过滤后的contig相关的文件)
TCR_BCR-2
| Column | Description |
|---|---|
| barcode | 该序列的cell barcode |
| is_cell | True 或 False 值，指示cell barcode是否是一个细胞 |
| contig_id | 此序列的唯一标识符 |
| high_confidence | True 或 False 值，指示该序列是否被称为高置信度（不太可能是嵌合序列或其他伪影） |
| length | 该序列的长度（以核苷酸为单位） |
| chain | 与此序列相关的链：TRA、TRB、IGK、IGL 或 IGH |
| v_gene | 得分最高的 V 片段 |
| d_gene | 得分最高的 D 片段 |
| j_gene | 得分最高的 J 片段 |
| c_gene | 得分最高的 C 片段 |
| full_length | True 或 False 值，指示该序列是否为全长 |
| productive | True 或 False 值，指示该序列是否为富有成效 |
| cdr3 | 预测的 CDR3 氨基酸序列 |
| cdr3_nt | 预测的 CDR3 核苷酸序列 |
| reads | 与该序列对齐的读取数 |
| umis | 与该序列对齐的不同 UMI 的数量 |
| raw_clonotype_id | 此cell barcode所分配到的克隆型 ID |
| raw_consensus_id | 此序列分配到的共有序列的 ID |
1.3 consensus_annotations.csv (一致性序列相关的文件)
TCR_BCR-3
| Column | Description |
|---|---|
| clonotype_id | 此共有序列所分配到的克隆型 ID |
| consensus_id | 该共识序列的 ID |
| v_start | 共有序列上 V 区起始位置的基于 0 的索引 |
| v_end | 共有序列上 V 区末端位置的基于 0 的索引 |
| v_end_ref | 参考上 V 基因结束位置的从 0 开始的索引 |
| j_start | 共有序列上 J 区起始位置的基于 0 的索引 |
| j_start_ref | 参考上 J 基因起始位置的基于 0 的索引 |
| j_end | 共有序列上 J 区末端位置的基于 0 的索引 |
| cdr3_start | 共有序列上 CDR3 区域起始位置的基于 0 的索引 |
| cdr3_end | 共有序列上 CDR3 区域结束位置的基于 0 的索引 |
2. 
由于 CellRanger 的输出是 TCRA 和 TCRB 链的量化，没有细胞barcode与clonotype直接对应的数据，因此不能直接整合到seurat中。scRepertoire分析的第一步是把上面的文件转换为barcode与clonotype对应的数据。这是使用执行的combineTCR（combineBCR）函数，其中输入的是 contig_list。还可以通过样品和 ID 信息重新标记条形码，以防止重复。
## 安装并加载所需的R包
# install.packages("immunarch")
# BiocManager::install("scRepertoire")
library(scRepertoire)
library(immunarch)
library(Seurat)
library(tibble)
library(tidyr)
library(RColorBrewer)
library(scales)
library(ggpubr)
library(ggplot2)
# 读入contig_annotations.csv文件
N11 <- read.csv('/data/N11_BCR_filtered_contig_annotations.csv')
N12 <- read.csv('/data/N12_BCR_filtered_contig_annotations.csv')
N5 <- read.csv('/data/N5_BCR_filtered_contig_annotations.csv')
N6 <- read.csv('/data/N6_BCR_filtered_contig_annotations.csv')
N7 <- read.csv('/data/N7_BCR_filtered_contig_annotations.csv')
N2 <- read.csv('/data/N2_BCR_filtered_contig_annotations.csv')
P17 <- read.csv('/data/P17_BCR_filtered_contig_annotations.csv')
P18 <- read.csv('/data/P18_BCR_filtered_contig_annotations.csv')
N9 <- read.csv('/data/N9_BCR_filtered_contig_annotations.csv')
P11 <- read.csv('/data/P11_BCR_filtered_contig_annotations.csv')
P13 <- read.csv('/data/P13_BCR_filtered_contig_annotations.csv')
P7 <- read.csv('/data/P7_BCR_filtered_contig_annotations.csv')
# scRepertoire需要将这些样本储存在list中，后续命名添加分组等等
BCR_list <- list(N11,N12,N5,N6,N7,N2,P17,P18,N9,P11,P13,P7)
# 如果barcode多，样本多的话，这一步执行还是需要花费一段时间的
data_bcr <- combineBCR(BCR_list,
ID = c("N11","N12","N5","N6","N7","N2","P17","P18","N9","P11","P13","P7"),
samples = c("HC","HC","HC","HC","HC","Mild","Mild","Mild","Severe","Severe","Severe"),
threshold = 0.85)
# 还有其他的分组信息添加，还可以继续使用addVariable
# data_bcr <- addVariable(data_bcr, name = "Age",
# variables = c("50", "40", "80", "60", "45", "77"))
# 也可在combineBCR()后利用subsetContig()删除特定的列表元素（子集化）
# subset <- subsetContig(data_bcr, name = "ID",
# variables = c("N11","N2"))
参数说明：
ꔷ cells：指定VDJ类型。"T-AB"代表α-β TCR；"T-GD"代表 γ-δ TCR
ꔷ call.related.clones: 使用核苷酸序列和V基因来调用相关克隆。 默认设置为 TRUE。 FALSE 将返回 CTstrict 或严格克隆型作为 V 基因 + 氨基酸序列
ꔷ threshold: 数字越高，用于聚类的序列相似度就越高
ꔷ removeNA: "TRUE"代表删除任何没有值的链；"FALSE"（ 默认设置）代表并合并具有NA 值的单元格
ꔷ removeMulti: "TRUE"代表删除任何具有 2 条以上免疫受体链的细胞barcode；"FALSE"（ 默认设置）代表并合并具有 > 2 个链的单元
ꔷ filterMulti: "TRUE"代表用多个链分离细胞barcode中前2个表达的链；"FALSE"（ 默认设置）代表允许选择最高表达的轻链和重链
3. 
3.1 量化Clonotypes（使用quantContig 探索每个样本的unique clone信息）
library(ggsci)
p1 <- quantContig(data_bcr,
cloneCall="strict",
scale = T, # 将数值转换为百分比
chain = "both") + # 柱形图具体的数值可以添加 exportTable = T函数导出
# 自行修改配色
# scale_fill_manual(values = c("#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd",
# "#8c564b", "#e377c2", "#7f7f7f", "#bcbd22", "#17becf",
# "#ff0000", "#0000ff")) +
scale_color_rickandmorty() +
theme(axis.text.x = element_text(angle = 60, hjust = 1))
# 只想可视化特异性的chain，那么在chain这里选择参数即可。
# TCR可以选择“TRA”, “TRB”, “TRD”, “TRG”， BCR则是“IGH” or “IGL”
p2 <- quantContig(data_bcr,
cloneCall="strict",
scale = T,
chain = "IGL") +
# 自行修改配色
# scale_fill_manual(values = c("#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd",
# "#8c564b", "#e377c2", "#7f7f7f", "#bcbd22", "#17becf",
# "#ff0000", "#0000ff")) +
scale_color_rickandmorty() +
theme(axis.text.x = element_text(angle = 60, hjust = 1))
TCR_BCR-4
ꔷ 这里cloneCall函数有四个参数：
1）"gene"：使用VDJC 基因
2）"nt"：使用 CDR3 区域的核苷酸序列
3）"aa"：使用 CDR3 区域的氨基酸序列
4）"strict" 使用包含VDJC 基因 + CDR3 区域的核苷酸序列
3.2 Clonotype丰度（按丰度检查克隆型的相对分布）
p3 <- abundanceContig(data_bcr, cloneCall = "strict", scale = F)
#也可以添加分组，直接用分组来展示
p4 <- abundanceContig(data_bcr, cloneCall = "strict", scale = F, group.by = 'sample')
TCR_BCR-5
3.3 Clonotype长度（查看CDR3序列的长度分布）
p5<- lengthContig(data_bcr,
cloneCall="aa", # 此处cloneCall只能是"nt"或"aa"
chain = "both")
TCR_BCR-6
3.4 比较Clonotype（查看样本之间的克隆型和动态变化）
p6 <- compareClonotypes(data_bcr,
numbers = 5, # 每组中要绘制图表的克隆型的最大数量
cloneCall="strict",
graph = "alluvial") + # "alluvial" or "area"
scale_colour_brewer(palette = "Greens")
# 看任意样本之间的比较
p7 <- compareClonotypes(data_bcr,
cloneCall="strict",
numbers = 3,
graph = "alluvial",
samples = c("HC_N11", "Severe_P7")) +
scale_colour_brewer(palette = "Greens")
TCR_BCR-7
3.5 TCR 或 BCR 基因的相对使用
p8 <- vizGenes(data_bcr,
gene = "V",
chain = "IGL",
plot = "bar", # "bar" or "heatmap"
order = "variance",
scale = TRUE) # "TRUE"代表按每个样本的基因数量缩放；"FALSE"代表原始数据
TCR_BCR-8
# 查看单个链中基因的使用情况。例如，对不同样本之间 IGL V 和 J 使用的差异感兴趣：
p9 <- vizGenes(data_bcr[c(1:5)],
gene = "V",
chain = "IGL",
y.axis = "J", # 沿 y 轴分隔计数的变量
plot = "heatmap",
scale = TRUE,
order = "gene") # "gene"代表按基因名称排序；"variance"代表按单独变量类别之间的方差排序
p10 <- vizGenes(data_bcr[c(6:8)],
gene = "V",
chain = "IGL",
y.axis = "J",
plot = "heatmap",
scale = TRUE,
order = "gene")
p11 <- vizGenes(data_bcr[c(9:12)],
gene = "V",
chain = "IGL",
y.axis = "J",
plot = "heatmap",
scale = TRUE,
order = "gene")
TCR_BCR-9
4. 
4.1 克隆空间稳态（查看特定比例的克隆所占据的样本全部的相对比例）
通过检查克隆空间，可以有效地查看特定比例的克隆所占据的相对空间。可以将整个免疫受体测序视为一个量杯。克隆空间稳态询问杯子中不同比例的克隆填充的百分比。比例分割点在函数中的cloneType变量下设置，并且可以在基线处进行调整
p12 <- clonalHomeostasis(data_bcr, cloneCall = "aa",
cloneTypes = c(Rare = 1e-04,
Small = 0.001,
Medium = 0.01,
Large = 0.1,
Hyperexpanded = 1))
#也可以按照分组展示，只需要添加group.by参数
p13 <- clonalHomeostasis(data_bcr, cloneCall = "aa",
cloneTypes = c(Rare = 1e-04,
Small = 0.001,
Medium = 0.01,
Large = 0.1,
Hyperexpanded = 1),
group.by = 'sample')
TCR_BCR-10
4.2 克隆比例
与克隆空间稳态一样，克隆比例的作用是将克隆放入单独的容器中。关键区别在于，该clonalProportion()函数不是查看克隆与总数的相对比例，而是按总数对克隆进行排名。分割代表按拷贝或出现频率对克隆型进行排名，这意味着 1:10 是每个样本中前 10 个克隆型。默认箱位于函数中的split变量下，可以调整。
p14 <- clonalProportion(data_bcr, cloneCall = "gene",
split = c(10, 100, 1000, 10000, 30000, 1e+05))
TCR_BCR-11
4.3 重叠分析
查看样本之间的相似性。目前可以执行四种方法：1）overlap，2）Morisita指数，3）Jaccard指数，4）raw。overlap着眼于较小样本中独特克隆型长度的克隆型重叠；Morisita指数更为复杂，衡量种群内个体分散程度的生态指标，纳入了种群规模。Jaccard相似度指数与overlap非常相似 - Jaccard指数的分母是两个比较的并集，而不是使用较小样本的长度，从而导致数字小得多。
p15 <- clonalOverlap(data_bcr,
cloneCall = "strict",
method = "morisita")
# 分析样本克隆大小分布的相似性
p16 <- clonesizeDistribution(data_bcr,
cloneCall = "strict",
method="ward.D2")
TCR_BCR-12
4.4 Diversity Analysis 多样性分析
多样性也可以通过样本或其他变量来衡量。使用五个指标计算多样性：1）Shannon，2）inverse Simpson， 3）Chao1，4）Abundance-based Coverage Estimator (ACE)，5）inverse Pielou’s 均匀度度量。前两者一般用于估计基线多样性，Chao/ACE指数用于估计样本的丰富度。该函数的新实现包括使用最少数量的独特克隆型对 100 个引导程序 (n.boots) 进行下采样，作为更稳健的多样性估计。
p17 <- clonalDiversity(data_bcr,
cloneCall = "aa",
n.boots = 1000,
x.axis = 'sample', # 这里的sample就是我们前面data_bcr的分组
group = 'ID') # ID就是每个样本
TCR_BCR-13
4.5 分散比较
计算相对克隆型，并生成比较两个样本的散点图
p18 <- scatterClonotype(data_bcr,
cloneCall ="gene",
x.axis = "HC_N12",
y.axis = "Severe_P11",
dot.size = "total",
graph = "proportion") # 绘制比例或原始克隆型计数图
TCR_BCR-14
5. 
5.1 数据预处理
head(data_bcr$HC_N11$barcode)
## [1] "HC_N11_AGTTAGCGCTCACACGC" "HC_N11_TATCTGTTACTAAACGC" "HC_N11_TCGTAGAATGCGCACTT" "HC_N11_GTACAGTATGACTGGGC" "HC_N11_CTACTATATGCTAGCCC" "HC_N11_CTACGGGATGGTGTGGC"
head(scRNA_harmony)
orig.ident nCount_RNA nFeature_RNA percent.mt percent.HB RNA_snn_res.0.8 seurat_clusters
## N11_ACAAGCTTACGGAGGTG N11 500 244 0.600000 0 4 4
## N11_ATAGGCTGCTGCACGCT N11 500 258 1.200000 0 4 4
## N11_ATTCAGGGCTATTCACG N11 500 242 1.200000 0 0 0
## N11_CTTTCAAATGCATGAAT N11 500 272 1.600000 0 0 0
## N11_GGATGTTTACACGATAC N11 500 261 3.600000 0 4 4
## N11_TCCTCCCGCTTTGCAGT N11 500 256 1.400000 0 0 0
## N11_TCCTCTTATGGATATAC N11 500 260 1.400000 0 4 4
## N11_AAACGAAGCTGAATTGA N11 501 259 1.397206 0 0 0
## N11_ATTCCATATGGTAGCCG N11 501 260 1.996008 0 4 4
## N11_CACTGAACGACTTGAAC N11 501 281 1.397206 0 0 0
# 因为barcode和seurat流程中的不一致，故需要修改（scRepertoire提供函数stripBarcode()可去除标签，这里因自己的优点不一样就自己修改了）
# 定义一个函数，按指定符号分割字符串，并合并前两个分割后的元素
split_and_merge <- function(string, delimiter) {
splitted <- strsplit(string, delimiter)[[1]]
return(paste(splitted[2:3], collapse = delimiter))
}
for (i in 1:length(data_bcr)) {
barcode_i <- data_bcr[[i]]$barcode
for (n in 1:length(barcode_i)) {
barcode_i[n] <- split_and_merge(barcode_i[n], "_")
}
data_bcr[[i]]$barcode <- barcode_i
}
for (i in 1:length(data_tcr)) {
barcode_i <- data_tcr[[i]]$barcode
for (n in 1:length(barcode_i)) {
barcode_i[n] <- split_and_merge(barcode_i[n], "_")
}
data_tcr[[i]]$barcode <- barcode_i
}
# 结合 TCR 和 BCR
list.receptors <- c(data_bcr, data_tcr)
seurat <- combineExpression(list.receptors,
scRNA_harmony,
cloneCall = "gene",
proportion = FALSE,
cloneTypes = c(Single = 1, Small = 5, Medium = 20, Large = 100, Hyperexpanded = 500))
DimPlot(seurat, group.by = "cloneType", raster=FALSE) +
scale_color_manual(values = colorblind_vector(5), na.value="grey") +
theme(plot.title = element_blank())
TCR_BCR-15
ꔷ 为了对频率进行分类，scRepertoire提供可变比例，如果 proportion 选择"TRUE" 允许相对比例；当 选择"FALSE" 时将使用绝对频率来定义克隆型组，cloneTypes充当放置标签的容器。默认情况下，cloneTypes设置为等于 cloneTypes = c(Rare = 1e-4、Small = 0.001、Medium = 0.01、Large = 0.1、Hyperexpanded = 1)
ꔷ 需要注意的是，如果有重复的cell barcode（如果细胞同时具有 Ig 和 TCR），则不会添加免疫受体信息
5.2 克隆覆盖（可视化克隆扩增细胞在降维图上的位置）
默认为"PCA"和freq.cutpoint以生成等高线图，可通过选择bins的数量或freq.cutpoint的数量来修改轮廓。clonalOverlay()可用于查看所有单元格或使用facet通过元数据变量进行分面。当进行分面时，将保持整体降维，而等高线图将根据分面变量进行调整。
clonalOverlay(seurat,
reduction = "UMAP",
freq.cutpoint = 30,
bins = 10,
facet = "Group") +
guides(color = "none")
TCR_BCR-16
5.2 克隆网络
与 clonalOverlay() 类似，可以使用 clonalNetwork() 来查看簇之间共享的克隆型的相互作用网络。此函数显示来自起始节点的克隆的相对比例，结束节点由箭头指示。
可以使用 3 种方法来过滤克隆：
filter.clones:
ꔷ 选择一个数字来分离包含前 X 个细胞的克隆 (filter.clones = 2000)
ꔷ 选择"min"以确保所有组都缩放到最小组的大小
filter.identity:
对于选择要可视化的身份，显示单个组的往返网络连接
filter.proportion:
删除组中克隆数量少于一定比例的克隆
# No Identity filter
clonalNetwork(seurat,
reduction = "UMAP",
identity = "seurat_clusters",
filter.clones = NULL,
filter.identity = NULL,
cloneCall = "aa")
# Examining Cluster 2 only
clonalNetwork(seurat,
reduction = "UMAP",
identity = "seurat_clusters",
filter.identity = "C2",
cloneCall = "aa")
TCR_BCR-17
5.3 突出显示克隆型
通过highlightClonotypes()查看特定序列的克隆型。为了突出显示克隆型，首先需要使用cloneCall我们将使用的序列类型，然后使用sequence特定序列本身。
seurat <- highlightClonotypes(seurat,
cloneCall= "aa",
sequence = c("CARDAGTITVVGPSGIDRW_CMQALQTPRTF",
"NA_CMQALQTPRTF"))
DimPlot(seurat, group.by = "highlight",raster=FALSE) +
theme(plot.title = element_blank())
TCR_BCR-18
5.4 占用的剧目
可以通过使用occupiedscRepertoire()函数，并选择x 轴来显示单细胞对象元数据中的簇或其他变量，按簇查看分配到特定频率范围的细胞计数。
ꔷ proportion 可用于查看相对级别分组
ꔷ label 将返回克隆型的绝对数量
occupiedscRepertoire(seurat, x.axis = "Group")
TCR_BCR-19
5.5 冲积克隆型
修改meta数据后，可以使用"alluvialClonotypes()"函数查看多个类别的克隆型。本质上我们能够使用这些图来检查分类变量的交换，因为此函数将生成一个图表，其中每个克隆型按所谓的分层排列（需要一些时间，具体取决于总细胞的大小）
alluvialClonotypes(seurat, cloneCall = "gene",
y.axes = c("orig.ident", "Group", "celltype"),
color = "orig.ident")
# 特定子集
alluvialClonotypes(seurat, cloneCall = "gene",
y.axes = c("orig.ident", "Group", "celltype"),
color = "TRAV12-2.TRAJ42.TRAC_TRBV20-1.TRBD2.TRBJ2-3.TRBC2") +
scale_fill_manual(values = c("grey", colorblind_vector(2)[2]))
TCR_BCR-20
5.6 圆环
与冲积图一样，可以使用 circlize R 包中的弦图来可视化簇的互连。
library(circlize)
library(scales)
circles <- getCirclize(seurat,
group.by = "seurat_clusters")
# 为每个簇分配颜色
grid.cols <- scales::hue_pal()(length(unique(seurat@active.ident)))
names(grid.cols) <- levels(seurat@active.ident)
# 绘图
circlize::chordDiagram(circles,
self.link = 1,
grid.col = grid.cols)
# 对单细胞对象进行子集化来探索B细胞
subset <- subset(seurat, celltype == "B_cell")
circles <- getCirclize(subset, group.by = "seurat_clusters")
grid.cols <- scales::hue_pal()(length(unique(subset@active.ident)))
names(grid.cols) <- levels(subset@active.ident)
chordDiagram(circles, self.link = 1,
grid.col = grid.cols, directional = 1,
direction.type = "arrows",
link.arr.type = "big.arrow")
TCR_BCR-21
