Setup the Basic R Document:

Setup the Library and determine the different columns present in gbm_data

library(TCGAbiolinks)
library(SummarizedExperiment)

gbm_data <- readRDS("gbm_data.rds")


#see the column names of our glioblastoma data
assayNames(gbm_data)
## [1] "unstranded"       "stranded_first"   "stranded_second"  "tpm_unstrand"    
## [5] "fpkm_unstrand"    "fpkm_uq_unstrand"

Previewing the data

expr_mat <- assay(gbm_data, "tpm_unstrand")
dim(expr_mat)
## [1] 60660   391
print(expr_mat[1:5, 1:5])
##                    TCGA-06-6390-01A-11R-A96S-41 TCGA-06-5411-01A-01R-1849-01
## ENSG00000000003.15                      35.2685                      31.9528
## ENSG00000000005.6                        0.3694                       0.3769
## ENSG00000000419.13                      34.4229                      94.3751
## ENSG00000000457.14                       7.5257                       5.4417
## ENSG00000000460.17                       4.0369                       3.5200
##                    TCGA-06-5411-01A-01R-A96S-41 TCGA-12-3648-01A-01R-A96T-41
## ENSG00000000003.15                      12.4084                      19.5305
## ENSG00000000005.6                        0.0738                       0.0496
## ENSG00000000419.13                      35.5003                      20.5027
## ENSG00000000457.14                       3.9462                       3.7070
## ENSG00000000460.17                       2.6398                       2.2320
##                    TCGA-06-A7TK-01A-21R-A96S-41
## ENSG00000000003.15                      37.2806
## ENSG00000000005.6                        0.3463
## ENSG00000000419.13                      64.0628
## ENSG00000000457.14                      11.2977
## ENSG00000000460.17                       9.5765

Log Transformation

log_mat <- log2(expr_mat + 1)
summary(as.vector(log_mat))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.260   1.263   1.948  19.646

Define Tumor vs Non-Tumor groups

colData(gbm_data)$sample_type
##   [1] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##   [4] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##   [7] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [10] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [13] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [16] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [19] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [22] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [25] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [28] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [31] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [34] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [37] "Recurrent Tumor"     "Primary Tumor"       "Primary Tumor"      
##  [40] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [43] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [46] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [49] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [52] "Recurrent Tumor"     "Primary Tumor"       "Primary Tumor"      
##  [55] "Primary Tumor"       "Primary Tumor"       "Recurrent Tumor"    
##  [58] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [61] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [64] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [67] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [70] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [73] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [76] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [79] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [82] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [85] "Solid Tissue Normal" "Primary Tumor"       "Primary Tumor"      
##  [88] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [91] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [94] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
##  [97] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [100] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [103] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [106] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [109] "Primary Tumor"       "Primary Tumor"       "Solid Tissue Normal"
## [112] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [115] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [118] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [121] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [124] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [127] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [130] "Primary Tumor"       "Primary Tumor"       "Recurrent Tumor"    
## [133] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [136] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [139] "Primary Tumor"       "Recurrent Tumor"     "Primary Tumor"      
## [142] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [145] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [148] "Recurrent Tumor"     "Primary Tumor"       "Primary Tumor"      
## [151] "Primary Tumor"       "Primary Tumor"       "Recurrent Tumor"    
## [154] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [157] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [160] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [163] "Recurrent Tumor"     "Recurrent Tumor"     "Primary Tumor"      
## [166] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [169] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [172] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [175] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [178] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [181] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [184] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [187] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [190] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [193] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [196] "Solid Tissue Normal" "Primary Tumor"       "Solid Tissue Normal"
## [199] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [202] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [205] "Solid Tissue Normal" "Primary Tumor"       "Primary Tumor"      
## [208] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [211] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [214] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [217] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [220] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [223] "Primary Tumor"       "Primary Tumor"       "Recurrent Tumor"    
## [226] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [229] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [232] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [235] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [238] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [241] "Primary Tumor"       "Recurrent Tumor"     "Primary Tumor"      
## [244] "Recurrent Tumor"     "Recurrent Tumor"     "Primary Tumor"      
## [247] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [250] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [253] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [256] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [259] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [262] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [265] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [268] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [271] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [274] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [277] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [280] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [283] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [286] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [289] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [292] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [295] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [298] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [301] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [304] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [307] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [310] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [313] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [316] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [319] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [322] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [325] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [328] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [331] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [334] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [337] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [340] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [343] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [346] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [349] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [352] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [355] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [358] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [361] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [364] "Recurrent Tumor"     "Primary Tumor"       "Primary Tumor"      
## [367] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [370] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [373] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [376] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [379] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [382] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [385] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [388] "Primary Tumor"       "Primary Tumor"       "Primary Tumor"      
## [391] "Primary Tumor"
table(colData(gbm_data)$sample_type)
## 
##       Primary Tumor     Recurrent Tumor Solid Tissue Normal 
##                 372                  14                   5

