For precision oncology, it is essential to identify which patients are likely to respond to a given therapy. This is particularly important for targeted therapies and immunotherapies, including immune checkpoint blockade, where molecular profiling of tumor samples can provide clinically relevant information. Here, we analyze bulk RNA-sequencing data from a colorectal cancer cohort from The Cancer Genome Atlas (TCGA-CRC) (n=621), integrated with clinical patient information including overall survival and consensus molecular subtype classification. We aim to determine which hallmark biological processes, immunological pathways, and cell-type signatures are enriched in individual patients and across consensus molecular subtypes. Furthermore, we investigate acquired and intrinsic MEK inhibitor resistance signatures for Selumetinib, together with mutations in the MAPK/ERK signaling pathway, to assess their potential predictive value for patient therapy response.

1. Load packages

Include all R and Bioconductor packages

library(conflicted)
library(readr)
library(tidyverse)
library(writexl)
#install.packages("remotes")
#remotes::install_github("omnideconv/immunedeconv")
library(immunedeconv)
#install.packages("devtools")
#devtools::install_github("bmbolstad/preprocessCore")
#devtools::install_github("Moonerss/CIBERSORT")
library(CIBERSORT)
library(GSVA)
library(GSEABase)
library(ComplexHeatmap)
library(circlize)
library(ggbeeswarm)
library(ggpubr)
library(survival)
library(survminer)
conflicts_prefer(dplyr::filter)
conflicts_prefer(dplyr::select)
conflicts_prefer(dplyr::slice_head)

2. Working directory and input files

Set your preferred working directory and ensure that the required input files, CRC_CLIN.txt, CRC_PRIM_TPM.txt, GENESETS.gmt, LM22.txt, and MEKi_selumetinib_CRC_signatures.gmt (from previous analyses) are stored in this directory.

  setwd("L:/2026/WM7/R") # use your preferred directory

3. Processing of bulk RNA sequencing and clinical information from TCGA-CRC cohort

Clinical information (Clinical_Pick_Tier1) and gene expression data (illuminahiseq_rnaseqv2-RSEM_genes; RSEM-processed RNA-sequencing V2 data) were downloaded from the COADREAD project via Firebrowse/Broad Institute. Normalized expression values were calculated as transcripts per million (TPM; scaled_estimate × 10^6) and log2-transformed after adding a pseudocount of 1. The final clinical annotation table was stored as CRC_CLIN.txt, and primary tumor expression data were stored as CRC_PRIM_TPM.txt, comprising a total of 621 patients.

Expression data (TPM) were loaded from CRC_PRIM_TPM.txt and log2 transformed and also converted into matrix fromat.

CRC_TPM <- read_tsv("CRC_PRIM_TPM.txt")
CRC_TPM_MAT <- column_to_rownames(CRC_TPM, var = "ID")
CRC_LOG2TPM1 <- CRC_TPM |>
  mutate(across(where(is.numeric), \(x) log2(x + 1)))
CRC_LOG2TPM1_MAT <- column_to_rownames(CRC_LOG2TPM1, var = "ID")

Expression data CRC_TPM has this format:

CRC_TPM |>
  select(1:8) |>
  slice_head(n = 20)

Clinical data were also loaded

CRC_CLIN <- read_tsv("CRC_CLIN.txt")
CRC_CLIN

4. Estimation of cell type fraction from the mixed bulk RNA sequencing using deconvolution approaches

For each tumor sample expression data (note these analyses need TPM as input) are used and two deconvolution methods (CIBERSORT and quanTIseq) are used. Information on the cell type fraction fro each patient are stored and visualized as stacked barplot. (Note for CIBERSORT the matrix LM22.txt is included)

CIBERSORT analyses:

sig_matrix <- system.file("extdata", "LM22.txt", package = "CIBERSORT")
results <- cibersort(sig_matrix, "CRC_PRIM_TPM.txt")
R<-as.data.frame(results)
write.table(R,"CIBERSORT_RESULTS.txt", sep="\t",quote=FALSE,row.names=TRUE, col.names=NA)
CB<-data.frame(Cell_type=rep(names(R[,1:22]),nrow(R)),PAT=rep(row.names(R), each=22), PRO=c(t(R[,1:22])))
cbp<-ggplot(CB, aes(x = PAT,y = PRO, fill = Cell_type)) +  
  geom_bar(stat = "identity") + 
  theme(axis.title.x=element_blank(),
       axis.text.x=element_blank(),
       axis.ticks.x=element_blank()) +
  theme(legend.text = element_text(size=10)) +
  theme(legend.title = element_text(size=10)) +
  theme(legend.key.size = unit(0.3, 'cm')) + 
  theme(plot.margin = margin(1,1,1,1, "cm")) +
  scale_fill_manual(values=rainbow(22)) +
  ggtitle("CIBERSORT")
