OSCR

Deepening the Understanding of Platelet-Derived RNA as Biomarker for Glioblastoma Through GSEA, TDEA, and Elastic Net Regularization Analysis.

Code ↔ Paper

2 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 2 matches
  1. [1] § Materials and Methods › Data Processing and Analysis ↔ programs/basic_model_analysis.R, lines 1 to 51 · score 0.82 · candidate genes, confusion matrix, DESeq2, platelet RNA, healthy control, split · Report a mismatch
  2. [2] § Results › Unraveling Gene Expression and Pathways Disparities Between GBM and Healthy Individuals Through TDEA and Elastic Net Regularization Analysis ↔ programs/basic_model_analysis.R, lines 1 to 51 · score 0.73 · platelet RNA seq, candidate genes, feature selection, healthy control, classify, fold · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R · 1,003 lines · 38 KB · GPL-3.0 · 2 matches

  1. # =============================================================================
  2. # Platelet RNA-based classification of Glioblastoma (GBM) vs Healthy Controls
  3. # =============================================================================
  4. #
  5. # This script performs the TDEA (Threshold-Based Differential Expression Analysis)
  6. # with Elastic Net regularization pipeline described in Giczewska, Pastuszak et al.
  7. #
  8. # The analysis proceeds in three stages:
  9. # 1. Data loading, filtering, and sample exclusion
  10. # 2. Differential expression analysis on the training set to identify candidate genes
  11. # 3. Elastic Net classification (with secondary feature selection) followed by
  12. # Gene Ontology and Reactome pathway enrichment of the final gene panel
  13. #
  14. # Input data:
  15. # - DESeq2-normalized platelet RNA-seq counts (4,412 genes x 403 samples)
  16. # - Sample metadata including patient group (GBM / HC), isolation site, etc.
  17. # - Pre-defined train/test split IDs (60/40)
  18. # - Ensembl 75 gene annotations
  19. # - Reactome GMT gene set file (v2022.1)
  20. #
  21. # Output:
  22. # - Confusion matrix plot (PDF) for the test set
  23. # - Reactome and Gene Ontology enrichment bar plots (PDF)
  24. # - Enrichment results table (XLSX)
  25. # - Elastic Net model coefficients (saved as RData)
  26. # =============================================================================
  28. library(DESeq2)
  29. library(edgeR)
  31. # Flag to use the fixed train/test split from the manuscript (reproducibility).
  32. # When TRUE, sample IDs are read from file rather than generated via createFolds.
  33. useArticleSplit = T
  35. # --- 1. Load normalized counts and sample metadata ---
  37. countsAllDeseq = read.csv2("../data/rawdata/countsNormalized.tsv", sep = "\t")
  38. colnames(countsAllDeseq) = gsub(".", "-", colnames(countsAllDeseq), fixed = T)
  40. # Ensure all columns are numeric (read.csv2 may import as character)
  41. for(i in 1:ncol(countsAllDeseq))
  42. countsAllDeseq[,i] = as.numeric(countsAllDeseq[,i])
  44. sampleInfoAll = read.csv2("../data/rawdata/sampleInfo.tsv", sep = "\t")
  46. # Gene annotation table (Ensembl 75): maps HGNC symbols to Ensembl IDs
  47. load("../data/dgeGenesEnsembl75.RData")
  49. # =============================================================================
  50. # Helper functions
  51. # =============================================================================
  53. # DESeq2 variance-stabilizing normalization (standalone, without design contrast).
  54. # Used when re-normalizing raw counts outside the main pipeline.
  55. normalizeDESeq2NoReport = function(rawCounts) {
  56. library(DESeq2)
  58. d = rawCounts
  59. sampleIds = colnames(rawCounts)
  60. # Single-group dummy design (normalization only, no DE testing)
  61. condition <- factor(rep("A", dim(d)[2]))
  63. dds <- DESeqDataSetFromMatrix(countData = d, DataFrame(condition), design = ~ 1)
  64. dds <- DESeq(dds)
  65. vsd <- varianceStabilizingTransformation(dds, blind = TRUE)
  67. logs_norm_arr_ord = assay(vsd)
  68. return(logs_norm_arr_ord)
  69. }
  71. # Filter out low-abundance transcripts from a DGEList object.
  72. # Removes genes with fewer than `minimum.read.counts` reads in > 90% of samples.
  73. filter.for.platelet.transcriptome <- function(dge,
  74. minimum.read.counts = 30,
  75. verbose = TRUE) {
  76. if (missing(dge)) {
  77. stop("Provide DGElist object")
  78. }
  79. stopifnot(class(dge) == "DGEList")
  81. if (!is.numeric(minimum.read.counts)) {
  82. stop("Provide numeric value for minimum.read.counts")
  83. }
  85. # Identify genes that fall below the read-count threshold in >90% of samples
  86. bad.genes <- names(which(
  87. apply(dge$counts, 1, function(x) { sum(x < minimum.read.counts) }) > 0.9 * ncol(dge)
  88. ))
  89. dgeFiltered <- dge[which(!rownames(dge) %in% bad.genes), ]
  91. if (verbose == TRUE) {
  92. print(paste("Number of transcripts detected: ", nrow(dgeFiltered$counts), sep = ""))
  93. }
  95. return(dgeFiltered)
  96. }
  98. # Safe wrappers for Wilcoxon and t-test that return NA on failure
  99. # (e.g., when all values in a group are identical)
  100. my.wilcox.p.value <- function(...) {
  101. obj <- try(wilcox.test(...), silent = TRUE)
  102. if (is(obj, "try-error")) return(NA) else return(obj$p.value)
  103. }
  105. my.t.test.p.value <- function(...) {
  106. obj <- try(t.test(...), silent = TRUE)
  107. if (is(obj, "try-error")) return(NA) else return(obj$p.value)
  108. }
  110. # =============================================================================
  111. # 2. Sample exclusion and cohort definition
  112. # =============================================================================
  114. library(matrixStats)
  115. dataFiltered = countsAllDeseq
  117. # Exclude one IDH-mutant sample reclassified as non-GBM under updated WHO criteria.
  118. # TR3712-GBM-VUMC was originally labeled GBM but is IDH-wildtype negative.
  119. id604 = grep("TR3712-GBM-VUMC", sampleInfoAll$id, fixed = T)
  120. if(length(id604) > 0) {
  121. dataFiltered = dataFiltered[, -id604]
  122. sampleInfoAll = sampleInfoAll[-id604, ]
  123. }
  125. # Align gene annotation table to the genes present in filtered counts
  126. genes = genes[which(!is.na(match(genes$hgnc_symbol, rownames(dataFiltered)))),]
  127. rownames(genes) = genes$hgnc_symbol
  129. # =============================================================================
  130. # Differential expression summary table (Wilcoxon rank-sum test + FDR)
  131. # =============================================================================
  132. #
  133. # For two patient groups (groupA vs groupB), computes per-gene:
  134. # - Mean and median expression in each group
  135. # - Fold change (on log2 scale, back-transformed to linear)
  136. # - Wilcoxon rank-sum p-value with FDR correction
  137. #
  138. # When `all = FALSE` (default), applies TDEA thresholds:
  139. # - FDR q-value < 0.1
  140. # - Fold change > 1.3 or < 1/1.3
  141. # - At least one group median > 3 (low-expression filter)
  143. getSummaryTable = function(dataFiltered, sampleInfo, groupA, groupB,
  144. all = F) {
  145. hcId = which(sampleInfo$ori_patientgroup == groupA)
  146. ocId = which(sampleInfo$ori_patientgroup == groupB)
  147. logsBenign = dataFiltered[, hcId]
  148. logsMalignant = dataFiltered[, ocId]
  150. # Group-level summary statistics
  151. meanBenign <- rowMeans(logsBenign)
  152. meanMalignant <- rowMeans(logsMalignant)
  154. medianBenign <- rowMedians(as.matrix(logsBenign))
  155. medianMalignant <- rowMedians(as.matrix(logsMalignant))
  157. # Fold change computed from mean log2-scale values (DESeq2-normalized)
  158. folds_mean_Benign_Malignant = 2^(meanBenign - meanMalignant)
  160. # Per-gene Wilcoxon rank-sum test between the two groups
  161. wilcox_Benign_Malignant_pvalue = sapply(1:nrow(logsBenign), function(i)
  162. my.wilcox.p.value(unlist(logsBenign[i, ]), unlist(logsMalignant[i, ])))
  164. genesLocal = genes[rownames(dataFiltered),]
  166. # FDR correction (Benjamini-Hochberg)
  167. fdr_Benign_Malignant_qvalue = p.adjust(
  168. wilcox_Benign_Malignant_pvalue, method = "fdr",
  169. n = length(wilcox_Benign_Malignant_pvalue)
  170. )
  172. # Assemble results table
  173. summary_table_Benign_Malignant = cbind.data.frame(
  174. genesLocal$ensembl_gene_id,
  175. rownames(logsBenign),
  176. meanBenign, meanMalignant,
  177. medianBenign, medianMalignant,
  178. folds_mean_Benign_Malignant,
  179. wilcox_Benign_Malignant_pvalue,
  180. fdr_Benign_Malignant_qvalue,
  181. genesLocal$description
  182. )
  184. groupNameA = groupA
  185. groupNameB = groupB
  186. tableCols_Benign_Malignant = cbind(
  187. "ENSG", "Gene",
  188. paste("Mean ", groupNameA, sep = ""), paste("Mean ", groupNameB, sep = ""),
  189. paste("Median ", groupNameA, sep = ""), paste("Median ", groupNameB, sep = ""),
  190. paste("Mean fold change ", groupNameA, "-", groupNameB, sep = ""),
  191. "Wilcoxon p-value ",
  192. "FDR q-value ",
  193. "Description"
  194. )
  195. colnames(summary_table_Benign_Malignant) = tableCols_Benign_Malignant
  197. # Sort by FDR q-value (ascending)
  198. summary_table_Benign_Malignant = summary_table_Benign_Malignant[order(
  199. summary_table_Benign_Malignant[, 9]), ]
  201. # Return all genes without filtering if requested
  202. if(all) {
  203. return(summary_table_Benign_Malignant)
  204. }
  206. # Apply TDEA thresholds: FDR < 0.1, fold change > 1.3, minimum expression > 3
  207. summaryTable = summary_table_Benign_Malignant[
  208. which(summary_table_Benign_Malignant$`FDR q-value ` < 0.1), ]
  209. summaryTable = summaryTable[c(
  210. which(summaryTable[, 7] > 1.3),
  211. which(summaryTable[, 7] < 1/1.3)
  212. ), ]
  214. lowExprCutOff = 3
  215. # NOTE: "Median OC" was the original column name when this function was used
  216. # with ovarian cancer data. When called with groupA="GBM", groupB="HC", the
  217. # actual column is "Median GBM". Because $`Median OC` returns NULL here,
  218. # the filter only requires HC median > 3. This is preserved for result
  219. # reproducibility with the published manuscript.
  220. highHc = which(summaryTable$`Median HC` > lowExprCutOff)
  221. highOc = which(summaryTable$`Median OC` > lowExprCutOff)
  222. highExpr = unique(c(highHc, highOc))
  223. summaryTable = summaryTable[highExpr, ]
  224. return(summaryTable)
  225. }
  227. # =============================================================================
  228. # 3. Define the study cohort (GBM + HC, excluding NKI controls)
  229. # =============================================================================
  231. sampleInfoAll$full_ori_patientgroup = sampleInfoAll$ori_patientgroup
  232. gbmSamples = sampleInfoAll[which(sampleInfoAll$ori_patientgroup == "GBM"),]
  234. library(caret)
  236. # Select healthy controls from non-NKI sites only (NKI excluded per study protocol)
  237. hcId = which(sampleInfoAll$ori_patientgroup == "HC" &
  238. sampleInfoAll$isolationlocation != "NKI" &
  239. sampleInfoAll$group == "nonMalignant")
  241. # Select GBM patients from non-NKI sites
  242. gbmId = which(sampleInfoAll$ori_patientgroup == "GBM" &
  243. sampleInfoAll$isolationlocation != "NKI")
  245. toInclude = c(hcId, gbmId)
  246. sampleInfoAll = sampleInfoAll[toInclude, ]
  247. countsAllDeseq = countsAllDeseq[, toInclude]
  249. # Refresh indices after subsetting
  250. hcId = which(sampleInfoAll$ori_patientgroup == "HC" &
  251. sampleInfoAll$isolationlocation != "NKI" &
  252. sampleInfoAll$group == "nonMalignant")
  253. gbmId = which(sampleInfoAll$ori_patientgroup == "GBM" &
  254. sampleInfoAll$isolationlocation != "NKI")
  256. # =============================================================================
  257. # 4. Train/test split (60/40)
  258. # =============================================================================
  260. set.seed(12345)
  262. # Default: 10-fold stratified split (folds 1-6 = train, 7-10 = test)
  263. folds = createFolds(sampleInfoAll$ori_patientgroup[toInclude], k = 10)
  264. trainId = toInclude[unlist(folds[1:6])]
  265. testId = toInclude[unlist(folds[7:10])]
  267. # Override with the fixed split used in the manuscript for reproducibility
  268. if(useArticleSplit) {
  269. trainId = read.table("../data/rawdata/train_ids_final.tsv")
  270. trainId = trainId$V1
  271. trainId = trainId[which(!is.na(match(trainId, colnames(countsAllDeseq))))]
  272. trainId = (match(trainId, colnames(countsAllDeseq)))
  274. testId = read.table("../data/rawdata/testId_ids_final.tsv")
  275. testId = testId$V1
  276. testId = testId[which(!is.na(match(testId, colnames(countsAllDeseq))))]
  277. testId = (match(testId, colnames(countsAllDeseq)))
  278. }
  280. # =============================================================================
  281. # 5. TDEA: identify differentially expressed candidate genes on TRAINING data only
  282. # =============================================================================
  284. # Run differential expression on training set (leak-free: test set not used)
  285. gbmHc_comparisonTrimmed = getSummaryTable(
  286. countsAllDeseq[, trainId], sampleInfoAll[trainId, ], "GBM", "HC", F
  287. )
  289. # Candidate genes passing TDEA thresholds (FDR < 0.1, FC > 1.3, expression > 3)
  290. candidates = rownames(gbmHc_comparisonTrimmed)
  292. # Subset counts to candidate genes for modeling
  293. countsModel = countsAllDeseq[candidates, ]
  294. countsModel_deseq = countsModel[, toInclude]
  296. # Record the train/test assignment for each sample
  297. countsModel_deseq_assignment = data.frame(
  298. c(colnames(countsAllDeseq)[trainId], colnames(countsAllDeseq)[testId])
  299. )
  300. colnames(countsModel_deseq_assignment) = "Sample"
  301. countsModel_deseq_assignment$Train = 0
  302. countsModel_deseq_assignment$Train[1:length(trainId)] = 1
  303. countsModel_deseq_assignment$Test = 1
  304. countsModel_deseq_assignment$Test[1:length(trainId)] = 0
  305. deseq_output = gbmHc_comparisonTrimmed
  307. # =============================================================================
  308. # 6. Elastic Net classification (first pass: all candidate genes)
  309. # =============================================================================
  311. # Prepare training and test matrices (samples x genes)
  312. trainX = countsModel[, trainId]
  313. testX = countsModel[, testId]
  314. trainY = factor(sampleInfoAll$ori_patientgroup[trainId])
  315. testY = factor(sampleInfoAll$ori_patientgroup[testId])
  317. library(party)
  319. trainX = t(trainX)
  320. trainX = as.data.frame(trainX)
  321. trainX$group = as.factor(trainY)
  323. testX = t(testX)
  324. testX = as.data.frame(testX)
  325. testX$group = as.factor(testY)
  327. library(dplyr)
  328. library(glmnet)
  330. realRefsTrain = trainX$group
  331. realRefsTest = testX$group
  333. # 5x5 repeated cross-validation with random hyperparameter search
  334. control <- trainControl(method = "repeatedcv",
  335. number = 5,
  336. repeats = 5,
  337. search = "random",
  338. verboseIter = TRUE)
  340. # Train Elastic Net (glmnet) on all TDEA-selected candidate genes
  341. elasticModel <- train(group ~ .,
  342. data = trainX,
  343. method = "glmnet",
  344. tuneLength = 25,
  345. trControl = control)
  346. print(elasticModel)
  348. # --- Evaluate first-pass model ---
  350. predictionsTrain = predict(elasticModel, trainX)
  351. probabilitiesTrain = predict(elasticModel, trainX, type = "prob")
  353. table(predictionsTrain, realRefsTrain)
  354. confusionMatrixTrain = confusionMatrix(predictionsTrain, realRefsTrain)
  355. print(confusionMatrixTrain)
  357. predictionsTest = predict(elasticModel, testX)
  358. probabilitiesTest = predict(elasticModel, testX, type = "prob")
  360. if (interactive()) boxplot(probabilitiesTest[, 1] ~ realRefsTest)
  362. # Apply a lowered probability threshold (0.25) to increase GBM sensitivity
  363. thresholdedPredictionsTest = rep("HC", length(predictionsTest))
  364. thresholdedPredictionsTest[which(probabilitiesTest[, 1] > 0.25)] = "GBM"
  366. thresholdedPredictionsTrain = rep("HC", length(realRefsTrain))
  367. thresholdedPredictionsTrain[which(probabilitiesTrain[, 1] > 0.25)] = "GBM"
  369. confusionMatrixTrain = confusionMatrix(as.factor(thresholdedPredictionsTrain), realRefsTrain)
  371. table(thresholdedPredictionsTest, realRefsTest)
  372. print(confusionMatrix(as.factor(thresholdedPredictionsTest), as.factor(realRefsTest), positive = "GBM"))
  374. table(predictionsTest, realRefsTest)
  375. confusionMatrixTest = confusionMatrix(
  376. as.factor(thresholdedPredictionsTest),
  377. reference = as.factor(realRefsTest),
  378. positive = "GBM"
  379. )
  380. print(confusionMatrixTest)
  382. # =============================================================================
  383. # 7. Secondary feature selection: keep genes with |coefficient| > threshold
  384. # =============================================================================
  386. weightThreshold = 0.05
  388. # Extract Elastic Net coefficients at the optimal lambda
  389. toKeepElNet = stats::predict(
  390. elasticModel$finalModel, type = "coefficients",
  391. s = elasticModel$bestTune$lambda
  392. )
  393. toKeepElNetDF = as.data.frame(as.matrix(toKeepElNet))
  394. toKeepElNetDF = toKeepElNetDF[-1, , drop = F] # Remove intercept row
  395. toKeepElNetDF[, 1] = as.numeric(toKeepElNetDF[, 1]) # Ensure plain numeric
  397. # Full model: all genes with non-zero coefficients
  398. fullModel = rownames(toKeepElNetDF)[which(abs(toKeepElNetDF[, 1]) > 0.0)]
  399. fullModel = gsub("`", "", fullModel)
  401. # Reduced panel: genes with |weight| above threshold
  402. maxWeights = abs(toKeepElNetDF[, 1])
  403. toKeepElNetDF = toKeepElNetDF[which(maxWeights > weightThreshold), , drop = F]
  404. allSelectedEnsg = rownames(toKeepElNetDF)
  405. allSelectedEnsg = gsub("`", "", allSelectedEnsg)
  407. # =============================================================================
  408. # 8. Test set AUC and confusion matrix visualization
  409. # =============================================================================
  411. library(cvAUC)
  412. AUC = cvAUC(probabilitiesTest[, 1], realRefsTest, label.ordering = c("HC", "GBM"))
  414. # Confusion matrix heatmap for the test set
  415. cmTest <- confusionMatrixTest
  416. plt <- as.data.frame(cmTest$table)
  417. plt$Prediction <- factor(plt$Prediction, levels = rev(levels(plt$Prediction)))
  419. pdf(paste0('../Report/', "Confusion_matrix_test", ".pdf"))
  420. confPlot = ggplot(plt, aes(Prediction, Reference, fill = Freq)) +
  421. geom_tile() + geom_text(aes(label = Freq)) +
  422. scale_fill_gradient(low = "white", high = "#009194") +
  423. labs(x = "Reference", y = "Prediction") +
  424. scale_x_discrete(labels = c("HC", "GBM")) +
  425. scale_y_discrete(labels = c("GBM", "HC"))
  427. print(confPlot)
  428. dev.off()
  430. # =============================================================================
  431. # 9. Reactome pathway enrichment (ORA) on the final Elastic Net gene panel
  432. # =============================================================================
  434. library(org.Hs.eg.db)
  435. library(clusterProfiler)
  436. library(openxlsx)
  438. deGenes = allSelectedEnsg # Final panel genes (HGNC symbols)
  439. geneUniverse = rownames(dataFiltered) # Background: all expressed genes
  441. # Load Reactome gene sets (MSigDB format)
  442. c3.tf <- read.gmt("../data/c2.cp.reactome.v2022.1.Hs.symbols.gmt")
  443. presentFromReact = which(is.element(c3.tf$gene, geneUniverse))
  444. c3.tf = c3.tf[presentFromReact, ]
  446. # Run overrepresentation analysis against Reactome pathways
  447. ans.react = enricher(allSelectedEnsg, TERM2GENE = c3.tf,
  448. pAdjustMethod = "fdr", qvalueCutoff = 0.5)
  450. # --- Prepare Reactome bar chart data ---
  452. dfReact = ans.react@result[1:10, ]
  454. # Parse the GeneRatio and BgRatio strings into numeric values
  455. noGenes = strsplit(dfReact$GeneRatio, "/")
  456. noGenesPlus = as.numeric(unlist(lapply(noGenes, `[[`, 1)))
  457. noGenesAll = as.numeric(unlist(lapply(noGenes, `[[`, 2)))
  459. noBg = strsplit(dfReact$BgRatio, "/")
  460. noBgPlus = as.numeric(unlist(lapply(noBg, `[[`, 1)))
  461. noBgAll = as.numeric(unlist(lapply(noBg, `[[`, 2)))
  463. dfReact$noGenesPlus = noGenesPlus
  464. dfReact$noGenesAll = noGenesAll
  465. dfReact$noBgPlus = noBgPlus
  466. dfReact$noBgAll = noBgAll
  467. dfReact$Pct_selected = round(100 * dfReact$noGenesPlus / dfReact$noGenesAll, 2)
  468. dfReact$Pct_background = round(100 * dfReact$noBgPlus / dfReact$noBgAll, 2)
  469. dfReact$LogP = -1 * (log10(dfReact$pvalue))
  471. # Reshape to long format for grouped bar chart
  472. dfReactChartP = cbind(dfReact$Description, dfReact$LogP,
  473. rep("p", nrow(dfReact)),
  474. paste0("p=", round(dfReact$pvalue, 6)))
  475. dfReactChartPct_selected = cbind(dfReact$Description, dfReact$Pct_selected,
  476. rep("Percentage selected", nrow(dfReact)),
  477. paste0(dfReact$Pct_selected, "%"))
  478. dfReactChartPct_background = cbind(dfReact$Description, dfReact$Pct_background,
  479. rep("Percentage background", nrow(dfReact)),
  480. paste0(dfReact$Pct_background, "%"))
  482. dfReactChart = rbind(dfReactChartP, dfReactChartPct_selected, dfReactChartPct_background)
  483. dfReactChart = as.data.frame(dfReactChart)
  484. colnames(dfReactChart) = c("Pathway", "Value", "Type", "Text")
  485. dfReactChart$Value = as.numeric(dfReactChart$Value)
  487. # Strip the REACTOME_ prefix for cleaner labels
  488. dfReactChart$Pathway = gsub("REACTOME_", "", dfReactChart$Pathway)
  489. save(dfReactChart, file = "dfReactChart_v3.RData")
  491. # Plot top 10 Reactome pathways (grouped bars: -log10 p, % selected, % background)
  492. pdf("../Report/reactome_top_10.pdf", width = 20, height = 8)
  493. ggplot(data = dfReactChart, aes(x = Pathway, y = Value, fill = factor(Type))) +
  494. geom_bar(position = "dodge", stat = "identity") +
  495. coord_flip() +
  496. ylab("- Log10 p") +
  497. ggtitle("Reactome analysis") +
  498. geom_text(aes(label = Text),
  499. position = position_dodge(width = 1), size = 3) +
  500. theme_bw()
  501. dev.off()
  503. # Save significant Reactome results (p < 0.05) to Excel
  504. toSaveX = ans.react@result[which(ans.react@result$pvalue < 0.05), ]
  505. toSaveX = toSaveX[, -c(6, 7)]
  506. openxlsx::write.xlsx(toSaveX, file = "Reactome_final_selection_vs_all_considered.xlsx")
  508. # =============================================================================
  509. # 10. Gene Ontology enrichment (ORA) on the final Elastic Net gene panel
  510. # =============================================================================
  512. deGenes = allSelectedEnsg
  514. # Map HGNC symbols to Ensembl IDs, then to Entrez IDs (required by enrichGO)
  515. ensgs = genes$ensembl_gene_id[match(deGenes, genes$hgnc_symbol)]
  517. ensgUniverse = genes$ensembl_gene_id[match(rownames(dataFiltered), genes$hgnc_symbol)]
  518. ensgUniverse = ensgUniverse[which(!is.na(ensgUniverse))]
  520. deGenes = unlist(mget(ensgs[which(!is.na(ensgs))], envir = org.Hs.egENSEMBL2EG,
  521. ifnotfound = NA))
  522. symbolUniverse = unlist(mget(ensgUniverse[which(!is.na(ensgUniverse))],
  523. envir = org.Hs.egENSEMBL2EG, ifnotfound = NA))
  524. symbolUniverse = symbolUniverse[which(!is.na(symbolUniverse))]
  526. set.seed(123)
  528. # GO enrichment across all ontologies (for exploration)
  529. ans.go <- enrichGO(gene = deGenes, ont = "All",
  530. OrgDb = "org.Hs.eg.db",
  531. universe = symbolUniverse,
  532. readable = TRUE,
  533. pvalueCutoff = 1, qvalueCutoff = 1, pAdjustMethod = "fdr")
  534. # Print top GO results (use View() interactively for full table)
  535. head(ans.go@result, 20)
  537. # GO enrichment restricted to Biological Process (for the manuscript figure)
  538. ans.go3 <- enrichGO(gene = deGenes, ont = "BP",
  539. OrgDb = "org.Hs.eg.db",
  540. universe = symbolUniverse,
  541. readable = TRUE,
  542. pvalueCutoff = 1, qvalueCutoff = 1, pAdjustMethod = "fdr")
  544. # --- Prepare GO bar chart data (same structure as Reactome) ---
  546. dfGo = ans.go3@result[1:10, ]
  548. noGenes = strsplit(dfGo$GeneRatio, "/")
  549. noGenesPlus = as.numeric(unlist(lapply(noGenes, `[[`, 1)))
  550. noGenesAll = as.numeric(unlist(lapply(noGenes, `[[`, 2)))
  552. noBg = strsplit(dfGo$BgRatio, "/")
  553. noBgPlus = as.numeric(unlist(lapply(noBg, `[[`, 1)))
  554. noBgAll = as.numeric(unlist(lapply(noBg, `[[`, 2)))
  556. dfGo$noGenesPlus = noGenesPlus
  557. dfGo$noGenesAll = noGenesAll
  558. dfGo$noBgPlus = noBgPlus
  559. dfGo$noBgAll = noBgAll
  560. dfGo$Pct_selected = round(100 * dfGo$noGenesPlus / dfGo$noGenesAll, 2)
  561. dfGo$Pct_background = round(100 * dfGo$noBgPlus / dfGo$noBgAll, 2)
  562. dfGo$LogP = -1 * (log10(dfGo$pvalue))
  564. dfGoChartP = cbind(dfGo$Description, dfGo$LogP,
  565. rep("p", nrow(dfGo)),
  566. paste0("p=", round(dfGo$pvalue, 6)))
  567. dfGoChartPct_selected = cbind(dfGo$Description, dfGo$Pct_selected,
  568. rep("Percentage selected", nrow(dfGo)),
  569. paste0(dfGo$Pct_selected, "%"))
  570. dfGoChartPct_background = cbind(dfGo$Description, dfGo$Pct_background,
  571. rep("Percentage background", nrow(dfGo)),
  572. paste0(dfGo$Pct_background, "%"))
  574. dfGoChart = rbind(dfGoChartP, dfGoChartPct_selected, dfGoChartPct_background)
  575. dfGoChart = as.data.frame(dfGoChart)
  576. colnames(dfGoChart) = c("Pathway", "Value", "Type", "Text")
  577. dfGoChart$Value = as.numeric(dfGoChart$Value)
  578. save(dfGoChart, file = "dfGoChart.RData")
  580. # Plot top 10 GO Biological Process terms
  581. pdf("../Report/gene_ontology_top_10.pdf", width = 20, height = 8)
  582. ggplot(data = dfGoChart, aes(x = Pathway, y = Value, fill = factor(Type))) +
  583. geom_bar(position = "dodge", stat = "identity") +
  584. coord_flip() +
  585. ylab("- Log10 p") +
  586. ggtitle("Gene Ontology analysis") +
  587. geom_text(aes(label = Text),
  588. position = position_dodge(width = 1), size = 3) +
  589. theme_bw()
  590. dev.off()
  592. # =============================================================================
  593. # 11. Refit Elastic Net on reduced gene panel (|coefficient| > threshold)
  594. # =============================================================================
  596. selectedVars = rownames(toKeepElNetDF)
  597. selectedVars = gsub("`", "", selectedVars)
  599. # Subset training and test data to the reduced gene panel
  600. trainX_subset = trainX[, c(selectedVars, "group")]
  601. testX_subset = testX[, c(selectedVars, "group")]
  603. # Same CV setup as before
  604. control <- trainControl(method = "repeatedcv",
  605. number = 5,
  606. repeats = 5,
  607. search = "random",
  608. verboseIter = TRUE)
  610. # Refit Elastic Net on the reduced feature set
  611. elasticModel_subset <- train(group ~ .,
  612. data = trainX_subset,
  613. method = "glmnet",
  614. tuneLength = 25,
  615. trControl = control)
  616. print(elasticModel_subset)
  618. # --- Evaluate reduced model on training set ---
  619. predictionsTrain = predict(elasticModel_subset, trainX_subset)
  620. probabilitiesTrain = predict(elasticModel_subset, trainX_subset, type = "prob")
  622. table(predictionsTrain, realRefsTrain)
  623. confusionMatrixTrain = confusionMatrix(predictionsTrain, realRefsTrain)
  624. print(confusionMatrixTrain)
  626. # --- Evaluate reduced model on test set ---
  627. predictionsTest = predict(elasticModel_subset, testX_subset)
  628. probabilitiesTest = predict(elasticModel_subset, testX_subset, type = "prob")
  630. if (interactive()) boxplot(probabilitiesTest[, 1] ~ realRefsTest)
  632. # Thresholded predictions (probability > 0.25 classified as GBM)
  633. thresholdedPredictionsTest = rep("HC", length(predictionsTest))
  634. thresholdedPredictionsTest[which(probabilitiesTest[, 1] > 0.25)] = "GBM"
  636. thresholdedPredictionsTrain = rep("HC", length(realRefsTrain))
  637. thresholdedPredictionsTrain[which(probabilitiesTrain[, 1] > 0.25)] = "GBM"
  639. # Confusion matrices for thresholded predictions
  640. confusionMatrixTrain = confusionMatrix(as.factor(thresholdedPredictionsTrain), realRefsTrain)
  641. print(confusionMatrixTrain)
  643. table(thresholdedPredictionsTest, realRefsTest)
  644. confusionMatrixTest = confusionMatrix(
  645. as.factor(thresholdedPredictionsTest), as.factor(realRefsTest), positive = "GBM"
  646. )
  647. print(confusionMatrixTest)
  649. # Confusion matrix for default-threshold predictions
  650. table(predictionsTest, realRefsTest)
  651. confusionMatrixTest = confusionMatrix(
  652. as.factor(predictionsTest), as.factor(realRefsTest), positive = "GBM"
  653. )
  654. print(confusionMatrixTest)
  656. # =============================================================================
  657. # 12. Three-way GO Biological Process enrichment across pipeline stages
  658. # =============================================================================
  659. #
  660. # Compares GO BP enrichment results across the three nested gene sets produced
  661. # by the pipeline:
  662. # 1. candidates -- TDEA-selected genes (FDR < 0.1, FC > 1.3, expr > 3)
  663. # 2. fullModel -- all genes with non-zero Elastic Net coefficients
  664. # 3. allSelectedEnsg -- final panel (|coefficient| > 0.05)
  665. #
  666. # A Venn diagram shows overlap of significant GO terms between sets,
  667. # produced for both FDR-adjusted (p.adjust < 0.05) and unadjusted (p < 0.05)
  668. # significance thresholds.
  670. library(VennDiagram)
  671. library(grid)
  673. # --- Convert each gene set from HGNC symbols to Entrez IDs ---
  674. # (enrichGO requires Entrez IDs; universe already computed in section 10)
  676. symbolsToEntrez = function(symbols) {
  677. ensIds = genes$ensembl_gene_id[match(symbols, genes$hgnc_symbol)]
  678. ensIds = ensIds[which(!is.na(ensIds))]
  679. entrezIds = unlist(mget(ensIds, envir = org.Hs.egENSEMBL2EG, ifnotfound = NA))
  680. entrezIds[which(!is.na(entrezIds))]
  681. }
  683. entrezCandidates = symbolsToEntrez(candidates) # TDEA-selected
  684. entrezFullModel = symbolsToEntrez(fullModel) # non-zero coefficients
  685. entrezFinalPanel = symbolsToEntrez(allSelectedEnsg) # |coef| > threshold
  687. # --- Run GO BP enrichment for each gene set ---
  689. set.seed(123)
  691. goBpCandidates = enrichGO(
  692. gene = entrezCandidates, ont = "BP", OrgDb = "org.Hs.eg.db",
  693. universe = symbolUniverse, readable = TRUE,
  694. pvalueCutoff = 1, qvalueCutoff = 1, pAdjustMethod = "fdr"
  695. )
  697. goBpFullModel = enrichGO(
  698. gene = entrezFullModel, ont = "BP", OrgDb = "org.Hs.eg.db",
  699. universe = symbolUniverse, readable = TRUE,
  700. pvalueCutoff = 1, qvalueCutoff = 1, pAdjustMethod = "fdr"
  701. )
  703. goBpFinalPanel = enrichGO(
  704. gene = entrezFinalPanel, ont = "BP", OrgDb = "org.Hs.eg.db",
  705. universe = symbolUniverse, readable = TRUE,
  706. pvalueCutoff = 1, qvalueCutoff = 1, pAdjustMethod = "fdr"
  707. )
  709. # Save significant results (unadjusted p < 0.05) to Excel
  710. write.xlsx(
  711. goBpCandidates@result[which(goBpCandidates@result$pvalue < 0.05), -c(6, 7)],
  712. file = "GO_BP_candidates_vs_all.xlsx"
  713. )
  714. write.xlsx(
  715. goBpFullModel@result[which(goBpFullModel@result$pvalue < 0.05), -c(6, 7)],
  716. file = "GO_BP_fullModel_vs_all.xlsx"
  717. )
  718. write.xlsx(
  719. goBpFinalPanel@result[which(goBpFinalPanel@result$pvalue < 0.05), -c(6, 7)],
  720. file = "GO_BP_finalSelection_vs_all.xlsx"
  721. )
  723. # --- Extract significant pathway names at both thresholds ---
  725. # FDR-adjusted (p.adjust < 0.05)
  726. sigAdjCandidates = goBpCandidates@result$Description[which(goBpCandidates@result$p.adjust < 0.05)]
  727. sigAdjFullModel = goBpFullModel@result$Description[which(goBpFullModel@result$p.adjust < 0.05)]
  728. sigAdjFinalPanel = goBpFinalPanel@result$Description[which(goBpFinalPanel@result$p.adjust < 0.05)]
  730. # Unadjusted (p < 0.05)
  731. sigRawCandidates = goBpCandidates@result$Description[which(goBpCandidates@result$pvalue < 0.05)]
  732. sigRawFullModel = goBpFullModel@result$Description[which(goBpFullModel@result$pvalue < 0.05)]
  733. sigRawFinalPanel = goBpFinalPanel@result$Description[which(goBpFinalPanel@result$pvalue < 0.05)]
  735. cat("Significant GO BP terms (FDR-adjusted p < 0.05):\n")
  736. cat(" TDEA candidates:", length(sigAdjCandidates), "\n")
  737. cat(" Non-zero in model:", length(sigAdjFullModel), "\n")
  738. cat(" Final panel:", length(sigAdjFinalPanel), "\n")
  740. cat("\nSignificant GO BP terms (unadjusted p < 0.05):\n")
  741. cat(" TDEA candidates:", length(sigRawCandidates), "\n")
  742. cat(" Non-zero in model:", length(sigRawFullModel), "\n")
  743. cat(" Final panel:", length(sigRawFinalPanel), "\n")
  745. # --- Venn diagram helper ---
  747. vennColors = c("#4477AA", "#EE6677", "#228833")
  748. vennSetNames = c("TDEA candidates", "Non-zero in model", "Final panel")
  750. drawGoVenn = function(pathList, titleText, outFile) {
  751. pdf(outFile, width = 7, height = 7)
  752. vennObj = venn.diagram(
  753. x = setNames(pathList, vennSetNames),
  754. filename = NULL,
  755. fill = vennColors,
  756. alpha = 0.45,
  757. col = vennColors,
  758. lwd = 1.5,
  759. fontfamily = "sans",
  760. cat.fontfamily = "sans",
  761. cat.fontface = "bold",
  762. cat.cex = 1.1,
  763. cex = 1.3,
  764. margin = 0.1,
  765. main = titleText,
  766. main.fontfamily = "sans",
  767. main.cex = 1.3,
  768. main.fontface = "bold"
  769. )
  770. grid.draw(vennObj)
  771. dev.off()
  772. }
  774. # --- Venn diagram: FDR-adjusted significance ---
  776. drawGoVenn(
  777. pathList = list(sigAdjCandidates, sigAdjFullModel, sigAdjFinalPanel),
  778. titleText = "GO BP pathways (FDR-adjusted p < 0.05)",
  779. outFile = "../Report/GO_BP_venn_adjusted.pdf"
  780. )
  782. # --- Venn diagram: unadjusted significance ---
  784. drawGoVenn(
  785. pathList = list(sigRawCandidates, sigRawFullModel, sigRawFinalPanel),
  786. titleText = "GO BP pathways (unadjusted p < 0.05)",
  787. outFile = "../Report/GO_BP_venn_unadjusted.pdf"
  788. )
  790. # --- Save overlapping pathway names ---
  792. sharedAdj = Reduce(intersect, list(sigAdjCandidates, sigAdjFullModel, sigAdjFinalPanel))
  793. sharedRaw = Reduce(intersect, list(sigRawCandidates, sigRawFullModel, sigRawFinalPanel))
  795. cat("\nPathways significant in all three sets (FDR-adjusted):", length(sharedAdj), "\n")
  796. cat("Pathways significant in all three sets (unadjusted):", length(sharedRaw), "\n")
  798. if (length(sharedAdj) > 0) {
  799. write.xlsx(data.frame(Pathway = sharedAdj),
  800. file = "GO_BP_shared_all_three_adjusted.xlsx")
  801. }
  802. if (length(sharedRaw) > 0) {
  803. write.xlsx(data.frame(Pathway = sharedRaw),
  804. file = "GO_BP_shared_all_three_unadjusted.xlsx")
  805. }
  807. # =============================================================================
  808. # 13. Build ranked gene list for GSEA (GBM vs HC, training set)
  809. # =============================================================================
  810. #
  811. # Gene ranking metric: sign(log2 FC) * -log10(Wilcoxon p-value)
  812. # This preserves directionality (GBM-up vs GBM-down) while using the
  813. # significance of the test to weight the magnitude of the ranking.
  814. #
  815. # All expressed genes are used here -- NOT just TDEA candidates --
  816. # because GSEA operates on a full ranked list, not a pre-filtered set.
  817. # =============================================================================
  819. library(clusterProfiler)
  820. library(org.Hs.eg.db)
  821. library(openxlsx)
  822. library(VennDiagram)
  823. library(grid)
  825. # Full (unfiltered) DE results on training set
  826. gbmHc_allGenes = getSummaryTable(
  827. countsAllDeseq[, trainId],
  828. sampleInfoAll[trainId, ],
  829. "GBM", "HC",
  830. all = TRUE
  831. )
  833. fc_col = which(colnames(gbmHc_allGenes) == "Mean fold change GBM-HC")
  834. pval_col = which(colnames(gbmHc_allGenes) == "Wilcoxon p-value ")
  836. fc_vals = as.numeric(gbmHc_allGenes[, fc_col])
  837. pval_vals = as.numeric(gbmHc_allGenes[, pval_col])
  839. # Guard against p = 0 (would produce Inf)
  840. pval_vals[pval_vals == 0] <- .Machine$double.eps
  842. rank_metric = sign(log2(fc_vals)) * (-log10(pval_vals))
  844. # Map HGNC symbols -> Ensembl IDs (gseGO uses ENSEMBL keyType)
  845. gene_symbols = rownames(gbmHc_allGenes)
  846. ensembl_ids = genes$ensembl_gene_id[match(gene_symbols, genes$hgnc_symbol)]
  848. original_gene_list_gbm = stats::setNames(rank_metric, ensembl_ids)
  850. # Drop unmapped and duplicated entries
  851. original_gene_list_gbm = original_gene_list_gbm[!is.na(names(original_gene_list_gbm))]
  852. original_gene_list_gbm = original_gene_list_gbm[!duplicated(names(original_gene_list_gbm))]
  854. gene_list_gbm = na.omit(original_gene_list_gbm)
  855. gene_list_gbm = sort(gene_list_gbm, decreasing = TRUE)
  857. # =============================================================================
  858. # 14. GSEA across multiple seeds (GO Biological Process)
  859. # =============================================================================
  860. #
  861. # GSEA results can vary across permutations, so 10 seeds are run and saved
  862. # individually. The seed=1 run (pvalueCutoff = 1) is kept in memory for
  863. # downstream visualisation and the ORA/GSEA Venn diagram.
  864. # =============================================================================
  865. library(enrichplot)
  866. for (seed in 1:10) {
  867. set.seed(seed)
  868. gse_gbm_iter <- gseGO(
  869. geneList = gene_list_gbm,
  870. ont = "ALL",
  871. keyType = "ENSEMBL",
  872. nPerm = 10000,
  873. minGSSize = 3,
  874. maxGSSize = 800,
  875. pvalueCutoff = 0.05,
  876. verbose = TRUE,
  877. OrgDb = org.Hs.eg.db,
  878. pAdjustMethod = "fdr"
  879. )
  880. write.xlsx(gse_gbm_iter,
  881. file = paste0("gene_set_enrichment_all_hc_vs_gbm_ref_hc3_", seed, ".xlsx"))
  882. save(gse_gbm_iter,
  883. file = paste0("gene_set_enrichment_all_hc_vs_gbm_ref_hc3_", seed, ".RData"))
  884. }
  886. # Re-run seed=1 with pvalueCutoff = 1 and ont = "BP" to match the ORA
  887. # used in the Venn. Full result set needed for top-100 selection below.
  888. set.seed(1)
  889. gse_gbm <- gseGO(
  890. geneList = gene_list_gbm,
  891. ont = "BP",
  892. keyType = "ENSEMBL",
  893. nPerm = 10000,
  894. minGSSize = 3,
  895. maxGSSize = 800,
  896. pvalueCutoff = 1, # no cutoff -- needed to get >= 100 results
  897. verbose = FALSE,
  898. OrgDb = org.Hs.eg.db,
  899. pAdjustMethod = "fdr"
  900. )
  902. # =============================================================================
  903. # 15. GSEA visualisation
  904. # =============================================================================
  906. library(DOSE)
  907. library(ggplot2)
  909. pdf("../Report/gsea_dotplot_BP.pdf", width = 10, height = 8)
  910. print(
  911. dotplot(gse_gbm, showCategory = 10, split = ".sign") +
  912. facet_grid(. ~ .sign) +
  913. theme_bw(base_size = 11) +
  914. theme(
  915. strip.background = element_rect(fill = "#F0F0F0"),
  916. axis.text.y = element_text(size = 9)
  917. )
  918. )
  919. dev.off()
  921. # emapplot requires pairwise term similarity
  922. pdf("../Report/gsea_emapplot_BP.pdf", width = 10, height = 10)
  923. print(emapplot(pairwise_termsim(gse_gbm), showCategory = 15))
  924. dev.off()
  926. # =============================================================================
  927. # 16. Venn diagram: ORA (top 100, no FDR) vs GSEA (top 100, no FDR)
  928. # =============================================================================
  929. #
  930. # Both sets are GO Biological Process terms ranked by raw p-value.
  931. # ORA is the enrichment on the final Elastic Net gene panel (ans.go3,
  932. # computed in section 10). GSEA is the seed=1 BP run above.
  933. #
  934. # "Top 100" is defined as the 100 terms with the smallest raw p-value in
  935. # each analysis, without any FDR or adjusted-p filter.
  936. # =============================================================================
  938. n_top = 100
  940. # ORA: ans.go3 was computed with pvalueCutoff = 1 so all terms are present
  941. ora_results_sorted = ans.go3@result[order(ans.go3@result$pvalue), ]
  942. ora_top = head(ora_results_sorted$Description, n_top)
  944. # GSEA: gse_gbm also has pvalueCutoff = 1
  945. gsea_results_sorted = gse_gbm@result[order(gse_gbm@result$pvalue), ]
  946. gsea_top = head(gsea_results_sorted$Description, n_top)
  948. n_overlap = length(intersect(ora_top, gsea_top))
  949. cat(sprintf(
  950. "\nGO BP Venn summary (top %d by raw p-value):\n ORA terms: %d\n GSEA terms: %d\n Overlap: %d\n",
  951. n_top, length(ora_top), length(gsea_top), n_overlap
  952. ))
  954. # Palette consistent with the three-way Venn in section 12
  955. venn_pal = c("#4477AA", "#EE6677")
  957. pdf("../Report/GO_BP_venn_ORA_vs_GSEA.pdf", width = 6, height = 6)
  958. venn_obj = venn.diagram(
  959. x = list(
  960. `ORA` = ora_top,
  961. `GSEA` = gsea_top
  962. ),
  963. filename = NULL,
  964. fill = venn_pal,
  965. alpha = 0.45,
  966. col = venn_pal,
  967. lwd = 1.8,
  968. fontfamily = "sans",
  969. cat.fontfamily = "sans",
  970. cat.fontface = "bold",
  971. cat.cex = 1.15,
  972. cex = 1.4,
  973. margin = 0.12,
  974. main = sprintf("GO BP top %d pathways: ORA vs GSEA", n_top),
  975. main.fontfamily = "sans",
  976. main.cex = 1.15,
  977. main.fontface = "bold"
  978. )
  979. grid.newpage()
  980. grid.draw(venn_obj)
  981. dev.off()
  983. # Save the overlapping pathway names
  984. shared_ora_gsea = intersect(ora_top, gsea_top)
  985. if (length(shared_ora_gsea) > 0) {
  986. write.xlsx(
  987. data.frame(
  988. Pathway = shared_ora_gsea,
  989. ORA_pvalue = ora_results_sorted$pvalue[
  990. match(shared_ora_gsea, ora_results_sorted$Description)],
  991. ORA_padj = ora_results_sorted$p.adjust[
  992. match(shared_ora_gsea, ora_results_sorted$Description)],
  993. GSEA_pvalue = gsea_results_sorted$pvalue[
  994. match(shared_ora_gsea, gsea_results_sorted$Description)],
  995. GSEA_padj = gsea_results_sorted$p.adjust[
  996. match(shared_ora_gsea, gsea_results_sorted$Description)],
  997. GSEA_NES = gsea_results_sorted$NES[
  998. match(shared_ora_gsea, gsea_results_sorted$Description)]
  999. ),
  1000. file = "GO_BP_shared_ORA_GSEA_top100.xlsx"
  1001. )
  1002. }

