【收藏贴】IOBR-让你轻松玩转多组学免疫分析！

##需要R版本大于等于3.6.3##安装IBOR以及一些分析依赖的包if (!requireNamespace("IOBR", quietly = TRUE)) remotes::install_github("IOBR/IOBR")BiocManager::install("maftools")if (!requireNamespace("EPIC", quietly = TRUE)) devtools::install_github("GfellerLab/EPIC", ref="master")if (!requireNamespace("estimate", quietly = TRUE)){ rforge <- "http://r-forge.r-project.org" install.packages("estimate", repos=rforge, dependencies=TRUE)}library(IOBR)library(EPIC)library(estimate) library(tidyverse)library(tidyHeatmap)library(maftools)library(ggpubr)library(ggplot2)

library(UCSCXenaTools)##下载TCGA胃癌RNA-seq数据eset_stad<-XenaGenerate(subset = XenaCohorts =="GDC TCGA Stomach Cancer (STAD)") %>%  XenaFilter(filterDatasets = "TCGA-STAD.htseq_counts.tsv") %>%  XenaQuery() %>% XenaDownload() %>%  XenaPrepare()##展示eset_stad[1:5,1:5]

#去掉版本号eset_stad$Ensembl_ID<-substring(eset_stad$Ensembl_ID, 1, 15)eset_stad<-column_to_rownames(eset_stad, var = "Ensembl_ID")# 去标准化eset_stad<-(2^eset_stad)-1eset_stad<-count2tpm(countMat = eset_stad, idType = "Ensembl", org="hsa")

• GEOquery")# 下载数据eset_geo<-getGEO(GEO     = "GSE100935",getGPL = F,destdir = "./")
  

eset <-eset_geo[[1]]eset <-exprs(eset)head(anno_hug133plus2)eset<-anno_eset(eset = eset,annotation = anno_hug133plus2,symbol = "symbol",probe = "probe_id",method = "mean")

data("eset_stad")##CIBERSORTcibersort<-deconvo_tme(eset = eset_stad, method = "cibersort", arrays = FALSE, perm = 200 )

data("imvigor210_sig")data("imvigor210_pdata")

pdata_group<-imvigor210_pdata[!imvigor210_pdata$BOR_binary=="NA",c("ID","BOR","BOR_binary")]##肿瘤特征tumor_signatureres<-iobr_cor_plot(pdata_group = pdata_group,id1 = "ID",feature_data = imvigor210_sig,id2 = "ID",target = NULL,group = "BOR_binary",is_target_continuous = FALSE,padj_cutoff = 1,index = 1,category = "signature",signature_group = sig_group[c(1,3,5)],ProjectID = "IMvigor210",palette_box = "paired1",palette_corplot = "pheatmap",palette_heatmap = 2,feature_limit = 26,character_limit = 30,show_heatmap_col_name = FALSE,show_col = FALSE,show_plot = TRUE,path = "1-BOR-relevant-signatures")

##CD_8_T_effectorres<-iobr_cor_plot(pdata_group = pdata_group,id1 = "ID",feature_data = imvigor210_eset,id2 = "ID",target = NULL,group = "BOR_binary",is_target_continuous = FALSE,padj_cutoff = 1,index = 1,category = "gene",signature_group = signature_collection[c(1:2,4)],ProjectID = "IMvigor210",palette_box = "paired1",palette_corplot = "pheatmap",palette_heatmap = 4,feature_limit = 26,character_limit = 30,show_heatmap_col_name = FALSE,show_col = FALSE,show_plot = TRUE,path = "2-BOR-relevant-genes")

##tumor_signaturepdata_group<-as.data.frame(pdata_group[,c("ID","HCP5","LINC00657")])head(pdata_group)res<-iobr_cor_plot(pdata_group = pdata_group,id1 = "ID",feature_data = imvigor210_sig,id2 = "ID",target = "HCP5",group = "group3",is_target_continuous = TRUE,padj_cutoff = 1,index = 1,category = "signature",signature_group = sig_group[1:3],ProjectID = "IMvigor210",palette_box = "set2",palette_corplot = "pheatmap",palette_heatmap = 2,feature_limit = 26,character_limit = 30,show_heatmap_col_name = FALSE,show_col = FALSE,show_plot = TRUE,path = "3-HCP5-relevant-signatures")

