eResearch nf-core-RNAseq pipeline

Owned by Roberto Barrero Gumiel
Last updated: Oct 16, 2023
5 min read

Aims

  • Run the nextflow nf-core/rnaseq pipeline in the HPC cluster. Exercises include:

    • Running a test to verify the execution of the pipeline

    • Running QC check to determine read trimming parameters

    • Running the full nf-core/rnaseq pipeline.

Work in the HPC

Work in the HPC

Before we start using the HPC, let’s start an interactive session:

qsub -I -S /bin/bash -l walltime=10:00:00 -l select=1:ncpus=1:mem=4gb

Get a copy of the scripts to be used in this module

Use the terminal to log into the HPC and create a /RNAseq/ folder to run the nf-core/rnaseq pipeline. For example:

mkdir -p $HOME/workshop/scripts cp /work/training/rnaseq/scripts/* $HOME/workshop/scripts/ ls -l $HOME/workshop/scripts/

  • Line 1: The -p indicates create 'parental directories as required. Thus the line 1 command creates both /workshop/ and the subfolder /workshop/scripts/

  • Line 2: Copies all files from /work/datasets/workshop/scripts/ as noted by an asterisk to the newly created folder $HOME/workshop/scripts/

Copy public data to your $HOME

mkdir -p $HOME/workshop/data cp /work/training/rnaseq/data/* $HOME/workshop/data/ # list the content of the $HOME/workshop/data/

  • Line 1: The first command creates the folder /scripts/

  • Line 2: Copies all files from /work/datasets/workshop/scripts/ folder as noted by an asterisk to newly created $HOME/workshop/scripts/ folder

  • Line 3: a quick challenge - see the previous section for hints

Create a folder for running the nf-RNA-seq pipeline

Let’s create an “RNAseq” folder to run the nf-core/rnaseq pipeline and move into it. For example:

  • Lines 1-4: create sub-folders for each exercise

  • Line 5: change the directory to the folder “run1_test”

  • Line 6: print the current working directory

Exercise 1: Running a test with nf-core sample data

First, let’s assess the execution of the nf-core/rnaseq pipeline by running a test using sample data.

Copy the launch_nf-core_RNAseq_test.pbs to the working directory

View the content of the script as follows:

#!/bin/bash -l

#PBS -N nfrnaseq_test

#PBS -l select=1:ncpus=2:mem=4gb

#PBS -l walltime=24:00:00

 

#work on current directory

cd $PBS_O_WORKDIR

 

#load java and set up memory settings to run nextflow

module load java

export NXF_OPTS='-Xms1g -Xmx4g'

 

nextflow run nf-core/rnaseq -r 3.12.0 -profile test,singularity --outdir results

#!/bin/bash -l

#PBS -N nfrnaseq_test

#PBS -l select=1:ncpus=2:mem=4gb

#PBS -l walltime=24:00:00

 

#work on current directory

cd $PBS_O_WORKDIR

 

#load java and set up memory settings to run nextflow

module load java

export NXF_OPTS='-Xms1g -Xmx4g'

 

nextflow run nf-core/rnaseq -r 3.12.0 -profile test,singularity --outdir results

  • nextflow command: nextflow run

  • pipeline name: nf-core/rnaseq

  • pipeline version: -r 3.12.0

  • container type and sample data: -profile test,singularity

  • output directory: --outdir results

Submitting the job

Submit the test job to the HPC cluster as follows:

Monitoring the Run

Exercise 2: Run RNA-seq QC check

The pipeline requires preparing at least 2 files:

  • Metadata file (samplesheet.csv) that specifies the name of the samples, location of FASTQ files ('Read 1' and ‘Read 2’), and strandedness (forward, reverse, or auto. Note: auto is used when the strandedness of the data is unknown)

  • PBS Pro script (launch_nf-core_RNAseq_QC.pbs) with instructions to run the pipeline

Create the metadata file (samplesheet.csv):

Change to the data folder directory:

Copy the bash script to the working folder

  • Note: you could replace ‘$HOME/workshop/data’ with “.” A dot indicates ‘current directory’ and will copy the file to the directory where you are currently located

View the content of the script:

Example for Single-End data (only ‘Read 1’ is available):

#!/bin/bash -l

 

#User defined variables

##########################################################

DIR='$HOME/workshop/data/'

INDEX='samplesheet.csv'

##########################################################

 

#load python module

module load python/3.10.8-gcccore-12.2.0

 

#fetch the script to create the sample metadata table

wget -L https://raw.githubusercontent.com/nf-core/rnaseq/master/bin/fastq_dir_to_samplesheet.py

chmod +x fastq_dir_to_samplesheet.py

 

#generate initial sample metadata file

./fastq_dir_to_samplesheet.py $DIR $INDEX \

        --strandedness auto \

        --read1_extension .fastq.gz

#!/bin/bash -l

 

#User defined variables

##########################################################

DIR='$HOME/workshop/data/'

INDEX='samplesheet.csv'

##########################################################

 

#load python module

module load python/3.10.8-gcccore-12.2.0

 

#fetch the script to create the sample metadata table

wget -L https://raw.githubusercontent.com/nf-core/rnaseq/master/bin/fastq_dir_to_samplesheet.py

chmod +x fastq_dir_to_samplesheet.py

 

#generate initial sample metadata file

./fastq_dir_to_samplesheet.py $DIR $INDEX \

        --strandedness auto \

        --read1_extension .fastq.gz

Let’s generate the metadata file by running the following command:

Check the newly created samplesheet.csv file:

 

sample,fastq_1,fastq_2,strandedness

CD49fmNGFRm_rep1,/work/eresearch_bio/training/data/rnaseq_data/mouse_PRJNA862107/SRR20622174_1.fastq.gz,,auto

CD49fmNGFRm_rep2,/work/eresearch_bio/training/data/rnaseq_data/mouse_PRJNA862107/SRR20622175_1.fastq.gz,,auto

CD49fmNGFRm_rep3,/work/eresearch_bio/training/data/rnaseq_data/mouse_PRJNA862107/SRR20622177_1.fastq.gz,,auto

CD49fpNGFRp_rep1,/work/eresearch_bio/training/data/rnaseq_data/mouse_PRJNA862107/SRR20622178_1.fastq.gz,,auto

CD49fpNGFRp_rep2,/work/eresearch_bio/training/data/rnaseq_data/mouse_PRJNA862107/SRR20622179_1.fastq.gz,,auto

CD49fpNGFRp_rep3,/work/eresearch_bio/training/data/rnaseq_data/mouse_PRJNA862107/SRR20622180_1.fastq.gz,,auto

MTEC_rep1,/work/eresearch_bio/training/data/rnaseq_data/mouse_PRJNA862107/SRR20622172_1.fastq.gz,,auto

MTEC_rep2,/work/eresearch_bio/training/data/rnaseq_data/mouse_PRJNA862107/SRR20622173_1.fastq.gz,,auto

MTEC_rep3,/work/eresearch_bio/training/data/rnaseq_data/mouse_PRJNA862107/SRR20622176_1.fastq.gz,,auto

sample,fastq_1,fastq_2,strandedness

CD49fmNGFRm_rep1,/work/eresearch_bio/training/data/rnaseq_data/mouse_PRJNA862107/SRR20622174_1.fastq.gz,,auto

CD49fmNGFRm_rep2,/work/eresearch_bio/training/data/rnaseq_data/mouse_PRJNA862107/SRR20622175_1.fastq.gz,,auto

CD49fmNGFRm_rep3,/work/eresearch_bio/training/data/rnaseq_data/mouse_PRJNA862107/SRR20622177_1.fastq.gz,,auto

CD49fpNGFRp_rep1,/work/eresearch_bio/training/data/rnaseq_data/mouse_PRJNA862107/SRR20622178_1.fastq.gz,,auto

CD49fpNGFRp_rep2,/work/eresearch_bio/training/data/rnaseq_data/mouse_PRJNA862107/SRR20622179_1.fastq.gz,,auto

CD49fpNGFRp_rep3,/work/eresearch_bio/training/data/rnaseq_data/mouse_PRJNA862107/SRR20622180_1.fastq.gz,,auto

MTEC_rep1,/work/eresearch_bio/training/data/rnaseq_data/mouse_PRJNA862107/SRR20622172_1.fastq.gz,,auto

MTEC_rep2,/work/eresearch_bio/training/data/rnaseq_data/mouse_PRJNA862107/SRR20622173_1.fastq.gz,,auto

MTEC_rep3,/work/eresearch_bio/training/data/rnaseq_data/mouse_PRJNA862107/SRR20622176_1.fastq.gz,,auto

Copy the PBS Pro script for QC (launch_nf-core_RNAseq_QC.pbs)

Copy and paste the code below to the terminal:

  • Line 1: Copy the samplesheet.csv file to the working directory

  • Line 2: move to the working directory

  • Line 3: copy the launch_nf-core_RNAseq_QC.pbs submission script to the working directory

View the content of the launch_nf-core_RNAseq_QC.pbs script:

#!/bin/bash -l

#PBS -N nfrnaseq_QC

#PBS -l select=1:ncpus=2:mem=4gb

#PBS -l walltime=24:00:00

 

#work on current directory

cd $PBS_O_WORKDIR

 

#load java and set up memory settings to run nextflow

module load java

export NXF_OPTS='-Xms1g -Xmx4g'

 

nextflow run nf-core/rnaseq \

      -profile singularity \

      -r 3.12.0 \

      --input samplesheet.csv \

      --outdir results \

      --genome GRCm38-local \

      --skip_trimming \

      --skip_alignment \

      --skip_pseudo_alignment

#!/bin/bash -l

#PBS -N nfrnaseq_QC

#PBS -l select=1:ncpus=2:mem=4gb

#PBS -l walltime=24:00:00

 

#work on current directory

cd $PBS_O_WORKDIR

 

#load java and set up memory settings to run nextflow

module load java

export NXF_OPTS='-Xms1g -Xmx4g'

 

nextflow run nf-core/rnaseq \

      -profile singularity \

      -r 3.12.0 \

      --input samplesheet.csv \

      --outdir results \

      --genome GRCm38-local \

      --skip_trimming \

      --skip_alignment \

      --skip_pseudo_alignment

  • We recommend running the nextflow nf-core/rnaseq pipeline once and then assessing the fastqc results folder to assess if sequence biases are present in the 5'-end and 3'-end ends of the sequences.

  • Version 3.12.0 allows running the pipeline to do quality assessment only, without any alignment, read counting, or trimming.

Submitting the job

Once you have created the folder for the run, the samplesheet.csv file, and launch.pbs, you are ready to submit the job to the HPC scheduler:

Monitoring the Run

to check on the jobs, you are running. Nextflow will launch additional jobs during the run.

You can also check the .nextflow.log file for details on what is going on.

Once the pipeline has finished running - Assess the QC report:

NOTE: To proceed, you need to be on QUT’s WiFi network or signed via VPN.

To browse the working folder in the HPC type in the file finder:

Windows PC

Mac

Evaluate the nucleotide distributions in the 5'-end and 3'-end of the sequenced reads (Read1 and Read2). Look into the “MultiQC” folder and open the provided HTML report.

Items to check:

  • The overall quality of the experiment and reads. Look at the “Sequence Quality Histogram” plot. For example, if Phred assigns a quality score of 30 to a base, the chances that this base is called incorrectly are 1 in 1000. Phred quality scores are logarithmically linked to error probabilities.

Phred Quality Score

Probability of incorrect base call

Base call accuracy

Phred Quality Score

Probability of incorrect base call

Base call accuracy

10

1 in 10

90%

20

1 in 100

99%

30

1 in 1000

99.9%

40

1 in 10,000

99.99%

50

1 in 100,000

99.999%

60

1 in 1,000,000

99.9999%

  • Assess QC reports (FastQC and MultiQC) to define how many nucleotides should be trimmed from the 5'-end and/or 3-end regions of the FASTQ reads (see Case 3 below).

Exercise 3: Run RNA-seq pipeline

Copy and paste the code below to the terminal:

  • Line 1: Copy the samplesheet.csv file to the working directory

  • Line 2: move to the working directory

  • Line 3: print working directory → verify the folder location

Copy the PBS Pro script to run the nf-core/rnaseq pipeline:

Adjusting the Trim Galore (read trimming) options

Print the content of the launch_RNAseq.pbs script:

#!/bin/bash -l

#PBS -N nfRNAseq

#PBS -l select=1:ncpus=2:mem=4gb

#PBS -l walltime=48:00:00

 

#work on current directory

cd $PBS_O_WORKDIR

 

#load java and set up memory settings to run nextflow

module load java

export NXF_OPTS='-Xms1g -Xmx4g'

 

nextflow run nf-core/rnaseq --input samplesheet.csv \

        --outdir results \

        -r 3.12.0 \

        --genome GRCm38-local \

        -profile singularity \

        --aligner star_salmon \

        --extra_trimgalore_args "--quality 30 --clip_r1 10 --clip_r2 10 --three_prime_clip_r1 1 --three_prime_clip_r2 1 "

#!/bin/bash -l

#PBS -N nfRNAseq

#PBS -l select=1:ncpus=2:mem=4gb

#PBS -l walltime=48:00:00

 

#work on current directory

cd $PBS_O_WORKDIR

 

#load java and set up memory settings to run nextflow

module load java

export NXF_OPTS='-Xms1g -Xmx4g'

 

nextflow run nf-core/rnaseq --input samplesheet.csv \

        --outdir results \

        -r 3.12.0 \

        --genome GRCm38-local \

        -profile singularity \

        --aligner star_salmon \

        --extra_trimgalore_args "--quality 30 --clip_r1 10 --clip_r2 10 --three_prime_clip_r1 1 --three_prime_clip_r2 1 "

Submitting the job

Monitoring the Run

Outputs

The pipeline will produce two folders, one called “work,” where all the processing is done, and another called “results,” where we can find the pipeline's outputs. The content of the results folder is as follows:

The quantification of the gene and transcript expressions can be found in the ‘star_salmon’ directory.

The following feature count tables are generated:

Tip: Read the help information for Nextflow pipelines

Information on how to run a nextflow pipeline and additional available parameters can be provided on the pipeline website (i.e., https://nf-co.re/rnaseq/3.12.0/docs/usage/ ). You can also run the following command to get help information:

Some pipelines may need file names, and others may want a CSV file with file names, the path to raw data files, and other information.