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

生信补给站

发布于 2023-08-25 10:15:33

3.7K0

发布于 2023-08-25 10:15:33

文章被收录于专栏：生信补给站

经过前面单细胞免疫组库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

  contigs     values total   scaled
1     795  t2_Normal  1221 65.11057
2    1672    t2_PBMC  1903 87.86127
3     759  t3_Normal  1044 72.70115
4    4462    t3_PBMC  6695 66.64675
5    2387  t3_Center  7615 31.34603
6    1159  t4_Normal  1575 73.58730
7    2661    t4_PBMC  3274 81.27673
8    2403  t4_Center  5395 44.54124
9    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               9 AAACCTGCAACACCTA-1
#t2_Normal_AAACCTGCACCTCGTT-1  t2_Normal       2194         1188  10.027347               7               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      13                       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

可以大概看出来CD8细胞中clone更多，符合文献的研究结果。

3.4 highlightClonotypes 高亮特定克隆型

使用highlightClonotypes函数展示特定序列的克隆型分布

seurat <- highlightClonotypes(seurat, 
                              cloneCall= "aa", 
                              sequence = c("CAVGGRSNSGYALNF;CAASRNNNDMRF_CASSSLGNEQFF", 
                                           "CAASRNNNDMRF_CASSSLGNEQFF"))
DimPlot(seurat, group.by = "highlight") + 
  theme(plot.title = element_blank())

3.5 occupiedscRepertoire cloneType分布

查看CD8T和 CD4T细胞中cloneType的分布情况

seurat_T <- subset(seurat, cluster_name %in% c("CD8A+ Exhausted","CD8A+ Exhausted IEG", 
               "CD8A+ NK-like","CD8A+ Proliferating","CD8A+ Tissue-resident",
               "CD4+ Activated IEG","CD4+ Effector" ,"CD4+ Naive","CD4+ Proliferating" ,
               "CD4+ Treg") )

occupiedscRepertoire(seurat_T, x.axis = "cluster_name")

3.6 alluvialClonotypes

使用alluvialClonotypes函数绘制T细胞种各个celltype之间的clonetype信息的桑葚图

alluvialClonotypes(seurat_T, cloneCall = "gene", 
                   y.axes = c("cluster_name", "cloneType"), 
                   color = "cluster_name")

3.7 getCirclize

使用getCirclize函数绘制circos图

library(circlize)
library(scales)
circles <- getCirclize(seurat, 
                       cloneCall = "gene+nt",
                       groupBy = "cluster_name"
                       )

四 scRepertoire基于celltype

第二部分中是根据sample计算的clone的一系列指标，如克隆型分布，长度，空间稳态等。那经过第三部分中和单细胞转录组结合后就可以根据celltype维度进行clone的上述所有统计，这里仅示例几个。

4.1 按照celltype拆分

首先使用expression2List函数将seurat按照cluster_name进行拆分，然后就可进行二中的所有分析。

注意这里不是官网的split.by ，而要使用group

combined2 <- expression2List(seurat_T, 
                             group = "cluster_name")
length(combined2) 
p01 <- clonalDiversity(combined2, 
                cloneCall = "nt")
p02 <- clonalHomeostasis(combined2, 
                  cloneCall = "nt")
p03 <- clonalOverlap(combined2, 
              cloneCall="aa", 
              method="overlap")
p01
p02
#保存数据(后台获取)，以备后用
save(combined,seurat_T,file = "combined_seurat_T.RData")