This vignette shows the relevant steps to QC RNA-seq data. The rationale of this workflow is to assess whether the overall experiment meets desired quality standards and to detect low quality samples.

Preparations

We start with loading relevant packages for the analysis.

library("RNAseqQC")
library("DESeq2")
library("ensembldb")
library("dplyr")
library("ggplot2")
library("purrr")
library("tidyr")
library("tibble")
library("magrittr")

Input data

The input to the workflow consists of a genes \(\times\) samples count matrix and data.frame of metadata annotating the samples, where the number of metadata rows must equal the number of matrix columns. In addition, the count matrix must have column names and the row names must be ENSEMBL gene IDs.

As an example, we are analyzing a data set in which T47D cells with different mutation statuses were treated with E2 (estradiol) or vehicle (Bahreini et al. 2017). The input data are stored in this package as a DESeqDataSet called T47D. For demonstration purposes, we extract the count matrix and sample metadata to show how to construct an annotated DESeqDataSet, which is required for running this package.

count_mat <- counts(T47D)
meta <- data.frame(colData(T47D))

# count matrix; rownames must be ENSEMBL gene IDs
count_mat[head(which(rowSums(count_mat) > 0)), 1:10]
#>                 veh_WT_rep1 veh_WT_rep2 veh_WT_rep3 veh_WT_rep4 veh_D538G_rep1
#> ENSG00000223972         772         386         487         929            568
#> ENSG00000278267         290         391         311         623            319
#> ENSG00000227232       11914       13539       14108       14427           9616
#> ENSG00000284332           0           0           0          63              0
#> ENSG00000243485          89         353         153         352              0
#> ENSG00000237613           0           0           0         150              0
#>                 veh_D538G_rep2 veh_D538G_rep4 veh_Y537S_rep1 veh_Y537S_rep2
#> ENSG00000223972            210            891            518            760
#> ENSG00000278267            318            267            354            283
#> ENSG00000227232          10344          10747          11011           9569
#> ENSG00000284332              0             52              0              0
#> ENSG00000243485            253            501            103            137
#> ENSG00000237613              0              0              0              0
#>                 veh_Y537S_rep3
#> ENSG00000223972            406
#> ENSG00000278267            324
#> ENSG00000227232          13232
#> ENSG00000284332              0
#> ENSG00000243485            118
#> ENSG00000237613              0

# metadata of the samples, where row i corresponds to column i in the count matrix
meta
#>                treatment mutation replicate
#> veh_WT_rep1          veh       WT      rep1
#> veh_WT_rep2          veh       WT      rep2
#> veh_WT_rep3          veh       WT      rep3
#> veh_WT_rep4          veh       WT      rep4
#> veh_D538G_rep1       veh    D538G      rep1
#> veh_D538G_rep2       veh    D538G      rep2
#> veh_D538G_rep4       veh    D538G      rep4
#> veh_Y537S_rep1       veh    Y537S      rep1
#> veh_Y537S_rep2       veh    Y537S      rep2
#> veh_Y537S_rep3       veh    Y537S      rep3
#> veh_Y537S_rep4       veh    Y537S      rep4
#> veh_D538G_rep3       veh    D538G      rep3
#> E2_WT_rep1            E2       WT      rep1
#> E2_WT_rep2            E2       WT      rep2
#> E2_WT_rep3            E2       WT      rep3
#> E2_WT_rep4            E2       WT      rep4
#> E2_Y537S_rep1         E2    Y537S      rep1
#> E2_Y537S_rep2         E2    Y537S      rep2
#> E2_Y537S_rep3         E2    Y537S      rep3
#> E2_Y537S_rep4         E2    Y537S      rep4
#> E2_D538G_rep1         E2    D538G      rep1
#> E2_D538G_rep2         E2    D538G      rep2
#> E2_D538G_rep3         E2    D538G      rep3
#> E2_D538G_rep4         E2    D538G      rep4

Create DESeqDataSet