Compute Difference of Means

tumor_samples  <- colnames(log_mat)[colData(gbm_data)$sample_type == "Primary Tumor"]
normal_samples <- colnames(log_mat)[colData(gbm_data)$sample_type == "Solid Tissue Normal"]

tumor_mean  <- rowMeans(log_mat[, tumor_samples])
normal_mean <- rowMeans(log_mat[, normal_samples])

DGE_df <- data.frame(
    gene = rownames(log_mat),
    normal_expression = normal_mean,
    tumor_expression = tumor_mean,
    tumor_minus_normal = tumor_mean - normal_mean,
    stringsAsFactors = FALSE
)

DGE_df$mean_expression <- rowMeans(DGE_df[, c("normal_expression", "tumor_expression")])

head(DGE_df)
##                                  gene normal_expression tumor_expression
## ENSG00000000003.15 ENSG00000000003.15         3.2196665        5.3780144
## ENSG00000000005.6   ENSG00000000005.6         0.2925813        0.8063988
## ENSG00000000419.13 ENSG00000000419.13         5.8170297        5.6922891
## ENSG00000000457.14 ENSG00000000457.14         2.6133246        2.7336163
## ENSG00000000460.17 ENSG00000000460.17         1.0722900        2.4622600
## ENSG00000000938.13 ENSG00000000938.13         2.7327035        3.3763050
##                    tumor_minus_normal mean_expression
## ENSG00000000003.15          2.1583480       4.2988405
## ENSG00000000005.6           0.5138176       0.5494901
## ENSG00000000419.13         -0.1247406       5.7546594
## ENSG00000000457.14          0.1202917       2.6734704
## ENSG00000000460.17          1.3899701       1.7672750
## ENSG00000000938.13          0.6436016       3.0545043

Create a histogram to visualize the expression differences between genes

hist(DGE_df$tumor_minus_normal, breaks = 60,
     main = "Differential expression: tumor minus normal",
     xlab = "log2 expression difference")
abline(v = 0, lwd = 2)
abline(v = c(-5, 5), lty = 2)

## Volcano Plot

DGE_df$expression_group <- "Similar expression"
DGE_df$expression_group[DGE_df$tumor_minus_normal > 5] <- "Higher in tumor"
DGE_df$expression_group[DGE_df$tumor_minus_normal < -5] <- "Higher in normal"

plot_colors <- c("Similar expression" = "gray70",
                  "Higher in tumor" = "firebrick",
                  "Higher in normal" = "royalblue")

plot(DGE_df$tumor_minus_normal, DGE_df$mean_expression,
     pch = 16, cex = 0.45,
     col = plot_colors[DGE_df$expression_group],
     xlab = "Expression difference: Tumor minus Normal",
     ylab = "Gene Expression",
     main = "Volcano plot: GBM tumor vs normal")
abline(v = 0, lwd = 2)
abline(v = c(-5, 5), lty = 2)
legend("topright", legend = names(plot_colors), col = plot_colors, pch = 16, cex = 0.8, bty = "n")

Order the genes

Clearly, there are some genes that are significantly higher in tumor and higher in normal, so we can determine those by ordering the dataset.

