In lm687/CompSign: Study of mutational signatures using compositional data analysis (CoDA) techniques
knitr::opts_chunk$set(
collapse = TRUE,
cache = FALSE,
comment = "#>"
)
CompSign is a toolkit for differential abundance analysis of mutational signatures using a mixed effects Dirichlet-multinominal model (or simpler variations). The compositional nature of mutational signature exposures has often been overlooked but has important implications, as the analyses must be done in relative terms.
paste0('Date of compilation: ', Sys.time())
Installation
CompSign can be installed as usual from github:
library(devtools)
devtools::install_github("lm687/CompSign")
These are the libraries used in this vignette:
library(CompSign)
library(gridExtra)
library(TMB)
library(HMP)
library(gridExtra)
library(ggplot2)
theme_set(theme_bw())
Datasets
Several datasets and objects can be loaded as follows using the data() function:
data(PancEndocrine_signaturesMSE)
data(ProstAdenoCA_chrom)
data(res_patient_lambda)
data(simplified_object)
data(ProstAdenoCA_trinucleotide)
data(obj_multilambda)
data(obj_multilambda_parameters)
The package contains the following datasets of exposures of mutational signatures and metadata of the corresponding samples. These datasets are:
PancEndocrine_signaturesMSE: Signature exposures for early and late mutations, in the PCAWG Panc-Endocrine cohortProstAdenoCA_chrom: Signature exposures for each chromosome, in the PCAWG Prost-AdenoCA cohortsimplified_object: Nucleotide changes in clonal and subclonal mutations in the Lung-AdenoCA cohortProstAdenoCA_trinucleotide: Trinucleotide abundandances for clonal and subclonal mutations in the PCAWG Prost-AdenoCA cohort
Minimal example workflow (Example 1)
This is a minimal example for running the model and extracting the estimated coefficients.
Input dataset
The input dataset is the argument object, which is a list with the following structure:
- x: covariate matrix (p x n)
- z: matrix of random effects indicating which patient-specific subsample belongs to which patient (n x N)
- Y (n x d)
- d
Fitting the model
The function wrapper_run_TMB() is used to run all variations of the model.
A minimal example is:
diagDM_no_small_sigs <- wrapper_run_TMB(object = simplified_object,
model = "diagRE_DM", use_nlminb=T, smart_init_vals=F)
in which simplified_object is a list containing simplified_object$x, simplified_object$z, simplified_object$Y. The object input for wrapper_run_TMB is always a list containing a matrix of covariates x, a matrix of observations Y and a matrix matching observations patients and the labels corresponding to the random effects. All three matrices should have the same number of rows - the number of observations. In this case, simplified_object corresponds to the nucleotides abundances in clonal and subclonal mutations in the Lung-AdenoCA cohort. simplified_object$x, contains two columns (a column of 1 for the intercept and a column of 0 and 1 indicating whether the observation corresponds to clonal (0) or subclonal (1) mutations. simplified_object$z contains a matrix where each observation is paired with its corresponding patient, simplified_object$Y is a matrix of counts for the nucleotide abundances for each patient and (clonal/subclonal) group.
Here we are using the model diagRE_DM, which is a mixed-effect Dirichlet-multinomial model with non-correlated random effects and two dispersion parameters (one per group). See below for the list of alternative models available.
Data included in CompSign
You can load simplified_object, as well as several other example datasets, using the data funtion:
data(package='CompSign')
data(simplified_object)
Testing for differential abundance
Depending on the model, the output of wrapper_run_TMB contains the estimated values of different coefficients. The parameters of interest are beta_0, beta_1 and lambda. Both beta are values in log-ratio space in with respect to the last category included in the input object for wrapper_run_TMB (i.e. the last colun of object$Y). beta_0 indicates the coefficient for the baseline - in the two-group example of simplified_object, this baseline reflects the abundance in the first (clonal) group. beta_1 indicates the cofficient for the second column - here representing the difference between clonal and subclonal mutations. If beta_1 is a vector of 0, this indicates that none of the log-ratios of abundances have changed between the two groups, and hence that there is no differential abundance.
To test for differential abundance, a generalised Wald test can be used with the function wald_TMB_wrapper, which gives a p-value as output:
wald_TMB_wrapper(diagDM_no_small_sigs)
Visualising the output
The estimated beta coefficients can be visualised as follows:
library(ggplot2)
plot_betas(diagDM_no_small_sigs, return_ggplot = T)
These beta coefficients are in log-ratio space with respect to a common category (the last column of object$Y). Therefore, it is difficult to determine which of the abundances changed in absolute terms. Assuming most categories do not change (i.e. assuming minimal perturbation), we can detect which categories deviate from the general beta_1 value. Without being able to draw any conclusion about the baseline category, as we do not have standard error values for it, we can detect for which categories in the numerator their beta_1 has extreme values: this can be done as follows with the plot_betas function:
plot_betas(diagDM_no_small_sigs, return_ggplot = T, add_median = T, line_zero = F, add_confint = T)
Workflow with mutational signatures (example 2)
In this example workflow we determine differential abundance between clonal and subclonal mutations in Pancreatic neuroendocrine tumors.
sign objects: how to load them
PancEndocrine_signaturesMSE is an S4 object of class sign. sign objects include all information about the name of the samples, the mutation signature exposures (of all signatures, or only of the subset of active signatures as determined by Alexandrov et al. 2020). However, the use of these S4 objects is not necessary for running the models of differential abundance, which instead take a list as argument for the input data.
PancEndocrine_signaturesMSE = load_PCAWG("../inst/extdata/roo/Panc-Endocrine_signaturesMSE_ROO.RDS",
read_directly = T,
typedata = "signaturesMSE", override_warning_X_Z = T)
PancEndocrine_signaturesMSE_v2 = load_PCAWG(ct = "Panc-Endocrine", typedata = "signaturesMSE", path_to_data = "../inst/extdata/", load_all_sigs = F, override_warning_X_Z = T)
Signature extraction
If signature extraction has not been performed, the wrapper function extract_sigs_TMB_obj outputs signature exposures from a matrix of trinucleotide abundances. We use the trinucleotide abundances from the object ProstAdenoCA_trinucleotide. In this wrapper function two methods for signature extraction are implemented: mutSigExtractor and quadratic programming (QP); this can be specified in the signature_fitting_method argument. The version of signature definitions matrix can be specified in the argument signature_version (which can be v2 or v3 of the COSMIC signatures. Alternatively, any signature definition matrix can be specified in the signature_definition, where the signature definition should be a F x K signature definition matrix, where K is the number of signatures and F the number of features (for COSMIC signatures, F=96). The active signature list can be specified in subset_signatures.
active_sigs_in_ct <- c("SBS1", "SBS2", "SBS3", "SBS5", "SBS8", "SBS13", "SBS18", "SBS37", "SBS40", "SBS41")
exposures_from_mutSigExtractor <- extract_sigs_TMB_obj(dataset_obj_trinucleotide=ProstAdenoCA_trinucleotide,
subset_signatures = active_sigs_in_ct,
signature_definition=NULL, signature_fitting_method = 'mutSigExtractor')
exposures_from_QP <- extract_sigs_TMB_obj(dataset_obj_trinucleotide=ProstAdenoCA_trinucleotide,
subset_signatures = active_sigs_in_ct,
signature_definition=NULL, signature_fitting_method = 'QP')
plot(log(as.vector(exposures_from_mutSigExtractor$Y)), log(as.vector(exposures_from_QP$Y)),
xlab='Exposures from mutSigExtractor', ylab='Exposures from quadratic programming')
abline(coef = c(0,1), lty='dashed')
Visualisation of mutational signature exposures
All samples - clonal and subclonal - sorted by increasing SBS3 exposure:
createBarplot(normalise_rw(non_duplicated_rows(PancEndocrine_signaturesMSE$Y)),
order_labels = names(sort(non_duplicated_rows(PancEndocrine_signaturesMSE$Y)[,'SBS3'],
decreasing = F)), remove_labels=T,
custom_color_palette=color_list)+ggtitle('Sorted by SBS3')
Subsetting the exposures
We create a simplified object containing exposures of fewer signatures (i.e. a subcomposition of the original signature vectors):
simplified_object <- give_subset_sigs_TMBobj(PancEndocrine_signaturesMSE,
sigs_to_remove = c('SBS13', 'SBS17a', 'SBS17b', 'SBS30'))
In which simplified_object is a list containing simplified_object$x, simplified_object$z, simplified_object$Y, simplified_object$d (see Input dataset section below)
The clonal and subclonal exposures are, respectively, the two barplots below:
do.call('grid.arrange', list(grobs=lapply(split_matrix_in_half(simplified_object$Y), function(i) createBarplot(normalise_rw(i), remove_labels = T, custom_color_palette=color_list)), nrow=1))
Variations of the model
There are several variations on the model, each with a particular combination of the following:
- Base model: Dirichlet-multinomial in most cases, but two simpler models -- multinomial models diagRE_M and fullRE_M are implemented for comparison.
- Random effects: correlated or not (e.g. uncorrelated in diagRE_DM, correlated in fullRE_DM).
- Number of dispersion parameters: generally one per group (e.g. in diagRE_DM), but possibly fewer or more if specified in the name (e.g. one single dipersion parameter in diagRE_DM_singlelambda, one dispersion parameter per patient in diagRE_DM_patientlambda).
The names of the model correspond to {random effects}_{base model}_{particularities of the dispersion parameter, if any}.
The first column is the <model> argument in the function wrapper_run_TMB().
| Name of model (for user) | Description and usage |
|---|---|
| FE_DM_singlelambda | Dirichlet-multinomial with no RE and a single lambda which can be used in a two-group comparison (if we consider the dipersion to be the same in both groups), a multi-group comparison, or any regression setting. |
| FE_DM | Dirichlet-multinomial with no random effects and two lambda, used for the comparison between two groups in which we want to account for group-specific dispersion. Model slightly more complex than FE_DM_singlelambda |
| diagRE_M | Multinomial with non-correlated multivariate effects. Simpler model than diagRE_DM, as no dispersion is included |
| diagRE_DM_singlelambda | Dirichlet-multinomial with non-correlated multivariate RE and one: equivalent of diagRE_DM but with a shared . It can be used in a two-group comparison (if we consider the dipersion to be the same in both groups), a multi-group comparison, or any regression setting. |
| singleRE_DM | Dirichlet-multinomial with a single RE intercept and two lambda: simple model in which random effects are not multivariate, but where group-specific dispersion is needed, in a two-group comparison. |
| diagRE_DM | Dirichlet-multinomial with independent RE and two lambda: used most commonly throughout the paper. The data are matched to warrant multivariate RE, but correlations between categories are not explicitly modelled. Faster than fullRE_DM with often comparable results for beta_1. |
| fullRE_DM_singlelambda | Dirichlet-multinomial with correlated RE and two lambdab: equivalent of diagR_DM_singlelambda that can be used if categories have strong correlations. |
| fullRE_M | Multinomial with correlated RE: assuming no overdispersion, it can be used in any regression setting |
| fullRE_DM | Dirichlet-multinomial with correlated RE and two lambda: model with multivariate RE with correlations modelled explicitly. As a (K-1)*(K-1) covariance matrix is estimated, this model is not recommended where the ratio of signatures to samples is high. |
| diagRE_DM_patientlambda | Dirichlet-multinomial with non-correlated RE and one lambda$ per patient: model more complex than diag_RE_DM that can be used when there are several observations per patient and we wish to include a patient-specific dispersion parameter |
| fullRE_DM_patientlambda | Dirichlet-multinomial with correlated RE and one per patient: equivalent of diagRE_DM_patientlambda, but modelling correlations between categories |
Using an alternative model, and parameter interpretation (example 3)
We run the model diagRE_DM_singlelambda with the dataset simplified_object:
diagDM_no_small_sigs <- wrapper_run_TMB(object = simplified_object,
model = "diagRE_DM_singlelambda", use_nlminb=T, smart_init_vals=F)
This is the resulting object with the estimates and their standard deviations:
diagDM_no_small_sigs
We can see how in this case we only have a single estimate for the dispersion parameter (log_lambda). The number of $\beta$ (20) corresonds to $2*(d-1)$ where $d$ is the number of mutations in the object (11, as seen in the bar plots above). The last signature or category, SBS39, is used as baseline signature in the ALR transformation. All beta parameters are in ALR-transformed space, using SBS39 as baseline. $\beta_0$ and $\beta_1$ results are intercalated, i.e. the $\beta$ in the first row corresponds to the $\beta_0$ of SBS1 wrt (with respect to) SBS39, the $\beta$ in the third row corresponds to the $\beta_0$ of SBS2 wrt (with respect to) SBS39, the $\beta$ in the first row corresponds to the $\beta_0$ of SBS1 wrt (with respect to) SBS39, etc.
These $\beta$ estimates can be plotted as follows:
plot_betas(diagDM_no_small_sigs)
Using an alternative model with patient-specific dispersion (example 4)
In a scenario where we have several samples for each patient, which can still be grouped in two or more groups (fixed effects), we might be interested in patient-specific dispersion parameters, besides the patient-specific intercepts that we had. This has been implemented:
Simulating data
n <- 40 ## patients
num_samples_per_patient = 30
stopifnot((num_samples_per_patient %% 2)==0)
z <- do.call('rbind', replicate(expr = diag(n), n= num_samples_per_patient, simplify = F))
x <- cbind(rep(1, n*num_samples_per_patient),
c(rep(0, (n*num_samples_per_patient)/2),
rep(1, (n*num_samples_per_patient)/2)))
d = 5 ## number of signatures
random_effects <- matrix(runif(n*(d-1)), nrow=n, ncol=d-1)
fixed_effects <- matrix(runif(2*(d-1)), nrow=2, ncol=d-1)
alpha_mat <- x %*% fixed_effects + z %*% random_effects
alpha_mat_softmax = softmax( cbind(alpha_mat, 0) )
## lambda per patient
log_lambda_vec = runif(n, 3, 4)
## get patient-specific overdispersions
log_lambda_vec_per_obs <- z %*% log_lambda_vec
## simulate counts
Nm = rep(1000, n*num_samples_per_patient) ## fixed number of mutations
y = matrix(NA, nrow = n*num_samples_per_patient, ncol=d)
for(l in 1:(n*num_samples_per_patient)){
## for each patient-specific subsample
alpha_l = alpha_mat_softmax[l,]*exp(log_lambda_vec_per_obs[l])
y[l,] = HMP::Dirichlet.multinomial(Nrs = Nm[l], shape = alpha_l)
}
##' create object with information about covariates, which observations correspond
##' to which patients, and the input (counts)
obj_multilambda <- list(Y = y, x=x, z = z, num_individuals = n)
Estimate parameters
res_patient_lambda <- wrapper_run_TMB(obj_multilambda, model = 'diagRE_DM_patientlambda')
data(res_patient_lambda)
data(obj_multilambda)
data(obj_multilambda_parameters)
Comparison of estimates with ground truth
The true parameters used for the simulation are found in the object obj_multilambda_parameters.
We can see how the parameters of interest are well recovered:
- $\beta_1$
- $\log(\lambda)$
and two additional sets of parameters suffer from a bias due to having used the uncorrelated version of the model (i.e. diagRE_DM_patientlambda instead of fullDMpatientlambda):
- $\beta_0$
- Random intercepts for each patient
Note that in the current setting the number of parameters to estimate is already large, and this is why
diagRE_DM_patientlambda is used.
appender <- function(string){
sapply(string, function(stringb){
if(stringb == 'Intercept'){
TeX(paste("$\\beta_0$"))
}else if (stringb == 'Slope'){
TeX(paste("$\\beta_1$"))
}
})
}
fixed_effects_estimate <- matrix(python_like_select_name(res_patient_lambda$par.fixed, 'beta'), nrow=2)
random_effects_estimate <- matrix(res_patient_lambda$par.random, ncol=ncol(obj_multilambda_parameters$fixed_effects))
log_lambda_vec_estimate <- python_like_select_name(res_patient_lambda$par.fixed, 'log_lambda')
list_comparisons_true_and_estimate = list(
all_betas = data.frame(x = as.vector(obj_multilambda_parameters$fixed_effects),
y = as.vector(fixed_effects_estimate)),
beta_1 = data.frame(x = obj_multilambda_parameters$fixed_effects[2,],
y = fixed_effects_estimate[2,]),
random_intercepts = data.frame(x = as.vector(obj_multilambda_parameters$random_effects),
y = as.vector(random_effects_estimate)),
log_lambda = data.frame(x = as.vector(obj_multilambda_parameters$log_lambda_vec),
y = as.vector(log_lambda_vec_estimate)))
plots_comparisons_true_and_estimate <- lapply(1:length(list_comparisons_true_and_estimate), function(i){
ggplot(list_comparisons_true_and_estimate[[i]], aes(x=x, y=y))+
geom_point()+
geom_abline(slope = 1, intercept = 0, lty='dashed')+
ggtitle(gsub('_', ' ', names(list_comparisons_true_and_estimate)[i]))+
labs(x='True value', y='Estimate')
})
do.call('grid.arrange', c(plots_comparisons_true_and_estimate, nrow=1))
Session info
sessionInfo()
lm687/CompSign documentation built on Jan. 16, 2025, 11:25 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, cache = FALSE, comment = "#>" )
CompSign is a toolkit for differential abundance analysis of mutational signatures using a mixed effects Dirichlet-multinominal model (or simpler variations). The compositional nature of mutational signature exposures has often been overlooked but has important implications, as the analyses must be done in relative terms.
paste0('Date of compilation: ', Sys.time())
Installation
CompSign can be installed as usual from github:
library(devtools) devtools::install_github("lm687/CompSign")
These are the libraries used in this vignette:
library(CompSign) library(gridExtra) library(TMB) library(HMP) library(gridExtra) library(ggplot2) theme_set(theme_bw())
Datasets
Several datasets and objects can be loaded as follows using the data() function:
data(PancEndocrine_signaturesMSE) data(ProstAdenoCA_chrom) data(res_patient_lambda) data(simplified_object) data(ProstAdenoCA_trinucleotide) data(obj_multilambda) data(obj_multilambda_parameters)
The package contains the following datasets of exposures of mutational signatures and metadata of the corresponding samples. These datasets are:
PancEndocrine_signaturesMSE: Signature exposures for early and late mutations, in the PCAWG Panc-Endocrine cohortProstAdenoCA_chrom: Signature exposures for each chromosome, in the PCAWG Prost-AdenoCA cohortsimplified_object: Nucleotide changes in clonal and subclonal mutations in the Lung-AdenoCA cohortProstAdenoCA_trinucleotide: Trinucleotide abundandances for clonal and subclonal mutations in the PCAWG Prost-AdenoCA cohort
Minimal example workflow (Example 1)
This is a minimal example for running the model and extracting the estimated coefficients.
Input dataset
The input dataset is the argument object, which is a list with the following structure:
- x: covariate matrix (p x n)
- z: matrix of random effects indicating which patient-specific subsample belongs to which patient (n x N)
- Y (n x d)
- d
Fitting the model
The function wrapper_run_TMB() is used to run all variations of the model.
A minimal example is:
diagDM_no_small_sigs <- wrapper_run_TMB(object = simplified_object, model = "diagRE_DM", use_nlminb=T, smart_init_vals=F)
in which simplified_object is a list containing simplified_object$x, simplified_object$z, simplified_object$Y. The object input for wrapper_run_TMB is always a list containing a matrix of covariates x, a matrix of observations Y and a matrix matching observations patients and the labels corresponding to the random effects. All three matrices should have the same number of rows - the number of observations. In this case, simplified_object corresponds to the nucleotides abundances in clonal and subclonal mutations in the Lung-AdenoCA cohort. simplified_object$x, contains two columns (a column of 1 for the intercept and a column of 0 and 1 indicating whether the observation corresponds to clonal (0) or subclonal (1) mutations. simplified_object$z contains a matrix where each observation is paired with its corresponding patient, simplified_object$Y is a matrix of counts for the nucleotide abundances for each patient and (clonal/subclonal) group.
Here we are using the model diagRE_DM, which is a mixed-effect Dirichlet-multinomial model with non-correlated random effects and two dispersion parameters (one per group). See below for the list of alternative models available.
Data included in CompSign
You can load simplified_object, as well as several other example datasets, using the data funtion:
data(package='CompSign') data(simplified_object)
Testing for differential abundance
Depending on the model, the output of wrapper_run_TMB contains the estimated values of different coefficients. The parameters of interest are beta_0, beta_1 and lambda. Both beta are values in log-ratio space in with respect to the last category included in the input object for wrapper_run_TMB (i.e. the last colun of object$Y). beta_0 indicates the coefficient for the baseline - in the two-group example of simplified_object, this baseline reflects the abundance in the first (clonal) group. beta_1 indicates the cofficient for the second column - here representing the difference between clonal and subclonal mutations. If beta_1 is a vector of 0, this indicates that none of the log-ratios of abundances have changed between the two groups, and hence that there is no differential abundance.
To test for differential abundance, a generalised Wald test can be used with the function wald_TMB_wrapper, which gives a p-value as output:
wald_TMB_wrapper(diagDM_no_small_sigs)
Visualising the output
The estimated beta coefficients can be visualised as follows:
library(ggplot2) plot_betas(diagDM_no_small_sigs, return_ggplot = T)
These beta coefficients are in log-ratio space with respect to a common category (the last column of object$Y). Therefore, it is difficult to determine which of the abundances changed in absolute terms. Assuming most categories do not change (i.e. assuming minimal perturbation), we can detect which categories deviate from the general beta_1 value. Without being able to draw any conclusion about the baseline category, as we do not have standard error values for it, we can detect for which categories in the numerator their beta_1 has extreme values: this can be done as follows with the plot_betas function:
plot_betas(diagDM_no_small_sigs, return_ggplot = T, add_median = T, line_zero = F, add_confint = T)
Workflow with mutational signatures (example 2)
In this example workflow we determine differential abundance between clonal and subclonal mutations in Pancreatic neuroendocrine tumors.
sign objects: how to load them
PancEndocrine_signaturesMSE is an S4 object of class sign. sign objects include all information about the name of the samples, the mutation signature exposures (of all signatures, or only of the subset of active signatures as determined by Alexandrov et al. 2020). However, the use of these S4 objects is not necessary for running the models of differential abundance, which instead take a list as argument for the input data.
PancEndocrine_signaturesMSE = load_PCAWG("../inst/extdata/roo/Panc-Endocrine_signaturesMSE_ROO.RDS", read_directly = T, typedata = "signaturesMSE", override_warning_X_Z = T) PancEndocrine_signaturesMSE_v2 = load_PCAWG(ct = "Panc-Endocrine", typedata = "signaturesMSE", path_to_data = "../inst/extdata/", load_all_sigs = F, override_warning_X_Z = T)
Signature extraction
If signature extraction has not been performed, the wrapper function extract_sigs_TMB_obj outputs signature exposures from a matrix of trinucleotide abundances. We use the trinucleotide abundances from the object ProstAdenoCA_trinucleotide. In this wrapper function two methods for signature extraction are implemented: mutSigExtractor and quadratic programming (QP); this can be specified in the signature_fitting_method argument. The version of signature definitions matrix can be specified in the argument signature_version (which can be v2 or v3 of the COSMIC signatures. Alternatively, any signature definition matrix can be specified in the signature_definition, where the signature definition should be a F x K signature definition matrix, where K is the number of signatures and F the number of features (for COSMIC signatures, F=96). The active signature list can be specified in subset_signatures.
active_sigs_in_ct <- c("SBS1", "SBS2", "SBS3", "SBS5", "SBS8", "SBS13", "SBS18", "SBS37", "SBS40", "SBS41") exposures_from_mutSigExtractor <- extract_sigs_TMB_obj(dataset_obj_trinucleotide=ProstAdenoCA_trinucleotide, subset_signatures = active_sigs_in_ct, signature_definition=NULL, signature_fitting_method = 'mutSigExtractor') exposures_from_QP <- extract_sigs_TMB_obj(dataset_obj_trinucleotide=ProstAdenoCA_trinucleotide, subset_signatures = active_sigs_in_ct, signature_definition=NULL, signature_fitting_method = 'QP') plot(log(as.vector(exposures_from_mutSigExtractor$Y)), log(as.vector(exposures_from_QP$Y)), xlab='Exposures from mutSigExtractor', ylab='Exposures from quadratic programming') abline(coef = c(0,1), lty='dashed')
Visualisation of mutational signature exposures
All samples - clonal and subclonal - sorted by increasing SBS3 exposure:
createBarplot(normalise_rw(non_duplicated_rows(PancEndocrine_signaturesMSE$Y)), order_labels = names(sort(non_duplicated_rows(PancEndocrine_signaturesMSE$Y)[,'SBS3'], decreasing = F)), remove_labels=T, custom_color_palette=color_list)+ggtitle('Sorted by SBS3')
Subsetting the exposures
We create a simplified object containing exposures of fewer signatures (i.e. a subcomposition of the original signature vectors):
simplified_object <- give_subset_sigs_TMBobj(PancEndocrine_signaturesMSE, sigs_to_remove = c('SBS13', 'SBS17a', 'SBS17b', 'SBS30'))
In which simplified_object is a list containing simplified_object$x, simplified_object$z, simplified_object$Y, simplified_object$d (see Input dataset section below)
The clonal and subclonal exposures are, respectively, the two barplots below:
do.call('grid.arrange', list(grobs=lapply(split_matrix_in_half(simplified_object$Y), function(i) createBarplot(normalise_rw(i), remove_labels = T, custom_color_palette=color_list)), nrow=1))
Variations of the model
There are several variations on the model, each with a particular combination of the following:
- Base model: Dirichlet-multinomial in most cases, but two simpler models -- multinomial models diagRE_M and fullRE_M are implemented for comparison.
- Random effects: correlated or not (e.g. uncorrelated in diagRE_DM, correlated in fullRE_DM).
- Number of dispersion parameters: generally one per group (e.g. in diagRE_DM), but possibly fewer or more if specified in the name (e.g. one single dipersion parameter in diagRE_DM_singlelambda, one dispersion parameter per patient in diagRE_DM_patientlambda).
The names of the model correspond to {random effects}_{base model}_{particularities of the dispersion parameter, if any}.
The first column is the <model> argument in the function wrapper_run_TMB().
| Name of model (for user) | Description and usage |
|---|---|
| FE_DM_singlelambda | Dirichlet-multinomial with no RE and a single lambda which can be used in a two-group comparison (if we consider the dipersion to be the same in both groups), a multi-group comparison, or any regression setting. |
| FE_DM | Dirichlet-multinomial with no random effects and two lambda, used for the comparison between two groups in which we want to account for group-specific dispersion. Model slightly more complex than FE_DM_singlelambda |
| diagRE_M | Multinomial with non-correlated multivariate effects. Simpler model than diagRE_DM, as no dispersion is included |
| diagRE_DM_singlelambda | Dirichlet-multinomial with non-correlated multivariate RE and one: equivalent of diagRE_DM but with a shared . It can be used in a two-group comparison (if we consider the dipersion to be the same in both groups), a multi-group comparison, or any regression setting. |
| singleRE_DM | Dirichlet-multinomial with a single RE intercept and two lambda: simple model in which random effects are not multivariate, but where group-specific dispersion is needed, in a two-group comparison. |
| diagRE_DM | Dirichlet-multinomial with independent RE and two lambda: used most commonly throughout the paper. The data are matched to warrant multivariate RE, but correlations between categories are not explicitly modelled. Faster than fullRE_DM with often comparable results for beta_1. |
| fullRE_DM_singlelambda | Dirichlet-multinomial with correlated RE and two lambdab: equivalent of diagR_DM_singlelambda that can be used if categories have strong correlations. |
| fullRE_M | Multinomial with correlated RE: assuming no overdispersion, it can be used in any regression setting |
| fullRE_DM | Dirichlet-multinomial with correlated RE and two lambda: model with multivariate RE with correlations modelled explicitly. As a (K-1)*(K-1) covariance matrix is estimated, this model is not recommended where the ratio of signatures to samples is high. |
| diagRE_DM_patientlambda | Dirichlet-multinomial with non-correlated RE and one lambda$ per patient: model more complex than diag_RE_DM that can be used when there are several observations per patient and we wish to include a patient-specific dispersion parameter |
| fullRE_DM_patientlambda | Dirichlet-multinomial with correlated RE and one per patient: equivalent of diagRE_DM_patientlambda, but modelling correlations between categories |
Using an alternative model, and parameter interpretation (example 3)
We run the model diagRE_DM_singlelambda with the dataset simplified_object:
diagDM_no_small_sigs <- wrapper_run_TMB(object = simplified_object, model = "diagRE_DM_singlelambda", use_nlminb=T, smart_init_vals=F)
This is the resulting object with the estimates and their standard deviations:
diagDM_no_small_sigs
We can see how in this case we only have a single estimate for the dispersion parameter (log_lambda). The number of $\beta$ (20) corresonds to $2*(d-1)$ where $d$ is the number of mutations in the object (11, as seen in the bar plots above). The last signature or category, SBS39, is used as baseline signature in the ALR transformation. All beta parameters are in ALR-transformed space, using SBS39 as baseline. $\beta_0$ and $\beta_1$ results are intercalated, i.e. the $\beta$ in the first row corresponds to the $\beta_0$ of SBS1 wrt (with respect to) SBS39, the $\beta$ in the third row corresponds to the $\beta_0$ of SBS2 wrt (with respect to) SBS39, the $\beta$ in the first row corresponds to the $\beta_0$ of SBS1 wrt (with respect to) SBS39, etc.
These $\beta$ estimates can be plotted as follows:
plot_betas(diagDM_no_small_sigs)
Using an alternative model with patient-specific dispersion (example 4)
In a scenario where we have several samples for each patient, which can still be grouped in two or more groups (fixed effects), we might be interested in patient-specific dispersion parameters, besides the patient-specific intercepts that we had. This has been implemented:
Simulating data
n <- 40 ## patients num_samples_per_patient = 30 stopifnot((num_samples_per_patient %% 2)==0) z <- do.call('rbind', replicate(expr = diag(n), n= num_samples_per_patient, simplify = F)) x <- cbind(rep(1, n*num_samples_per_patient), c(rep(0, (n*num_samples_per_patient)/2), rep(1, (n*num_samples_per_patient)/2))) d = 5 ## number of signatures random_effects <- matrix(runif(n*(d-1)), nrow=n, ncol=d-1) fixed_effects <- matrix(runif(2*(d-1)), nrow=2, ncol=d-1) alpha_mat <- x %*% fixed_effects + z %*% random_effects alpha_mat_softmax = softmax( cbind(alpha_mat, 0) ) ## lambda per patient log_lambda_vec = runif(n, 3, 4) ## get patient-specific overdispersions log_lambda_vec_per_obs <- z %*% log_lambda_vec ## simulate counts Nm = rep(1000, n*num_samples_per_patient) ## fixed number of mutations y = matrix(NA, nrow = n*num_samples_per_patient, ncol=d) for(l in 1:(n*num_samples_per_patient)){ ## for each patient-specific subsample alpha_l = alpha_mat_softmax[l,]*exp(log_lambda_vec_per_obs[l]) y[l,] = HMP::Dirichlet.multinomial(Nrs = Nm[l], shape = alpha_l) }
##' create object with information about covariates, which observations correspond ##' to which patients, and the input (counts) obj_multilambda <- list(Y = y, x=x, z = z, num_individuals = n)
Estimate parameters
res_patient_lambda <- wrapper_run_TMB(obj_multilambda, model = 'diagRE_DM_patientlambda')
data(res_patient_lambda) data(obj_multilambda) data(obj_multilambda_parameters)
Comparison of estimates with ground truth
The true parameters used for the simulation are found in the object obj_multilambda_parameters.
We can see how the parameters of interest are well recovered:
- $\beta_1$
- $\log(\lambda)$
and two additional sets of parameters suffer from a bias due to having used the uncorrelated version of the model (i.e. diagRE_DM_patientlambda instead of fullDMpatientlambda):
- $\beta_0$
- Random intercepts for each patient
Note that in the current setting the number of parameters to estimate is already large, and this is why
diagRE_DM_patientlambdais used.
appender <- function(string){ sapply(string, function(stringb){ if(stringb == 'Intercept'){ TeX(paste("$\\beta_0$")) }else if (stringb == 'Slope'){ TeX(paste("$\\beta_1$")) } }) }
fixed_effects_estimate <- matrix(python_like_select_name(res_patient_lambda$par.fixed, 'beta'), nrow=2) random_effects_estimate <- matrix(res_patient_lambda$par.random, ncol=ncol(obj_multilambda_parameters$fixed_effects)) log_lambda_vec_estimate <- python_like_select_name(res_patient_lambda$par.fixed, 'log_lambda') list_comparisons_true_and_estimate = list( all_betas = data.frame(x = as.vector(obj_multilambda_parameters$fixed_effects), y = as.vector(fixed_effects_estimate)), beta_1 = data.frame(x = obj_multilambda_parameters$fixed_effects[2,], y = fixed_effects_estimate[2,]), random_intercepts = data.frame(x = as.vector(obj_multilambda_parameters$random_effects), y = as.vector(random_effects_estimate)), log_lambda = data.frame(x = as.vector(obj_multilambda_parameters$log_lambda_vec), y = as.vector(log_lambda_vec_estimate))) plots_comparisons_true_and_estimate <- lapply(1:length(list_comparisons_true_and_estimate), function(i){ ggplot(list_comparisons_true_and_estimate[[i]], aes(x=x, y=y))+ geom_point()+ geom_abline(slope = 1, intercept = 0, lty='dashed')+ ggtitle(gsub('_', ' ', names(list_comparisons_true_and_estimate)[i]))+ labs(x='True value', y='Estimate') }) do.call('grid.arrange', c(plots_comparisons_true_and_estimate, nrow=1))
Session info
sessionInfo()
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.