Overview of today’s session:
During Session 2 a basic generic RNA-seq pipeline has been run, without specifying additional parameters that can ensure the removal of sequence biases to have a more precise estimation of gene expression (feature counts). In this session, we will do the following tasks:
inspect the results from Session 2
run an advanced RNA-seq pipeline to measure the expression of genes
(optional) run statistical analysis to identify differentially expressed genes
Task 1: Evaluation of RNA-seq results using a basic (generic) nextflow pipeline
The nextflow/RNA-seq pipeline automatically generates two output folders:
results - contains the main outputs generated by each of the pipeline steps. Most users only need to look at the files contained in this folder.
work - this folder contains all the files generated during the running of the pipeline including intermediate files. Most users do not need to look at the content in this folder, unless the pipeline did not run properly.
To view the results from the completed pipeline, enter the run folder (i.e., run1_star_salmon)
#access the run folder for your samples. For example: cd run1_star_salmon #then access the results folder cd results
In the results folder you will find the following sub-folders:
fastqc/ trimgalore/ multiqc/ star_salmon/ pipeline_info/
FASTQC Report - assessing the quality of input reads
Connect to the work folder via HPC-FS (See session 2). Browse to the fastqc output folder: run1_star_salmon → results → fastqc. Then click on the HTML reports for each file to assess the quality of raw data. You may also copy the files to your laptop by simply drag-and-drop to a relevant folder.
The main items to verify are denoted below.
Per base sequence quality:
Inspect the overall quality of the generated data per nucleotide position.
Reads with a quality score above 20 (Q20) are 90.0% accurate, and those with >= Q30 are 99.9% accurate.
For most applications, it is recommended to set a quality trimming score of 30. Note, by default the pipeline will remove poor quality reads and bases below Q20.
Per base sequence content:
Determine if biases in the distribution of A, T, C, and G nucleotides are present on either the 5'-end and 3'-end of the reads
Recommendation: remove the first 10 nucleotides from the 5'-end (hexamer primer bias during PCR amplification) and 2 nucleotides from the 3'-end of reads (these bases can interfere with the proper mapping of reads onto reference genomes/transcriptomes).
Check other items reported in the FASTQC report such as level of duplication, highly abundant sequences, and presence of adapter sequences.
MultiQC Report - provides an overview of the quality, trimming, mapping, PCA, and many informative statistics of all files in the experiment in a single report.
Connect to the work folder via HPC-FS (See session 2). Browse to the fastqc output folder: run1_star_salmon → results → multiqc.
Task 2: Run the nextflow nf-core/rnaseq pipeline by including advanced filtering parameters
Requirements:
index.csv → a file that provides a list of sample IDs and their associated FASTQ files (read 1 and read 2)
launch.pbs → a script to submit the job to the HPC cluster
Example index.csv file for nf-core/rnaseq version 3.3:
group,fastq_1,fastq_2,strandedness control_r1,/work/kenna_team/data/raw_data/SRR1039508_1.fastq.gz,/work/kenna_team/data/raw_data/SRR1039508_2.fastq.gz,unstranded dex_r1,/work/kenna_team/data/raw_data/SRR1039509_1.fastq.gz,/work/kenna_team/data/raw_data/SRR1039509_2.fastq.gz,unstranded control_r2,/work/kenna_team/data/raw_data/SRR1039512_1.fastq.gz,/work/kenna_team/data/raw_data/SRR1039512_2.fastq.gz,unstranded dex_r2,/work/kenna_team/data/raw_data/SRR1039513_1.fastq.gz,/work/kenna_team/data/raw_data/SRR1039513_2.fastq.gz,unstranded control_r3,/work/kenna_team/data/raw_data/SRR1039516_1.fastq.gz,/work/kenna_team/data/raw_data/SRR1039516_2.fastq.gz,unstranded dex_r3,/work/kenna_team/data/raw_data/SRR1039517_1.fastq.gz,/work/kenna_team/data/raw_data/SRR1039517_2.fastq.gz,unstranded control_r4,/work/kenna_team/data/raw_data/SRR1039520_1.fastq.gz,/work/kenna_team/data/raw_data/SRR1039520_2.fastq.gz,unstranded dex_r4,/work/kenna_team/data/raw_data/SRR1039521_1.fastq.gz,/work/kenna_team/data/raw_data/SRR1039521_2.fastq.gz,unstranded
Example of a launch.pbs script with advanced parameter options:
#!/bin/bash -l #PBS -N nfrnaseq #PBS -l select=1:ncpus=2:mem=4gb #PBS -l walltime=24:00:00 #Use the current directory to run the workflow cd $PBS_O_WORKDIR module load java NXF_OPTS='-Xms1g -Xmx4g' #run the nextflow RNA-seq pipeline: nextflow run nf-core/rnaseq -profile singularity -r 3.3 --aligner star_salmon --input index.csv --genome GRCh38 --clip_r1 10 --clip_r2 10 --three_prime_clip_r1 2 --three_prime_clip_r2 2 --save_trimmed -resume
Session 3 exercises:
Run the nf-core/rnaseq pipeline using the Airway smooth muscle public data (PMID: 24926665. GEO: GSE52778) - aligner option set to ‘star_salmon’
Same as above but aligner option set to ‘star_rsem’
Create a new working folder:
mkdir session3 cd session3 mkdir run1_star_salmon cd run1_star_salmon cd ..
Copy index.csv and launch.pbs files to the newly created folder
cp /work/kenna_team/scripts/star_salmon/session3/* .
Check that files were copied into the new working folder
ls -a ./ ../ index.csv launch.pbs #verify the content of index.csv cat index.csv #also check the PBS Pro submission script cat launch.pbs
Run the workflow:
qsub launch.pbs
Monitor the progress of the workflow:
qjobs
or
qstats -u userID
Repeat the above process for ‘star_rsem’
The only variation is copying the index.csv and launch.pbs script. As follows:
cp /work/kenna_team/scripts/star_rsem/session3/* .