Nextflow RNAseq pipeline

Simple RNA-Seq workflow (from Seqera training materials)

To demonstrate a real-world biomedical scenario, we will implement a proof-of-concept RNA-seq workflow that indexes a transcriptome, checks read quality, quantifies expression, and produces a summary report.

We will build it up across a series of seven scripts (script1.nf – script7.nf), each one adding a step to the previous. The workflow uses three third-party tools that may be new to you:

• Salmon — quantifies transcripts from RNA-seq data.
• FastQC — performs quality control on high-throughput sequencing data.
• MultiQC — searches a directory for analysis logs and compiles them into a single HTML report.

These may not be the exact tools you will use in your own field, but they can be swapped for any comparable tool — that is the power of Nextflow.

The pipeline at a glance

Before we build it piece by piece, here is the whole workflow. It takes a reference transcriptome and paired-end RNA-seq reads, and produces per-sample expression estimates plus a single quality report. Each box is one step; the arrows show what feeds into what.

%%{init: {'themeVariables': {'lineColor': '#0f7040'}}}%%
flowchart TD
    T([transcriptome.fa]) --> IDX["INDEX — salmon index"]
    IDX --> SI([salmon_index])
    R([paired-end reads]) -->|fromFilePairs| RP([read pairs, per sample])
    SI --> Q["QUANTIFICATION — salmon quant"]
    RP --> Q
    RP --> FQ["FASTQC — quality control"]
    Q --> M["MULTIQC — summary report"]
    FQ --> M
    Q --> QR([quantification results])
    M --> REP([multiqc_report.html])

    linkStyle default stroke-width:4px

Step by step — and why each tool is used when:

1. Index the transcriptome reference — salmon index. Salmon cannot quantify reads against a plain FASTA; it first needs to turn the reference into a fast, searchable index. We do this once at the start, and every sample reuses the same index.

2. Quantify expression — salmon quant. For each sample, Salmon maps the read pair against the index and estimates how much of each transcript is present. Here Salmon turns raw reads into numbers to analyse.

3. Check read quality — FastQC. In parallel, FastQC inspects the raw reads of each sample (base quality, adapter content, and so on). We run it alongside quantification, not after, because quality control is independent of the counting.

4. Summarise everything — MultiQC. Finally, MultiQC gathers the logs from all the samples (Salmon and FastQC) into one HTML report, so you read a single summary instead of opening dozens of files. It waits for the all earlier steps to finish before it runs.

Note

Notice the two branches: quantification and quality control don’t depend on each other, so Nextflow runs them at the same time. They only come back together at MultiQC, which needs both before it can build the report. Keeping independent work parallel — and joining it only where necessary — is a large part of what makes a workflow manager worth using.

Preparations

If you haven’t done so yet, copy the directory “nf-training” from the common_data directory into your personal directory. Change directory into it.

Set the executor

So far we have used the local executor for our scripts, which is a big no no on clusters since that uses the log-in node for the execution. The pipeline we are building up now is a bit more complicated and requires more computing power than string manipulation, so it is time to set slurm as the executor. This can be done in a file called nextflow.config. I have prepared that file for you: locate the file nextflow.config in your folder and replace its content with the following:

process{
    executor = "slurm"
    time = '2.h'  // Set default time limit for all processes
    
    // Other settings
    cpus = 4
}

resume = true
singularity.enabled = true

executor {
    account = 'your_project_account'
}

• Here, we set the executor for every process to slurm.
• We set resume to true, so it is automatically used for all our runs (which means we do not have to specify this in our nextflow run command anymore).
• We enable singularity and set some singularity run options.
• Lastly, we specify the account name for the slurm execution, replace the placeholder with the actual compute allocation ID.

This nextflow.config file is implicitly called when we execute nextflow from the folder it is in. We do not have to explicitely name it in the nextflow run command.

Define the workflow parameters

Parameters are inputs and options that can be changed when the workflow is run.

The script script1.nf defines the workflow input parameters:

#!/usr/bin/env nextflow

params.reads = "$projectDir/data/ggal/gut_{1,2}.fq"
params.transcriptome_file = "$projectDir/data/ggal/transcriptome.fa"
params.multiqc = "$projectDir/multiqc"

workflow {
    println "reads: $params.reads"
}

Run it by using the following command:

pixi run nextflow run script1.nf

Note

Note that we do not need to specify that we want to use the nextflow.config file. Since it is in the same directory it will be used automatically.

ImportantTo do for you

Modify script1.nf by adding a fourth parameter named outdir and set it to a default path that will be used as the workflow output directory.

Create a transcriptome index file

Nextflow allows the execution of any command or script by using a process definition.

A process is defined by providing three main declarations: the process input, output, and the command script.

To add a transcriptome INDEX processing step, add the following code blocks to your script1.nf. Alternatively, these code blocks have already been added to script2.nf.

process INDEX {
    input:
    path transcriptome

    output:
    path 'salmon_index'

    script:
    """
    salmon index --threads $task.cpus -t $transcriptome -i salmon_index
    """
}

Additionally, replace the existing workflow scope with an input channel definition and the index process:

workflow {
    index_ch = INDEX(file(params.transcriptome_file, checkIfExists: true))
}

