DgeaHeatmap Vignette

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 with:

# 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:

library(DgeaHeatmap)
x <- 1
matrixCounts <- build_matrix(input_data, x)
No_Template_Control_NA_NA_NA_NA disease3_DKD_glomerulus_Geometric_Segment_7 disease3_DKD_glomerulus_Geometric_Segment_8 disease4_DKD_tubule_PanCK_1
A2M 0 1183.26647 983.54729 87.92538
NAT2 0 36.59587 67.57195 74.04243
ACADM 0 75.63146 56.30996 92.55303
ACADS 0 65.87257 63.81795 106.43599
ACAT1 0 70.75202 90.09594 236.01024
ACVRL1 0 126.86568 82.58794 64.78712 
parameter1 <- "DKD"
parameter2 <- "glomerulus"
factors_for_individual_matrix <- list(parameter1, parameter2)
indiMatrix <- individual_matrix(factors_for_individual_matrix, matrixCounts)
disease3_DKD_glomerulus_Geometric_Segment_7 disease3_DKD_glomerulus_Geometric_Segment_8 disease1B_DKD_glomerulus_Geometric_Segment_14 disease2B_DKD_glomerulus_WT1_1
A2M 1183.26647 983.54729 684.32136 102.5730
NAT2 36.59587 67.57195 68.43214 0.0000
ACADM 75.63146 56.30996 34.21607 0.0000
ACADS 65.87257 63.81795 91.24285 0.0000
ACAT1 70.75202 90.09594 68.43214 0.0000
ACVRL1 126.86568 82.58794 125.45892 153.8595

Filtering the matrix for only an x amount of most variable expressed genes:

top_number_of_genes <- 20
varGenesMatrix <- filtering_for_top_exprGenes(indiMatrix, top_number_of_genes)
print(nrow(varGenesMatrix))
#> [1] 20

Using Z-score scaling to scale the values of the matrix:

scaled_counts <- scale_counts(varGenesMatrix)

How to visualize the data distribution of the scaled counts:

show_data_distribution(scaled_counts)

Generating an elbow plot to choose k for k-mean generation:

seed <- 1 # setting a seed for a reproducible outcome
set.seed(seed) # set seed for reproducibility
elbow_plot(scaled_counts)
#> Warning: empty cluster: try a better set of initial centers
#> Warning: empty cluster: try a better set of initial centers
#> Warning: empty cluster: try a better set of initial centers
#> Warning: empty cluster: try a better set of initial centers
#> Warning: empty cluster: try a better set of initial centers
#> Warning: empty cluster: try a better set of initial centers
#> Warning: empty cluster: try a better set of initial centers
#> Warning: empty cluster: try a better set of initial centers
#> Warning: empty cluster: try a better set of initial centers
#> Warning: empty cluster: try a better set of initial centers
#> Warning: empty cluster: try a better set of initial centers
#> Warning: empty cluster: try a better set of initial centers

This can also be conducted with the transposed matrix to get the number of clusters the columns can be divided into:

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.

probes <- list("disease3_DKD_glomerulus_Geometric_Segment", "disease1B_DKD_glomerulus_Geometric_Segment", "disease2B_DKD_glomerulus_WT")
sumBioRepsMatrix <- summarise_bio_replicates(scaled_counts, probes)
disease3_DKD_glomerulus_Geometric_Segment disease1B_DKD_glomerulus_Geometric_Segment disease2B_DKD_glomerulus_WT
A2M 0.7326170 -0.1149016 -1.3503324
C3 -0.2893549 1.4431712 -0.8644613
ALDOB -0.8564917 0.6752631 1.0377204
SERPING1 0.1467075 1.0335593 -1.3269742
APOE 0.8517822 -0.6496811 -1.0538832
ASS1 0.2686383 0.8501191 -1.3873957

