Heatmap

[aston9090]

Hello,

Thanks for the package! I've encountered some issues while trying to annotate the KO pathways and plot a heatmap. Below are the details:

1. Annotation of KO Pathways:

> daa_results_df_ALDE <- pathway_daa(abundance = kegg_abundance, metadata = metadata1, group = "Group", daa_method = "ALDEx2", select = NULL, p.adjust = "BH", reference = NULL)
> daa_results_df_ALDE_sig <- daa_results_df_ALDE %>% filter(p_adjust < 0.05)
> daa_results_df_ALDE_sig_sub <- daa_results_df_ALDE_sig[daa_results_df_ALDE_sig$method == "ALDEx2_Kruskal-Wallace test", ]
> daa_annotated_results_df_ALDE_sig_sub <- pathway_annotation(pathway = "KO", daa_results_df = daa_results_df_ALDE_sig_sub, ko_to_kegg = TRUE)

Error Received:

The number of statistically significant pathways exceeds the database's query limit. Please consider breaking down the analysis into smaller queries or selecting a subset of pathways for further investigation.

nrow(daa_annotated_results_df_ALDE_sig_sub)
204

Solution Attempt:
I divided daa_annotated_results_df_ALDE_sig_sub into dataframes with fewer rows and managed to annotate them.

2. Plotting Heatmap:

> daa_results_df_ALDE1 <- daa_results_df_ALDE_sig_sub[c(1:100),]
> daa_results_df_ALDE2 <- daa_results_df_ALDE_sig_sub[c(101:200),]
> daa_results_df_ALDE3 <- daa_results_df_ALDE_sig_sub[c(201:204),]

> daa_annotated_results_df_ALDE1 <- pathway_annotation(pathway = "KO", daa_results_df = daa_results_df_ALDE1, ko_to_kegg = TRUE)
> daa_annotated_results_df_ALDE2 <- pathway_annotation(pathway = "KO", daa_results_df = daa_results_df_ALDE2, ko_to_kegg = TRUE)
> daa_annotated_results_df_ALDE3 <- pathway_annotation(pathway = "KO", daa_results_df = daa_results_df_ALDE3, ko_to_kegg = TRUE)

> daa_annotated_results_df_ALDE <- rbind(daa_annotated_results_df_ALDE1, daa_annotated_results_df_ALDE2, daa_annotated_results_df_ALDE3)

> pathway_heatmap(
    abundance = daa_annotated_results_df_ALDE,
    metadata = metadata1,
    group = "Group"
)

Error Received:

Error in pull():
! Can't extract columns past the end.
ℹ Location 1 doesn't exist.
ℹ There are only 0 columns.

---

Any assistance would be much appreciated. Thank you!

[cafferychen777]

Hello @aston9090 ,

It appears there's a misunderstanding about the data structure required by the pathway_heatmap() function. Here's a reference code that might clarify the expected input and steps:

Load necessary data
data("metacyc_abundance")
data("metadata")

Perform differential abundance analysis
metacyc_daa_results_df <- pathway_daa(abundance = metacyc_abundance %>%
column_to_rownames("pathway"),
metadata = metadata, group = "Environment", daa_method = "LinDA")

Annotate the results
annotated_metacyc_daa_results_df <- pathway_annotation(pathway = "MetaCyc",
daa_results_df = metacyc_daa_results_df, ko_to_kegg = FALSE)

Filter pathways with p-adjusted value less than 0.05
feature_with_p_0.05 <- metacyc_daa_results_df %>% filter(p_adjust < 0.05)

Prepare the abundance data for heatmap
pathway_heatmap(abundance = metacyc_abundance %>%
right_join(annotated_metacyc_daa_results_df %>%
select(all_of(c("feature","description"))), by = c("pathway" = "feature")) %>%
filter(pathway %in% feature_with_p_0.05$feature) %>%
select(-"pathway") %>%
column_to_rownames("description"), metadata = metadata, group = "Environment")

I recommend aligning your data structure and steps with the ones presented in this example. Please let me know if you have further questions or need assistance.

Best wishes!

[aston9090]

Thank you!
On another issue, may I know how I can change the colours of heatmap and the points in pathway_pca? Thanks!

[cafferychen777]

Dear @aston9090 ,

Sorry for the delayed reply,
Regarding the customization of color schemes:
I've taken the liberty to provide adjusted functions with customizable color schemes. You can directly replace the color values within the colors_choices and colors variables to suit your manuscript's theme.

In the pathway_pca function:

colors_choices <- c("#d93c3e", "#3685bc", ...)[1:levels]

