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

生信补给站

发布于 2023-08-25 10:17:53

1.7K0

发布于 2023-08-25 10:17:53

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

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

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

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

一准备R包，数据

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

#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 pipeline
out <- 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相关信息（barcode ID，clone ID，clone 状态（是否是clone）等），以及单细胞细胞类型注释的信息（CD4,CD8 ，亚型）。

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

二 VDJ数据处理

2.1 VDJ数据合并

首先将上篇推文单细胞免疫组库VDJ| 从零开始scRepertoire分析，解决真实场景中可能的问题中提到的所有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)

2.2 VDJ数据过滤

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

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，注意进行对应的修改

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)

按照Cell_name合并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结果一致的唯一Cell_name列，然后进行合并。

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 ID
tmp = 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
#  <chr>                              <int> <chr>   
#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）根据celltype拆分出CD4和CD8；

（2）Clone_NUM 大于1，即为Clonal

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")

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

四 Startrac分析

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

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

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

4.1 cluster level index by patients

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))

4.2 cluster level index of all 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())

4.3 transition index between two major clusters

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)] <- 0
yrange <- 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)