## ----setup, include = FALSE--------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", crop = NULL, # Related to https://stat.ethz.ch/pipermail/bioc-devel/ # 2020-April/016656.html dpi = 96, fig.width = 6, fig.height = 4, fig.keep = "all", max.height = "300px" ) local({ hook_output <- knitr::knit_hooks$get("output") knitr::knit_hooks$set(output = function(x, options) { if (!is.null(options$max.height)) { options$attr.output <- c( options$attr.output, sprintf('style="max-height: %s;"', options$max.height) ) } hook_output(x, options) }) }) ## ----install, eval = FALSE---------------------------------------------------- # if (!requireNamespace("BiocManager", quietly = TRUE)) { # install.packages("BiocManager") # } # # BiocManager::install("TiDEomics") ## ----load-packages------------------------------------------------------------ library(TiDEomics) library(SummarizedExperiment) library(org.Mm.eg.db) ## ----load-example------------------------------------------------------------- data("tutorial_sample_info", package = "TiDEomics") data("tutorial_data", package = "TiDEomics") ## ----create-input------------------------------------------------------------- data_obj <- create_input( data = tutorial_data, sample_ann = tutorial_sample_info) data_obj ## ----view-data---------------------------------------------------------------- assays(data_obj)[["orig"]] |> as.data.frame() |> utils::head() ## ----coldata------------------------------------------------------------------ colData(data_obj) # sample annotation ## ----palette------------------------------------------------------------------ custom_palette <- c( "untreated" = "#1b9e77", "IFNbeta" = "#d95f02", "IFNgamma" = "#7570b3", "LPS" = "#e7298a" ) set_custom_palette(custom_palette) ## ----abundance---------------------------------------------------------------- plot_distribution(data_obj, facet_by = "Group") ## ----norm0-------------------------------------------------------------------- data_obj <- normalise_to_start(data_obj) ## ----split-------------------------------------------------------------------- data_obj_list <- split_groups(data_obj) ## ----merge-rep---------------------------------------------------------------- data_obj_merged_list <- merge_replicates(data_obj_list) ## ----merge-groups------------------------------------------------------------- data_obj_merged <- merge_groups(data_obj_merged_list) ## ----cormat------------------------------------------------------------------- plot_cor_matrix(data_obj, method = "spearman", label_rep = TRUE, label_batch = TRUE, cellwidth = 2, cellheight = 2 ) ## ----pca---------------------------------------------------------------------- PC <- plot_pca(data_obj, # pc1 = 1, pc2 = 2, # default to plot PC1 and PC2 plot_screeplot = TRUE, plot_loadings = FALSE, plot_morepc = TRUE, circle = FALSE, morepc = 1:5 ) PC$p_list ## ----pca-3D------------------------------------------------------------------- plot_pca_3D(PC$pca, pcs = 1:3) ## ----pca-eigencor------------------------------------------------------------- PCAtools::eigencorplot(PC$pca, metavars = c("Group", "Time"), components = paste0("PC", 1:5), col = colorRampPalette(c("#3C5488FF", "white", "#E64B35FF"))(100), colCorval = "black" ) ## ----umap--------------------------------------------------------------------- umap <- plot_umap(data_obj, seed = 1234) umap$p_list ## ----pca-umap-group, fig.width = 8, fig.height = 2.5-------------------------- plot_pca_by_group(data_obj, circle = TRUE, arrow = TRUE, legend_pos = "top") # also accepts a list: plot_pca_by_group(data_obj_list) plot_umap_by_group(data_obj, seed = 1234, legend_pos = "top") # also accepts a list: plot_umap_by_group(data_obj_list) ## ----de-between-time, fig.width = 8, fig.height = 2, warning = FALSE---------- DE_between_time_out <- DE_between_time(data_obj, assay = 1, filter = 1, trend = FALSE) ## ----de-between-time-all------------------------------------------------------ DE_between_time_out$all_list$IFNbeta$`t2-t0` |> utils::head() ## ----de-between-time-de------------------------------------------------------- DE_between_time_out$de_list$IFNbeta$`t2-t0` |> utils::head() ## ----de-between-time-plot, fig.width = 7, fig.height = 2---------------------- # Filtered with thresholds in `DE_between_time()` plot_DE_between_time(DE_between_time_out, fontsize = 8, value = FALSE, nrow = 1, heatmap_width = 3 ) # Re-filtering with new thresholds plot_DE_between_time(DE_between_time_out, fontsize = 8, value = FALSE, nrow = 1, heatmap_width = 3, adjP_thres = 0.01, logFC_thres = 1 ) ## ----de-between-group--------------------------------------------------------- DE_between_group_out <- DE_between_group(data_obj, assay = 2, filter = 1, trend = TRUE) ## ----de-between-group-plot, fig.width = 3, fig.height = 2, out.width = "70%"---- # Filtered with thresholds in `DE_between_group()` plot_DE_between_group(DE_between_group_out) # Re-filtering with new thresholds plot_DE_between_group(DE_between_group_out, adjP_thres = 0.01, logFC_thres = 1) ## ----de-between-group-all----------------------------------------------------- DE_between_group_out$all_list$`IFNgamma-untreated`$`24` |> utils::head() ## ----de-between-group-de------------------------------------------------------ DE_between_group_out$de_list$`IFNgamma-untreated` |> utils::head() ## ----volcano------------------------------------------------------------------ plot_volcano(DE_between_group_out, group1 = "untreated", group2 = "IFNgamma", time = 24, logFC_thres = 0.5, adjP_thres = 0.05, label = TRUE) ## ----property----------------------------------------------------------------- data_obj_merged_list <- calc_feature_property(data_obj_merged_list, threshold = 0) property_tb <- summarise_feature_property(data_obj_merged_list) utils::head(property_tb) ## ----group-specific----------------------------------------------------------- group_specific_features(property_tb, groups = c("untreated"), genename = FALSE, GO = FALSE ) ## ----impute------------------------------------------------------------------- data_obj_merged_imp_list <- impute_groups(data_obj_merged_list) ## ----trendy-run--------------------------------------------------------------- set.seed(1234) random_features <- sample(rownames(data_obj_merged_imp_list[[1]]), 50) example_res_list <- run_Trendy(data_obj_merged_imp_list, feature = random_features, minExp = 0.5, maxK = 1, minNumInSeg = 2, meanCut = 0, NCores = 2 ) ## ----trendy-plot, fig.height = 2.5, fig.width = 8----------------------------- plot_segments(data_obj_merged_imp_list, example_res_list, feature = c("Slc25a51", "Aunip"), # example features ylab = "Log2(CPM + 1)" ) ## ----trendy-summary, fig.height = 2, fig.width = 4---------------------------- plot_breakpoints(example_res_list) trendy_summary <- summarise_Trendy(example_res_list) trendy_summary |> utils::head() ## ----trendy-list-------------------------------------------------------------- trendy_list <- extract_segment_trends(trendy_summary) trendy_list$IFNbeta ## ----variance-decomposition, fig.height = 3----------------------------------- # filter genes for variance decomposition: # at least 50% values > 0 in at least 2 groups decomp_filter_genes <- group_specific_features(property_tb, filter_ratio = 0.5, group_pct = 2 / 4, GO = FALSE, genename = FALSE )$features var_decomp <- decomp_variance(data_obj, features = decomp_filter_genes, assay = 1, core = 2 ) plot_variance(var_decomp, rank = "Time", top_n = 20) plot_variance(var_decomp, rank = "Group", top_n = 20) ## ----wgcna-prepare------------------------------------------------------------ # Example filtering by residual variance < Q3 var_res_q3 <- stats::quantile(var_decomp$Residual, 0.75, na.rm = TRUE) filter_wgcna <- var_decomp |> dplyr::filter(Residual < var_res_q3) |> dplyr::pull(Feature) data_obj_wgcna <- data_obj[filter_wgcna, ] ## ----choose-power------------------------------------------------------------- wgcna_input <- prepare_WGCNA(data_obj_wgcna, assay = 2, powers = seq(1, 20), networkType = "signed", RsquaredCut = 0.8 ) wgcna_input$plot picked_power <- wgcna_input$powerEstimate picked_power ## ----run-wgcna---------------------------------------------------------------- net <- run_WGCNA(wgcna_input, power = picked_power, # corType = "pearson", # other option is "bicor" numericLabels = TRUE ) plot_WGCNA(net, fontsize = 8) ## ----extract-wgcna-modules---------------------------------------------------- gene_module <- WGCNA_module(net, exclude_grey = TRUE) gene_module |> dplyr::group_by(Module) |> dplyr::summarise(n = dplyr::n()) ## ----module-metrics----------------------------------------------------------- summarise_module_metrics(net) ## ----plot-wgcna-modules, fig.width = 6, fig.height = 6------------------------ plot_modules_v(gene_module, data_obj_merged, scale = TRUE, ylabel = "Z-score of log2(CPM + 1)", height_ratio = 2 ) ## ----go_rank, fig.width=5, fig.height=3--------------------------------------- gse_group <- enrichGO_rank(var_decomp, gene_rank_by = "Group", OrgDb = org.Mm.eg.db, keyType = "SYMBOL", category = "BP", go_rank_by = "p.adjust") enrichplot::gseaplot2(gse_group[["BP"]], geneSetID = 1:3, base_size = 8) ## ----go-example, fig.width = 6------------------------------------------------ background_wgcna <- colnames(wgcna_input$data) go_list <- enrichGO_list( gene_list = gene_module, OrgDb = org.Mm.eg.db, universe = background_wgcna, pvalueCutoff = 0.9, # get more results for demonstration qvalueCutoff = 0.9, simplify = FALSE, keyType = "SYMBOL" ) plot_GO(go_list$all, plot_dotplot = TRUE, plot_emapplot = FALSE, plot_cnetplot = FALSE, showCategory_dotplot = 3) ## ----enrich-msigdb------------------------------------------------------------ # Mouse Hallmark gene sets hallmark_msigdb <- enrich_msigdb(gene_module, universe = background_wgcna, minGSSize = 5, category = "MH", species = "Mus musculus", db_species = "MM") # Human C2 (chemical and genetic perturbations) perturb_msigdb <- enrich_msigdb(gene_module, universe = background_wgcna, category = "C2", subcategory = "CGP", species = "Mus musculus") # Specific gene sets (category not required) ifn_msigdb <- enrich_msigdb(gene_module, universe = background_wgcna, gene_sets = c("HALLMARK_INTERFERON_ALPHA_RESPONSE", "HALLMARK_INTERFERON_GAMMA_RESPONSE"), minGSSize = 5, species = "Mus musculus", db_species = "MM") ## ----plot-wgcna-modules-go, fig.width = 8, fig.height = 5--------------------- plot_modules_h(gene_module, data_obj_merged, scale = TRUE, ylabel = "Z-score of log2(CPM + 1)", enrich_list = go_list$all, enrich_category = "BP", heatmap_width = 6, heatmap_height = 8 ) ## ----plot-wgcna-modules-multi, fig.width = 15, fig.height = 5----------------- hub_features <- extract_hubs(net, top_n = 3) plot_modules_h(gene_module, data_obj_merged, scale = TRUE, ylabel = "Z-score of log2(CPM + 1)", enrich_list = c(go_list$all, hallmark_msigdb), enrich_category = c("BP", "CC", "MH"), enrich_p_threshold = c(0.05, 0.05, 0.05), mark_features = hub_features, heatmap_width = 6, heatmap_height = 8 ) ## ----plot-feature------------------------------------------------------------- plot_trend(data_obj, assay = 1, features = c("Abtb1", "Dram1", "Ifi27", "Nufip1"), title = "Example features") # Or pre-calculate mean and sd with calc_mean_sd() table_mean_sd_orig <- calc_mean_sd(data_obj)$orig plot_trend(table_mean_sd_orig, features = c("Abtb1", "Dram1", "Ifi27", "Nufip1"), title = "Example features") ## ----plot-residual------------------------------------------------------------ plot_trend(data_obj, assay = 1, features = var_decomp |> dplyr::arrange(Residual) |> utils::head(12) |> dplyr::pull(Feature), title = "Features with lowest residual variance" ) plot_trend(data_obj, assay = 1, features = var_decomp |> dplyr::arrange(dplyr::desc(Residual)) |> utils::head(12) |> dplyr::pull(Feature), title = "Features with highest residual variance" ) ## ----prepare-tide------------------------------------------------------------- tide <- prepare_tide( data = tutorial_data, sample_ann = tutorial_sample_info, filter_ratio = 0.5, min_groups = 2, keep = "below_quantile", residual_threshold = 0.75 ) tide$filter_summary ## ----tide-de------------------------------------------------------------------ tide$DE <- list( between_group = DE_between_group_out, between_time = DE_between_time_out ) tide$WGCNA <- net tide$enrichment$GO_modules <- go_list tide$enrichment$MSigDB <- hallmark_msigdb tide$enrichment$GO_group <- gse_group ## ----isee, eval = FALSE------------------------------------------------------- # iSEE::iSEE(data_obj) ## ----deedee------------------------------------------------------------------- # Flatten nested results de_flat <- flatten_DE(tide$DE) enrich_flat <- flatten_enrich(tide$enrichment) # Build DeeDeeExperiment dde <- DeeDeeExperiment::DeeDeeExperiment( sce = tide$se, de_results = de_flat, enrich_results = enrich_flat ) ## ----session-info------------------------------------------------------------- sessionInfo()