nf-core/riboseq

Edit

Introduction

This document describes the output produced by the pipeline. Most of the plots are taken from the MultiQC report generated from the full-sized test dataset for the pipeline using a command similar to the one below:

nextflow run nf-core/riboseq -profile test_full,<docker/singularity/institute>

The directories listed below will be created in the results directory after the pipeline has finished. All paths are relative to the top-level results directory.

Pipeline overview

The pipeline is built using Nextflow and processes data using the following steps:

• nf-core/riboseq: Output
  • Introduction
  • Pipeline overview
  • Preprocessing
    • cat
    • FastQC
    • UMI-tools extract
    • TrimGalore
    • fastp
    • BBSplit
    • rRNA removal
  • Alignment
    • STAR
  • UMI-tools deduplication
  • Riboseq-specific QC
    • Ribo-TISH quality
    • Ribotricer detect-orfs QC outputs
  • ORF predictions
    • Ribo-TISH predict
    • Ribotricer detect-orfs
  • P-site identification
    • riboWaltz
  • Quantification
  • Translational efficiency
    • MultiQC
  • Workflow reporting and genomes
    • Reference genome files
    • Pipeline information

Preprocessing

cat

If multiple libraries/runs have been provided for the same sample in the input samplesheet (e.g. to increase sequencing depth) then these will be merged at the very beginning of the pipeline in order to have consistent sample naming throughout the pipeline. Please refer to the usage documentation to see how to specify these samples in the input samplesheet.

fq lint

fq lint runs several checks on input FASTQ files. It will fail with a non-zero error code when issues are found, which will terminate the workflow execution. In the absence of this, the successful linting produces the logs you will find here.

FastQC

FastQC gives general quality metrics about your sequenced reads. It provides information about the quality score distribution across your reads, per base sequence content (%A/T/G/C), adapter contamination and overrepresented sequences. For further reading and documentation see the FastQC help pages.

UMI-tools extract

UMI-tools and UMICollapse deduplicate reads based on unique molecular identifiers (UMIs) to address PCR-bias. Firstly, the UMI-tools extract command removes the UMI barcode information from the read sequence and adds it to the read name. Secondly, reads are deduplicated based on UMI identifier after mapping as highlighted in the UMI dedup section.

To facilitate processing of input data which has the UMI barcode already embedded in the read name from the start, --skip_umi_extract can be specified in conjunction with --with_umi.

TrimGalore

Trim Galore! is a wrapper tool around Cutadapt and FastQC to peform quality and adapter trimming on FastQ files. Trim Galore! will automatically detect and trim the appropriate adapter sequence. It is the default trimming tool used by this pipeline, however you can use fastp instead by specifying the --trimmer fastp parameter. You can specify additional options for Trim Galore! via the --extra_trimgalore_args parameters.

NB: TrimGalore! will only run using multiple cores if you are able to use more than > 5 and > 6 CPUs for single- and paired-end data, respectively. The total cores available to TrimGalore! will also be capped at 4 (7 and 8 CPUs in total for single- and paired-end data, respectively) because there is no longer a run-time benefit. See release notes and discussion whilst adding this logic to the nf-core/atacseq pipeline.

fastp

fastp is a tool designed to provide fast, all-in-one preprocessing for FastQ files. It has been developed in C++ with multithreading support to achieve higher performance. fastp can be used in this pipeline for standard adapter trimming and quality filtering by setting the --trimmer fastp parameter. You can specify additional options for fastp via the --extra_fastp_args parameter.

BBSplit

BBSplit is a tool that bins reads by mapping to multiple references simultaneously, using BBMap. The reads go to the bin of the reference they map to best. There are also disambiguation options, such that reads that map to multiple references can be binned with all of them, none of them, one of them, or put in a special “ambiguous” file for each of them.

This functionality would be especially useful, for example, if you have mouse PDX samples that contain a mixture of human and mouse genomic DNA/RNA and you would like to filter out any mouse derived reads.

