Hi all,

In the tutorial (https://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#data-quality-assessment-by-sample-clustering-and-visualization), the selected the top genes based on normalized counts, but then plot VST values (APPROACH 1).

### APPROACH 1
# Select top genes by rowMeans of normalized counts
select <- order(rowMeans(counts(dds, normalized=TRUE), 
                         na.rm=TRUE), 
                decreasing=TRUE)[1:100]

# Transform data and use 'select' indices for plotting
vsd <- vst(dds, blind=TRUE)
pheatmap(t(assay(vsd)[select, ]), 
         cluster_rows=TRUE, 
         cluster_cols=TRUE,
         ...)

If I want to plot the top highly expressed genes in a group of biological replicates (e.g., within the same cell type and same condition just across different donors), should I instead select top genes based on VST-transformed values, then plot VST (APPROACH 2)?

Thank you for your input!

### APPROACH 2
# 1) Transform data first
vsd <- vst(dds, blind=TRUE)

# 2) Select top genes by rowMeans of VST
gene_means_vst <- rowMeans(assay(vsd))
select <- order(gene_means_vst, decreasing=TRUE)[1:100]

# 3) Plot using the same VST data
pheatmap(t(assay(vsd)[select, ]), 
         cluster_rows=TRUE, 
         cluster_cols=TRUE,
         ...)

The choice of variable genes is typically part of a QC process. You select genes that in an unbiased fashion (unbiased because variance calculation does not know about sample-to-group relationship) separate samples and then feed these into clustering or PCA to identify outliers or batch effects.

You can plot "highly-expressed" genes, but I am not sure what this will give you. Absolute expression level (to me) does not mean much in RNA-seq since it is influenced by a lot of technical factors such as mappability, GC and amplification bias of the gene and the overall composition of the library.