We make a DESeqDataSet from the input count matrix and metadata. This step also requires to specify a record identifier from the AnnotationHub package, that allows to download an EnsDb object to annotate gene IDs and store them in the rowData slot. We choose the latest homo sapiens AnnotationHub record.

dds <- make_dds(counts = count_mat, metadata = meta, ah_record = "AH89426")

Records for other species can be found by running the following code. Ideally, the AnnotationHub record should match the genome build against which the sequencing reads were aligned.

mcols(AnnotationHub::AnnotationHub()) %>%
  as_tibble(rownames = "record_id") %>%
  dplyr::filter(rdataclass == "EnsDb")

QC plots on raw count matrix

Total sample counts

A good metric to start quality control is to look at the total number of counts for each sample. This is also referred to as library size and we typically expect all samples to have total counts within the same order of magnitude.

plot_total_counts(dds)

Library complexity

For each sample, it is shown what fraction of counts is taken up by what fraction of genes. Samples showing a different library complexity than the rest might be considered low quality. In our case, all samples have very similar complexity.

plot_library_complexity(dds)

Gene detection

Another interesting quantity is the number of detected genes for each sample. This number depends on the detection threshold, i.e. the minimum count required to call a gene detected and also on the used protocol (e.g. RNA-seq read counts vs. DRUG-seq UMI counts). Thus, different thresholds are of interest. Note that genes in this context refers to ENSEMBL genes and includes protein coding genes, lncRNAs, and various other biotypes.

plot_gene_detection(dds)

Gene biotypes

It is also helpful to stratify the total gene count by the different gene biotypes. The plot below shows for each sample the total (library size) normalized count for the major gene biotypes. In a high quality data set, we expect this plot to show almost horizontal lines. Deviations from the horizontal pattern may correspond to different biological groups in the data set, especially in case of different cell lines, but could also indicate batch effects or low quality samples.

plot_biotypes(dds)

Gene filtering

Before continuing with quality control, a useful intermediate step is to remove genes with low counts. This often substantially reduces the total number of genes, and thus the overall size of the data set and computation time. A good strategy is to determine the size of a biological group, i.e. all samples under the same biological condition, and filter (keep) genes with a certain count at least as often as the size of the biological group. In our case, each group is a combination of treatment and mutation status with 4 samples and we choose 5 counts as threshold.

dds <- filter_genes(dds, min_count = 5, min_rep = 4)

Variance stabilization

For downstream tasks like clustering or PCA, a transformation is necessary so that the variance of a gene does not depend on it’s mean, i.e. we want genes with low and high mean counts to have similar variances. A transformation that achieves this is called variance stabilizing and we use the vst function of DESeq2 for this task. To check if the variance is indeed stabilized, we plot for each gene the rank of the mean count versus the standard deviation. In the resulting plot, the red trend line should be relatively flat with respect to the scale on the y-axis, as can be seen for our example data.

vsd <- vst(dds)
mean_sd_plot(vsd)

Chromosomal expression

Chromosome expression plots can be useful to detect expression patterns related to the chromatin. For each chromosome, we make a heatmap where genes are arranged by chromosomal position along the x-axis and each row is a sample. The columns of the heatmap are scaled such that different genes can be compared to each other. Long red (blue) stretches along the x-axis indicate regions with increased (decreased) expression.

There can be several reasons for differential expression of a chromosomal region, for example duplication of a chromosome band or changes in chromatin accessibility. Also if the data comprise different cell lines or primary cancer samples, differential expression of a region or a whole chromosome could be observed.

In our example, chromosome 1 shows no obvious differentially expressed region, chromosome 5 shows differential expression related to the interaction of treatment and mutation status, and chromosome 14 is down regulated in samples with the D538G mutation.

map(c("1", "5", "14"), ~plot_chromosome(vsd, .x))
#> [[1]]

#> 
#> [[2]]

#> 
#> [[3]]

Replicate variability

