R - Bonus Content

Additional visualizations for gene/sample level QC assessment

Count boxplots

To evaluate the difference between the count distributions before and after normalization, let’s extract the raw counts and the rlog normalized counts, and coerce them to tibbles.

raw_counts = as_tibble(counts(dds), rownames = 'id')
norm_counts = as_tibble(assay(rld), rownames = 'id')

raw_counts

# A tibble: 55,492 × 7
   id                 sample_A sample_B sample_C sample_D sample_E sample_F
   <chr>                 <int>    <int>    <int>    <int>    <int>    <int>
 1 ENSMUSG00000000001     1041      905     1296     3481     1283     1921
 2 ENSMUSG00000000003        0        0        0        0        0        0
 3 ENSMUSG00000000028     1043     1232     1664     2690     1825     2720
 4 ENSMUSG00000000031     1819      914     1618     8618     1350     1222
 5 ENSMUSG00000000037       19       11       18       48       37       29
 6 ENSMUSG00000000049       18        1        4       24        1        1
 7 ENSMUSG00000000056    14972    14768    21026    22962    22263    23622
 8 ENSMUSG00000000058        1        0        0        0        0        0
 9 ENSMUSG00000000078      888      607      911     1689      738     1180
10 ENSMUSG00000000085      402      483      744      898      811     1261
# ℹ 55,482 more rows

norm_counts

# A tibble: 16,249 × 7
   id                 sample_A sample_B sample_C sample_D sample_E sample_F
   <chr>                 <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
 1 ENSMUSG00000000001    10.5     10.4     10.4     10.8     10.4     10.6 
 2 ENSMUSG00000000028    10.6     10.7     10.7     10.7     10.8     11.0 
 3 ENSMUSG00000000031    11.2     10.5     10.7     11.9     10.6     10.3 
 4 ENSMUSG00000000037     4.64     4.56     4.58     4.66     4.72     4.62
 5 ENSMUSG00000000049     3.02     2.82     2.84     2.94     2.81     2.81
 6 ENSMUSG00000000056    14.3     14.3     14.3     13.9     14.4     14.2 
 7 ENSMUSG00000000078    10.1      9.77     9.86     9.93     9.68     9.90
 8 ENSMUSG00000000085     9.28     9.41     9.53     9.23     9.60     9.80
 9 ENSMUSG00000000088    10.4     10.6     10.7     10.5     10.5     10.2 
10 ENSMUSG00000000093     3.88     3.89     3.93     3.91     3.94     3.99
# ℹ 16,239 more rows

Next, this is the perfect opportunity to use tidyr::pivot_longer() because we want to use ggplot() for the bar plots, but the data is currently in the wide form, not the tidy form!

tidy_raw = tidyr::pivot_longer(raw_counts, -id, names_to = 'sample', values_to = 'counts')
tidy_norm = tidyr::pivot_longer(norm_counts, -id, names_to = 'sample', values_to = 'counts')

We should also join in the sample metadata so that we can color on the sample groups.

samplesheet_tbl = as_tibble(colData(dds), rownames = 'sample')

tidy_raw = tidy_raw %>% left_join(samplesheet_tbl, by = 'sample')
tidy_norm = tidy_norm %>% left_join(samplesheet_tbl, by = 'sample')

raw_boxplot = ggplot(tidy_raw, aes(x = sample, y = log2(counts), fill = condition)) +
    geom_boxplot(notch = TRUE) +
    labs(
        title = 'Raw Counts',
        x = 'Sample',
        y = 'log2 Counts') +
    theme_bw() + theme(axis.text.x = element_text(angle = 90))
raw_boxplot

Warning: Removed 226538 rows containing non-finite outside the scale range
(`stat_boxplot()`).

Question

Why did we get that warning about non-finite values?

