ProTICS 揭示了不同分子亚型中肿瘤浸润免疫细胞的预后影响

生信技能树

发布于 2022-07-26 10:19:16

5730

发布于 2022-07-26 10:19:16

文章被收录于专栏：生信技能树

今天分享的是复旦大学和西北民族大学小伙伴合作的笔记

下面是笔记原文

1. 引入

这是生信技能树知识整理工作的第2个文档，前面的是：基于支持向量机模型的TNBC的分子亚型预测

在看到这个名字的时候，我本能的以为ProjTICS[1]应该也是个R包，但在进行诸多尝试之后才发现其实不然。

ProTICS 揭示了不同分子亚型中肿瘤浸润免疫细胞的预后影响

Github:https://github.com/liu-shuhui/ProTICS

[1] Shuhui Liu, Yupei Zhang, Xuequn Shang, Zhaolei Zhang, ProTICS reveals prognostic impact of tumor infiltrating immune cells in different molecular subtypes, Briefings in Bioinformatics, Volume 22, Issue 6, November 2021, bbab164, https://doi.org/10.1093/bib/bbab164

2. ProTICS 背景

同一种癌症的不同亚型往往表现出不同的基因组特征，需要有针对性的治疗。不同亚型肿瘤微环境的细胞和分子水平的差异对肿瘤的发病机制和预后有重要影响。虽然有关肿瘤浸润性淋巴细胞在选定组织学亚型中的预后关系的研究颇多，但很少有研究系统地报道如何通过多组学数据集使用机器学习方法量化免疫细胞在分子亚型中对预后的影响。

ProTICS （Prognostic association between Tumor-Infiltrating lymphocytes and Cancer Subtypes）描述了一种新的计算框架，用于量化肿瘤微环境中免疫细胞比例的差异，并估计其在不同亚型中的预后效应。首先，作者应用非负张量因子分解技术（nonnegative tensor factorization technique），根据基因表达和甲基化谱将患者分为不同的分子亚型。然后作者使用基于 mRNA 的反卷积（deconvolution）方法定量每个标本中细胞类型的比例。对于每个亚型的肿瘤，作者通过应用 Cox 比例风险回归（ Cox proportional hazard regression）估计免疫细胞类型的预后影响。在分子水平上，我们还预测了每个亚型的特征基因的预后。最后，作者对 ProTICS 在三个 TCGA 数据集和另一个独立的 METABRIC 数据集上的性能进行了基准测试。总的来说，ProTICS 成功地将肿瘤分为不同的分子亚型，表现为不同的总体存活率。

3. ProTICS的下载

一般的Github R 包的学习策略是直接访问Github 官网，ProTICS的使用方法也是直接写在了Github上，浏览其页面发现，ProTICS并没有相关的安装R包的代码，标准的安装github R包的方式对其均无效，这表明，ProTICS可能封装的并不完全。后续在其提供的说明文档中发现，ProTICS 调用函数的方式是 source，而非library R包后直接使用，这点也提示常规的方法在这里不适用。

source("./R/functions/normalization.R")
source("./R/functions/NTD_subtyping.R")

基于前述的探索，我们发现，ProTICS关键的部分在于其内置的几个函数，因此，将ProTICS从Github 上下载，并在环境下按提示使用Source 调用函数即可。

4. ProTICS 使用

ProTICS 是一个由三个部分组成的严格的pipeline，这三个部分有各自的目标。后面部分的实现依赖于前面部分的

4.1 ProTICS pipeline：Part 1

这里使用作者提供的演示数据进行初步的探索，输入数据需要提供甲基化水平和基因表达的数据结果。首先按作者提示配置好R包环境。

4.1.1 ProTICS pipeline 环境配置

#请在当前环境下加载如下R包
suppressMessages({
  library(data.table)
  library(dplyr)
  library(rTensor)
  library(nnTensor)
  library(survival)
  library(survminer) 
})

4.1.2 输入数据

# read gene expression data

