Help with plotting DEG-specific heatmap

Question

[marina.wakid]

Hi all,

So this code previously worked but for some reason it isn't anymore and I don't know why. I want to plot a heatmap of DEGs only (side note I know vst is superior, I just want to compare rlog and vst). Can someone please explain what's wrong and how to fix it? Many thanks!

colData table looks like this, with more columns of data and up to 48 rows (48 samples). I have two experimental groups, CTRL and SA:

res_order_FDR table (these are my DEGs):

# Heatmaps - (Hierarchical clustering)
# generating Heatmaps (Hierarchical clustering) with just DEGs

# read in table with only genes reaching FDR
res_order_FDR <- read.csv("DEGS_FDR0.1_tableforheatmap.csv", row.names=NULL, fill=TRUE)
res_order_FDR <- as.data.frame(res_order_FDR)

# take the regularized log
rld <- rlog(dds, blind = FALSE)

#isolate samples from groups of interest
ind_to_keep <- c(which(colData(rld)$Group=="CTRL"), which(colData(rld)$Group=="SA"))

# set up gene expression matrix
mat1 <- assay(rld)[rownames(res_order_FDR_05), ind_to_keep]

This is where I get the problem I cannot fix (nor really understand)

Error in assay(rld)[rownames(res_order_FDR), ind_to_keep] : subscript out of bounds

# scale matrix by each col. values
mat_scaled = t(apply(mat1, 1, scale))

# set up colors for heatmap
col = colorRamp2(c(-3, 0, 3), c("blue", "white", "red"))
cols1 <- brewer.pal(11, "Paired")

# subset coldata for samples in CTRL and SA groups
colData_sub <- colData(dds)[ind_to_keep, ]

# set up annotation bar for samples
ha1 = HeatmapAnnotation(Group = colData_sub$Group,
                        col = list(Group = c("CTRL" = cols1[1], "SA" = cols1[2])),
                        show_legend = TRUE)

# set up column annotation labels (samples)
ha = columnAnnotation(x = anno_text(colData_sub$SRR,
                                    which="column", rot = 45,
                                    gp = gpar(fontsize = 10)))

# generate heatmap object
ht1 = Heatmap(mat_scaled, name = "Expression", col = col,
              top_annotation = c(ha1),
              bottom_annotation = c(ha),
              show_row_names = FALSE,
              show_column_names = FALSE,
              cluster_columns = FALSE)
# plot the heatmap
draw(ht1, row_title = "Genes", column_title = "Hierachical clustering of DEGs (padj<0.1)")

As an alternative, I have also tried:

# set up gene expression matrix
mat1 <- assay(rld)[rownames(res_order_FDR), rownames(colData)]

I end up with the same error as before

[Trivas]

This type of error is usually because your indexing rownames(res_order_FDR_05) or ind_to_keep is larger than your object assay(rld). Use a combination of the functions dim, nrow, ncol, and length to see which part of your indexing code is out of bounds.

As a note, your code does not generate the object res_order_FDR_05.

[marina.wakid]

Thank you Trivas!

Answer 1

[jared.andrews07]

library(dittoSeq)

# take the regularized log
rld <- rlog(dds, blind = FALSE)
assay(dds, "vst") <- assay(rld)

dittoHeatmap(dds, genes = rownames(res_order_FDR_05), assay = "vst", 
             cells.use = colData(dds)$Group %in% c("CTRL", "SA"), annot.by = c("Group", "SRR"))

Should get you pretty close.

Your original code has a bunch of issues. I'd go through it a line at a time to ensure you're actually getting what you expect for command.

[marina.wakid]

Thank you Jared!