inst/doc/TNBC.CMS.R

In TNBC.CMS: TNBC.CMS: Prediction of TNBC Consensus Molecular Subtypes

## ----include=FALSE------------------------------------------------------------
library(knitr)
opts_chunk$set(concordance=TRUE)
options(scipen = 1, digits = 2, warn = -1)

## ----loadPackagesAndData,message=FALSE----------------------------------------
library(TNBC.CMS)
data("GSE25055")
dim(assays(GSE25055)[[1]])
assays(GSE25055)[[1]][1:5, 1:5]

## ----prediction---------------------------------------------------------------
predictions <- predictCMS(GSE25055)
table(predictions)
head(attr(predictions, "probabilities"))

## ----include=FALSE------------------------------------------------------------
resultGES <- t(as.matrix(colData(GSE25055)[,3:9]))

## ----eval=FALSE---------------------------------------------------------------
#  resultGES <- computeGES(expr = GSE25055, pred = predictions,
#                          rnaseq = FALSE)

## ----eval=TRUE----------------------------------------------------------------
resultGES[,1:4]

## ----label=gesplot,echo=FALSE,fig.height=4,fig.cap="Gene expression signature scores",fig.pos="H",message=FALSE----
library(grid)
library(ggpubr)
CMS.palette <- c("MSL" = "brown2", "IM" = "gold2",
                   "LAR" = "yellowgreen", "SL" = "midnightblue")
pval1 <- 8.1e-7
pval2 <- 3.6e-6
pval3 <- 0.016
pval4 <- 0.00034
TITLE_SIZE <- 12
SUBTITLE_SIZE <- 10

sigdat <- data.frame(t(resultGES))
sigdat$CMS <- predictions

sub1 <- bquote(paste("Wilcoxon (MSL vs. others) ",
                         italic(p), " = ", .(format(pval1, digits = 2))))
sub2 <- bquote(paste("Wilcoxon (MSL vs. others) ",
                         italic(p), " = ", .(format(pval2, digits = 2))))
sub3 <- bquote(paste("Wilcoxon (MSL vs. others) ",
                         italic(p), " = ", .(format(pval3, digits = 2))))
sub4 <- bquote(paste("Wilcoxon (MSL vs. others) ",
                         italic(p), " = ", .(format(pval4, digits = 2))))

p1 <- ggboxplot(sigdat, x = "CMS", y = "Stromal",
                fill = "CMS", palette = CMS.palette) +
  theme_bw() + labs(title = "Stromal", subtitle = sub1) +
  theme(legend.position = "none", axis.title = element_blank(),
        axis.text.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid = element_blank(),
        plot.title = element_text(size = TITLE_SIZE,
                                  hjust = 0.5, face = "bold"),
        plot.subtitle = element_text(size=  SUBTITLE_SIZE, hjust = 0.5))

p2 <- ggboxplot(sigdat, x = "CMS", y = "Immune",
                fill = "CMS", palette = CMS.palette) +
  theme_bw() + labs(title = "Immune", subtitle = sub2) +
  theme(legend.position = "none", axis.title = element_blank(),
        axis.text.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid = element_blank(),
        plot.title = element_text(size = TITLE_SIZE,
                                  hjust = 0.5, face = "bold"),
        plot.subtitle = element_text(size=  SUBTITLE_SIZE, hjust = 0.5))

p3 <- ggboxplot(sigdat, x = "CMS", y = "Hormone",
                fill = "CMS", palette = CMS.palette) +
  theme_bw() + labs(title = "Hormone", subtitle = sub3) +
  theme(legend.position = "none", axis.title = element_blank(),
        axis.text.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid = element_blank(),
        plot.title = element_text(size = TITLE_SIZE,
                                  hjust = 0.5, face = "bold"),
        plot.subtitle = element_text(size=  SUBTITLE_SIZE, hjust = 0.5))

p4 <- ggboxplot(sigdat, x = "CMS", y = "Stemness",
                fill = "CMS", palette = CMS.palette) +
  theme_bw() + labs(title = "Stem", subtitle = sub4) +
  theme(legend.position = "none", axis.title = element_blank(),
        axis.text.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid = element_blank(),
        plot.title = element_text(size = TITLE_SIZE,
                                  hjust = 0.5, face = "bold"),
        plot.subtitle = element_text(size=  SUBTITLE_SIZE, hjust = 0.5))

g1 <- ggplotGrob(p1)
g2 <- ggplotGrob(p2)
g3 <- ggplotGrob(p3)
g4 <- ggplotGrob(p4)
g <- cbind(rbind(g1, g3, size = "first"),
           rbind(g2,  g4, size = "first"), size = "first")
grid.newpage()
grid.draw(g)

## ----label=gsa,fig.height=6,fig.cap="GSVA enrichment scores",fig.pos="H",fig.show="hold"----
resultGSVA <- performGSVA(expr = GSE25055, pred = predictions,
                          gene.set = NULL)
head(resultGSVA[,1:4])

## ----label=os,fig.height=3,fig.width=5,fig.cap="Overall survival",fig.pos="H",fig.show="hold",fig.align="center"----
time <- colData(GSE25055)$DFS.month
event <- colData(GSE25055)$DFS.status
plotKM(pred = predictions, time = time, event = event)

## ----label=hr,fig.height=3,fig.cap="Forest plot of hazard ratios",fig.pos="H",fig.show="hold"----
library(survival)

#Test for difference of survival between low and high expression groups
surv <- Surv(time, event)
GSE25055.exprs <- assays(GSE25055)[[1]]
chisq <- apply(GSE25055.exprs, 1, function(x) survdiff(surv ~ (x > median(x)))$chisq)
pval <- 1 - pchisq(chisq, 1)

#Select 10 genes with lowest p-values for the log-rank test
gs <- names(sort(pval)[1:10])
gs
plotHR(expr = GSE25055, gene.symbol = gs, pred = predictions, time = time,
       event = event, by.subtype = FALSE)

## ----label=hrsub,fig.height=4,fig.cap="Forest plot of subtype-specific hazard ratios",fig.pos="H",fig.show="hold"----
plotHR(expr = GSE25055, gene.symbol = gs[1:4], pred = predictions,
       time = time, event = event, by.subtype = TRUE)

## ----label=ds,fig.height=3,fig.cap="Drug signature scores",fig.pos="H",fig.show="hold"----
resultDS <- computeDS(expr = GSE25055, pred = predictions)
head(resultDS[,1:4])

## ----saveResults--------------------------------------------------------------
dfCMS <- data.frame(row.names = colnames(GSE25055.exprs), CMS = predictions, t(resultGES),
                    t(resultDS), stringsAsFactors = FALSE)
head(dfCMS)
write.table(dfCMS, file = "GSE25055_CMS.txt")

## ----sessionInfo--------------------------------------------------------------
sessionInfo()