Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
30 changes: 24 additions & 6 deletions R/calZscore.R
Original file line number Diff line number Diff line change
Expand Up @@ -20,12 +20,18 @@ calZscore <- function(...) UseMethod("calZscore")

calZscore.default <- function(x, mad = T) {
if (mad) {
dev = mad(x, na.rm = T)
dev <- mad(x, na.rm = T)
} else {
dev = sd(x, na.rm = T)
dev <- sd(x, na.rm = T)
}
x = round((x-median(x, na.rm = T))/dev, 3)
x
# Avoid division by zero
if (dev == 0) {
warning("Standard deviation is zero. Returning original values.")
return(x)
}

x <- round((x-median(x, na.rm = T))/dev, 3)
return(x)
}

#' Calculate z-scores of QC metrics of a SampleDataset
Expand All @@ -52,6 +58,7 @@ calZscore.sampleDataset <- function (
object$zscore = .zscoreColNames(qcMetrics)
object$df[object$zscore] = NA


if (is.null(strat)) {
# calculate z-score across all qcMetrics
object$df[, object$zscore] = apply(object$df[, qcMetrics], 2, calZscore)
Expand All @@ -75,15 +82,26 @@ calZscore.sampleDataset <- function (
return(object)
}

#' Calculate maximum vcf z-score and update the object with a maxVcfzscore
#' column
#' Calculate maximum VCF z-score and update the object with a maxVcfZscore column
#'
#' @param object a sampleDataset object
#' @param qcMetrics a vector of names of QC metrics
#' @return updated sampleDataset object with a new maxVcfZscore column
#'
#' @export
.calMaxVcfZ <- function(object, qcMetrics) {
vcfQcMetrInUse = .zscoreColNames(intersect(object$vcfQcMetr, qcMetrics))
object$df$maxVcfzscore = apply(abs(object$df[vcfQcMetrInUse]), 1, max)
return(object)
}


#' Get z-score column names
#'
#' @param qcMetrics a vector of QC metric names
#' @return a vector of z-score column names
#'
#' @export
.zscoreColNames <- function (qcMetrics) {
return(paste(qcMetrics, 'Zscore', sep = ''))
}
Expand Down
12 changes: 0 additions & 12 deletions R/examples.R

This file was deleted.

6 changes: 3 additions & 3 deletions R/flagSamples.R
Original file line number Diff line number Diff line change
Expand Up @@ -66,12 +66,12 @@ flagSamples.sampleDataset <- function(sds, cutoffs, zscore = NULL){
stop("None of metrics in the cutoff table is in the SampleDataset,
please double check your input data.")
}
if (!is.null(zscore) && !('zscore' %in% attributes(sds)$names)) {
if (!is.null(zscore) && !('zscore' %in% names(attributes(sds)))) {
stop("Input dataset has no z-score, please calculate z-score before
applying z-score filters.")
}

sds$df$flaggedReason = ''
sds$df$flaggedReason <- ''
for (i in 1:dim(cutoffs)[1]) {
# to do: re-write with sapply
sds$df = flagSamples(sds$df, cutoffs$qcMetrics[i], cutoffs$value[i],
Expand All @@ -84,6 +84,6 @@ flagSamples.sampleDataset <- function(sds, cutoffs, zscore = NULL){
sds$df = flagSamples(sds$df, i, -zscore, greater = F)
}
}
sds['flaggedReason'] = 'flaggedReason'
sds['flaggedReason'] <- 'flaggedReason'
return(sds)
}
61 changes: 29 additions & 32 deletions R/plot.R
Original file line number Diff line number Diff line change
Expand Up @@ -19,27 +19,24 @@ sampleQcPlot.default <- function(
main = 'QC', geom = c('scatter', 'violin', 'hist')
) {
geom <- match.arg(geom)
if (!is.null(annot)) {
if (geom == 'scatter') {
# scatter plot stratified by sample
plt = scatter(data = data, x = 'index', y = qcMetric, strat=annot,
xlab = 'samples', legend = legend, main = main,
outliers = outliers, primaryID = sds$primaryID)
}
if (geom == 'violin') {
data[[annot]] = factor(.toSameLength(data[[annot]]))
plt <- (ggplot2::ggplot(data, ggplot2::aes_string(annot, qcMetric,
if (!is.null(annot) && geom == 'scatter') {
# scatter plot stratified by sample
plt <- scatter(data = data, x = 'index', y = qcMetric, strat=annot,
xlab = 'samples', legend = legend, main = main,
outliers = outliers, primaryID = primaryID)
}else if(!is.null(annot) && geom == 'violin'){
data[[annot]] = factor(.toSameLength(data[[annot]]))
plt <- (ggplot2::ggplot(data, ggplot2::aes_string(annot, qcMetric,
color = annot))
+ ggplot2::geom_violin()
+ ggplot2::geom_jitter(height = 0, width = 0.3))
}
} else {
if (geom == 'hist') {
plt = ggplot2::qplot(data[[qcMetric]], geom = 'histogram', bins = 100,
xlab = qcMetric)
}
+ ggplot2::geom_violin()
+ ggplot2::geom_jitter(height = 0, width = 0.3))
}else if(geom == 'hist'){
plt <- ggplot2::qplot(data[[qcMetric]], geom = 'histogram', bins = 100,
xlab = qcMetric)
}
if(!is.null(plt)){
plt <- plt + ggplot2::ggtitle(main)
}
plt = plt + ggplot2::ggtitle(main)
return(plt)
}

Expand All @@ -65,17 +62,17 @@ sampleQcPlot.sampleDataset <- function(
ncols = 5, show = FALSE, sort = TRUE
) {
if(length(qcMetrics) == 1) { ncols = 1 }
if (!is.null(outliers)) {
if(!all(outliers %in% sds$df[[sds$primaryID]])) {
stop("Not all outliers are in the Sample Dataset. Please double check.")
}
if (!is.null(outliers) && length(setdiff(outliers, sds$df[[sds$primaryID]])) > 0) {
stop("Not all outliers are in the Sample Dataset. Please double check.")
}
geom <- match.arg(geom)
position <- match.arg(position)
if(is.null(qcMetrics)) {
qcMetrics = sds$qcMetrics
qcMetrics <- sds$qcMetrics
}
if (sort) sds = sort(sds, by = annot)
if (sort) sds <- sort(sds, by = annot)

geom <- match.arg(geom)
position <- match.arg(position)

plots = sapply(
qcMetrics,
function(x) sampleQcPlot(
Expand All @@ -102,29 +99,29 @@ outlierPlots <- function(...) UseMethod('outlierPlots')
outlierPlots.default <- function(tab, qcMetrics, strat, main, outliers,
primaryID=NULL, type='violin'){
plots <- list()
plots[[1]] = scatter(tab, x = 'sample', y = qcMetrics, strat = strat,
plots[[1]] <- scatter(tab, x = 'sample', y = qcMetrics, strat = strat,
xlab = 'sample', outliers = outliers, main = main,
legend = T, primaryID = primaryID)
if (type == 'density') {
plots[[2]] = (ggplot2::qplot(tab[, qcMetrics], color = tab[, strat],
plots[[2]] <- (ggplot2::qplot(tab[, qcMetrics], color = tab[, strat],
geom = "density", main= main, xlab = qcMetrics)
+ geom_vline(linetype="dashed",
xintercept = tab[outlier.index, qcMetrics]))
} else {
plots[[2]] = (ggplot2::ggplot(data = tab,
plots[[2]] <- (ggplot2::ggplot(data = tab,
ggplot2::aes_string(x=strat, y=qcMetrics,
color = strat))
+ ggplot2::geom_violin()
+ ggplot2::geom_jitter(height = 0, width = 0.3, alpha = 0.3)
+ ggplot2::ggtitle(main)
+ ggplot2::theme(axis.text.x = ggplot2::element_blank()))
if(!is.null(outliers)) {
plots[[2]] = (plots[[2]]
plots[[2]] <- (plots[[2]]
+ ggplot2::geom_point(data=tab[tab$sampleId %in% outliers,],
colour = 'black', size = 3))
}
}
multiplot(plotlist = plots, ncols = 2)
return(multiplot(plotlist = plots, ncols = 2))
}

#' Produce outlier plots for a sampleDataset
Expand Down
2 changes: 1 addition & 1 deletion R/runExplorer.R
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
runQcExplorer <- function () {
appDir <- system.file("apps", "qcExplorer", package = "samplyzer")
if (appDir == "") {
stop("Could not find example directory. Try re-installing `samplyzer`.",
stop("Could not find qcExplorer directory. Try re-installing `samplyzer`.",
call. = FALSE)
}
shiny::runApp(appDir, display.mode = "normal")
Expand Down
38 changes: 20 additions & 18 deletions R/sampleDataset.R
Original file line number Diff line number Diff line change
@@ -1,21 +1,22 @@
.loadInput <- function (input) {
if (is.data.frame(input)) {
df = input
df <- input
} else if (is.character(input)) {
df = read.csv(input, sep='\t')
df <- read.csv(input, sep='\t')
} else {
stop("Input is not a tsv file or a data.frame")
}
return(df)
}

.mergeMetr <- function(df, metrDf, primaryID) {
df = merge(df, metrDf, by = primaryID)
df <- merge(df, metrDf, by = primaryID)
return(df)
}

.getNames <- function(df, primaryID) {
df = names(df)[-which(names(df) == primaryID)]
return(names(df)[-which(names(df) == primaryID)])
#df = names(df)[-which(names(df) == primaryID)]
}

#' Sample dataset object
Expand Down Expand Up @@ -53,37 +54,38 @@ sampleDataset <- function(

# import annotation
if (is.null(df)) {
df = .loadInput(annotInput)
annotNames = .getNames(df, primaryID)
df <- .loadInput(annotInput)
annotNames <- .getNames(df, primaryID)

if (!all(stratify %in% annotNames)) {
stop("Not all stratified factors are in the annotation file,
please double check your input stratify option.")
}
bamQcMetrName = NULL
vcfQcMetrName = NULL

bamQcMetrName <- NULL
vcfQcMetrName <- NULL
if (!is.null(bamQcInput)) {
bamQcMetr = .loadInput(bamQcInput)
df = .mergeMetr(df, bamQcMetr, primaryID)
bamQcMetrName = .getNames(bamQcMetr, primaryID)
bamQcMetr <- .loadInput(bamQcInput)
df <- .mergeMetr(df, bamQcMetr, primaryID)
bamQcMetrName <- .getNames(bamQcMetr, primaryID)
}
if (!is.null(vcfQcInput)) {
vcfQcMetr = .loadInput(vcfQcInput)
df = .mergeMetr(df, vcfQcMetr, primaryID)
vcfQcMetrName = .getNames(vcfQcMetr, primaryID)
vcfQcMetr <- .loadInput(vcfQcInput)
df <- .mergeMetr(df, vcfQcMetr, primaryID)
vcfQcMetrName <- .getNames(vcfQcMetr, primaryID)
}
} else {
# directly create SDS from an already defined data.frame
bamQcMetrName = bamQcInput
vcfQcMetrName = vcfQcInput
annot = annotInput
bamQcMetrName <- bamQcInput
vcfQcMetrName <- vcfQcInput
annot <- annotInput
}

sds = list(df = df, qcMetrics = c(bamQcMetrName, vcfQcMetrName),
bamQcMetr = bamQcMetrName, vcfQcMetr = vcfQcMetrName,
annot = annotNames, primaryID = primaryID)

sds$df$index = 1:length(df[[primaryID]])
sds$df$index <- 1:length(df[[primaryID]])
class(sds) <- 'sampleDataset'
return(sds)
}