ggsave("CIBERSORT.png",width = 24,height = 8, units = "cm", dpi = 300)
Estimated immune cell fractions inferred using CIBERSORT.
Estimated immune cell fractions inferred using CIBERSORT.

QUANTISEQ analyses:

Q<-deconvolute(CRC_TPM_MAT, "quantiseq")
write.table(Q, file="QUANTISEQ_RESULTS.txt", sep="\t",quote=FALSE,row.names=FALSE)
D<-data.frame(Cell_type=rep(Q$cell_type,(ncol(Q)-1)),PAT=rep(names(Q)[2:ncol(Q)], each=nrow(Q)),PRO=as.numeric(data.frame(col =unlist(Q))[(nrow(Q)+1):(nrow(Q)*ncol(Q)),1]))
gp<-ggplot(D, aes(x = PAT,y = PRO, fill = Cell_type)) +  
  geom_bar(stat = "identity") + 
  theme(axis.title.x=element_blank(),
       axis.text.x=element_blank(),
       axis.ticks.x=element_blank()) +
 theme(legend.text = element_text(size=10)) +
 theme(legend.title = element_text(size=10)) +
 theme(legend.key.size = unit(0.4, 'cm')) + 
 theme(plot.margin = margin(1,1,1,1, "cm")) +
 scale_fill_manual(values=rainbow(11)) +
 ggtitle("quanTIseq")
ggsave("QUANTISEQ.png",width = 16,height = 8, units = "cm", dpi = 300)

Estimated immune cell fractions inferred using QUANTISEQ. ## 5. Gene Set Variation Analysis (GSVA) and single sample Gene Set Enrichment Analyses (ssGSEA)

Hallmark gene signatures were downloaded from teh Human Molecular Signatures Database (MSigDB) and complemented with specific immunological gene signatures and stored in GMT format (GENESETS.gmt). The GSVA package is used for sample wise GSVA and ssGSEA on gene expression files in matrix format (CRC_LOG2TPM1_MAT).The GSEABase package is used to read gmt format of signatures into a list object (GENE_SETS). For visualization purpose scores are z-transformed.

GMT <- getGmt("GENESETS.gmt")
GENE_SETS <- geneIds(GMT)
PARAM <- gsvaParam(exprData = as.matrix(CRC_LOG2TPM1_MAT), geneSets = GENE_SETS)
SCORES<- gsva(PARAM)
GSVA_SCORES<-as.data.frame(SCORES)
write.table(GSVA_SCORES,file="GSVA_CRC.txt",sep="\t",quote=FALSE, row.names=TRUE,col.names=NA)
GSVA_MAT <- GSVA_SCORES |>
  as.matrix()
storage.mode(GSVA_MAT) <- "numeric"
GSVA_Z <- t(scale(t(GSVA_MAT)))

