单细胞免疫组库VDJ|和Nature学STARTRAC，定量T细胞动态变化

是发表于2018年的 文章（ of T cells in ）中的分析方法，可以应用于单细胞免疫组库数据来揭示T细胞动态变化的分析。原理假设认为克隆型一致的细胞来源一致，可以定量刻画T细胞的组织分布、克隆扩增、组织迁移和状态变化等。

5'/>

上图中不同颜色的圆球代表不同的T细胞类型，圆球上不同颜色的“Y”代表了不同的TCR克隆型，右边给出了简单的算法。

其中指不同T细胞在某个细胞分群中的克隆程度；指相同克隆型的T细胞在不同组织间的扩散程度；指相同克隆型的T细胞在不同细胞类型之间的共享程度。

简单的了解一下原理以及指标的含义，实现的话就相对比较简单了。

一 准备R包，数据

首先上加载R包和示例数据，然后将我们自己的数据整理成示例数据的格式，然后运行的话只需要一行代码即可。

#install.packages("devtools")#devtools::install_github("Japrin/STARTRAC")
library("Startrac")library("tictoc")library("tidyverse")library("Seurat")library("data.table")library("ggpubr")library("ComplexHeatmap")library("RColorBrewer")library("circlize")
dat.file <- system.file("extdata/example.cloneDat.Zhang2018.txt",package = "Startrac")in.dat <- read.table(dat.file,stringsAsFactors = F,head=T)#run the STARTRAC pipelineout <- Startrac.run(in.dat, proj="CRC", cores=NULL,verbose=F)#查看示例数据head(in.dat,2)#Cell_Name            clone.id clone.status patient sampleType stype majorCluster loc#1  TTH36-20180123 CRC.P0123_C000002:9       Clonal   P0123        TTH   CD4 CD4_C07-GZMK   T#2 TP7170-20180123 CRC.P0123_C000002:9       Clonal   P0123        TP7   CD4 CD4_C07-GZMK   T

可以看到包含 样本的基本信息（名称，类型，位置），clone相关信息（ ID，clone ID，clone 状态（是否是clone）等），以及单细胞细胞类型注释的信息（CD4,CD8 ，亚型）。

下面就需要将我们自己的VDJ数据 + 单细胞数据整理成这样的格式，其中样本信息（已知），细胞注释信息（单细胞免疫组库VDJ| 从零开始分析，解决真实场景中可能的问题）有，现在需要解决clone的ID 和 状态即可。

二 VDJ数据处理

2.1 VDJ数据合并

首先将上篇推文单细胞免疫组库VDJ| 从零开始分析，解决真实场景中可能的问题中提到的所有VDJ文件 合并在一起，可以linux中cat ，可以excel 中复制粘贴，可以R中一个个读入然后rbind ，也可以循环合并（注意保留样本名），最终效果如下

#添加file 标签read_tcr <- function(tcrfile){  p3_n <- read.csv(tcrfile)  p3_n$file <- sub('.filtered_contig_annotations.csv','',sub('^.*/','',tcrfile))  return(p3_n)}
tcrfiles <- list.files('./','.filtered_contig_annotations.csv',full.names = T)tcrfiles
if (all(file.exists(tcrfiles))){  tcr_list = list()  for (i in 1:length(tcrfiles)){    print(i)    tcr_list[[i]] = read_tcr(tcrfile = tcrfiles[i])  }}lapply(tcr_list,  dim)
vdj <- do.call(rbind, tcr_list) ; dim(vdj)head(vdj,2)table(vdj$file)

5__细胞免疫测定'/>

2.2 VDJ数据过滤

使用 为true，且为true的TRA TRB的序列 ，通过合并样本名+构建唯一

vdj <- vdj %>%   dplyr::filter(high_confidence =="true" &                   chain %in% c("TRA","TRB") &                  productive =="true")vdj$Cell_name <- paste0(vdj$file,'_',vdj$barcode)head(vdj,2)

注：true这里可能是True 也可能是TRUE，注意进行对应的修改

5_细胞免疫测定_'/>

2.3 拆分/合并 TRA ，TRB

前面也提到了clone一般是结合TRA 和 TRB的cdr3序列 ，因此这里先拆分TRA 和 TRB ，以备后面合并使用

vdj_a <- vdj %>% filter(chain =="TRA") %>% dplyr::arrange(desc(umis), desc(reads)) vdj_b<-vdj%>%filter(chain=="TRB")%>%dplyr::arrange(desc(umis),desc(reads))#### Get the best TRA or TRB test <- vdj_a %>%   dplyr::group_by(Cell_name) %>%   dplyr::summarise(reads=max(reads), umis=max(umis)) head(test)vdj_a <- data.frame(inner_join(vdj_a, test)) #Joining, by = c("reads", "umis", "Cell.name") 按照3列 join ，所以是最大的dim(vdj_a)
test <- vdj_b %>% group_by(Cell_name) %>%  dplyr::summarise(reads = max(reads), umis=max(umis) )vdj_b <- data.frame(inner_join(vdj_b, test))dim(vdj_b)

按照合并TRA 和 TRB

### merge TRA or TRB  final_vdj = dplyr::full_join(x = vdj_a, y=vdj_b, by = c("Cell_name"), suffix = c(".TRA",".TRB"))dim(final_vdj)head(final_vdj,2)save(final_vdj,file = 'final_vdj.rda')

三结合单细胞转录组数据

