## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", message = FALSE, fig.height = 10, fig.width = 10 ) ## ----setup-------------------------------------------------------------------- library(qmtools) library(SummarizedExperiment) library(vsn) library(pls) library(ggplot2) library(patchwork) set.seed(1e8) data(faahko_se) ## Only keep the first assay for the vignette assays(faahko_se)[2:4] <- NULL faahko_se ## ----filtering---------------------------------------------------------------- dim(faahko_se) # 206 features table(colData(faahko_se)$sample_group) ## Missing value filter based on 2 groups. dim(removeFeatures(faahko_se, i = "raw", group = colData(faahko_se)$sample_group, cut = 0.80)) # nothing removed dim(removeFeatures(faahko_se, i = "raw", group = colData(faahko_se)$sample_group, cut = 0.85)) # removed 65 features ## based on "WT" only dim(removeFeatures(faahko_se, i = "raw", group = colData(faahko_se)$sample_group, levels = "WT", cut = 0.85)) ## ----plotMiss, fig.wide = TRUE, fig.height = 5-------------------------------- ## Sample group information g <- factor(colData(faahko_se)$sample_group, levels = c("WT", "KO")) ## Visualization of missing values plotMiss(faahko_se, i = "raw", group = g) ## ----knn, fig.wide = TRUE, fig.height = 5------------------------------------- se <- imputeIntensity(faahko_se, i = "raw", name = "knn", method = "knn") se # The result was stored in assays slot: "knn" ## Standardization of input does not influence the result m <- assay(faahko_se, "raw") knn_scaled <- as.data.frame( imputeIntensity(scale(m), method = "knn") # Can accept matrix as an input ) knn_unscaled <- as.data.frame(assay(se, "knn")) idx <- which(is.na(m[, 1]) | is.na(m[, 2])) # indices for missing values p1 <- ggplot(knn_unscaled[idx, ], aes(x = ko15.CDF, y = ko16.CDF)) + geom_point() + theme_bw() p2 <- ggplot(knn_scaled[idx, ], aes(x = ko15.CDF, y = ko16.CDF)) + geom_point() + theme_bw() p1 + p2 + plot_annotation(title = "Imputed values: unscaled vs scaled") ## ----vsn, fig.wide = TRUE, fig.height = 5------------------------------------- se <- normalizeIntensity(se, i = "knn", name = "knn_vsn", method = "vsn") se # The result was stored in assays slot: "knn_vsn" p1 <- plotBox(se, i = "knn", group = g, log2 = TRUE) # before normalization p2 <- plotBox(se, i = "knn_vsn", group = g) # after normalization p1 + p2 + plot_annotation(title = "Before vs After normalization") ## ----PCA---------------------------------------------------------------------- ## PCA m_pca <- reduceFeatures(se, i = "knn_vsn", method = "pca", ncomp = 2) summary(m_pca) ## ----PLSDA-------------------------------------------------------------------- ## PLS-DA (requires information about each sample's group) m_plsda <- reduceFeatures(se, i = "knn_vsn", method = "plsda", ncomp = 2, y = g) summary(m_plsda) ## ----plotReduced, fig.wide = TRUE, fig.height = 5----------------------------- p_pca <- plotReduced(m_pca, group = g) p_plsda <- plotReduced(m_plsda, label = TRUE, ellipse = TRUE) p_pca + p_plsda + plot_annotation(title = "PCA and PLS-DA") ## ----clusterFeatures---------------------------------------------------------- se <- clusterFeatures(se, i = "knn_vsn", rtime_var = "rtmed", rt_cut = 10, cor_cut = 0.7) rowData(se)[, c("rtmed", "rtime_group", "feature_group")] ## ----rtime hclust, fig.wide = TRUE, fig.height = 5---------------------------- rts <- rowData(se)$rtmed rt_cut <- 10 fit <- hclust(dist(rts, "manhattan")) plot(as.dendrogram(fit), leaflab = "none") rect.hclust(fit, h = rt_cut) ## ----connected, fig.small = TRUE---------------------------------------------- se_connected <- clusterFeatures(se, i = "knn_vsn", rtime_var = "rtmed", rt_cut = 10, cor_cut = 0.7, cor_grouping = "connected") plotRTgroup(se_connected, i = "knn_vsn", group = "FG.22") ## ----louvain, fig.small = TRUE------------------------------------------------ se_louvain <- clusterFeatures(se, i = "knn_vsn", rtime_var = "rtmed", rt_cut = 10, cor_cut = 0.7, cor_grouping = "louvain") plotRTgroup(se_louvain, i = "knn_vsn", group = "FG.22") ## ----pairs, fig.wide = FALSE-------------------------------------------------- plotRTgroup(se_louvain, i = "knn_vsn", group = "FG.22", type = "pairs") ## ----compareSamples----------------------------------------------------------- ## Compute statisticis for the contrast: KO - WT fit <- compareSamples(se, i = "knn_vsn", group = "sample_group", class1 = "WT", class2 = "KO") ## List top 5 features head(fit, 5) ## ----compareSamples covariates------------------------------------------------ ## Include covariates colData(se)$covar <- c(rep(c("A", "B"), 6)) compareSamples(se, i = "knn_vsn", group = "sample_group", covariates = "covar", class1 = "WT", class2 = "KO", number = 5) ## ----session info------------------------------------------------------------- sessionInfo()