恶性肿瘤细胞可通过多种方式干扰肿瘤微环境的免疫细胞以减弱对肿瘤细胞的杀伤作用并诱导免疫细胞免疫耐受。而肿瘤微环境固有的复杂性，多细胞的动态特性，对获取免疫反应生物标志物和预测免疫治疗效果提出了巨大挑战。因此，作者期望使用 bulk RNA-seq数据与不同的先验知识来源（肿瘤浸润细胞，细胞内信号传导，TFs活性,细胞间通讯强弱）相结合，以得出基于系统的肿瘤微环境特征，量化免疫细胞组成以及细胞内和细胞间的通讯。作者通过将多任务学习应用于这些特征，预测免疫反应的不同特征，并基于可解释的Biomarker推导出癌症类型特异性模型。并将该模型应用于来自接受 PD-1/PD-L1 抑制剂治疗的癌症患者的独立 RNA-seq 数据，以证明EaSIeR 的方法可以准确预测治疗结果。
下面代码示例数据来自IMvigor210CoreBiologies包，参考：IMvigor210CoreBiologies包官方下载链接失效问题
代码语言：javascript复制BiocManager::install("easier")
library(IMvigor210CoreBiologies)
data(cds)
head(counts(cds))
head(fData(cds))
head(pData(cds))
meta <- pData(cds)
dataset_mariathasan <- get_Mariathasan2018_PDL1_treatment()
dataset_mariathasan
#获取PD-1治疗反应信息
# patient_ICBresponse <- dataset_mariathasan@colData@listData[["BOR"]]
# names(patient_ICBresponse) <- dataset_mariathasan@colData@listData[["pat_id"]]
library(SummarizedExperiment)
#获取PD-1治疗反应信息
patient_ICBresponse <- colData(dataset_mariathasan)[["BOR"]]
names(patient_ICBresponse) <- colData(dataset_mariathasan)[["pat_id"]]
table(meta[colData(dataset_mariathasan)[["pat_id"]],1])
# 获取肿瘤突变负荷信息
TMB <- colData(dataset_mariathasan)[["TMB"]]
names(TMB) <- colData(dataset_mariathasan)[["pat_id"]]
# 获取肿瘤类型
cancer_type <- metadata(dataset_mariathasan)[["cancertype"]]
# 获取基因表达矩阵（counts 和tpm）
RNA_counts <- assays(dataset_mariathasan)[["counts"]]
RNA_tpm <- assays(dataset_mariathasan)[["tpm"]]
# 为节约时间，随机抽样30例，设定随机数，以保证结果可重复
set.seed(2019)
pat_subset <- sample(names(patient_ICBresponse), size = 30)
patient_ICBresponse <- patient_ICBresponse[pat_subset]
TMB <- TMB[pat_subset]
RNA_counts <- RNA_counts[, pat_subset]
RNA_tpm <- RNA_tpm[, pat_subset]
# 考虑到部分基因有多个对应关系，需要进一步处理（保留作者给定的gene symbol）
genes_info <- easier:::reannotate_genes(cur_genes = rownames(RNA_tpm))
## 去除不支持的基因symbol
non_na <- !is.na(genes_info$new_names)
RNA_tpm <- RNA_tpm[non_na, ]
genes_info <- genes_info[non_na, ]
## 去除 entries that are withdrawn
RNA_tpm <- RNA_tpm[-which(genes_info$new_names == "entry withdrawn"), ]
genes_info <- genes_info[-which(genes_info$new_names == "entry withdrawn"), ]
## 找出重复基因
newnames_dup <- unique(genes_info$new_names[duplicated(genes_info$new_names)])
newnames_dup_ind <- do.call(c, lapply(newnames_dup, function(X) which(genes_info$new_names == X)))
newnames_dup <- genes_info$new_names[newnames_dup_ind]
## 检索重复基因的数据
tmp <- RNA_tpm[genes_info$old_names[genes_info$new_names %in% newnames_dup],]
## 删除重复基因的数据
RNA_tpm <- RNA_tpm[-which(rownames(RNA_tpm) %in% rownames(tmp)),]
## 整合重复基因的数据
dup_genes <- genes_info$new_names[which(genes_info$new_names %in% newnames_dup)]
names(dup_genes) <- rownames(tmp)
if (anyDuplicated(newnames_dup)){
tmp2 <- stats::aggregate(tmp, by = list(dup_genes), FUN = "mean")
rownames(tmp2) <- tmp2$Group.1
tmp2$Group.1 <- NULL
# 整理归纳到一个表达矩阵
RNA_tpm <- rbind(RNA_tpm, tmp2)
}
hallmarks_of_immune_response <- c("CYT", "Roh_IS", "chemokines", "Davoli_IS", "IFNy", "Ayers_expIS", "Tcell_inflamed", "RIR", "TLS")
immune_response_scores <- compute_scores_immune_response(RNA_tpm = RNA_tpm,
selected_scores = hallmarks_of_immune_response)
head(immune_response_scores)
#quanTIseq
cell_fractions <- compute_cell_fractions(RNA_tpm = RNA_tpm)
# 考虑到部分通路相关基因可能被纳入计算免疫反应评分，
# 因此参数remove_sig_genes_immune_response 设置为True， 去除这部分重复基因进行计算，
# 初次使用可以尝试使用效果。
pathway_activities <- compute_pathway_activity(RNA_counts = RNA_counts,
remove_sig_genes_immune_response = TRUE)
##====转录因子分析
tf_activities <- compute_TF_activity(RNA_tpm = RNA_tpm)
#> lrpair_weights 是个data.frame
lrpair_weights <- compute_LR_pairs(RNA_tpm = RNA_tpm,
cancer_type = "pancan")
#> LR signature genes found in data set: 629/644 (97.7%)
#> Ligand-Receptor pair weights computed
head(lrpair_weights)[,1:5]
ccpair_scores <- compute_CC_pairs(lrpairs = lrpair_weights,
cancer_type = "pancan")
head(ccpair_scores)
# 这个部分 cancer_type需要特别指定，需要注意的是，前面cancer_type 默认值都是 'pancan' == ‘泛癌’
# 但在这个部分 'pancan' 却不能作为其中一个参数，而需要特别指定'BRCA'（膀胱癌样本）
predictions <- predict_immune_response(pathways = pathway_activities,
immunecells = cell_fractions,
tfs = tf_activities,
lrpairs = lrpair_weights,
ccpairs = ccpair_scores,
cancer_type = cancer_type,
verbose = TRUE)
predict_immune_response
#有反应数据预测
output_eval_with_resp <- assess_immune_response(predictions_immune_response = predictions,
patient_response = patient_ICBresponse,
RNA_tpm = RNA_tpm,
TMB_values = TMB,
easier_with_TMB = "weighted_average",
weight_penalty = 0.5)
# inspect output
output_eval_with_resp[[1]]
output_eval_with_resp[[2]]
output_eval_with_resp[[3]]
#无反应数据预测
output_eval_no_resp <- assess_immune_response(predictions_immune_response = predictions,
TMB_values = TMB,
easier_with_TMB = "weighted_average",
weight_penalty = 0.5)
# 图片输出
output_eval_no_resp[[1]]
output_eval_no_resp[[2]]
#获得免疫反应评分
output_biomarkers <- explore_biomarkers(pathways = pathway_activities,
immunecells = cell_fractions,
tfs = tf_activities,
lrpairs = lrpair_weights,
ccpairs = ccpair_scores,
cancer_type = 'BLCA',
patient_response = patient_ICBresponse)
# 输出图片
output_biomarkers[[1]]
output_biomarkers[[2]]
output_biomarkers[[3]]
output_biomarkers[[4]]
output_biomarkers[[5]]
output_biomarkers[[6]]
