Versions Compared

Old Version 22

changes.mady.by.user Roberto Barrero Gumiel

Saved on

compared with

New Version Current

changes.mady.by.user Roberto Barrero Gumiel

Saved on

Key

  • This line was added.
  • This line was removed.
  • Formatting was changed.

Prior running the nf-core/sarek pipeline with real data, we will first run a test with sample data to make sure the pipeline runs properly.

...

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

You should be in your home directory, if unsure you can run the following command:

Code Block
cd ~

List the existing files and folders:

Code Block
ls -l

Let’s create a folder for the workshop

Code Block
mkdir -p $HOME/workshop/sarek

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

Use the terminal to log into the HPC and Now let’s create a /RNAseq/ folder to run the nf-core/rnaseq pipeline. For example‘scripts’ folder and copy all scripts that we will using in the session:

Code Block
mkdir -p $HOME/workshop/sarek/scripts
cp /work/training/sarek/scripts/* $HOME/workshop/sarek/scripts/
ls -l $HOME/workshop/sarek/scripts/

  • Line 1: The -p indicates create 'parental directories as required. Thus the line 1 command creates both /workshop/ and the subfolder /workshop/sarek/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

Code Block
mkdir -p $HOME/workshop/sarek/data/WES/trio
mkdir -p $HOME/workshop/sarek/data/WES/liver
cp /work/training/sarek/data/WES/trio/* $HOME/workshop/sarek/data/WES/trio
cp /work/training/sarek/data/WES/liver/* $HOME/workshop/sarek/data/WES/liver

  • Lines 1 -2: Command creates the folders to copy data

  • Line 3: Copies all files from /work/datasets/workshop/sarek/data/WES/trio folder as noted by an asterisk to newly created $HOME/workshop/sarek/data/WES/trio folder.

  • Line 4: Copies all files from /work/datasets/workshop/sarek/data/WES/liver folder as noted by an asterisk to newly created $HOME/workshop/sarek/data/WES/liver folder.sarek/scripts/

Create folders for running the nf-core/sarek pipeline

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

Code Block
mkdir -p $HOME/workshop/sarek
mkdir $HOME/workshop/sarek/run1_test
mkdir $HOME/workshop/sarek/run2_trio
mkdir $HOME/workshop/sarek/run3_liver
cd $HOME/workshop/sarek

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

  • Line 54: 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_RNAseqsarek_test.pbs to the working directory

...

Code Block
cat launch_nf-core_sarek_test.pbs

#!/bin/bash -l

#PBS -N nfsarek_run1_test

#PBS -l walltime=48:00:00

#PBS -l select=1:ncpus=1:mem=5gb

cd $PBS_O_WORKDIR

NXF_OPTS='-Xms1g -Xmx4g'

module load java

#specify the nextflow version to use to run the workflow

export NXF_VER=23.10.1

#run the sarek pipeline

nextflow run nf-core/sarek \

        -r 3.3.2 \

        -profile test,singularity \

        --outdir ./results

  • nextflow command: nextflow run

  • pipeline name: nf-core/sarek

  • pipeline version: -r 3.3.2

  • 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:

...

Code Block
results/
├── csv
│   ├── markduplicates.csv
│   ├── markduplicates_no_table.csv
│   ├── recalibrated.csv
│   └── variantcalled.csv
├── multiqc
│   ├── multiqc_data
│   ├── multiqc_plots
│   └── multiqc_report.html
├── pipeline_info
│   ├── execution_report_2024-05-08_15-28-38.html
│   ├── execution_timeline_2024-05-08_15-28-38.html
│   ├── execution_trace_2024-05-08_15-28-38.txt
│   ├── params_2024-05-08_15-41-30.json
│   ├── pipeline_dag_2024-05-08_15-28-38.html
│   └── software_versions.yml
├── preprocessing
│   ├── markduplicates
│   ├── recalibrated
│   └── recal_table
├── reports
│   ├── bcftools
│   ├── fastqc
│   ├── markduplicates
│   ├── mosdepth
│   ├── samtools
│   └── vcftools
├── tabix
│   ├── genome.bed.gz
│   └── genome.bed.gz.tbi
└── variant_calling
    └── strelka

Exercise 2: Run nf-core/sarek using trio data

The pipeline requires preparing at least 2 files:

...

 Code Block
patient,sample,lane,fastq_1,fastq_2
ID1,S1,L002,/full/path/to/ID1_S1_L002_R1_001.fastq.gz,/full/path/to/ID1_S1_L002_R2_001.fastq.gz
  • PBS Pro script (launch_nf-core_sarek_trio.pbs) with instructions to run the pipeline
Create the metadata file (samplesheet.csv):
Change to the data folder directory:
 Code Block
cd $HOME/workshop/sarek/data/trio
pwd
Copy the python script “create_samplesheet_nf-core_sarek.py" to the working folder
 Code Block
cp /work/training/sarek/scripts/create_samplesheet_nf-core_sarek.py $HOME/workshop/sarek/data/trio
  • Note: you could replace ‘$HOME/workshop/sarek/data’ with “.” A dot indicates ‘current directory’ and will copy the file to the directory where you are currently located
Check help option on how to run the script:
 Code Block
python create_samplesheet_nf-core_sarek.py --help
 Code Block
python create_samplesheet_nf-core_sarek.py -h
usage: create_samplesheet_nf-core_sarek.py [-h] [--dir DIR] [--read1_extension READ1_EXTENSION] [--read2_extension READ2_EXTENSION] [--out OUT]
Extract metadata from fastq files in a directory.
optional arguments:
  -h, --help            show this help message and exit
  --dir DIR             Directory to search for files (default: current directory)
  --read1_extension READ1_EXTENSION
                        Extension for fastq_1 files (default: R1_001.fastq.gz)
  --read2_extension READ2_EXTENSION
                        Extension for fastq_2 files (default: R2_001.fastq.gz)
  --out OUT             Output metadata CSV file
Let’s generate the metadata file by running the following command:
 Code Block
python create_samplesheet_nf-core_sarek.py --dir $HOME/workshop/sarek/data/trio \
  --read1_extension R1.fastq.gz \
  --read2_extension R2.fastq.gz \
  --out samplesheet.csv
Check the newly created samplesheet.csv file:
 Code Block
ls -l
cat samplesheet.cvs
patient,sample,lane,fastq_1,fastq_2
SRR14724455,NA12892a,L001,/work/training/sarek/data/WES/trio/SRR14724455_NA12892a_L001_R1.fastq.gz,/work/training/sarek/data/WES/trio/SRR14724455_NA12892a_L001_R2.fastq.gz
SRR14724456,NA12891a,L001,/work/training/sarek/data/WES/trio/SRR14724456_NA12891a_L001_R1.fastq.gz,/work/training/sarek/data/WES/trio/SRR14724456_NA12891a_L001_R2.fastq.gz
SRR14724463,NA12878a,L001,/work/training/sarek/data/WES/trio/SRR14724463_NA12878a_L001_R1.fastq.gz,/work/training/sarek/data/WES/trio/SRR14724463_NA12878a_L001_R2.fastq.gz
SRR14724474,NA12892b,L001,/work/training/sarek/data/WES/trio/SRR14724474_NA12892b_L001_R1.fastq.gz,/work/training/sarek/data/WES/trio/SRR14724474_NA12892b_L001_R2.fastq.gz
SRR14724475,NA12891b,L001,/work/training/sarek/data/WES/trio/SRR14724475_NA12891b_L001_R1.fastq.gz,/work/training/sarek/data/WES/trio/SRR14724475_NA12891b_L001_R2.fastq.gz
SRR14724483,NA12878b,L001,/work/training/sarek/data/WES/trio/SRR14724483_NA12878b_L001_R1.fastq.gz,/work/training/sarek/data/WES/trio/SRR14724483_NA12878b_L001_R2.fastq.gz
Copy the PBS Pro script for running the nf-core/sarek pipeline (launch_nf-core_sarek_trio.pbs) 
Copy and paste the code below to the terminal:
 Code Block
cp $HOME/workshop/sarek/data/WES/trio/samplesheet.csv $HOME/workshop/sarek/runs/run2_sarek_trio
cp $HOME/workshop/sarek/scripts/launch_nf-core_sarek_trio.pbs $HOME/workshop/sarek/runs/run2_trio
cd $HOME/workshop/sarek/runs/run2_trio
  • Line 1: Copy the samplesheet.csv file generated above to the working directory
  • Line 2: copy the launch_nf-core_sarek_trio.pbs submission script to the working directory
  • Line 3: move to the working directory
View the content of the launch_nf-core_RNAseq_QC.pbs script:
 Code Block
cat launch_nf-core_RNAseq_QC.pbs
#!/bin/bash -l
#PBS -N nfsarek_run2_trio
#PBS -l walltime=48:00:00
#PBS -l select=1:ncpus=1:mem=5gb
cd $PBS_O_WORKDIR
NXF_OPTS='-Xms1g -Xmx4g'
module load java
#specify the nextflow version to use to run the workflow
export NXF_VER=23.10.1
#run the sarek pipeline
nextflow run nf-core/sarek \
        -r 3.3.2 \
        -profile singularity \
        --genome GATK.GRCh38 \
        --input samplesheet.csv \
        --wes \
        --outdir ./results \
        --step mapping \
        --tools haplotypecaller,snpeff,vep \
        --snpeff_cache /work/training/sarek/NXF_SINGULARITY_CACHEDIR/snpeff_cache \
        --vep_cache /work/training/sarek/NXF_SINGULARITY_CACHEDIR/vep_cache \
        -resume
  • The above script will screen for germline (inherited) mutations using GATK’s haplotypecaller and then annotate the identified variants using snpeff and VEP. 
  • Version 3.3.2 allows running the pipeline to do quality assessment only, without any alignment, read counting, or trimming.
  • The pipeline enables use to start at distinct stages, we are commencing from the start “--step mapping”
Submitting the job
Once 
...
 the 
...
 Code Block
qsub launch_nf-core_sarek_trio.pbs
Monitoring the Run 
 Code Block
qjobs
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.
...
 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:
...
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:
...
.
...
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: Healthy vs. disease liver samples
Copy and paste the code below to the terminal:
 Code Block
cp $HOME/workshop/data/samplesheet.csv $HOME/workshop/RNAseq/run3_RNAseq
cd $HOME/workshop/RNAseq/run3_RNAseq
pwd
  • 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:
 Code Block
cp $HOME/workshop/scripts/launch_nf-core_RNAseq_pipeline.pbs $HOME/workshop/RNAseq/run3_RNAseq
Adjusting the Trim Galore (read trimming) options
Print the content of the launch_RNAseq.pbs script:
 Code Block
cat launch_nf-core_RNAseq_pipeline.pbs
#!/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
 Code Block
qsub launch_nf-core_RNAseq_pipeline.pbs
Monitoring the Run
 Code Block
qjobs
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: