DatasetExperiment objects

The DatasetExperiment object is an extension of the SummarizedExperiment class used by the Bioconductor community. It contains three main parts:

1. data A data frame containing the measured data for each sample.
2. sample_meta A data frame of additional information related to the samples e.g. group labels.
3. variable_meta A data frame of additional information related to the variables (features) e.g. annotations

Like all struct objects it also contains name and description fields (called “slots” is R language).

A key difference between DatasetExperiment and SummarizedExperiment objects is that the data is transposed. i.e. for DatasetExperiment objects the samples are in rows and the features are in columns, while the opposite is true for SummarizedExperiment objects.

All slots are accessible using dollar notation.

# show some data
head(D$data[,1:4])

##   Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1          5.1         3.5          1.4         0.2
## 2          4.9         3.0          1.4         0.2
## 3          4.7         3.2          1.3         0.2
## 4          4.6         3.1          1.5         0.2
## 5          5.0         3.6          1.4         0.2
## 6          5.4         3.9          1.7         0.4

Statistical models

Before we can apply e.g. PCA we first need to create a PCA object. This object contains all the inputs, outputs and methods needed to apply PCA. We can set parameters such as the number of components when the PCA model is created, but we can also use dollar notation to change/view it later.

P = PCA(number_components=15)
P$number_components=5
P$number_components

## [1] 5

The inputs for a model can be listed using param_ids(object):

param_ids(P)

## [1] "number_components"

Or a summary of the object can be printed to the console:

P

## A "PCA" object
## --------------
## name:          Principal Component Analysis (PCA)
## description:   PCA is a multivariate data reduction technique. It summarises the data in a smaller number of
##                  Principal Components that maximise variance.
## input params:  number_components 
## outputs:       scores, loadings, eigenvalues, ssx, correlation, that 
## predicted:     that
## seq_in:        data

Model sequences

Unless you have good reason not to, it is usually sensible to mean centre the columns of the data before PCA. Using the STRUCT framework we can create a model sequence that will mean centre and then apply PCA to the mean centred data.

M = mean_centre() + PCA(number_components = 4)

In structToolbox mean centring and PCA are both model objects, and joining them using “+” creates a model_sequence object. In a model_sequence the outputs of the first object (mean centring) are automatically passed to the inputs of the second object (PCA), which allows you chain together modelling steps in order to build a workflow.

The objects in the model_sequence can be accessed by indexing, and we can combine this with dollar notation. For example, the PCA object is the second object in our sequence and we can access the number of components as follows:

M[2]$number_components

## [1] 4

Training/testing models

Model and model_sequence objects need to be trained using data in the form of a DatasetExperiment object. For example, the PCA model sequence we created (M) can be trained using the iris DatasetExperiment object (‘D’).

M = model_train(M,D)

This model sequence has now mean centred the original data and calculated the PCA scores and loadings.

Model objects can be used to generate predictions for test datasets. For the PCA model sequence this involves mean centring the test data using the mean of training data, and the projecting the centred test data onto the PCA model using the loadings. The outputs are all stored in the model sequence and can be accessed using dollar notation. For this example we will just use the training data again (sometimes called autoprediction), which for PCA allows us to explore the training data in more detail.

M = model_predict(M,D)

Sometimes models don’t make use the training/test approach e.g. univariate statsitics, filtering etc. For these models the model_apply method can be used instead. For models that do provide training/test methods, model_apply applies autoprediction by default i.e. it is a short-cut for applying model_train and model_predict to the same data.

M = model_apply(M,D)

The available outputs for an object can be listed and accessed like input params, using dollar notation:

output_ids(M[2])

## [1] "scores"      "loadings"    "eigenvalues" "ssx"         "correlation"
## [6] "that"

M[2]$scores

## A "DatasetExperiment" object
## ----------------------------
## name:          
## description:   
## data:          150 rows x 4 columns
## sample_meta:   150 rows x 2 columns
## variable_meta: 4 rows x 1 columns

Model charts

The struct framework includes chart objects. Charts associated with a model object can be listed.

chart_names(M[2])