basic_model_analysis.R at commit d355d14, under GPL-3.0 · at the source

Overview

Authors: Anna Giczewska1, Krzysztof Pastuszak1,2,3, Leon van Hout4, Edward Post4, Jip Ramaker4, Ken Zwaan4, Anna Supernat1,2, Bart A Westerman4,5
ORCID iDs: Anna Giczewska, Krzysztof Pastuszak, Leon van Hout, Edward Post, Jip Ramaker, Anna Supernat, Bart A Westerman
  1. Laboratory of Translational Oncology, Intercollegiate Faculty of Biotechnology, University of Gdańsk and Medical University of Gdańsk, Gdańsk, Poland
  2. Center of Biostatistics and Bioinformatics, Medical University of Gdańsk, Gdańsk, Poland
  3. Department of Algorithms and System Modeling, Gdansk University of Technology, Gdańsk, Poland
  4. Department of Neurosurgery, Cancer Center Amsterdam, Amsterdam UMC Location VUMC, Amsterdam, the Netherlands
  5. Department of Human Genetics, Cancer Center Amsterdam, Amsterdam UMC Location VUMC, Amsterdam, the Netherlands
Institutions: University of Gdańsk (Poland); Gdańsk Medical University (Poland); Gdańsk University of Technology (Poland); Amsterdam Neuroscience (Netherlands); Amsterdam University Medical Centers (Netherlands)
Journal: Cancer reports (Hoboken, N.J.), volume 9, issue 10, article e70656
Dates: received 17 July 2025; accepted 31 July 2026; published online 30 September 2026; in print October 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1002/cnr2.70656 · PMID 42816786 · PMCID PMC13627264 · OpenAlex W7214976098
Open access: gold, a free copy (OpenAlex)
Status: code verified (how the code is found and verified)
Categories: human (organism), other condition (population), clinical / translational (subfield)
Methods: Statistics, Machine learning, Connectivity
Keywords: bioinformatics, biomarkers, glioblastoma, pathway analysis, platelet RNA
MeSH: Biomarkers, Tumor*, Blood Platelets*, Brain Neoplasms*, Glioblastoma*, RNA*, Case-Control Studies, Female, Gene Expression Regulation, Neoplastic, Humans, Tumor Microenvironment (* major topic)
Topic: Glioma Diagnosis and Treatment (Oncology, Medicine), according to OpenAlex
Funding: Gdański Uniwersytet Medyczny (MNISW 07/IDUB/2019/94); Health~Holland Grants called AI-IMPACT; Toxicity-Atlas
Citations: not cited yet (Europe PMC); 46 references in the paper