The BBSplit index will have to be built at least once with this pipeline by providing --bbsplit_fasta_list which has to be a file containing 2 columns: short name and full path to reference genome(s):

mm10,/path/to/mm10.fa
ecoli,/path/to/ecoli.fa
sarscov2,/path/to/sarscov2.fa

You can save the index by using the --save_reference parameter and then provide it via --bbsplit_index for future runs. As described in the Output files dropdown box above the FastQ files relative to the main reference genome will always be called *primary*.fastq.gz.

rRNA removal

When --remove_ribo_rna is specified (default: true), the pipeline removes ribosomal RNA reads using one of three tools, selected via the --ribo_removal_tool parameter.

SortMeRNA (default)

When --ribo_removal_tool sortmerna is used (the default), the pipeline uses SortMeRNA for the removal of ribosomal RNA. SortMeRNA uses k-mer matching against rRNA databases. By default, rRNA databases defined in the SortMeRNA GitHub repo are used. You can see an example in the pipeline GitHub repository in assets/rrna-db-defaults.txt which is used by default via the --ribo_database_manifest parameter. Please note that commercial/non-academic entities require licensing for SILVA for these default databases.

Bowtie2

When --ribo_removal_tool bowtie2 is specified, the pipeline uses Bowtie2 for alignment-based rRNA removal. Reads are aligned against rRNA reference sequences specified via --ribo_database_manifest, and reads that align to rRNA are filtered out. The unaligned reads (non-rRNA) are kept for downstream analysis.

RiboDetector

When --ribo_removal_tool ribodetector is specified, the pipeline uses RiboDetector for machine learning-based rRNA detection. RiboDetector uses a pre-trained neural network to identify rRNA reads without requiring a reference database. This makes it particularly useful for organisms that lack well-characterized rRNA sequences, or when you want to avoid database licensing requirements. The tool automatically determines the appropriate read length from your data.

Alignment

STAR

STAR is a read aligner designed for splice aware mapping typical of RNA sequencing data. STAR stands for Spliced Transcripts Alignment to a Reference, and has been shown to have high accuracy and outperforms other aligners by more than a factor of 50 in mapping speed, but it is memory intensive. Using --aligner star_salmon is the default alignment and quantification option.

The STAR section of the MultiQC report shows a bar plot with alignment rates: good samples should have most reads as Uniquely mapped and few Unmapped reads.

The original BAM files generated by the selected alignment algorithm are further processed with SAMtools to sort them by coordinate, for indexing, as well as to generate read mapping statistics.

UMI dedup

After extracting the UMI information from the read sequence (see UMI-tools extract), the second step in the removal of UMI barcodes involves deduplicating the reads based on both mapping and UMI barcode information. UMI deduplication can be carried out either with UMI-tools or UMICollapse, set via the umi_dedup_tool parameter. The output BAM files are the same, though UMI-tools has some additional outputs, as described above. Either method will generate a filtered BAM file after the removal of PCR duplicates.

Riboseq-specific QC

Read distribution metrics around annotated protein coding regions or based on alignments alone, plus related metrics.

Ribo-TISH quality

Ribo-TISH quality control analyzes the distribution of ribosome-protected fragment (RPF) reads around annotated coding sequences. Key metrics include read length distribution and reading frame periodicity, which are essential for assessing ribosome profiling data quality.

Ribotricer detect-orfs QC outputs

ORF predictions

Ribo-TISH predict

Ribotricer detect-orfs

P-site identification

riboWaltz

P-sites are identified by passing transcriptome-level alignment BAM files to riboWaltz, producing the following outputs:

Quantification

Quantification is done by passing transcriptome-level alignment BAM files to Salmon, producing the following outputs:

Raw outputs from Salmon are available for each sample:

Translational efficiency

The pipeline supports two methods for translational efficiency analysis: anota2seq (default) and deltaTE. The method used depends on the --translational_efficiency_method parameter.

anota2seq outputs

anota2seq produces the following outputs:

deltaTE outputs

When using --translational_efficiency_method deltate, the pipeline produces the following outputs:

