Key Takeaways

• This end-to-end RNA-Seq analysis workflow runs from raw FASTQ files to interpretable biomarker results - using the nf-core pipeline on AWS HealthOmics and Sonrai Discovery for downstream analysis
• All data is real, publicly available data from a peer-reviewed publication - every step shown here is fully reproducible
• The complete workflow requires no coding and can be completed by a researcher in under 20 minutes on the Sonrai Discovery platform

Jump to section

Objective

The purpose of this document is to describe a workflow detailing the end-to-end RNA-Seq analysis workflow from raw data (FASTQ files) to interpretable results on the Sonrai Discovery platform. This guide will first look at the analysis via the NextFlow RNA-Seq pipeline to process raw data, followed by a no-code exploration and visualization of this data via the Sonrai Discovery Storyboards.

Introduction

Bulk RNA-Seq can be an incredibly powerful tool to help understand gene expression patterns, discover new genes, transcripts and variants, in an incredible variety of cells, tissues and organisms. Despite well established laboratory protocols for sequencing this data, there is often a lack of gold standard for bulk RNA-Seq data analysis [1]. Historically, both academia and industry have had to produce and maintain their own bespoke pipelines which, despite yielding results, may be esoteric, unreproducible when using another pipeline, become unreproducible if versions of software are changed internally, or become outdated with lack of maintenance.

The use of the NextFlow-Core RNA-Seq pipeline [2], and nf-core community pipelines broadly, have been a community-wide effort to tackle the issues stated above by producing easy-to-use, gold-standard pipelines, that are version-controlled, highly-documented, and flexible with their input and analysis depending on the exact outputs required (Fig. 1). In addition, the nf-core pipelines can leverage the power of cloud computing, allowing us to split our workloads across multiple computers, improving the speed at which pipelines run from days on a single computer to hours or minutes using cloud technology.

Our analysis today will leverage the green route shown on the nf-core RNA-Seq pipeline outlined in Fig. 1 - aligning our reads to the most recent version of the human reference genome at the time of writing (the T2T-CHM13v2.0 [3]) using the STAR aligner and Salmon Quantifier [4]. Moreover, Sonrai's nf-core pipeline will leverage Amazon Web Services (AWS) HealthOmics [5], Amazon's platform for hosting pipelines to optimize pipeline execution for both cost and performance.

Figure 1. An adapted version of the nf-core RNA-Seq "tube map" showing the different steps within the pipeline and the routes that can be taken. This demonstration takes the green route (Aligner: STAR, Quantification: Salmon). Gene counts produced by running the pipeline on Sonrai's cloud infrastructure via AWS HealthOmics are loaded into Sonrai Discovery for downstream analysis.

Although the nextflow pipeline produces high-quality outputs containing the raw counts of different genes for different samples, and a detailed multiQC report for understanding the quality of sequencing data, it does not produce immediate and effective visualizations that help extract biological meaning from this data. For a multi-modal approach combining RNA-seq with proteomics and flow cytometry, see our MOFA+ clinical trial case study.

This is where the Sonrai Discovery Platform [6] and its large variety of statistical and bioinformatic tools can be used to ingest and analyze the nf-core RNA-Seq outputs. For this workflow, we will:

1. Understand the inherent quality and structure of the data using principal component analysis (PCA)
2. Perform differential gene expression analysis using a volcano plot to select biomarkers
3. Visualize the pattern of expression via a heatmap to further refine our biomarker search
4. Perform classical statistical analysis via a boxplot and confirm our discovery by revisiting the PCA

This workflow requires no coding or bioinformatician, and can be done by a researcher or non-technical person in less than 20 minutes.

Dataset Analyzed

The data used in this workflow is taken from the paper "Functional Analysis and Transcriptomic Profiling of iPSC-Derived Macrophages and Their Application in Modeling Mendelian Disease" by Zhang et al., 2015 [7].

We will be using 6 samples of macrophages derived from human monocytes (HMDMs). Three samples have been treated with an endotoxin and interferon-gamma for 18-20 hours to induce an inflammatory response (M1), whereas the other 3 have been left untreated (M0). We can use our RNA-Seq pipeline followed by Sonrai Discovery to interrogate the quality of our sequencing data and investigate differences in gene expression between these two treatment groups.

Table 1. Metadata for the 6 samples used in this workflow. The description column (M0-HMDM vs M1-HMDM) provides the binary comparison used throughout. A1BG and A2M are shown as illustrative examples of the bias-corrected count values across 42,000+ genes.

Reproducibility note

All data used in this workflow is publicly available via NCBI SRA (accession numbers SRR1182374-SRR2910664). Every step shown here is fully reproducible - and the same workflow can be applied directly to your proprietary data inside Sonrai's governed research environment.

Quality Control - MultiQC

A key output of the pipeline prior to analysis on Sonrai Discovery is its interactive multiQC report [8], which compiles analysis logs from numerous tools throughout the pipeline with over 40 different metrics from 10 different tools. Despite the large amount of data, the multiQC report also provides overall summaries coalescing these data into summary tables, and will throw warnings or failure messages if a particular metric has not met the threshold, meaning an in-depth knowledge of every plot and table is not needed.

