Versions Compared

Key

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

Exercise 3: Run nf-core/sarek using a liver samples 

What is NAFLD?

Non-alcoholic fatty liver disease (NAFLD) is a condition characterized by an accumulation of fat in liver cells (hepatocytes). Excess fat in the liver can lead to significant damage over the years. There are two types of NAFLD:

...

NASH (Non-alcoholic steatohepatitis) – It is a malignant condition and in this fat is accumulated in the liver can cause scarring of liver, liver fibrosis and liver damage further it can cause cirrhosis of liver (also called as end-stage liver disease), further cirrhosis of the liver may lead to liver cancer.

...

 

The pipeline requires preparing at least 2 files:

...

Location of raw data:

Code Block
/work/training/sarek/data/WES/liver
Code Block
├── liver
│   ├── Control1_C1_L001_R1_001.fastq.gz
│   ├── Control1_C1_L001_R2_001.fastq.gz
│   ├── Control2_C2_L001_R1_001.fastq.gz
│   ├── Control2_C2_L001_R2_001.fastq.gz
│   ├── Control3_C3_L001_R1_001.fastq.gz
│   ├── Control3_C3_L001_R2_001.fastq.gz
│   ├── Control4_C4_L001_R1_001.fastq.gz
│   ├── Control4_C4_L001_R2_001.fastq.gz
│   ├── NAFLD1_P1_L001_R1_001.fastq.gz
│   ├── NAFLD1_P1_L001_R2_001.fastq.gz
│   ├── NAFLD2_P2_L001_R1_001.fastq.gz
│   ├── NAFLD2_P2_L001_R2_001.fastq.gz
│   ├── NAFLD3_P3_L001_R1_001.fastq.gz,/full/path/to/ID1_S1_L002
│   ├── NAFLD3_P3_L001_R2_001.fastq.gz
  • PBS Pro script (launch_nf-core_sarek_trio.pbs) with instructions to run the pipeline


│   ├── NAFLD4_P4_L001_R1_001.fastq.gz
│   ├── NAFLD4_P4_L001_R2_001.fastq.gz
│   └── samplesheet.csv

STEP1: Create the metadata file (samplesheet.csv):

Change to the data folder directory:

Code Block
cd $HOME/workshop/sarek/dataruns/triorun3_liver
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/triorun3_liver
  • Note: you could replace ‘$HOME/workshop/sarek/runs/data’ liver’ 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 -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/trioliver \
  --read1_extension R1.fastq.gz \
  --read2_extension R2.fastq.gz \
  --out samplesheet.csv

...

Code Block
ls -l
cat samplesheet.cvs

patient,sample,lane,fastq_1,fastq_2

SRR14724455,NA12892a

Control1,C1,L001,/sarek/data/WES/liver/Control1_C1_L001_R1_001.fastq.gz,/sarek/data/WES/liver/Control1_C1_L001_R2_001.fastq.gz

Control2,C2,L001,/sarek/data/WES/

trio

liver/

SRR14724455

Control2_

NA12892a

C2_L001_R1_001.fastq.gz,/sarek/data/WES/

trio

liver/

SRR14724455

Control2_

NA12892a

C2_L001_R2_001.fastq.gz

SRR14724456

Control3,

NA12891a

C3,L001,/sarek/data/WES/

trio

liver/

SRR14724456

Control3_

NA12891a

C3_L001_R1_001.fastq.gz,/sarek/data/WES/

trio

liver/

SRR14724456

Control3_

NA12891a

C3_L001_R2_001.fastq.gz

SRR14724463

Control4,

NA12878a

C4,L001,/sarek/data/WES/

trio

liver/

SRR14724463

Control4_

NA12878a

C4_L001_R1_001.fastq.gz,/sarek/data/WES/

trio

liver/

SRR14724463

Control4_

NA12878a

C4_L001_R2_001.fastq.gz

SRR14724474

NAFLD1,

NA12892b

P1,L001,/sarek/data/WES/

trio

liver/

SRR14724474

NAFLD1_

NA12892b

P1_L001_R1_001.fastq.gz,/sarek/data/WES/

trio

liver/

SRR14724474

NAFLD1_

NA12892b

P1_L001_R2_001.fastq.gz

SRR14724475

NAFLD2,

NA12891b

P2,L001,/sarek/data/WES/

trio

liver/

SRR14724475

NAFLD2_

NA12891b

P2_L001_R1_001.fastq.gz,/sarek/data/WES/

trio

liver/

SRR14724475

NAFLD2_

NA12891b

P2_L001_R2_001.fastq.gz

SRR14724483

NAFLD3,

NA12878b

P3,L001,/sarek/data/WES/

trio

liver/

SRR14724483

NAFLD3_

NA12878b

P3_L001_R1_001.fastq.gz,/sarek/data/WES/

trio

liver/

SRR14724483

NAFLD3_

NA12878b

P3_L001_R2_001.fastq.gz

Copy the PBS Pro script for running the nf-core/sarek pipeline (launch_nf-core_sarek_trio.pbs)

...

NAFLD4,P4,L001,/sarek/data/WES/liver/NAFLD4_P4_L001_R1_001.fastq.gz,/sarek/data/WES/liver/NAFLD4_P4_L001_R2_001.fastq.gz

Alternatively copy the samplesheet.csv file:

Code Block
cp $HOME/work/workshoptraining/sarek/data/WES/trioliver/samplesheet.csv $HOME/workshop/sarek/runs/run2_sarek_trio
.

STEP2 - Run the nf-core/sarek pipeline

Copy and paste the code below to the terminal:

Code Block
cp $HOME/workshop/sarek/scripts/launch_nf-core_sarek_trioliver.pbs $HOME/workshop/sarek/runs/run2run3_trioliver
cd $HOME/workshop/sarek/runs/run2run3_trioliver
  • 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_RNAseqsarek_QCliver.pbs

#!/bin/bash -l

#PBS -N nfsarek_

run2_trio

liver

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

...

Once the pipeline has finished running - Assess the results as follows:

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.

During execution of the workflow two output folders are generated:

  • work - where all intermediate results and tasks are run

  • results - where all final results for all stages of the pipeline are copied

Let’s browse the results of the pipeline:

Code Block
results/
├── annotation
│   └── haplotypecaller
├── 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-11_20-17-13.html
│   ├── execution_timeline_2024-05-11_20-17-13.html
│   ├── execution_trace_2024-05-11_20-17-13.txt
│   ├── params_2024-05-11_23-31-28.json
│   ├── pipeline_dag_2024-05-11_20-17-13.html
│   └── software_versions.yml
├── preprocessing
│   ├── markduplicates
│   ├── recalibrated
│   └── recal_table
├── reports
│   ├── bcftools
│   ├── EnsemblVEP
│   ├── fastp
│   ├── fastqc
│   ├── markduplicates
│   ├── mosdepth
│   ├── samtools
│   ├── snpeff
│   └── vcftools
├── tabix
│   ├── wgs_calling_regions_noseconds.hg38.bed.gz
│   └── wgs_calling_regions_noseconds.hg38.bed.gz.tbi
└── variant_calling
    └── haplotypecaller