DGE_ordered <- DGE_df[order(DGE_df$tumor_minus_normal, decreasing = TRUE), ]

# Top genes higher in tumor
top_up <- head(DGE_ordered$gene, 15)
top_up
##  [1] "ENSG00000263740.2" "ENSG00000265735.2" "ENSG00000200488.1"
##  [4] "ENSG00000200312.1" "ENSG00000252010.1" "ENSG00000201428.1"
##  [7] "ENSG00000202058.1" "ENSG00000212232.1" "ENSG00000286522.2"
## [10] "ENSG00000197061.5" "ENSG00000271394.1" "ENSG00000184357.5"
## [13] "ENSG00000200087.1" "ENSG00000239899.3" "ENSG00000276168.1"
# Top genes higher in normal
top_down <- tail(DGE_ordered$gene, 15)
top_down
##  [1] "ENSG00000198695.2"  "ENSG00000074317.11" "ENSG00000166448.15"
##  [4] "ENSG00000130540.14" "ENSG00000181418.8"  "ENSG00000006116.4" 
##  [7] "ENSG00000163032.12" "ENSG00000008056.14" "ENSG00000176884.16"
## [10] "ENSG00000124507.11" "ENSG00000126583.11" "ENSG00000104722.14"
## [13] "ENSG00000154146.13" "ENSG00000157005.4"  "ENSG00000104888.10"

From the data, it seems that genes ENSG00000263740.2, ENSG00000265735.2, and ENSG00000200488.1 are significantly more expressed in tumors. On the other hand, genes ENSG00000198695.2, ENSG00000074317.11, and ENSG00000166448.15 are expressed more in normal brain tissue.

Create a Heatmap

heatmap_genes <- c(top_up, top_down)

set.seed(1)
sample_subset <- c(normal_samples, sample(tumor_samples, 30)) # using only a random 30 columns of samples since there is too much for the heatmap

sum(heatmap_genes %in% rownames(log_mat))   
## [1] 30
length(heatmap_genes)                        
## [1] 30
missing_genes <- heatmap_genes[!heatmap_genes %in% rownames(log_mat)]

expr_sub <- log_mat[heatmap_genes, sample_subset]



expr_sub_scaled <- t(scale(t(expr_sub)))

heat_colors <- colorRampPalette(c("blue", "white", "red"))(100)

heatmap(
  expr_sub_scaled,
  labRow = heatmap_genes,   
  labCol = "",
  margins = c(3, 8),
  xlab = "samples (tumor + normal)",
  ylab = "possible biomarker genes",
  col = heat_colors,
  zlim = c(-2, 2),
  main = "Biomarker genes: GBM tumor vs normal"
)

Clearly, through these 15 genes, there is a separation of tumor and non-tumor samples.

We can combine all of these genes to create a biomarker.

But now, we need to test the feasability and accuracy in detecting whether a sample is a tumor or not.

Final Stage: Testing our Biomarker

# first, we create an overall z-score for our biomarker

all_samples <- c(tumor_samples, normal_samples)
expr_all <- log_mat[heatmap_genes, all_samples]
expr_all_scaled <- t(scale(t(expr_all)))

# average z-score of tumor-up genes - average z-score of normal-up genes = biomarker score
up_score   <- colMeans(expr_all_scaled[top_up, ])
down_score <- colMeans(expr_all_scaled[top_down, ])
biomarker_score <- up_score - down_score
# create a new dataframe with the corresponding score 
score_df <- data.frame(
    sample = all_samples,
    group  = colData(gbm_data)$sample_type[match(all_samples, colnames(log_mat))],
    score  = biomarker_score
)

Next, we can create a boxplot to visualize whether the biomarker accurately separates the boxes enough.

Boxplot

# creates 2 boxplots of the dataframe
boxplot(score ~ group, data = score_df,
        main = "Composite biomarker score: Tumor vs Normal",
        ylab = "Biomarker score", col = c("lightblue", "salmon"))

Clearly, there is a statistically significant difference between the primary tumor group and hte solid tissue normal group, demonstrating that our biomarker is viable.