## [1] "pca_biplot"           "pca_correlation_plot" "pca_dstat_plot"      
## [4] "pca_loadings_plot"    "pca_scores_plot"      "pca_scree_plot"

Like model objects, chart objects need to be created before they can be used. Here we will plot the PCA scores plot for our mean centred PCA model.

C = pca_scores_plot(factor_name='class') # colour by class
chart_plot(C,M[2])

Note that indexing the PCA model is required because the pca_scores_plot object requires a PCA object as input, not a model_sequence.

If we make changes to the input parameters of the chart, chart_plot must be called again to see the effects.

# add petal width to meta data of pca scores
M[2]$scores$sample_meta$example=D$data[,1]
# update plot
C$factor_name='example'
chart_plot(C,M[2])

The chart_plot method returns a ggplot object so that you can easily combine it with other plots using the gridExtra or cowplot packages for example.

# scores plot
C1 = pca_scores_plot(factor_name='class') # colour by class
g1 = chart_plot(C1,M[2])

# scree plot
C2 = pca_scree_plot()
g2 = chart_plot(C2,M[2])

# arange in grid
grid.arrange(grobs=list(g1,g2),nrow=1)

Ontology

Within the struct framework (and structToolbox) an ontology slot is provided to allow for standardardised definitions for objects and its inputs and outputs using the Ontology Lookup Service (OLS).