3.1 合并单细胞数据

单细胞数据同样需要构建与VDJ结果一致的唯一列，然后进行合并。

subT <- get(load("E:/bioinformation/scTCR_BCR/seurat_T.RData") )subT@meta.data <- subT@meta.data %>%   mutate(Cell_name = rownames(subT@meta.data)) %>%   inner_join(final_vdj, by = "Cell_name")
head(subT@meta.data)

3.2 计算Clone信息

结合TRA 和TRB的cdr3序列构建clone ，并统计每种clone的个数

subT@meta.data$Clone_AA = paste(subT@meta.data$cdr3.TRA, subT@meta.data$cdr3.TRB, sep="_")
subT@meta.data = subset(subT@meta.data, productive.TRA == "true" & productive.TRB == "true"  ) ; dim(subT@meta.data)subT@meta.data = subT@meta.data %>% arrange(., Clone_AA)
### calculate clone number and clone IDtmp = subT@meta.data %>%   group_by(Clone_AA) %>%  summarize(Clone_NUM = n()) %>%  mutate(Clone_ID = paste0("Clone_",rownames(.)))head(tmp)
# A tibble: 6 × 3#  Clone_AA                       Clone_NUM Clone_ID#                                    #1 CAAAAAGKSTF_CASSQGDSSYEQYF             1 Clone_1 #2 CAAAAAGRRALTF_CSARGGWGGITGELFF         1 Clone_2 #3 CAAAANYGGATNKLIF_CASSLEYNEQFF          2 Clone_3 #4 CAAADGQKLLF_CASSYNSNQPQHF              1 Clone_4 #5 CAAADNYGQNFVF_CASSESSPEQFF             1 Clone_5 #6 CAAADSGGSEKLVF_CASSGLMNTGELFF          1 Clone_6

subT@meta.data = merge.data.frame(subT@meta.data, tmp) head(subT@meta.data,2)

3.3 根据示例数据筛选列

subT@meta.data中有很多信息，根据示例数据筛选出来对应的信息，并修改列名字 。

（1）根据拆分出CD4和CD8；

（2） 大于1，即为

subT.meta <- subT@meta.data %>%   select(Cell_name,Clone_ID,Clone_NUM,orig.ident,Sample,type,cluster,cluster_name,pos)head(subT.meta)
subT.meta$stype <- ifelse(subT.meta$cluster_name %in% c("CD4+ Activated IEG","CD4+ Effector","CD4+ Naive","CD4+ Proliferating","CD4+ Treg"),"CD4","CD8")subT.meta$clone.status <- ifelse(subT.meta$Clone_NUM >1 ,"Clonal","NoClonal")
subT.meta <- subT.meta %>%   select(Cell_name, Clone_ID ,clone.status, orig.ident ,Sample   ,stype , cluster_name , pos )names(subT.meta) <- names(in.dat)save(subT.meta,file = "subT.meta.Rdata")

5_细胞免疫测定'/>

保存结果，后台回复即可获取.rda 和 subT.meta.Rdata文件。

四 分析

准备好了subT.meta文件，分析就是一行代码的事情

tic("Startrac.run")out2 <- Startrac.run(subT.meta, proj="CRC",verbose=F)#plot(out2,index.type="cluster.all",byPatient=T)

可以输出结果，但是在按照官网文档使用plot的相关函数时候会报错。影响不大，可以自己提取数据绘制 或者 直接参考官网的函数 。可以先str(out2) 看一下数据结构，，和的结果可以对应的进行提取。

5_'/>

4. level

ggboxplot(as.data.table(out2@cluster.sig.data)[,][order(majorCluster),],          x="majorCluster",y="value",palette = "npg",          color = "index", add = "point", outlier.colour=NULL) +  facet_wrap(~index,ncol=1,scales = "free_y") +  theme(axis.text.x=element_text(angle = 60,hjust = 1))

5_细胞免疫测定_'/>

4. level index ofall data

dat.plot <- as.data.table(out2@cluster.sig.data)[aid==out2@proj,]ggbarplot(dat.plot[order(majorCluster),],               x="majorCluster",y="value",palette = "npg",fill = "index") +  facet_wrap(~index,ncol=1,scales = "free_y") +  coord_cartesian(clip="off") +theme(axis.text.x=element_text(angle=60,hjust=1),strip.background=element_blank())

5'/>

4.3 index two major

dat.plot <- as.matrix(subset(out2@pIndex.tran,aid==out2@proj)[,c(-1,-2,-3)])rownames(dat.plot) <- subset(out2@pIndex.tran,aid==out2@proj)[,3]dat.plot[is.na(dat.plot)] <- 0yrange <- pretty(dat.plot)col.heat <- colorRamp2(seq(0,max(yrange),length=15),                       colorRampPalette(rev(brewer.pal(n=7,name="RdBu")))(15),                       space = "LAB")Heatmap(dat.plot,name="pIndex.tran",col = col.heat)

5_细胞免疫测定_'/>

当时使用的还比较少，而TCR的定量刻画又很有意义，你确定不在文章中试试？

后面会分享一下发表在2021年 的Pan- -cell of tumor- T cells文章中使用的相关指数与 “目标指数”之间的相关分析内容。

5'/>

参考资料

://///blob///.html

◆◆◆ ◆◆

精心整理（含图PLUS版）|R语言生信分析，可视化（R统计，绘图，生信图形可视化汇总）