These are fully self-contained, downloadable, customisable prior to running the pipeline and afterwards, and provide the version of each and every tool. Most importantly, they communicate the quality of our data from the sequencing to the final, analysis-ready tables - forming the foundation of a reproducible translational analytics workflow.

Figure 2. Initial opening view of the multiQC report. This summary table shows key sequencing metrics including duplication, reads mapped, GC content, and strand bias. The full report video can be found here.

From an initial scan of the multiQC report, the majority of metrics throw no errors, and the measurements which do are due to particular modules that anticipate different data modalities [9] (e.g. whole genome sequencing as opposed to RNA-Seq) and so will inevitably fail. Given the overwhelming number of metrics that have passed, this RNA-Seq data is amenable to downstream analysis.

Principal Components Analysis (PCA)

Our first port-of-call for a high-dimensional dataset (where each dimension equals a gene or transcript) may be dimensionality reduction, which allows us to detect biologically meaningful clusters and establish outliers in our dataset by collapsing all of our data into 2 or 3 visual dimensions. (Read our article on K-Means clustering for further analysis.) A well-established and commonly used technique, particularly for -omics data is principal components analysis [10].

We can do this easily through our inbuilt PCA functionality by selecting the number of principal components we wish to calculate (up to 100), as well as what we wish to color our samples by. In this case, we are interested in whether the structure of the data reveals differences between our untreated (M0) and inflamed (M1) HMDMs, and so color the points accordingly.

Figure 3. PCA of RNA-Seq expression data. Left: 3D representation of the first 3 principal components, coloured by M0 or M1 macrophage status. Right: Scree plot showing cumulative variation explained per principal component - all variation is captured within the first 5 components. For a video of this app in action: ipsc_pca.webm.

As we can see from the PCA analysis, there is a distinct separation (if not tight clustering) of the control (M0) and treated (M1) groups, indicating inherently different structures in the gene expression data between the control and inflamed macrophages. Although PCA does not provide immediate insight into specific biological differences, it tells us there are underlying differences in our data that can be explored later in our workflow.

Volcano Plot - Differential Gene Expression

Having visualized the underlying structure of this data, we want to try and uncover what genes in particular are significantly differentially expressed. To this end, we can deploy a volcano plot [11] using the python implementation of the DESeq2 algorithm [12] to understand both fold change (how much a gene is up or down regulated), and how significant this change is.

In our plot, we choose which two comparisons we would like to make (in this case, between M0 and M1), and produce our volcano as simply as that. Given the interactive nature of this plot, we can immediately pick out which genes have been both significantly, and most greatly up and downregulated. As we input M1 and M0 as "Contrast variable 1" and "Contrast Variable 2" respectively, we can read the log fold change as "compared to control".

Figure 4. Volcano plot showing the upregulation (green) and downregulation (red) of 42,000+ genes in M0 with respect to M1 via the DESeq2 algorithm. Cut-offs used: log2 fold change of +/-2 (4x up or downregulated) and p-value <0.05. With these filters, 3,084 genes are significantly differentially regulated between M0 and M1. For a video of this app in action: ipsc_volcanoplot_browsing.webm.

Figure 5. Top: The top 5 most differentially expressed genes in M1 (inflamed) macrophages compared to M0 (control). Bottom: Top differentially expressed genes in M0 vs M1. Columns show gene name, mean expression, log2 fold change, standard error, T-statistic, p-value, adjusted p-value (Benjamini-Hochberg), and direction of regulation.

Heatmap

We can begin to break down expressions for the 10 most differentially expressed genes by leveraging a heatmap, which will not only show us the expression for each sample for all 10 genes, but how they cluster together - helping us understand if these genes can be used to explain the difference between the M0 and M1 groups.

Figure 6. Heatmap of the 10 most differentially expressed genes from the volcano plot analysis. The heatmap reveals similarity between expression profiles for M0 and M1 macrophages, clustering between specific genes (e.g. LINC01961 and LINC03025), and variability in expression between samples (e.g. LOC124904143 in M0).

The heatmap reveals several interesting features. Firstly, it shows that these samples will cluster according to the values of these 10 genes, indicating that M0 and M1 macrophage populations can be distinguished from each other based on their expression patterns. Interestingly, LINC02930 and LOC124904143 show high expression in HMDM_MAC_Rep3 (an M0 sample) but not in the other two M0 samples, suggesting it could be anomalous. Conversely, we see a consistent increase in expression across all M1 samples for CSAG3 and LOC101929800, suggesting these may be more robust biomarker candidates for inflammation.

Boxplot

Digging deeper into these 10 differentially expressed genes, we can now perform some classical statistics to compare their variation more directly. We can select the Boxplot tab, click "Add graph" to have two total graphs, and then select the genes of interest (5 most up and down regulated in each respective plot) to compare the raw expression between M0 and M1, as well as understand the variation in expression via the interquartile ranges.

