Step 1: Preparation of input data

Input data should consist of lists of expression matrices, and should be pre-normalised, batch-corrected, and standardised if required. Care should be taken to ensure that genes of interest are present in the dataset and not reported primarily as NA values. Gene signatures must be annotated in a manner consistent with the input data.

dataset_1 <- replicate(1000, rnorm(500,mean=0,sd=1))#random matrix - 1000 genes x 500 samples
row.names(dataset_1) <- as.character(1:(dim(dataset_1)[1])) #set the gene names
colnames(dataset_1) <- paste0('d1_',as.character(1:(dim(dataset_1)[2]))) #set the sample IDs

dataset_2 <- replicate(1000, rnorm(250,mean=0,sd=2))#random matrix - 1000 genes x 250 samples
row.names(dataset_2) <- as.character(1:(dim(dataset_2)[1])) #set the gene names
colnames(dataset_2) <- paste0('d2_',as.character(1:(dim(dataset_2)[2]))) #set the sample IDs

#gene sets that we are going to use for our gene signatures
gene_sig_1 <- c('1','2','3','5','8','13','21','34','55','89','144','233','377','610','987') 
gene_sig_2 <- c('2','4','8','16','32','64','128','256','512')

Any specific expression thresholds for expression (other than global median) should be computed, as this is the default the package uses as an expression cutoff.

Step 2: Creation of input variables

Next, we load the input datasets and gene signatures into list format.

mRNA_expr_matrix = list()
mRNA_expr_matrix[["dataset_1"]] = dataset_1
mRNA_expr_matrix[["dataset_2"]] = dataset_2

gene_sigs_list = list()
gene_sigs_list[['gene_sig_1']] = as.matrix(gene_sig_1)
gene_sigs_list[['gene_sig_2']] = as.matrix(gene_sig_2)

We then set the other non-default valued input variables as needed for our analysis:

showResults <- FALSE # we do not want to show the reuslts in R graphics windows
doNegativeControl <- FALSE # we do not want to compute the negative or permutation controls for time purposes

Step 3: Running of sigQC package

library("sigQC")

make_all_plots(gene_sigs_list, mRNA_expr_matrix, showResults = showResults, doNegativeControl = doNegativeControl)