data1<-fread(file = "./Data/data1.txt",header = T)
dim(data1)
## [1] 2000  201
head(data1)[1:5,1:5]
##    symbol TCGA-3C-AALJ TCGA-5L-AAT1 TCGA-5T-A9QA TCGA-A1-A0SP
## 1:   DPM1     46.28950     31.09528     52.98086     34.70152
## 2:  FUCA2     14.94863     20.04903     27.75997     43.72134
## 3: SEMA3F     36.01887     14.26484     48.44247     17.15032
## 4:    BAD     28.23475     22.86731     13.66222     17.29690
## 5:   LAP3     20.44311     77.64164     33.18052     33.60387
# read DNA methylation data
data2<-fread(file = "./Data/data2.txt",header = T)
dim(data2)
## [1] 2000  201
head(data2)[1:5,1:5]
##   Gene_Symbol TCGA-3C-AALJ TCGA-5L-AAT1 TCGA-5T-A9QA TCGA-A1-A0SP
##1:        DPM1    0.9397264    0.9248537    0.9620376    0.9112450
##2:       FUCA2    0.8777380    0.9305087    0.9422915    0.9562860
##3:      SEMA3F    0.6188444    0.6711587    0.4135552    0.8241335
##4:         BAD    0.7404509    0.7059423    0.7164417    0.8291710
##5:        LAP3    0.8147975    0.4907758    0.5874193    0.6480680
table(colnames(data1)== colnames(data2))

## FALSE  TRUE 
##     1   200 
clinicdata<-fread(file ="./Data/clinic_Data.txt",header = T)
colnames(clinicdata)<-c("patient_id", "death", "survival")
dim(clinicdata)
## [1] 200   3
head(clinicdata)
##     patient_id death survival
##1: TCGA-3C-AALJ     0     1474
##2: TCGA-5L-AAT1     0     1471
##3: TCGA-5T-A9QA     0      303
##4: TCGA-A1-A0SP     0      584
##5: TCGA-A2-A0CK     0     4159
##6: TCGA-A2-A0CR     0     3283
table(clinicdata$patient_id == colnames(data2)[2:ncol(data2)])
## TRUE 
## 200

需要注意的是：

这里的行名不是gene symbol 也不是patient id，和常规的认知需要相互区分。
三个输入数据具有的样本量、样本名及其顺序应是一致的

4.1.3 加载函数

source("./R/functions/normalization.R")
source("./R/functions/NTD_subtyping.R")

4.1.4 NTD (Nonnegative Tensor Decomposition,非负张量分解)运算及生存分析可视化

Subtype= NTD_subtyping(data1,data2,k=2, n=100)
# 简要查看NTD_subtyping函数code
NTD_subtyping
function(data1,data2,k,n){
  ### define a 3-mode tensor
  
  arr <- array(0,dim = c(dim(data1[,-1]),2)) # rows: genes; columns: patiens
  arrT <- as.tensor(arr)
  
  arrT[,,1] <- unlist(normalization(data1[,-1]))
  arrT[,,2] <- unlist(normalization(data2[,-1]))

  ##k:the number of subtypes; n:The number of interation step (Default: 100)
  
  output <- NTD(arrT, rank=c(k, k, k),num.iter=n)  # NTD function from nnTensor

  matrice_B<-t(output$A[[2]]) # matrice_B saved the latent factors information of patiens
  group<-max.col(matrice_B) # subtypes information
  return(group)
}
<bytecode: 0x00000182227e1920>
class(Subtype)
## [1] "integer"
table(Subtype)
## Subtype
##   1   2 
## 106  94

到这里，是完成了NTD非负张量分解的运算。输入参数 data1 和data2是前期准备的数据，k 对应后面 nnTensor R包的核心函数NTD的rank 参数， n 对应后续的 num.iter参数。在以往的学习中，我们经常看到的是 NMF，即非负矩阵分解。这里的NTD 是NMF在高维数组中的进一步推广[2]。这是一个很深层次的数学问题，临床专业的我实在是搞不定，放弃深入理解。

[2] Hao, N., Horesh, L., Kilmer, M.E. (2014). Nonnegative Tensor Decomposition. In: Carmi, A., Mihaylova, L., Godsill, S. (eds) Compressed Sensing & Sparse Filtering. Signals and Communication Technology. Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-642-38398-4_5

在NTD运算后，Subtype 即是NTD运算的分组信息，具体的分组数有k值决定。因此，作者在文献中提及

