Sample workflow
The premise behind sigQC is to input primarily two things -
the gene signatures and the datasets to test them upon. Once these
inputs are prepared in the correct format, the
make_all_plots function can be run, with settings as per
user needs, generating the quality control plots needed to ascertain how
each signature performs on each dataset considered.
Below, we consider the case of two random datasets and random gene signatures, and show how to customise each of the inputs.
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:
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.
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.
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.
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.
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.’