Generating K-means for clustering in the heatmap:

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)
#>          disease3_DKD_glomerulus_Geometric_Segment
#> ALDOB                                  -0.85649173
#> AGL                                    -0.48956476
#> A2M                                     0.73261699
#> C3                                     -0.28935494
#> SERPING1                                0.14670747
#> APOE                                    0.85178217
#> ASS1                                    0.26863831
#> ASL                                     0.55092494
#> ACAT1                                   0.58597573
#> ACADS                                   0.24723832
#> SGCA                                   -0.27603565
#> ATP7B                                   0.26014710
#> APOC3                                   0.30221751
#> ADRB2                                   0.83176413
#> ARSB                                    0.71899892
#> BTD                                     0.46960879
#> BTK                                     0.08823472
#> ACADM                                   0.75282225
#> NAT2                                    0.27612361
#> APC                                     0.55773089
#>          disease1B_DKD_glomerulus_Geometric_Segment disease2B_DKD_glomerulus_WT
#> ALDOB                                    0.67526308                   1.0377204
#> AGL                                     -0.51545613                   1.4945857
#> A2M                                     -0.11490160                  -1.3503324
#> C3                                       1.44317122                  -0.8644613
#> SERPING1                                 1.03355929                  -1.3269742
#> APOE                                    -0.64968111                  -1.0538832
#> ASS1                                     0.85011909                  -1.3873957
#> ASL                                      0.36697061                  -1.4688205
#> ACAT1                                    0.28183183                  -1.4537833
#> ACADS                                    0.92623864                  -1.4207153
#> SGCA                                     1.42822535                  -0.8761540
#> ATP7B                                    0.87700526                  -1.3972995
#> APOC3                                    0.83738667                  -1.4418217
#> ADRB2                                   -0.50003065                  -1.1634976
#> ARSB                                    -0.04670777                  -1.3912901
#> BTD                                      0.55755977                  -1.4967774
#> BTK                                      1.12997367                  -1.3064431
#> ACADM                                   -0.22565881                  -1.2799857
#> NAT2                                     0.78140353                  -1.3336507
#> APC                                      0.35054379                  -1.4660056
#>           
#> ALDOB    1
#> AGL      1
#> A2M      2
#> C3       2
#> SERPING1 2
#> APOE     2
#> ASS1     2
#> ASL      2
#> ACAT1    2
#> ACADS    2
#> SGCA     2
#> ATP7B    2
#> APOC3    2
#> ADRB2    2
#> ARSB     2
#> BTD      2
#> BTK      2
#> ACADM    2
#> NAT2     2
#> APC      2
disease3_DKD_glomerulus_Geometric_Segment disease1B_DKD_glomerulus_Geometric_Segment disease2B_DKD_glomerulus_WT
ALDOB -0.8564917 0.6752631 1.0377204 1
AGL -0.4895648 -0.5154561 1.4945857 1
A2M 0.7326170 -0.1149016 -1.3503324 2
C3 -0.2893549 1.4431712 -0.8644613 2
SERPING1 0.1467075 1.0335593 -1.3269742 2
APOE 0.8517822 -0.6496811 -1.0538832 2

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.

number_of_annotations_per_cluster <- 2
k_clusters <- 2
mostVarGeneslist <- most_variable_genes(K_meanTable, number_of_annotations_per_cluster, k_clusters)
print(mostVarGeneslist)
#> [[1]]
#> [1] "AGL"
#> 
#> [[2]]
#> [1] "ALDOB"
#> 
#> [[3]]
#> [1] "BTK"
#> 
#> [[4]]
#> [1] "ACADS"

To set the annotation for a heatmap the following function can be used:

number_of_annotations_per_cluster <- 2
fontsize_rowAnnotation <- 8
annotation_for_heatmap <- set_annotation(sumBioRepsMatrix, mostVarGeneslist, fontsize_rowAnnotation)
print(annotation_for_heatmap)
#> An AnnotationFunction object
#>   function: anno_mark()
#>   position: row 
#>   items: -1 
#>   width: 14.525mm 
#>   height: 1npc 
#>   imported variable: .pos, labels_gp, extend, side, at2index, labels_rot, at, link_height, link_width, link_gp, at2labels, .scale, padding 
#>   subsettable variable: at 
#>   this object is subsettable

“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.

k_clusters <- 2
split_heatmap_clusters <- performing_kMeans(sumBioRepsMatrix, k_clusters)
colorPalette <- "RdBu"
color_setting(colorPalette)
#>  [1] "#053061" "#0A3B70" "#10467F" "#16518E" "#1B5C9E" "#2166AC" "#2870B1"
#>  [8] "#2F79B5" "#3682BA" "#3D8BBF" "#4695C4" "#569FC9" "#66A9CF" "#76B3D4"
#> [15] "#86BDDA" "#95C6DF" "#A2CDE2" "#AFD4E6" "#BCDAEA" "#C9E1ED" "#D4E6F0"
#> [22] "#DBEAF2" "#E3EDF3" "#EBF1F4" "#F3F5F6" "#F7F4F2" "#F8EEE8" "#FAE8DE"
#> [29] "#FBE3D4" "#FCDDCA" "#FBD4BE" "#FAC9B0" "#F8BEA2" "#F6B394" "#F4A886"
#> [36] "#EF9B7A" "#E98D6F" "#E37E64" "#DD7059" "#D7624F" "#D05447" "#C84540"
#> [43] "#C13639" "#BA2832" "#B2192B" "#A41328" "#940E26" "#850923" "#760421"
#> [50] "#67001F"

Finally, a heatmap with clusters can be generated as in this example:

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.

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:

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.

# 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.

# 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))
#> [1] TRUE