Here, the params.transcriptome_file parameter is used as the input for the INDEX process. The INDEX process (using the salmon tool) creates salmon_index, an indexed transcriptome that is passed as an output to the index_ch channel.

Run with:

pixi run nextflow run script1.nf 

Caution

Why doesn’t this work? What have we forgotten?

EITHER add Salmon to your environment

pixi add salmon

OR add Salmon in a container to the process

Adding a container to the process and running the tool via the container is more reproducible than adding it to your environment. Use the container directive and add the URL to the container to the process scope, above the input block:

container 'oras://community.wave.seqera.io/library/salmon:2.7.0--cb05356ce9eb88b7'

Note

You can also specify a process specific output directory in the process block:

publishDir "$params.outdir/salmon"

This command line would typically be located just below the container directive.

By default publishDir links the published files back into the work/ directory rather than copying them. This is fast and saves space, but it means the files disappear if you clean work/. Add mode: 'copy' when you want an independent copy that survives a cleanup, for example a result you intend to keep or share.

Define run time environment

Remember how Nextflow is big on execution abstraction? This is where we will put this concept into practice. In the nextflow.config add the following just before the closing squiggly bracket from the process definitions:

    withName:'INDEX'{
        time = 15.m
        cpus = 2
    }

With this we specify the alloted time and cpus for this specific process.

To specify memory add mem = xx.GB - check with your cluster to see how much memory is available and how to specify it. In our case we don’t need to because the memory implicit in the 2 cpus is enough for our purposes.

pixi run nextflow run script1.nf 

Interlude: Visualize tuples

Salmon needs both reads of a pair together to quantify a sample. Our data folder holds six files — three samples, two files each:

gut_1.fq    gut_2.fq
liver_1.fq  liver_2.fq
lung_1.fq   lung_2.fq

If we simply collected every .fq into one channel, Nextflow would see six loose files with nothing saying that gut_1 and gut_2 belong together. We need them grouped by sample.

That is what the fromFilePairs channel factory does, and we can have a look at it in script3.nf. The {1,2} in the pattern marks the part that differs between the two files of a pair; Nextflow groups the files that share everything else:

read_pairs_ch = channel.fromFilePairs(params.reads)

To see what it produces, view the channel from inside the workflow:

workflow {
    read_pairs_ch = channel.fromFilePairs(params.reads)
    read_pairs_ch.view()
}

Run it:

pixi run nextflow run script3.nf

It prints one item per sample:

[gut, [/.../data/ggal/gut_1.fq, /.../data/ggal/gut_2.fq]]

Read that line piece by piece — each item is a tuple with two parts:

• gut — the sample ID, taken from the shared prefix of the file names, and
• [gut_1.fq, gut_2.fq] — a list of the two files that belong to that sample.

Tip

Think of it like matching socks: fromFilePairs puts _1 and _2 together by their shared name, and the {1,2} glob tells Nextflow which part is allowed to differ.

Note

params.reads points at gut only, so you see one pair. Point it at all three samples to watch them grouped side by side:

pixi run nextflow run script3.nf --reads 'data/ggal/*_{1,2}.fq'

Because each item carries the ID and both files together, Salmon can process each sample in a single task — and Nextflow runs all three in parallel. That is the behaviour you will use in the next step.

Expression quantification

Look at the file “script4.nf”. This script adds a gene expression QUANTIFICATION process and a call to it within the workflow scope. Quantification requires the index transcriptome and RNA-Seq read pair fastq files.

In the workflow scope, note how the index_ch channel is assigned as output in the INDEX process.

Next, note that the first input channel for the QUANTIFICATION process is the previously declared index_ch, which contains the path to the salmon_index.

Also, note that the second input channel for the QUANTIFICATION process, is the read_pair_ch we just created. This being a tuple** composed of two elements (a value: “sample_id” and a list of paths to the fastq reads: “reads”) in order to match the structure of the items emitted by the fromFilePairs channel factory.

**tuple: a data structure consisting of multiple parts.

pixi run nextflow run script4.nf

Execute the same script again with more read files, as shown below:

pixi run nextflow run script4.nf --reads 'data/ggal/*_{1,2}.fq'

You will notice that the QUANTIFICATION process is executed multiple times.

Nextflow parallelizes the execution of your workflow simply by providing multiple sets of input data to your script.

ImportantTo do for you

Add the process specific runtime environment definition for QUANTIFICATION to nextflow.config.

Quality control

Now we want to add another process, using our beloved FastQC to check the samples. The input is the same as the read pairs used in the QUANTIFICATION step.

ImportantTo do for you

Sit together with your neighbour and first think 5 minutes what you need to do to implement this new process.

We will then discuss in class.

MultiQC

This last step collects the outputs from the QUANTIFICATION and FASTQC processes to create a final report using MultiQC.

Again, what do you need to consider when adding this process to the pipeline?

Tip

We only want one task of MultiQC to be executed to produce one report. For this, it needs to wait for the report generating processes to stop. We can use the mix channel operator to combine the two channels followed by the collect operator, to return the complete channel contents as a single element.

So the new line in the workflow scope looks like this:

MULTIQC(quant_ch.mix(fastqc_ch).collect())