‘’To obtain the optimal number of subtypes, we traverse the number of subtypes from 2 to 6‘’ [1]

意即分别设置k=2,3,4,5,6，而后根据总体生存率最低值确定合适的亚型数。这里意在探索R包，简要使用k =2 完成后续代码分析。

下面比较两个亚群的生存结局

# 合并至临床数据
survivaldata<-cbind(clinicdata,Subtype)
# 输出保存
write.table(survivaldata, file = "./output/overallsurvival_subtypes.txt",
            sep = "\t", col.names = T, quote = F, row.names = F)
#生存分析
survdiff(Surv(survival,death)~Subtype, data=survivaldata)

survival_out<-survfit(Surv(survival,death)~Subtype, data=survivaldata)
#画图
ggsurvplot(survival_out, data = survivaldata, risk.table = T,
           xlab="Survival time/day", ylab="Survival rate")

4.2 ProTICS pipeline：Part 2

Part 1 主要完成对肿瘤数据的亚型进行区分，这里Part 2主要完成差异基因的分析和展示

4.2.1 加载所需R包

suppressMessages({
library(edgeR)
library(limma)
library(Glimma)
library(gplots)
library(org.Mm.eg.db)
library(grDevices)
library(data.table)
library(dplyr)
})

4.2.2 输入数据及数据处理

##读入count 矩阵
sig_expr <- fread("./Data/signature_count.txt",sep = "\t",
                  header = TRUE) ##行是signature genes
dim(sig_expr)  ##发现这里只有177个gene，确实是经过筛选的
## [1] 177 201  
survival_data <- fread("./output/overallsurvival_subtypes.txt", 
                       sep = "\t",header = TRUE) ### subtype_label can be saved from part1.R
subtypes<-survival_data$Subtype

ID<- which(subtypes==1 | subtypes==2)
Surv<-survival_data[ID,]
seqd<-select(sig_expr,c(colnames(sig_expr)[1],Surv$patient_id))
## Error in (function (classes, fdef, mtable)  : 
##  unable to find an inherited method for function ‘select’ for signature ‘"data.table"’

# 在这里出现一个经典报错，可以的解决方案是 dplyr::select
seqd<-dplyr::select(sig_expr,c(colnames(sig_expr)[1],Surv$patient_id))

不管是说明文档还是demo code 里面，都没有对这个signature gene 进行说明，但在作者的原文献中提及，TIP方法进行浸润分析用到的是178个预定义 signature genes，作者提供在github上的数据中少了一个gene。浏览了一下TIP （http://biocc.hrbmu.edu.cn/TIP/）的网页，目前的gene set 更新到273个了，但这并不影响目前这个探索结果。

In addition to CIBERSORT, another tool, TIP [32], adopts the same principle but is able to work with RNA-seq data and has compiled a set of 178 pre-curated signature genes‘’[2]

4.2.3 差异分析及可视化

source("./R/functions/subtypes_DEA.R")
GS<-subtypes_DEA(Surv,seqd)

# 差异基因热图展示

sig_expr<-sig_expr[is.element(sig_expr$symbol,GS),]

IDD<-c(which(subtypes==1),which(subtypes==2))
survd_new<-survival_data[IDD,]

sigdata<-dplyr::select(sig_expr,c(colnames(sig_expr)[1],survd_new$patient_id))

anno_c<-data.frame(Types = factor(survd_new$Subtype,c("1","2"),c("Sub1","Sub2")))
colnames(anno_c)<-c("  ")
row.names(anno_c)<-survd_new$patient_id

subtypes_DEA函数将limma 分析的流程封装在一步中完成，并了选取meaningful 的基因，选取标准是abs(logFC)>=1 & (adj.P.Val < 1e-2))，只保留 ≤ 20个 genes。后续是数据处理，标准化及热图展示，不多赘述。

source("./R/functions/normalization.R")
data<-normalization(log2(sigdata[,-1]+1))

library(pheatmap)
rownames(data)<-sigdata$symbol
pheatmap(data,cluster_rows=T,
         color = colorRampPalette(c( "#0077FF","#FFEEFF","#FF7700"))(1000),
         cluster_cols=F,show_rownames = TRUE,show_colnames=F,
         annotation=anno_c,annotation_legend=TRUE,main="dataset")