A good quality metric on a per-sample basis is to check whether biological and/or technical replicates are closely correlated. To do this, we define groups of samples by specifying a column of the colData slot as grouping variable, and define all samples that belong to the same level as group. Subsequently, a reference sample is computed for each group by taking for each gene the median over all samples in the group. Finally, we make an MA-plot where each sample is plotted against its group reference.

In our example, we define a new grouping variable as the combination of treatment and mutation status and show replicate variability for the sample groups veh_WT and veh_Y537S. All shown replicates have high concordance with their respective reference. The only exception is sample veh_Y537S_rep3, which has an increased number down-regulated genes with high mean count.

# define new grouping variable
colData(vsd)$trt_mut <- paste0(colData(vsd)$treatment, "_", colData(vsd)$mutation)

ma_plots <- plot_sample_MAs(vsd, group = "trt_mut")
cowplot::plot_grid(plotlist = ma_plots[17:24], ncol = 2)

Clustering

We can also have a look at the clustering of samples by a distance metric such as euclidean distance or pearson correlation distance. It is useful to define annotation variables for the distance heatmap and relate appearing clusters to the annotation. In our example, samples cluster according to the biological groups defined in the previous section, highlighting the good quality of the data set.

# set seed to control random annotation colors
set.seed(1)
plot_sample_clustering(vsd, anno_vars = c("treatment", "mutation", "replicate"), distance = "euclidean")

PCA

The package also provides wrappers for principal component analysis (PCA) to check whether the main variation in the data is driven by biological or batch effects.

plot_pca(vsd, PC_x = 1, PC_y = 2, color_by = "treatment", shape_by = "mutation")

In addition to defining the PCs to plot and variables of interest for color or shape, we can also plot the PC loadings, i.e. the genes driving a principal axis. The plot_loadings function offers several options to annotate loadings, e.g. by a column of the rowData slot, by highlighting genes of interest or the top genes of a principal axis. It is also helpful to investigate if principal axes are driven by the genes GC content, indicating a potentially sequencing related effect.

pca_res <- plot_pca(vsd, show_plot = FALSE)
plot_loadings(pca_res, PC = 1, annotate_top_n = 5)

plot_loadings(pca_res, PC = 1, highlight_genes = c("CD34", "FLT1", "MAPT"))

plot_loadings(pca_res, PC = 4, color_by = "gene_biotype", show_plot = F)$plot +
  theme(legend.position = "bottom")

plot_loadings(pca_res, PC = 2, color_by = "gc_content")

To investigate multiple PCs simultaneously, we can also make a matrix of scatter plots, to plot all pairs of principal components against each other.

plot_pca_scatters(vsd, n_PCs = 5, color_by = "treatment", shape_by = "mutation")

Differential expression

Differential expression analysis is one of the most common tasks for investigating RNA-seq data and can be also used for quality assessment. Before doing differential testing, we start by actually plotting the expression pattern of interesting genes.

Plot a gene

Sometimes one has prior knowledge of which genes should be differentially expressed in a certain condition or one wants to get a better understanding of the genes driving a principal component. For such cases, this package offers methods to plot a gene of interest. In order to make counts between different samples comparable, sample-specific normalization factors need be computed prior to plotting. This is either done automatically when the DESeq function has already been run on the DESeqDataSet object, or alternatively directly via estimateSizeFacors.

dds <- estimateSizeFactors(dds)
plot_gene("CLEC2B", dds, x_var = "mutation", color_by = "treatment")

# modify the plot
plot_gene("CLEC2B", dds, x_var = "mutation", color_by = "treatment", show_plot = F)$plot + ggsci::scale_color_jco()

# a custom plot type
plot_data <- plot_gene("CLEC2B", dds, show_plot = F)$data
plot_data %>%
  ggplot(aes(treatment, count, color = mutation)) +
  geom_point() +
  geom_path(
    aes(x = treatment, y = mean_count, color = mutation, group = mutation),
    data = plot_data %>%
      group_by(treatment, mutation) %>%
      summarize(mean_count = mean(count), .groups = "drop")
  ) +
  labs(y = "log2(norm. count)", title="CLEC2B") +
  cowplot::theme_cowplot()

