## to install MOVICS
### note: there are many dependencies; you may get an error
### for a missing package; download and try again
### it may take several times
### you can refer to their IMPORTS file from their DESCRIPTION file on Github
### for a list of dependencies:
### https://github.com/xlucpu/MOVICS/blob/master/DESCRIPTION
# if (!requireNamespace("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
# if (!require("devtools"))
# install.packages("devtools")
# devtools::install_github("xlucpu/MOVICS")
install.packages("~/Downloads/SNFtool_2.3.0.tar.gz", repos = NULL, type="source")
library(MOVICS)
set.seed(4444)
library("jpeg")
jj <- readJPEG("MOVICS_pipeline.jpeg",native=TRUE)
plot(0:1,0:1,type="n",ann=FALSE,axes=FALSE)
rasterImage(jj,0,0,1,1)
jj <- readJPEG("methods_comparison.jpeg",native=TRUE)
plot(0:1,0:1,type="n",ann=FALSE,axes=FALSE)
rasterImage(jj,0,0,1,1)
# install.packages("Rfssa") if you want to download through github
library(Rfssa)
url <- "https://github.com/KechrisLab/ASAShortCourse-MultiOmics/blob/main/Lecture%202/brca_dat.Rdata"
load_github_data(url)
# load("brca_dat.Rdata")
# let's get a quick look at our data
names(brca_dat)
paste("dim of clinical data:", dim(brca_dat[["clinical"]]))
head(brca_dat[["clinical"]])
# check sample names all match
identical(brca_dat[["clinical"]]$bcr_patient_barcode, colnames(brca_dat[["MO"]][["Expression"]]))
identical(brca_dat[["clinical"]]$bcr_patient_barcode, colnames(brca_dat[["MO"]][["Methylation"]]))
identical(brca_dat[["clinical"]]$bcr_patient_barcode, colnames(brca_dat[["MO"]][["miRNA"]]))
identical(colnames(brca_dat[["MO"]][["Expression"]]), colnames(brca_dat[["MO"]][["Methylation"]]))
identical(colnames(brca_dat[["MO"]][["Expression"]]), colnames(brca_dat[["MO"]][["miRNA"]]))
identical(colnames(brca_dat[["MO"]][["Methylation"]]), colnames(brca_dat[["MO"]][["miRNA"]]))
# should all be TRUE (6)
paste("names of MO data:", names(brca_dat[["MO"]]))
paste("dim of mRNA data:", dim(brca_dat[["MO"]][["Expression"]]))
paste("dim of methylation data:", dim(brca_dat[["MO"]][["Methylation"]]))
paste("dim of miRNA data:", dim(brca_dat[["MO"]][["miRNA"]]))
# data checking -- are there any missing values?
sum(is.na(brca_dat[["MO"]][["Expression"]]))
sum(is.na(brca_dat[["MO"]][["Methylation"]]))
sum(is.na(brca_dat[["MO"]][["miRNA"]]))
# exp
range(brca_dat[["MO"]][["Expression"]])
plot(density(brca_dat[["MO"]][["Expression"]]), main = "Expression")
# methyl
range(brca_dat[["MO"]][["Methylation"]])
plot(density(brca_dat[["MO"]][["Methylation"]]), main = "Methylation")
# miRNA
range(brca_dat[["MO"]][["miRNA"]])
plot(density(brca_dat[["MO"]][["miRNA"]]), main = "miRNA")
jj <- readJPEG("ex_breastcancer_pic.jpeg")
plot(0:1,0:1,type="n",ann=FALSE,axes=FALSE)
rasterImage(jj,0,0,1,1)
# https://doi.org/10.1016/B978-0-12-800886-7.00021-2
optk.brca <- getClustNum(data = brca_dat[["MO"]],
is.binary = c(F,F,F), # all omics data is continuous (not binary)
try.N.clust = 2:8, # try cluster number from 2 to 8
fig.name = "CLUSTER NUMBER OF TCGA-BRCA")
optk.brca
# what if we use the suggested k=3 clusters?
# you don't need to run this during lab; I'm just presenting it as an example
mo_rslts_3 <- getMOIC(data = brca_dat[["MO"]],
methodslist = list("SNF", "PINSPlus", "NEMO", "LRAcluster", "IntNMF"),
N.clust = 3,
type = c("gaussian", "gaussian", "gaussian"))
cmoic.brca_3 <- getConsensusMOIC(moic.res.list = mo_rslts_3,
fig.name = "CONSENSUS HEATMAP - 3 Clusters",
distance = "euclidean",
linkage = "average")
getSilhouette(sil = cmoic.brca_3$sil, # a sil object returned by getConsensusMOIC()
fig.path = getwd(),
fig.name = "SILHOUETTE",
height = 5.5,
width = 5)
mo_rslts <- getMOIC(data = brca_dat[["MO"]],
methodslist = list("SNF", "PINSPlus", "NEMO", "LRAcluster", "IntNMF"),
N.clust = 4, # set number of clusters
type = c("gaussian", "gaussian", "gaussian")) # what is the distribution of the datasets in MO list (same order)
cmoic.brca <- getConsensusMOIC(moic.res.list = mo_rslts,
fig.name = "CONSENSUS HEATMAP - 4 Clusters",
distance = "euclidean",
linkage = "ward.D")
getSilhouette(sil = cmoic.brca$sil, # a sil object returned by getConsensusMOIC()
fig.path = getwd(),
fig.name = "SILHOUETTE",
height = 5.5,
width = 5)
# convert beta value to M value for stronger signal
std_dat <- brca_dat[["MO"]]
std_dat[["Methylation"]] <- log2(std_dat[["Methylation"]] / (1 - std_dat[["Methylation"]]))
# data normalization for heatmap
plotdata <- getStdiz(data = std_dat,
halfwidth = c(2,2,2), # no truncation for mutation
centerFlag = c(T,T,T), # no center for mutation
scaleFlag = c(T,T,T)) # no scale for mutation
mRNA.col <- c("#00FF00", "#008000", "#000000", "#800000", "#FF0000")
meth.col <- c("#0074FE", "#96EBF9", "#FEE900", "#F00003")
miRNA.col <- c("#6699CC", "white" , "#FF3C38")
col.list <- list(mRNA.col, meth.col, miRNA.col)
# extract PAM50, pathologic stage for sample annotation
annCol <- brca_dat[["clinical"]][,c("BRCA_Subtype_PAM50"), drop = FALSE]
# generate corresponding colors for sample annotation
annColors <- list(
BRCA_Subtype_PAM50 = c("Basal" = "blue",
"Her2" = "red",
"LumA" = "yellow",
"LumB" = "green",
"Normal" = "black")
)
# comprehensive heatmap
getMoHeatmap(data = plotdata,
row.title = names(std_dat),
is.binary = c(F,F,F), # we don't have any binary omics data (ex mutation)
legend.name = c("mRNA","M value","miRNA"),
clust.res = mo_rslts$SNF$clust.res, # cluster results for SNF
color = col.list,
# annCol = annCol, # annotation for samples (if you want to show PAM50 classes too)
# annColors = annColors, # annotation color
width = 10, # width of each subheatmap
height = 5, # height of each subheatmap
fig.name = "COMPREHENSIVE HEATMAP OF SNF")
getMoHeatmap(data = plotdata,
row.title = names(std_dat),
is.binary = c(F,F,F), # all data is continuous
legend.name = c("mRNA","M value","miRNA"),
clust.res = mo_rslts$PINSPlus$clust.res, # cluster results for PINSPlus
color = col.list,
width = 10, # width of each subheatmap
height = 5, # height of each subheatmap
fig.name = "COMPREHENSIVE HEATMAP OF PINSPlus")
# comprehensive heatmap (may take a while)
getMoHeatmap(data = plotdata,
row.title = names(std_dat),
is.binary = c(F,F,F),
legend.name = c("mRNA","M value","miRNA"),
clust.res = mo_rslts$NEMO$clust.res, # cluster results for NEMO
color = col.list,
width = 10, # width of each subheatmap
height = 5, # height of each subheatmap
fig.name = "COMPREHENSIVE HEATMAP OF PINSPlus")
# comprehensive heatmap (may take a while)
getMoHeatmap(data = plotdata,
row.title = names(std_dat),
is.binary = c(F,F,F),
legend.name = c("mRNA","M value","miRNA"),
clust.res = mo_rslts$LRAcluster$clust.res, # cluster results for LRAcluster
color = col.list,
width = 10,
height = 5,
fig.name = "COMPREHENSIVE HEATMAP OF PINSPlus")
# comprehensive heatmap (may take a while)
getMoHeatmap(data = plotdata,
row.title = names(std_dat),
is.binary = c(F,F,F),
legend.name = c("mRNA","M value","miRNA"),
clust.res = mo_rslts$IntNMF$clust.res, # cluster results for intNMF
color = col.list,
width = 10, # width of each subheatmap
height = 5, # height of each subheatmap
fig.name = "COMPREHENSIVE HEATMAP OF IntNMF")
getMoHeatmap(data = plotdata,
row.title = names(plotdata),
is.binary = c(F,F,F), # no binary omics data
legend.name = c("mRNA","M value","miRNA"),
clust.res = cmoic.brca$clust.res, # consensusMOIC results
clust.dend = NULL, # show no dendrogram for samples
show.colnames = FALSE, # show no sample names
show.row.dend = c(T,T,T), # show dendrogram for features
annRow = NULL, # no selected features
color = col.list,
annCol = annCol, # annotation for samples
annColors = annColors, # annotation color
width = 10, # width of each subheatmap
height = 5, # height of each subheatmap
fig.name = "COMPREHENSIVE HEATMAP OF CONSENSUSMOIC")
clust_rslts_SNF_df <- data.frame(mo_rslts$SNF$clust.res)
colnames(clust_rslts_SNF_df) <- c("samID", "SNF")
clust_rslts_CM_df <- data.frame(cmoic.brca$clust.res)
colnames(clust_rslts_CM_df) <- c("samID", "Consensus")
clust_rslts_df <- merge(clust_rslts_SNF_df, clust_rslts_CM_df, by="samID")
# head(clust_rslts_df)
table(clust_rslts_df$SNF, clust_rslts_df$Consensus)
# survival comparison
brca_dat[["clinical"]]$futime = as.numeric(brca_dat[["clinical"]]$futime)
head(brca_dat[["clinical"]])
surv.brca <- compSurv(moic.res = cmoic.brca,
surv.info = brca_dat[["clinical"]],
convt.time = "m", # convert day unit to month
surv.median.line = "h", # draw horizontal line at median survival
xyrs.est = c(5,10), # estimate 5 and 10-year survival
fig.name = "KAPLAN-MEIER CURVE OF CONSENSUSMOIC")
print(surv.brca)
# clinVars_df <- brca_dat[["clinical"]][,c("ajcc_pathologic_stage", "age_at_diagnosis","ajcc_pathologic_t", "ajcc_pathologic_n","ajcc_pathologic_m")]
clinVars_df <- brca_dat[["clinical"]]
rownames(clinVars_df) <- clinVars_df$bcr_patient_barcode
clinVars_df <- clinVars_df[,c( "ajcc_pathologic_stage", "ajcc_pathologic_t", "ajcc_pathologic_n", "ajcc_pathologic_m", "age_at_diagnosis")]
head(clinVars_df)
clin.brca <- compClinvar(moic.res = cmoic.brca,
var2comp = clinVars_df, # data.frame needs to summarize (must has row names of samples)
strata = "Subtype", # stratifying variable (e.g., Subtype in this example)
# factorVars = c("ajcc_pathologic_stage"), # features that are considered categorical variables
factorVars = c("ajcc_pathologic_stage", "ajcc_pathologic_t", "ajcc_pathologic_n", "ajcc_pathologic_m"), # features that are considered categorical variables
nonnormalVars = "age_at_diagnosis", # feature(s) that are considered using nonparametric test
exactVars = NULL, # feature(s) that are considered using exact test
doWord = FALSE, # generate .docx file in local path
tab.name = "SUMMARIZATION OF CLINICAL FEATURES")
clin.brca
# compare agreement with other subtypes
sub_df = data.frame(
BRCA_Subtype_PAM50 = brca_dat[["clinical"]][,c("BRCA_Subtype_PAM50")])
rownames(sub_df) = brca_dat[["clinical"]][,c("bcr_patient_barcode")]
head(sub_df)
# agreement comparison (support up to 6 classifications include current subtype)
agree.brca <- compAgree(moic.res = cmoic.brca,
subt2comp = sub_df,
doPlot = TRUE,
box.width = 0.2,
fig.name = "AGREEMENT OF CONSENSUSMOIC WITH PAM50 Subtype")
# run DEA with limma
runDEA(dea.method = "limma",
expr = brca_dat[["MO"]][["Expression"]], # normalized expression data
moic.res = cmoic.brca,
overwt = T,
res.path = getwd(), # path to save marker files
prefix = "de_TCGA-BRCA")
# choose limma result to identify subtype-specific DOWN-regulated biomarkers
marker.dn <- runMarker(moic.res = cmoic.brca,
dea.method = "limma",
prefix = "de_TCGA-BRCA",
dirct = "down",
dat.path = getwd(),
res.path = getwd(),
p.cutoff = 0.05, # p cutoff to identify significant DEGs
p.adj.cutoff = 0.05, # padj cutoff to identify significant DEGs
n.marker = 200, # number of biomarkers for each subtype
doplot = T,
annCol = annCol,
annColors = annColors,
norm.expr = brca_dat[["MO"]][["Expression"]],
fig.name = "UPREGULATED BIOMARKER HEATMAP")
# subtype-specific UP-regulated biomarkers
marker.up <- runMarker(moic.res = cmoic.brca,
dea.method = "limma",
prefix = "de_TCGA-BRCA",
dirct = "up",
dat.path = getwd(),
res.path = getwd(),
p.cutoff = 0.05, # p cutoff to identify significant DEGs
p.adj.cutoff = 0.05, # padj cutoff to identify significant DEGs
n.marker = 200, # number of biomarkers for each subtype
doplot = T,
annCol = annCol,
annColors = annColors,
norm.expr = brca_dat[["MO"]][["Expression"]],
fig.name = "UPREGULATED BIOMARKER HEATMAP")
# MUST locate ABSOLUTE path of msigdb file
MSIGDB.FILE <- system.file("extdata", "c5.bp.v7.1.symbols.xls", package = "MOVICS", mustWork = TRUE)
# # run GSEA to identify DOWN-regulated GO pathways using results from edgeR
gsea.down <- runGSEA(moic.res = cmoic.brca,
dea.method = "limma", # name of DEA method
prefix = "de_TCGA-BRCA", # MUST be the same of argument in runDEA()
dat.path = getwd(), # path of DEA files
res.path = getwd(), # path to save GSEA files
msigdb.path = MSIGDB.FILE, # MUST be the ABSOLUTE path of msigdb file
norm.expr = brca_dat[["MO"]][["Expression"]], # use normalized expression to calculate enrichment score
dirct = "down", # direction of dysregulation in pathway
p.cutoff = 0.05, # p cutoff to identify significant pathways
p.adj.cutoff = 0.25, # padj cutoff to identify significant pathways
gsva.method = "gsva", # method to calculate single sample enrichment score
norm.method = "mean", # normalization method to calculate subtype-specific enrichment score
fig.name = "DOWNREGULATED PATHWAY HEATMAP")
data.frame(gsea.down$gsea.list$CS2[1:6,3:6])
head(round(gsea.down$grouped.es,3))
# # run GSEA to identify up-regulated GO pathways using results from limma
gsea.up <- runGSEA(moic.res = cmoic.brca,
dea.method = "limma", # name of DEA method
prefix = "detesting_TCGA-BRCA", # MUST be the same of argument in runDEA()
dat.path = getwd(), # path of DEA files
res.path = getwd(), # path to save GSEA files
msigdb.path = MSIGDB.FILE, # MUST be the ABSOLUTE path of msigdb file
norm.expr = brca_dat[["MO"]][["Expression"]], # use normalized expression to calculate enrichment score
dirct = "up", # direction of dysregulation in pathway
p.cutoff = 0.05, # p cutoff to identify significant pathways
p.adj.cutoff = 0.25, # padj cutoff to identify significant pathways
gsva.method = "gsva", # method to calculate single sample enrichment score
norm.method = "mean", # normalization method to calculate subtype-specific enrichment score
fig.name = "UPREGULATED PATHWAY HEATMAP")
# MUST locate ABSOLUTE path of gene set file
GSET.FILE <-
system.file("extdata", "gene sets of interest.gmt", package = "MOVICS", mustWork = TRUE)
# run GSVA to estimate single sample enrichment score based on given gene set of interest
gsva.res <-
runGSVA(moic.res = cmoic.brca,
norm.expr = brca_dat[["MO"]][["Expression"]],
gset.gmt.path = GSET.FILE, # ABSOLUTE path of gene set file
gsva.method = "gsva", # method to calculate single sample enrichment score
annCol = annCol,
annColors = annColors,
fig.path = getwd(),
fig.name = "GENE SETS OF INTEREST HEATMAP",
height = 5,
width = 10)
message("check raw enrichment score")
gsva.res$raw.es[1:3,1:3]
message("check z-scored and truncated enrichment score")
gsva.res$scaled.es[1:3,1:3]
Datasets
Models
