The data set consists of N=200 samples from a cohort of Chronic Lymphocytic Leukemia (CLL) patients where four omics were profiled from blood samples: DNA methylation (450K Illumina microarray), bulk RNA-seq, somatic mutations and drug response data from blood for N=200 patients with. The data set was introduced in detail in this article and can be downloaded here. The MOFA analysis was originally published here
Make sure that you have installed the MOFA2 and the MOFAdata package.
library(data.table)
library(ggplot2)
library(tidyverse)
library(MOFA2)
library(MOFAdata)
library(randomForest)
library(survival)
library(survminer)
Data is stored as a list of matrices. Features are stored in the rows and samples in the columns.
utils::data("CLL_data")
lapply(CLL_data,dim)
## $Drugs
## [1] 310 200
##
## $Methylation
## [1] 4248 200
##
## $mRNA
## [1] 5000 200
##
## $Mutations
## [1] 69 200
The mRNA expression has been normalised by library size, followed by a variance stabilizing transformation using DESeq2:
hist(CLL_data$mRNA)
DNA methylation is calculated for every CpG site using the M-value, which provides a better summary statistic for downstream analysis. For the MOFA analysis we selected the top 1% (N=4248) most variable sites.
hist(CLL_data$Methylation)
In this study the authors have measured the effect of multiple drugs ex vivo using a high-throughput platform. For each drug they have measured 5 different concentrations. The value reported is the viability score (0=all cells died, 1=no cells died).
hist(CLL_data$Drugs)
Mutations are assessed using a panel of common cancer mutations and are summarised in a binary format (0=no mutation, 1=mutation):
table(CLL_data$Mutations)
##
## 0 1
## 8474 667
Load sample metadata as a data.frame. Important columns are:
CLL_metadata <- fread("data/sample_metadata.txt")
head(CLL_metadata)
## sample Gender age TTT TTD treatedAfter died IGHV trisomy12
## 1: H005 m 75.26575 0.57494867 2.625599 TRUE FALSE 1 0
## 2: H006 m NA NA NA NA NA NA NA
## 3: H007 f NA NA NA NA NA NA NA
## 4: H008 m NA NA NA NA NA NA NA
## 5: H010 f 72.78082 2.93223819 2.932238 FALSE FALSE 0 0
## 6: H011 f 72.99452 0.01916496 2.951403 TRUE FALSE 1 0
Create the MOFA object.
MOFAobject <- create_mofa(CLL_data)
MOFAobject
## Untrained MOFA model with the following characteristics:
## Number of views: 4
## Views names: Drugs Methylation mRNA Mutations
## Number of features (per view): 310 4248 5000 69
## Number of groups: 1
## Groups names: group1
## Number of samples (per group): 200
##
Visualise the number of views (rows) and the number of groups (columns), what are their corresponding dimensionalities, and how many missing samples they have (grey bars).
plot_data_overview(MOFAobject)
Two important options:
model_opts <- get_default_model_options(MOFAobject)
model_opts$likelihoods["Mutations"] <- "bernoulli"
model_opts$num_factors <- 15
model_opts
## $likelihoods
## Drugs Methylation mRNA Mutations
## "gaussian" "gaussian" "gaussian" "bernoulli"
##
## $num_factors
## [1] 15
##
## $spikeslab_factors
## [1] FALSE
##
## $spikeslab_weights
## [1] FALSE
##
## $ard_factors
## [1] FALSE
##
## $ard_weights
## [1] TRUE
Prepare the MOFA object.
MOFAobject <- prepare_mofa(MOFAobject,
model_options = model_opts
)
Train the model: this should take ~5min, but DO NOT RUN THIS. The software has evolved since the original publication and the results will not be 100% identical to the original publication. Please load the pretrained model below for reproducibility purposes.
MOFAobject <- run_mofa(MOFAobject)
Load precomputed model.
MOFAobject <- readRDS("data/MOFA2_CLL.rds")
The sample metadata must be provided as a data.frame and it must
contain a column sample
with the sample IDs. Make sure that
the samples in the metadata match the samples in the model.
stopifnot(CLL_metadata$sample %in% samples_metadata(MOFAobject)$sample)
samples_metadata(MOFAobject) <- CLL_metadata
In this step, we load metadata from additional files and use it to replace the uninformative drug and gene IDs with descriptive names.
Keep the model with the original variable names for the gene set enrichment analysis section.
MOFAobject.ensembl <- MOFAobject
updated_features_names <- features_names(MOFAobject)
# Rename drug IDs (i.e. D_001) to drug names (i.e. navitoclax)
drug_metadata <- fread("data/drugs.txt.gz")
tmp <- drug_metadata$name; names(tmp) <- drug_metadata$drug_id
updated_features_names[["Drugs"]] <- stringr::str_replace_all(features_names(MOFAobject)[["Drugs"]], tmp)
# Rename mRNA from ENSEMBLE IDs (i.e. ENSG00000223972) to gene names (i.e. DDX11L1)
gene_metadata <- fread("data/Hsapiens_genes_BioMart.87.txt.gz")
gene_metadata[,symbol:=ifelse(symbol=="",ens_id,symbol)]
tmp <- gene_metadata$symbol; names(tmp) <- gene_metadata$ens_id
# avoid duplicated names with the Mutations view
tmp[tmp%in%features_names(MOFAobject)[["Mutations"]]] <- paste0(tmp[tmp%in%features_names(MOFAobject)[["Mutations"]]],"_mRNA")
updated_features_names[["mRNA"]] <- stringr::str_replace_all(features_names(MOFAobject)[["mRNA"]], tmp)
# Update features names in model
features_names(MOFAobject) <- updated_features_names
The most important insight that MOFA generates is the variance decomposition analysis. This plot shows the percentage of variance explained by each factor across each data modality.
plot_variance_explained(MOFAobject, max_r2=10)
What insights from the data can we learn just from inspecting this plot?
(Q) Based on the MOFA output, if you were to profile just one molecular layer, which one would you choose to maximise the amount of sources of variation captured?
There are a few systematic strategies to characterise the molecular signal that underlies each MOFA Factor and to relate them to existent sample covariates:
Let’s test the association between the MOFA factors versus Gender and age:
correlate_factors_with_covariates(MOFAobject,
covariates = c("Gender","age","died"),
plot = "log_pval"
)
## Warning in correlate_factors_with_covariates(MOFAobject, covariates =
## c("Gender", : There are non-numeric values in the covariates data.frame,
## converting to numeric...
Most Factors don’t have a clear association with any of the covariates. However Factor 1 and Factor 14 appear to have some association with survival outcome. We will explore association with clinical measurements later in the vignette.
How do we interpret the factor values?
Each factor captures a different source of variability in the data.
Mathematically, each Factor is defined by a linear combination of the
input features. Each Factor ordinates cells along a one-dimensional axis
that is centered at zero. Samples with different signs manifest opposite
phenotypes along the inferred axis of variation, with higher absolute
value indicating a stronger effect.
Note that the interpretation of MOFA factors is analogous to the
interpretation of the principal components in PCA.
plot_factors(MOFAobject,
factors = c(1,3),
dot_size = 2.5
)
How do we interpret the feature weights?
The weights provide a score for each feature on each factor. Features
with no association with the corresponding factor are expected to have
values close to zero, whereas features with strong association with the
factor are expected to have large absolute values. The sign of the
weights indicates the direction of the effect: a positive weight
indicates that the feature has higher levels in the cells with positive
factor values, and vice-versa.
By looking at the variance explained plot, we saw that Factor 1 captures variation in all data modalities. Out of all omics, the somatic mutation data is a good place to start, as somatic mutations are very sparse, easy to interpret and any change in the DNA is likely to have downstream consequences to all other molecular layers. Let’s plot the weights:
plot_weights(MOFAobject,
view = "Mutations",
factor = 1,
nfeatures = 10, # Top number of features to highlight
scale = T # Scale weights from -1 to 1
)
Notice that most features lie at zero, indicating that most features have no association with Factor 1. There is however one gene that clearly stands out: IGHV (immunoglobulin heavy chain variable region). This is the main clinical marker for CLL.
An alternative visualisation to the full distribution of weights is to do a line plot that displays only the top features with the corresponding weight sign on the right:
plot_top_weights(MOFAobject,
view = "Mutations",
factor = 1,
nfeatures = 10, # Top number of features to highlight
scale = TRUE # Scale weights from -1 to 1
)
IGHV has a positve weight. This means that samples with positive Factor 1 values have the IGHV mutation whereas samples with negative Factor 1 values do not have the IGHV mutation. To confirm this, let’s plot the Factor values and colour the IGHV mutation status.
plot_factor(MOFAobject,
factors = 1,
color_by = "IGHV",
add_violin = TRUE,
dodge = TRUE,
show_missing = FALSE
)
We can also plot Factor values coloured by other covariates, for
example Gender
. As concluded from the association analysis
above, this variable has no association with Factor 1:
plot_factor(MOFAobject,
factors = 1,
color_by = "Gender",
dodge = TRUE,
add_violin = TRUE
)
From the variance explained plot we know that Factor 1 drives variation across all data modalities. Let’s visualise the mRNA expression changes that are associated with Factor 1:
plot_weights(MOFAobject,
view = "mRNA",
factor = 1,
nfeatures = 10
)
In this case we have a large amount of genes with high positive and
high negative weights. Genes with large positive values will be more
expressed in the samples with IGHV mutation, whereas genes with large
negative values will be more expressed in the samples without the IGHV
mutation. Let’s verify this. The function plot_data_scatter
generates a scatterplot of Factor 1 values (x-axis) versus expression
values (y-axis) for the top 4 genes with largest positive weight.
Samples are coloured by IGHV status:
plot_data_scatter(MOFAobject,
view = "mRNA",
factor = 1,
features = 4,
sign = "positive",
color_by = "IGHV"
) + labs(y="RNA expression")
This function generates a scatterplot of Factor 1 values (x-axis) versus expression values (y-axis) for the top 4 genes with largest negative weight. Samples are coloured by IGHV status:
plot_data_scatter(MOFAobject,
view = "mRNA",
factor = 1,
features = 4,
sign = "negative",
color_by = "IGHV"
) + labs(y="RNA expression")
An alternative visualisation is to use a heatmap.
plot_data_heatmap(MOFAobject,
view = "mRNA",
factor = 1,
features = 25,
cluster_rows = FALSE, cluster_cols = FALSE,
show_rownames = TRUE, show_colnames = FALSE,
scale = "row"
)
(Q) Can you suggest new RNA expression and DNA methylation
markers for personalised treatment recommendations according to Factor 1
(the IGHV status)?
First explore the MOFA weights, then go back to the input data and do
boxplots for the chosen markers (x-axis being the IGHV status and y-axis
being the marker’s expression or DNA methylation values).
Hints:
- The section Customized analysis may
be helpful to extract the weights and the data in a long data.frame
format - the IGHV status for each sample can be fetched from the
CLL_metadata
object - the ggpubr
package provides very useful ggplot-based visualisations, including boxplots
with p-values. I highly recommend it!
(Q) Your task is to provide a characterisation for Factor 3.
Try a similar pipeline as for Factor 1 and answer the following
questions:
- Which mutation underlies Factor 3?
- Can you identify mRNA markers?
- Do a (small) bibliographical search to check if your predictions make
sense
Now that we have characterised the etiology of the two main Factors, let’s do a scatterplot colouring each patient by their somatic mutation profile:
p <- plot_factors(MOFAobject,
factors = c(1,3),
color_by = "IGHV",
shape_by = "trisomy12",
dot_size = 2.5,
show_missing = T
)
p <- p +
geom_hline(yintercept=-1, linetype="dashed") +
geom_vline(xintercept=(-0.5), linetype="dashed")
print(p)
## Warning: Removed 27 rows containing missing values (`geom_point()`).
The scatterplot of Factor 1 vs Factor 3 reveals that a few samples are missing the somatic mutation status. In this case, the doctors were not able to classify patients into their clinical subgroups. But we can now use MOFA to exploit the molecular profiles and attempt to impute the IGHV and trisomy12 status.
# Prepare data
df <- as.data.frame(get_factors(MOFAobject, factors=c(1,2))[[1]])
# Train the model for IGHV after removing missing observations
df$IGHV <- as.factor(samples_metadata(MOFAobject)$IGHV)
model.ighv <- randomForest(IGHV ~ ., data=df[!is.na(df$IGHV),])
# Do predictions
samples_metadata(MOFAobject)$IGHV.pred <- stats::predict(model.ighv, df)
# Prepare data
df <- as.data.frame(get_factors(MOFAobject, factors=c(1,2))[[1]])
# Train the model for Trisomy12 after removing missing observations
df$trisomy12 <- as.factor(samples_metadata(MOFAobject)$trisomy12)
model.trisomy12 <- randomForest(trisomy12 ~ ., data=df[!is.na(df$trisomy12),])
samples_metadata(MOFAobject)$trisomy12.pred <- stats::predict(model.trisomy12, df)
Plot predictions for IGHV.
samples_metadata(MOFAobject)$IGHV.pred_logical <- c("True","Predicted")[as.numeric(is.na(samples_metadata(MOFAobject)$IGHV))+1]
p <- plot_factors(MOFAobject,
factors = c(1,3),
color_by = "IGHV.pred",
shape_by = "IGHV.pred_logical",
dot_size = 2.5,
show_missing = T
)
p <- p +
geom_hline(yintercept=-1, linetype="dashed") +
geom_vline(xintercept=(-0.5), linetype="dashed")
print(p)
In addition to exploring the individual weights for each factor, we can use enrichment analysis to look for signiificant associations of factors to genesets. Here, we use the Reactome genesets for illustrations, which is contained in the MOFAdata package. For more details on how the GSEA works we encourage the users to read the GSEA vignette
Gene set annotations are provided as a binary membership matrix. Genes are stored in the rows, pathways are stored in the columns. A value of 1 indicates that gene \(j\) belongs to the pathway \(i\).
utils::data(reactomeGS) # from MOFAdata
head(colnames(reactomeGS))
## [1] "ENSG00000187634" "ENSG00000188976" "ENSG00000187961" "ENSG00000187583"
## [5] "ENSG00000187642" "ENSG00000188290"
head(rownames(reactomeGS))
## [1] "Interleukin-6 signaling"
## [2] "Apoptosis"
## [3] "Hemostasis"
## [4] "Intrinsic Pathway for Apoptosis"
## [5] "Cleavage of Growing Transcript in the Termination Region "
## [6] "PKB-mediated events"
These are the steps for doing Gene Set Enrichment Analysis (GSEA) with MOFA:
j
belongs to
pathway i
. A value of 0 indicates elsewise.mean.diff
(difference in the average weight
between foreground and background genes) or rank.sum
(difference in the sum of ranks between foreground and background
genes).parametric
(a simple and very liberal parametric t-test),
cor.adj.parametric
(parametric t-test adjusted by the
correlation between features), permutation
(unparametric,
the null distribution is created by permuting the weights. This option
is computationally expensive, but it preserves the correlation structure
between features in the data.)enrichment.results <- run_enrichment(
object = MOFAobject.ensembl,
view = "mRNA",
feature.sets = reactomeGS,
set.statistic = "mean.diff",
statistical.test = "parametric"
)
The enrichment analysis returns a list of 5 elements:
names(enrichment.results)
## [1] "feature.sets" "pval" "pval.adj"
## [4] "feature.statistics" "set.statistics" "sigPathways"
Plot an overview of the number of significant pathways per
factor.
It seems that most of the Factors do not have clear gene set signatures.
An exception is Factor 8, which has a very strong enrichment for genes
with positive weights.
plot_enrichment_heatmap(enrichment.results)
(Q) Can you characterise Factor 8 based on the GSEA results?
Which genes are driving the top enriched pathways?
Hint: use the function plot_enrichment
.
(Q) Which drugs are associated with Factor 8? What is their
target pathway? Do they make biological sense?
Hint: use the drug_metadata
object.
For customized exploration of weights and factors, you can directly
fetch the variables from the model using get_weights
and
get_factors
:
weights <- get_weights(MOFAobject,
views = "all",
factors = "all",
as.data.frame = TRUE
)
head(weights)
## feature factor value view
## 1 navitoclax_1 Factor1 0.0005978113 Drugs
## 2 navitoclax_2 Factor1 -0.0028452079 Drugs
## 3 navitoclax_3 Factor1 -0.0294330772 Drugs
## 4 navitoclax_4 Factor1 -0.0268737264 Drugs
## 5 navitoclax_5 Factor1 -0.0162836239 Drugs
## 6 ibrutinib_1 Factor1 0.0014393722 Drugs
factors <- get_factors(MOFAobject,
factors = "all",
as.data.frame = TRUE
)
head(factors)
## sample factor value group
## 1 H045 Factor1 -2.287864 group1
## 2 H109 Factor1 -2.253917 group1
## 3 H024 Factor1 2.537418 group1
## 4 H056 Factor1 2.012213 group1
## 5 H079 Factor1 -1.965762 group1
## 6 H164 Factor1 -2.562547 group1
The factors inferred by MOFA can be related to clinical outcomes such as time to treatment or survival times. As this type of data is censored (not for all samples we have already observed the event) we will use Cox models for this purpose. In a Cox proportional hazards model we model the hazard of an event occurring (e.g. death or treatment) as a function of some covariates (here the factors). If a factor has a influence on the survival time or time to treatment it will receive a high absolute coefficient in the Cox model. In particular:
To fit these models we will use the coxph
function in
the survival
package. The survival data is stored in a
survival object that contains both the time a sample has been followed
up and whether the event has occurred (as 0,1).
Let’s take time to treatment as an example here. The sample metadata
contains the follow-up times per sample in years in the column
TTT
, and the column treatedAfter
indicated
whether a treatment occurred.
SurvObject <- Surv(samples_metadata(MOFAobject)$TTT, samples_metadata(MOFAobject)$treatedAfter)
Z <- get_factors(MOFAobject)[[1]]
fit <- coxph(SurvObject ~ Z)
fit
## Call:
## coxph(formula = SurvObject ~ Z)
##
## coef exp(coef) se(coef) z p
## ZFactor1 -0.35121 0.70383 0.06753 -5.201 1.99e-07
## ZFactor2 0.53641 1.70986 0.14997 3.577 0.000348
## ZFactor3 0.08833 1.09235 0.06681 1.322 0.186124
## ZFactor4 -0.31316 0.73113 0.10120 -3.095 0.001971
## ZFactor5 -0.10717 0.89838 0.08275 -1.295 0.195279
## ZFactor6 0.01357 1.01366 0.09228 0.147 0.883122
## ZFactor7 0.24953 1.28343 0.08331 2.995 0.002742
## ZFactor8 0.43294 1.54179 0.10291 4.207 2.59e-05
## ZFactor9 0.17822 1.19508 0.09171 1.943 0.051981
## ZFactor10 0.29090 1.33763 0.10322 2.818 0.004830
## ZFactor11 -0.52174 0.59348 0.10001 -5.217 1.82e-07
## ZFactor12 0.16270 1.17668 0.11105 1.465 0.142895
## ZFactor13 -0.33730 0.71369 0.09403 -3.587 0.000334
## ZFactor14 0.04600 1.04708 0.07343 0.627 0.530968
## ZFactor15 -0.18496 0.83114 0.10575 -1.749 0.080301
##
## Likelihood ratio test=112.9 on 15 df, p=< 2.2e-16
## n= 174, number of events= 96
## (26 observations deleted due to missingness)
We can see that several factors have a significant association to time to treatment. For example, Factor 1 has a negative coefficient. Samples with low factor values have an increased hazard compared to samples with a large factor values.
(Q) Which Factors are associated with the clinical covariate (time to next treatment)?
Extract p-values and Cox model coefficients (i.e. hazard ratios).
s <- summary(fit)
coef <- s[["coefficients"]]
df <- data.frame(
factor = factor(rownames(coef), levels = rev(rownames(coef))),
p = coef[,"Pr(>|z|)"],
coef = coef[,"exp(coef)"],
lower = s[["conf.int"]][,"lower .95"],
higher = s[["conf.int"]][,"upper .95"]
)
Plot the Hazard ratio per factor, together with 95% confidence intervals.
ggplot(df, aes(x=factor, y=coef, ymin=lower, ymax=higher, color=p<0.01)) +
geom_pointrange() +
coord_flip() +
scale_x_discrete() +
labs(y="Hazard Ratio", x="") +
geom_hline(aes(yintercept=1), linetype="dotted") +
theme_bw()
By inspecting the variance explained plot, we can see that the RNA expression is capturing most of the sources of variation in this data set. There are only a few Factors that cannot be captured using RNA expression ( for example Factors 2 and 14). The ex vivo drug response assay also captures a lot of variability, but it is much harder to obtain from patient cohorts than RNA expression data. The other two data modalities are less informative: DNA methylation data is only active in Factor 1, 7, and 9, and Somatic mutations are only associated with Factor 1.
If were were to profile just one molecular layer for a large number of patients in a cost-effective way, we would need to compare the feasibility and costs of Drug response assays and RNA sequencing. The latter is much cheaper and more standardized in the community.
We first collect the genes with the largest weights for Factor 1. Then we can do boxplots stratifying samples by IGHV status, followed by statistical testing with the null hypothesis that the average gene expression does not differ between groups (a simple t-test should work for a first exploration).
Extract mRNA weights from the MOFA object.
rna.weights <- get_weights(MOFAobject,
views = "mRNA",
factors = 1,
abs = TRUE, # we do not distinguish between direction of effect
as.data.frame = TRUE
)
# Extract top N genes
top.genes <- rna.weights %>%
.[order(rna.weights$value, decreasing = T),] %>%
head(n=9) %>% .$feature %>% as.character
top.genes
## [1] "ZNF667" "SOWAHC" "SEPT10" "ZNF471" "LPL" "PLD1" "ADAM29"
## [8] "LDOC1" "L3MBTL4"
Fetch mRNA data from the MOFAobject for the top genes.
rna.data <- get_data(MOFAobject,
views = "mRNA",
as.data.frame = TRUE,
features = list("mRNA"=top.genes)
)
head(rna.data)
## view group feature sample value
## 1 mRNA group1 ZNF667 H045 10.323723
## 2 mRNA group1 SOWAHC H045 10.134042
## 3 mRNA group1 SEPT10 H045 11.473791
## 4 mRNA group1 ZNF471 H045 9.680574
## 5 mRNA group1 LPL H045 12.686457
## 6 mRNA group1 PLD1 H045 8.453543
Add IGHV status from the sample metadata.
to.plot <- rna.data %>%
merge(CLL_metadata[,c("sample","IGHV")], by="sample")
colnames(to.plot)
## [1] "sample" "view" "group" "feature" "value" "IGHV"
(Optional) Remove samples with unknown IGHV status.
to.plot <- to.plot[!is.na(to.plot$IGHV),]
Box plots with statistical comparison of means.
ggpubr::ggboxplot(to.plot, x = "IGHV", y = "value", fill = "IGHV", facet="feature") +
stat_compare_means() +
labs(x="IGHV status", y="mRNA expression") +
theme(legend.position="none")
Plotting the GSEA results for Factor 4 reveals that this Factor is capturing differences in the stress response of the blood cells. We have significantly enriched pathways such as cellular response to stress or heat shock response.
plot_enrichment(enrichment.results, factor = 4, max.pathways = 15)
The top genes that are driving this pathway are Heat Shock Proteins and some inflammatory markers such as TNF.
plot_top_weights(MOFAobject,
view = "mRNA",
factor = 4,
nfeatures = 15
)
Almost all of these genes have a negative weight, with the exception of EML6, which means that they have higher levels in the samples with negative Factor 4 values:
plot_data_scatter(MOFAobject,
factor = 4,
view = "mRNA",
features = 6,
add_lm = TRUE
) + labs(y="RNA expression")
plot_top_weights(MOFAobject,
view = "Drugs",
factor = 4,
nfeatures = 15
)
plot_data_scatter(MOFAobject,
view = "Drugs",
factor = 4,
features = 6,
add_lm = TRUE
) + labs(y="Viability")
Out of the 5 top drugs, the target category for 3 of them is Reactive Oxygen Species, which are closely related to stress response mechanisms.
drug_metadata[grep("SD51|BAY|SD07|MIS-43|doxorubicine",drug_metadata$name),]
## drug_id name main_targets
## 1: D_041 BAY 11-7085 NFkB
## 2: D_127 SD07 ROS
## 3: D_141 SD51 ROS
## 4: D_149 MIS-43 ROS
## 5: D_159 doxorubicine DNA intercalation, Topoisomerase II
## target_category pathway
## 1: NFkB Other
## 2: Reactive oxygen species Reactive oxygen species
## 3: Reactive oxygen species Reactive oxygen species
## 4: Reactive oxygen species Reactive oxygen species
## 5: DNA damage response Other
sessionInfo()
## R version 4.3.1 (2023-06-16 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] survminer_0.4.9 ggpubr_0.6.0 survival_3.5-5
## [4] randomForest_4.7-1.1 MOFAdata_1.16.1 MOFA2_1.10.0
## [7] lubridate_1.9.3 forcats_1.0.0 stringr_1.5.0
## [10] dplyr_1.1.3 purrr_1.0.2 readr_2.1.4
## [13] tidyr_1.3.0 tibble_3.2.1 tidyverse_2.0.0
## [16] ggplot2_3.4.4 data.table_1.14.8
##
## loaded via a namespace (and not attached):
## [1] mnormt_2.1.1 gridExtra_2.3 rlang_1.1.1
## [4] magrittr_2.0.3 matrixStats_1.0.0 compiler_4.3.1
## [7] mgcv_1.8-42 dir.expiry_1.8.0 png_0.1-8
## [10] vctrs_0.6.4 reshape2_1.4.4 pkgconfig_2.0.3
## [13] crayon_1.5.2 fastmap_1.1.1 backports_1.4.1
## [16] labeling_0.4.3 KMsurv_0.1-5 utf8_1.2.4
## [19] rmarkdown_2.25 tzdb_0.4.0 xfun_0.40
## [22] cachem_1.0.8 jsonlite_1.8.7 rhdf5filters_1.12.1
## [25] DelayedArray_0.26.7 Rhdf5lib_1.22.1 psych_2.3.9
## [28] broom_1.0.5 parallel_4.3.1 R6_2.5.1
## [31] bslib_0.5.1 stringi_1.7.12 RColorBrewer_1.1-3
## [34] reticulate_1.32.0.9001 car_3.1-2 jquerylib_0.1.4
## [37] Rcpp_1.0.11 knitr_1.44 zoo_1.8-12
## [40] R.utils_2.12.2 IRanges_2.34.1 Matrix_1.6-1.1
## [43] splines_4.3.1 timechange_0.2.0 tidyselect_1.2.0
## [46] rstudioapi_0.15.0 abind_1.4-5 yaml_2.3.7
## [49] lattice_0.21-8 plyr_1.8.9 basilisk.utils_1.12.1
## [52] withr_2.5.1 evaluate_0.22 Rtsne_0.16
## [55] survMisc_0.5.6 pillar_1.9.0 filelock_1.0.2
## [58] MatrixGenerics_1.12.3 carData_3.0-5 corrplot_0.92
## [61] stats4_4.3.1 generics_0.1.3 S4Vectors_0.38.1
## [64] hms_1.1.3 munsell_0.5.0 scales_1.2.1
## [67] xtable_1.8-4 glue_1.6.2 pheatmap_1.0.12
## [70] tools_4.3.1 ggsignif_0.6.4 cowplot_1.1.1
## [73] rhdf5_2.44.0 grid_4.3.1 colorspace_2.1-0
## [76] nlme_3.1-162 basilisk_1.12.1 HDF5Array_1.28.1
## [79] cli_3.6.1 km.ci_0.5-6 fansi_1.0.5
## [82] S4Arrays_1.0.6 uwot_0.1.16 gtable_0.3.4
## [85] R.methodsS3_1.8.2 rstatix_0.7.2 sass_0.4.7
## [88] digest_0.6.33 BiocGenerics_0.46.0 ggrepel_0.9.3
## [91] farver_2.1.1 htmltools_0.5.6.1 R.oo_1.25.0
## [94] lifecycle_1.0.3