inst/doc/MOFA_example_scMT.R

In MOFA: Multi-Omics Factor Analysis (MOFA)

## ---- warning=FALSE, message=FALSE--------------------------------------------
library(MultiAssayExperiment)
library(MOFA)
library(MOFAdata)
library(ggplot2)

## -----------------------------------------------------------------------------
data("scMT_data")
scMT_data

## -----------------------------------------------------------------------------
MOFAobject <- createMOFAobject(scMT_data)
MOFAobject

## -----------------------------------------------------------------------------
plotDataOverview(MOFAobject, colors=c("#31A354","#377EB8","#377EB8","#377EB8"))

## -----------------------------------------------------------------------------
DataOptions <- getDefaultDataOptions()
DataOptions

## -----------------------------------------------------------------------------
ModelOptions <- getDefaultModelOptions(MOFAobject)
ModelOptions

## -----------------------------------------------------------------------------
TrainOptions <- getDefaultTrainOptions()
TrainOptions$seed <- 2018

# Automatically drop factors that explain less than 1% of variance in all omics
TrainOptions$DropFactorThreshold <- 0.01

TrainOptions

## -----------------------------------------------------------------------------
MOFAobject <- prepareMOFA(
  MOFAobject, 
  DataOptions = DataOptions,
  ModelOptions = ModelOptions,
  TrainOptions = TrainOptions
)

## ---- eval=FALSE--------------------------------------------------------------
#  MOFAobject <- runMOFA(MOFAobject)

## -----------------------------------------------------------------------------
filepath <- system.file("extdata", "scMT_model.hdf5", package = "MOFAdata")
MOFAobject <- loadModel(filepath, MOFAobject)
MOFAobject

## -----------------------------------------------------------------------------
# Calculate the variance explained (R2) per factor in each view 
r2 <- calculateVarianceExplained(MOFAobject)
r2$R2Total

# Variance explained by each factor in each view
head(r2$R2PerFactor)

# Plot it
plotVarianceExplained(MOFAobject)

## -----------------------------------------------------------------------------
# Plot all weights and highlight specific gene markers
plotWeights(
  object = MOFAobject, 
  view = "RNA expression", 
  factor = 1, 
  nfeatures = 0, 
  abs = FALSE, 
  manual = list(c("Zfp42","Esrrb","Morc1","Fbxo15",
                  "Jam2","Klf4","Tcl1","Tbx3","Tex19.1"))
)

# Plot top 10 genes
plotTopWeights(
  object = MOFAobject,
  view = "RNA expression", 
  factor = 1, 
  nfeatures = 10
)

## -----------------------------------------------------------------------------
# Add metadata to the plot
factor1 <- sort(getFactors(MOFAobject,"LF1")[,1])
order_samples <- names(factor1)
df <- data.frame(
  row.names = order_samples,
  culture = getCovariates(MOFAobject, "culture")[order_samples],
  factor = factor1
)

plotDataHeatmap(
  object = MOFAobject, 
  view = "RNA expression", 
  factor = "LF1", 
  features = 20, 
  transpose = TRUE, 
  show_colnames=FALSE, show_rownames=TRUE, # pheatmap options
  cluster_cols = FALSE, annotation_col=df # pheatmap options
)

## -----------------------------------------------------------------------------
plotWeights(
  object = MOFAobject,
  view = "Met Enhancers", 
  factor = 1, 
  nfeatures = 0,
  abs = FALSE,
  scale = FALSE
)

## -----------------------------------------------------------------------------
plotDataHeatmap(
  object = MOFAobject, 
  view = "Met Enhancers", 
  factor = 1, 
  features = 500,
  transpose = TRUE,
  cluster_rows=FALSE, cluster_cols=FALSE, # pheatmap options
  show_rownames=FALSE, show_colnames=FALSE, # pheatmap options
  annotation_col=df  # pheatmap options
)

## -----------------------------------------------------------------------------
# Plot all weights and highlight specific gene markers
plotWeights(
  object = MOFAobject, 
  view="RNA expression", 
  factor = 2, 
  nfeatures = 0, 
  manual = list(c("Krt8","Cald1","Anxa5","Tagln",
                  "Ahnak","Dsp","Anxa3","Krt19")), 
  scale = FALSE,
  abs = FALSE
)

# Plot top 10 genes
plotTopWeights(
  object = MOFAobject, 
  view="RNA expression", 
  factor = 2, 
  nfeatures = 10
)

## -----------------------------------------------------------------------------
plotWeights(
  object = MOFAobject, 
  view = "Met Enhancers", 
  factor = 2, 
  nfeatures = 0, 
  abs = FALSE,
  scale = FALSE
)

## ---- message=FALSE-----------------------------------------------------------
factor2 <- sort(getFactors(MOFAobject,"LF2")[,1])
order_samples <- names(factor2)
df <- data.frame(
  row.names = order_samples,
  culture =  getCovariates(MOFAobject, "culture")[order_samples],
  factor = factor2
)
plotDataHeatmap(
  object = MOFAobject, 
  view = "Met Enhancers", 
  factor = 2, 
  features = 500,
  transpose = TRUE,
  cluster_rows=FALSE, cluster_cols=FALSE,  # pheatmap options
  show_rownames=FALSE, show_colnames=FALSE,  # pheatmap options
  annotation_col=df  # pheatmap options
)

## ---- message=FALSE-----------------------------------------------------------
p <- plotFactorScatter(
  object = MOFAobject, 
  factors=1:2, 
  color_by = "culture")

p + scale_color_manual(values=c("lightsalmon","orangered3"))

## ---- message=FALSE-----------------------------------------------------------
plotFactorBeeswarm(
  object = MOFAobject,
  factors = 3, 
  color_by = "cellular_detection_rate"
)

## ---- message=FALSE-----------------------------------------------------------
cdr <- colMeans(getTrainData(MOFAobject)$`RNA expression`>0,na.rm=TRUE)
factor3 <- getFactors(MOFAobject,
                      factors=3)

foo <- data.frame(factor = as.numeric(factor3), cdr = cdr)
ggplot(foo, aes_string(x = "factor", y = "cdr")) + 
  geom_point() + xlab("Factor 3") +
  ylab("Cellular Detection Rate") +
  stat_smooth(method="lm") +
  theme_bw()

## -----------------------------------------------------------------------------
MOFAobject <- impute(MOFAobject)
nonImputedMethylation <- getTrainData(MOFAobject,
                                      view="Met CpG Islands")[[1]]
imputedMethylation <- getImputedData(MOFAobject,
                                     view="Met CpG Islands")[[1]]

## -----------------------------------------------------------------------------
# non-imputed data
pheatmap::pheatmap(nonImputedMethylation[1:100,1:20],
         cluster_rows = FALSE, cluster_cols = FALSE,
         show_rownames = FALSE, show_colnames = FALSE)

# imputed data
pheatmap::pheatmap(imputedMethylation[1:100,1:20],
         cluster_rows = FALSE, cluster_cols = FALSE,
         show_rownames = FALSE, show_colnames = FALSE)

## ---- message=FALSE-----------------------------------------------------------
# kmeans clustering with K=3 using Factor 1
set.seed(1234)
clusters <- clusterSamples(MOFAobject, k=3, factors=1)

# Scatterplot colored by the predicted cluster labels,
# and shaped by the true culture conditions
plotFactorScatter(
  object = MOFAobject,
  factors = 1:2,
  color_by = clusters,
  shape_by = "culture"
)

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