Introduction
RolDE is a composite method, consisting of three independent modules
with different approaches to detecting longitudinal differential
expression. The combination of these diverse modules allows RolDE to
robustly detect varying differences in longitudinal trends and
expression levels in diverse data types and experimental settings.
The RegROTS module merges the power of regression modelling with the
power of the established differential expression method Reproducibility
Optimized Test Statistic (ROTS) (Elo et al., Suomi et al.). A polynomial
regression model for protein expression over time is fitted separately
for each replicate (individual) in each condition. Differential
expression between two replicates (individuals) in different conditions
is examined by comparing the coefficients of the replicate-specific
regression models. If all coefficient differences are zero, no
longitudinal differential expression between the two replicates in
different conditions exist. For a thorough exploration of differential
expression between the conditions, all possible combinations of
replicates in different conditions are examined.
In the DiffROTS module the expression of replicates (individuals) in
different conditions are directly compared at all time points. Again, if
the expression level differences at all time points are zero, no
differential expression between the examined replicates (individuals) in
different conditions exist. Similarly to the RegROTS module,
differential expression is examined between all possible combinations of
replicates (individuals) in the different conditions. In non-aligned
time point data, the overall expression level differences between the
conditions is examined when accounting for time-associated trends of
varying complexity in the data. More specifically, the expression level
differences between the conditions are examined when adjusting for
increasingly complex time-related expression trends of polynomial
degrees d=0,1,.,d where d is the maximum degree for the polynomial and
the same degree as is used for the PolyReg module.
In the PolyReg module, polynomial regression modelling is used to
detect longitudinal differential expression. Condition is included as a
categorical factor within the models and by investigating the condition
related intercept and longitudinal-term related coefficients at
different levels of the condition factor, average differences in
expression levels as well as differences in longitudinal expression
patterns between the conditions can be examined.
Finally, to conclusively detect any differential expression, the
findings from the different modules are combined using the rank product.
For more details about the method, see the original RolDE publication
(Valikangas et al.).
Applying RolDE in example datasets with aligned time points
First, load the dataset and the design matrix to be used:
data(data1)
data("des_matrix1")
To understand the structure of the design matrix and what a design
matrix for RolDE should look like, let’s explore the data and the design
matrix:
Dimensions of the data:
dim(data1)
#> [1] 1045 30
Column names of the data:
head(colnames(data1))
#> [1] "Condition1_timepoint_1_Sample_1_r1" "Condition1_timepoint_1_Sample_1_r2"
#> [3] "Condition1_timepoint_1_Sample_1_r3" "Condition1_timepoint_2_Sample_2_r1"
#> [5] "Condition1_timepoint_2_Sample_2_r2" "Condition1_timepoint_2_Sample_2_r3"
Dimensions of the design matrix:
dim(des_matrix1)
#> [1] 30 4
Column names of the design matrix:
colnames(des_matrix1)
#> [1] "Sample Names" "Condition" "Time point"
#> [4] "Replicate/Individual"
Contents of the design matrix:
| Condition1_timepoint_1_Sample_1_r1 |
Condition1 |
1 |
1 |
| Condition1_timepoint_1_Sample_1_r2 |
Condition1 |
1 |
2 |
| Condition1_timepoint_1_Sample_1_r3 |
Condition1 |
1 |
3 |
| Condition1_timepoint_2_Sample_2_r1 |
Condition1 |
2 |
1 |
| Condition1_timepoint_2_Sample_2_r2 |
Condition1 |
2 |
2 |
| Condition1_timepoint_2_Sample_2_r3 |
Condition1 |
2 |
3 |
| [25,] |
Condition2_timepoint_4_Sample_4_r1 |
Condition2 |
4 |
4 |
| [26,] |
Condition2_timepoint_4_Sample_4_r2 |
Condition2 |
4 |
5 |
| [27,] |
Condition2_timepoint_4_Sample_4_r3 |
Condition2 |
4 |
6 |
| [28,] |
Condition2_timepoint_5_Sample_5_r1 |
Condition2 |
5 |
4 |
| [29,] |
Condition2_timepoint_5_Sample_5_r2 |
Condition2 |
5 |
5 |
| [30,] |
Condition2_timepoint_5_Sample_5_r3 |
Condition2 |
5 |
6 |
The first column of the design matrix should contain the sample names
of the data (column names of the data). The second column should
indicate the condition/group factor status for each sample to be
explored (e.g. sick/healthy, control/case). The third column should
indicate the time point for each sample (data with aligned time points)
or time value (or equivalent, e.g. age, temperature, disease
progression) information in data with non-aligned time points. The
fourth and final column of the design matrix should indicate the
replicate / individual each sample came from. Let’s look at the
exemplary design matrix a little bit more:
Unique samples:
unique(des_matrix1[,1])
#> [1] "Condition1_timepoint_1_Sample_1_r1" "Condition1_timepoint_1_Sample_1_r2"
#> [3] "Condition1_timepoint_1_Sample_1_r3" "Condition1_timepoint_2_Sample_2_r1"
#> [5] "Condition1_timepoint_2_Sample_2_r2" "Condition1_timepoint_2_Sample_2_r3"
#> [7] "Condition1_timepoint_3_Sample_3_r1" "Condition1_timepoint_3_Sample_3_r2"
#> [9] "Condition1_timepoint_3_Sample_3_r3" "Condition1_timepoint_4_Sample_4_r1"
#> [11] "Condition1_timepoint_4_Sample_4_r2" "Condition1_timepoint_4_Sample_4_r3"
#> [13] "Condition1_timepoint_5_Sample_5_r1" "Condition1_timepoint_5_Sample_5_r2"
#> [15] "Condition1_timepoint_5_Sample_5_r3" "Condition2_timepoint_1_Sample_1_r1"
#> [17] "Condition2_timepoint_1_Sample_1_r2" "Condition2_timepoint_1_Sample_1_r3"
#> [19] "Condition2_timepoint_2_Sample_2_r1" "Condition2_timepoint_2_Sample_2_r2"
#> [21] "Condition2_timepoint_2_Sample_2_r3" "Condition2_timepoint_3_Sample_3_r1"
#> [23] "Condition2_timepoint_3_Sample_3_r2" "Condition2_timepoint_3_Sample_3_r3"
#> [25] "Condition2_timepoint_4_Sample_4_r1" "Condition2_timepoint_4_Sample_4_r2"
#> [27] "Condition2_timepoint_4_Sample_4_r3" "Condition2_timepoint_5_Sample_5_r1"
#> [29] "Condition2_timepoint_5_Sample_5_r2" "Condition2_timepoint_5_Sample_5_r3"
Conditions:
unique(des_matrix1[,2])
#> [1] "Condition1" "Condition2"
table(des_matrix1[,2])
Time points:
unique(des_matrix1[,3])
#> [1] "1" "2" "3" "4" "5"
table(des_matrix1[,3])
Replicates:
unique(des_matrix1[,4])
#> [1] "1" "2" "3" "4" "5" "6"
table(des_matrix1[,4])
In this example, we have a dataset of 30 samples, 2 conditions, 6
replicates each with 5 time points. The timepoints in the data are
aligned.
This is how a design matrix for a dataset for RolDE should look like.
This gives all the essential information for RolDE it needs to determine
longitudinal differential expression between the conditions.
By bare minimum, RolDE needs the data and the design matrix. By
default RolDE assumes that the time points in the data are aligned. If
not defined otherwise, RolDE will by default use sequential processing.
However, using parallel processing and multiple threads will greatly
reduce the computational time required by RolDE.
Please use set.seed() for reproducibility.
Running RolDE using parallel processing and 3 threads:
set.seed(1)
data1.res<-RolDE(data=data1, des_matrix=des_matrix1, n_cores=3)
RolDE supports the SummarizedExperiment data structure and the data
and design matrix can be provided to RolDE as a SummarizedExperiment
object. In this case, the data argument for RolDE must be a
SummarizedExperiment object, where the data matrix is included as a list
in the assays argument and the design matrix in the
colData argument. The format of the data matrix and the design
matrix within the SummarizedExperiment object must be the same as when
provided separately.
Constructing a SummarizedExperiment object from data1 and the
associated design matrix for RolDE:
SE_object_for_RolDE = SummarizedExperiment(assays=list(data1),
colData = des_matrix1)
Running RolDE using a SummarizedExperiment object including the data
and the metadata:
set.seed(1)
data1.res<-RolDE(data=SE_object_for_RolDE, n_cores=3)
The results of RolDE are returned in a list with lot of
information:
names(data1.res)
#> [1] "RolDE_Results" "RegROTS_Results" "RegROTS_P_Values"
#> [4] "DiffROTS_Results" "DiffROTS_P_Values" "PolyReg_Results"
#> [7] "PolyReg_P_Values" "ROTS_Runs" "Method_Degrees"
#> [10] "Input"
The main results of RolDE are located in the first element of the
provided result list. Elements 2,4 and 6 provide the result for the
different modules of RolDE separately (the RegROTS, DiffROTS and PolyReg
modules).
Elements 3 and 5 provide the significance values for the different
ROTS runs of the RegROTS and DiffROTS modules, respectively. The used
ROTS runs within those modules are given in element 8. The comparisons
between the replicates in the different conditions are divided into
different runs so that each sample is used only once within each run to
preserve the proper degrees of freedom for statistical testing.
All the condition related significance values of the polynomial
regression models in the PolyReg module are given in element 7. The used
polynomial degrees for the regression models in the RegROTS and the
PolyReg modules are given in element 9. And in element 10, all the
inputs used by RolDE (both given by the user and those determined
automatically by the method) are given.
Typically, the main (only) interest of the user is in element 1,
where the main results of RolDE are given:
RolDE.data1<-data1.res$RolDE_Results
dim(RolDE.data1)
#> [1] 1045 4
colnames(RolDE.data1)
#> [1] "Feature ID"
#> [2] "RolDE Rank Product"
#> [3] "Estimated Significance Value"
#> [4] "Adjusted Estimated Significance Value"
By default, the features in the result data frame are given in the
same order as entered.
Let’s order the results based on the strength of longitudinal
differential expression detected by RolDE:
RolDE.data1<-RolDE.data1[order(as.numeric(RolDE.data1[,2])),]
head(RolDE.data1, 5)
| sp|P32263|P5CR_YEAST |
sp|P32263|P5CR_YEAST |
4.447960 |
0.0010988 |
0.9983582 |
| P01344ups|IGF2_HUMAN_UPS |
P01344ups|IGF2_HUMAN_UPS |
6.073178 |
0.0028508 |
0.9983582 |
| sp|Q04935|COX20_YEAST |
sp|Q04935|COX20_YEAST |
7.989570 |
0.0044291 |
0.9983582 |
| sp|P40008|FMP52_YEAST |
sp|P40008|FMP52_YEAST |
9.356095 |
0.0053161 |
0.9983582 |
| sp|P15180|SYKC_YEAST |
sp|P15180|SYKC_YEAST |
9.619002 |
0.0061172 |
0.9983582 |
Explore the distribution of the estimated significance values:
hist(RolDE.data1$`Estimated Significance Value`, main = "Estimated significance values", xlab="")
As can be observed, the distribution of the estimated significance
values is approximately uniform. Data1 is random data; the values for
the features have been randomly assigned from a normal distribution with
a mean of 22 and a standard deviation of 1.5. Overall, the null
hypothesis between Condition1 and Condition2 is true; the conditions are
not differentially expressed and a uniform significance value
distribution is expected.
After correcting for the simultaneous testing of multiple hypothesis
by using the Bejamini-Hochberg (FDR) correction, no differentially
epxressed feature remains with the commonly used FDR level of 0.05:
length(which(RolDE.data1$`Adjusted Estimated Significance Value`<=0.05))
#> [1] 0
The calculation of significance values in RolDE is controlled via the
parameter sigValSampN which control if the significance values
should be calculated and how many permutations should be used when
determining the significance values. Computational time required by
RolDE can be reduced by not calculating the significance values by
setting sigValSampN to 0. By increasing sigValSampN
from the default number, the significance values can be estimated more
accurately but the required computational time will also increase. If
parallel processing for RolDE is enabled (n_cores > 1), it
will also be utilized when estimating the significance values.
The estimated significance values can be adjusted by any method
included in the p.adjust method in the stats package. Alternatively,
q-values as defined by Storey et al. in the Bioconductor package qvalue
can be used. Valid values for sig_adj_meth are then: “holm”,
“hochberg”, “hommel”, “bonferroni”, “BH”, “BY”,“fdr”,“none”,
“qvalue”.
Semi-simulated spike-in dataset
In addition to the random data already explored, a semi-simulated
spike-in dataset with differential expression between the conditions is
also included in the RolDE package:
data("data3")
data("des_matrix3")
?data3
A semi-simulated UPS1 spike-in dataset with differences in longitudinal expression for the spike-in proteins - no missing values.
Description
A longitudinal proteomics dataset with five timepoints, two conditions and
three replicates for each sample at each timepoint in each condition. The expression
values for each protein in each sample have been generated by using the rnorm
function and the means and standard deviations of the corresponding proteins and samples
in the original experimental UPS1 spike-in data. In this manner, the expression for each
protein in each sample could be replicated with some random variation multiple times when
necessary. The pattern of missing values was directly copied from the used samples
in the original UPS1 spike-in dataset. All proteins with missing values have been filtered out.
Linear-like trends for the spike-in proteins for Condition 1 were generated using the 4,4,10,25 and 50 fmol samples of the original
UPS1 spike-in data. Linear-like trends for the spike-in proteins for Condition 2 were generated using the 2,4,4,10 and 25 fmol samples
of the original UPS1 spike-in data. For more information about the generation of the semi-simulated spike-in datasets,
see the original RolDE publication. In the spike-in datasets, the UPS1
spike-in proteins are expected to differ between conditions while the expression of the rest
of the proteins (the background proteins) are expected to remain stable between the conditions,
excluding experimental noise.
Usage
data3
Format
A matrix with 1033 rows and 30 variables.
Source
https://www.ebi.ac.uk/pride/archive/projects/PXD002099
Let’s run RolDE for data3 using the Bonferroni procedure for multiple
hypothesis adjustement:
set.seed(1)
data3.res<-RolDE(data=data3, des_matrix=des_matrix3, n_cores=3, sig_adj_meth = "bonferroni")
Retrieve the final RolDE results and organize the results based on
strength of differential expression:
RolDE.data3<-data3.res$RolDE_Results
RolDE.data3<-RolDE.data3[order(as.numeric(RolDE.data3[,2])),]
head(RolDE.data3, 3)
| P41159ups|LEP_HUMAN_UPS |
P41159ups|LEP_HUMAN_UPS |
5.235409 |
0 |
0 |
| P01133ups|EGF_HUMAN_UPS |
P01133ups|EGF_HUMAN_UPS |
5.473704 |
0 |
0 |
| P01112ups|RASH_HUMAN_UPS |
P01112ups|RASH_HUMAN_UPS |
7.894447 |
0 |
0 |
Explore the distribution of the estimated significance values:
hist(RolDE.data3$`Estimated Significance Value`, main = "Estimated significance values", xlab="")
We can now observe, that the distribution of signficance values is
different than for the random data of data1, where the null hypothesis
was true. In data3, true differential expression between the conditions
exists. The spike-in proteins are expected to change between the
conditions, while most of the background proteins are expected to remain
stable between the conditions. However, in reality this is not always
the case - some background proteins might be changing too due to
variations in the experimental conditions during the preparation of the
dataset.
Let’s see how RolDE has ranked the spike-in proteins known to be
changing between the conditions:
grep("ups", RolDE.data3[,1])
#> [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 21 22 23 24 25 26 27
#> [26] 29 30 31 32 33 34 35 36 37 38 39 41 42 43 44 46 47 48 54 57 60 67
Most of the spike-in proteins are located near the top of the result
list, as expected.
How many of the proteins remain signifcant at the alpha level of 0.05
after the Bonferroni procedure to adjust for multiple hypothesis
testing:
length(which(RolDE.data3[,4]<=0.05))
#> [1] 74
Findings from the DE analysis can be plotted with the
plotFindings function included in the RolDE package:
| plotFindings | R Documentation |
Plot RolDE results
Description
Plot the findings from longitudinal differential expression analysis with RolDE.
Usage
plotFindings(file_name = NULL, RolDE_res, top_n, col1 = "blue", col2 = "red")
Arguments
file_name |
a string indicating the file name in which the results should be plotted. Should have a ".pdf" extension. Default is NULL, no file is created.
|
RolDE_res |
the RolDE result object.
|
top_n |
an integer or a vector of integers indicating what top differentially expressed features should be plotted. If top_n is a single number, the top_n most
differentially expressed feature will be plotted (e.g top_n=1 will plot the most differentially expressed feature). If top_n is a vector of numbers,
the differentially expressed features corresponding to top detections within the given range will be plotted (e.g. top_n=seq(1:50) will plot the top 50 differentially expressed features).
If more than one feature will be plotted, it is advisable to define a suitable file name in file_name.
|
col1 |
a string indicating which color should be used for Individuals / Replicates in condition 1. The default is blue.
|
col2 |
a string indicating which color should be used for Individuals / Replicates in condition 2. The default is red.
|
Details
The function plots the longitudinal expression of the top RolDE findings. The function can plot either the expression of a single finding
or multiple top findings as indicated by the top_n. The findings can be plotted into a pdf file as indicated by the file_name.
The given file_name should have a ".pdf" extension. If the plottable feature has missing values, a mean value over the feature values will
be imputted for visualization purposes. The missing / imputed value will be indicated with an empty circle symbol.
Value
plotFindings Plots the results from the RolDE object.
Examples
data("res3")
#Plotting the most DE finding. DE results are in the res3 object.
plotFindings(file_name = NULL, RolDE_res = res3, top_n = 1)
Changing the settings for RolDE
While RolDE performs typically very well with the default parameter
values, the user might sometimes wish to alter some of these values for
a specific kind of analysis. Some of the most important parameters
include:
Parameter min_comm_diff controls how many common time points
must two replicates (individuals) have in different conditions to be
compared. The first value controls the number of common time points for
the RegROTS module, while the second one controls the number of common
time points for the DiffROTS module. If min_comm_diff is set to
“auto”, RolDE will use a value of 3 for the RegROTS module and a value
of 1 for the DiffROTS module. In the case of data with non-aligned time
points, the first value of min_comm_diff controls how many time
values (or similar, e.g. age, temperature) must both replicates
(individuals) in different conditions have in the common time interval
to be compared.
Defining larger values for the minimum number of common time points
to be used in RolDE with data1:
set.seed(1)
data1.res<-RolDE(data=data1, des_matrix=des_matrix1, doPar=T, n_cores=3, min_comm_diff = c(4,3))
Using the above values for min_comm_diff requires the
replicates in different conditions to have 4 common time points to be
compared in the RegROTS module and 3 common time points to be compared
in the DiffROTS module.
min_feat_obs controls the number of non-missing values a
feature must have for a replicate (an individual) in a condition to be
compared in the RegROTS module and the DiffROTS module (in data with
aligned time points). A feature is required to have at least
min_feat_obs non-missing values for both replicates in the
different conditions to be compared. If lowered, more missing values are
allowed but the analysis may become less accurate. If increased, only
replicates (individuals) with less missing values for a feature are
compared.
The user can control the degree of polynomials used by the RegROTS
and the PolyReg modules via the degtree_RegROTS and the
degree_PolyReg parameters. If left to “auto”, RolDE will
determine the suitable degrees automatically as described in (Valikangas
et al.)
Using RolDE with non default user given polynomial degrees for the
RegROTS and PolyReg modules:
set.seed(1)
data1.res<-RolDE_Main(data=data1, des_matrix=des_matrix1, n_cores=3, degree_RegROTS = 2, degree_PolyReg = 3)
By default, RolDE uses fixed effects only regression with a common
intercept and slope for the replicates (individuals) when time points in
the data are aligned and mixed effects models with a random effect for
the individual baseline (intercept) if the time points are non aligned
for the PolyReg and the DiffROTS (only in data with non aligned time
points) modules. This behaviour is controlled with the parameter
model_type and the default behaviour is induced when
model_type is allowed to be “auto”. However, the user can
choose to use mixed effects regression modelling when appropriate by
setting the parameter model_type as “mixed0” for random effects
for the individual baseline and setting model_type as “mixed1”
for an individual baseline and slope. Fixed effects only models can be
chosen to be used by setting as “fixed”. Valid inputs for are “auto”
(the default), “mixed0”, “mixed1” and “fixed”.
Analyzing data1 using mixed effects modelling with random intercepts
for the replicates in the PolyReg module of RolDE:
set.seed(1)
data1.res<-RolDE_Main(data=data1, des_matrix=des_matrix1, n_cores=3, model_type="mixed0")
Analyzing data1 using mixed effects modelling with random intercepts
and also random linear slopes for the replicates in the PolyReg AND the
DiffROTS modules of RolDE:
set.seed(1)
data1.res<-RolDE(data=data1, des_matrix=des_matrix1, n_cores=3, model_type="mixed1")
Altering the model_type parameter has an effect for the
DiffROTS module only when the time points in the data are non-aligned.
In non-aligned time point data, the expression level differences between
the conditions in the DiffROTS module is examined when accounting for
time-associated trends of varying complexity in the data.