PaperlayerView exact revision on arXiv

Source-linked AI summary

Count-based differential expression analysis of RNA sequencing data using R and Bioconductor

Simon Anders, Davis J. McCarthy, Yunshen Chen, Michal Okoniewski, Gordon K. Smyth, Wolfgang Huber, Mark D. Robinson

arXiv:1302.3685v3q-bio.GN

TL;DR

RNA-seq differential expression requires careful counting, variability modeling, quality control, and experimental-design specification, especially when replication is limited. This protocol organizes those tasks into an R/Bioconductor workflow using DESeq and edgeR, with support for common two-group comparisons and blocking factors. It provides a practical standard workflow but leaves more complex designs and isoform-level questions to further adaptation and specialized methods.

  • Problem

    RNA-seq differential expression needs guidance on read counting, biological variability, quality control, and statistical modeling across experimental conditions and covariates.

  • Method

    The protocol presents a count-based R/Bioconductor workflow using DESeq and edgeR, including negative-binomial modeling and generalized linear models for blocking factors.

  • Results

    The workflow provides a standard RNA-seq analysis for two-group comparisons while modeling blocking factors through generalized linear models when needed.

  • Takeaways & Limitations

    Users can apply the workflow to common count-based differential analyses while adapting its modular stages and design specification to their experiments.

  • Takeaways & Limitations

    Gene-level count-based methods may miss true differential expression when isoforms of one gene change in opposite directions.

Abstract

from arXiv · show

RNA sequencing (RNA-seq) has been rapidly adopted for the profiling of transcriptomes in many areas of biology, including studies into gene regulation, development and disease. Of particular interest is the discovery of differentially expressed genes across different conditions (e.g., tissues, perturbations), while optionally adjusting for other systematic factors that affect the data collection process. There are a number of subtle yet critical aspects of these analyses, such as read counting, appropriate treatment of biological variability, quality control checks and appropriate setup of statistical modeling. Several variations have been presented in the literature, and there is a need for guidance on current best practices. This protocol presents a "state-of-the-art" computational and statistical RNA-seq differential expression analysis workflow largely based on the free open-source R language and Bioconductor software and in particular, two widely-used tools DESeq and edgeR. Hands-on time for typical small experiments (e.g., 4-10 samples) can be <1 hour, with computation time <1 day using a standard desktop PC.

INTRODUCTION

This protocol presents an accessible, modular workflow for count-based RNA-seq differential expression, centered on two-group comparisons while accommodating blocking factors and extensions. It emphasizes quality assessment, statistical modeling for modest replication, and adaptation to more complex or isoform-aware analyses.

  • Workflow: Reads are counted at annotated features, and statistical analysis of the count table identifies expression changes between experimental groups.The direct counting strategy also applies to other count-based datasets, including ChIP-seq, proteomic spectral counts, and metagenomics.
  • Statistical framework: Negative binomial models are used because biological replication is often modest and Poisson variation captures technical rather than biological variability.The protocol focuses on DESeq and edgeR, which implement NB-based differential analyses with different dispersion estimation approaches.
  • Software and modularity: R and Bioconductor integrate statistical algorithms, visualization, annotation resources, and interactive quality checks across the workflow.The workflow can incorporate mature tools outside Bioconductor and alternative aligners or feature-counting strategies.
  • Scope and extensions: Gene-level counting can miss differential expression when isoforms from the same gene change in opposite directions.Exon-level, transcript-level, and isoform-deconvolution methods address related limitations, but isoform-specific estimation remains difficult in several settings.
  • Scope of this protocol: The workflow targets RNA-seq experiments comparing two treatment groups and can account for a blocking factor.The protocol is designed as a concise standard workflow for new users and R newcomers.
  • Experimental design considerations: More complicated designs require an appropriately specified design matrix representing biological factors and covariates such as batch effects.The protocol states that users must adapt the pipeline and consult broader documentation for tasks beyond the standard workflow.

MATERIALS

The materials support a Unix-like, R/Bioconductor-based RNA-seq workflow with splice-aware alignment, read counting, visualization, and example data. Software versions and input formats require attention because recommendations and compatibility can change.

  • Operating system: The protocol assumes Linux or macOS with a bash-like shell, although Windows use through Cygwin is possible but not covered.Commands are intended for execution in a terminal window.
  • Software: Required software includes a splice-aware aligner, R, Bioconductor packages, samtools, and HTSeq; IGV or Savant are optional visualization tools.The illustrated aligner is TopHat2, but several alternatives are listed.
  • Versioning: Users should install recent stable software versions and consult version-specific documentation because package features and recommendations change over time.R and Bioconductor versions are matched during package installation.
  • Input files: The workflow starts generally from compressed FASTQ files, while modifications for reads from other sequencing platforms are not discussed.FASTQ is identified as the common starting format for Illumina sequencing data.
  • Example data: The example dataset contains seven Drosophila S2-cell RNA-seq samples: three pasilla knockdowns and four untreated controls.The analysis goal is to identify genes changing between Knockdown and Control.
  • Experimental design: Library type can be modeled as a blocking factor, turning the example from a simple two-group comparison into a complex design analyzed with generalized linear models.Both edgeR and DESeq provide functions for fitting these GLMs.

