Two complementary approaches were used to derive MEK inhibitor resistance gene signatures. First, to identify genes associated with intrinsic resistance, baseline expression data from colorectal cancer cell lines in the Cancer Cell Line Encyclopedia (CCLE) were correlated with Selumetinib IC50 values using linear regression. Second, to identify genes associated with acquired resistance, differential expression analysis was performed using Gene Expression Omnibus dataset GSE126109, which includes five colorectal cancer cell lines and their matched Selumetinib-resistant derivatives generated by long-term drug exposure.
Include all R and Bioconductor packages
library(conflicted)
library(kableExtra)
library(readr)
library(AnnotationDbi)
library(org.Hs.eg.db)
library(tidyverse)
library(writexl)
library(DESeq2)
conflicts_prefer(dplyr::filter)
conflicts_prefer(dplyr::select)
conflicts_prefer(dplyr::slice_head)Set your preferred working directory and ensure that the required
input files, CCLE_TPM.txt, GDSC1.txt,
GSE126109_SAMPLES.txt, and
GSE126109_COUNT_MATRIX.txt, are stored in this
directory.
setwd("L:/2026/WM7/R") # use your preferred directoryFrom the Cancer Cell Line Encyclopedia (CCLE) downloaded a file including gene expression values (TPM) for > 1000 cancer cell lines https://depmap.org/portal/data_page/?tab=allData&releasename=CCLE%202019&filename=CCLE_RNAseq_rsem_genes_tpm_20180929.txt.gz the final file was renamed to
CCLE_TPM.txtWe can use the read_tsv() function from the
readr package to read the file and assign it to
ccle_expr.
CCLE_EXPR <- read_tsv("CCLE_TPM.txt")
CCLE_EXPR |>
select(1:8) |>
slice_head(n = 20) Transfrom data: Select only colon cancer cell lines (LARGE_INTESTINE), log2-transform TPM adding a pseudocount of 1, adding gene symbol column, filter for genes with no expression in more than 20 cancer cell lines
COL_ind <- grep("LARGE_INTESTINE", names(CCLE_EXPR))
CCLE <- CCLE_EXPR |>
select(all_of(COL_ind)) |>
mutate(across(everything(), \(x) log2(x + 1))) |>
rename_with(\(x) gsub("_LARGE_INTESTINE", "", x))
CCLE_ANNO <- CCLE |>
mutate(
ENSEMBL = gsub("\\..*", "", CCLE_EXPR$gene_id),
.before = 1
) |>
left_join(
tibble(ENSEMBL = unique(gsub("\\..*", "", CCLE_EXPR$gene_id))) |>
mutate(
SYMBOL = AnnotationDbi::mapIds(
org.Hs.eg.db,
keys = ENSEMBL,
column = "SYMBOL",
keytype = "ENSEMBL",
multiVals = "first"
)
),
by = "ENSEMBL"
) |>
relocate(SYMBOL, .after = ENSEMBL)
CCLE_ANNO_FILT <- CCLE_ANNO |>
mutate(n_low = rowSums(across(where(is.numeric), \(x) x < 0.001))) |>
filter(!is.na(SYMBOL),n_low < 20) |>
select(-n_low)The final CCLE table will look like:
CCLE_ANNO_FILT |>
select(ENSEMBL, SYMBOL, 3:8) |>
slice_head(n = 20)The final CCLE table will be formatted and saved as text file:
CCLEP <- CCLE_ANNO_FILT |>
mutate(across(where(is.numeric),\(x) sprintf("%.4f", x)))
write_tsv(CCLEP,"CCLE_LOG2TPM1_CRC_FILT.txt") The Genomics of Drug Sensitivity in Cancer (GDSC) database was used
to obtain drug sensitivity data for cancer cell lines, which were saved
in the file GDSC1.txt. For this analysis, the primary
information used was the IC50 values (or log-transformed IC50; ln IC50)
for the MEK inhibitor Selumetinib in colorectal cancer cell lines. The
GDSC1 (legacy) release was used and remains available through the GDSC
website or can be downloaded directly as an Excel file from GDSC1_fitted_dose_response.
GDSC<-read_tsv("GDSC1.txt")
GDSC<- GDSC |>
mutate(CELL_LINE_NAME = toupper(gsub("-", "", CELL_LINE_NAME)))|>
filter(TCGA_DESC=="COREAD",DRUG_NAME=="Selumetinib", DRUG_ID=="1062")
write_tsv(GDSC,"GDSC_CRC_SELUMETINIB_IC50.txt")
GDSC1<-GDSC |>
column_to_rownames("CELL_LINE_NAME")The final IC50 table from GDSC GDSC1
GDSC1CL <- base::intersect(GDSC$CELL_LINE_NAME, names(CCLE_ANNO_FILT))
IC50<-GDSC1[CL,"LN_IC50"]
CCLEX <- CCLE_ANNO_FILT |>
select(all_of(c("SYMBOL",CL)))
CCLEXU <- CCLEX |>
distinct(SYMBOL, .keep_all = TRUE)This section defines the plotting function
plotlinreg(gene_name, IC50, expression), which generates a
linear regression plot for a given gene. The Selumetinib ln(IC50) value
is shown on the x-axis, and gene expression, expressed as log2(TPM + 1),
is shown on the y-axis. The plot includes the fitted regression line
with confidence intervals, together with Pearson’s correlation
coefficient and the unadjusted p-value.
plotlinreg<-function (name,y,x) {
plot_df <- tibble(
IC50 = y,
expr = x
)
di<-max(x)-min(x)
fit1<-lm(IC50 ~ x,data=plot_df)
s1 <- summary(fit1)
p1 <- s1$coefficients[2,4]
r21<-r2 <- s1$r.squared
r1 <- sign(coef(fit1)[2]) * sqrt(r21)
g <- ggplot(plot_df, aes(IC50, expr)) +
geom_point() +
geom_smooth(method = "lm") +
labs(
x = "Selumetinib ln(IC50)",
y = paste0(name, " log2(TPM+1)")
) +
theme_classic(base_size = 16) +
theme(
plot.margin = margin(1.5,1.5,1.5,1.5, "cm"),
axis.title.y = element_text(vjust = +4),
axis.title.x = element_text(vjust = -2),
axis.text.y=element_text(size=16),
axis.text.x=element_text(size=16),
) +
annotate("text", x=0.75*max(IC50), y=max(x)+0.2*di, label= paste0("R = ",round(r1, 3)), size=5) +
annotate("text", x=0.75*max(IC50), y=max(x)+0.1*di, label= paste0("p = ",signif(p1, 3)), size=5)
ggsave(
paste0(name, "_REG_PLOT.pdf"),
g,
width = 14,
height = 12,
units = "cm"
)
}The expression object CCLEXU is split by gene symbol and
processed using map_dfr() from the purrr package. For each
gene, a linear regression model is fitted between gene expression and
Selumetinib ln(IC50) values. Regression statistics, including slope
(coefficient), Pearson’s correlation coefficient (R), coefficient of
determination (R²), and p-value, were extracted and combined into a
single results table RESULTS. Only genes meeting the
predefined significance criteria (P < 0.01 and R² > 0.1) are
retained. For genes with p < 0.002 a regresion plot is saved in pdf
format using the defined plotting function
plotlinreg(gene_name, IC50, expression).
RESULTS <- CCLEXU |>
split(CCLEXU$SYMBOL) |>
map_dfr(\(gene) {
x <- unlist(gene[-1])
fit <- lm(IC50 ~ x)
s <- summary(fit)
p <- s$coefficients[2,4]
r2 <- s$r.squared
r <- sign(coef(fit)[2]) * sqrt(r2)
if (p < 0.002) {
plotlinreg(gene$SYMBOL[1],IC50,x)
}
if (p < 0.01 && r2 > 0.1) {
tibble(
symbol = gene$SYMBOL[1],
slope = coef(fit)[2],
R = r,
R2 = r2,
P = p
)
} else {
tibble()
}
})
RESULTS<-RESULTS |>
arrange(desc(R))The final RESULTS table:
RESULTSWrite results to text file and Excel file
write_tsv(RESULTS,"SELUMETINIB_CRC_RESISTANT_SENSITIVE_GENES.txt")
write_xlsx(RESULTS,"SELUMETINIB_CRC_RESISTANT_SENSITIVE_GENES.xlsx") Raw sequencing data for GSE126109 were downloaded from the Gene Expression Omnibus (GEO) as the archive GSE126109_RAW.tar (https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE126109&format=file). The following samples were included in the analysis: HCT116R, HCT116R_6244R, HKH2, HKH2_6244R, LoVo, LoVo_6244R, SW480, SW480_6244R, SW620, and SW620_6244R, representing five colorectal cancer cell lines and their corresponding Selumetinib (6244R)-resistant derivatives. For each sample, gene-level count tables were extracted from the archive and combined into a single count matrix (GSE126109_COUNT_MATRIX.txt). Sample annotations, including cell line and treatment status, were compiled into a corresponding sample information table (GSE126109_SAMPLES.txt).
SAMPLES<- read_tsv("GSE126109_SAMPLES.txt")
COUNTS<- read_tsv("GSE126109_COUNT_MATRIX.txt")The sample information
SAMPLESThe final counts table
COUNTS |>
select(1:8) |>
slice_head(n = 20)As differences between cell lines are expected to be greater than
those associated with acquired resistance, as indicated by the PCA
analysis, differential expression analysis was performed using a paired
design in DESeq2 design = ~ CELL + COND).
count_mat <- COUNTS |> column_to_rownames("gene") |> round()
sample_info <- SAMPLES |> select(-FILE)
dds <- DESeqDataSetFromMatrix(countData = count_mat,colData = sample_info,design= ~ CELL + COND)
dds$COND = relevel(dds$COND, "CTRL")
dds <- DESeq(dds)
res <- results(dds)
vsd <- vst(dds, blind=FALSE)
res_significant <- as.data.frame(res) |>
filter(!is.na(padj), padj < 0.1, baseMean>10) |>
arrange(desc(log2FoldChange))
res_sig_tbl <- res_significant |> rownames_to_column("gene") |> tibble::as_tibble()The final significant results are stored in
res_sig_tbl
res_sig_tbl |>
select(gene, baseMean, log2FoldChange, padj) |>
slice_head(n = 50)Write results to a text file and an Excel file
write_tsv(res_sig_tbl,"SELUMETINIB_CRC_ACQUIRED_RESISTANCE_GENES.txt")
write_xlsx(res_sig_tbl,"SELUMETINIB_CRC_ACQUIRED_RESISTANCE_GENES.xlsx")Generate a list gene signatures for both acquired resistance (sensitivity) and intrinsic resistance (sensitivity) and save in GMT format as used for Gene Set Enrichment Analysis (GSEA).
L<-list()
L$acquired_up_genes <- pull(res_sig_tbl |> filter(log2FoldChange > 0.3) |> select(gene), gene)
L$acquired_down_genes <- pull(res_sig_tbl |> filter(log2FoldChange < -0.3) |> select(gene),gene)
L$IC50_positive_genes <- as.vector(pull(RESULTS |> filter(R > 0.3) |> select(symbol),symbol))
L$IC50_negative_genes <- as.vector(pull(RESULTS |> filter(R < -0.3) |> select(symbol),symbol))
write_gmt <- function(gene_sets, file) {
con <- file(file, open = "wt")
on.exit(close(con))
for (set_name in names(gene_sets)) {
genes <- unique(na.omit(gene_sets[[set_name]]))
n_genes <- length(genes)
line <- paste(c(set_name, n_genes, genes),collapse = "\t")
writeLines(line, con)
}
}
write_gmt(L, "MEKi_selumetinib_CRC_signatures.gmt")To assess the likelihood of a patient’s response to therapy, we defined two Selumetinib (MEK inhibitor, MEKi) resistance scores representing intrinsic and acquired resistance. Both scores were calculated using single-sample gene set enrichment analysis (ssGSEA) implemented in the GSVA package on bulk RNA sequencing data from the tumor. The intrinsic resistance score is defined as:
MEKi_score_intrinsic = ssGSEA(IC50-positive genes) − ssGSEA(IC50-negative genes), whereas the acquired resistance score is defined as:
MEKi_score_acquired = ssGSEA(acquired-upregulated genes) − ssGSEA(acquired-downregulated genes).
Higher scores indicate a gene expression profile associated with increased resistance to Selumetinib and, consequently, a lower likelihood of therapeutic response.