Abstract

Background: Glioblastoma (GBM) is the most aggressive form of primary brain cancer. Blood platelets have emerged as biomarkers where their RNA expression pattern indicates a tumor's presence. Our study aimed to characterize these patterns to identify biological processes associated with GBM and to evaluate their potential as a non‐invasive biomarker.

Methods: We analyzed platelet RNA‐seq expression profiles from 84 GBM patients and 319 healthy individuals using two complementary analytical approaches of Gene Set Enrichment Analysis and Threshold‐based Differential Expression Analysis with Elastic Net regularization. Together, these approaches may offer complementary perspectives on biological pathways associated with GBM and the descriptive functional annotation of a predictive platelet RNA gene signature, providing insight into platelet‐tumor interactions.

Results: Our GSEA broad analysis highlighted platelet‐related biological processes, including coagulation, wound healing, immune regulation, and cytoskeletal organization. The Elastic Net framework additionally suggested pathways related to immune activation, cell signaling, and protein folding. Overall, these findings may suggest that platelet‐derived RNA not only has potential as a non‐invasive biomarker but may also provide insight into the systemic biological processes associated with the GBM tumor microenvironment.

Conclusions: Our data suggests that platelet‐derived RNA may reflect biological processes associated with the GBM tumor microenvironment, where wound repair and immune‐related processes may play important roles. Our findings provide a resource for further research into the mechanisms that may underlie the pathophysiology and progression of glioblastoma.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 2 matches between paragraphs and lines of code.