This will produce, in the desired output directory, a number of plots in PDF files which may be analysed as described in the subsequent steps. The package also creates an output file `log.log’ in the output directory, a text file, which summarises the run, and reports any errors that may have occurred if they are not printed to the console. This should be consulted if any issues are encountered in the running of this principal function and for troubleshooting purposes.

Step 4: Analysis of expression

Expression of signature genes should be evaluated across samples in all datasets, and this is done by analysis of the plots sig_expr_*.pdf.

These plots describe the proportion of samples with supra-threshold expression of each signature gene (either for the median threshold, or user-specified as above), and the proportion of samples with non-NA values, identifying non-expressed signature components.

Barchart of expression of signature genes above threshold

Step 5: Analysis of variability

To check signature gene variability, load the file ‘sig_mean_vs_sd.pdf.’ These plots describe the mean and standard deviation of expression of all genes reported (in grey) versus all signature genes (in red), with corresponding dashed lines over the plots describing the 10th, 25th, 50th, 75th and 90th percentiles of both mean and standard deviation.

This facilitates the easy identification of those signature genes, which are not variable or expressed sufficiently among the samples, as well as a global evaluation of signature behaviour across samples of a dataset.

Plot of mean vs. standard deviation of signature genes’ expression (red) and all other genes (grey)

Step 6: Analysis of co-correlation of scoring metrics

Next, loading the files called ’sig_compare_metrics_*.pdf’ it is possible to analyse the co-correlation of mean, median and first principal component (PCA1) as scoring metrics across the samples for each signature across each dataset.

Also shown in the fourth row of panels of these plots is a principal components analysis (PCA) scree plot, which describes the proportion of the variance attributable to each principal component, which may reflect whether the first PCA represents a reasonable scoring summary metric for a particular gene signature.

Plot of co-correlations and distributions of scoring metrics for gene signature 1 across datasets

In addition, loading the files called ’sig_compare_ES_metrics_*.pdf’, we provide analogous co-correlation plots for each of GSVA, ssGSEA, and PLAGE scores, also showing the distribution of these scoring metrics. These scores are based on the techniques implemented by the GSVA R package, which is a dependency of sigQC, and the original package and publications by Hanzelman et al., 2013, Barbie et al., 2009, and Tomfohr et al., 2005 should be consulted for details of the calculation of these metrics. As these are enrichment socre-based metrics, these are best used in concert with the quality control for signatures based on enrichment analysis results, such as those from MSigDB.

Also, loading the files called ’scoring_metrics_corr_*.pdf’, one can observe the co-correlation between all of the six scoring metrics considered (mean, median, PCA1, GSVA, ssGSEA, and PLAGE) as a heatmap of Spearman correlation coefficients. Again, this is of particular use for signatures derived based on enrichment analyses.

Analysis of normality of scoring metrics

In the files called ’sig_qq_plots_*.pdf.’ we show QQ plots against the normal distribution, to depict how close to normally the signature scores are distributed across the datasets, for the mean, median, and first principal component.

Plot of QQ plots for distributions of scoring metrics for gene signature 1 across datasets

Analysis of modality of scoring metrics

Next, to understand the modality of the distributions of the scoring metrics across the datasets, we use the mclust package to compute the Bayes Information Criterion (BIC) for the likelihood of different Gaussian mixture models. The files called ’sig_gaussian_mixture_model_*.pdf’ provide plots of the modality on the x axis for each type of considered Gaussian mixture model, and the corresponding BIC on the y-axis.

Note that we have also provided a text output interpreting the modality for each scoring metric in each dataset/signature combination in the file ‘mixture_model_out.txt.’

Lastly, we save each of the raw ouptuts from the mclust call in the file ‘mixture_models_raw_out.rda’ in a list variable called mixture_models, with the following structure.

mixture_models[[signature_name]][[dataset_name]][[metric_type]] will contain the mclust output for the Gaussian mixture model fitting for the metric_type (one of ‘median’,‘mean’,or ‘pca1’, referring to the scoring criteria) for the dataset called ‘dataset_name’ and the signature called ‘signature_name.’

Plot of BIC vs modality for Gaussian mixture models of scoring metric distributions across datasets

Step 7: Analysis of data standardisation effects

To study data standardisation effects, load the output file called ‘sig_standardisation_comp.pdf.’. This plot provides the comparison of median of gene signature expression on the raw data provided versus the median of the gene signature expression on the z-transformed (standardised to zero mean and unit variance) dataset, for each sample in each dataset and each gene signature under consideration.

Plot of median vs. Z-transformed median score for gene signature and dataset combinations

Step 8: Analysis of signature compactness

The files produced in the output directory called ‘sig_autocor_hmaps.pdf’ and ‘sig_autocor_dens.pdf’ provide heatmap and kernel density estimate plots of the correlation of each signature genes’ expression with the expression of every other signature gene, providing an analysis of signature compactness.

In addition to these plots produced for the analysis of intra-signature correlation, the files ’sig_autocor_rankProd_*.pdf’ are produced when there is more than one dataset present for a given gene signature (not shown for this analysis). These plots represent the output of the BioConductor RankProduct package in the evaluation of signature genes whose median intra-signature correlation with all other genes consistently ranks low with the other signature genes. This facilitates the process for finding genes with consistent poor intra-signature correlation with the other genes of the gene signature, particularly when refining a given signature for optimal performance across a number of different datasets (e.g. multiple clinical cohorts, or clinical data and cell line data).

Plot of intra-signature correlation heatmaps for gene signature and dataset combinations

Step 9: Analysis of signature structure

Firstly, signature structure can be evaluated by hierarchical clustering on the provided expression values of the signature elements over all samples, in conjunction with annotations for the samples, if they are provided. These plots are present in the output directory and are named ’sig_eval_struct_clustering_*.pdf,’ which are clustered based on each dataset in turn, and run over each signature and each dataset present.

Plot of signature gene expression across samples for gene signature 1

In addition to hierarchical clustering on patient samples and signature genes, in the output directory of {}, biclustering results can be found if there is a case of biclustering found among the sample groups and signature genes are found (not shown for this analysis). The particular details of the biclustering outputs and algorithm used may be found in the R Package documentation. If no biclusters are found, the files ‘sig_eval_bivariate_clustering.pdf’ are not produced.

Step 10: Optional: Comparison of multiple signatures

The file produced, named ‘sig_radarplot.pdf’ describes each signature applied to each dataset in a holistic, radar chart format. This plot evaluates the gene signature across a number of QC metrics, many of which are summary metrics for those in steps 4-9 of this procedure, and these are described in detail in the Table below.

Radar plot of signature metrics summarising results of QC for each dataset/signature combination

Step 11: Optional: Analysis of null distributions of sigQC metrics (resulting in longer runtimes)

The file produced called ‘boxplot_metrics.pdf,’ (not produced in this example) in the negative_control subfolder of the results output shows the distributions of each of the fourteen summary metrics as reported in the radar plot for each signature and dataset combination. These distributions are generated for the number of repeats as specified by the input parameter numResampling, with default set to 50. The values for the gene signature and dataset combination in question are shown in red overlaid with the other points in grey, giving a sense of significance of each of the metrics.

In addition, the file produced called ‘boxplot_metrics.pdf’ (not produced in this example) in the permutation_control subfolder of the results output shows the distributions of each of the fourteen summary metrics for each signature and dataset combination as well. These distributions are produced by permuting the gene labels for each sample, for each of the signature genes present in each dataset. The number of times this process is done is given by the input parameter numResampling, with default set to 50. As with the random gene signatures, the values for the gene signature and dataset combination in question are shown in red overlaid with the other points in grey. These metrics help provide an interpretation of the gene signature’s significance with respect to resampling among the same genes, to generate another set of null distributions that can be used in the evaluation of a gene signature.

Imagemagick dependency

Issues may be experienced if the ‘ImageMagick’ dependency is not installed on the user’s system (particularly for Windows systems). To install this dependency, please follow directions at: https://imagemagick.org/script/download.php.

RStudio plot size issue

There is a known issue with RStudio plot size being too small causing errors with the plots being produced, even if they are not displayed in RStudio. The error in this case will appear in the log file as: Error in graphics::plot.new(): figure margins too large.

In this case, expand the plot size in the RStudio editor, and also click the maximise button on the plot window of RStudio, and re-run the package.

Input data specifications

• Gene signatures must be formatted as a list of matrices, of dimension \(k\) rows by 1 column, for a signature of length \(k\) genes. Inputting a single list as a vector will cause an error to the program.
• Datasets must also be formatted as a list of matrices, such that genes are the rownames of the dataset, and samples are organised by columns of the dataset.
• Gene signatures and datasets must be annotated in the same way, as if the names of the genes of a signature are not found in a dataset, the computation will not continue.
• Care must be taken to ensure that NA valued genes are removed as optimally as possible, as if there are too many values in the expression matrix for the gene signature are NA, calculations dependent upon singular value decomposition (e.g. principal components analysis) cannot be carried out.