Calculation of z-scores

zscores.from.metrics takes quality control data for each sample and calculates the z-score across each metric. Each row of qc.data should correspond to a sample, and each column to a test metric. An example input dataframe can be seen below.

zscores <- zscores.from.metrics(qc.data = example.qc.dataframe);

The function returns a dataframe containing the z-scores for each sample and test metric. The example data in this package would return the dataframe below.

Adjusting metric directionality

This package is designed to be independent of test metrics used, but since some test metrics are better if they are larger and some are better if they are smaller, the sign of a poor z-score will differ between the metrics. We must therefore adjust the metrics such that negative z-scores is considered a bad measurement across all tests. The function correct.zscore.signs requires input that states if a positive or a negative z-score is good for each metric. This input takes the form of a dataframe with a column for metric and a column for sign which have the potential values ‘pos’ or ‘neg’ for each metric used.

The function correct.zscore.signs also sets all non-negative z-scores to zero to make sure positive and negative values do not cancel each other out when we calculate an accumulated score.

Apart from the zscores, and signs.data, the function also takes the names of the columns containing the metric name and the sign instructions in sign.data.

zscores.corrected <- correct.zscore.signs(
  zscores = zscores,
  signs.data = sign.correction,
  metric.col.name = 'Metric',
  signs.col.name = 'Sign'
  );

The output of the function, zscores.corrected, can now be used to calculate accumulated scores for each sample.

Calculating total quality score

accumulate.zscores takes the zscores.corrected dataframe and calculates a total quality score for each sample. It does this by summing over all negative z-scores. It then orders the data by the magnitude of the quality score and returns the result in a dataframe.

quality.scores <- accumulate.zscores(zscores.corrected = zscores.corrected);

The resulting quality.scores is shown below:

Finding the best fitting distribution

The function fit.and.evaluate takes the quality.scores and evaluates how well they fit to some common distributions. It returns a Bayesian Information Criterion score as well as a Kolmogorov–Smirnov test result for each distribution.

fit.results <- fit.and.evaluate(
    quality.scores = quality.scores,
    trim.factor = 0.15
    );

fit.and.evaluate returns the dataframe printed below:

The iterative method for outlier nomination using cosine similarity

The function cosine.similarity.iterative takes quality.scores, trims the proportion of data indicated by trim.factor from each extreme, and fits it to the distribution given. It then tests the largest datapoint compared a null distribution of size no.simulations. If the largest datapoint has a significant p-value, it moves onto the 2nd largest datapoint and so on, until it reaches a datapoint whose p-value is insignificant.

outlier.detect.iterative.res <- cosine.similarity.iterative(
    quality.scores = quality.scores,
    distribution = 'lnorm',
    no.simulations = 1000,
    trim.factor = 0.15,
    alpha.significant = 0.05
    );

Number of outliers found:

print(outlier.detect.iterative.res$no.outliers);
#> [1] 11

The sample labels of the outliers:

print(outlier.detect.iterative.res$outlier.labels);
#>  [1] "CPCG0375" "CPCG0382" "CPCG0498" "CPCG0464" "CPCG0437" "CPCG0235"
#>  [7] "CPCG0486" "CPCG0263" "CPCG0266" "CPCG0256" "CPCG0502"

The cutoff method for outlier nomination using cosine similarity

The function cosine.similarity.cutoff takes quality.scores, trims the proportion of data indicated by trim.factor from each extreme, and fits it to the distribution given. It then simulates as many datasets as stated by no.simulations, and computes the cosine similarity of each dataset against theoretical distribution. Using this simulatd dataset, it estimates the cutoff threshold for cosine similarity which would correspond to a statistically significant p-value (determined by alpha.significant) and nominates outliers that surpass this threshold.

outlier.detect.cutoff.res <- cosine.similarity.cutoff(
    quality.scores = quality.scores,
    distribution = 'lnorm',
    no.simulations = 1000,
    trim.factor = 0.15,
    alpha.significant = 0.05
    );

Quality score cutoff:

print(outlier.detect.cutoff.res$cutoff);
#>      5% 
#> 24.8339

Number of outliers found:

print(outlier.detect.cutoff.res$no.outliers);
#> [1] 6

The sample labels of the outliers:

print(outlier.detect.cutoff.res$outlier.labels);
#> [1] "CPCG0235" "CPCG0437" "CPCG0464" "CPCG0498" "CPCG0382" "CPCG0375"

Quality score barplot

The function get.qc.barplot takes the accumulated quality scores (quality.scores) and returns a barplot displaying the scores from lowest to highest. It can also optionally depict a cut-off for outlier nomination. If filename is not NULL, it is saved to file. Otherwise, the trellis object is returned.

qc.barplot <- get.qc.barplot(
    quality.scores = quality.scores,
    abline.h = - outlier.detect.cutoff.res$cutoff
    );

Z-score heatmap

The function get.qc.heatmap takes the dataframe with the corrected z-scores for each sample (zscores.corrected) and metric as well as a vector of labels (ylabels) for the y-axis. The labels should be the full names of the metrics in the same order as they are in the dataframe. The function also takes quality.scores to make sure the samples are in the correct order on the heatmap. The function returns a standardized heatmap. If filename is not NULL, it is saved to file. Otherwise, the trellis object is returned.

qc.heatmap <- get.qc.heatmap(
  zscores = zscores.corrected,
  quality.scores = quality.scores,
  yaxis.lab = ylabels
  );

Aggregating the plots

Now that both the heat map and barplot has been generated they can be aggregated to a multipanelplot. get.qc.multipanelplot takes the barplot, heatmap and a filename and concatenates the plots in a standardized format and saves the multipanelplot to file.

qc.multipanel <- get.qc.multipanelplot(
  barplot = qc.barplot,
  heatmap = qc.heatmap
  );