Zehra Alibasic

RNA sequencing (RNA-seq) has become an indispensable tool for studying transcriptomes. By quantifying RNA molecules across conditions, it enables researchers to explore gene expression patterns, discover novel transcripts, and investigate molecular mechanisms underlying various biological phenomena, including diseases. RNA-seq starts with extracting RNA from cells, sequencing it, and producing millions of short reads. These reads are then aligned to a reference genome to produce a matrix of counts representing the number of reads per gene. However, raw counts alone are not sufficient for meaningful biological insights due to technical and biological variations. This is where Differential Expression Analysis (DEA) comes into play, quantifying and statistically validating differences in gene expression between experimental groups. There are many tools and bioinformatics programs used for the RNA-seq differential expression analysis. Why Use Bioconductor for RNA-Seq Analysis? Bioconductor is an R-based ecosystem specifically designed for bioinformatics. It includes packages like DESeq2, which handle the complexities of RNA-seq analysis, such as normalization, model fitting, and visualization. Bioconductor emphasizes reproducibility, a cornerstone of good scientific practice, while providing flexibility for diverse experimental designs. The statistical rigor of its tools ensures that results are both accurate and interpretable, making it an ideal choice for RNA-seq workflows. Introduction to the RNA-Seq Workflow Before diving into the differential expression analysis, it is essential to understand the broader RNA-Seq workflow and its foundational steps. RNA-Seq is not just about analyzing gene expression differences; it begins with careful experimental design and follows a structured pipeline to ensure reliable results. A robust experimental design is crucial, as errors at this stage can compromise the entire analysis. Biological replicates, typically three or more, are essential for capturing natural variability and improving the robustness of statistical tests. Minimizing batch effects and avoiding confounding variables, such as sex or age differences, is equally important to ensure that observed changes in gene expression reflect true biological differences rather than experimental artifacts. The workflow proceeds with sample preparation, where RNA is isolated, and non-essential components like ribosomal RNA are removed. The RNA is then reverse-transcribed into cDNA, fragmented, and prepared into libraries for sequencing. Sequencing generates millions of short reads, which are output as FASTQ files, representing the raw data for analysis. Quality control is performed on these files to assess sequencing accuracy and check for contamination or technical issues. Next, the reads are aligned to a reference genome, mapping them to their most likely origin in the genome. This alignment is particularly complex in RNA-Seq, as mRNA consists only of exons, and many reads span intronic boundaries, necessitating tools capable of aligning across introns. Once aligned, reads are quantified to produce a count matrix, a table where rows represent genes and columns represent samples. This count matrix forms the input for statistical analyses, where the primary goal is to determine whether gene expression differences between conditions are statistically significant, accounting for natural variability within sample groups. The results of the differential expression analysis include log2 fold changes and adjusted p-values for each gene, identifying those with significant expression differences. This process not only pinpoints individual genes of interest but also highlights broader patterns and pathways that may drive the observed biological phenomena. Bioconductor workflow Step 1: Data Preparation and Preprocessing Let’s try to perform DEA with Bioconductor using a publicly available RNA-seq dataset for fibrosis in Mus musculus (house mouse)*. Fibrosis is a pathological condition in which connective tissue replaces normal parenchymal tissue with connective tissue during wound healing, which ultimately reduces organ function and can cause failure and death**. The dataset used for this analysis contains raw RNA counts for each gene in 4 samples with induced fibrosis condition in mice through overexpression of the SMOC2 gene, which has been found to play one of the key roles in fibrosis condition, and 3 normal samples without the pathology.  The first step would be installing all of the packages needed for this analysis and then loading our raw count matrix. {r} # Install BiocManager (for Bioconductor packages like DESeq2) install.packages("BiocManager") # Install DESeq2 BiocManager::install("DESeq2") A raw count matrix from RNA-Seq is a table that contains the raw, unnormalized counts of sequencing reads that map to genomic features (e.g., genes, transcripts, or exons). These counts represent the initial data obtained after read alignment and quantification. After loading our raw data, it should look something like this: Step 2: Metadata Creation and Alignment Metadata describes the experimental conditions for each sample. In this case, metadata includes information about treatment groups “fibrosis” and “normal.” Proper metadata creation ensures each sample is correctly linked to its condition, forming the basis for statistical comparisons. {r} metadata <- data.frame( treatment = c("smoc2", "smoc2", "smoc2", "smoc2", "smoc2", "smoc2", "smoc2"), group =c("fibrosis", "fibrosis", "fibrosis", "fibrosis", "normal", "normal", "normal") ) rownames(metadata) <- colnames(data) The resulting metadata should look something like this: The alignment between the metadata and the count matrix is verified to prevent mismatches, which could lead to invalid results. This step reflects careful experimental planning, as the metadata directly determines how the statistical model interprets the data. Step 3: Creating the DESeq2 Dataset Using the DESeq2 package, we encapsulate the count data, metadata, and experimental design into a DESeq2 dataset object. The design formula specifies the condition of interest, such as the "group" variable in this case, and accounts for sources of variation. Output from this step can be observed below. {r} dds <- DESeqDataSetFromMatrix( countData = data, colData = metadata, design = ~ group ) This step sets the stage for differential expression analysis by modeling RNA-seq data using the negative binomial distribution. This distribution accommodates the variability typical of RNA-seq experiments, which often exceeds what simpler distributions like the Poisson can handle. Step 4: Data Normalization Raw RNA-seq counts are affected by sequencing depth and library size differences. To address this, DESeq2 calculates size factors for each sample, normalizing the data to make gene expression levels comparable across samples. {r} dds <- estimateSizeFactors(dds) normalized_counts <- counts(dds, normalized = TRUE) Normalization adjusts for biases and ensures that observed differences in gene expression are due to biological, rather than technical, factors. Skipping normalization could lead to inaccurate identification of differentially expressed genes, as systematic biases would confound the analysis. Step 5: Quality Assessment with VST, Heatmaps, and PCA After normalization, it is crucial to assess the data quality and explore relationships between samples. This ensures that biological replicates cluster together and conditions separate as expected. Variance Stabilizing Transformation (VST) RNA-seq data often exhibits a dependence of variance on mean expression, which can obscure patterns in the data. To address this, DESeq2 applies a Variance Stabilizing Transformation (VST) that makes variance approximately constant across expression levels. This transformation enhances the interpretability of downstream clustering and visualizations. {r} vsd <- vst(dds, blind = TRUE) The blind = TRUE argument ensures the transformation is not biased by experimental design, making it ideal for quality assessment. Heatmaps for Sample Clustering Heatmaps provide a visual representation of sample similarity based on pairwise correlation of gene expression. Samples are hierarchically clustered, revealing whether replicates group together and conditions separate clearly. A correlation heatmap shows the pairwise correlation coefficients between samples based on their VST-transformed expression data. These coefficients range from -1 (completely dissimilar) to 1 (perfectly similar). This heatmap reveals how closely related each sample is to the others in terms of overall gene expression. {r} pheatmap(cor(assay(vsd)), main = "Sample Correlation Heatmap") If replicates cluster tightly and conditions are well-separated, like it is the case in our heatmap above, it indicates good experimental reproducibility. Unexpected clustering patterns or outliers might require further investigation. Distance Heatmap A distance heatmap calculates the pairwise Euclidean distances between samples based on their VST-transformed data. It provides an alternative perspective on sample relationships, focusing on absolute differences in expression profiles rather than relative similarity. {r} pheatmap(as.matrix(dist(t(assay(vsd)))), main = "Sample Distance Heatmap") In the distance heatmap, smaller distances (dark red) indicate greater similarity between samples. Biological replicates should form tight clusters with small distances, while samples from different conditions should exhibit larger distances (dark blue). This visualization is particularly sensitive to outliers, as samples with unusual expression profiles will appear far from all others. Principal Component Analysis (PCA) PCA is another powerful technique to examine sample similarity. It reduces data dimensionality, with the first principal component (PC1) capturing the largest source of variation. In RNA-seq, PCA often reveals whether the experimental condition drives variation across samples. {r} plotPCA(vsd, intgroup = "group") In the PCA plot, biological replicates should cluster, while different conditions should separate along one of the principal components. Any failure to see separation might indicate confounding factors or subtle biological differences. We can see that in our sample groups, normal and fibrosis, separate well on PC1. This means that our condition corresponds to PC1, which represents 93% of the variance in our data, while 2% is explained by PC2. This is great since it seems that a lot of the variation in gene expression in the dataset can likely be explained by the differences between sample groups.  Step 6: Differential Expression Analysis The differential expression analysis aims to identify genes with significant expression differences between conditions, such as fibrosis and normal. DESeq2 fits the count data to a negative binomial model and performs statistical tests to compute log2 fold changes and adjusted p-values for each gene using the DESeq2()function in Bioconductor. {r} dds <- DESeq(dds) results <- results(dds, contrast = c("group", "fibrosis", "normal"), alpha = 0.05) The alpha parameter represents the significance threshold for determining which genes are considered "differentially expressed" (DE). An alpha of 0.05 means we are allowing up to a 5% chance of identifying a gene as DE when it is actually not DE (a false positive). This threshold corresponds to a p-value cutoff of 0.05. Genes with adjusted p-values (padj) below 0.05 are considered statistically significant. Exploring Gene Variability with Dispersion Plot The dispersion plot is an essential diagnostic tool for RNA-seq analysis. It shows how variability (dispersion) in gene expression relates to mean expression levels. {r} plotDispEsts(dds) Black Dots: Raw dispersion estimates for each gene. Red Line: Expected dispersion trend based on the mean expression. Blue Dots: Shrunken estimates used in the model for better accuracy. A good fit is indicated by points clustering around the red line. Large deviations might signal poor data quality or unaccounted variability. Shrinking dispersion values ensures more reliable differential expression results, especially for genes with low counts. MA Plot Another way to visually provide insight into gene behavior and ensure the reliability of differential expression analysis is an MA plot. The MA plot visualizes the relationship between the mean expression and the log2 fold change for all genes, helping to identify differentially expressed genes. {r} plotMA(results, ylim = c(-8, 8)) Genes with significant expression changes are highlighted, but genes with low counts may show large, unreliable fold changes. In this case, we can apply the shrinkage method. Shrinking reduces exaggerated fold changes for genes with low counts, improving reliability without affecting statistical significance. {r} results_shrunken <- lfcShrink( dds, contrast = c("group", "fibrosis", "normal"), type = "normal") plotMA(results_shrunken, ylim = c(-8, 8)) The difference between normal and shrunken MA plots can be observed below. Refining Results After identifying differentially expressed genes (DEGs) with DESeq2, we apply thresholds to focus on biologically significant findings to ensure selected genes show both meaningful expression changes and statistical significance. The lfcThreshold specifies the minimum log2 fold change (LFC) that is considered biologically meaningful. A value of 0.32 corresponds to an absolute fold change of approximately 1.25. The alpha parameter sets the significance level for the adjusted p-value (Benjamini-Hochberg False Discovery Rate), and a value of 0.05 means we accept a 5% false discovery rate. These thresholds can be adjusted based on the study's requirements. {r} results_lfc <- results( dds, contrast = c("group", "fibrosis", "normal"), lfcThreshold = 0.32, alpha = 0.05 ) To add biological context, Ensembl gene IDs are annotated using a very convenient Bioconductor package: the biomaRt: {r} mart <- useMart("ensembl", dataset = "mmusculus_gene_ensembl") annotations <- getBM( attributes =c("ensembl_gene_id", "external_gene_name", "description"), filters = "ensembl_gene_id", values =significant_results_lfc$ensgene, mart = mart ) Our significant differentially expressed genes with their biological descriptions and roles can be clearly observed now: Step 7: Visualization of Results Visualizations like volcano plots and heatmaps make it easier to interpret differential expression results. The expression heatmap visualizes the normalized counts of significant genes (in this case, the top 20 were chosen), highlighting patterns of gene expression across samples. {r} pheatmap(sig_norm_counts, color = heat_colors, scale ="row") Genes and samples with similar expression profiles group together, and it can be seen how significant genes are expressed differently between sample groups. A Z-score standardizes the gene expression values for each gene (row) to highlight relative expression patterns. Negative values (e.g., -2) show that expression is lower than average for that gene. Positive values (e.g., 2) show that expression is higher than average for that gene. Another way to visualize the expression levels of the most significant genes across all samples is the expression plot. {r} ggplot(top_20, aes(x = ensgene, y = normalized_counts, color = group)) + geom_point(size = 3, alpha = 0.8) + scale_y_log10() + labs( title = "Expression of Top 20 Genes", x = "Gene ID", y = "Normalized Counts (log10 scale)", color = "Condition" ) + theme_minimal() + theme( axis.text.x = element_text(angle = 45, hjust = 1), plot.title =element_text(hjust = 0.5) ) The plot highlights the variation in the expression of these genes between the sample groups, and it is easy to observe upregulated or downregulated genes between the conditions. Conclusion and Next Steps These steps show the basic workflow for RNA-seq differential expression analysis using DESeq2, from preprocessing data to identifying significant differentially expressed genes (DEGs). While these results are valuable, they are just the beginning. The next steps involve exploring and validating the findings: Research Genes of Interest: Focusing on the most significant DEGs. Investigating their known functions, roles in relevant pathways, or associations with specific conditions using databases like Ensembl or Gene Ontology. Experimental Validation: Confirming the computational findings through laboratory experiments such as qPCR or western blotting. Functional Analysis: Performing pathway enrichment or gene ontology analysis to uncover biological processes and pathways associated with the significant genes. Bioconductor provides tools like clusterProfiler for this. Hypothesis Generation: Using results to form new hypotheses about the biological mechanisms underlying the condition of interest. By following these steps, RNA-seq data can be turned into meaningful biological insights, driving further research and discovery. References *Oakley, F., Gee, L. M., Sheerin, N. S., & Borthwick, L. A. (2019). Implementation of pre-clinical methodologies to study fibrosis and test anti-fibrotic therapy. Current Opinion in Pharmacology, 49, 95–101. https://doi.org/10.1016/j.coph.2019.10.004  **Oakley, F., Gee, L. M., Sheerin, N. S., & Borthwick, L. A. (2019). Implementation of pre-clinical methodologies to study fibrosis and test anti-fibrotic therapy. Current Opinion in Pharmacology, 49, 95–101. https://doi.org/10.1016/j.coph.2019.10.004 Work was inspired by Datacamp’s RNA-Seq with Bioconductor in R course. "Using Bioconductor for RNA-seq Differential Expression Analysis" Tech Bite was brought to you by Zehra Alibašić, Junior Data Analyst at Atlantbh. (more…)