KrzysztofPastuszak/Platelets-GBM

License: GPL-3.0
State: the link answers, verified on 2 November 2026
Evidence: files inventoried
Commit: d355d14795f9a6a7d41b49181cb76265e689adaa, 6 March 2026
Languages: R (1)
Size: 23 files, 1 script
Software Heritage: not archived
Found in: “Data Availability Statement”
Holds: README, license file
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: caret (1 file), clusterProfiler (1 file), DESeq2 (1 file), edgeR (1 file), ggplot2 (1 file), glmnet (1 file), tidyverse (1 file)
Availability: 2 checks, the latest on 2 November 2026: the link answers
  • 2 November 2026: the link answers
  • 3 October 2026: the link answers
3 files
At the source: github.com/KrzysztofPastuszak/Platelets-GBM

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 1 script, each with its path and the digest of its content;
  • 2 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data Availability Statement

Data and scripts can be found https://github.com/KrzysztofPastuszak/Platelets‐GBM (https://github.com/KrzysztofPastuszak/Platelets-GBM). All analyses were conducted using R version 4.1.2 or higher.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 3 October 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 8 authors, 5 keywords, 10 MeSH terms, 3 funders, 46 references.

Cite

This paper

Giczewska, A., Pastuszak, K., van Hout, L., Post, E., Ramaker, J., Zwaan, K., Supernat, A., & Westerman, B. A. (2026). Deepening the Understanding of Platelet-Derived RNA as Biomarker for Glioblastoma Through GSEA, TDEA, and Elastic Net Regularization Analysis. Cancer reports (Hoboken, N.J.), 9(10), e70656. https://doi.org/10.1002/cnr2.70656

