orthos 1.2.0
RNAseq-based differential expression analysis upon cellular perturbations, such as gene knockouts, RNA knockdowns or compound treatment experiments, is the most commonly used tool for probing molecular mechanisms of action due to its simplicity and low cost.
However, interpretation of such gene expression contrasts is confounded by the complex and nuanced impacts of experimental treatments on cellular processes.
For example, knockout or over-expression of a transcription factor will not only alter the transcription of its direct target genes, but also cause many secondary expression changes. In addition, treatments or treatment delivery agents typically elicit a variety of unintended, systemic responses (such as immune, toxic, metabolic) that cannot be well-controlled for by the design of the study.
The final experimentally measured gene expression changes are a hard to assess convolution of specific and non-specific secondary and lateral treatment effects.
orthos
is a generative modelling-based approach that disentangles the
experiment-specific from the non-specific effects of perturbations on gene
expression. It is trained on a large corpus of gene expression contrasts
(per organism >60K annotated, >0.5M augmented), compiled from the ARCHS4
database of uniformly processed RNAseq experiments (Lachmann et al. (2018)).
It accurately captures and isolates non-specific effects (effects
that are observed across multiple treatments) while accounting for context
(tissue or cell-line experimental background).
The residual specific component obtained from this decomposition offers a more unequivocal experimental signature and is more closely related to the direct molecular effects of the perturbation when compared to the raw signal.
In addition to providing a clearer understanding of the effects of experimental
treatments on gene expression, orthos
also enables researchers to query the
contrast database with arbitrary contrasts and identify experiments with similar
specific effects, ultimately helping to map treatments to mechanisms of action.
orthos
can be installed from from Bioconductor using BiocManager::install()
:
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("orthos")
# or also...
BiocManager::install("orthos", dependencies = TRUE)
After installation, the package can be loaded with:
library(orthos)
A typical analysis involves two steps:
Decomposing one or several contrasts into their corresponding specific and
non-specific components using the decomposeVar()
function and
Performing queries with the original and decomposed specific and non-specific
contrasts against the contrast database using the queryWithContrasts()
function.
To demonstrate the functionality of orthos
we use a dataset from the the GEO
series GSE215150.
This series was not part of the orthos
training or the orthosData
contrast
database; it was only publicly released on January 1st 2023 after the freeze of
the training data to the ARCHS4 v2.1.2 database.
The performed experiment involves over-expression (OE) of the MKL/megakaryoblastic leukemia 1 gene (also termed MRTFA/myocardin related transcription factor A) and a constitutively active mutant MKL1 (caMKL1, described in Hu et al. (2019)). Both OE experiments were performed in mouse LM2 and human 4T1 tumor-derived breast cancer cell lines. In addition to the MKL1/caMKL1 OE samples, the series also contains no-treatment controls for each of the two cell lines.
For simplicity the three biological replicates of each available condition have been collapsed in the data provided in the package.
In the provided form each of the two datasets (Mouse, Human) contains raw counts for over 55,000 genes identified by gene symbols in three conditions: Control (Ctrl), MKL1 OE (MKL1) and constitutively-active MKL1 OE (caMKL1).
Load the human dataset:
MKL1_human <- readRDS(system.file("extdata", "GSE215150_MKL1_Human.rds",
package = "orthos"))
head(MKL1_human)
#> Ctrl MKL1 caMKL1
#> DDX11L1 6 3 4
#> WASH7P 47 57 41
#> MIR6859-1 10 10 6
#> MIR1302-2HG 0 6 1
#> MIR1302-2 0 0 0
#> FAM138A 1 1 5
dim(MKL1_human)
#> [1] 59453 3
Load the mouse dataset:
MKL1_mouse <- readRDS(system.file("extdata", "GSE215150_MKL1_Mouse.rds",
package = "orthos"))
head(MKL1_mouse)
#> Ctrl MKL1 caMKL1
#> 4933401J01Rik 0 0 0
#> Gm26206 0 0 0
#> Xkr4 0 0 0
#> Gm18956 0 0 0
#> Gm37180 0 0 0
#> Gm37363 0 0 0
dim(MKL1_mouse)
#> [1] 55367 3
decomposeVar()
orthos
modelsThe workhorse behind orthos
are organism-specific conditional variational
autoencoder (cVAE) models that break down the variance of a given differential
expression experiment into a non-specific and an experiment-specific
component.
The non-specific component corresponds to gene variance that has been observed across multiple instances during training, while the experiment-specific variance is fairly unique to the experiment.
The inputs to the models are gene counts in the form of log2-transformed counts per million (LCPMs) that are used to encode the context of the performed experiment as well as the actual gene expression contrasts in the form of log2 fold-changes (LFCs), i.e log2-transformed CPM ratios.
As we will see, calculation of those inputs is by default performed internally given only raw gene counts and a specification of the contrasted conditions.
Given these inputs, the model maps the contrast to a concise latent representation (zD) which retains its recurring -and therefore compressible- traits. The compressed latent representation is then used to reconstruct a decoded version of the contrast*.
The decoded contrast corresponds directly to the non-specific component of the observed effects; it subsumes regularities i.e gene variance that the model can account for because it has been repeatedly encountered, in some form, during training.
The residual obtained after removing the decoded contrast from the original one is the specific component; this encompasses the gene variance that the model cannot account for (experiment-specific biological effects + random noise).
From now on the terms decoded and non-specific will be used in conjunction or interchangeably. Ditto for the terms residual and specific.
* Notice that both the latent encoding and the decoded output are conditioned on context
(i.e they are context-specific). This means that decomposing a contrast on a different context
will produce a different output. An interesting ancillary application of this conditioning is that
one can “morph” decoded contrasts to in-silico evaluate non-specific effects in new contexts.
In essence, we can infer what the non-specific effects would look like had the experiment been performed in
e.g different cell-lines, tissues, batches, patients or using different library-preparation protocols.
The inferred decodings will reflect the new context in multiple ways. For example, genes that
were non-detected/non-expressed in the original contrast (and therefore had neutral LFCs )
will produce new (possibly non-neutral) decoded LFCs if present in the new context and vice versa.
This “out-of-context” type of inference is limited to the decoded contrast, as by definition, the residual is not part of
the model’s generative capacity.
Mechanically, out-of-context inference of non-specific effects is a simple as evaluating the same contrast using multiple contexts when
calling
decomposeVar()
(see section below).
decomposeVar()
The function that performs this contrast decomposition into the non-specific and
specific components is decomposeVar()
.
There are two available modes in which the required inputs can be fed into the
function:
In the first mode the user passes the matrix M
of raw counts (genes in rows,
conditions in columns) and two vectors treatm
and cntr
of equal length
specifying the column indices in M
corresponding to the “treatments” and
their respective “controls”. The same column indices can be repeated multiple
times in these vectors, for example in the case where multiple treatments are
paired to the same control.
In the second mode the user passes the matrix M
of raw counts and a second
matrix MD
that contains pre-calculated log2 fold-changes for the contrasts to
be analyzed. In this mode, M
specifies the contexts for the corresponding
columns of MD
, and thus the two matrices need to have the same
dimensionality and identical row- and column-names.
This would be the mode of choice if e.g one wishes to produce the LFCs independently
or if one wants to evaluate the decoding of the same contrast(s) in multiple contexts
(e.g for “out-of-context” inference of non-specific effects described above). In the
latter case, copies of the same contrast in columns of MD
will be paired with
columns of M
specifying the different contexts.
In both modes the rownames of M
(and MD
if specified) need to correspond to
valid gene identifiers (currently orthos
supports Entrez gene identifiers,
ENSEMBL gene identifiers, gene symbols or ARCHS4 gene identifiers).
By default the type of gene identifier is detected automatically.
The first time that decomposeVar
is executed for a particular organism, the
models required for inference will be automatically downloaded from
ExperimentHub
and cached in the user ExperimentHub directory
(see ExperimentHub::getExperimentHubOption("CACHE")
) using the orthosData
companion package.
For the MKL1 data, which are stored as raw counts, it is more natural to
call decomposeVar()
using the first mode.
Decomposing the human contrasts:
#Decompose MKL1-vs-Cntrl and caMKL1-vs-Cntrl contrasts for human:
dec_MKL1_human <- decomposeVar(M = MKL1_human, treatm = c(2, 3), cntr = c(1, 1),
organism = "Human", verbose = FALSE)
#> 1/1 - 0s - 422ms/epoch - 422ms/step
#> 1/1 - 0s - 125ms/epoch - 125ms/step
#> 1/1 - 0s - 79ms/epoch - 79ms/step
dec_MKL1_human
#> class: SummarizedExperiment
#> dim: 20411 2
#> metadata(0):
#> assays(4): INPUT_CONTRASTS DECODED_CONTRASTS RESIDUAL_CONTRASTS CONTEXT
#> rownames(20411): A1BG A1BG-AS1 ... WFDC5 XAGE2
#> rowData names(1): User_provided_IDs
#> colnames(2): MKL1 caMKL1
#> colData names(3): ACCOUNTED_VARIANCE.DECODED
#> ACCOUNTED_VARIANCE.RESIDUAL ACCOUNTED_VARIANCE.COMMON
Decomposing the mouse contrasts:
#Decompose MKL1-vs-Cntrl and caMKL1-vs-Cntrl contrasts for mouse:
dec_MKL1_mouse <- decomposeVar(M = MKL1_mouse, treatm = c(2, 3), cntr = c(1, 1),
organism = "Mouse", verbose = FALSE)
#> 1/1 - 0s - 94ms/epoch - 94ms/step
#> 1/1 - 0s - 110ms/epoch - 110ms/step
#> 1/1 - 0s - 94ms/epoch - 94ms/step
dec_MKL1_mouse
#> class: SummarizedExperiment
#> dim: 20339 2
#> metadata(0):
#> assays(4): INPUT_CONTRASTS DECODED_CONTRASTS RESIDUAL_CONTRASTS CONTEXT
#> rownames(20339): 0610005C13Rik 0610009B22Rik ... n-R5s146 n-R5s149
#> rowData names(1): User_provided_IDs
#> colnames(2): MKL1 caMKL1
#> colData names(3): ACCOUNTED_VARIANCE.DECODED
#> ACCOUNTED_VARIANCE.RESIDUAL ACCOUNTED_VARIANCE.COMMON
The output of decomposeVar()
is a SummarizedExperiment
object with
dimensions N
x M
, where N
is the number of orthos
genes* for that organism and M
is the number of contrasts specified during input.
The SummarizedExperiment
output also has 4 assay slots corresponding to
the input contrasts, decoded (non-specific), and residual (specific) components,
as well as the gene context.
Contrasts are represented as log2 fold-changes (LFCs) and context is represented
as log2-transformed counts per million (log2 CPM).
We can use the returned object to produce an MA plot for the original contrast or to check how the input and decomposed contrasts are related to each other.
For example for the mouse caMKL1 contrast:
suppressPackageStartupMessages({
library(ggplot2)
library(SummarizedExperiment)
})
assays(dec_MKL1_mouse)
#> List of length 4
#> names(4): INPUT_CONTRASTS DECODED_CONTRASTS RESIDUAL_CONTRASTS CONTEXT
#MA plot of for the input contrasts:
DF <- data.frame(L2CPM= assay(dec_MKL1_mouse,"CONTEXT")[,2],
L2FC_INPUT=assay(dec_MKL1_mouse,"INPUT_CONTRASTS")[,2],
L2FC_DECODED=assay(dec_MKL1_mouse,"DECODED_CONTRASTS")[,2],
L2FC_RESIDUAL=assay(dec_MKL1_mouse,"RESIDUAL_CONTRASTS")[,2]
)
#MA plot of for the input contrast
P1 <- ggplot(data=DF, aes(x=L2CPM, y=L2FC_INPUT)) +
geom_point(alpha=0.4, size=1.8) +
geom_hline(aes(yintercept = 0), colour = "darkgray", linewidth = 0.5) +
xlab("Expression (Log2 CPMs)") +
ylab("Log2 Fold Change")
#Delta-delta plots for the input and decomposed contrast fractions
P2 <- ggplot(data=DF, aes(x=L2FC_INPUT, y=L2FC_DECODED)) +
geom_point(alpha=0.4, size=1.8) +
geom_hline(aes(yintercept = 0), colour = "darkgray", linewidth = 0.5) +
xlab("Log2 Fold Change INPUT") +
ylab("Log2 Fold Change DECODED")
P3 <- ggplot(data=DF, aes(x=L2FC_INPUT, y=L2FC_RESIDUAL)) +
geom_point(alpha=0.4, size=1.8) +
geom_hline(aes(yintercept = 0), colour = "darkgray", linewidth = 0.5) +
xlab("Log2 Fold Change INPUT") +
ylab("Log2 Fold Change RESIDUAL")
P4 <- ggplot(data=DF, aes(x=L2FC_DECODED, y=L2FC_RESIDUAL)) +
geom_point(alpha=0.4, size=1.8) +
geom_hline(aes(yintercept = 0), colour = "darkgray", linewidth = 0.5) +
xlab("Log2 Fold Change DECODED") +
ylab("Log2 Fold Change RESIDUAL")
cowplot::plot_grid(P1,P2,P3,P4)
As expected, both the decoded and residual components are correlated to the input contrast. However, the residual and decoded components are largely uncorrelated.
The colData
of the object summarizes the proportion of variance accounted for
in each decomposed component:
colData(dec_MKL1_human)
#> DataFrame with 2 rows and 3 columns
#> ACCOUNTED_VARIANCE.DECODED ACCOUNTED_VARIANCE.RESIDUAL
#> <numeric> <numeric>
#> MKL1 NA NA
#> caMKL1 NA NA
#> ACCOUNTED_VARIANCE.COMMON
#> <numeric>
#> MKL1 NA
#> caMKL1 NA
colData(dec_MKL1_mouse)
#> DataFrame with 2 rows and 3 columns
#> ACCOUNTED_VARIANCE.DECODED ACCOUNTED_VARIANCE.RESIDUAL
#> <numeric> <numeric>
#> MKL1 NA NA
#> caMKL1 NA NA
#> ACCOUNTED_VARIANCE.COMMON
#> <numeric>
#> MKL1 NA
#> caMKL1 NA
* Notice that, of the total gene features present in the input (over 55,000), only
~20,000 genes are part of the orthos
model () and the decomposeVar()
output.
These ~20,000 orthos
genes are “sanctioned” according to several criteria
(located on canonical chromosomes, no pseudogenes, no ribosomal protein genes,
detected in at least a small fraction of the ARCHS4 database).
The model is highly robust to small fractions of orthos
genes not being part
of the user input, even if those genes are expressed in the context under
consideration. That being noted, it is safer to feed-in inputs that are as
comprehensive as possible, i.e not filtered in any way, in terms of gene
features.
queryWithContrasts()
Typically, the next step of the analysis involves querying the contrasts
database (orthosData
) to identify public experiments similar to the one(s)
under investigation, either in terms of the original or decomposed decoded (non-specific)
and residual (specific) contrasts. As we will see in the following examples the results’
of these queries can guide the interpretation of of the different contrast fractions.
orthosData
contains over 100,000 differential gene expression experiments
compiled from the ARCHS4 database of publicly available expression data (Lachmann et al. (2018)).
Each entry in orthosData
corresponds to a pair of RNAseq samples contrasting
a treatment vs a control condition.
A combination of metadata, semantic and quantitative analyses was used to
determine the proper assignment of samples to such pairs in orthosData
.
The function that performs the queries against orthosData
is
queryWithContrasts()
. The input to this function is the SummarizedExperiment
object obtained in the previous step from decomposeVar()
, either the complete
object or one that has been column-subsetted, allowing to query the
contrast database with only a subset of the decomposed contrasts.
As was the case for the orthos
models, a database will be automatically downloaded from
ExperimentHub
and cached in the user ExperimentHub directory
(see ExperimentHub::getExperimentHubOption("CACHE")
) using the orthosData
companion package, the first time queryWithContrasts()
is called for that
database or the first time the user attempts to access the database directly
with loadContrastDatabase()
(see Accessing the contrast database) .
The queryWithContrasts()
function returns a list with three elements per query contrast:
orthosData
entries for a specific
component (input, decoded/non-specific, residual/specific).orthosData
metadata for each of the
top detailTopn
hits per component (default 10).In the following examples, please note that the queries are run using
mode = "DEMO"
in order to keep computations short. For actual analyses,
the default mode = "ANALYSIS"
should be used.
Examples queries using the decomposed human MKL1 data:
# parallelization parameters:
params <- BiocParallel::MulticoreParam(workers = 2)
# for demonstration purposes (for actual analyses, use 'mode = "ANALYSIS"'):
query.res.human <- queryWithContrasts(dec_MKL1_human, organism = "Human",
BPPARAM = params, verbose = FALSE,
mode = "DEMO")
names(query.res.human)
#> [1] "pearson.rhos" "zscores" "TopHits"
names(query.res.human$zscores)
#> [1] "INPUT_CONTRASTS" "DECODED_CONTRASTS" "RESIDUAL_CONTRASTS"
# query contrasts in rows, `orthosData` entries in columns:
dim(query.res.human$zscores$RESIDUAL_CONTRASTS)
#> [1] 2 1000
summary(t(query.res.human$zscores$RESIDUAL_CONTRASTS))
#> MKL1 caMKL1
#> Min. :-3.67488 Min. :-3.67488
#> 1st Qu.:-0.64686 1st Qu.:-0.64686
#> Median :-0.06019 Median :-0.06019
#> Mean : 0.00000 Mean : 0.00000
#> 3rd Qu.: 0.62312 3rd Qu.: 0.62312
#> Max. : 3.66232 Max. : 3.66232
#Information on the top hits of the query using the residual human MKL1/caMKL1 contrasts:
query.res.human$TopHits$RESIDUAL_CONTRASTS
#> $MKL1
#> DataFrame with 10 rows and 8 columns
#> Zscore TREATM_geo_accession title
#> <numeric> <character> <character>
#> GSM4296923 3.66232 GSM4296923 TEAD4 knock Down HTS..
#> GSM5750634 3.33309 GSM5750634 HLTE_A20_IAV infection
#> GSM5022425 3.20173 GSM5022425 D7_hiPSC1_NAM_2
#> GSM2754776 3.12605 GSM2754776 143_09_0278_DAY_0_T_..
#> GSM3976595 3.07313 GSM3976595 FSK_IL-1beta_3
#> GSM3976593 3.02627 GSM3976593 FSK_IL-1beta_1
#> GSM3976594 3.02548 GSM3976594 FSK_IL-1beta_2
#> GSM3291117 2.97296 GSM3291117 TGFb and JQ1 treated..
#> GSM4451228 2.90588 GSM4451228 22PcEP_FBS_NRG1_1
#> GSM4648524 2.58702 GSM4648524 RNAseq_RBMX_KD2_1
#> characteristics_ch1 series_id source_name_ch1
#> <character> <character> <character>
#> GSM4296923 tissue: Placenta GSE144809,GSE144812 cell line
#> GSM5750634 infection group: Inf.. GSE192528 lung
#> GSM5022425 tissue: Retinal Orga.. GSE164884 Retinal Organoid
#> GSM2754776 cell type: Tcells,st.. GSE103147 T cells
#> GSM3976595 cell type: hTERT-HM^.. GSE134896 myometrium
#> GSM3976593 cell type: hTERT-HM^.. GSE134896 myometrium
#> GSM3976594 cell type: hTERT-HM^.. GSE134896 myometrium
#> GSM3291117 tissue: liver,cell t.. GSE117329 primary hepatic stel..
#> GSM4451228 cell line: 22Pc-EP,p.. GSE147976 Prostate cancer cell..
#> GSM4648524 cell line: MOLM13 ce.. GSE153637 MOLM13 cell line
#> CNT_geo_accession corr_TREATM_CNT
#> <character> <numeric>
#> GSM4296923 GSM4296920 0.984434
#> GSM5750634 GSM5750635 0.984432
#> GSM5022425 GSM5022423 0.988207
#> GSM2754776 GSM2754789 0.993094
#> GSM3976595 GSM3976579 0.987181
#> GSM3976593 GSM3976579 0.987186
#> GSM3976594 GSM3976579 0.985832
#> GSM3291117 GSM3291096 0.966951
#> GSM4451228 GSM4451226 0.985999
#> GSM4648524 GSM4648516 0.989333
#>
#> $caMKL1
#> DataFrame with 10 rows and 8 columns
#> Zscore TREATM_geo_accession title
#> <numeric> <character> <character>
#> GSM4296923 3.66232 GSM4296923 TEAD4 knock Down HTS..
#> GSM5750634 3.33309 GSM5750634 HLTE_A20_IAV infection
#> GSM5022425 3.20173 GSM5022425 D7_hiPSC1_NAM_2
#> GSM2754776 3.12605 GSM2754776 143_09_0278_DAY_0_T_..
#> GSM3976595 3.07313 GSM3976595 FSK_IL-1beta_3
#> GSM3976593 3.02627 GSM3976593 FSK_IL-1beta_1
#> GSM3976594 3.02548 GSM3976594 FSK_IL-1beta_2
#> GSM3291117 2.97296 GSM3291117 TGFb and JQ1 treated..
#> GSM4451228 2.90588 GSM4451228 22PcEP_FBS_NRG1_1
#> GSM4648524 2.58702 GSM4648524 RNAseq_RBMX_KD2_1
#> characteristics_ch1 series_id source_name_ch1
#> <character> <character> <character>
#> GSM4296923 tissue: Placenta GSE144809,GSE144812 cell line
#> GSM5750634 infection group: Inf.. GSE192528 lung
#> GSM5022425 tissue: Retinal Orga.. GSE164884 Retinal Organoid
#> GSM2754776 cell type: Tcells,st.. GSE103147 T cells
#> GSM3976595 cell type: hTERT-HM^.. GSE134896 myometrium
#> GSM3976593 cell type: hTERT-HM^.. GSE134896 myometrium
#> GSM3976594 cell type: hTERT-HM^.. GSE134896 myometrium
#> GSM3291117 tissue: liver,cell t.. GSE117329 primary hepatic stel..
#> GSM4451228 cell line: 22Pc-EP,p.. GSE147976 Prostate cancer cell..
#> GSM4648524 cell line: MOLM13 ce.. GSE153637 MOLM13 cell line
#> CNT_geo_accession corr_TREATM_CNT
#> <character> <numeric>
#> GSM4296923 GSM4296920 0.984434
#> GSM5750634 GSM5750635 0.984432
#> GSM5022425 GSM5022423 0.988207
#> GSM2754776 GSM2754789 0.993094
#> GSM3976595 GSM3976579 0.987181
#> GSM3976593 GSM3976579 0.987186
#> GSM3976594 GSM3976579 0.985832
#> GSM3291117 GSM3291096 0.966951
#> GSM4451228 GSM4451226 0.985999
#> GSM4648524 GSM4648516 0.989333
Example queries using the decomposed mouse MKL1 data:
# query the database using only the "caMKL1" mouse contrast, suppress plotting:
# for demonstration purposes (for actual analyses, use 'mode = "ANALYSIS"'):
query.res.mouse <- queryWithContrasts(dec_MKL1_mouse[, "caMKL1"], organism = "Mouse",
BPPARAM = params, verbose = FALSE,
plotType = "none", mode = "DEMO")
# plot results for individual contrasts using violin plots:
ViolinPlots_mouse <- plotQueryResultsViolin(query.res.mouse, doPlot = FALSE)
ViolinPlots_mouse[["caMKL1"]]
# plot results for individual contrasts using composite Manhattan/Density plots:
ManhDensPlots_mouse <- plotQueryResultsManh(query.res.mouse, doPlot = FALSE)
ManhDensPlots_mouse[["caMKL1"]]
#Information on the top hits of the query using the residual mouse caMKL1 contrasts:
query.res.mouse$TopHits$RESIDUAL_CONTRASTS
#> $caMKL1
#> DataFrame with 10 rows and 8 columns
#> Zscore TREATM_geo_accession title
#> <numeric> <character> <character>
#> GSM5170833 4.03706 GSM5170833 P2-1
#> GSM1958354 3.59669 GSM1958354 7dayC57CS3
#> GSM1902647 3.40848 GSM1902647 16_RightAtriumAHR_WT..
#> GSM5823264 2.98148 GSM5823264 ONCM, biological rep 1
#> GSM5825381 2.85630 GSM5825381 MYCS62A_Splenomegaly..
#> GSM2535506 2.84680 GSM2535506 48h +T Rep 2
#> GSM3939086 2.51029 GSM3939086 Mm PNC p53NULL_2
#> GSM5073873 2.49571 GSM5073873 WT-DTX.3
#> GSM3188993 2.44023 GSM3188993 Female_Het_D08_Rep1
#> GSM3942315 2.38631 GSM3942315 945351-2
#> characteristics_ch1 series_id source_name_ch1
#> <character> <character> <character>
#> GSM5170833 cell type: an insuli.. GSE168825 β-TC6 cells
#> GSM1958354 strain: C57BL6,genot.. GSE75513 Alveolar macrophages
#> GSM1902647 cell type: cardiomyo.. GSE73787 mouse heart tissue
#> GSM5823264 cell type: cultured .. GSE193893 mouse DRG tissue
#> GSM5825381 cell type: Lymphocyt.. GSE193995 Thymic tissue
#> GSM2535506 treatment: Tamoxifen.. GSE96555,GSE96595 tatsf1floxCreERT2 em..
#> GSM3939086 tissue: Mouse TG neu.. GSE134220 Primary neuron culture
#> GSM5073873 strain: BALB/c,cell .. GSE166536 4T1 murine breast ca..
#> GSM3188993 Sex: Female,tissue: .. GSE115766,GSE115768 Liver
#> GSM3942315 tissue: Kidney glome.. GSE134327 Kidney glomerulus
#> CNT_geo_accession corr_TREATM_CNT
#> <character> <numeric>
#> GSM5170833 GSM5170831 0.989282
#> GSM1958354 GSM1958350 0.995726
#> GSM1902647 GSM1902627 0.989253
#> GSM5823264 GSM5823263 0.989100
#> GSM5825381 GSM5825372 0.853545
#> GSM2535506 GSM2535498 0.995861
#> GSM3939086 GSM3939082 0.909963
#> GSM5073873 GSM5073868 0.991759
#> GSM3188993 GSM3188995 0.989233
#> GSM3942315 GSM3942321 0.991093
The top hits obtained for the residual (specific) fractions of either MKL1 or caMKL1 contrasts both in human and mouse are more clearly separated from the background compared to those obtained from the input or decoded (non-specific) fractions.
More importantly closer inspection of those top hits for the residual contrasts in both experiments (e.g hits from series GSE77120, GSE112277 or GSE140898 in human or GSE164860 in mouse) reveal that they correspond to treatments involving either MKL/MRTFA overexpression or overexpression of the MKL related transcription factor MYOCD. Of note, these treatments were performed in various cell contexts, different from the ones of the MKL study under consideration (LM2 and 4T1 cell lines for mouse and human respectively).
In general, as in this example, the residual specific fraction of a DGE profile will be a better query proxy for molecularly and mechanistically related treatments as it is largely stripped of nuisance variance present in the original contrast.
On the other hand the decoded non-specific fraction and its corresponding query hits can also be of interest in some applications as they provide information on the extent and type of downstream/secondary or lateral treatment effects.
The orthos
package provides functionality for direct access to contrast databases of orthosData
with the
loadContrastDatabase()
function.
This can be used to retrieve contrast values for all or subsets of genes or metadata for specific datasets, e.g for hits identified with queryWithContrasts()
.
The organism-specific databases are compiled as HDF5SummarizedExperiment objects.
As was the case for the orthos
models, a database will be automatically downloaded from
ExperimentHub
and cached in the user ExperimentHub directory
(see ExperimentHub::getExperimentHubOption("CACHE")
) using the orthosData
companion package, the first time loadContrastDatabase()
is called for that
database either directly or via queryWithContrasts()
.
The HDF5SummarizedExperiment object contains pre-calculated INPUT, RESIDUAL and DECODED log2 fold change contrasts as well as
the corresponding expression CONTEXT in log2 CPM representation for all the datasets in orthosData
.
Extensive gene and contrast annotation is available in the object’s rowData
and colData
respectively.
organism <- "Mouse"
orthosDB <- loadContrastDatabase(organism = "Mouse", mode = "DEMO")
orthosDB
#> class: SummarizedExperiment
#> dim: 20339 988
#> metadata(0):
#> assays(4): INPUT_CONTRASTS DECODED_CONTRASTS RESIDUAL_CONTRASTS CONTEXT
#> rownames(20339): 0610005C13Rik 0610009B22Rik ... n-R5s146 n-R5s149
#> rowData names(15): seqnames start ... ENTREZ_GENE_ID ARCHS4_ID
#> colnames(988): GSM1053282 GSM1061179 ... GSM970471 GSM984548
#> colData names(33): aligned_reads channel_count ... HasAssignedCNT
#> CNTname
#Available contrast annotations:
colnames(colData(orthosDB))
#> [1] "aligned_reads" "channel_count" "characteristics_ch1"
#> [4] "contact_address" "contact_city" "contact_country"
#> [7] "contact_institute" "contact_name" "contact_zip"
#> [10] "data_processing" "extract_protocol_ch1" "geo_accession"
#> [13] "instrument_model" "last_update_date" "library_selection"
#> [16] "library_source" "library_strategy" "molecule_ch1"
#> [19] "organism_ch1" "platform_id" "relation"
#> [22] "series_id" "singlecellprobability" "source_name_ch1"
#> [25] "sra_id" "status" "submission_date"
#> [28] "taxid_ch1" "title" "type"
#> [31] "Cor2CNT" "HasAssignedCNT" "CNTname"
#Available gene annotations:
colnames(rowData(orthosDB))
#> [1] "seqnames" "start"
#> [3] "end" "width"
#> [5] "strand" "ENSEMBL_GENE_ID"
#> [7] "gene_name" "gene_biotype"
#> [9] "seq_coord_system" "description"
#> [11] "ENSEMBL_GENE_ID_VERSION" "ENSEMBL_CANONICAL_TRANSCRIPT"
#> [13] "GENE_SYMBOL" "ENTREZ_GENE_ID"
#> [15] "ARCHS4_ID"
#Retrieve partial annotation for a specific contrast
#returned as a top-hit in the mouse caMKL1 query above:
colData(orthosDB)["GSM5021181", c("title", "series_id", "CNTname")]
#> DataFrame with 1 row and 3 columns
#> title series_id CNTname
#> <character> <character> <character>
#> GSM5021181 E0771_pRetrox_MRTFA_.. GSE164860 GSM5021179
# Compare context and individual contrast fractions between
# the mouse caMKL1 contrast under consideration and the "GSM5021181"
# query hit:
par(mfrow = c(2, 2))
queryID <- "GSM5021181"
for (contrast in names(assays(dec_MKL1_mouse))[c(4, 1, 2, 3)]) {
unit <- "L2FC"
if (contrast == "CONTEXT") {unit <- "L2CPM"}
plot(assays(dec_MKL1_mouse)[[contrast]][, "caMKL1"],
assays(orthosDB)[[contrast]][, queryID],
pch = 16, cex = 0.5, col = "darkslategrey", main = contrast,
xlab = paste0(unit, " caMKL1"), ylab = paste0(unit, " ", queryID))
abline(0, 1, col = "darkred", lwd = 0.8, lty = 2)
}
As (1) typical orthos
use cases do not require direct access to the models and
(2) use of the models requires loading of a conda environment via basilisk
this functionality
is by default not exposed to the user and is carried out transparently
by the non-exported functions .predictEncoder()
and .predictEncoderD()
.
However, as we envision cases where directly accessing the models might be of interest we provide here a brief overview and examples for direct calls to these functions.
The orthos
models are implemented in Keras
. For each organism there are two types of models:
orthos
models).As noted previously the first time these models are requested either by decomposeVar
or directly
by .predictEncoder()
and .predictEncoderD()
they are downloaded and
cached in the user ExperimentHub directory (see ExperimentHub::getExperimentHubOption("CACHE")
)
using the orthosData
companion package.
When calling the .predictEncoder()
and .predictEncoderD()
methods directly be attentive to the following:
orthos:::orthosenv
needs to be activated using basilisk::basiliskStart()
basilisk::basiliskStop()
.We now demonstrate calls to the context encoder for generating a latent embedding of a specific context and to the contrast conditional variational autoencoder for producing a contrast latent embedding and decoding.
# mouse MKL1 context and contrast with the appropriate shape and representation.
#
# Shape of models input is M x N,
# where M is the number of conditions,
# N the number of features -i.e orthos Genes
#
# Representation is L2CPMs for contexts and L2FCs for contrasts.
#
CONTEXT <- t(assay(dec_MKL1_mouse,"CONTEXT")[,1])
CONTRAST <- t(assay(dec_MKL1_mouse,"INPUT_CONTRASTS")[,1])
# Activate the `basilisk` environment:
library(basilisk)
cl <- basiliskStart(orthos:::orthosenv,
testload = "tensorflow")
# Produce a latent embedding for the context with .predictEncoder:
LATC <- basilisk::basiliskRun(proc = cl,
fun = orthos:::.predictEncoder,
organism = "Mouse",
gene_input = CONTEXT)
#> 1/1 - 0s - 94ms/epoch - 94ms/step
# Produce a latent embedding and decoding for the contrast with .predictEncoderD:
res <- basilisk::basiliskRun(proc = cl,
fun = orthos:::.predictEncoderD,
organism = "Mouse",
delta_input = CONTRAST, context = LATC)
#> 1/1 - 0s - 125ms/epoch - 125ms/step
#> 1/1 - 0s - 94ms/epoch - 94ms/step
# Deactivate the `basilisk` environment:
basilisk::basiliskStop(cl)
# Access the contrast latent embedding and decoding from the .predictEncoderD returned result:
LATD <- res$LATD
DEC <- res$DEC
Calls similar to the ones above are carried out under the hood when
decomposeVar()
is called.
sessionInfo()
#> R version 4.4.0 beta (2024-04-15 r86425 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows Server 2022 x64 (build 20348)
#>
#> Matrix products: default
#>
#>
#> locale:
#> [1] LC_COLLATE=C
#> [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] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] basilisk_1.16.0 reticulate_1.36.1
#> [3] ggplot2_3.5.1 keras_2.15.0
#> [5] orthosData_1.1.0 orthos_1.2.0
#> [7] SummarizedExperiment_1.34.0 Biobase_2.64.0
#> [9] GenomicRanges_1.56.0 GenomeInfoDb_1.40.0
#> [11] IRanges_2.38.0 S4Vectors_0.42.0
#> [13] BiocGenerics_0.50.0 MatrixGenerics_1.16.0
#> [15] matrixStats_1.3.0 knitr_1.46
#> [17] BiocStyle_2.32.0
#>
#> loaded via a namespace (and not attached):
#> [1] DBI_1.2.2 rlang_1.1.3 magrittr_2.0.3
#> [4] compiler_4.4.0 RSQLite_2.3.6 dir.expiry_1.12.0
#> [7] png_0.1-8 vctrs_0.6.5 stringr_1.5.1
#> [10] pkgconfig_2.0.3 crayon_1.5.2 fastmap_1.1.1
#> [13] backports_1.4.1 dbplyr_2.5.0 XVector_0.44.0
#> [16] labeling_0.4.3 fontawesome_0.5.2 utf8_1.2.4
#> [19] rmarkdown_2.26 UCSC.utils_1.0.0 purrr_1.0.2
#> [22] bit_4.0.5 xfun_0.43 zlibbioc_1.50.0
#> [25] cachem_1.0.8 jsonlite_1.8.8 blob_1.2.4
#> [28] highr_0.10 rhdf5filters_1.16.0 DelayedArray_0.30.0
#> [31] Rhdf5lib_1.26.0 BiocParallel_1.38.0 tensorflow_2.16.0
#> [34] broom_1.0.5 parallel_4.4.0 R6_2.5.1
#> [37] stringi_1.8.3 bslib_0.7.0 car_3.1-2
#> [40] jquerylib_0.1.4 Rcpp_1.0.12 bookdown_0.39
#> [43] base64enc_0.1-3 Matrix_1.7-0 tidyselect_1.2.1
#> [46] abind_1.4-5 yaml_2.3.8 codetools_0.2-20
#> [49] curl_5.2.1 plyr_1.8.9 lattice_0.22-6
#> [52] tibble_3.2.1 withr_3.0.0 basilisk.utils_1.16.0
#> [55] KEGGREST_1.44.0 evaluate_0.23 BiocFileCache_2.12.0
#> [58] ExperimentHub_2.12.0 Biostrings_2.72.0 pillar_1.9.0
#> [61] BiocManager_1.30.22 ggpubr_0.6.0 filelock_1.0.3
#> [64] carData_3.0-5 whisker_0.4.1 generics_0.1.3
#> [67] BiocVersion_3.19.1 munsell_0.5.1 scales_1.3.0
#> [70] glue_1.7.0 tools_4.4.0 AnnotationHub_3.12.0
#> [73] ggsignif_0.6.4 cowplot_1.1.3 rhdf5_2.48.0
#> [76] grid_4.4.0 tidyr_1.3.1 AnnotationDbi_1.66.0
#> [79] colorspace_2.1-0 GenomeInfoDbData_1.2.12 HDF5Array_1.32.0
#> [82] cli_3.6.2 rappdirs_0.3.3 tfruns_1.5.3
#> [85] fansi_1.0.6 viridisLite_0.4.2 S4Arrays_1.4.0
#> [88] dplyr_1.1.4 gtable_0.3.5 ggsci_3.0.3
#> [91] rstatix_0.7.2 zeallot_0.1.0 sass_0.4.9
#> [94] digest_0.6.35 SparseArray_1.4.0 ggrepel_0.9.5
#> [97] farver_2.1.1 memoise_2.0.1 htmltools_0.5.8.1
#> [100] lifecycle_1.0.4 httr_1.4.7 mime_0.12
#> [103] bit64_4.0.5
Hu, Xiao, Zongzhi Z Liu, Xinyue Chen, Vincent P Schulz, Abhishek Kumar, Amaleah A Hartman, Jason Weinstein, et al. 2019. “MKL1-Actin Pathway Restricts Chromatin Accessibility and Prevents Mature Pluripotency Activation.” Nature Communications 10 (1): 1695.
Lachmann, Alexander, Denis Torre, Alexandra B Keenan, Kathleen M Jagodnik, Hoyjin J Lee, Lily Wang, Moshe C Silverstein, and Avi Ma’ayan. 2018. “Massive Mining of Publicly Available Rna-Seq Data from Human and Mouse.” Nature Communications 9 (1): 1366.