EaSIeR包预测免疫治疗反应

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]]