Download the example data

The example-data procedure retrieves selected SRA samples, converts them to FASTQ, and prepares reference-genome resources for alignment. Metadata-driven commands automate sample selection and downloads, while paired-end reads require explicit splitting.

  • Sample selection: Only the four Control and three Knockdown samples from the example SRA dataset are selected for download.The dataset is identified by GEO accession GSE18508 and corresponding SRA records.
  • Metadata-driven download: SraRunInfo.csv supplies metadata used in R to select matching libraries and construct download paths.The example filters library names containing CG8144 or Untreated-.
  • Downloading: R can automate downloading the selected SRA files, while wget, curl, or a web browser provide alternatives.The download procedure checks that files exist before conversion.
  • SRA conversion: The SRA Toolkit converts compressed SRA files to FASTQ using fastq-dump, with R optionally constructing and executing the shell commands.The protocol illustrates processing the downloaded files in a loop.
  • SRA conversion: The --split-3 option separates mate-pair reads into distinct files, producing two FASTQ files for paired-end data and one for single-end data.Paired-end outputs receive 1.FASTQ and 2.FASTQ suffixes.
  • Reference genome: Reference genomes should be downloaded in genomic FASTA form rather than cDNA, without repeat masking, for genome alignment.The protocol lists Ensembl, Ensembl Plants, Ensembl Metazoa, and UCSC as reference sources.

Get gene model annotations

The workflow begins by obtaining a gene-model annotation and matching it to the reference sequence, then constructs and verifies sample metadata before downstream analysis.

  • Download a GTF file containing gene models for the organism of interest.
  • Use matching assembly coordinates for the reference FASTA and GTF annotation, preferably downloading both from the same provider.
  • A UCSC Table Browser-generated GTF can be used, but its gene IDs may require manual correction before downstream tools work correctly.
  • Preprocess the reference FASTA into an aligner-specific index before read alignment.
  • Construct one metadata row per sample, add experimental groupings, and verify filenames, labels, conditions, and populated columns before statistical analysis.

Align the reads (using tophat2) to reference genome

Reads are aligned to the indexed reference with tophat2 using annotation and sample-specific FASTQ inputs, after which mapped files are organized for downstream counting and visualization.

  • Construct tophat2 commands from the metadata table, specifying the annotation GTF, thread count, output directory, reference index, and sample FASTQ files.
  • For paired-end data, provide the two FASTQ argument lists in matching mate order, with filenames concatenated by commas without spaces.
  • The annotation supplied with -G facilitates mapping reads across exon-exon junctions, including some junctions found de novo.
  • Organize accepted_hits.bam files by creating name-sorted files for htseq-count and position-sorted, indexed BAM files for IGV.
  • Sorting and conversion operations can consume substantial disk space because multiple SAM and BAM derivatives are created for each sample.

Inspect alignments with IGV

The workflow uses IGV to inspect mapped reads against the annotation and known controls before generating count files for differential-expression analysis.

  • Load the position-sorted BAM files and GTF annotation into IGV after selecting the appropriate reference genome.
  • Inspect expressed transcripts at read-level resolution to check alignment at and across exon-exon junctions expected from the annotation.
  • Use known positive and negative controls, when available, to verify that relative read density differs in the expected regions.
  • Run htseq-count on the name-sorted SAM files with the selected GTF annotation, recording each resulting COUNT file in the metadata table.
  • Set htseq-count options according to the experiment: -s specifies strand protocol status, while -a sets a minimum alignment-quality score.
  • Choose edgeR or DESeq procedures according to whether the experimental design is simple or complex.

A. edgeR - simple design

For a simple two-group design, edgeR imports and filters count data, normalizes and examines sample relationships, estimates dispersion, and tests differential expression.

  • Read the HTSeq count files into edgeR, remove non-informative features, and create a DGEList container for the retained counts.
  • In edgeR, retain features with at least 1 read per million in the number of samples equal to the smallest replicate group.
  • Estimate normalization factors and inspect relationships among samples using a multidimensional scaling plot.
  • For simple designs, estimate tagwise dispersion and visualize the mean-variance relationship with plotMeanVar and plotBCV.
  • Test for differential expression, summarize results with topTags, inspect normalized expression for leading genes, and visualize selected genes in an M-versus-A plot.
  • For designs including additional factors, specify them in a design matrix, estimate design-adjusted dispersion, fit feature-wise GLMs, and perform a likelihood-ratio test.

