预测癌症免疫治疗反应-EaSleR学习和知识整理

EaSIeR通过RNA-seq数据与其他不同来源的先验知识/数据相结合，得出系统性的肿瘤微环境特征以及量化了免疫细胞组成和细胞内和细胞间的通讯。开发团队对这种特征应用了多任务的学习方式赖预测免疫反应的不同标志物，并根据可解释的生物标志物得出特异性的癌症模型。

这种模型可以应用于接受 PD-1/PD-L1 抑制剂治疗的患者的RNA-seq数据。

输入数据：

1、Count数据；2、TPM数据；3、TMB数据；4、既往免疫治疗反应的数据(若存在的话)

分析步骤

1、导入数据和加载R包

注：示例数据不是RNA-seq数据， 仅用于实操练习！

rm(list = ls())
proj = "combat"
library("easier")
load("~/Desktop/data.Rdata")
head(exp)[1:4,1:4]
# GSM1020099 GSM1020100 GSM1020101 GSM1020102
# RFC2     7.727199   5.109329   8.021387   7.880099
# HSPA6    7.601476   6.538490   6.340193   6.016021
# PAX8     5.910106   5.848390   6.124806   5.894079
# GUCA1A   3.577948   3.591067   3.566542   3.553448

# TIDE预测得到了免疫治疗反应情况
meta_TIDE <- read.csv("~/Desktop/meta_TIDE.csv",
                      row.names = 1)
identical(colnames(exp),rownames(meta_TIDE))
# TRUE

2、获取免疫特征评分

Hallmarks <- c("CYT", "Roh_IS", "chemokines", "Davoli_IS", 
               "IFNy", "Ayers_expIS", "Tcell_inflamed", "RIR", "TLS")
scores <- compute_scores_immune_response(RNA_tpm = exp, 
                                         selected_scores = Hallmarks)
colnames(scores)
scores <- scores[,-c(5)]

# 对数据进行标准化
data_scaled <- t(scale(t(scores)))
pheatmap::pheatmap(data_scaled, 
                   cluster_rows = TRUE, 
                   cluster_cols = TRUE, 
                   show_rownames = TRUE, 
                   show_colnames = TRUE,
                   fontsize_row = 8, 
                   fontsize_col = 8)

3、对TME特征进行定量为后续分析做准备

# quanTIseq法对免疫细胞比例进行定量
cell_res <- compute_cell_fractions(RNA_tpm = exp)

# PROGENy法推断14条通路活性
# 部分通路相关基因可能被纳入计算免疫反应评分
# remove_sig_genes_immune_response设置为True，去除这部分重复基因进行计算，
# 这里是应该使用count数据
pathway_res <- compute_pathway_activity(RNA_counts = exp,
                                        remove_sig_genes_immune_response = TRUE)

# DoRothEA法进行转录因子分析
tf_res <- compute_TF_activity(RNA_tpm = exp)

# 肿瘤特异性细胞间通讯网络分析
# "BLCA"   "BRCA"   "CESC"   "CRC"    "GBM"    "HNSC"   "KIRC"   "LIHC"   
# "LUAD" # "LUSC"   "OV"     "PAAD"   "pancan" "PRAD"   "SKCM"  
# "STAD"   "THCA"   "UCEC"
# 其中pancan是指泛癌

lrpair_weights <- compute_LR_pairs(RNA_tpm = exp,
                                  cancer_type = "HNSC")
ccpair_scores <- compute_CC_pairs(lrpairs = lrpair_weights, 
                                   cancer_type = "HNSC")                                 

4、预测患者的免疫反应

# cancer_type需要特别指定
predictions <- predict_immune_response(pathways = pathway_activities,
                                       immunecells = cell_fractions,
                                       tfs = tf_activities,
                                       lrpairs = lrpair_weights,
                                       ccpairs = ccpair_scores,
                                       cancer_type = "HNSC", 
                                       verbose = TRUE)
# 若已知病人的免疫治疗结果，那么可以判断predictions的准确性，若未知免疫治疗效果，则预测的结果无从判断。

