本文首发于“生信补给站”公众号 https://mp.weixin.qq.com/s/dqxGfDRS2jg8jV_y1BBonA
经过前面单细胞免疫组库VDJ|从数据下载开始完成cellranger vdj分析（1）中的cellranger count 和 cellranger vdj的分析后，得到了单细胞转录组 和 单细胞TCR的结果数据 。考虑到样本数太多，选择既有RNA又有TCR样本的数据进行后续分析。重点解决以下三个小问题
1，只看TCR时（二），不同样本，不同分组情况下的clone差异以及变化；
2，TCR结合转录组数据（三+四），展示clone的分布可视化 以及 不同celltype中clone的差异以及变化；
3，一些可能的报错：（1）结合单细胞数据时；（2）官网链接中一些group相关的函数，按照推文或者使用??查看报错函数的具体说明
一 准备R包，数据集
安装scRepertoire包，如果出现如下的Error，按照提示安装gsl后再安装scRepertoire。
BiocManager::install("scRepertoire") suppressMessages(library(scRepertoire)) #Error: package or namespace load failed for ‘scRepertoire’ in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]): # 不存在叫‘gsl’这个名字的程辑包 #install.packages("gsl") suppressMessages(library(Seurat)) library(ggsci)
二 scRepertoire VDJ分析
1，处理VDJ数据
使用combineTCR函数结合所有样本的filtered_contig_annotations.csv数据，使用samples 和 ID 函数对样本进行标注和分组
S1 <- read.csv("./scTCR_data/t2_Normal.filtered_contig_annotations.csv") S2 <- read.csv("./scTCR_data/t2_PBMC.filtered_contig_annotations.csv") S3 <- read.csv("./scTCR_data/t3_Normal.filtered_contig_annotations.csv") S4 <- read.csv("./scTCR_data/t3_PBMC.filtered_contig_annotations.csv") S5 <- read.csv("./scTCR_data/t3_Center.filtered_contig_annotations.csv") S6 <- read.csv("./scTCR_data/t4_Normal.filtered_contig_annotations.csv") S7 <- read.csv("./scTCR_data/t4_PBMC.filtered_contig_annotations.csv") S8 <- read.csv("./scTCR_data/t4_Center.filtered_contig_annotations.csv") S9 <- read.csv("./scTCR_data/UT1_Normal.filtered_contig_annotations.csv") contig_list <- list(S1, S2, S3, S4, S5, S6, S7,S8,S9) combined <- combineTCR(contig_list, samples = c("t2", "t2", "t3", "t3","t3","t4", "t4", "t4", "UT1"), ID = c("Normal", "PBMC", "Normal", "PBMC","Center", "Normal", "PBMC", "Center","Normal"), cells ="T-AB")
也可以使用addVariable函数添加其他信息，如tissue ，age ，分期等，这个很实用，后续可以使用这些分组进行 clone 的各类比较
example <- addVariable(combined, name = "Tissue", variables = c("Tumor", "Tumor", "Tumor", "Tumor","Tumor", "Tumor", "Tumor", "Tumor","Normal")) example[[1]][1:5,]
2，VCJ分析以及可视化
2.1 Quantify Clonotypes 探索克隆型
使用quantContig 探索每个样本的unique clone信息（分别为比例 以及 个数）
p1 <- quantContig(combined, cloneCall="gene+nt", scale = T)+ theme(axis.text.x = element_text(angle = 30,vjust = 0.85,hjust = 0.75))#X坐标加点斜体，出图好看 p2 <- quantContig(combined, cloneCall="gene+nt", scale = F)+ theme(axis.text.x = element_text(angle = 30,vjust = 0.85,hjust = 0.75)) p1 + p2
注：这里cloneCall函数有四个参数
(1) “gene” ：使用包含 TCR/Ig 的 VDJC 基因
(2) “nt”：使用 CDR3 区域的核苷酸序列
(3) “aa” ：使用 CDR3 区域的氨基酸序列
(4) “gene+nt” 使用包含 TCR/Ig + CDR3 区域的核苷酸序列的 VDJC 基因
柱形图具体的数值可以添加 exportTable = T函数导出
quantContig_output <- quantContig(combined, cloneCall="gene+nt", scale = T, exportTable = T) quantContig_output
2.2 Clonotype Abundance克隆型丰度
使用abundanceContig函数 查看各个样本的克隆型总数的折线图，也可以添加exportTable = T函数输出折线图的具体结果
abundanceContig(combined, cloneCall = "gene", scale = F)
2.3 Length of Clonotypes克隆型长度
使用lengtheContig函数查看 CDR3 序列的长度分布。chain 函数可以选择是如左图的全部展示，还是如右图的TRA ，TRB分别展示
p11 <- lengthContig(combined, cloneCall="aa", chain = "combined")#左图p22 <- lengthContig(combined, cloneCall="aa", chain = "single") #右图p11 + p22
2.4 Compare Clonotypes比较克隆型
使用compareClonotypes函数查看样本之间的克隆型的比例和动态变化。cloneCall可以选择“gene+nt”
compareClonotypes(combined, numbers = 10, samples = c("t3_PBMC","t3_Normal","t3_Center"), #samples = c("t2_Normal", "t2_PBMC","t3_PBMC","t3_Normal","t3_Center"), cloneCall="aa", graph = "alluvial")
samples可以选择任意样本，但是如果top number clonotype sequences之间没有共享的话，则没有连线。
2.5 Clonal Space Homeostasis克隆空间稳态
可以通过clonalHomeostasis函数查看，各特定比例的克隆型（cloneTypes：Rare，，，Hyperexpanded）所占据该样本的相对比例
clonalHomeostasis(combined, cloneCall = "gene",cloneTypes = c(Rare = 1e-04,Small = 0.001,Medium = 0.01,Large = 0.1,Hyperexpanded = 1))
2.6 Clonal Proportion 克隆比例
可以通过clonalProportion函数查看克隆型的比例，按克隆型的出现频率将其进行排名，1:10表示每个样本中的前10个克隆型。
clonalProportion(combined, cloneCall = "nt")
2.7 Overlap Analysis 重叠分析
使用clonalOverlap函数分析样本之间的相似性，使用 clonesizeDistribution函数对样本进行聚类。
clonalOverlap(combined, cloneCall = "gene+nt", method = "overlap") clonesizeDistribution(combined, cloneCall = "gene+nt", method="ward.D2")
2.8 Diversity Analysis多样性分析
使用clonalDiversity函数进行多样性分析，会通过以下4各指标计算：
1)Shannon, 2) inverse Simpson, 3) Chao1和4)based Coverage Estimator (ACE)。
前两者一般用于估计基线多样性，Chao/ACE 指数用于估计样本的丰富度。
clonalDiversity(combined, cloneCall = "gene+nt", n.boots = 100) clonalDiversity(combined, cloneCall = "gene+nt", group = "ID", n.boots = 100)
注意此处不要用group.by 而是要用group 。