4.3 ProTICS pipeline：Part 3

Part 3 部分是这个pipeline 的一个小亮点，将不同subtype的免疫细胞与预后相关联

4.3.1 加载所需R包

library(forestplot)
library(data.table)
library(survival)
library("survminer")
library(dplyr)

4.3.2 输入数据及数据简要处理

survdata <- fread("./output/overallsurvival_subtypes.txt", sep = "\t",header = TRUE)
cell<-fread(file = "./Data/CellProportion.txt", sep = "\t",header = T)
head(cell)[1:5,1:5]
##        Mixture B cells CD4 Naive CD4 Memory CD8 Naive
##  1: TCGA-A1-A0SP       0     0.059      0.090         0
##  2: TCGA-EW-A1IZ       0     0.045      0.063         0
##  3: TCGA-UL-AAZ6       0     0.031      0.076         0
##  4: TCGA-A2-A3XU       0     0.062      0.109         0
##  5: TCGA-A7-A6VY       0     0.028      0.090         0
## clolumns[16:18] 不是细胞类型评分，需要去除.
cell<-cell[,-c(16:18)]
###  去除方差非常小的列
id=which(apply(cell[,-1],2,var)>1e-05)+1 
cell_new<-dplyr::select(cell,c(colnames(cell)[c(1,id)]))

### colnames of cell types
covariates<-c("`CD4 Naive`","`CD4 Memory`","`CD8 Memory`",
              "`CD8 Effector`", "`Th cell`", "`Monocytes CD16`",
              "`Monocytes CD14`","DC","pDC","Plasma")

# 这里也可以考虑 covariates <- colnames(cell_new)[2:ncol(cell_new)]

这里的细胞评分是通过TIP网页评分完成的。根据文献的描述，和TIP网页的提示，这里只需要讲178个signature genes 的log2 transform 的表达值输入即能完成。后续导入R中，完成后续的可视化步骤。

4.3.3 免疫细胞分布可视化

##### 1. 绘制10种免疫细胞在不同分子亚型中的比例分布
`Cell types` = c(rep(covariates, each=length(which(survdata$Subtype==1))),
                 rep(covariates, each=length(which(survdata$Subtype==2))))
`Patient type` = c(rep(c("Subtyp1"),each=length(which(survdata$Subtype==1))*10),
                   rep(c("Subtyp2"),each=length(which(survdata$Subtype==2))*10))

# 整理不同Subtype的 patientid
ID1<-sapply(survdata$patient_id[which(survdata$Subtype==1)],
                 function(x) which(cell_new$Mixture==x))
ID2<-sapply(survdata$patient_id[which(survdata$Subtype==2)],
                 function(x) which(cell_new$Mixture==x))

`Relative proportions of the 10 immune cell types` <-c(as.vector(as.matrix(cell_new[ID1,-1])),
                                                  as.vector(as.matrix(cell_new[ID2,-1])))

data<-data.frame(`Cell types`,`Patient type`,`Relative proportions of the 10 immune cell types`)

实际上,前面这几步，也可以考虑数据合并后，通过reshape::melt()完成

data$Cell.type <- factor(data$Cell.type,levels=covariates,ordered = TRUE)

            ggplot(data, aes(`Cell types`, y=`Relative proportions of the 10 immune cell types`, color=`Patient type`)) +
  theme(
    panel.background = element_rect(linetype = 1, colour = "white", size = 1,fill = "lightblue"),
    axis.text.x = element_text(angle = 20, hjust = 0.6,vjust = 0.75),
    plot.title = element_text(colour = "black",face = "bold",size = 12, vjust = 1),
    plot.margin = unit(c(0.2, 0.2, 0.2, 0.2), "inches")
  )+
  stat_boxplot(geom ='errorbar', width = 0.8) +
  geom_boxplot(width = 0.8) ##  +   
# facet_grid(.~Cell.type, scales = "free_x") # facet_grid 可以设置分面

4.3.4 单因素分析

##### 2. 免疫细胞预后森林图

surv_sub<-survdata[which(survdata$Subtype==1),] ## This runs for subtype 1;
surv_sub$survival<-scale(surv_sub$survival,center = FALSE, scale = TRUE)