Sometimes it is also useful to inspect a large number of genes simultaneously to get a better understanding of the data. As example, we make a list of plots for the first 100 genes in the dataset and save them as a 4-page PDF file with 25 plots per page.

plots <- rownames(dds)[1:100] %>% 
  map(~plot_gene(.x, dds, x_var="mutation", color_by="treatment", show_plot = FALSE)$plot)
save_plots_to_pdf(plots, file="genes.pdf", ncol=5, nrow=5)

Differential testing

Differential testing allows to uncover genes that are differentially expressed for a certain comparison. We use DESeq2 for this purpose and refer the reader to the DESeq2 vignette on how to define a design formula and specify contrasts of interest. A good diagnostic plot for assessing the validity of differential testing assumptions is to plot the per-gene dispersion estimates versus the mean expression and check whether the red trend line shows an appropriate fit.

# design variables need to be factors
# since we update the design of the dds
# object, we have to do it manually
dds$mutation <- as.factor(dds$mutation)
dds$treatment <- as.factor(dds$treatment)
design(dds) <- ~ mutation + treatment

dds <- DESeq(dds, parallel = T)
plotDispEsts(dds)

Plot a testing result

Differential testing results can be visualized in various ways. A good diagnostic is the MA-plot that plots the mean normalized count versus the log₂-fold change. Ideally, the plot should look rather symmetric around the horizontal axis. In our example, there are genes with zero counts in one condition and relatively high counts in another, leading to very high log₂-fold changes. We thus use lfcShrink instead of the results function, which would be the default way to compute differential testing results.

There are three different colors in the plot: Light-gray points correspond to genes with very low counts that are not considered for differential testing by the independent filtering step. Points in different shades of gray are not significant and simultaneously display the number of neighboring points to get a better understand where the majority of points are located. Red points indicate differentially expressed genes and we also place labels for the top 5 up- and down-regulated genes.

# get testing results
de_res <- lfcShrink(dds, coef = "mutation_WT_vs_D538G", lfcThreshold = log2(1.5), type = "normal", parallel = TRUE)

# MA plot
plot_ma(de_res, dds, annotate_top_n = 5)

In addition, the plot_ma function allows to highlight genes, and they are colored by their respective category, i.e. excluded by independent filtering, not differentially expressed, and differentially expressed.

plot_ma(de_res, dds, highlight_genes = c("CLEC2B", "PAGE5", "GAPDH"))

SessionInfo

