...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
Aim:
Assess the quality of raw datasets
Define quality trimming parameters prior RNAseq gene profiling
Prepare Working Directory Space for Session 4 |
|---|
Before we start using the HPC, let’s start an interactive session:
| Code Block |
|---|
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:
| Code Block |
|---|
mkdir -p $HOME/workshop/2024/rnaseq/scripts
cp /work/training/2024/rnaseq/scripts/* $HOME/workshop/2024/rnaseq/scripts/
ls -l $HOME/workshop/2024/rnaseq/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
| Code Block |
|---|
mkdir -p $HOME/workshop/2024/rnaseq/data
cp /work/training/2024/rnaseq/data/* $HOME/workshop/2024/rnaseq/data/
# list the content of the $HOME/workshop/2024/rnaseq/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 “runs” folder in the ~/workshop/2024/rnaseq folder to run the nf-core/rnaseq pipeline. For example:
| Code Block |
|---|
mkdir -p $HOME/workshop/2024/rnaseq/runs
mkdir $HOME/workshop/2024/rnaseq/runs/run1_test
mkdir $HOME/workshop/2024/rnaseq/runs/run2_QC
mkdir $HOME/workshop/2024/rnaseq/runs/run3_RNAseq
cd $HOME/workshop/2024/rnaseq/runs |
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
...
running the complete RNAseq gene profiling pipeline
Run RNA-seq QC check
The pipeline requires preparing at least 2 files:
Metadata file (samplesheet.csv) thatspecifies 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:
| Code Block |
|---|
cd $HOME/workshop/2024-2/rnaseqsession4_RNAseq/data/mouse |
Copy the bash script to the working folder
| Code Block |
|---|
cp /work/training/2024/rnaseq/scripts/create_samplesheet_nf-core_RNAseq_SEdata.sh $HOME/workshop/2024-2/rnaseqsession4_RNAseq/data/mouse |
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:
| Code Block |
|---|
cat create_samplesheet_nf-core_RNAseq_SEdata.sh |
Example for Single-End data (when 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
NOTE: when paired end data is available (R1 and R2 FASTQ files for each sample). Then use:
--read1_extension R1.fastq.gz \
--read2_extension R2.fastq.gz
...
NOTE: modify ‘read1_extension’ as appropriate for your data. For example: _1.fastq.gz or _R1_001.fastq.gz or _R1.fq.gz
...
, etc
Let’s generate the metadata file by running the following command:
| Code Block |
|---|
sh create_samplesheet_nf-core_RNAseq_samplesheetSEdata.sh $HOME/workshop/2024-2/session4_RNAseq/data/mouse |
Check the newly created samplesheet.csv file:
| Code Block |
|---|
ls -l cat samplesheet.cvs |
...
csv |
sample,fastq_1,fastq_2,strandedness SRR20622172,/work/training/rnaseq/data/SRR20622172.fastq.gz,,auto SRR20622173,/work/training/rnaseq/data/SRR20622173.fastq.gz,,auto SRR20622174,/work/training/rnaseq/data/SRR20622174.fastq.gz,,auto SRR20622175,/work/training/rnaseq/data/SRR20622175.fastq.gz,,auto SRR20622176,/work/training/rnaseq/data/SRR20622176.fastq.gz,,auto SRR20622177,/work/training/rnaseq/data/SRR20622177.fastq.gz,,auto SRR20622178,/work/training/rnaseq/data/SRR20622178.fastq.gz,,auto SRR20622179,/work/training/rnaseq/data/SRR20622179.fastq.gz,,auto SRR20622180,/work/training/rnaseq/data/SRR20622180.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:
| Code Block |
|---|
cp $HOME/workshop/2024-2/rnaseqsession4_RNAseq/data/mouse/samplesheet.csv $HOME/workshop/2024-2/rnaseqsession4_RNAseq/runs/run2run1_QC cp $HOME/workshop/2024-2/session4_RNAseq/scripts/launch_nf-core_RNAseq_QC.pbs $HOME/workshop/2024-2/rnaseqsession4_RNAseq/runs/run2run1_QC cd $HOME/workshop/2024-2/rnaseqsession4_RNAseq/runs/run2run1_QC |
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:
| Code Block |
|---|
cat launch_nf-core_RNAseq_QC.pbs |
#!/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.14.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:
...
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:
...
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% |
...
Inspect precomputed QC reports
FastQC:
location: /work/training/2024/rnaseq/runs/run2_QC/results/fastqc
Quality control metrics to assess - an individual html report is produced for each sample:
Per Base Sequence Quality: Measures the quality score across all bases at each position in the reads, providing a visual summary of how the quality varies across the length of the reads.
Per Sequence Quality Scores: Shows the distribution of mean quality scores for all reads, highlighting if the dataset contains many low-quality reads.
Per Base Sequence Content: Evaluates the proportion of each nucleotide (A, T, C, G) at every base position, revealing any biases in nucleotide composition
Per Base GC Content: Measures the percentage of GC content per base, checking for consistency with the expected GC content for the species or sample type.
Per Sequence GC Content: Examines the overall distribution of GC content across all reads, helping to detect unusual GC bias or contamination
Per Base N Content: Reports the proportion of uncalled bases (“N”s) at each position in the reads, indicating sequencing or base-calling issues
Sequence Length Distribution: Displays the distribution of read lengths, which can help assess whether the sequencing run generated reads of the expected length.
Sequence Duplication Levels: Identifies the level of duplicate sequences in the dataset, which can indicate issues such as overamplification in PCR or technical biases.
Overrepresented Sequences: Detects sequences that appear more frequently than expected, which might indicate contamination, adapters, or other technical artifacts.
Adapter Content: Checks for the presence of sequencing adapters that were not properly removed during preprocessing, which can affect downstream analysis.
K-mer Content: Identifies enriched k-mers (short nucleotide sequences) that occur more frequently than expected, possibly indicating contamination or other biases.
MultiQC:
Combined FastQC results for all samples:
Per Base Sequence Quality: Overview of quality scores across bases.
Per Sequence Quality Scores: Mean quality scores of reads.
Per Base Sequence Content: Distribution of nucleotide composition across reads.
Per Sequence GC Content: GC content distribution.
Sequence Length Distribution: The length of reads in the data.
Adapter Content: Detection of sequencing adapters.
Hit: define how many nucleotides should be trimmed from the 5'-end and/or 3-end regions of the FASTQ reads
...
.