Fold change plots

Both methods produce fold change plots that visualize the relationship between RNA-seq and Ribo-seq changes. The examples below were generated using the test_full profile on the same dataset, allowing direct comparison between methods.

anota2seq fold change plot:

The anota2seq fold change plot shows total mRNA log2FC on the X-axis and translated mRNA log2FC on the Y-axis, with genes colored by their regulatory classification. The solid diagonal line represents equal changes in both data types.

In the test dataset: The plot reveals 623 genes with significant translational regulation across several categories:

• Translation down (357) / up (225): Genes showing pure translational regulation without corresponding RNA changes, appearing as vertical deviations from the diagonal
• mRNA abundance down (16) / up (25): Genes with coordinated RNA and Ribo-seq changes, appearing along the diagonal
• Buffered (mRNA down) (0) / (mRNA up) (0): Genes where translation dampens RNA changes (none detected in this dataset)

The predominance of translation-regulated genes (582 of 623) indicates that translational control is the dominant regulatory mechanism in this experimental system.

deltaTE scatter plot:

The deltaTE method produces a similar visualization showing total mRNA log2FC (RNA-seq) on the X-axis and translated mRNA log2FC (Ribo-seq) on the Y-axis, with genes colored by their classification. The diagonal dashed line represents equal changes in both data types.

In the test dataset: The plot reveals 1,694 DTEGs across several categories:

• Translation down (895) / up (710): Genes showing translational regulation without corresponding RNA changes, appearing as vertical deviations from the diagonal along the y-axis
• mRNA abundance down (4) / up (30): Genes with coordinated RNA and Ribo-seq changes, appearing along the diagonal
• Buffering down (14) / up (16): Genes where translation dampens RNA changes, appearing below the diagonal
• Intensified (25): Genes where translation amplifies RNA changes, appearing above the diagonal

The striking vertical distribution of points (translation category) demonstrates that translational regulation is the dominant mode of gene expression control in this experimental system, with most changes occurring in Ribo-seq independently of RNA-seq.

Diagnostic plots

Both methods provide diagnostic plots for quality control and model assessment:

anota2seq diagnostic plots

anota2seq produces several QC plots for assessing model fit and statistical assumptions.

Residual distribution summary:

This plot summarizes the distribution of regression residuals across quantiles of the standard normal distribution. It helps assess whether the residuals follow the expected normal distribution, which is an assumption of the linear model.

In the test dataset: The obtained outlier rate (1.588%) is close to the expected 1%, indicating the model assumptions are reasonably well met.

Residuals vs fitted values:

This classic diagnostic plot shows residuals plotted against fitted values. Ideally, residuals should be randomly scattered around zero with constant variance (homoscedasticity).

In the test dataset: The residuals show a characteristic funnel shape with greater spread at lower fitted values, which is typical for count-based expression data. This heteroscedasticity is expected and accounted for by the RVM approach.

RVM fit for interactions:

The Random Variance Model (RVM) fit plot shows how well the empirical variance distribution matches the theoretical F-distribution. The left panel shows a Q-Q plot of empirical vs theoretical quantiles, and the right panel shows the cumulative distribution functions with a Kolmogorov-Smirnov test p-value.

In the test dataset: The KS p-value of 4.24e-09 indicates significant deviation from the theoretical distribution, suggesting the variance structure in this dataset differs from standard assumptions. This is not uncommon in real experimental data and the RVM approach helps account for such deviations.

deltaTE visualizations

The deltaTE method produces diagnostic plots for quality control and interpretation:

• Fold change plot: The main visualization showing RNA-seq vs Ribo-seq log2 fold changes, with genes colored by their deltaTE classification (mRNA_abundance, translation, buffering, intensified). This provides direct visual comparison with anota2seq’s equivalent plot.
• Interaction p-value distribution: Histogram showing the distribution of adjusted p-values from the DESeq2 interaction term analysis.
• PCA plots: Principal component analysis plots for both Ribo-seq and RNA-seq samples to assess sample relationships and identify potential batch effects.
• Heatmap: Clustered heatmap of the most significantly regulated genes, showing expression patterns across samples.

