首页站长学院单细胞免疫组库VDJ| 从零开始scRepertoire分析，解决真实场景中可能的问题

单细胞免疫组库VDJ| 从零开始scRepertoire分析，解决真实场景中可能的问题

站长网2023-07-27 10:55:413阅

选择既有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 contigs

values total scaled

795 t2_Normal 1221 65.11057

1672

t2_PBMC 1903 87.86127

759 t3_Normal 1044 72.70115

4462

t3_PBMC 6695 66.64675

2387 t3_Center 7615 31.34603

1159 t4_Normal 1575 73.58730

2661

t4_PBMC 3274 81.27673

2403 t4_Center 5395 44.54124

1618 UT1_Normal 2918 55.44894

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 。

三 scRepertoire-VDJ RNA

以上即为单细胞免疫组库的常规分析，我们更多的还是与scRNA-seq数据结合，就可以在umap上查看clone的分布情况以及基于cluster/celltype展示clonotype的分布情况。

3.1 Seurat标准流程

folders的样本文件夹内为cellranger的3个标准文件，的简单的走一遍标准流程单细胞工具箱|Seurat官网标准流程，此处仅为示例就不考虑批次问题了。

folders=list.files('./')folders#[1] "t2_Normal" "t2_PBMC"

"t3_Center" "t3_Normal" "t3_PBMC"

"t4_Center" "t4_Normal" "t4_PBMC"

"UT1_Normal"

#批量读取10X单细胞转录组数据文件夹sceList = lapply(folders,function(folder){ CreateSeuratObject(counts = Read10X(folder),

project = folder )})sce.all <- merge(sceList[[1]],

y = c(sceList[[2]],sceList[[3]],sceList[[4]],

sceList[[5]],sceList[[6]],sceList[[7]],sceList[[8]],sceList[[9]]),

add.cell.ids = folders, #添加样本名

project = "scRNA")sce.all #查看#线粒体比例sce.all[["percent.mt"]] <- PercentageFeatureSet(sce.all, pattern = "^MT-")#过滤sce <- subset(sce.all, subset = nFeature_RNA > 500 &

nFeature_RNA < 5000 &

percent.mt < 20)VlnPlot(sce, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)#标准化sce <- NormalizeData(sce)#高变基因sce <- FindVariableFeatures(sce, selection.method = "vst", nfeatures = 2000)all.genes <- rownames(sce)#归一化sce <- ScaleData(sce, features = all.genes)#降维聚类sce <- RunPCA(sce, npcs = 50)sce <- FindNeighbors(sce, dims = 1:30)sce <- FindClusters(sce, resolution = 0.5)sce <- RunUMAP(sce, dims = 1:30)sce <- RunTSNE(sce, dims = 1:30)p1 <- DimPlot(sce, reduction = "umap", pt.size=0.5, label = F,repel = TRUE)p2 <- DimPlot(sce, reduction = "umap",group.by = "orig.ident", pt.size=0.5, label = F,repel = TRUE)p1 p2

(1) 添加meta，处理barcode

读取文章提供的meta信息，如下图所示，红框中的不一致的部分需要tidyverse处理一下Tidyverse|数据列的分分合合，一分多，多合一

meta <- read.table("ccRCC_6pat_cell_annotations.txt",

sep = "\t",header = T)meta2 <- meta %>% separate(cell,into = c( "CB","sample","pos"),sep = "_") %>% mutate(cell = paste(sample,pos,CB,sep = "_"))

sce@meta.data <- sce@meta.data %>% rownames_to_column("cell") %>% inner_join(meta2,by = "cell") %>% column_to_rownames("cell")head(sce@meta.data)#

orig.ident nCount_RNA nFeature_RNA percent.mt RNA_snn_res.0.5 seurat_clusters

CB#t2_Normal_AAACCTGCAACACCTA-1 t2_Normal

7863

2244 3.700878

9 AAACCTGCAACACCTA-1#t2_Normal_AAACCTGCACCTCGTT-1 t2_Normal

2194

1188 10.027347

7 AAACCTGCACCTCGTT-1#

sample

pos type region Sample Sample2 cluster

cluster_name

UMAP1

UMAP2#t2_Normal_AAACCTGCAACACCTA-1

t2 Normal Normal Normal

t2 t2_Normal

cDC2 -9.836647 0.07936505#t2_Normal_AAACCTGCACCTCGTT-1

t2 Normal Normal Normal

t2 t2_Normal

11 CD45- Vascular Endothelium -7.601253 9.00771236#

Sample_name#t2_Normal_AAACCTGCAACACCTA-1 Ipi/Nivo resistant#t2_Normal_AAACCTGCACCTCGTT-1 Ipi/Nivo resistant

(2) sce的细胞数和 meta信息不一致

因为前面Seurat流程使用的自己的参数，而meta信息是文章给的。

解决：使用subset函数提取 meta.data中的cell

sce@meta.data <- sce@meta.data %>% rownames_to_column("Cell") %>% mutate(Cell2 = Cell) %>% column_to_rownames("Cell2")sce <- subset(sce,Cell %in% sce@meta.data$Cell)

3.2 VDJ结果结合RNA

使用combineExpression函数将VDJ结果和 RNA结果结合。

注：combined的barcode 和 sce的rownames 一定要一致，有出入的使用tidyverse包进行数据处理。

seurat <- combineExpression(combined, sce,

cloneCall="gene",

# group.by = "Sample",

proportion = FALSE,

cloneTypes=c(Single=1, Small=5, Medium=20, Large=100, Hyperexpanded=500))head(seurat)

因为很多barcode 并没有TCR clone ，因此cloneTypes 含有一些NA是正常的。

3.3 clonoType umap分布

可以使用DimPlot函数绘制 clonetype在UUMAP中的分布

colorblind_vector <- colorRampPalette(rev(c("#0D0887FF", "#47039FFF",

"#7301A8FF", "#9C179EFF", "#BD3786FF", "#D8576BFF",

"#ED7953FF","#FA9E3BFF", "#FDC926FF", "#F0F921FF")))p111 <- DimPlot(seurat, group.by = "cluster_name",label = T)p222 <- DimPlot(seurat, group.by = "cloneType",label = T) #NoLegend() scale_color_manual(values=colorblind_vector(5)) theme(plot.title = element_blank())p111 p222