Hi all,
I am working with RNA-Seq data and need to control for a variable which can be measured via HOX gene expression. The study is a simple case/control study with a treated and untreated condition.
I have performed some clustering on the dds object (code below) and generated a correlation heatmap of HOX gene expression. To my mind, there are roughly 4 clusters. The labels aren't clear, but all of my controls seem to fall into a single cluster (they have shorter names than the cases).
When I assign cluster membership to these samples in the metadata, all of the untreated controls fall into cluster 1, but the treated controls seem to be split across clusters 2-3. This doesn't seem to bear any resemblance to what is shown on the heatmap, nor does it really make biological sense.
Would you have any insights? Is there a better way to extract clusters from the dendrogram?
Many thanks,
Rob.
Code Below (All code written in Rstudio):
dds <- DESeqDataSetFromMatrix(counts, meta, design = ~condition)
dds <- estimateSizeFactors(dds)
Performing HOX Clustering (HOXA10-HOXA9 removed as there were 0 counts)
vst_dds <- vst(dds, blind=TRUE)
vst_mat <- assay(vst_dds)
vst_mat_hox <- vst_mat[grepl("^HOX", row.names(vst_mat)), ]
vst_mat_hox <- vst_mat_hox[row.names(vst_mat_hox) != "HOXA10-HOXA9", ]
vst_cor_hox <- cor(vst_mat_hox)
Generate the clustering heatmap
library(pheatmap)
hox_heatmap <- pheatmap(vst_cor_hox)
Assign clusters to the metadata
meta$hox_group <- cutree(as.hclust(hox_heatmap$tree_row), k = 4)
meta$hox_group <- as.factor(meta$hox_group)
Example of HOX gene group membership of treated controls
