單細胞免疫組庫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統計，繪圖，生信圖形可視化彙總）