# 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)
#> [1] "3_DKD_glomerulus_Geometric_S"  "1B_DKD_glomerulus_Geometric_S"
#> [3] "2B_DKD_glomerulus_WT"
print(sample_metadata)
#>                                                   Group....group_assignment
#> disease3_DKD_glomerulus_Geometric_Segment_7    3_DKD_glomerulus_Geometric_S
#> disease3_DKD_glomerulus_Geometric_Segment_8    3_DKD_glomerulus_Geometric_S
#> disease1B_DKD_glomerulus_Geometric_Segment_14 1B_DKD_glomerulus_Geometric_S
#> disease2B_DKD_glomerulus_WT1_1                         2B_DKD_glomerulus_WT
#>                                                                   row.names....sample_names
#> disease3_DKD_glomerulus_Geometric_Segment_7     disease3_DKD_glomerulus_Geometric_Segment_7
#> disease3_DKD_glomerulus_Geometric_Segment_8     disease3_DKD_glomerulus_Geometric_Segment_8
#> disease1B_DKD_glomerulus_Geometric_Segment_14 disease1B_DKD_glomerulus_Geometric_Segment_14
#> disease2B_DKD_glomerulus_WT1_1                               disease2B_DKD_glomerulus_WT1_1

The heatmap parameters can then be changed to generate a more advanced and custom 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:

# 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
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”.

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.

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)
X disease3_DKD_glomerulus_Geometric_Segment_7 disease3_DKD_glomerulus_Geometric_Segment_8 disease3_DKD_glomerulus_Geometric_Segment_9 disease4_DKD_tubule_PanCK_1
A2M 485 262 225 19
NAT2 15 18 5 16
ACADM 31 15 15 20
ACADS 27 17 24 23
ACAT1 29 24 21 51
ACVRL1 52 22 24 14

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:

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:

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))
disease3_DKD_glomerulus_Geometric_Segment_7 disease3_DKD_glomerulus_Geometric_Segment_8 disease3_DKD_glomerulus_Geometric_Segment_9 disease4_DKD_tubule_PanCK_1 disease4_DKD_tubule_neg_1
A2M 485 262 225 19 367
NAT2 15 18 5 16 59
ACADM 31 15 15 20 99
ACADS 27 17 24 23 127
ACAT1 29 24 21 51 272
ACVRL1 52 22 24 14 65

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.

all(colnames(copy_df_Expr) %in% rownames(coldata_2)) # check if all column names of data are in rownames of metadata #check if TRUE
#> [1] TRUE
class(coldata_2)
#> [1] "matrix" "array"

all(colnames(copy_df_Expr) == rownames(coldata_2)) # check if the order of the data column names == order of metadata rownames
#> [1] TRUE

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.

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:

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")
)

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.

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:

library(ggvenn)
up_genes <- extractDEGenes(results_list_d, contrasts, only_up = TRUE)
ggvenn(up_genes, fill_color = c("red", "green", "blue"))

Venn Diagram for down regulated genes extracted from the DESeq2 analysis:

down_genes <- extractDEGenes(results_list_d, contrasts, only_down = TRUE)
ggvenn(down_genes, fill_color = c("red", "green", "blue"))

Venn Diagram for significant genes extracted from the DESeq2 analysis:

sig_genes <- extractDEGenes(results_list_d, contrasts, only_sig = TRUE)
ggvenn(sig_genes, fill_color = c("red", "green", "blue"))

Venn Diagram for up and down regulated genes extracted from the DESeq2 analysis:

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")
)

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.

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:

# 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")
)

Normalized Counts can for example be extracted as in this example:

#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] DgeaHeatmap_0.99.3 BiocStyle_2.41.0  
#> 
#> loaded via a namespace (and not attached):
#>  [1] ComplexHeatmap_2.29.0 jsonlite_2.0.0        rjson_0.2.23         
#>  [4] compiler_4.6.1        BiocManager_1.30.27   crayon_1.5.3         
#>  [7] parallel_4.6.1        cluster_2.1.8.3       jquerylib_0.1.4      
#> [10] IRanges_2.47.5        png_0.1-9             yaml_2.3.12          
#> [13] fastmap_1.2.0         R6_2.6.1              generics_0.1.4       
#> [16] shape_1.4.6.1         knitr_1.52            BiocGenerics_0.59.12 
#> [19] iterators_1.0.14      GetoptLong_1.1.1      circlize_0.4.18      
#> [22] maketools_1.3.2       RColorBrewer_1.1-3    bslib_0.12.0         
#> [25] rlang_1.3.0           cachem_1.1.0          xfun_0.61            
#> [28] GlobalOptions_0.1.4   sass_0.4.10           sys_3.4.3            
#> [31] otel_0.2.0            doParallel_1.0.17     cli_3.6.6            
#> [34] magrittr_2.0.5        digest_0.6.39         foreach_1.5.2        
#> [37] grid_4.6.1            clue_0.3-68           lifecycle_1.0.5      
#> [40] S4Vectors_0.51.10     vctrs_0.7.3           evaluate_1.0.5       
#> [43] codetools_0.2-20      buildtools_1.0.0      stats4_4.6.1         
#> [46] colorspace_2.1-4      rmarkdown_2.32        purrr_1.2.2          
#> [49] matrixStats_1.5.0     tools_4.6.1           htmltools_0.5.9