For example, STATO is a general purpose STATistics Ontology (http://stato-ontology.org). From the webpage:

Its aim is to provide coverage for processes such as statistical tests, their conditions of application, and information needed or resulting from statistical methods, such as probability distributions, variables, spread and variation metrics. STATO also covers aspects of experimental design and description of plots and graphical representations commonly used to provide visual cues of data distribution or layout and to assist review of the results.

The ontology for an object can be set by assigning the ontology term identifier to the ontology slot of any struct_class object at design time. The ids can be listed using $ notation:

# create an example PCA object
P=PCA()

# ontology for the PCA object
P$ontology

## [1] "OBI:0200051"

The ontology method can be used obtain more detailed ontology information. When cache = NULL the struct package will automatically attempt to use the OLS API (via the rols package) to obtain a name and description for the provided identifiers. Here we used cached versions of the ontology definitions provided in the structToolbox package to prevent issues connecting to the OLS API when building the package.

ontology(P,cache = ontology_cache()) # set cache = NULL (default) for online use

## [[1]]
## An object of class "ontology_list"
## Slot "terms":
## [[1]]
## term id:       OBI:0200051
## ontology:      obi
## label:         principal components analysis dimensionality reduction
## description:   A principal components analysis dimensionality reduction is a dimensionality reduction
##                  achieved by applying principal components analysis and by keeping low-order principal
##                  components and excluding higher-order ones.
## iri:           http://purl.obolibrary.org/obo/OBI_0200051
## 
## 
## 
## [[2]]
## An object of class "ontology_list"
## Slot "terms":
## [[1]]
## term id:       STATO:0000555
## ontology:      stato
## label:         number of predictive components
## description:   number of predictive components is a count used as input to the principle component analysis
##                  (PCA)
## iri:           http://purl.obolibrary.org/obo/STATO_0000555

Note that the ontology method returns definitions for the object (PCA) and the inputs/outputs (number_of_components).

Cross-validation

Cross validation is a common technique for assessing the performance of classification models. For this example we will use a Partial least squares-discriminant analysis (PLS-DA) model. Data should be mean centred prior to PLS, so we will build a model sequence first.

M = mean_centre() + PLSDA(number_components=2,factor_name='class')
M

## A model_seq object containing:
## 
## [1]
## A "mean_centre" object
## ----------------------
## name:          Mean centre
## description:   The mean sample is subtracted from all samples in the data matrix. The features in the centred
##                  matrix all have zero mean.
## input params:  mode 
## outputs:       centred, mean_data, mean_sample_meta 
## predicted:     centred
## seq_in:        data
## 
## [2]
## A "PLSDA" object
## ----------------
## name:          Partial least squares discriminant analysis
## description:   PLS is a multivariate regression technique that extracts latent variables maximising
##                  covariance between the input data and the response. The Discriminant Analysis
##                  variant uses group labels in the response variable. For >2 groups a 1-vs-all
##                  approach is used. Group membership can be predicted for test samples based on
##                  a probability estimate of group membership, or the estimated y-value.
## input params:  number_components, factor_name, pred_method 
## outputs:       scores, loadings, yhat, design_matrix, y, reg_coeff, probability, vip, pls_model, pred, threshold, sr, sr_pvalue 
## predicted:     pred
## seq_in:        data

iterator objects like the k-fold cross-validation object (kfold_xval) can be created just like any other struct object. Parameters can be set at creation using the equals sign, and accessed or changed later using dollar notation.

# create object
XCV = kfold_xval(folds=5,factor_name='class')

# change the number of folds
XCV$folds=10
XCV$folds

## [1] 10

The model to be cross-validated can be set/accessed using the models method.

models(XCV)=M
models(XCV)

## A model_seq object containing:
## 
## [1]
## A "mean_centre" object
## ----------------------
## name:          Mean centre
## description:   The mean sample is subtracted from all samples in the data matrix. The features in the centred
##                  matrix all have zero mean.
## input params:  mode 
## outputs:       centred, mean_data, mean_sample_meta 
## predicted:     centred
## seq_in:        data
## 
## [2]
## A "PLSDA" object
## ----------------
## name:          Partial least squares discriminant analysis
## description:   PLS is a multivariate regression technique that extracts latent variables maximising
##                  covariance between the input data and the response. The Discriminant Analysis
##                  variant uses group labels in the response variable. For >2 groups a 1-vs-all
##                  approach is used. Group membership can be predicted for test samples based on
##                  a probability estimate of group membership, or the estimated y-value.
## input params:  number_components, factor_name, pred_method 
## outputs:       scores, loadings, yhat, design_matrix, y, reg_coeff, probability, vip, pls_model, pred, threshold, sr, sr_pvalue 
## predicted:     pred
## seq_in:        data

Alternatively, iterators can be combined with models using the multiplication symbol was shorthand for the models assignement method:

# cross validation of a mean centred PLSDA model
XCV = kfold_xval(
        folds=5,
        method='venetian',
        factor_name='class') * 
      (mean_centre() + PLSDA(factor_name='class'))

The run method can be used with any iterator object. The iterator will then run the set model or model sequence multiple times.

In this case we run cross-validation 5 times, splitting the data into different training and test sets each time.

The run method also needs a metric to be specified, which is another type of struct object. This metric may be calculated once after all iterations, or after each iteration, depending on the iterator type (resampling, permutation etc). For cross-validation we will calculate “balanced accuracy” after all iterations.

XCV = run(XCV,D,balanced_accuracy())
XCV$metric

##              metric mean sd
## 1 balanced_accuracy 0.11 NA

Note The balanced_accuracy metric actually reports 1-accuracy, so a value of 0 indicates perfect performance. The standard deviation “sd” is NA in this example because there is only one permutation.

Like other struct objects, iterators can have chart objects associated with them. The chart_names function will list them for an object.

chart_names(XCV)

## [1] "kfoldxcv_grid"   "kfoldxcv_metric"

Charts for iterator objects can be plotted in the same way as charts for any other object.

C = kfoldxcv_grid(
  factor_name='class',
  level=levels(D$sample_meta$class)[2]) # first level
chart_plot(C,XCV)

It is possible to combine multiple iterators by using the multiplication symbol. This is equivalent to nesting one iterator inside the other. For example, we can repeat our cross-validation multiple times by permuting the sample order.

# permute sample order 10 times and run cross-validation
P = permute_sample_order(number_of_permutations = 10) * 
    kfold_xval(folds=5,factor_name='class')*
    (mean_centre() + PLSDA(factor_name='class',number_components=2))
P = run(P,D,balanced_accuracy())
P$metric

##              metric   mean          sd
## 1 balanced_accuracy 0.1095 0.004972145

Training and Test sets

Splitting data into training and test sets is an important aspect of machine learning. In structToolbox this is implemented using the split_data object for random subsampling across the whole dataset, and stratified_split for splitting based on group sizes, which is the approach used by Mendez et al.

# prepare model
M = stratified_split(p_train=0.75,factor_name='Class')
# apply to filtered data
M = model_apply(M,filtered)
# get data from object
train = M$training
train

## A "DatasetExperiment" object
## ----------------------------
## name:          Gastric cancer (NMR)
##                (Training set)
## description:   • 1H-NMR urinary metabolomic profiling for diagnosis of gastric cancer
##                • A subset of the data has been selected as a training set
## data:          62 rows x 53 columns
## sample_meta:   62 rows x 5 columns
## variable_meta: 53 rows x 1 columns

cat('\n')

test = M$testing
test

## A "DatasetExperiment" object
## ----------------------------
## name:          Gastric cancer (NMR)
##                (Testing set)
## description:   • 1H-NMR urinary metabolomic profiling for diagnosis of gastric cancer
##                • A subset of the data has been selected as a test set
## data:          21 rows x 53 columns
## sample_meta:   21 rows x 5 columns
## variable_meta: 53 rows x 1 columns

Optimal number of PLS components

In Mendez et al a k-fold cross-validation is used to determine the optimal number of PLS components. 100 bootstrap iterations are used to generate confidence intervals. In strucToolbox these are implemented using “iterator” objects, that can be combined with model objects. R2 is used as the metric for optimisation, so the PLSR model in structToolbox will be used. For speed only 10 bootstrap iterations are used here.

# scale/transform training data
M = log_transform(base = 10) +
    autoscale() + 
    knn_impute(neighbours = 3,by='samples')

# apply model
M = model_apply(M,train)

# get scaled/transformed training data
train_st = predicted(M)

# prepare model sequence
MS = grid_search_1d(
        param_to_optimise = 'number_components',
        search_values = as.numeric(c(1:6)),
        model_index = 2,
        factor_name = 'Class_num',
        max_min = 'max') *
     permute_sample_order(
        number_of_permutations = 10) *
     kfold_xval(
        folds = 5,
        factor_name = 'Class_num') * 
     (mean_centre(mode='sample_meta')+
      PLSR(factor_name='Class_num'))

# run the validation
MS = struct::run(MS,train_st,r_squared())

#
C = gs_line()
chart_plot(C,MS)

The chart plotted shows Q2, which is comparable with Figure 13 of Mendez et al . Two components were selected by Mendez et al, so we will use that here.

PLS model evalutation

To evaluate the model for discriminant analysis in structToolbox the PLSDA model is appropriate.

# prepare the discriminant model
P = PLSDA(number_components = 2, factor_name='Class')

# apply the model
P = model_apply(P,train_st)

# charts
C = plsda_predicted_plot(factor_name='Class',style='boxplot')
g1 = chart_plot(C,P)

C = plsda_predicted_plot(factor_name='Class',style='density')
g2 = chart_plot(C,P)+xlim(c(-2,2))

C = plsda_roc_plot(factor_name='Class')
g3 = chart_plot(C,P)

plot_grid(g1,g2,g3,align='vh',axis='tblr',nrow=1, labels=c('A','B','C'))

# AUC for comparison with Mendez et al
MET = calculate(AUC(),P$y$Class,P$yhat[,1])
MET

## A "AUC" object
## --------------
## name:          Area under ROC curve
## description:   The area under the ROC curve of a classifier is estimated using the trapezoid method.
## value:         0.9739583

Note that the default cutoff in A and B of the figure above for the PLS models in structToolbox is 0, because groups are encoded as +/-1. This has no impact on the overall performance of the model.

Permutation test

A permutation test can be used to assess how likely the observed result is to have occurred by chance. In structToolbox permutation_test is an iterator object that can be combined with other iterators and models.

# model sequence
MS = permutation_test(number_of_permutations = 20,factor_name = 'Class_num') * 
     kfold_xval(folds = 5,factor_name = 'Class_num') *
     (mean_centre(mode='sample_meta') + PLSR(factor_name='Class_num', number_components = 2))

# run iterator
MS = struct::run(MS,train_st,r_squared())

# chart
C = permutation_test_plot(style = 'density') 
chart_plot(C,MS) + xlim(c(-1,1)) + xlab('R Squared')

This plot is comparable to the bottom half of Figure 17 in Mendez et. al.. The unpermuted (true) Q2 values are consistently better than the permuted (null) models. i.e. the model is reliable.

PLS projection plots

PLS can also be used to visualise the model and interpret the latent variables.

# prepare the discriminant model
P = PLSDA(number_components = 2, factor_name='Class')

# apply the model
P = model_apply(P,train_st)

C = pls_scores_plot(components=c(1,2),factor_name = 'Class')
chart_plot(C,P)

PLS feature importance

Regression coefficients and VIP scores can be used to estimate the importance of individual features to the PLS model. In Mendez et al bootstrapping is used to estimate the confidence intervals, but for brevity here we will skip this.

# prepare chart
C = pls_vip_plot(ycol = 'HE')
g1 = chart_plot(C,P)

C = pls_regcoeff_plot(ycol='HE')
g2 = chart_plot(C,P)

plot_grid(g1,g2,align='hv',axis='tblr',nrow=2)

Data transformation

The data is log2 transformed, then scaled such that the mean of the medians is equal for all conditions. These steps are available in structToolbox using log_transform and mean_of_medians objects.

# prepare model sequence
M = log_transform(
  base=2) + 
  mean_of_medians(
    factor_name = 'condition')
# apply model sequence
M = model_apply(M,DS) 

# get transformed data
DST = predicted(M)

The Reporter genes are plotted again for comparison.

# chart object
C = feature_boxplot(feature_to_plot=Ldha,factor_name='time',label_outliers=FALSE)
g1=chart_plot(C,DST)+ggtitle('Ldha')+ylab('log2(expression)')
C = feature_boxplot(feature_to_plot=Hk2,factor_name='time',label_outliers=FALSE)
g2=chart_plot(C,DST)+ggtitle('Hk2')+ylab('log2(expression)')

plot_grid(g1,g2,nrow=1,align='vh',axis='tblr')

Missing value filtering

Missing value filtering involves removing any feature (gene) where there are at least 3 missing values per group in at least 11 samples.

This specific filter is not in structToolbox at this time, but can be achieved by combining filter_na_count and filter_by_name objects.

Specifically, the default output of filter_na_count is changed to return a matrix of NA counts per class. This output is then connected to the ‘names’ input of filter_by_names and converted to TRUE/FALSE using the ‘seq_fcn’ input.

The ‘seq_fcn’ function processes the NA counts before they are used as inputs for filter_by_names. When data is passed along the model sequence it passes unchanged through the filter_na_count object becuase the default output has been changed, so the filter_na_count and filter_by_name objects are working together as a single filter.

# build model sequence
M2 = filter_na_count(
  threshold=2,
  factor_name='condition', 
  predicted='na_count') +    # override the default output
  filter_by_name(
    mode='exclude',
    dimension='variable',
    names='place_holder',
    seq_in='names',
    seq_fcn=function(x) {      # convert NA count pre group to true/false
      x=x>2                  # more the two missing per group
      x=rowSums(x)>10        # in more than 10 groups
      return(x)
    }
  )
# apply to transformed data
M2 = model_apply(M2,DST)

# get the filtered data
DSTF = predicted(M2)

Missing value imputation

STATegra uses two imputation methods that are not available as struct objects, so we create temporary STATegra_impute objects to do this using some functions from the struct package.

The first imputation method imputes missing values for any treatment where values are missing for all samples using a “random value below discovery”. We create a new struct object using set_struct_obj in the global environment, and a “method_apply” method that implements the imputation.

# create new imputation object
set_struct_obj(
  class_name = 'STATegra_impute1',
  struct_obj = 'model',
  params=c(factor_sd='character',factor_name='character'),
  outputs=c(imputed='DatasetExperiment'),
  prototype = list(
    name = 'STATegra imputation 1',
    description = 'If missing values are present for all one group then they are replaced with min/2 + "random value below discovery".',
    predicted = 'imputed'
  ) 
)

# create method_apply for imputation method 1
set_obj_method(
  class_name='STATegra_impute1',
  method_name='model_apply',
  definition=function(M,D) {
    
    # for each feature count NA within each level
    na = apply(D$data,2,function(x){
      tapply(x,D$sample_meta[[M$factor_name]],function(y){
        sum(is.na(y))
      })
    })
    # count number of samples in each group
    count=summary(D$sample_meta[[M$factor_name]])
    # standard deviation of features within levels of factor_sd
    sd = apply(D$data,2,function(x) {tapply(x,D$sample_meta[[M$factor_sd]],sd,na.rm=TRUE)})
    sd = median(sd,na.rm=TRUE)
    
    # impute or not
    check=na == matrix(count,nrow=2,ncol=ncol(D)) # all missing in one class
    
    # impute matrix
    mi = D$data
    for (j in 1:nrow(mi)) {
      # index of group for this sample
      g = which(levels(D$sample_meta[[M$factor_name]])==D$sample_meta[[M$factor_name]][j])
      iv=rnorm(ncol(D),min(D$data[j,],na.rm=TRUE)/2,sd)
      mi[j,is.na(mi[j,]) & check[g,]] = iv[is.na(mi[j,]) & check[g,]]
    }
    D$data = mi
    M$imputed=D
    return(M)
  }
)

The second imputation method replacing missing values in any condition with exactly 1 missing value with the mean of the values for that condition. Again we create a new struct object and corresponding method for the the new object to implement the filter.

# create new imputation object
set_struct_obj(
  class_name = 'STATegra_impute2',
  struct_obj = 'model',
  params=c(factor_name='character'),
  outputs=c(imputed='DatasetExperiment'),
  prototype = list(
    name = 'STATegra imputation 2',
    description = 'For those conditions with only 1 NA impute with the mean of the condition.',
    predicted = 'imputed'
  ) 
)

# create method_apply for imputation method 2
set_obj_method(
  class_name='STATegra_impute2',
  method_name='model_apply',
  definition=function(M,D) {
    # levels in condition
    L = levels(D$sample_meta[[M$factor_name]])
    
    # for each feature count NA within each level
    na = apply(D$data,2,function(x){
      tapply(x,D$sample_meta[[M$factor_name]],function(y){
        sum(is.na(y))
      })
    })
    
    # standard deviation of features within levels of factor_sd
    sd = apply(D$data,2,function(x) {tapply(x,D$sample_meta[[M$factor_name]],sd,na.rm=TRUE)})
    sd = median(sd,na.rm=TRUE)
    
    # impute or not
    check=na == 1 # only one missing for a condition
    
    # index of samples for each condition
    IDX = list()
    for (k in L) {
      IDX[[k]]=which(D$sample_meta[[M$factor_name]]==k)
    }
    
    ## impute
    # for each feature
    for (k in 1:ncol(D)) {
      # for each condition
      for (j in 1:length(L)) {
        # if passes test
        if (check[j,k]) {
          # mean of samples in group
          m = mean(D$data[IDX[[j]],k],na.rm=TRUE)
          # imputed value
          im = rnorm(1,m,sd)
          # replace NA with imputed
          D$data[is.na(D$data[,k]) & D$sample_meta[[M$factor_name]]==L[j],k]=im
        }
      }
    }
    M$imputed=D
    return(M)
  }
)

The new STATegra imputation objects can now be used in model sequences like any other struct object. A final filter is added to remove any feature that has missing values after imputation.

# model sequence
M3 = STATegra_impute1(factor_name='treatment',factor_sd='condition') +
     STATegra_impute2(factor_name = 'condition') + 
     filter_na_count(threshold = 3, factor_name='condition')
# apply model
M3 = model_apply(M3,DSTF)
# get imputed data
DSTFI = predicted(M3)
DSTFI

## A "DatasetExperiment" object
## ----------------------------
## name:          STATegra Proteomics
## description:   downloaded from: https://github.com/STATegraData/STATegraData/
## data:          36 rows x 864 columns
## sample_meta:   36 rows x 4 columns
## variable_meta: 864 rows x 3 columns

Exploratory analysis

It is often useful to visualise the distribution of values across samples to verify that the transformations/normalisation/filtering etc have been effective.

The values are no longer skewed and show an approximately normal distribution. The boxplots are comparable in width with very few outliers indicated, so the transformations etc have had an overall positive effect.

PCA is used to provide a graphical representation of the data. For comparison with the outputs from STATegra a filter is included to reduce the data to include only the treated samples (IKA)

# model sequence
P = filter_smeta(mode='include',factor_name='treatment',levels='IKA') +
    mean_centre() + 
    PCA(number_components = 2)
# apply model
P = model_apply(P,DSTFI)

# scores plots coloured by factors
g = list()
for (k in c('batch','time')) {
  C = pca_scores_plot(factor_name=k,ellipse='none')
  g[[k]]=chart_plot(C,P[3])
}

plot_grid(plotlist = g,nrow=1)

There does not appear to be a strong batch effect. PC1 is dominated by time point “24” and some potentially outlying points from time points “2” and “0”.

Data preprocessing

In the STATegra project the LCMS data was combined with GCMS data and multiblock analysis was conducted. Here only the LCMS will be explored, so the data will be processed differently in comparison to Gomez-Cabrero et al. Some basic processing steps will be applied in order to generate a valid PCA plot from the biological and QC samples.

# prepare model sequence
MS = filter_smeta(mode = 'include', levels='QC', factor_name = 'sample_type') +
     knn_impute(neighbours=5) +
     vec_norm() + 
     log_transform(base = 10) 

# apply model sequence
MS = model_apply(MS, DE)

## Warning in knnimp(x, k, maxmiss = rowmax, maxp = maxp): 3 rows with more than 50 % entries missing;
##  mean imputation used for these rows

Exploratory analysis

First we will use PCA to look at the QC samples in order to make an assessment of the data quality.

# pca model sequence
M =  mean_centre() +
     PCA(number_components = 3)

# apply model
M = model_apply(M,predicted(MS))


# PCA scores plot
C = pca_scores_plot(factor_name = 'sample_type',label_factor = 'order',points_to_label = 'all')
# plot
chart_plot(C,M[2])

The QC labelled “36” is clearly very different to the other QCs. In STATegra this QC was removed, so we will exclude it here as well. This corresponds to QC H1. STATegra also excluded QC samples measured immediately after a blank, which we will also do here.

# prepare model sequence
MS = filter_smeta(
      mode = 'include', 
      levels='QC', 
      factor_name = 'sample_type') +
  
     filter_by_name(
      mode = 'exclude', 
      dimension='sample',
      names = c('1358BZU_0001QC_H1','1358BZU_0001QC_A1','1358BZU_0001QC_G1')) +
  
     knn_impute(
      neighbours=5) +
  
     vec_norm() + 
  
     log_transform(
       base = 10) + 
  
     mean_centre() +
  
     PCA(
       number_components = 3)

# apply model sequence
MS = model_apply(MS, DE)

## Warning in knnimp(x, k, maxmiss = rowmax, maxp = maxp): 4 rows with more than 50 % entries missing;
##  mean imputation used for these rows

# PCA scores plot
C = pca_scores_plot(factor_name = 'sample_type',label_factor = 'order',points_to_label = 'all')
# plot
chart_plot(C,MS[7])

Now we will plot the QC samples in context with the samples. There are several possible approaches, and we will apply the approach of applying PCA to the full dataset including the QCs. We will exclude the blanks as it is likely that they will dominate the plot if not removed. All samples from batch 12 were excluded from STATegra and we will replicate that here.

# prepare model sequence
MS = filter_smeta(
      mode = 'exclude', 
      levels='Blank', 
      factor_name = 'sample_type') +
  
     filter_smeta(
      mode = 'exclude', 
      levels='12', 
      factor_name = 'biol.batch') +
  
     filter_by_name(
      mode = 'exclude', 
      dimension='sample',
      names = c('1358BZU_0001QC_H1',
                '1358BZU_0001QC_A1',
                '1358BZU_0001QC_G1')) +
  
     knn_impute(
      neighbours=5) +
  
     vec_norm() + 
  
     log_transform(
       base = 10) + 
  
     mean_centre() +
  
     PCA(
       number_components = 3)

# apply model sequence
MS = model_apply(MS, DE)

## Warning in knnimp(x, k, maxmiss = rowmax, maxp = maxp): 2 rows with more than 50 % entries missing;
##  mean imputation used for these rows

# PCA scores plots
C = pca_scores_plot(factor_name = 'sample_type')
# plot
chart_plot(C,MS[8])

The QCs appear to representative of the samples, but there are strong clusters in the data, including the QC samples which have no biological variation. There is likely to be a number of ‘low quality’ features that should be excluded, so we will do that now, and use more sophisticated normalisation (PQN) and scaling methods (glog).

MS =  filter_smeta(
       mode = 'exclude', 
       levels = '12', 
       factor_name = 'biol.batch') +
  
      filter_by_name(
       mode = 'exclude', 
       dimension='sample',
       names = c('1358BZU_0001QC_H1',
                 '1358BZU_0001QC_A1',
                 '1358BZU_0001QC_G1')) +

      blank_filter(
       fold_change = 20,
       qc_label = 'QC',
       factor_name = 'sample_type') +

      filter_smeta(
       mode='exclude',
       levels='Blank',
       factor_name='sample_type') +
  
      mv_feature_filter(
       threshold = 80, 
       qc_label = 'QC', 
       factor_name = 'sample_type', 
       method = 'QC') +
     
      mv_feature_filter(
        threshold = 50, 
        factor_name = 'sample_type', 
        method='across') +
  
     rsd_filter(
       rsd_threshold=20, 
       qc_label='QC',
       factor_name='sample_type') +
  
     mv_sample_filter(
       mv_threshold = 50) +
     
     pqn_norm(
       qc_label='QC',
       factor_name='sample_type') +
     
     knn_impute(
       neighbours=5, 
       by='samples') +
     
     glog_transform(
       qc_label = 'QC',
       factor_name = 'sample_type') +
     
     mean_centre() + 
     
     PCA(
       number_components = 10)

# apply model sequence
MS = model_apply(MS, DE)


# PCA plots using different factors
g=list()
for (k in c('order','biol.batch','time.point','condition')) {
  C = pca_scores_plot(factor_name = k,ellipse='none')
  # plot
  g[[k]]=chart_plot(C,MS[length(MS)])
}

plot_grid(plotlist = g,align='vh',axis='tblr',nrow=2,labels=c('A','B','C','D'))

We can see now that the QCs are tightly clustered. This indicates that the biological variance of the remaining high quality features is much greater than the technical variance represented by the QCs.

There does not appear to be a trend by measurement order (A), which is an important indicator that instrument drift throughout the run is not a large source of variation in this dataset.

There does not appear to be strong clustering related to biological batch (B).

There does not appear to be a strong trend with time (C) but this is likely to be a more subtle variation and might be masked by other sources of variance at this stage.

There is some clustering related to condition (D) but with overlap.

To further explore any trends with time, we will split the data by the condition factor and only explore the Ikaros group. Removing the condition factor variation will potentially make it easier to spot any more subtle trends. We will extract the glog transformed matrix from the previous model sequence and continue from there.

# get the glog scaled data
GL = predicted(MS[11])

# extract the Ikaros group and apply PCA
IK = filter_smeta(
      mode='include',
      factor_name='condition',
      levels='Ikaros') +
     mean_centre() + 
     PCA(number_components = 5)

# apply the model sequence to glog transformed data
IK = model_apply(IK,GL)

# plot the PCA scores
C = pca_scores_plot(factor_name='time.point',ellipse = 'sample')
g1=chart_plot(C,IK[3])
g2=g1 + scale_color_viridis_d() # add continuous scale colouring

plot_grid(g1,g2,nrow=2,align='vh',axis = 'tblr',labels=c('A','B'))

Colouring by groups (A) makes the time point trend difficult to see, but by adding a ggplot continuous colour scale “viridis” (B) the trend with time along PC1 becomes much clearer.