norm_boxplot = ggplot(tidy_norm, aes(x = sample, y = counts, fill = condition)) +
    geom_boxplot(notch = TRUE) +
    labs(
        title = 'rlog Normalized Counts',
        x = 'Sample',
        y = 'rlog Counts') +
    theme_bw() + theme(axis.text.x = element_text(angle = 90))
norm_boxplot

Observe that the normalized plots truly do appear normalized; their means are more uniform. Let’s go ahead and save these plots in our outputs/figures folder.

ggsave(filename = 'outputs/figures/Boxplot_raw_condition.pdf', plot = raw_boxplot, height = 6, width = 6)
ggsave(filename = 'outputs/figures/Boxplot_rlog_condition.pdf', plot = norm_boxplot, height = 6, width = 6)

Heatmaps

Let’s create a heatmap based on the distance between pair-wise sample expression profiles. This is another vantage point of how similar and dissimilar the samples are from one another. To get started, we actually want a plain matrix, with out a column for the gene IDs, as we did for the boxplot.

norm_mat = assay(rld)
head(norm_mat)

                    sample_A  sample_B  sample_C  sample_D  sample_E  sample_F
ENSMUSG00000000001 10.514813 10.366709 10.419463 10.840373 10.410449 10.578771
ENSMUSG00000000028 10.604461 10.734506 10.735026 10.682714 10.820938 10.990999
ENSMUSG00000000031 11.160276 10.498747 10.742763 11.861617 10.578156 10.298022
ENSMUSG00000000037  4.642433  4.555500  4.578934  4.656950  4.719012  4.620204
ENSMUSG00000000049  3.017478  2.820198  2.843894  2.936455  2.814833  2.811619
ENSMUSG00000000056 14.321672 14.284652 14.337197 13.904238 14.393912 14.235971

Next, we’ll use the dist() function on the transpose of norm_mat to compute the distance. We will use the default Euclidean distance.

dist_mat = dist(t(norm_mat), upper = TRUE)

Next, we’ll plot a very simple pheatmap() using this matrix. But let’s check the class of dist_mat first, because the input to pheatmap() needs to be a matrix.

class(dist_mat)

[1] "dist"

# Have to coerce it
dist_mat = as.matrix(dist_mat)

dist_heatmap = pheatmap(
    mat = dist_mat,
    cluster_rows = TRUE,
    cluster_cols = TRUE,
    show_rownames = TRUE,
    show_colnames = TRUE
)
dist_heatmap

This is nice, but there are a few tweaks we might consider:

1. The color scale used by default is divergent, but our distances are values in [0, Inf), so a linear color scale would be more appropriate here.
2. We might like for the sample metadata to be included so that we can easily tell if the samples cluster by their condition.

To accomplish the first change, we’ll use the colorRampPalette() function from the RColorBrewer package. This is the package to use for creating color scales of all varieties. We highly recommend exploring the package website and documentation.

colors = colorRampPalette(brewer.pal(9, 'Blues'))(50)

Now let’s use that in the pheatmap() call using the color parameter:

dist_heatmap = pheatmap(
    mat = dist_mat,
    color = colors,
    cluster_rows = TRUE,
    cluster_cols = TRUE,
    show_rownames = TRUE,
    show_colnames = TRUE
)
dist_heatmap

This is a lot more reasonable. More distant samples are a deeper blue, while more similar samples are closer to white. Next, let’s add some annotation data.

annotation_tbl = samplesheet_tbl %>%
    dplyr::select(sample, condition) %>%
    tibble::column_to_rownames(var = 'sample')

dist_heatmap = pheatmap(
    mat = dist_mat,
    color = colors,
    cluster_rows = TRUE,
    cluster_cols = TRUE,
    show_rownames = TRUE,
    show_colnames = TRUE,
    annotation_col = annotation_tbl
)
dist_heatmap

And now we have our conditions as a colored annotation bar along the columns.

Sources

