本文首发于“生信补给站”公众号 https://mp.weixin.qq.com/s/Lnba8cBT7RpA1bhwx5ppkQ
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)
当时使用Startrac的还比较少，而TCR的定量刻画又很有意义，你确定不在文章中试试？
后面会分享一下发表在2021年Science 的Pan-cancer single-cell landscape of tumor-infiltrating T cells文章中使用Startrac的相关指数与 “目标指数”之间的相关分析内容。