ID<-sapply(surv_sub$patient_id, function(x) which(cell_new$Mixture==x))
cell_new<-cell_new[ID,-1]
#cell_new<-logcell<-log2(cell_new+1)
cutoff<-as.matrix(apply(cell_new,2,median))

tem<-t(replicate(dim(cell_new)[1],cutoff[,1]))
mat_bip<-as.matrix(cell_new>tem)
mat_bip[mat_bip==TRUE]<-1

data1<-cbind(surv_sub,mat_bip)
head(data1)[1:5,1:5]
##     patient_id death  survival Subtype CD4 Naive
## 1: TCGA-3C-AALJ     0 0.9465150       1         1
## 2: TCGA-5L-AAT1     0 0.9445885       1         0
## 3: TCGA-5T-A9QA     0 0.1945685       1         1
## 4: TCGA-A2-A0CK     0 2.6706619       1         0
## 5: TCGA-A2-A0SW     1 0.8765217       1         1

这里为了使用cox 分析，survival参数使用scale函数将其转变为连续型变量，后续，截断值也是取median，后续大于cutoff的定为1，否则为0。而后以suvival 需要注意的是，这里的数据处理只包含了 subtype = =1。

## 单因素回归
source("./R/functions/uni_cox.R")
result<-uni_cox(covariates,data1)
res1<-result[[1]]
res2<-result[[2]]
# forest plot
forestplot(res1, mean = res2$HR, lower = res2$lower, upper = res2$upper,
           graph.pos = 2,graphwidth = unit(18,"mm"),
           hrzl_lines = list("2" = gpar(lty=2,columns=1:4)),
           is.summary = c(TRUE,rep(FALSE,10)),
           txt_gp = fpTxtGp(ticks = gpar(cex=0.8),summary = gpar(cex=0.8),cex = 0.8),
           boxsize = 0.2,
           line.margin = unit(6,"mm"),
           lineheight = unit(6,"mm"),
           col=fpColors(box="blue",line="blue",summary="blue"),
           clip = c(0,5),
           xticks = c(0, 0.5, 1, 2,3,4,5),
           lwd.ci=2, ci.vertices=TRUE, ci.vertices.height = 0.12,
           colgap = unit(2,"mm"),zero = 1,
           title = "Subtype 1")

4.3.5 多因素分析

#### 多因素回归
source("./R/functions/multi_cox.R")
result<-multi_cox(covariates,data1)
res1<-result[[1]]
res2<-result[[2]]
### 绘图
forestplot(res1, mean = res2$HR, lower = res2$lower, upper = res2$upper,
           graph.pos = 2,graphwidth = unit(18,"mm"),
           hrzl_lines = list("2" = gpar(lty=2,columns=1:4)),
           is.summary = c(TRUE,rep(FALSE,10)),
           txt_gp = fpTxtGp(ticks = gpar(cex=0.8),summary = gpar(cex=0.8),cex = 0.8),
           boxsize = 0.2,
           line.margin = unit(6,"mm"),
           lineheight = unit(6,"mm"),
           col=fpColors(box="blue",line="blue",summary="blue"),
           clip = c(0,5),
           xticks = c(0, 0.5, 1, 2,3,4,5),
           lwd.ci=2, ci.vertices=TRUE, ci.vertices.height = 0.12,
           colgap = unit(2,"mm"),zero = 1,
           title = "Subtype 1")

5. 探索总结

前期第一个写的乳腺癌的分类包TNBC.CMS是通过预定义的基因集将乳腺癌分为四个亚群，输入数据是单一的转录组数据。这里的ProTICS相较于前者，输入数据是包含转录组数据、甲基化表达谱和临床预后数据，使用相较新的算法NTD（非负张量分解）将肿瘤患者的样本分成多个分子亚型。因此，相较于前者，其探索出的结果具有更多样性，但作者在文献中仅探讨了免疫细胞评分的差异。另外，作者在多个肿瘤样本中探索pipeline的结果可信性，表明这个pipeline在其他肿瘤中的可行性，总的来说，确实是一个很棒的分类包。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2022-06-17，如有侵权请联系 cloudcommunity@tencent.com 删除

tcp/ip