5、根据免疫治疗反应数据评估EaSIeR的预测

# 准备肿瘤突变负荷数据
# 由于这个数据集中是没有TMB数据的
# 自己分析的时候需要把TMB数据转化成一行，内容是TMB的值，列名是样本名
# TMB <- colData(dataset_mariathasan)[["TMB"]]
# names(TMB) <- colData(dataset_mariathasan)[["pat_id"]]

############
# 有免疫治疗数据评估免疫反应
output_have <- assess_immune_response(predictions_immune_response = predictions,
                                      patient_response = meta_TIDE$Responder,
                                      RNA_tpm = exp,
                                      TMB_values = TMB,
                                      easier_with_TMB = "weighted_average",
                                      weight_penalty = 0.5)
# Error: object 'TMB' not found
# In addition: Warning message:
# In assess_immune_response(predictions_immune_response = predictions,  :
#   Gene expression (TPM) data is missing, skipping computation of scores of immune response..

###################
# 无免疫治疗数据评估免疫反应
output_none <- assess_immune_response(predictions_immune_response = predictions,
                                      TMB_values = TMB,
                                      easier_with_TMB = "weighted_average",
                                      weight_penalty = 0.5)
# Error: object 'TMB' not found
# In addition: Warning message:
# In assess_immune_response(predictions_immune_response = predictions,  :
#   Gene expression (TPM) data is missing, skipping computation of scores of immune response..

################
# 获得免疫反应评分
derived_scores <- retrieve_easier_score(predictions_immune_response = predictions,
                                        TMB_values = TMB,
                                        easier_with_TMB = c("weighted_average", 
                                                             "penalized_score"),
                                        weight_penalty = 0.5)

该数据集中没有TMB数据，所以这里的很多步没法执行

6、可解释的Biomarker预测/展示免疫疗法的反应

identical(colnames(exp),rownames(meta_TIDE))
#meta_TIDE$threshold = factor(ifelse(meta_TIDE$TIDE > 0,'NR','R'),levels=c('R','NR'))
ICB_response = ifelse(meta_TIDE$TIDE > 0,'NR','R')
ICB_response <- t(ICB_response)
names(ICB_response) <- colnames(exp)
output_biomarkers <- explore_biomarkers(pathways = pathway_res,
                                        immunecells = cell_res,
                                        tfs = tf_res,
                                        lrpairs = lrpair_weights,
                                        ccpairs = ccpair_scores,
                                        cancer_type = 'HNSC',
                                        patient_label = ICB_response)
# 输出图片
output_biomarkers[[1]]
output_biomarkers[[2]]
output_biomarkers[[3]]
output_biomarkers[[4]]
output_biomarkers[[5]]
output_biomarkers[[6]]

重要结果解读

第一张图

个人理解：蓝色代表正相关，红色代表负相关，数值代表不同类型细胞在肿瘤中浸润的程度。Z-score就是展示同一种细胞亚型对不同分组的患者的影响分数。

第二张图

个人理解：左边是没有反应患者，右边是有反应患者情况。中间的灰色或者蓝色点的大小代表权重。形状代表正负相关性。里面的文字标识代表了不同的细胞-细胞，细胞-分子，分子-分子之间的作用情况/方向。整体从左往右看的话可以观察患者无反应到有反应过程中的分子调控的变化。

更详细的结果解释可读一下原文~

参考资料：

1、EaSleR: https://bioconductor.org/packages/release/bioc/vignettes/easier/inst/doc/easier_user_manual.html

2、Interpretable systems biomarkers predict response to immune-checkpoint inhibitors. Patterns (N Y). 2021 Jun 30;2(8):100293

2、生信菜鸟团： https://mp.weixin.qq.com/s/_Z5mWDNi8hzAc_mLivcjsQ https://mp.weixin.qq.com/s/AHJCgY_341ZY5pmq6vSQCw

致谢：感谢曾老师、小洁老师以及生信技能树团队全体成员。

注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。更多内容可关注公众号：生信方舟

- END -