C. DESeq - simple design

The simple-design workflow prepares count data and metadata, assesses sample relationships and dispersion, then tests differential expression and reviews the results with visual and statistical checks.

  • DESeq creates a CountDataSet from count tables and corresponding sample metadata.
  • The workflow estimates normalization factors, inspects sample relationships with variance stabilization and PCA, and examines estimated dispersions.
  • Differential expression is tested with nbinomTest, followed by MA-plot visualization of log-fold-changes against expression strength.
  • Genes with adjusted P-values below 0.1 are inspected by fold-change direction, and their number is counted.
  • Results are saved persistently and checked with a histogram of unadjusted P-values.

D. DESeq - complex design

The complex-design workflow estimates design-relative dispersions, compares full and reduced negative-binomial models, adjusts P-values, and inspects significant genes and their genomic read evidence.

  • For complex designs, dispersion estimates are calculated relative to the specified factors using the CR adjusted profile likelihood.
  • Differential expression is tested by fitting a full model and a reduced model with the factor of interest removed.
  • Likelihood-ratio statistics and associated P-values are computed from the two fitted models, then P-values are adjusted for multiple testing.
  • The result table combines full-model fits with raw and adjusted P-values, and significant genes are inspected at FDR below 10%.
  • Top differentially expressed genes are spot-checked in IGV using GTF and BAM files to assess whether counting and statistics are appropriately represented.

TIMING

The protocol’s runtime is concentrated in quality checks, alignment, and feature counting, while differential analysis is comparatively quick; downstream checks assess data quality and modeling behavior.

  • TIMING: Approximately 10 hours are required for the full workflow on eight cores and 8 GB of RAM, with computation dominated by quality checks, alignment, and feature counting.
  • TIMING: Quality checks take ∼2 h, alignment ∼6 h, feature counting ∼3 h, and differential analysis is often <20 min.
  • Quality checks: Sequencing-quality reports should be inspected for persistent low-quality scores and adapter overrepresentation before removing samples, trimming reads, or modifying alignment parameters.
  • Feature counting: Count tables contain integer read-overlap counts, and independent filtering of uninformative or low-count features is typically beneficial for subsequent differential-expression power.
  • Normalization: Normalization scales counts for comparable modeling without altering raw read counts, while sample-specific factors are incorporated into the negative-binomial model.
  • Sample relations: MDS or PCA sample-relationship checks can reveal systematic library-type effects and biological outliers that influence dispersion and differential-expression discovery.
  • Sanity checks: A useful P-value histogram has a left peak for strong differential expression, an approximately uniform [0, 1] floor for null genes, and sometimes an upper-end peak from low counts.

FIGURE LEGENDS

The paper’s figures and supporting passages depict a modular count-based RNA-seq differential-expression workflow, its statistical foundations, and diagnostic plots for sample relationships, dispersion, expression changes, and testing.

  • Statistical model: The negative binomial model accommodates biological variability and supports generalized linear models for complex designs and covariate-adjusted expression testing.Its dispersion parameter represents over-dispersion relative to the Poisson distribution, while the biological coefficient of variation is its square root.
  • Statistical model: Stable gene-wise dispersion estimates require information sharing across genes, through averaging, global trends, or compromises between gene-specific and global estimates.The paper emphasizes dispersion estimation as critical for reliable statistical testing when sample numbers are limited.
  • Workflow: The workflow moves from read sequences through feature counting to differential-expression discovery, with quality checks throughout and alternative entry points for alignments or count tables.The protocol supports simple and complex experimental designs, while edgeR and DESeq share many steps but use different commands.
  • Diagnostic plots: Diagnostic figures visualize sample relations with MDS or PCA, mean-variance and dispersion patterns, log-fold changes against expression strength, and gene-level P-value distributions.The plots use count-specific distances, vst-transformed counts, normalized expression measures, and log-average read counts to support quality assessment.
  • Counting and normalization: Accurate inference requires raw read counts, correct handling of paired-end fragments and ambiguous assignments, and normalization factors used within the statistical model rather than altering raw counts.The counting and normalization choices are presented as prerequisites for valid negative-binomial inference.

SUPPLEMENTARY MATERIAL

The supplementary material provides example analysis files and records acknowledgements, funding, competing interests, and author contributions.

  • Supplementary files: Supplementary File 1 contains intermediate COUNT files, the count table, metadata table, and the original SraRunInfo CSV downloaded from the Short Read Archive.These files support the example workflow and statistical analysis.
  • Credits and disclosures: The paper acknowledges contributors and funding sources, declares no competing financial interests, and identifies package authors and protocol contributors.The author-contribution statement distinguishes contributors associated with DESeq, edgeR, and protocol development.
Loading 1302.3685v3…