In lujunyan1118/mofaCLL: Multi-omics reveals clinically relevant proliferative drive associated with mTOR-MYC-OXPHOS activity in chronic lymphocytic leukemia

Prepare MOFA model

plotDir = ifelse(exists(".standalone"), "", "part1/")
if(plotDir!="") if(!file.exists(plotDir)) dir.create(plotDir)
knitr::opts_chunk$set(fig.path=plotDir, dev=c("png", "pdf"))


library(mofaCLL)
library(reticulate)
library(DESeq2)
library(sva)
library(MultiAssayExperiment)
library(MOFA)
library(tidyverse)
options(stringsAsFactors=FALSE)

#set the global ggplot theme
theme_set(theme_bw() + theme(axis.text = element_text(size=12), 
                             axis.title = element_text(size=14),
                             plot.title = element_text(size = 15, hjust =0.5, face="bold")))

Prepare each single omic dataset

Load datasets

data("drug","gene","rna")

Download methylation dataset from server to a temporary folder

methPath <- file.path(tempdir(),"meth.RData")
if(!file.exists(methPath)) {
  download.file("https://www.huber.embl.de/users/jlu/data/meth.RData",
                methPath)
  load(methPath)
} else {
  load(methPath)
}

Drug screening

viabMat <- mutate(drug, id = paste0(Drug,"_",concIndex)) %>%
  select(patientID, id, normVal) %>%
  spread(key = patientID, value = normVal) %>%
  data.frame() %>% column_to_rownames("id") %>% 
  as.matrix()

RNA sequencing

rna.vst<-DESeq2::varianceStabilizingTransformation(rna)
exprMat <- assay(rna.vst)
nTop = 5000
sds <- genefilter::rowSds(exprMat)
exprMat <- exprMat[order(sds, decreasing = T)[1:nTop],]

DNA methylation array

methData <- assays(meth)[["beta"]]
nTop = 5000
sds <- genefilter::rowSds(methData)
methData <- methData[order(sds, decreasing = T)[1:nTop],]

Genetics

gene <- t(gene)

Assemble multiAssayExperiment object

Create object

mofaData <- list(Drugs = viabMat, 
                mRNA = exprMat,
                Mutations = gene, 
                Methylation = methData)

# Create MultiAssayExperiment object 
mofaData <- MultiAssayExperiment::MultiAssayExperiment(
  experiments = mofaData
)

Only keep samples that have at least three assays

useSamples <- MultiAssayExperiment::sampleMap(mofaData) %>%
  as_tibble() %>% group_by(primary) %>% summarise(n= length(assay)) %>%
  filter(n >= 3) %>% pull(primary)
mofaData <- mofaData[,useSamples]

Dimensions for each dataset

experiments(mofaData)

How many samples have the complete datasets

table(table(sampleMap(mofaData)$primary))

Build MOFA object

Build MOFA object from multiAssayExperiment object

# Build the MOFA object
MOFAobject <- createMOFAobject(mofaData)
MOFAobject

Overview of data

p <- plotDataOverview.m(MOFAobject)
p$plot

Setup MOFA training parameters

Define data options

DataOptions <- getDefaultDataOptions()

Define model options

ModelOptions <- getDefaultModelOptions(MOFAobject)
ModelOptions$numFactors <- 30
ModelOptions$sparsity <- TRUE
ModelOptions

Define training options

TrainOptions <- getDefaultTrainOptions()

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

Run MOFA model

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

MOFAobject <- runMOFA(MOFAobject)

lujunyan1118/mofaCLL documentation built on Dec. 21, 2021, 12:42 p.m.