pathway_pca <- function(abundance, metadata, group){
# due to NSE notes in R CMD check
PC1 = PC2 = Group = NULL
# Perform PCA on the abundance data, keeping the first two principal components
pca_axis <- stats::prcomp(t(abundance), center = TRUE, scale = TRUE)$x[,1:2]

# Calculate the proportion of total variation explained by each PC
pca_proportion <- stats::prcomp(t(abundance), center = TRUE, scale = TRUE)$sdev[1:2]/sum(stats::prcomp(t(abundance), center = TRUE, scale = TRUE)$sdev)*100

# Combine the PCA results with the metadata information
pca <- cbind(pca_axis, metadata %>% select(all_of(c(group))))
pca$Group <- pca[,group]

levels <- length(levels(factor(pca$Group)))



colors_choices <- c("#d93c3e", "#3685bc","#208A42","#89288F","#F47D2B",
          "#FEE500","#8A9FD1","#C06CAB","#E6C2DC","#90D5E4",
          "#89C75F","#F37B7D","#9983BD","#D24B27","#3BBCA8",
          "#6E4B9E","#0C727C", "#7E1416","#D8A767","#3D3D3D")[1:levels]


# Create a ggplot object for the PCA scatter plot
Fig1a.taxa.pca <- ggplot2::ggplot(pca,ggplot2::aes(PC1,PC2))+
 ggplot2::geom_point(size=2,ggplot2::aes(color=Group),show.legend = T)+
 ggplot2::scale_color_manual(values=colors_choices)+
 ggplot2::stat_ellipse(ggplot2::aes(color = Group),fill="white",geom = "polygon",
              level=0.95,alpha = 0.01,show.legend = F)+
 ggplot2::labs(x=paste0("PC1(",round(pca_proportion[1],1),"%)"),y=paste0("PC2(",round(pca_proportion[2],1),"%)"),color = group)+
 ggplot2::theme_classic()+
 ggplot2::theme(axis.line=ggplot2::element_line(colour = "black"),
       axis.title=ggplot2::element_text(color="black",face = "bold"),
       panel.grid.major = ggplot2::element_blank(),
       panel.grid.minor = ggplot2::element_blank(),
       panel.background = ggplot2::element_blank(),
       axis.text = ggplot2::element_text(color="black",size=10,face = "bold"),
       legend.text = ggplot2::element_text(size = 16, face = "bold"),
       legend.title = ggplot2::element_text(size = 16, face = "bold"))+
 ggplot2::geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
 ggplot2::geom_hline(yintercept = 0, linetype = "dashed", color = "black")

# Create a ggplot object for the density plot of PC1
# Plot the density of PC1 colored by Group
Fig1a.taxa.pc1.density <-
 ggplot2::ggplot(pca) +
 ggplot2::geom_density(ggplot2::aes(x=PC1, group=Group, fill=Group), # Plot the density of PC1
              color="black", alpha=1, position = 'identity',show.legend = F) +
 ggplot2::scale_fill_manual(values=colors_choices) + # Manually set the fill color for each group
 ggplot2::theme_classic() + # Set the classic theme
 ggplot2::scale_y_discrete(expand = c(0,0.001)) + # Scale the y-axis with a small expansion to improve appearance
 ggplot2::labs(x=NULL, y=NULL) + # Remove x and y labels
 ggplot2::theme(axis.text.x=ggplot2::element_blank(), # Remove x axis text
       axis.ticks.x = ggplot2::element_blank()) # Remove x axis ticks

# Plot the density of PC2 colored by Group
Fig1a.taxa.pc2.density <-
 ggplot2::ggplot(pca) +
 ggplot2::geom_density(ggplot2::aes(x=PC2, group=Group, fill=Group), # Plot the density of PC2
              color="black", alpha=1, position = 'identity',show.legend = F) +
 ggplot2::scale_fill_manual(values=colors_choices) + # Manually set the fill color for each group
 ggplot2::theme_classic() + # Set the classic theme
 ggplot2::scale_y_discrete(expand = c(0,0.001)) + # Scale the y-axis with a small expansion to improve appearance
 ggplot2::labs(x=NULL, y=NULL) + # Remove x and y labels
 ggplot2::theme(axis.text.y = ggplot2::element_blank(), # Remove y axis text
       axis.ticks.y = ggplot2::element_blank()) + # Remove y axis ticks
 ggplot2::coord_flip() # Flip the plot to change it from a horizontal to a vertical orientation

# Combine the two plots into a single plot with the PC1 plot on top and the PC2 plot on the right
Fig1a.taxa.pca %>%
 aplot::insert_top(Fig1a.taxa.pc1.density, height = 0.3) %>% # Insert the PC1 plot on top with a height of 0.3
 aplot::insert_right(Fig1a.taxa.pc2.density, width=0.3) %>% # Insert the PC2 plot on the right with a width of 0.3
 ggplotify::as.ggplot() # Convert the combined plot to a ggplot object
}

And in the pathway_heatmap function:

colors <- c("#d93c3e", "#3685bc", ...)