Figure 7. Left: Genes downregulated in M1 macrophages compared to M0. Right: Genes upregulated in M1 compared to M0. The greater variation of upregulated M0 genes (LHS) suggests some upregulation may be partially anomalous. Conversely, larger and more consistent upregulation is seen in M1 macrophages (RHS), indicating significant differences in the expression of these genes.

The left plot shows that the genes upregulated in the uninflamed M0 samples are much more variable, which could indicate that this upregulation is partially by chance, with one sample having an abnormally high expression. Conversely, we see three highly significantly upregulated genes in the M1 macrophages (LINC03025, LOC101929800, CSAG3) and with the exception of LINC03025, we see much tighter distributions for upregulated genes in M1 - notably for LOC101929800 and CSAG3, confirming what we saw in our heatmap and reaffirming our suspicion these are candidate biomarkers for inflammation. The consistency of expression across a treatment group is exactly the kind of signal that underpins biomarker-guided patient stratification.

As we are primarily focused on understanding genes upregulated during the inflammatory response, we might be particularly interested in CSAG3, as it is not only significantly different in M1, but also consistently upregulated. As an additional confirmatory step, we can recolour our initial PCA by the expression of CSAG3.

Figure 8. Comparison between the initial PCA coloured by sample (left) and coloured by CSAG3 expression (right). The M0 and M1 groups show the same separation when coloured by CSAG3 expression as when coloured by their M0/M1 status - strongly suggesting CSAG3 as a candidate biomarker of inflammation.

When we recolor our samples by their CSAG3 expression, we reproduce the distinction between the M0 and M1 groups in the same way coloring by the sample in the first place might do. This strongly suggests that CSAG3 could be a biomarker of inflammation.

Summary

To summarize, we have taken the outputs of the nf-core RNA-Seq analysis workflow and demonstrated a complete path from raw data to candidate biomarker identification. Using the generated multiQC report we confirmed the quality of our sequencing prior to downstream analysis. Through Sonrai Discovery, we used PCA to show the data has patterns distinguishing M0 and M1 macrophages, a volcano plot to extract statistically significant differentially expressed genes, a heatmap to understand clustering and inter-sample variation, and boxplots to confirm the statistical significance and consistency of our top candidates.

This led us to select CSAG3 as a candidate biomarker of inflammation - a finding confirmed by returning to the PCA and recolouring by CSAG3 expression. Lastly, this demonstrates the power of the nf-core RNA-Seq pipeline in combination with the Sonrai Discovery platform to quickly and accurately convert raw data into decision-ready evidence for biomarker discovery.

Further Reading

• AI-Powered Single-Cell Analysis - how single-cell RNA sequencing complements bulk RNA-Seq analysis
• K-Means Clustering - how K-Means clustering can be applied to analyze high-dimensional data
• What is t-SNE? - how t-SNE is employed for high-dimensional data visualization
• What is XGBoost? - how XGBoost is applied in predictive modeling and biomarker discovery
• MOFA+ Clinical Trial Case Study - combining RNA-seq with proteomics and flow cytometry in a multi-modal workflow

References

1. A Beginner's Guide to RNA-Seq Analysis - Novogene webinar. Also: A written guide to RNA-Seq analysis by Matthew N. Bernstein.
2. NextFlow-Core RNAseq pipeline - documentation, parameters, outputs and results.
3. T2T-CHM13v2.0 - the reference genome used for RNA-Seq alignment.
4. STAR aligner and Salmon Quantifier - tools used for aligning and quantifying RNA-Seq reads.
5. AWS HealthOmics - Amazon's platform for hosting pipelines to reduce cost and increase execution time.
6. Sonrai Discovery Platform - the AI-driven platform used to analyze the raw RNA-Seq data.
7. Zhang et al., 2015 - "Functional Analysis and Transcriptomic Profiling of iPSC-Derived Macrophages and Their Application in Modeling Mendelian Disease."
8. MultiQC - reporting tool for aggregating results from multiple tools and samples.
9. Bioinformatics StackExchange - why particular modules may fail in multiQC with different data modalities.
10. Principal Components Analysis - an introduction to the dimensionality reduction technique.
11. Volcano plots - an explanation of volcano plots for visualizing differential gene expression.
12. DESeq2 algorithm - an explanation of the DESeq2 algorithm used to calculate differential expression.

No cost. No risk.

Want to run this workflow on your own data?

Challenge Sonrai to find what others haven't. We'll apply our RNA-Seq analysis workflow to your dataset - at no cost.

Explore the free showcase

Use Case

AI-Powered Single-Cell Analysis

Use Case

AI Biomarker Discovery Workflow

Case Study

MOFA+ Clinical Trial Case Study

Talk to our team

Sonrai acts as an extension of your team, combining scientific expertise, transparent and reproducible analysis, and a governed research environment to generate decision-ready evidence that can be revisited, reproduced, and built upon throughout the life of a programme.

Book a demo