The results from GSVA are stored in the dataframe CGSVA_SCORESand as text fileGSVA_CRC.txt`

GSVA_SCORES

Similarly, acquired and intrinsic MEK inhibitor resistance signatures were analyzed using sample-wise ssGSEA. MEK inhibitor resistance scores were calculated as:

MEKi_score_intrinsic = ssGSEA(IC50-positive genes) − ssGSEA(IC50-negative genes), MEKi_score_acquired = ssGSEA(acquired-upregulated genes) − ssGSEA(acquired-downregulated genes)

The resulting resistance scores were stored in the data frame MEKI_SCORE and exported as a text file.

MEKI<- getGmt("MEKi_selumetinib_CRC_signatures.gmt")
MEKI_SETS <- geneIds(MEKI)
PARAM_MEKI <- ssgseaParam(exprData = as.matrix(CRC_LOG2TPM1_MAT), geneSets = MEKI_SETS)
SCORES_MEKI<- gsva(PARAM_MEKI)
SSGSEA_SCORES<-as.data.frame(SCORES_MEKI)
MEKI_SCORE <- SSGSEA_SCORES[c(TRUE, FALSE), ] - SSGSEA_SCORES[c(FALSE, TRUE), ]
row.names(MEKI_SCORE)<- c("MEKi_Resist_Acqu","MEKi_Resist_Intr")
write.table(MEKI_SCORE,file="MEKI_RESISTANCE_SCORES.txt",sep="\t",quote=FALSE, row.names=TRUE,col.names=NA)

The results from ssGSEA for MEKi resistance signature are stored in the dataframe MEKI_SCORE and as text file MEKI_RESISTANCE_SCORES.txt

GSVA_SCORES

6. Heatmap indicating profile of enrichment scores of selected signatures

Enrichment profiles of selected gene signatures were calculated across all patients and visualized as z-scores. The heatmap included additional tumor-level annotations, including MEK inhibitor resistance scores, KRAS mutation status, and BRAF mutation status. Patients were grouped according to consensus molecular subtype (CMS), as defined by Guinney et al. (Nature Medicine, 2015), and hierarchical clustering was applied within the heatmap visualization using ComplexHeatmap.

meki_a <- as.numeric(MEKI_SCORE["MEKi_Resist_Acqu", ])
meki_i <- as.numeric(MEKI_SCORE["MEKi_Resist_Intr", ])
kras_status <- ifelse(CRC_CLIN$KRAS >= 1, "Mut", "WT")
braf_status <- ifelse(CRC_CLIN$BRAF >= 1, "Mut", "WT")
ha <- HeatmapAnnotation(
  CMS = CRC_CLIN$CMS,
  KRAS = kras_status,
  BRAF = braf_status,
  MEKi_A = meki_a,
  MEKi_I = meki_i,
  col = list(
    CMS = c(
      CMS1 = "orange",
      CMS2 = "blue",
      CMS3 = "magenta",
      CMS4 = "aquamarine4"
    ),
    KRAS = c(
      WT = "grey80",
      Mut = "black"
    ),
    BRAF = c(
      WT = "grey80",
      Mut = "darkblue"
    ),
    MEKi_A = circlize::colorRamp2(
      c(min(meki_a, na.rm = TRUE),
        median(meki_a, na.rm = TRUE),
        max(meki_a, na.rm = TRUE)),
      c("lightgrey", "white", "darkred")
    ),
    MEKi_I = circlize::colorRamp2(
      c(min(meki_i, na.rm = TRUE),
        median(meki_i, na.rm = TRUE),
        max(meki_i, na.rm = TRUE)),
      c("lightgrey", "white", "darkorange")
    )
  )
)
png("GSVA_HEATMAP.png",width = 4800,height = 2400,res = 300)
Heatmap(
  GSVA_Z,
  top_annotation = ha,
  column_split = CRC_CLIN$CMS,
  row_names_gp = gpar(fontsize = 10),
  row_names_max_width = unit(12, "cm"),
  show_column_names = FALSE,
  col = circlize::colorRamp2(
    c(-2, 0, 2),
    c("blue", "white", "red")
  )
)
dev.off()
## png 
##   2
Heatmap of hallmark gene signatures.
Heatmap of hallmark gene signatures.

7. Boxplots for signature scores or gene expression across different CMS

Define variable (gene signature) to show scores across consensus molecular subtypes in a boxplot.

VAR<-"IMMUNE_SIGNATURE" # define signature of interest       

signature_df <- tibble(
     PAT = colnames(GSVA_SCORES),
     !!VAR := as.numeric(GSVA_SCORES[VAR, ])
     )
CRC_CLIN_2 <- CRC_CLIN |>
      left_join(signature_df, by = "PAT")
CRC_CLIN_FILT <- CRC_CLIN_2 |>
      filter(!is.na(CMS))
df <- CRC_CLIN_FILT |>
  mutate(
    CMS = factor(CMS, levels = c("CMS1", "CMS2", "CMS3", "CMS4")),
    VAR = VAR 
  )

The final clinical information including the selected signatures is stored in CRC_CLIN_FILT

CRC_CLIN_FILT |>
  slice_head(n = 20)

Boxplots of gene expression or gene signature scores across consensus molecular subtypes (CMS) are generated using ggplot2. Individual patient values were overlaid using quasirandom point placement implemented in the ggbeeswarm package to visualize the underlying data distribution and sample density. Pairwise differences between CMS groups were assessed using the Wilcoxon rank-sum test, and significance levels were displayed on the plots.

plo <- ggplot(df, aes(x = CMS, y = .data[[VAR]])) +
  geom_boxplot(
    fill = "grey85",
    colour = "black",
    width = 0.7,
    outlier.shape = NA,
    alpha = 0.8,
    linewidth = 0.3 
  ) +
  ggbeeswarm::geom_quasirandom(
    aes(fill = CMS),
    shape = 21,
    colour = "black",
    stroke = 0.4,
    size = 2.5,
    width = 0.25,
    alpha = 0.9
  ) +
  ggpubr::stat_compare_means(
    comparisons = list(
      c("CMS1", "CMS2"),
      c("CMS1", "CMS3"),
      c("CMS1", "CMS4"),
      c("CMS2", "CMS3"),
      c("CMS2", "CMS4"),
      c("CMS3", "CMS4")
    ),
    method = "wilcox.test",
    label = "p.signif",
    size = 5,
    bracket.size = 0.4,
    step.increase = 0.1,
    vjust = 0.05
  ) +
  scale_fill_manual(values = c(
    CMS1 = "orange",
    CMS2 = "blue",
    CMS3 = "magenta",
    CMS4 = "aquamarine4"
  )) +
  labs(
    x = NULL,
    y = VAR
  ) +
  theme_classic(base_size = 16) +
  theme(
    axis.line = element_line(
      linewidth = 0.4,
      colour = "black"
    )
  ) +
  theme(
    legend.position = "none",
    axis.text.x = element_text(colour = "black", size = 16),
    axis.text.y = element_text(colour = "black")
  )
ggsave(
  paste0(VAR,"_CMS_BOXPLOT.png"),
  plo,
  width = 6,
  height = 5,
  dpi = 300
)
Boxplot of signature score across CMS
Boxplot of signature score across CMS

8. Survival analysis and Kaplan-Meier plots

Patient are dichotomized on median values from the selected signature (gene) and overal survival are visualized as Kaplan-Meier curves using packages survival and survminer.

CRC_CLIN_OS <- CRC_CLIN_2 |>
  mutate(
    OS_TIME = coalesce(DAYS_TO_DEATH, DAYS_TO_LAST_FOLLOWUP)* 12 / 365.25,
    OS_EVENT = ifelse(VITAL_STATUS == "1", 1, 0),
    VAR_GROUP = ifelse(
      .data[[VAR]] >= median(.data[[VAR]], na.rm = TRUE),
      "High",
      "Low"
    )
  ) |>
  filter(
    !is.na(OS_TIME),
    !is.na(OS_EVENT),
    !is.na(.data[[VAR]])
  )
CRC_CLIN_OS <- CRC_CLIN_OS |>
  mutate(
    VAR_GROUP = factor(VAR_GROUP, levels = c("Low", "High"))
  )

fit <- survfit(Surv(OS_TIME, OS_EVENT) ~ VAR_GROUP,data = CRC_CLIN_OS)
cox <- coxph(Surv(OS_TIME, OS_EVENT) ~ VAR_GROUP,data = CRC_CLIN_OS)
s <- summary(cox)
HR <- round(s$conf.int[1, "exp(coef)"], 2)
LCL <- round(s$conf.int[1, "lower .95"], 2)
UCL <- round(s$conf.int[1, "upper .95"], 2)
PVAL <- signif(s$coefficients[1, "Pr(>|z|)"], 3)
HR_LABEL <- paste0("HR = ", HR," (95% CI ", LCL, "-", UCL, ")\n","P = ", PVAL)

pos <- ggsurvplot(
  fit,
  data = CRC_CLIN_OS,
  pval = FALSE,
  risk.table = "absolute",
  break.time.by = 25,
  risk.table.height = 0.25,
  risk.table.y.text = FALSE,
  conf.int = TRUE,
  palette = c("blue", "red"),
  legend.title = "",
  legend.labs = c("Low","High"),
  xlab = "Months",
  ylab = "Overall survival probability",
)
pos$plot <- pos$plot +
  ggtitle(VAR) +
  theme(
    plot.title = element_text(
      face = "bold",
      hjust = 0.5,
      size = 16
    )
  )
pos$plot <- pos$plot +
  annotate(
    "text",
    x = 0,
    y = 0.1,
    label = HR_LABEL,
    hjust = 0,
    size = 5
  )
png(paste0(VAR, "_KM.png"),  width = 2100,height = 2100,res = 300)
print(pos, newpage = FALSE)
dev.off()
## png 
##   2
Kaplan-Meier curves
Kaplan-Meier curves