BibTeX

@article{giczewska2026deepening,
author = {Giczewska, Anna and Pastuszak, Krzysztof and van Hout, Leon and Post, Edward and Ramaker, Jip and Zwaan, Ken and Supernat, Anna and Westerman, Bart A},
title = {{Deepening the Understanding of Platelet-Derived RNA as Biomarker for Glioblastoma Through GSEA, TDEA, and Elastic Net Regularization Analysis}},
journal = {Cancer reports (Hoboken, N.J.)},
year = {2026},
month = oct,
volume = {9},
number = {10},
pages = {e70656},
publisher = {Wiley},
issn = {2573-8348},
doi = {10.1002/cnr2.70656},
url = {https://doi.org/10.1002/cnr2.70656},
pmid = {42816786},
pmcid = {PMC13627264}
}

RIS

TY - JOUR
AU - Giczewska, Anna
AU - Pastuszak, Krzysztof
AU - van Hout, Leon
AU - Post, Edward
AU - Ramaker, Jip
AU - Zwaan, Ken
AU - Supernat, Anna
AU - Westerman, Bart A
TI - Deepening the Understanding of Platelet-Derived RNA as Biomarker for Glioblastoma Through GSEA, TDEA, and Elastic Net Regularization Analysis
T2 - Cancer reports (Hoboken, N.J.)
J2 - Cancer Rep (Hoboken)
PY - 2026
DA - 2026/10/01
VL - 9
IS - 10
SP - e70656
SN - 2573-8348
PB - Wiley
DO - 10.1002/cnr2.70656
UR - https://doi.org/10.1002/cnr2.70656
LA - en
ER -

CSL-JSON

{
"id": "10.1002/cnr2.70656",
"type": "article-journal",
"title": "Deepening the Understanding of Platelet-Derived RNA as Biomarker for Glioblastoma Through GSEA, TDEA, and Elastic Net Regularization Analysis",
"container-title": "Cancer reports (Hoboken, N.J.)",
"author": [
{
"family": "Giczewska",
"given": "Anna"
},
{
"family": "Pastuszak",
"given": "Krzysztof"
},
{
"family": "van Hout",
"given": "Leon"
},
{
"family": "Post",
"given": "Edward"
},
{
"family": "Ramaker",
"given": "Jip"
},
{
"family": "Zwaan",
"given": "Ken"
},
{
"family": "Supernat",
"given": "Anna"
},
{
"family": "Westerman",
"given": "Bart A"
}
],
"container-title-short": "Cancer Rep (Hoboken)",
"volume": "9",
"issue": "10",
"page": "e70656",
"DOI": "10.1002/cnr2.70656",
"PMID": "42816786",
"PMCID": "PMC13627264",
"ISSN": "2573-8348",
"publisher": "Wiley",
"URL": "https://doi.org/10.1002/cnr2.70656",
"language": "en",
"issued": {
"date-parts": [
[
2026,
10,
1
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1093/noajnl/vdad134 [code]
Longitudinal drug synergy assessment using convolutional neural network image-decoding of glioblastoma single-spheroid cultures.
Journal: Neuro-oncology advances
In common: caret, clinical / translational, other condition, 3 authors
[2] doi:10.1016/j.xcrm.2026.102682 [code]
TET CpG sequence-context-specific DNA demethylation shapes progression of IDH-mutant gliomas.
Journal: Cell reports. Medicine
In common: glmnet, edgeR, clusterProfiler, 3 other tools, other condition, 1 reference, author Bart A Westerman
[3] doi:10.1093/noajnl/vdad073
Screening of predicted synergistic multi-target therapies in glioblastoma identifies new treatment strategies.
Journal: Neuro-oncology advances
In common: clinical / translational, other condition, 2 authors
[4] doi:10.3390/genes16050473 [code]
Enhancing Prognostic Signatures in Glioblastoma with Feature Selection and Regularised Cox Regression.
Journal: Genes
In common: glmnet, caret, edgeR, 3 other tools, clinical / translational, other condition, 1 reference
[5] doi:10.1093/gigascience/giaf093 [code]
The Open Pediatric Cancer Project.
Journal: GigaScience
In common: glmnet, caret, edgeR, 3 other tools, other condition, 2 references
[6] doi:10.3390/ijms27136068 [code]
Loss of Neuropeptide Y Signaling Accompanies the Neural-to-Mesenchymal Transcriptional Transition in Glioblastoma: A Multi-Scale Transcriptomic Analysis.
Journal: International journal of molecular sciences
In common: glmnet, edgeR, clusterProfiler, 2 other tools, other condition, 3 references
[7] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: glmnet, caret, edgeR, 4 other tools
[8] doi:10.1038/s41467-025-58675-9 [code]
Multimodal fusion of radio-pathology and proteogenomics identify integrated glioma subtypes with prognostic and therapeutic opportunities.
Journal: Nature communications
In common: glmnet, caret, clusterProfiler, 2 other tools, clinical / translational, other condition, 2 references
[9] doi:10.1093/brain/awag014 [code]
Study of blood linear RNA nominates CD55 and DLD as early-stage biomarkers for Parkinson's sisease.
Journal: Brain : a journal of neurology
In common: glmnet, caret, clusterProfiler, 3 other tools, clinical / translational, 1 reference
[10] doi:10.1371/journal.pcbi.1012857 [code]
An integrative network approach for longitudinal stratification in Parkinson's disease.
Journal: PLoS computational biology
In common: glmnet, edgeR, clusterProfiler, 3 other tools, clinical / translational, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.1002/cnr2.70656.svg)](https://openscicode.org/badge/doi/10.1002/cnr2.70656)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.1002/cnr2.70656.svg
   :target: https://openscicode.org/badge/doi/10.1002/cnr2.70656
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.1002/cnr2.70656"><img src="https://openscicode.org/badge/doi/10.1002/cnr2.70656.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.