pathway_heatmap <- function(abundance, metadata, group) {
# Heatmaps use color changes to visualize changes in values. However, if the
# data for plotting the heat map are too different, for example, if the heat
# map is plotted using gene expression data, gene1 is expressed above 1000 in
# all samples and gene2 is expressed between 1-10 in all samples, it is
# difficult to plot the heat map with small changes in the expression of two
# genes in different samples by the colors to reflect. Therefore, when
# plotting a heat map, we usually normalize the gene expression data, that
# is, we subtract the mean value of each gene expression from the expression
# of this gene in all samples and divide it by its standard deviation, and
# this normalization is called standard normalization or Z-score processing.
# The processed values are reduced equally, and the expression of each gene
# in all samples becomes a set of values with a mean of 0 and a standard
# deviation of 1. At this point, the plotted heat map gives a good indication
# of the variation in expression of all genes across samples.

# Check that 'group' is a column in 'metadata'
if (!group %in% colnames(metadata)) {
 stop(paste("group:", group, "must be a column in metadata"))
}

# Find the column in metadata that matches the column names of abundance
sample_name_col <- colnames(metadata)[sapply(colnames(metadata), function(x) all(colnames(abundance) %in% metadata[[x]]))]
metadata$sample_name <- metadata %>% select(all_of(c(sample_name_col))) %>% pull()

if (!all(colnames(abundance) %in% metadata$sample_name)) {
 stop("Samples in abundance and metadata must match")
}

# Now sample_name_col contains the column name in metadata that stores the sample names

z_abundance <- t(apply(abundance, 1, scale))
colnames(z_abundance) <- colnames(abundance)

# Convert the abundance matrix to a data frame
z_df <- as.data.frame(z_abundance)

metadata <- metadata %>% as.data.frame()

# Order the samples based on the environment information
ordered_metadata <- metadata[order(metadata[, group]),]
ordered_sample_names <- ordered_metadata$sample_name
order <- ordered_metadata$sample_name
ordered_group_levels <- ordered_metadata %>% select(all_of(c(group))) %>% pull()


# Convert the abundance data frame to a long format
long_df <- z_df %>%
 tibble::rownames_to_column() %>%
 tidyr::pivot_longer(cols = -rowname,
                     names_to = "Sample",
                     values_to = "Value") %>% left_join(metadata %>% select(all_of(c("sample_name",group))), by = c("Sample" = "sample_name"))

# Set the order of the samples in the heatmap
long_df$Sample <- factor(long_df$Sample, levels = order)

# Compute breaks from the data
breaks <- range(long_df$Value, na.rm = TRUE)

colors <- c("#d93c3e", "#3685bc", "#6faa3e", "#e8a825", "#c973e6", "#ee6b3d", "#2db0a7", "#f25292")

# Create the heatmap using ggplot
p <-
 ggplot2::ggplot(data = long_df,
                 mapping = ggplot2::aes(x = Sample, y = rowname, fill = Value)) +
 ggplot2::geom_tile() +
 ggplot2::scale_fill_gradient2(low = "#0571b0", mid = "white", high = "#ca0020", midpoint = 0) +
 ggplot2::labs(x = NULL, y = NULL) +
 ggplot2::scale_y_discrete(expand = c(0, 0), position = "left") +
 ggplot2::scale_x_discrete(expand = c(0, 0)) +
 # Customize the appearance of the heatmap
 ggplot2::theme(
   axis.text.x = ggplot2::element_blank(),
   axis.text.y = ggplot2::element_text(size = 12, color = "black"),
   axis.ticks = ggplot2::element_blank(),
   axis.text = ggplot2::element_text(
     color = "black",
     size = 10,
     face = "bold"
   ),
   panel.spacing = unit(0, "lines"),
   legend.title = ggplot2::element_text(size = 12, color = "black",face = "bold"),
   legend.text = ggplot2::element_text(size = 12, color = "black",face = "bold"),
   panel.background = ggplot2::element_blank(),
   legend.margin = ggplot2::margin(l = 0, unit = "cm"),
   strip.text = element_text(size = 12, face = "bold")
 ) +
 # Add a color bar to the heatmap
 ggplot2::guides(
   fill = ggplot2::guide_colorbar(
     direction = "vertical",
     reverse = F,
     barwidth = unit(0.6, "cm"),
     barheight = unit(9, "cm"),
     title = "Z Score",
     title.position = "top",
     title.hjust = -1,
     ticks = TRUE,
     label = TRUE
   )
 ) + ggh4x::facet_nested(cols = vars(!!sym(group)), space = "free", scale = "free", switch = "x", strip =strip_nested(background_x = elem_list_rect(fill = colors)))

# Print the ordered sample names and group levels
cat("The Sample Names in order from left to right are:\n")
cat(ordered_sample_names, sep = ", ")
cat("\n")

cat("The Group Levels in order from left to right are:\n")
cat(ordered_group_levels, sep = ", ")
cat("\n")

return(p)
}

You can replace the color codes inside the c(...) with your desired colors. The rest of the function logic remains the same and should adapt to the new color scheme you've chosen.

I hope these solutions address your concerns. If you have further questions or need more clarifications, please don't hesitate to ask.

Best regards,
Chen YANG