PCA plots:

PCA plots help assess sample clustering and identify potential batch effects or outliers. In well-designed experiments, samples should cluster primarily by treatment condition rather than technical factors.

In the test dataset: Both PCA plots show clear separation between control (coral) and treated (teal) samples along PC1, which captures 57% of variance for Ribo-seq and 39% for RNA-seq. This indicates that the treatment effect is the dominant source of variation in both data types. The treated samples cluster together on the right side of PC1 in both plots, while control samples are distributed on the left, with one control sample (SRX11780879 in RNA-seq, SRX11780886 in Ribo-seq) showing some separation from other controls along PC2.

Expression heatmap:

This clustered heatmap shows the expression patterns of the top 50 differentially translated genes (DTEGs) across all samples, with rows representing genes and columns representing samples. The color scale indicates normalized expression levels (Z-scores). Samples are annotated by condition (control/treated) and sequencing type (rnaseq/riboseq).

In the test dataset: The heatmap reveals distinct expression patterns between control and treated samples. Samples cluster hierarchically first by condition (control vs treated), with Ribo-seq control samples on the far left showing predominantly high expression (red), and treated samples showing lower expression (blue/purple) for most genes. The top 50 DTEGs show strong coordinated downregulation in treated samples compared to controls, particularly evident in the Ribo-seq data, reflecting the predominance of translation-down regulation in this dataset.

Interpreting deltaTE diagnostic plots

Interaction p-value distribution:

This histogram shows the distribution of adjusted p-values from the interaction term in the DESeq2 model, which tests for differences in translational efficiency between conditions. For well-powered experiments with genuine translational regulation:

• You should see a mixture of p-values with enrichment near 0 (significant genes) and relatively uniform distribution across higher p-values
• A completely uniform distribution suggests no translational regulation effects
• Heavy skewing toward 1.0 may indicate underpowered analysis or technical issues

In the test dataset: The distribution shows strong enrichment of small adjusted p-values (large peak near 0), indicating widespread translational regulation in this experimental system. The dashed red line marks the significance threshold (α = 0.05). The analysis identified 1,694 differentially translated genes (DTEGs), with translational regulation being the dominant mode of gene expression control.

The deltaTE approach focuses on the key fold change visualization that directly shows the relationship between RNA and ribosome changes, making it easy to interpret translational regulation patterns. The intensified category (genes where translation amplifies RNA changes in the same direction) is unique to deltaTE and provides additional biological insights not available in anota2seq.

MultiQC

MultiQC is a visualization tool that generates a single HTML report summarising all samples in your project. Most of the pipeline QC results are visualised in the report and further statistics are available in the report data directory.

Results generated by MultiQC collate pipeline QC from supported tools i.e. FastQC, Cutadapt, SortMeRNA, STAR, RSEM, HISAT2, Salmon, SAMtools, Picard, RSeQC, Qualimap, Preseq and featureCounts. Additionally, various custom content has been added to the report to assess the output of dupRadar, DESeq2 and featureCounts biotypes, and to highlight samples failing a mimimum mapping threshold or those that failed to match the strand-specificity provided in the input samplesheet. The pipeline has special steps which also allow the software versions to be reported in the MultiQC output for future traceability. For more information about how to use MultiQC reports, see http://multiqc.info.

Workflow reporting and genomes

Reference genome files

A number of genome-specific files are generated by the pipeline because they are required for the downstream processing of the results. If the --save_reference parameter is provided then these will be saved in the genome/ directory. It is recommended to use the --save_reference parameter if you are using the pipeline to build new indices so that you can save them somewhere locally. The index building step can be quite a time-consuming process and it permits their reuse for future runs of the pipeline to save disk space.

Pipeline information

Nextflow provides excellent functionality for generating various reports relevant to the running and execution of the pipeline. This will allow you to troubleshoot errors with the running of the pipeline, and also provide you with other information such as launch commands, run times and resource usage.