--- title: "DgeaHeatmap Vignette" author: "Lancelle, Leonie Johanna" date: "`r Sys.Date()`" output: BiocStyle::html_document: toc: true vignette: > %\VignetteIndexEntry{DgeaHeatmap Vignette} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ```{r setup} library(DgeaHeatmap) ``` # DgeaHeatmap The goal of DgeaHeatmap is to enable R users to generate heatmaps more easily and to help with preprocessing read counts. Furthermore, the package is aimed to simplify the extraction of raw read counts from .dcc and .pkc files generated through Nanostring GeoMx DSP. ## Installation You can install the development version of DgeaHeatmap from [GitLab](https://gitlab.ub.uni-bielefeld.de/spittaulab/Dgea_Heatmap_Package) with: ``` r # install.packages("pak") devtools::install_gitlab(repo = "https://gitlab.ub.uni-bielefeld.de/spittaulab/Dgea_Heatmap_Package.git") ``` ## Usage This is a basic example which shows you how to solve a common problem: ```{r inputnormalizedReadCountsExample, echo=FALSE, results='hide'} input_data <- read.csv(system.file("extdata/testfile_counts.csv", package = "DgeaHeatmap")) ``` ```{r exampleBuildMatrix, echo=TRUE, results='hide'} library(DgeaHeatmap) x <- 1 matrixCounts <- build_matrix(input_data, x) ``` ```{r showTable, echo=FALSE, message=FALSE, warning=FALSE, include=TRUE, out.height="60%", out.width="60%"} knitr::kable(head(matrixCounts[, seq_len(4)]), "pipe") ``` ```{r exampleIndividualMatrix, echo=TRUE} parameter1 <- "DKD" parameter2 <- "glomerulus" factors_for_individual_matrix <- list(parameter1, parameter2) indiMatrix <- individual_matrix(factors_for_individual_matrix, matrixCounts) ``` ```{r showIndividualMatrix, echo=FALSE} knitr::kable(head(indiMatrix[, seq_len(4)]), "pipe") ``` Filtering the matrix for only an x amount of most variable expressed genes: ```{r exampleFilteringForTopExprGenes, echo=TRUE} top_number_of_genes <- 20 varGenesMatrix <- filtering_for_top_exprGenes(indiMatrix, top_number_of_genes) print(nrow(varGenesMatrix)) ``` Using Z-score scaling to scale the values of the matrix: ```{r exampleZScoreScaling, echo=TRUE, results='hide'} scaled_counts <- scale_counts(varGenesMatrix) ``` How to visualize the data distribution of the scaled counts: ```{r exampleShowDataDistribution, echo=TRUE} show_data_distribution(scaled_counts) ``` Generating an elbow plot to choose k for k-mean generation: ```{r exampleElbowPlot, echo=TRUE} seed <- 1 # setting a seed for a reproducible outcome set.seed(seed) # set seed for reproducibility elbow_plot(scaled_counts) ``` This can also be conducted with the transposed matrix to get the number of clusters the columns can be divided into: ```{r exampleElbowPlot_columns, echo=TRUE, warning=FALSE} maxK <- 4 # only 4 samples therefore a max of 4 clusters possible seed <- 1 # setting a seed for a reproducible outcome transposed_matrix <- t(scaled_counts) # transposes matrix set.seed(seed) # set seed for reproducibility elbow_plot(transposed_matrix, maxK) ``` As desired, the samples can further be summarized as biological replicates to help generate more clearly arranged heatmaps and giving a better overview. ```{r exampleSummariseBioReplicates, echo=TRUE, results='hide'} probes <- list("disease3_DKD_glomerulus_Geometric_Segment", "disease1B_DKD_glomerulus_Geometric_Segment", "disease2B_DKD_glomerulus_WT") sumBioRepsMatrix <- summarise_bio_replicates(scaled_counts, probes) ``` ```{r showSummedBioReplicates, echo=FALSE} knitr::kable(head(sumBioRepsMatrix[, seq_len(3)]), "pipe") ``` Generating K-means for clustering in the heatmap: ```{r exampleKmeanGeneration, echo=TRUE} seed <- 4 # set a seed for reproducibility, but try different seeds first k_clusters <- 2 set.seed(seed) K_meanTable <- Kmean_generation(sumBioRepsMatrix, k_clusters) print(K_meanTable) ``` ```{r showGeneratedKmeans, echo=FALSE} knitr::kable(head(K_meanTable[, seq_len(4)]), "pipe") ``` The most variable genes of each cluster are selected through the following function. A new object is created with the purpose of saving the most variable genes for each cluster, while the matrix with included cluster for each gene is divided into individual cluster matrices. These matrices are each sorted by their variance. ```{r exampleMostVariableGenes, echo=TRUE} number_of_annotations_per_cluster <- 2 k_clusters <- 2 mostVarGeneslist <- most_variable_genes(K_meanTable, number_of_annotations_per_cluster, k_clusters) print(mostVarGeneslist) ``` To set the annotation for a heatmap the following function can be used: ```{r exampleSetAnnotation, echo=TRUE} number_of_annotations_per_cluster <- 2 fontsize_rowAnnotation <- 8 annotation_for_heatmap <- set_annotation(sumBioRepsMatrix, mostVarGeneslist, fontsize_rowAnnotation) print(annotation_for_heatmap) ``` "performing_kMeans" is a function specifically written to perform k-mean clustering outside of the general heatmap function. This allows a more reliable splitting of the heatmap by their assigned clusters. ```{r examplePerformingKmeansClustering, echo=TRUE} k_clusters <- 2 split_heatmap_clusters <- performing_kMeans(sumBioRepsMatrix, k_clusters) ``` ```{r} colorPalette <- "RdBu" color_setting(colorPalette) ``` Finally, a heatmap with clusters can be generated as in this example: ```{r} seed <- 1 title <- "Heatmap DKD Glomerulus" fontsize_columnNames <- 6 fontsize_rowNames <- 4 title_heatmap_legend <- "Expression" WidthNum <- 4.5 HeightNum <- 3 UnitSize <- "cm" color_Palette <- "RdBu" annotation_for_heatmap <- set_annotation(sumBioRepsMatrix, number_of_annotations_per_cluster, fontsize_rowAnnotation) set.seed(seed) print_heatmap(sumBioRepsMatrix, title, split_heatmap_clusters, annotation_for_heatmap, fontsize_columnNames, fontsize_rowNames, title_heatmap_legend, WidthNum, HeightNum, UnitSize, color_Palette) ``` With "function_complexHeatmap_var" a heatmap can be created with automatic annotation of x most variable genes from each cluster. As input a matrix with the scaled counts of the most varible genes from a dataset can be used. The function is able to summarize biological replicates. ```{r exampleFunctionComplexHeatmapVar, echo=TRUE} title <- "Heatmap Cortex Iba1 positive" fontsize_columnNames <- 6 fontsize_rowNames <- 4 title_heatmap_legend <- "Expression" WidthNum <- 4.5 HeightNum <- 3 UnitSize <- "cm" colorPalette <- "RdBu" set.seed(seed) hm <- function_complexHeatmap_var(sumBioRepsMatrix, probes, number_of_annotations_per_cluster, k_clusters, title, fontsize_rowAnnotation, fontsize_columnNames, fontsize_rowNames, title_heatmap_legend, WidthNum, HeightNum, UnitSize, colorPalette) ``` A heatmap can further be generated with annotation of specific genes, as in this example: ```{r exampleHeatmapSpecificAnnotation} probes <- list("disease3_DKD_glomerulus_Geometric_Segment", "disease1B_DKD_glomerulus_Geometric_Segment", "disease2B_DKD_glomerulus_WT") sumBioRepsMatrix <- summarise_bio_replicates(scaled_counts, probes) seed <- 1 k_clusters <- 2 set.seed(seed) K_meanTable <- Kmean_generation(sumBioRepsMatrix, k_clusters) anno_specific_genes <- list("C3", "ALDOB", "ACADS", "APOE") annotation_for_heatmap <- set_annotation(sumBioRepsMatrix, anno_specific_genes, fontsize_rowAnnotation) split_heatmap_clusters <- performing_kMeans(sumBioRepsMatrix, k_clusters) title <- "Heatmap with specific annotation" k_clusters <- 2 fontsize_columnNames <- 6 fontsize_rowNames <- 4 title_heatmap_legend <- "Expression" WidthNum <- 4.5 HeightNum <- 3 UnitSize <- "cm" colorPalette <- "RdBu" set.seed(seed) print_heatmap(sumBioRepsMatrix, title, split_heatmap_clusters, annotation_for_heatmap, fontsize_columnNames, fontsize_rowNames, title_heatmap_legend, WidthNum, HeightNum, UnitSize, colorPalette) ``` **Generate heatmaps with advanced customization:** With the function "adv_Heatmap" users can generate heatmaps using a lot of different settings. The options include: - changeable heatmap title - changeable color palettes in the heatmap - choice between the cluster methods: hierarchical and k-means - choice between distancing methods for hierarchical clustering (e.g. euclidean, correlation) - optional clustering of rows - optional clustering of columns - optional splitting of rows through k-means - optional splitting of columns through k-means - optional annotation of the columns based on an input metadata file - changeable color schemes for the column annotation - changeable side of the column annotation - optional display of row names - optional display of column names - optional annotation of some row names: either automatic through choosing the most variable rows or specific using an input list of the rownames - changeable number of automatic row annotations - changeable font size of the row annotation - changeable font size of the column names - changeable font size of the row names - changeable title of the heatmap legend - changeable width size of the heatmap - changeable height size of the heatmap - changeable unit used for the heatmap sizes The advanced heatmap can then be generated usind adv_Heatmap() and choosing parameter options by specifically setting them or alternatively use the parameters default. The only input necessary to directly generate a heatmap is the input matrix as shown below. ```{r exampleadv_Heatmap_default} # default settings of the parameters ncounts_matrix <- scaled_counts # input matrix column_name <- "Heatmap" # name for heatmap, default = "Heatmap" colorPalette <- NULL # available color palettes from RColorBrewer (), default = NULL cluster_method <- "hierarchical" # cluster methods, either "hierarchical" or "kmeans", default = "hierarchical" distance_method <- "euclidean" # method for creating the distance matrix for hierarchical clustering, default = "euclidean" cluster_rows <- TRUE # optional clustering of rows, default = TRUE cluster_columns <- FALSE # optional clustering of columns, default = FALSE k_row <- NULL # splitting of rows in heatmaps using k-means, would be an integer, default = NULL k_col <- NULL # splitting of columns in heatmaps using k-means, would be an integer, default = NULL sample_metadata <- NULL # dataframe containing the metadata information of the file, default = NULL annotation_colors <- NULL # list containing the column groups and the choosen colors for the column annotation per group, default = NULL annotation_name_side <- "right" # optional change: side of annotation name, default = "right" show_row_names <- FALSE # optional change: show of rownames on = TRUE & off = FALSE, default = FALSE show_column_names <- TRUE # optional change: show of column names on = TRUE & off = FALSE, default = TRUE row_annotation <- FALSE # optional row annotation, default = FALSE row_annotation_method <- "auto" # set if row_annotation = TRUE: options are "auto" & "specific", default = "auto" row_anno_names <- NULL # set if row_annotation = TRUE & row_annotation_method = "specific: input list of specific genes for the row annotation, default = NULL row_anno_number <- 5 # optional change: number of automatic row annotations per cluster, default = 5 fontsize_title <- 15 # optional change: fontsize the heatmap title, default = 15 fontsize_rowAnnotation <- 8 # optional change: fontsize of the optional row annotation, default = 8 fontsize_columnNames <- 6 # optional change: fontsize of the optional column names, default = 6 fontsize_rowNames <- 4 # optional change: fontsize of the optional row names, default = 4 fontsize_cluster_labels <- 8 # optional change: fontsize of the cluster labels, default = 8 fontsize_group_annotation <- 8 # optional change: font size of the group annotation title, default = 8. fontsize_group_annotation_legend <- 10 # optional change: fontsize of optional group annotation legend title, default = 10 fontsize_group_annotation_labels <- 8 # optional change: fontsize of optional group annotation legend labels, default = 8 fontsize_heatmap_legend <- 8 # optional change: fontsize of heatmap legend, default = 10 fontsize_heatmap_legend_labels <- 8 # optional change: fontsize of heatmap legend labels, default = 8 title_heatmapLegend <- "Expression" # changeable title of the legend, default "Expression" WidthNum <- 4.5 # optional change of heatmap width, default = 4.5 HeightNum <- 3 # optional change of heatmap height, default = 3 UnitSize <- "cm" # optional change of heatmap unit for sizes, default = "cm" set.seed(1) # Using the adv_Heatmap() function with only the input matrix file. All other parameters use their default option. hm <- adv_Heatmap(scaled_counts) ``` All the optional parameters can then be changed to the taste and specifications of the user. For example, the color scheme can be changed, the row names can be shown in the heatmap and a group annotation can be pictured using a metadata file as input. This metadata file can either be loaded into R as an excel or csv file, or alternatively be produced within R as shown below. ```{r generatingMetaFile} # set a list of the groups (names should be contained within the sample names) groups <- c("3_DKD_glomerulus_Geometric_S", "1B_DKD_glomerulus_Geometric_S", "2B_DKD_glomerulus_WT") # get the sample names from the used data set sample_names <- c(colnames(scaled_counts)) # Match each sample name to the correct group group_assignment <- sapply(sample_names, function(sample) { matched <- groups[sapply(groups, function(g) grepl(g, sample))] if (length(matched) > 0) matched[1] else NA }) # Ensure same length stopifnot(length(sample_names) == length(group_assignment)) # Create a dataframe with the sample names and their respective group assignments sample_metadata <- data.frame(Group <- group_assignment, row.names <- sample_names) # confirm matrix column names match the metadata rownames all(colnames(scaled_counts) == rownames(sample_metadata)) # set a list with the Groups and choose colors for them group_colors <- list(Group = c("3_DKD_glomerulus_Geometric_S" = "#1b9e77", "1B_DKD_glomerulus_Geometric_S" = "#7570b3", "2B_DKD_glomerulus_WT" = "#e7298a")) names(group_colors$Group) print(sample_metadata) ``` The heatmap parameters can then be changed to generate a more advanced and custom heatmap: ```{r exampleadv_heatmap} # parameters and their options in adv_Heatmap() ncounts_matrix <- scaled_counts # input matrix column_name <- "Heatmap 1" # name for heatmap, default = "Heatmap" colorPalette <- "RdBu" # available color palettes from RColorBrewer (), default = NULL cluster_method <- "hierarchical" # cluster methods, either "hierarchical" or "kmeans", default = "hierarchical" distance_method <- "euclidean" # method for creating the distance matrix for hierarchical clustering, default = "euclidean" cluster_rows <- TRUE # optional clustering of rows, default = TRUE cluster_columns <- FALSE # optional clustering of columns, default = FALSE k_row <- NULL # splitting of rows in heatmaps using k-means, would be an integer, default = NULL k_col <- NULL # splitting of columns in heatmaps using k-means, would be an integer, default = NULL sample_metadata <- sample_metadata # dataframe containing the metadata information of the file, default = NULL annotation_colors <- group_colors # list containing the column groups and the choosen colors for the column annotation per group, default = NULL annotation_name_side <- "right" # optional change: side of annotation name, default = "right" show_row_names <- TRUE # optional change: show of rownames on = TRUE & off = FALSE, default = FALSE show_column_names <- TRUE # optional change: show of column names on = TRUE & off = FALSE, default = TRUE row_annotation <- FALSE # optional row annotation, default = FALSE row_annotation_method <- "auto" # set if row_annotation = TRUE: options are "auto" & "specific", default = "auto" row_anno_names <- NULL # set if row_annotation = TRUE & row_annotation_method = "specific: input list of specific genes for the row annotation, default = NULL row_anno_number <- 5 # optional change: number of automatic row annotations per cluster, default = 5 fontsize_title <- 15 # optional change: fontsize the heatmap title, default = 15 fontsize_rowAnnotation <- 8 # optional change: fontsize of the optional row annotation, default = 8 fontsize_columnNames <- 6 # optional change: fontsize of the optional column names, default = 6 fontsize_rowNames <- 4 # optional change: fontsize of the optional row names, default = 4 fontsize_cluster_labels <- 8 # optional change: fontsize of the cluster labels, default = 8 fontsize_group_annotation <- 9 # optional change: font size of the group annotation title, default = 8. fontsize_group_annotation_legend <- 9 # optional change: fontsize of optional group annotation legend title, default = 10 fontsize_group_annotation_labels <- 8 # optional change: fontsize of optional group annotation legend labels, default = 8 fontsize_heatmap_legend <- 9 # optional change: fontsize of heatmap legend, default = 10 fontsize_heatmap_legend_labels <- 8 # optional change: fontsize of heatmap legend labels, default = 8 title_heatmapLegend <- "Expression" # changeable title of the legend, default "Expression" WidthNum <- 4.5 # optional change of heatmap width, default = 4.5 HeightNum <- 3 # optional change of heatmap height, default = 3 UnitSize <- "cm" # optional change of heatmap unit for sizes, default = "cm" set.seed(1) hm <- adv_Heatmap(ncounts_matrix, column_name = column_name, colorPalette = colorPalette, cluster_method = cluster_method, cluster_rows = cluster_rows, cluster_columns = cluster_columns, k_row = k_row, k_col = k_col, sample_metadata = sample_metadata, annotation_colors = annotation_colors, annotation_name_side = annotation_name_side, show_row_names = show_row_names, show_column_names = show_column_names, row_annotation = row_annotation, row_annotation_method = row_annotation_method, row_anno_names = row_anno_names, row_anno_number = row_anno_number, fontsize_title = fontsize_title, fontsize_rowAnnotation = fontsize_rowAnnotation, fontsize_columnNames = fontsize_columnNames, fontsize_rowNames = fontsize_rowNames, fontsize_cluster_labels = fontsize_cluster_labels, fontsize_group_annotation = fontsize_group_annotation, fontsize_group_annotation_legend = fontsize_group_annotation_legend, fontsize_group_annotation_labels = fontsize_group_annotation_labels, fontsize_heatmap_legend = fontsize_heatmap_legend, fontsize_heatmap_legend_labels = fontsize_heatmap_legend_labels, title_heatmapLegend = title_heatmapLegend, WidthNum = WidthNum, HeightNum = HeightNum, UnitSize = UnitSize) ``` Alternatively all these options can also be conducted using k-means Clustering: ```{r K_meansadv_heatmap} # parameters and their options in adv_heatmap() ncounts_matrix <- scaled_counts # input matrix column_name <- "Heatmap 2" # name for heatmap, default = "Heatmap" colorPalette <- "RdBu" # available color palettes from RColorBrewer (), default = NULL cluster_method <- "kmeans" # cluster methods, either "hierarchical" or "kmeans", default = "hierarchical" distance_method <- "euclidean" # method for creating the distance matrix for hierarchical clustering, default = "euclidean" cluster_rows <- TRUE # optional clustering of rows, default = TRUE cluster_columns <- TRUE # optional clustering of columns, default = FALSE k_row <- 3 # splitting of rows in heatmaps using k-means, would be an integer, default = NULL k_col <- 2 # splitting of columns in heatmaps using k-means, would be an integer, default = NULL sample_metadata <- sample_metadata # dataframe containing the metadata information of the file, default = NULL annotation_colors <- group_colors # list containing the column groups and the choosen colors for the column annotation per group, default = NULL annotation_name_side <- "right" # optional change: side of annotation name, default = "right" show_row_names <- FALSE # optional change: show of rownames on = TRUE & off = FALSE, default = FALSE show_column_names <- FALSE # optional change: show of column names on = TRUE & off = FALSE, default = TRUE row_annotation <- TRUE # optional row annotation, default = FALSE row_annotation_method <- "auto" # set if row_annotation = TRUE: options are "auto" & "specific", default = "auto" row_anno_names <- NULL # set if row_annotation = TRUE & row_annotation_method = "specific: input list of specific genes for the row annotation, default = NULL row_anno_number <- 5 # optional change: number of automatic row annotations per cluster, default = 5 fontsize_title <- 15 # optional change: fontsize the heatmap title, default = 15 fontsize_rowAnnotation <- 8 # optional change: fontsize of the optional row annotation, default = 8 fontsize_columnNames <- 6 # optional change: fontsize of the optional column names, default = 6 fontsize_rowNames <- 4 # optional change: fontsize of the optional row names, default = 4 fontsize_cluster_labels <- 12 # optional change: fontsize of the cluster labels, default = 8 fontsize_group_annotation <- 10 # optional change: font size of the group annotation title, default = 8. fontsize_group_annotation_legend <- 10 # optional change: fontsize of optional group annotation legend title, default = 10 fontsize_group_annotation_labels <- 8 # optional change: fontsize of optional group annotation legend labels, default = 8 fontsize_heatmap_legend <- 10 # optional change: fontsize of heatmap legend, default = 10 fontsize_heatmap_legend_labels <- 8 # optional change: fontsize of heatmap legend labels, default = 8 title_heatmapLegend <- "Expression" # changeable title of the legend, default "Expression" WidthNum <- 5 # optional change of heatmap width, default = 4.5 HeightNum <- 7 # optional change of heatmap height, default = 3 UnitSize <- "cm" # optional change of heatmap unit for sizes, default = "cm" set.seed(1) hm <- adv_Heatmap(ncounts_matrix, column_name = column_name, colorPalette = colorPalette, cluster_method = cluster_method, cluster_rows = cluster_rows, cluster_columns = cluster_columns, k_row = k_row, k_col = k_col, sample_metadata = sample_metadata, annotation_colors = annotation_colors, annotation_name_side = annotation_name_side, show_row_names = show_row_names, show_column_names = show_column_names, row_annotation = row_annotation, row_annotation_method = row_annotation_method, row_anno_names = row_anno_names, row_anno_number = row_anno_number, fontsize_title = fontsize_title, fontsize_rowAnnotation = fontsize_rowAnnotation, fontsize_columnNames = fontsize_columnNames, fontsize_rowNames = fontsize_rowNames, fontsize_cluster_labels = fontsize_cluster_labels, fontsize_group_annotation = fontsize_group_annotation, fontsize_group_annotation_legend = fontsize_group_annotation_legend, fontsize_group_annotation_labels = fontsize_group_annotation_labels, fontsize_heatmap_legend = fontsize_heatmap_legend, fontsize_heatmap_legend_labels = fontsize_heatmap_legend_labels, title_heatmapLegend = title_heatmapLegend, WidthNum = WidthNum, HeightNum = HeightNum, UnitSize = UnitSize) ``` ### Differential Gene Expression Analysis Using Nanostring Data **Functions to analyze Nanostring GeoMx DSP data:** ```{{r loadingPackagesNanostringData, eval=requireNamespace("GeoMxWorkflows", quietly = TRUE), warning=FALSE, include=FALSE, results='hide'} library(NanoStringNCTools) library(GeomxTools) library("GeoMxWorkflows") library(EnvStats) library(ggiraph) if (packageVersion("GeomxTools") < "2.1" && packageVersion("GeoMxWorkflows") >= "1.0.1") { stop("GeomxTools and Workflow versions do not match. Please use the same version. Rhis workflow is meant to be used with most current version of packages.If you are using an older version of Bioconductor please reinstall GeoMxWorkflows and use vignette(GeoMxWorkflows) instead") } if (packageVersion("GeomxTools") > "2.1" && packageVersion("GeoMxWorkflows") <= "1.0.1") { stop("GeomxTools and Workflow versions do not match. Please use the same version, see install instructions above.") } ``` ```{r readInNanostringRawFiles, eval=FALSE, echo=TRUE, results='hide', warning=FALSE} datadir <- system.file("extdata/WTA_NGS_Example", package = "GeoMxWorkflows") DCCFiles <- dir(file.path(datadir, "dccs"), pattern = ".dcc$", full.names = TRUE, recursive = TRUE ) PKCFiles <- unzip(zipfile = dir(file.path(datadir, "pkcs"), pattern = ".zip$", full.names = TRUE, recursive = TRUE )) # has to be in file directory # SampleAnnotationFile <- file.path(datadir, "kidney_AOI_Annotations_all_vignette.xlsx") SampleAnnotationFile <- dir(file.path(datadir, "annotation"), pattern = ".xlsx$", full.names = TRUE, recursive = TRUE ) ``` The raw Nanostring GeoMx DSP files are read in and then united into an instance of class "NanostringGeoMxSet". ```{r ReadData, eval=FALSE, echo=TRUE, results='hide', warning=FALSE} rawDataObject <- suppressWarnings(readNanoStringGeoMxSet( dccFiles = subset_dccs, pkcFiles = PKCFiles, phenoDataFile = SampleAnnotationFile, phenoDataSheet = "Template", phenoDataDccColName = "Sample_ID", protocolDataColNames = c("aoi", "roi"), configFile = NULL, analyte = "RNA", phenoDataColPrefix = "", experimentDataColNames = NULL )) ``` Anon Analyzing GeoMx-NGS RNA Expression Data with GeomxTools. Available from: https://bioconductor.statistik.tu-dortmund.de/packages/3.18/workflows/vignettes/GeoMxWorkflows/inst/doc/GeomxTools_RNA-NGS_Analysis.html (September 21, 2023b). Anon Developer Introduction to the NanoStringGeoMxSet. Available from: https://www.bioconductor.org/packages/devel/bioc/vignettes/GeomxTools/inst/doc/Developer_Introduction_to_the_NanoStringGeoMxSet.html (September 21, 2023c). To preprocess the data, the logarithm of the count matrix is computed and added to the demoData object. Furthermore, the data is summarized by splitting it by a chosen column and calculating the mean. ```{r examplePreProcess, eval=FALSE, echo=TRUE, results='hide'} PrePro_rawDataObject <- add_demoElem(rawDataObject) class(PrePro_rawDataObject) # loop over the features(1) or samples(2) of the assayData element and get the mean assayDataApply(PrePro_rawDataObject, MARGIN = 1, FUN = mean, elt = "demoElem")[seq_len(5)] # split the data by group column with feature, pheno or protrocol data then get the mean VGroup <- "aoi" elt <- "demoElem" PrePro_rawDataObject <- split_data_by_column(PrePro_rawDataObject, VGroup, elt) # Quality Control QCPassed <- aExprsDataQC(PrePro_rawDataObject, "QCFlags") df_Exp <- genRawReadCountTable(PrePro_rawDataObject) ``` ```{r showExpressionMatrix, echo=FALSE, warning=FALSE} file_path <- system.file("extdata", "RawDataExamplePackageNanostring.csv", package = "DgeaHeatmap") df_Exp <- read.csv(file_path) copy_df_Expr <- build_matrix(df_Exp, 1) knitr::kable(head(df_Exp[, seq_len(5)]), "pipe") path_metadata <- system.file("extdata", "MetaDataPackageNanostring.rds", package = "DgeaHeatmap") coldata_2 <- readRDS(path_metadata) coldata_2 <- build_matrix(coldata_2, 1) ``` Anon Analyzing GeoMx-NGS RNA Expression Data with GeomxTools. Available from: https://bioconductor.statistik.tu-dortmund.de/packages/3.18/workflows/vignettes/GeoMxWorkflows/inst/doc/GeomxTools_RNA-NGS_Analysis.html (September 21, 2023b). Anon Developer Introduction to the NanoStringGeoMxSet. Available from: https://www.bioconductor.org/packages/devel/bioc/vignettes/GeomxTools/inst/doc/Developer_Introduction_to_the_NanoStringGeoMxSet.html (September 21, 2023c). Next, an usable matrix of column data is extracted from the sample data: ```{r ColdataDEG, eval=FALSE, echo=TRUE, results='hide', warning=FALSE} annotation_matrix <- sData(rawDataObject) annotation_matrix$Samplename <- paste(annotation_matrix$`slide name`, annotation_matrix$class, annotation_matrix$region, annotation_matrix$segment, annotation_matrix$roi, sep = "_") class(annotation_matrix) library(tibble) annotation_matrix <- tibble::rownames_to_column(annotation_matrix, "Samplenummern") rownames(annotation_matrix) annotation_matrix1 <- data.frame(annotation_matrix$Samplenummern, annotation_matrix$Samplename) coldata_2 <- data.frame(annotation_matrix$Samplename, annotation_matrix$`slide name`, annotation_matrix$class, annotation_matrix$region, annotation_matrix$segment, annotation_matrix$roi, annotation_matrix$Samplenummern) # Renaming the columns names(coldata_2)[names(coldata_2) == "annotation_matrix.Samplename"] <- "Samplename" names(coldata_2)[names(coldata_2) == "annotation_matrix..slide.name"] <- "slide_name" names(coldata_2)[names(coldata_2) == "annotation_matrix.class"] <- "class" names(coldata_2)[names(coldata_2) == "annotation_matrix.segment"] <- "segment" names(coldata_2)[names(coldata_2) == "annotation_matrix.sample"] <- "sample" names(coldata_2)[names(coldata_2) == "annotation_matrix.region"] <- "region" names(coldata_2)[names(coldata_2) == "annotation_matrix.Samplenummern"] <- "Samplenummern" coldata_2$Samplename <- gsub(" ", "_", coldata_2$Samplename) list_rownames <- as.list(coldata_2[, 1]) # Combines sample describing columns into one coldata_2$comp <- paste(coldata_2$segment, coldata_2$region, coldata_2$class, coldata_2$slide_name, sep = "_") annotationMatrix <- data.frame(coldata_2$Samplenummern, coldata_2$Samplename) coldata_2 <- coldata_2[, -1] coldata_2 <- as.matrix(coldata_2) list_rownames <- gsub("\\.", "_", list_rownames) uniq_rownames <- make.names(list_rownames, unique = TRUE) # makes all names in the list_rownames unique to later use them as rownames print(uniq_rownames) uniq_rownames <- gsub("\\.", "_", uniq_rownames) rownames(coldata_2) <- uniq_rownames # Creates a new dataframe wich only includes the necessary information coldata_df <- as.data.frame(coldata_2) class(coldata_df) # exchange of "." with "_" in the rownames of the column data rownames(coldata_2) <- gsub("\\.", "_", rownames(coldata_2)) ``` The columns in the counts table can now be replaced with the corresponding sample names from the table containing the column data. To do this, the following code can be used: ```{r exampleExchangeColumnNamesWMatchingSampleNames, eval=FALSE, echo=TRUE, results='hide'} list_columnNames <- as.list(colnames(df_Exp)) # list of column names in counts table list_rownames <- as.list(annotationMatrix[, 2]) annotationMatrix <- as.matrix(annotationMatrix) rownames(annotationMatrix) <- list_rownames x <- match(rownames(annotationMatrix), colnames(df_Exp)) print(x) copy_df_Expr <- df_Exp # creates copy of counts table annotationMatrix <- as.data.frame(annotationMatrix) for (i in list_columnNames) { matchingRow <- which(annotationMatrix$coldata_2.Samplenummern == i) matchingColumn <- match(i, names(copy_df_Expr)) print(matchingColumn) print(matchingRow) colnames(copy_df_Expr)[matchingColumn] <- rownames(annotationMatrix)[matchingRow] print(rownames(annotationMatrix)[matchingRow]) } # exchange of "." with "_" in the column names of copy_df_Expr colnames(copy_df_Expr) <- gsub("\\.", "_", colnames(copy_df_Expr)) ``` ```{r exampleColumnNameExchange, echo=FALSE} knitr::kable(head(copy_df_Expr[, seq_len(5)]), "pipe") ``` **Differential gene expression analysis** The process of performing a differential gene expression analysis starts by checking whether all the columns of the raw data are present in the metadata and if the column names appear in the same order in both the raw data and metadata. ```{r EqualOrder, echo=TRUE, warning=FALSE} all(colnames(copy_df_Expr) %in% rownames(coldata_2)) # check if all column names of data are in rownames of metadata #check if TRUE class(coldata_2) all(colnames(copy_df_Expr) == rownames(coldata_2)) # check if the order of the data column names == order of metadata rownames ``` Love M, Ahlmann-Eltze C, Forbes K, Anders S, Huber W, FP7 RE, NHGRI N, CZI (2023) DESeq2: Differential gene expression analysis based on the negative binomial distribution. https://doi.org/10.18129/B9.bioc.DESeq2 The differential gene expression analysis can be conducted based on limma voom, DESeq2, and edgeR. **DEA Limma** Utilizing limma voom, a matrix is created in which each coefficient represents the average of the samples from one group. Afterwards, the contrasts are set to direct testing between chosen groups in the analysis. To scale the raw library sizes, the normalization factors are calculated. The count data is then transformed to log2-counts per million, and the mean-variance relationship is estimated. By fitting a linear model to each gene, both fold changes and standard errors are estimated. The fitted model objects are then reorientated from the coefficients of the design matrix to any set of contrasts of the original coefficients. The standard errors are smoothed by empirical Bayes. Finally, the results of the differential gene expression analysis are classified as up, down, or not significantly differently expressed. ```{r PerformingDEGLimmaVoomfct,echo=TRUE, include=TRUE, error=TRUE, warning=FALSE, eval=FALSE} coldata2 <- as.data.frame((coldata_2)) comparisons <- list( "Geometric_Segment_glomerulus_DKD_disease3_vs_PanCK_tubule_DKD_disease4" = c("Geometric_Segment_glomerulus_DKD_disease3", "PanCK_tubule_DKD_disease4"), "Geometric_Segment_glomerulus_DKD_disease3_vs_neg_tubule_DKD_disease4" = c("Geometric_Segment_glomerulus_DKD_disease3", "neg_tubule_DKD_disease4"), "PanCK_tubule_DKD_disease4_vs_neg_tubule_DKD_disease4" = c("PanCK_tubule_DKD_disease4", "neg_tubule_DKD_disease4") ) groupingColumns <- c("segment", "region", "class", "slide_name") DGEA_results <- DGEALimma(copy_df_Expr, coldata2, grouping_columns = groupingColumns, comparisons = comparisons) results_all_DEG <- limma::decideTests(DGEA_results$fit) summary(results_all_DEG) result1 <- topTable(fit2, coef = "Geometric_Segment_glomerulus_DKD_disease3 - PanCK_tubule_DKD_disease4", number = Inf, adjust.method = "fdr") %>% as.data.frame() # differentially expressed genes are obtained by topTreat() function # filter results to get significantly differentially expressed genes topUp <- result1[which(result1$logFC > 0), ][seq_len(100), ] # up reg top 100 topDown <- result1[which(result1$logFC < 0), ][seq_len(100), ] # down reg top 100 class(topUp) NamesUpReg <- row.names(topUp) NamesDownReg <- row.names(topDown) print(NamesDownReg) result1 <- result1 %>% dplyr::mutate(isSignificant = case_when( adj.P.Val < 0.05 & abs(logFC) > 1 ~ TRUE, TRUE ~ FALSE # If condictions in the line above are not met, gene is not DE. )) sigDEresults <- result1 %>% dplyr::filter(isSignificant == TRUE) ``` inesdesantiago (2020) 10 Tips & Tricks for complex model.matrix designs in DGE analysis. Sequencing QC and data analysis blog. Available from: Bibliography 96 https://seqqc.wordpress.com/2020/11/28/10-tips-tricks-for-complex-model-matrix-designs-in-dge-analysis/ (September 21, 2023). Anon Developer Introduction to the NanoStringGeoMxSet. Available from: https://www.bioconductor.org/packages/devel/bioc/vignettes/GeomxTools/inst/doc/Developer_Introduction_to_the_NanoStringGeoMxSet.html (September 21, 2023c). **Visualisation with a Venn Diagramm:** The results from the differential gene expression analysis are visualized using a Venn diagramm to represent the up- and down- regulated genes between groups. This is set up in the following way: ```{r VennDiagramm, echo=TRUE, error=TRUE, warning=FALSE, eval=FALSE} library(VennDiagram) library(limma) venn.plot <- vennDiagram(results_all_DEG, imagetype = "tiff", include = c("up", "down"), mar = rep(1, 4), cex = c(1, 1, 0.7), lwd = 1, counts.col = c("red", "blue"), circle.col = c("red", "blue", "green3") ) ``` ```{r ImageVennLimmaupdownreg, echo=FALSE, out.width="100%"} knitr::include_graphics( system.file("extdata", "VennLimmaupdownreg.png", package = "DgeaHeatmap") ) ``` Chen, H., Boutros, P.C. VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R. BMC Bioinformatics 12, 35 (2011). https://doi.org/10.1186/1471-2105-12-35 **DEA DESeq2** The raw counts are first structured into a DESeqDataSEt and comparisons are used to test specific differences between groups. Next, the size factors are estimated for normalization and dispersion values are estimated for each gene. A GLM (generalized linear model) is fit to the data through negative niomial distribution. A Wald test is performed for statistical testing and p-values are adjusted for multiple testing through Benjamin-Hochberg procedure. The data is then classified into up and down reguated and not significantly expressed genes. ```{r PerformingDEGDESeq2fct,echo=TRUE, include=TRUE, error=TRUE, warning=FALSE, eval=FALSE} comparisons <- list( "disease1b_vs_disease2B" = c("Geometric_Segment_glomerulus_DKD_disease1B", "Geometric_Segment_glomerulus_DKD_disease2B"), "disease1b_vs_disease3" = c("Geometric_Segment_glomerulus_DKD_disease1B", "Geometric_Segment_glomerulus_DKD_disease3"), "disease2B_vs_disease3" = c("Geometric_Segment_glomerulus_DKD_disease2B", "Geometric_Segment_glomerulus_DKD_disease3") ) int_counts <- round(copy_df_Expr) groupingColumns <- c("segment", "region", "class", "slide_name") results_DESeq2 <- DGEADESeq2(int_counts, coldata2, groupingColumns, comparisons) results_list_d <- results_DESeq2$results ``` Venn Diagram for up regulated genes extracted from the DESeq2 analysis: ```{r DESeq2_Venn_upreg, warning=FALSE, eval=FALSE} library(ggvenn) up_genes <- extractDEGenes(results_list_d, contrasts, only_up = TRUE) ggvenn(up_genes, fill_color = c("red", "green", "blue")) ``` ```{r ImageVennDESeq2upreg, echo=FALSE, out.width="100%"} knitr::include_graphics( system.file("extdata", "VennDESeq2upreg.png", package = "DgeaHeatmap") ) ``` Venn Diagram for down regulated genes extracted from the DESeq2 analysis: ```{r DESeq2_Venn_downreg, eval=FALSE} down_genes <- extractDEGenes(results_list_d, contrasts, only_down = TRUE) ggvenn(down_genes, fill_color = c("red", "green", "blue")) ``` ```{r ImageVennDESeq2downreg, echo=FALSE, out.width="100%"} knitr::include_graphics( system.file("extdata", "VennDESeq2downreg.png", package = "DgeaHeatmap") ) ``` Venn Diagram for significant genes extracted from the DESeq2 analysis: ```{r DESeq2_Venn_sigreg, eval=FALSE} sig_genes <- extractDEGenes(results_list_d, contrasts, only_sig = TRUE) ggvenn(sig_genes, fill_color = c("red", "green", "blue")) ``` ```{r ImageVennDESeq2sigreg, echo=FALSE, out.width="100%"} knitr::include_graphics( system.file("extdata", "VennDESeq2sigreg.png", package = "DgeaHeatmap") ) ``` Venn Diagram for up and down regulated genes extracted from the DESeq2 analysis: ```{r DESeq2_Venn_updownreg, eval=FALSE} state_mat <- extractDEGenes(results_list_d, contrasts, up_down = TRUE) head(state_mat) venn.plot <- vennDiagram(state_mat, imagetype = "tiff", include = c("up", "down"), mar = rep(1, 4), cex = c(1.5, 1, 0.7), lwd = 1, counts.col = c("red", "blue"), circle.col = c("red", "blue", "green3") ) ``` ```{r ImageVennDESeq2updownreg, echo=FALSE, out.width="100%"} knitr::include_graphics( system.file("extdata", "VennDESeq2updownreg.png", package = "DgeaHeatmap") ) ``` **DEA edgeR** The raw counts are first used to construct a DGEList object and grouping factors are identified. Then genes with low expression are filtered before the normalization factors are calculated using trimmed mean of m-values (TMM). A linear model is fitted to each gene and fold changes, as well as standard errors are estimated. The model is then reorientated using contrasts (set as comparison) of the original coefficients. The results of the differential gene expression analysis are classified as up, down, or not significantly differently expressed. ```{r DEA_edgeR, eval=FALSE} comparisons <- list( "GeometricSegment_glomerulus_Dkd_disease1b_vs_disease2B" = c("Geometric_Segment_glomerulus_DKD_disease1B ", "Geometric_Segment_glomerulus_DKD_disease2B "), "GeometricSegmentglomerulus_DKD_disease1B_vs_disease3" = c("Geometric_Segment_glomerulus_DKD_disease1B ", "Geometric_Segment_glomerulus_DKD_disease3"), "GeometricSegmentglomerulus_DKD_disease2B_vs_disease3" = c("Geometric_Segment_glomerulus_DKD_disease2B ", "Geometric_Segment_glomerulus_DKD_disease3") ) groupingColumns <- c("segment", "region", "class", "slide_name") results_edgeR <- DGEAedgeR(int_counts, coldata2, groupingColumns, comparisons, prefix = "DEA") edgeR_summary <- summarize_edgeR_DEA(results_edgeR) ``` Venn Diagram for up and down regulated genes extracted from the edgeR analysis: ```{r DEA_Venn_edgeR, warning=FALSE, eval=FALSE} # Suppose `edgeR_summary` is the output of summarize_edgeR_DEA() state_matrix <- edgeR_summary$classified_results %>% lapply(function(df) df$decision) %>% # extract the –1/0/+1 vector for each contrast do.call(cbind, .) # bind into a matrix colnames(state_matrix) <- names(edgeR_summary$classified_results) venn.plot <- vennDiagram(state_matrix, imagetype = "tiff", include = c("up", "down"), mar = rep(1, 4), cex = c(1, 0.8, 0.7), lwd = 1, counts.col = c("red", "blue"), circle.col = c("red", "blue", "green3") ) ``` ```{r ImageVennedgeRupdownreg, echo=FALSE, out.width="100%"} knitr::include_graphics( system.file("extdata", "VennedgeRupdownreg.png", package = "DgeaHeatmap") ) ``` Normalized Counts can for example be extracted as in this example: ```{r extractionNormalizedCounts, echo=FALSE, warning=FALSE, results='hide', eval=FALSE} library(data.table) y <- DGEA_results$normFactors class(y) head(y) CPM <- edgeR::cpm(y) head(CPM) class(CPM) df_normCounts <- data.frame("genes" = row.names(CPM), CPM) class(df_normCounts) write.csv(df_normCounts, "/Users/llancell/Documents/gitLab_SpittauLab/Dgea_Heatmap_Package/tests/testthat/NormCounts_ExampleDataNanostring.csv", row.names = FALSE) ``` ```{r sessionInfo, echo=FALSE} sessionInfo() ```