pdata_group<-as.data.frame(imvigor210_pdata[,c("ID","Pan_F_TBRs")])pdata_group$Pan_F_TBRs<-scale(as.numeric(pdata_group$Pan_F_TBRs))head(pdata_group)res<-iobr_cor_plot(pdata_group = pdata_group,id1 = "ID",feature_data = imvigor210_sig,id2 = "ID",target = "Pan_F_TBRs",group = "group3",is_target_continuous = TRUE,padj_cutoff = 1,index = 5,category = "signature",signature_group = sig_group[1:2],ProjectID = "IMvigor210",palette_box = "set2",palette_corplot = "pheatmap",palette_heatmap = 2,feature_limit = 26,character_limit = 30,show_heatmap_col_name = FALSE,show_col = FALSE,show_plot = TRUE,path = "5-Pan_F_TBRs-relevant-signatures")

res<-iobr_cor_plot(pdata_group = pdata_group,id1 = "ID",feature_data = imvigor210_sig,id2 = "ID",target = "Pan_F_TBRs",group = "group3",is_target_continuous = TRUE,padj_cutoff = 1,index = 6,category = "signature",signature_group = sig_group[20:24],ProjectID = "IMvigor210",palette_box = "jco",palette_corplot = "pheatmap",palette_heatmap = 3,feature_limit = 26,character_limit = 30,show_heatmap_col_name = FALSE,show_col = FALSE,show_plot = TRUE,path = "6-Pan_F_TBRs-relevant-TME-cell")

#下载突变数据maf_file<-"TCGA.STAD.mutect.c06465a3-50e7-46f7-b2dd-7bd654ca206b.DR-10.0.somatic.maf"mut_list<-make_mut_matrix(maf = maf_file, isTCGA = T,                           category = "multi")var_stad<-XenaGenerate(subset = XenaCohorts =="GDC TCGA Stomach Cancer (STAD)") %>%  XenaFilter(filterDatasets = "TCGA-STAD.mutect2_snv.tsv") %>%  XenaQuery() %>% XenaDownload() %>%   XenaPrepare() mut_list2<-make_mut_matrix(mut_data = var_stad, category = "multi", Tumor_Sample_Barcode = "Sample_ID", Hugo_Symbol = "gene", Variant_Classification = "effect", Variant_Type = "Variant_Type")

mut<-mut_list$snpres<-find_mutations(mutation_matrix = mut, signature_matrix = tcga_stad_sig,id_signature_matrix = "ID",signature = "CD_8_T_effector",min_mut_freq = 0.01,plot = TRUE,method = "Wilcoxon",save_path = paste0("7-CD_8_T_effector-relevant-mutations"),palette = "jco",show_plot = T,width = 8, height = 4,oncoprint_group_by = "mean",oncoprint_col = "#224444",gene_counts = 10)

data("imvigor210_sig")data("imvigor210_pdata")# 用于分析二元变量input<-imvigor210_pdata %>%  dplyr::select(ID,BOR_binary) %>%  inner_join(.,imvigor210_sig,by="ID") %>% filter(!is.na(.$BOR_binary)) %>% filter(!.$BOR_binary=="NA")
# Feature engineeringres<-batch_wilcoxon(data = as.data.frame(input), target = "BOR_binary", group_names = c("NR","R"), feature = colnames(input)[3:ncol(input)])head(res)
model_feas<-as.character(res[res$p.value<0.05,]$sig_names)
input<-as.data.frame(imvigor210_sig)feas<-colnames(input)[colnames(input)%in%model_feas]input<-input[, c("ID",feas)]
# 目标数据pdata_group <- imvigor210_pdata[!imvigor210_pdata$BOR_binary=="NA",c("ID","BOR_binary")]pdata_group$BOR_binary <- ifelse(pdata_group$BOR_binary == "R", 1, 0)
#特征选择binomial_result <- BinomialModel(x = input,  y = pdata_group,  seed = "123456",  scale = TRUE, train_ratio = 0.7,  nfold = 10,  plot = T)

plot(binomial_result$lasso_result$model)