sessionInfo()
#> R version 4.2.0 (2022-04-22)
#> Platform: x86_64-apple-darwin17.0 (64-bit)
#> Running under: macOS Big Sur/Monterey 10.16
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
#> 
#> locale:
#> [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] magrittr_2.0.3              tibble_3.1.7               
#>  [3] tidyr_1.2.0                 purrr_0.3.4                
#>  [5] ggplot2_3.3.6               dplyr_1.0.9                
#>  [7] ensembldb_2.20.1            AnnotationFilter_1.20.0    
#>  [9] GenomicFeatures_1.48.3      AnnotationDbi_1.58.0       
#> [11] DESeq2_1.36.0               SummarizedExperiment_1.26.1
#> [13] Biobase_2.56.0              MatrixGenerics_1.8.0       
#> [15] matrixStats_0.62.0          GenomicRanges_1.48.0       
#> [17] GenomeInfoDb_1.32.2         IRanges_2.30.0             
#> [19] S4Vectors_0.34.0            BiocGenerics_0.42.0        
#> [21] RNAseqQC_0.1.4             
#> 
#> loaded via a namespace (and not attached):
#>   [1] colorspace_2.0-3         rjson_0.2.21             ellipsis_0.3.2          
#>   [4] gghighlight_0.3.3        circlize_0.4.15          XVector_0.36.0          
#>   [7] GlobalOptions_0.1.2      clue_0.3-61              rstudioapi_0.13         
#>  [10] hexbin_1.28.2            farver_2.1.0             affyio_1.66.0           
#>  [13] ggrepel_0.9.1            bit64_4.0.5              fansi_1.0.3             
#>  [16] xml2_1.3.3               codetools_0.2-18         splines_4.2.0           
#>  [19] doParallel_1.0.17        cachem_1.0.6             geneplotter_1.74.0      
#>  [22] knitr_1.39               jsonlite_1.8.0           Cairo_1.5-15            
#>  [25] Rsamtools_2.12.0         annotate_1.74.0          cluster_2.1.3           
#>  [28] vsn_3.64.0               dbplyr_2.2.0             png_0.1-7               
#>  [31] BiocManager_1.30.18      compiler_4.2.0           httr_1.4.3              
#>  [34] assertthat_0.2.1         Matrix_1.4-1             fastmap_1.1.0           
#>  [37] lazyeval_0.2.2           limma_3.52.1             cli_3.3.0               
#>  [40] htmltools_0.5.2          prettyunits_1.1.1        tools_4.2.0             
#>  [43] gtable_0.3.0             glue_1.6.2               GenomeInfoDbData_1.2.8  
#>  [46] affy_1.74.0              rappdirs_0.3.3           Rcpp_1.0.8.3            
#>  [49] jquerylib_0.1.4          vctrs_0.4.1              Biostrings_2.64.0       
#>  [52] preprocessCore_1.58.0    rtracklayer_1.56.0       iterators_1.0.14        
#>  [55] xfun_0.31                stringr_1.4.0            lifecycle_1.0.1         
#>  [58] restfulr_0.0.14          XML_3.99-0.9             zlibbioc_1.42.0         
#>  [61] scales_1.2.0             hms_1.1.1                ProtGenerics_1.28.0     
#>  [64] parallel_4.2.0           RColorBrewer_1.1-3       ComplexHeatmap_2.12.0   
#>  [67] yaml_2.3.5               curl_4.3.2               memoise_2.0.1           
#>  [70] sass_0.4.1               biomaRt_2.52.0           stringi_1.7.6           
#>  [73] RSQLite_2.2.14           highr_0.9                genefilter_1.78.0       
#>  [76] BiocIO_1.6.0             foreach_1.5.2            filelock_1.0.2          
#>  [79] BiocParallel_1.30.3      shape_1.4.6              rlang_1.0.2             
#>  [82] pkgconfig_2.0.3          bitops_1.0-7             evaluate_0.15           
#>  [85] lattice_0.20-45          patchwork_1.1.1          GenomicAlignments_1.32.0
#>  [88] labeling_0.4.2           cowplot_1.1.1            bit_4.0.4               
#>  [91] tidyselect_1.1.2         ggsci_2.9                R6_2.5.1                
#>  [94] generics_0.1.2           DelayedArray_0.22.0      DBI_1.1.2               
#>  [97] pillar_1.7.0             withr_2.5.0              survival_3.3-1          
#> [100] KEGGREST_1.36.0          RCurl_1.98-1.6           crayon_1.5.1            
#> [103] utf8_1.2.2               ggpointdensity_0.1.0     BiocFileCache_2.4.0     
#> [106] rmarkdown_2.14           GetoptLong_1.0.5         progress_1.2.2          
#> [109] locfit_1.5-9.5           grid_4.2.0               blob_1.2.3              
#> [112] digest_0.6.29            xtable_1.8-4             munsell_0.5.0           
#> [115] viridisLite_0.4.0        bslib_0.3.1

Bahreini, Amir, Zheqi Li, Peilu Wang, Kevin M Levine, Nilgun Tasdemir, Lan Cao, Hazel M Weir, et al. 2017. “Mutation Site and Context Dependent Effects of Esr1 Mutation in Genome-Edited Breast Cancer Cell Models.” Breast Cancer Research 19 (1): 1–10.

RNA-seq Quality Control

Frederik Ziebell