• HBC QC tutorial: https://hbctraining.github.io/DGE_workshop/lessons/03_DGE_QC_analysis.html
• Detailed Heatmap tutorial from Galaxy: https://training.galaxyproject.org/training-material/topics/transcriptomics/tutorials/rna-seq-viz-with-heatmap2/tutorial.html

---

Session Info

sessionInfo()

R version 4.4.0 (2024-04-24)
Platform: x86_64-apple-darwin20
Running under: macOS Sonoma 14.4.1

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/Detroit
tzcode source: internal

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] biomaRt_2.60.1              data.table_1.15.4          
 [3] RColorBrewer_1.1-3          pheatmap_1.0.12            
 [5] ggrepel_0.9.5               lubridate_1.9.3            
 [7] forcats_1.0.0               stringr_1.5.1              
 [9] dplyr_1.1.4                 purrr_1.0.2                
[11] readr_2.1.5                 tidyr_1.3.1                
[13] tibble_3.2.1                ggplot2_3.5.1              
[15] tidyverse_2.0.0             DESeq2_1.44.0              
[17] SummarizedExperiment_1.34.0 Biobase_2.64.0             
[19] MatrixGenerics_1.16.0       matrixStats_1.3.0          
[21] GenomicRanges_1.56.1        GenomeInfoDb_1.40.1        
[23] IRanges_2.38.1              S4Vectors_0.42.1           
[25] BiocGenerics_0.50.0         knitr_1.47                 
[27] rmarkdown_2.27             

loaded via a namespace (and not attached):
 [1] DBI_1.2.3               httr2_1.0.2             rlang_1.1.4            
 [4] magrittr_2.0.3          compiler_4.4.0          RSQLite_2.3.7          
 [7] png_0.1-8               vctrs_0.6.5             pkgconfig_2.0.3        
[10] crayon_1.5.3            fastmap_1.2.0           dbplyr_2.5.0           
[13] XVector_0.44.0          labeling_0.4.3          utf8_1.2.4             
[16] tzdb_0.4.0              UCSC.utils_1.0.0        bit_4.0.5              
[19] xfun_0.44               zlibbioc_1.50.0         cachem_1.1.0           
[22] jsonlite_1.8.8          progress_1.2.3          blob_1.2.4             
[25] highr_0.11              DelayedArray_0.30.1     BiocParallel_1.38.0    
[28] parallel_4.4.0          prettyunits_1.2.0       R6_2.5.1               
[31] bslib_0.7.0             stringi_1.8.4           jquerylib_0.1.4        
[34] Rcpp_1.0.13             Matrix_1.7-0            timechange_0.3.0       
[37] tidyselect_1.2.1        rstudioapi_0.16.0       abind_1.4-5            
[40] yaml_2.3.8              codetools_0.2-20        curl_5.2.1             
[43] lattice_0.22-6          withr_3.0.1             KEGGREST_1.44.1        
[46] evaluate_0.23           BiocFileCache_2.12.0    xml2_1.3.6             
[49] Biostrings_2.72.1       filelock_1.0.3          pillar_1.9.0           
[52] BiocManager_1.30.23     generics_0.1.3          hms_1.1.3              
[55] munsell_0.5.1           scales_1.3.0            glue_1.7.0             
[58] tools_4.4.0             locfit_1.5-9.10         grid_4.4.0             
[61] AnnotationDbi_1.66.0    colorspace_2.1-1        GenomeInfoDbData_1.2.12
[64] cli_3.6.2               rappdirs_0.3.3          fansi_1.0.6            
[67] S4Arrays_1.4.1          gtable_0.3.5            sass_0.4.9             
[70] digest_0.6.35           SparseArray_1.4.8       farver_2.1.2           
[73] memoise_2.0.1           htmltools_0.5.8.1       lifecycle_1.0.4        
[76] httr_1.4.7              bit64_4.0.5            

---

These materials have been adapted and extended from materials listed above. These are open access materials distributed under the terms of the Creative Commons Attribution license (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

---

Previous lesson	Top of this lesson	To wrap-up