# How to Use ggplot2 for RNA-seq Visualization: A Practical Guide for Custom Publication-Ready Plots
## Key Takeaways
- **Data Reshaping is Foundational:** RNA-seq data, typically in a wide format (genes as rows, samples as columns), must be transformed into a long format using functions like `pivot_longer` from `tidyr` for effective visualization with `ggplot2`. This long format enables `ggplot2` to map variables like gene expression to aesthetics such as color or size, facilitating group comparisons and faceting.
- **Grammar of Graphics for Control:** `ggplot2`'s layered approach, based on the grammar of graphics, allows incremental construction of complex plots. This enables precise control over elements like aesthetic mappings (e.g., `color = condition`), geometric objects (`geom_point` for PCA, `geom_tile` for heatmaps), scales, and themes, crucial for tailoring figures to publication standards.
- **PCA for Global Variation and QC:** Principal Component Analysis (PCA) on normalized and variance-stabilized counts is essential for visualizing global sample variation and identifying potential batch effects or outliers. Plotting PC1 vs. PC2 with color-coded sample groups and confidence ellipses, as demonstrated with `stat_ellipse`, is a standard method for assessing sample relationships.
- **Heatmaps for Expression Patterns:** Heatmaps, constructed using `geom_tile` with scaled expression values, are vital for visualizing co-expression patterns across selected genes and samples. Hierarchical clustering of genes and samples, achieved by reordering factor levels based on `hclust` output, reveals biologically meaningful clusters and relationships.
- **Volcano and MA Plots for Differential Expression:** Volcano plots (log2 fold change vs. -log10 adjusted p-value) and MA plots (log2 mean expression vs. log2 fold change) are critical for visualizing differential gene expression results. Threshold lines and color-coding by significance status (e.g., upregulated, downregulated) highlight key genes, while `ggrepel` aids in labeling specific gene symbols.
- **Custom Themes and Reproducibility:** Publication-ready figures require consistent styling, achievable through custom `ggplot2` theme functions (e.g., `theme_publication`) that control fonts, grid lines, and panel appearance. Integrating visualization scripts into reproducible workflows using version control (Git) and workflow managers (Snakemake, Nextflow) ensures that figures can be regenerated and validated.
---
RNA sequencing produces large, multidimensional datasets that require careful visualization for interpretation and publication. The R package ggplot2 provides a flexible, layered system for creating custom plots from RNA-seq data, but many researchers struggle to move beyond default templates. This guide addresses the specific problem of building publication-ready visualizations from RNA-seq outputs using a structured grammar-of-graphics approach. You will learn how to reshape raw count data, construct principal component analysis (PCA) plots, generate heatmaps, customize themes for journal submission, and implement quality control visualizations that support reproducible analysis workflows.
## Understanding the Role of Visualization in RNA-seq Analysis
RNA-seq experiments generate count matrices that capture transcript abundance across thousands of genes and multiple biological conditions. These data require multiple visualization layers to support different analytical goals. Quality control plots reveal sample outliers and batch effects before downstream analysis. Dimensionality reduction plots such as PCA summarize global variation across samples. Differential expression visualizations such as volcano plots and MA plots highlight genes that change significantly between conditions. Heatmaps display expression patterns across gene sets or sample groups. Each plot type answers a distinct biological question and requires specific data structures and ggplot2 mappings.
The grammar of graphics, which underlies ggplot2, treats every plot as a combination of data, aesthetic mappings, geometric objects, scales, coordinate systems, and themes. This layered approach allows researchers to build complex figures incrementally and maintain full control over every visual element. For RNA-seq data, this means you can start with a basic scatter plot for PCA scores and progressively add group colors, confidence ellipses, axis labels, and publication themes without rewriting the entire plotting code.
Published RNA-seq studies commonly use ggplot2 for visualization across diverse biological contexts. A study on ADAM19 in systemic sclerosis used R packages including edgeR, limma, clusterProfiler, ggplot2, gseaplot2, and complexheatmap for data analysis and visualization of bulk RNA-seq data [1]. Similarly, research on PROS1 in glioma performed gene ontology enrichment and KEGG pathway analyses using the clusterProfiler package and visualized results with ggplot2 [2]. A preeclampsia study used the Limma package for differential expression analysis and illustrated gene set enrichment analysis results with ClusterProfiler and ggplot2 [3]. A cataract study visualized differentially expressed gene enrichment results using the R ggplot2 package [4]. These examples demonstrate that ggplot2 serves as a standard visualization layer within broader RNA-seq analysis pipelines.
## Core Principles of the Grammar of Graphics for Genomic Data
The grammar of graphics organizes plot construction around seven fundamental components. Data specifies the dataset to be plotted. Aesthetics map variables to visual properties such as x and y coordinates, color, size, and shape. Geometries define the type of plot, including points, lines, bars, and tiles. Scales control how data values map to visual properties, including axis limits, color gradients, and legend breaks. Coordinate systems determine the spatial arrangement of the plot. Faceting splits data into multiple panels. Themes control non-data elements such as fonts, grid lines, and background colors.
For RNA-seq visualization, the most frequently used geometries include `geom_point` for PCA scores and volcano plots, `geom_tile` for heatmaps, `geom_boxplot` for expression distributions, `geom_histogram` for count distributions, and `geom_line` for trajectory or enrichment plots. Aesthetic mappings typically include x and y coordinates for sample or gene positions, color for experimental groups, and fill for density or heatmap values.
The key to effective RNA-seq visualization is understanding how to transform raw analytical outputs into the tidy data format that ggplot2 expects. Most RNA-seq analysis packages produce matrices where genes are rows and samples are columns. ggplot2 requires data in long format where each row represents one observation. The `tidyr` package provides the `pivot_longer` function to reshape wide count matrices into long format. This reshaping step is essential before any ggplot2 plotting can occur.
Consider a typical differential expression output from DESeq2 or edgeR. The result table contains columns for gene identifiers, log2 fold changes, p-values, and adjusted p-values. This table is already in a suitable format for volcano plots and MA plots because each row represents one gene with all relevant statistics. However, for plotting expression values across samples, you must reshape the count matrix to long format and join it with metadata.
## Preparing RNA-seq Data for ggplot2
### Required Input Data Structures
RNA-seq visualization requires three primary data structures. The count matrix contains raw or normalized read counts with genes as rows and samples as columns. The sample metadata table contains experimental information such as condition, treatment, time point, or batch for each sample. The differential expression results table contains gene-level statistics including log2 fold changes and adjusted p-values.
Count matrices typically come from quantification tools such as RSEM, which provides gene and isoform abundance estimates from single-end or paired-end RNA-seq data [5]. The output includes estimated counts, transcript per million values, and 95% credibility intervals. For ggplot2 visualization, you will usually work with normalized counts instead of raw counts to account for sequencing depth differences across samples.
Sample metadata should include unique sample identifiers that match the column names in the count matrix. This matching is critical for joining expression data with experimental conditions. Common metadata columns include sample name, biological condition, replicate number, sequencing batch, and any covariates such as sex or age.
### Data Reshaping with tidyr
The first practical step is converting the wide count matrix into long format. The `pivot_longer` function transforms the matrix so that each row contains a gene-sample pair with its expression value. This long format enables grouping and faceting by sample or condition in ggplot2.
```r
library(tidyr)
library(dplyr)
## Assume counts_matrix has genes as rows and samples as columns
## Assume sample_metadata has columns: sample_id, condition
counts_long <- counts_matrix %>%
rownames_to_column(var = "gene") %>%
pivot_longer(cols = -gene,
names_to = "sample_id",
values_to = "expression") %>%
left_join(sample_metadata, by = "sample_id")
```
This reshaped data frame now supports plotting expression distributions by condition using `geom_boxplot` or `geom_violin`. It also enables faceting by gene to compare expression patterns across multiple genes of interest.
### Joining Differential Expression Results with Annotations
For volcano plots and MA plots, you need a data frame that combines differential expression statistics with gene annotations. The `left_join` function merges the results table with a gene annotation table containing gene symbols, biotypes, or chromosome locations.
```r
## Assume de_results has columns: gene_id, log2FoldChange, pvalue, padj
## Assume gene_annotations has columns: gene_id, gene_symbol, description
de_annotated <- de_results %>%
left_join(gene_annotations, by = "gene_id") %>%
mutate(regulation = case_when(
padj < 0.05 & log2FoldChange > 1 ~ "upregulated",
padj < 0.05 & log2FoldChange < -1 ~ "downregulated",
TRUE ~ "not significant"
))
```
The `regulation` column enables coloring points by significance and direction in volcano plots. This annotation step also supports labeling specific genes of interest using `geom_text` or `ggrepel`.
## Constructing Principal Component Analysis Plots
### Computing PCA from Normalized Counts
PCA reduces the dimensionality of the expression matrix to capture the major sources of variation across samples. The `prcomp` function in base R performs PCA on a matrix where rows are genes and columns are samples. Before running PCA, you should filter low-expression genes and apply a variance-stabilizing transformation or regularized log transformation to normalize the data.
```r
## Assume normalized_counts has genes as rows and samples as columns
## Filter genes with low expression
keep_genes <- rowSums(normalized_counts > 10) >= 3
filtered_counts <- normalized_counts[keep_genes, ]
## Transpose for PCA (samples as rows, genes as columns)
pca_result <- prcomp(t(filtered_counts), scale. = TRUE)
## Extract variance explained
percent_var <- summary(pca_result)$importance[6] * 100
## Create data frame for plotting
pca_scores <- as.data.frame(pca_result$x) %>%
rownames_to_column(var = "sample_id") %>%
left_join(sample_metadata, by = "sample_id")
```
The `scale. = TRUE` argument standardizes each gene to have unit variance, which prevents genes with high expression from dominating the PCA. This is particularly important for RNA-seq data where expression values span several orders of magnitude.
### Building a Publication-Ready PCA Plot
A publication-ready PCA plot includes principal component scores as points, color coding for experimental groups, confidence ellipses, and axis labels with variance percentages.
```r
library(ggplot2)
pca_plot <- ggplot(pca_scores, aes(x = PC1, y = PC2, color = condition)) +
geom_point(size = 3, alpha = 0.8) +
stat_ellipse(aes(fill = condition), geom = "polygon", alpha = 0.2) +
labs(x = paste0("PC1: ", round(percent_var[7], 1), "% variance"),
y = paste0("PC2: ", round(percent_var[6], 1), "% variance"),
color = "Condition", fill = "Condition") +
theme_bw(base_size = 14) +
theme(legend.position = "right")
```
The `stat_ellipse` function adds confidence ellipses around each group, which helps visualize group separation. The axis labels include the percentage of variance explained by each principal component, a standard requirement for publication figures.
### Adding Sample Labels and Customizing Point Shapes
For small sample sizes, adding sample labels directly to the PCA plot improves readability. The `geom_text` function with `hjust` and `vjust` adjustments places labels near each point.
```r
pca_plot_labeled <- pca_plot +
geom_text(aes(label = sample_id), hjust = -0.3, vjust = 0.5, size = 3)
```
When samples come from multiple batches or conditions, you can map different aesthetics to different metadata variables. For example, color can represent biological condition while shape represents sequencing batch. This dual mapping helps identify batch effects that might confound biological variation.
## Creating Heatmaps for Gene Expression Patterns
### Preparing Data for Heatmap Visualization
Heatmaps display expression values across genes and samples using color intensity. While the `pheatmap` and `ComplexHeatmap` packages provide dedicated heatmap functions, ggplot2 offers flexibility for custom heatmap designs using `geom_tile`.
The data must be in long format with columns for gene, sample, and expression value. For publication, you typically display a subset of genes such as differentially expressed genes, pathway members, or marker genes.
```r
## Assume de_genes is a vector of gene identifiers of interest
## Assume normalized_counts has genes as rows and samples as columns
heatmap_data <- normalized_counts[de_genes, ] %>%
rownames_to_column(var = "gene") %>%
pivot_longer(cols = -gene,
names_to = "sample_id",
values_to = "expression") %>%
left_join(sample_metadata, by = "sample_id")
```
For heatmaps, you should scale expression values per gene so that the color scale is comparable across genes with different baseline expression levels. The `scale` function in base R centers and scales each gene row.
```r
scaled_matrix <- t(scale(t(normalized_counts[de_genes, ])))
```
### Building the Heatmap with geom_tile
The `geom_tile` geometry creates a grid of colored rectangles. The x aesthetic maps to samples, the y aesthetic maps to genes, and the fill aesthetic maps to scaled expression values.
```r
heatmap_plot <- ggplot(heatmap_data, aes(x = sample_id, y = gene, fill = scaled_expression)) +
geom_tile() +
scale_fill_gradient2(low = "blue", mid = "white", high = "red",
midpoint = 0, name = "Z-score") +
theme_minimal(base_size = 12) +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
axis.text.y = element_text(size = 8))
```
The `scale_fill_gradient2` function creates a diverging color scale with blue for low expression, white for medium expression, and red for high expression. This color scheme is standard for heatmaps because it highlights both upregulation and downregulation relative to the mean.
### Ordering Genes and Samples for Pattern Recognition
Heatmap interpretation improves when genes and samples are ordered by similarity instead of alphabetical order. Hierarchical clustering provides this ordering. The `hclust` function computes a dendrogram, and the `order` function extracts the leaf order.
```r
## Cluster genes
gene_dist <- dist(scaled_matrix)
gene_clust <- hclust(gene_dist)
gene_order <- gene_clust$labels[gene_clust$order]
## Cluster samples
sample_dist <- dist(t(scaled_matrix))
sample_clust <- hclust(sample_dist)
sample_order <- sample_clust$labels[sample_clust$order]
## Reorder factor levels for plotting
heatmap_data$gene <- factor(heatmap_data$gene, levels = gene_order)
heatmap_data$sample_id <- factor(heatmap_data$sample_id, levels = sample_order)
```
Setting factor levels based on clustering order ensures that `geom_tile` displays the heatmap with clustered rows and columns. This ordering reveals blocks of co-expressed genes and sample groups with similar expression profiles.
## Customizing Themes for Publication Standards
### Understanding ggplot2 Theme Components
The theme system controls all non-data elements of a plot, including text, lines, rectangles, and legend appearance. Journal-specific formatting often requires specific font sizes, font families, and panel backgrounds. The `theme` function modifies individual elements, while complete themes such as `theme_bw`, `theme_minimal`, and `theme_classic` provide starting points.
Common publication requirements include sans-serif fonts such as Arial or Helvetica, font sizes between 8 and 12 points for figure text, and minimal grid lines. Journals often specify exact dimensions in inches or centimeters for single-column or double-column figures.
### Creating a Reusable Publication Theme
Defining a custom theme function ensures consistency across all figures in a manuscript. This function can be saved in an R script and sourced at the beginning of each analysis session.
```r
theme_publication <- function(base_size = 10, base_family = "Arial") {
theme_bw(base_size = base_size, base_family = base_family) +
theme(
plot.title = element_text(face = "bold", size = rel(1.2), hjust = 0.5),
axis.title = element_text(face = "bold"),
axis.text = element_text(color = "black"),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "black", fill = NA, linewidth = 0.8),
legend.title = element_text(face = "bold"),
legend.background = element_rect(fill = "white", color = "black", linewidth = 0.3),
strip.background = element_rect(fill = "grey90", color = "black"),
strip.text = element_text(face = "bold")
)
}
```
This theme removes minor grid lines, sets black axis text, adds a black panel border, and formats legend and facet titles. Applying this theme to all plots creates a cohesive visual style across the manuscript.
### Adjusting Fonts and Exporting High-Resolution Figures
For journal submission, figures must be exported at sufficient resolution. The `ggsave` function exports plots to various formats including PDF, TIFF, and PNG. The `dpi` argument controls resolution for raster formats, and `width` and `height` control physical dimensions.
```r
ggsave("pca_plot.pdf", pca_plot + theme_publication(),
width = 6, height = 4, units = "in", dpi = 300)
```
For journals that require specific font embedding, the `cairo_pdf` device provides better font support than the default PDF device. The `showtext` package enables custom font families that may not be installed on the system.
## Visualizing Differential Expression Results
### Volcano Plots for Identifying Significant Genes
Volcano plots display log2 fold change on the x-axis and negative log10 adjusted p-value on the y-axis. Each point represents one gene. Genes with large fold changes and significant p-values appear in the upper left and upper right corners.
```r
volcano_plot <- ggplot(de_annotated, aes(x = log2FoldChange, y = -log10(padj))) +
geom_point(aes(color = regulation), size = 1.5, alpha = 0.6) +
scale_color_manual(values = c("upregulated" = "red",
"downregulated" = "blue",
"not significant" = "grey70")) +
geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "grey40") +
geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "grey40") +
labs(x = "Log2 fold change", y = "-Log10 adjusted p-value") +
theme_publication()
```
The threshold lines at adjusted p-value of 0.05 and log2 fold change of plus or minus 1 provide visual reference for significance criteria. These thresholds should match the criteria used in the differential expression analysis.
### Labeling Key Genes in Volcano Plots
For manuscripts, you often need to label specific genes of interest. The `ggrepel` package provides labels that avoid overlapping points.
```r
library(ggrepel)
## Identify genes to label
genes_to_label <- de_annotated %>%
filter(padj < 0.01 & abs(log2FoldChange) > 2) %>%
slice_max(order_by = abs(log2FoldChange), n = 20)
volcano_plot_labeled <- volcano_plot +
geom_text_repel(data = genes_to_label,
aes(label = gene_symbol),
size = 3, max.overlaps = 20)
```
The `geom_text_repel` function automatically adjusts label positions to minimize overlap, which is essential when many genes meet the labeling criteria.
### MA Plots for Assessing Expression-Dependent Changes
MA plots display the log2 mean expression on the x-axis and the log2 fold change on the y-axis. These plots reveal whether fold changes depend on expression level, which can indicate normalization issues.
```r
ma_plot <- ggplot(de_annotated, aes(x = log2(baseMean), y = log2FoldChange)) +
geom_point(aes(color = regulation), size = 1.5, alpha = 0.5) +
geom_hline(yintercept = 0, linetype = "solid", color = "black") +
scale_color_manual(values = c("upregulated" = "red",
"downregulated" = "blue",
"not significant" = "grey70")) +
labs(x = "Log2 mean expression", y = "Log2 fold change") +
theme_publication()
```
The `baseMean` column from DESeq2 results provides the mean normalized count across all samples. MA plots help identify systematic biases where low-expression genes show inflated fold changes.
## Quality Control Visualizations for RNA-seq Data
### Assessing Sample Similarity with Correlation Heatmaps
Before downstream analysis, you should verify that biological replicates cluster together and that no sample is an outlier. A sample-to-sample correlation heatmap provides this assessment.
```r
## Compute correlation matrix
cor_matrix <- cor(normalized_counts, method = "spearman")
## Reshape for ggplot2
cor_long <- cor_matrix %>%
as.data.frame() %>%
rownames_to_column(var = "sample1") %>%
pivot_longer(cols = -sample1, names_to = "sample2", values_to = "correlation")
cor_heatmap <- ggplot(cor_long, aes(x = sample1, y = sample2, fill = correlation)) +
geom_tile() +
scale_fill_gradient2(low = "blue", mid = "white", high = "red",
midpoint = 0.9, limits = c(0.8, 1)) +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
```
Samples with low correlation to all other samples may indicate failed library preparation, sample contamination, or mislabeling. These samples should be investigated before proceeding with differential expression analysis.
### Visualizing Count Distributions Across Samples
Boxplots of log-transformed counts reveal differences in sequencing depth and distribution shape across samples. Samples with substantially different distributions may require normalization or exclusion.
```r
count_distribution <- ggplot(counts_long, aes(x = sample_id, y = log2(expression + 1))) +
geom_boxplot(aes(fill = condition), outlier.size = 0.5) +
labs(x = "Sample", y = "Log2(count + 1)") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
```
The `+1` offset avoids taking the logarithm of zero. This visualization helps identify samples with unusually low total counts or abnormal distribution shapes.
### Detecting Batch Effects with PCA Before and After Correction
PCA plots generated before and after batch correction provide visual evidence of batch effect removal. If samples cluster by batch instead of biological condition in the uncorrected PCA, batch correction is necessary.
```r
## Before correction
pca_before <- plot_pca(normalized_counts, sample_metadata, color_by = "batch")
## After correction using ComBat-seq or similar
corrected_counts <- apply_batch_correction(normalized_counts, batch = sample_metadata$batch)
pca_after <- plot_pca(corrected_counts, sample_metadata, color_by = "batch")
```
Comparing these two plots helps determine whether batch correction successfully removed technical variation while preserving biological signal. The Crescendo method provides batch correction specifically for single-cell spatial transcriptomics count data, improving visualization and detection of spatial gene patterns [6].
## Visualizing Gene Set Enrichment Results
### Preparing Enrichment Results for Plotting
Gene set enrichment analysis produces results tables with pathway names, enrichment scores, p-values, and gene counts. The clusterProfiler package generates these results, and ggplot2 visualizes them effectively [2][3].
```r
## Assume enrich_results has columns: Description, GeneRatio, p.adjust, Count
## Assume enrich_results is ordered by p.adjust
top_pathways <- enrich_results %>%
slice_min(order_by = p.adjust, n = 20) %>%
mutate(Description = factor(Description, levels = rev(Description)))
```
The factor reordering ensures that pathways display in order of significance when plotted.
### Dot Plots for Enrichment Results
Dot plots display pathways on the y-axis, gene ratio on the x-axis, point size for gene count, and point color for adjusted p-value.
```r
dot_plot <- ggplot(top_pathways, aes(x = GeneRatio, y = Description)) +
geom_point(aes(size = Count, color = p.adjust)) +
scale_color_gradient(low = "red", high = "blue", name = "Adjusted p-value") +
scale_size_continuous(name = "Gene count") +
labs(x = "Gene ratio", y = "") +
theme_publication() +
theme(axis.text.y = element_text(size = 8))
```
This visualization communicates both the significance and the magnitude of enrichment for each pathway. The gene ratio represents the proportion of genes in the gene set that are present in the pathway.
### Bar Plots for Enrichment Scores
For gene set enrichment analysis using the GSEA algorithm, bar plots display normalized enrichment scores for each pathway.
```r
## Assume gsea_results has columns: Description, NES, p.adjust
gsea_bar <- ggplot(gsea_results, aes(x = reorder(Description, NES), y = NES, fill = p.adjust)) +
geom_bar(stat = "identity") +
coord_flip() +
scale_fill_gradient(low = "red", high = "blue", name = "Adjusted p-value") +
labs(x = "", y = "Normalized enrichment score") +
theme_publication()
```
The `coord_flip` function rotates the plot so pathway names display horizontally, which improves readability for long pathway names.
## Visualizing Single-Cell RNA-seq Data with ggplot2
### Adapting ggplot2 for Single-Cell Visualizations
Single-cell RNA sequencing produces data at the individual cell level, requiring specialized visualization approaches. Standard workflows generate UMAP or t-SNE embeddings where each cell is a point in two-dimensional space. The Seurat package provides built-in plotting functions, but ggplot2 offers greater customization control [7].
```r
## Assume seurat_obj is a Seurat object with UMAP coordinates
## Extract cell metadata and embeddings
umap_data <- data.frame(
UMAP1 = seurat_obj@reductions$umap@cell.embeddings[7],
UMAP2 = seurat_obj@reductions$umap@cell.embeddings[6],
cell_type = seurat_obj$cell_type,
condition = seurat_obj$condition
)
umap_plot <- ggplot(umap_data, aes(x = UMAP1, y = UMAP2, color = cell_type)) +
geom_point(size = 0.5, alpha = 0.6) +
theme_publication() +
theme(legend.position = "right")
```
The SCpubr package provides concise function calls for generating publication-ready visualizations commonly used in single-cell transcriptome analyses, building on the working knowledge of Seurat and ggplot2 that many experimental biologists already possess [7].
### Feature Plots for Gene Expression in Single-Cell Data
Feature plots display the expression level of a specific gene across all cells in the embedding. These plots help identify cell populations that express particular marker genes.
```r
## Extract gene expression for a specific gene
gene_expression <- GetAssayData(seurat_obj, assay = "RNA", layer = "data")["GENE_NAME", ]
umap_data$gene_expression <- gene_expression
feature_plot <- ggplot(umap_data, aes(x = UMAP1, y = UMAP2, color = gene_expression)) +
geom_point(size = 0.5, alpha = 0.8) +
scale_color_gradient(low = "grey90", high = "red", name = "Expression") +
theme_publication()
```
The color gradient from grey to red highlights cells with high expression while keeping low-expressing cells visually subtle.
### Visualizing Cell-Cell Communication Results
Cell-cell communication analysis produces results that describe ligand-receptor interactions between cell types. These results can be visualized as dot plots or heatmaps using ggplot2.
```r
## Assume communication_results has columns: source, target, ligand, receptor, score
communication_plot <- ggplot(communication_results,
aes(x = source, y = target, size = score)) +
geom_point(aes(color = score)) +
scale_size_continuous(range = c(1, 10)) +
scale_color_gradient(low = "blue", high = "red") +
labs(x = "Source cell type", y = "Target cell type") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
```
This visualization reveals which cell types communicate most strongly and which ligand-receptor pairs mediate the interactions.
## Working with Genome Coverage and Annotation Tracks
### Visualizing Coverage Data with ggplot2
Genome coverage visualization is important for inspecting and interpreting next-generation sequencing data, including RNA-seq alignments [8]. The ggcoverage package provides an R-based approach to visualize and annotate genome coverage from BAM, BigWig, BedGraph, and TSV formats [8]. While ggcoverage handles the specialized track layouts, ggplot2 can create custom coverage plots for specific genomic regions.
```r
## Assume coverage_data has columns: position, sample, coverage
coverage_plot <- ggplot(coverage_data, aes(x = position, y = coverage, color = sample)) +
geom_line(linewidth = 0.8) +
labs(x = "Genomic position", y = "Read coverage") +
theme_publication()
```
Coverage plots help verify that alignments are distributed as expected across exons and that no regions show abnormal pileup.
### Adding Annotation Tracks
Genome annotations such as gene models, exons, and transcription start sites provide context for coverage plots. These annotations can be added as rectangles or lines using ggplot2 geometries.
```r
## Assume annotation_data has columns: start, end, feature_type, gene_name
annotation_plot <- ggplot(annotation_data, aes(xmin = start, xmax = end, y = feature_type)) +
geom_rect(aes(fill = gene_name), alpha = 0.5) +
labs(x = "Genomic position", y = "") +
theme_publication()
```
Combining coverage and annotation tracks in a single figure requires aligning the x-axis scales. The `facet_grid` function can stack coverage and annotation panels with shared x-axis coordinates.
## Reproducibility and Workflow Integration
### Structuring Analysis Scripts for Reproducibility
Reproducible visualization requires structured scripts that document every step from raw data to final figure. The Carpentries lessons provide foundational training in computing, data handling, shell, Git, and programming that supports reproducible analysis workflows [9]. Key practices include using relative file paths, setting random seeds, and recording package versions.
```r
## Record session information
sessionInfo()
## Set random seed for reproducible results
set.seed(42)
```
The `sessionInfo` function records R version, platform, and package versions. Including this output in analysis documentation ensures that others can recreate the analysis environment.
### Using Workflow Managers for Pipeline Integration
Workflow managers such as Snakemake and Nextflow automate the execution of analysis pipelines, ensuring that each step runs in the correct order with proper inputs and outputs. The nf-core documentation describes community pipeline standards, usage, configuration, and reproducible workflow context [10]. The RAGER platform integrates popular bioinformatics tools in an automated thread for joint mining of RNA-seq and ATAC-seq data, implemented using Snakemake [11].
For visualization, workflow integration means that plotting scripts run automatically after differential expression analysis completes. The workflow manager passes result files to the plotting script and collects output figures.
### Version Control for Plotting Code
Version control with Git tracks changes to analysis scripts and enables collaboration. The Carpentries lessons include Git training that covers committing changes, branching, and merging [9]. For visualization code, version control ensures that figure generation is reproducible and that changes are documented.
```r
## Example Git workflow
## git init
## git add analysis_scripts/
## git commit -m "Add PCA and heatmap visualization scripts"
```
Each commit should correspond to a logical unit of work, such as adding a new plot type or fixing a visualization bug.
## Common Failure Patterns in ggplot2 RNA-seq Visualization
### Factor Level Ordering Issues
A frequent problem occurs when ggplot2 orders categorical variables alphabetically instead of in the desired biological order. This affects the order of samples on the x-axis of heatmaps, the order of pathways in dot plots, and the order of groups in boxplots.
The solution is to explicitly set factor levels using the `factor` function with the `levels` argument. This should be done after data reshaping and before plotting.
```r
sample_metadata$condition <- factor(sample_metadata$condition,
levels = c("control", "treatment"))
```
### Overlapping Points in Dense Plots
Volcano plots and single-cell UMAP plots often contain thousands of points that overlap, obscuring data density. Solutions include reducing point size, increasing transparency with the `alpha` argument, and using density-based visualization.
```r
## Use transparency to reveal overlapping points
geom_point(size = 1, alpha = 0.3)
## Use hexagonal binning for dense regions
geom_hex(bins = 50)
```
The `geom_hex` geometry divides the plot area into hexagonal bins and colors each bin by point count, which effectively visualizes density in large datasets.
### Legend and Axis Label Formatting Problems
Journals often require specific legend titles, axis labels, and font sizes. Default ggplot2 output uses variable names as labels, which are rarely publication-ready. The `labs` function overrides all labels in one call.
```r
labs(x = "Log2 fold change",
y = "-Log10 adjusted p-value",
color = "Regulation status",
title = "Differential expression in treated vs control")
```
### Color Blindness and Accessibility Issues
Default ggplot2 color scales may not be distinguishable for color-blind readers. The `viridis` package provides color scales that are perceptually uniform and accessible.
```r
scale_color_viridis_d() # for discrete variables
scale_fill_viridis_c() # for continuous variables
```
These scales maintain color distinction across the full range of values and are increasingly required by journals.
## Limitations and Interpretation Boundaries
### Statistical Limitations of Visualization
Visualizations reveal patterns but do not establish statistical significance. A PCA plot showing group separation does not prove that the separation is statistically meaningful. Permutation tests or multivariate statistical methods provide formal significance testing.
Similarly, volcano plots display adjusted p-values and fold changes but do not account for biological relevance. A gene with a large fold change and significant p-value may still have low biological importance if its absolute expression is negligible.
### Technical Limitations of ggplot2 for Large Datasets
ggplot2 can become slow or memory-intensive with very large datasets. Single-cell datasets with hundreds of thousands of cells may require downsampling or rasterization for efficient plotting. The `geom_point` function with `raster = TRUE` renders points as a raster layer, reducing file size and rendering time.
```r
geom_point(size = 0.5, alpha = 0.5, raster = TRUE)
```
For genome-wide coverage data, plotting every position is impractical. Aggregating coverage into bins or plotting only regions of interest reduces the data volume.
### Interpretation Caveats for Normalized Data
Visualizations based on normalized counts depend on the normalization method used. Different normalization approaches can produce different visual patterns, particularly for PCA and heatmaps. Researchers should document the normalization method and consider how it affects interpretation.
The choice between variance-stabilizing transformation, regularized log transformation, and TPM values affects the visual appearance of plots. Each method has different properties regarding the handling of low-count genes and library size differences.
## Professional Escalation Criteria
### When to Seek Bioinformatics Support
Researchers should escalate to bioinformatics specialists when visualization reveals patterns that suggest technical artifacts. These include samples that cluster by batch instead of condition, correlation heatmaps showing one sample with uniformly low correlation, and PCA plots where technical replicates do not cluster together.
Bioinformatics support is also warranted when the analysis requires methods beyond standard workflows, such as integrating multiple data types, applying advanced batch correction, or implementing custom statistical models.
### When to Consult Statistics Experts
Statistical consultation is appropriate when differential expression results show unexpected patterns, such as an unusually high number of significant genes, fold changes that correlate strongly with expression level, or inconsistent results across analysis methods.
A study on PCOS with insulin resistance identified 339 common differentially expressed genes using both DESeq2 and Limma, demonstrating that different tools can produce overlapping but not identical results [12]. When tools disagree substantially, statistical expertise helps determine the source of discrepancy.
### When to Involve Domain Experts
Biological interpretation of visualization results benefits from domain expertise. If enrichment analysis reveals unexpected pathways, or if marker genes do not align with known biology, consultation with a domain expert helps determine whether the results reflect genuine biology or technical artifacts.
The ADAM19 study in systemic sclerosis used publicly available transcriptome datasets to assess gene expression, then validated findings through real-time PCR, western blot, and immunostaining [1]. This validation approach demonstrates the importance of confirming computational findings with independent experimental methods.
## Records and Documentation for Visualization Workflows
### Maintaining Analysis Logs
Documentation of visualization workflows should include the input data files, R scripts, package versions, and output figures. This documentation enables others to reproduce the figures and understand the analytical decisions made.
```r
## Example analysis log entry
## Date: 2025-01-15
## Input: counts_matrix.csv, sample_metadata.csv
## Script: visualization_workflow.R
## Packages: ggplot2 3.4.4, tidyr 1.3.0, dplyr 1.1.3
## Output: pca_plot.pdf, heatmap_plot.pdf, volcano_plot.pdf
```
### Storing Figure Generation Code
Each figure should have a corresponding R script that generates it from the analysis results. These scripts should be stored in a version-controlled directory structure.
```text
figures/
figure1_pca.R
figure2_heatmap.R
figure3_volcano.R
figure4_enrichment.R
```
This structure ensures that figures can be regenerated when analysis parameters change or when journals request format modifications.
### Recording Parameter Choices
Visualization parameters such as color scales, point sizes, and threshold values should be documented. These choices affect interpretation and must be consistent across figures in a manuscript.
```r
## Document threshold choices
significance_threshold <- 0.05
fold_change_threshold <- 1
```
Recording these parameters in the analysis script makes the visualization workflow transparent and reproducible.
## Practical Implementation Steps
### Step 1: Install Required Packages
Start by installing the packages needed for RNA-seq visualization. The Bioconductor project provides official package, workflow, installation, and reproducible genomic-analysis documentation [13].
```r
install.packages(c("ggplot2", "tidyr", "dplyr", "ggrepel", "viridis"))
BiocManager::install(c("DESeq2", "edgeR", "clusterProfiler", "ComplexHeatmap"))
```
### Step 2: Load and Inspect Data
Load the count matrix, sample metadata, and differential expression results. Inspect the structure of each data frame to confirm that sample identifiers match across files.
```r
counts <- read.csv("counts_matrix.csv", row.names = 1)
metadata <- read.csv("sample_metadata.csv")
de_results <- read.csv("de_results.csv")
```
### Step 3: Reshape Data for ggplot2
Convert the count matrix to long format and join with metadata. This step is required for most ggplot2 visualizations.
```r
counts_long <- counts %>%
rownames_to_column(var = "gene") %>%
pivot_longer(cols = -gene, names_to = "sample_id", values_to = "count") %>%
left_join(metadata, by = "sample_id")
```
### Step 4: Generate Quality Control Plots
Create boxplots of count distributions, correlation heatmaps, and PCA plots to assess data quality before downstream analysis.
### Step 5: Create Publication Figures
Generate PCA plots, heatmaps, volcano plots, and enrichment visualizations using the custom publication theme.
### Step 6: Export Figures
Export each figure in the format and resolution required by the target journal.
## At a Glance: ggplot2 Visualization Workflow for RNA-seq
| Analysis Stage | Recommended Plot Type | Key ggplot2 Geometry | Primary Decision |
| --- | --- | --- | --- |
| Quality control | Sample correlation heatmap | `geom_tile` | Identify outlier samples or batch effects |
| Global variation | PCA score plot | `geom_point` with `stat_ellipse` | Assess group separation and clustering |
| Differential expression | Volcano plot | `geom_point` with color mapping | Identify significant up and downregulated genes |
| Expression patterns | Heatmap of scaled counts | `geom_tile` with `scale_fill_gradient2` | Visualize co-expression clusters |
| Pathway enrichment | Dot plot | `geom_point` with size and color | Rank pathways by significance and gene ratio |
| Single-cell data | UMAP feature plot | `geom_point` with color gradient | Visualize cell types and gene expression |
## Common Failure Patterns and Solutions
| Failure Pattern | Observable Symptom | Likely Cause | Corrective Action |
| --- | --- | --- | --- |
| Samples cluster by batch in PCA | PCA plot shows batch grouping | Technical batch effects | Apply batch correction method |
| Heatmap shows no clear pattern | Uniform colors across all tiles | Genes not scaled per row | Apply z-score scaling per gene |
| Volcano plot points overlap heavily | Dense point cloud obscures structure | Too many genes plotted | Reduce point size and increase transparency |
| Pathway names truncated in dot plot | Long names cut off at plot edge | Insufficient plot width | Use `coord_flip` and adjust margins |
| Legend labels show variable names | Legend displays "condition" instead of "Treatment" | Default label mapping | Use `labs` to set descriptive labels |
| Figure text too small for journal | Font size below journal minimum | Default base_size too small | Set `base_size = 10` or higher in theme |
## Safety and Ethical Considerations in Data Visualization
### Avoiding Misleading Visualizations
Visualizations must accurately represent the underlying data without exaggeration or distortion. Manipulating axis limits, color scales, or point sizes to exaggerate differences constitutes scientific misconduct. The choice of color scale should reflect the data distribution, and axis limits should be stated clearly.
### Data Sharing and Reproducibility Requirements
Many journals and funding agencies require that raw sequencing data be deposited in public databases. The NCBI provides official descriptions of databases, search systems, sequence resources, and analysis services for data sharing [14]. The EMBL-EBI Training program offers bioinformatics learning pathways, data-resource training, and practical analysis education [15].
### Transparency in Analysis Methods
Manuscripts should describe the visualization methods used, including package versions and parameter choices. This transparency enables others to reproduce the figures and assess the validity of the visual interpretations.
## Building a Visualization Decision Framework for RNA-seq Plot Selection
Choosing the wrong plot type wastes analysis time and can obscure biological signals in RNA-seq data. A structured decision framework helps you match the analytical question, data structure, and audience expectations to the correct ggplot2 visualization. This section provides a practical framework for selecting plot types based on your specific analytical stage and the decisions you need to make from the data.
### Defining the Decision Criteria for Plot Selection
Before writing any ggplot2 code, evaluate three criteria that determine the appropriate visualization. First, identify the analytical question you need to answer. Second, assess the dimensionality of the data you are working with. Third, determine the target audience for the figure, whether that is a lab meeting, a manuscript reviewer, or a clinical collaborator.
The analytical question drives the plot type. Quality assessment questions require correlation heatmaps and count distribution boxplots. Global structure questions require PCA plots. Differential expression questions require volcano plots and MA plots. Pattern discovery questions require clustered heatmaps. Enrichment questions require dot plots and bar plots. Each question maps to a specific set of ggplot2 geometries and data transformations.
Data dimensionality determines whether you need faceting, color mapping, or multiple geometries. Gene-level comparisons across two conditions need simple boxplots or violin plots. Genome-wide comparisons need dimensionality reduction before plotting. Single-cell data with hundreds of thousands of cells requires downsampling or rasterization strategies before any ggplot2 geometry can render efficiently.
Audience expectations shape formatting choices. A lab meeting figure can use default themes and exploratory color scales. A manuscript figure requires the publication theme, specific font sizes, and accessible color palettes. A clinical collaborator may need simplified figures that emphasize interpretability over statistical detail.
### A Practical Plot Selection Matrix
The following decision matrix organizes plot selection by analytical stage and the specific decision you need to make. Use this matrix when you are uncertain which visualization will best support your analysis step.
| Analytical Stage | Decision Required | Recommended Plot | Data Structure Needed | Key ggplot2 Components |
| --- | --- | --- | --- | --- |
| Pre-analysis QC | Are there outlier samples? | Sample correlation heatmap | Normalized count matrix | `geom_tile`, `scale_fill_gradient2` |
| Pre-analysis QC | Are count distributions comparable? | Boxplot of log counts | Long-format counts with metadata | `geom_boxplot`, `facet_wrap` |
| Global structure | Do biological groups separate? | PCA score plot | PCA coordinates with metadata | `geom_point`, `stat_ellipse` |
| Global structure | Are there batch effects? | PCA colored by batch | PCA coordinates with batch variable | `geom_point`, `aes(color = batch)` |
| Differential expression | Which genes are significant? | Volcano plot | DE results with regulation column | `geom_point`, `scale_color_manual` |
| Differential expression | Is there expression-dependent bias? | MA plot | DE results with baseMean column | `geom_point`, `geom_hline` |
| Pattern discovery | What genes co-express? | Clustered heatmap | Scaled expression matrix | `geom_tile`, factor level ordering |
| Enrichment | Which pathways are enriched? | Dot plot | Enrichment results table | `geom_point`, `scale_size`, `scale_color` |
| Single-cell | What cell types exist? | UMAP colored by cell type | UMAP coordinates with metadata | `geom_point`, `aes(color = cell_type)` |
| Single-cell | Where is a gene expressed? | Feature plot | UMAP coordinates with gene expression | `geom_point`, `scale_color_gradient` |
### Implementing the Decision Framework in Practice
Apply the framework by working through a structured sequence of questions before you open R. Start with the analytical stage and the specific decision you need to make. Then confirm that your data structure supports the required plot type. Finally, verify that the plot will communicate the necessary information to your audience.
For example, if you are at the differential expression stage and need to identify significant genes, the volcano plot is the correct choice. Your data must include log2 fold change, adjusted p-value, and a regulation status column. The plot will communicate significance and direction of change to any reader familiar with genomic analysis.
If you are at the pattern discovery stage and need to identify co-expressed gene modules, the clustered heatmap is the correct choice. Your data must be scaled per gene and ordered by hierarchical clustering. The plot will reveal blocks of co-expressed genes that may share regulatory mechanisms or biological functions.
### Recording Visualization Decisions for Reproducibility
Document every plot selection decision in your analysis log. Record the analytical question, the plot type selected, the data transformation applied, and the rationale for the choice. This documentation supports reproducibility and helps collaborators understand why specific visualizations were chosen over alternatives.
```r
## Example decision log entry
## Date: 2025-03-10
## Analytical stage: Differential expression
## Decision: Volcano plot for identifying significant genes
## Data: DESeq2 results with 18,542 genes
## Transformation: Added regulation column with padj < 0.05 and |log2FC| > 1
## Rationale: Need to communicate significance and direction simultaneously
## Alternative considered: MA plot rejected because expression-dependent bias not the primary question
```
This record system ensures that figure choices are transparent and can be revisited if analysis parameters change. The Galaxy Training Network provides accessible workflow training and analysis tutorials that emphasize documentation and reproducibility as core practices [16].
### Troubleshooting Plot Selection Errors
When a visualization fails to communicate the intended message, the problem often lies in the plot selection instead of the code. A PCA plot that shows no group separation may indicate that the analytical question requires a different visualization, such as a heatmap of specific marker genes or a supervised analysis method.
A volcano plot with too many significant genes may not be the right choice for communicating the most important biological changes. In this case, a heatmap of the top differentially expressed genes or a pathway enrichment dot plot may better serve the analytical goal.
A heatmap with no visible pattern may indicate that the gene selection was too broad or that the scaling method obscured the biological signal. Consider whether a different gene set or a different normalization approach would better reveal the pattern of interest.
### Comparing Visualization Approaches for the Same Data
Different plot types can answer different questions from the same dataset. A study on PROS1 in glioma used multiple visualization approaches to characterize the role of this gene in the tumor immune microenvironment, including enrichment analyses visualized with ggplot2 and survival curves [2]. The choice of visualization depended on the specific question being asked about PROS1 expression and its clinical significance.
Similarly, a study on ADAM19 in systemic sclerosis used bulk RNA-seq data analyzed and visualized with multiple R packages including ggplot2 to examine the role of this gene in extracellular matrix remodeling and fibroblast activation [1]. The researchers selected different plot types for different analytical questions, demonstrating the importance of matching visualization to the specific biological question.
### Building a Visualization Checklist for Each Figure
Create a checklist that you apply to every figure before considering it complete. The checklist should confirm that the plot type matches the analytical question, that the data structure supports the geometry, that the aesthetic mappings communicate the intended variables, and that the formatting meets the target journal requirements.
```text
Figure checklist:
1. Does the plot type answer the analytical question?
2. Is the data in the correct format for the geometry?
3. Are all aesthetic mappings intentional and labeled?
4. Are axis labels descriptive and include units where relevant?
5. Is the color palette accessible to color-blind readers?
6. Are thresholds and reference lines clearly marked?
7. Does the figure meet journal formatting requirements?
8. Is the figure reproducible from the analysis script?
```
This checklist prevents common visualization errors and ensures that each figure serves its intended purpose. The EMBL-EBI Training program offers bioinformatics learning pathways that include practical guidance on data visualization and analysis education [15].
### When to Escalate Visualization Problems
If a plot type consistently fails to reveal expected patterns, the problem may lie in the upstream analysis instead of the visualization. Samples that do not cluster by condition in PCA may indicate a batch effect that requires correction before any visualization will be meaningful. Genes that show no differential expression may indicate a problem with the experimental design or the statistical model.
Escalate to a bioinformatics specialist when you cannot determine whether a visualization problem reflects a technical artifact or a genuine biological finding. Escalate to a statistics expert when different visualization approaches produce conflicting interpretations of the same data. Escalate to a domain expert when the visualization results contradict established biological knowledge in your field.
The nf-core documentation describes community pipeline standards and reproducible workflow context that can help identify whether visualization problems stem from upstream analysis steps [10]. The Carpentries lessons provide foundational training in computing and data handling that supports troubleshooting skills [9].
## Frequently Asked Questions
### What is the minimum data preparation needed before using ggplot2 for RNA-seq data?
You need three data structures: a count matrix with genes as rows and samples as columns, a sample metadata table with experimental conditions, and differential expression results if you plan to create volcano or MA plots. The count matrix must be reshaped to long format using `pivot_longer` for most ggplot2 visualizations. Sample identifiers in the metadata must exactly match column names in the count matrix.
### How do I choose between different normalization methods before visualization?
The choice depends on your analytical goal. Variance-stabilizing transformation and regularized log transformation from DESeq2 are appropriate for PCA and heatmaps because they stabilize variance across expression levels. TPM values work for comparing expression across genes within a sample. Raw counts require library size normalization before any visualization that compares across samples.
### Why do my PCA plots show samples clustering by batch instead of condition?
Batch effects occur when technical variation from library preparation, sequencing runs, or sample processing dominates biological variation. Visualize the PCA with points colored by batch to confirm this pattern. Apply batch correction methods such as ComBat-seq or include batch as a covariate in the differential expression model. After correction, regenerate the PCA to verify that samples cluster by biological condition.
### How can I make my ggplot2 figures meet journal formatting requirements?
Create a custom theme function that sets font family, font size, panel borders, and legend formatting. Export figures using `ggsave` with the dimensions and resolution specified by the target journal. Most journals require 300 dpi for raster images and specify figure widths in inches or centimeters. Use sans-serif fonts such as Arial or Helvetica at sizes between 8 and 12 points.
### What is the best way to visualize hundreds of differentially expressed genes in a heatmap?
Filter to the most significant genes based on adjusted p-value and fold change, typically the top 50 to 100 genes. Scale expression values per gene using z-scores so that the color scale is comparable across genes. Order genes and samples by hierarchical clustering to reveal co-expression patterns. Consider splitting the heatmap into multiple panels if the gene list is very large.
### How do I add gene labels to a volcano plot without overlapping text?
Use the `ggrepel` package with `geom_text_repel`. This function automatically adjusts label positions to minimize overlap. Select a subset of genes to label, such as the top 20 genes by fold change or specific genes of biological interest. Set `max.overlaps` to control how many labels are placed in dense regions.
### Can ggplot2 handle single-cell RNA-seq visualization, or do I need specialized packages?
ggplot2 can create publication-ready single-cell visualizations when you extract the dimensionality reduction coordinates and metadata from Seurat or similar objects. The SCpubr package provides concise function calls for generating high-quality visualizations commonly used in single-cell transcriptome analyses, building on working knowledge of Seurat and ggplot2 [7]. For complex single-cell visualizations, specialized packages may save time, but ggplot2 offers maximum customization.
### How do I visualize gene set enrichment results with ggplot2?
The clusterProfiler package generates enrichment results that can be plotted directly with ggplot2. For dot plots, map gene ratio to the x-axis, pathway names to the y-axis, gene count to point size, and adjusted p-value to point color. For GSEA results, create bar plots of normalized enrichment scores with pathways ordered by score. The PROS1 study in glioma used clusterProfiler for GO and KEGG analyses and visualized results with ggplot2 [2].
## Related Bioinformatics Guides
- [RNA-Seq Visualization: Volcano Plots, Heatmaps, and PCA](/knowledge/bioinformatics/rna-seq-visualization-volcano-plots-heatmaps-and-pca)
- [Single-Cell RNA Sequencing Quality Control: A Practical Guide to Filtering and Metrics](/knowledge/bioinformatics/single-cell-rna-sequencing-quality-control-a-practical-guide-to-filtering-and-metrics)
- [RNA-Seq Databases: Accessing and Using Public RNA-Seq Data](/knowledge/bioinformatics/rna-seq-databases-accessing-and-using-public-rna-seq-data)
- [RNA-Seq Data Analysis in Galaxy: A User-Friendly Platform](/knowledge/bioinformatics/rna-seq-data-analysis-in-galaxy-a-user-friendly-platform)
- [RNA-Seq Data Analysis Workflow: From Raw Reads to Insights](/knowledge/bioinformatics/rna-seq-data-analysis-workflow-from-raw-reads-to-insights)
## Related Clinical & Scientific Guides
* [A Practical Guide to Detecting Antimicrobial Resistance Genes in Shotgun Metagenomic Data](/knowledge/bioinformatics/a-practical-guide-to-detecting-antimicrobial-resistance-genes-in-shotgun-metagenomic-data)
* [Computational Immunology: Modeling the Immune System](/knowledge/bioinformatics/computational-immunology-modeling-the-immune-system)
* [How to Set Hard Filters for Germline Variant Calling: A Practical Guide to GATK Best Practices](/knowledge/bioinformatics/how-to-set-hard-filters-for-germline-variant-calling-a-practical-guide-to-gatk-best-practices)
## References and Further Reading
[1] [ADAM Metallopeptidase domain 19 promotes skin fibrosis in systemic sclerosis via neuregulin-1.](https://pubmed.ncbi.nlm.nih.gov/39716051). Molecular medicine (Cambridge, Mass.), 2024.
[2] [PROS1 shapes the immune-suppressive tumor microenvironment and predicts poor prognosis in glioma.](https://pubmed.ncbi.nlm.nih.gov/36685506). Frontiers in immunology, 2022.
[3] [Identification of Key circRNAs/lncRNAs/miRNAs/mRNAs and Pathways in Preeclampsia Using Bioinformatics Analysis.](https://pubmed.ncbi.nlm.nih.gov/30833538). Medical science monitor : international medical journal of experimental and clinical research, 2019.
[4] [Transcriptome RNA sequencing reveals the global molecular responses and circRNA-miRNA-lncRNA interaction network in cataract.](https://pubmed.ncbi.nlm.nih.gov/41128953). Japanese journal of ophthalmology, 2026.
[5] [RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome](https://doi.org/10.1186/1471-2105-12-323). BMC Bioinformatics, 2011.
[6] [Batch correcting single-cell spatial transcriptomics count data with Crescendo improves visualization and detection of spatial gene patterns](https://doi.org/10.1186/s13059-025-03479-9). Genome Biology, 2025.
[7] [SCpubr: a user-friendly R-package for generating publication-ready visualizations of single-cell transcriptome analyses.](https://pubmed.ncbi.nlm.nih.gov/42327689). Bioinformatics advances, 2026.
[8] [ggcoverage: an R package to visualize and annotate genome coverage for various NGS data.](https://pubmed.ncbi.nlm.nih.gov/37559015). BMC bioinformatics, 2023.
[9] [The Carpentries Lessons](https://carpentries.org/lessons). The Carpentries.
[10] [nf-core Documentation](https://nf-co.re/docs). nf-core.
[11] [RAGER: A user-friendly computational platform for integrated analysis of RNA-Seq and ATAC-seq data.](https://doi.org/10.1371/journal.pone.0349941). 2026.
[12] [Molecular Characteristics of the Endometrium in Polycystic Ovary Syndrome with Insulin Resistance.](https://doi.org/10.1007/s43032-026-02069-9). 2026.
[13] [Bioconductor](https://bioconductor.org/). Bioconductor Project.
[14] [NCBI Data Resources](https://www.ncbi.nlm.nih.gov/). National Center for Biotechnology Information.
[15] [EMBL-EBI Training](https://www.ebi.ac.uk/training). European Bioinformatics Institute.
[16] [Galaxy Training Network](https://training.galaxyproject.org/). Galaxy Project.
> This article is educational and does not replace validated analysis plans, institutional policy, clinical interpretation, or specialist review.