My report for analysis X (project title)

Nina Dombrowski

2023-12-18

Description

The goal of this workflow is to analyse sequencing data received on 23.01.2023 from Data Science Sequencing Centre. The sequencing data consists of 4 samples, with a forward (R1) and a reverse (R2) read file for each sample. The code below documents the quality checks and the cleaning of the sequences.

Software

All analyses were run on the UvA Crunchomics HPC with a CentOS 7 operating system.

Dependencies:

• SeqKit v2.7.0 (pre-installed)
• fastp v1.4.0 (conda environment fastp_1.4.0)
• …

Installation notes for software that was not pre-installed:

mamba create -n fastp_1.4.0 -c bioconda fastp=1.4.0

Analysis workflow

Settings (change these before you run the workflow)

# Path to working directory
wdir="$HOME/personal/data_analysis"

# Custom variables
download_link="https://raw.githubusercontent.com/ndombrowski/cli_workshop/main/data/seq_project.tar.gz"

Getting the data

# Go to project dir
cd $wdir 

# Setup base folders
mkdir data scripts logs results

# Download data using link provided by the sequencing center on 23.01.2023
wget -P data $download_link
tar -xvf data/seq_project.tar.gz -C data
rm data/*.tar.gz

# Generate a list for the files we work with 
ls data/seq_project/*/*.gz | cut -f4 -d "/" > fastq_files.txt

# Generate a list of sample names (used in scripts/fastp.sh)
ls data/seq_project/*/*_R1.fastq.gz | cut -f4 -d "/" | cut -f1 -d "_" > samples.txt

Looking at the data

# We work with 8 files from 4 samples
wc -l fastq_files.txt
wc -l samples.txt

# We have as many R1 as we have R2 files
ls data/seq_project/*/*R1* | wc -l
ls data/seq_project/*/*R2* | wc -l

# The sequence header looks as expected
zcat data/seq_project/*/*gz | head

# R1 and R2 reads have the same number of reads
mkdir -p results/counts

for fastq_file in $(cat fastq_files.txt); do
    echo "$fastq_file: $(( $(zcat data/seq_project/*/$fastq_file | wc -l) /4 ))"
done > results/counts/read_counts.txt

# We expect 8 lines, one per fastq file
wc -l results/counts/read_counts.txt
cat results/counts/read_counts.txt

We work with so many reads per fastq file:

Sample-DUMMY1_R1.fastq.gz: 44846
Sample-DUMMY1_R2.fastq.gz: 44846
Sample-DUMMY2_R1.fastq.gz: 19733
Sample-DUMMY2_R2.fastq.gz: 19733
Sample-DUMMY3_R1.fastq.gz: 19784
Sample-DUMMY3_R2.fastq.gz: 19784
Sample-DUMMY4_R1.fastq.gz: 24380
Sample-DUMMY4_R2.fastq.gz: 24380

Read statistics with seqkit

Run seqkit stats

mkdir -p results/seqkit

# Submitted batch job 380625
sbatch scripts/seqkit.sh

# We expect one table with a header line and one line per fastq file (9 lines)
ls -l results/seqkit
wc -l results/seqkit/stats_raw.txt
cat results/seqkit/stats_raw.txt

# Get some basic summary statistics with awk:
# -F'\t': the columns in the table are separated by tabs
# NR>1: skip the first line (the header)
# reads+=$4; bp+=$5: add up column 4 (num_seqs) and column 5 (sum_len); n++ counts the files
# END {...}: after the last line, print the totals and the mean per fastq file
awk -F'\t' 'NR>1 {reads+=$4; bp+=$5; n++}
    END {print "files:", n; print "reads:", reads; print "bp:", bp; printf "mean reads per file: %.1f\n", reads/n; printf "mean bp per file: %.1f\n", bp/n}' results/seqkit/stats_raw.txt

Comments on the read statistics from seqkit stats:

• We work with a total of
  • 217,486 reads (avg 27,186 read pairs per sample)
  • 54,493,987 bp (avg 6,811,748 bp per fastq file)
• The number of reads (num_seqs) matches the read counts from the sanity checks above, and R1 and R2 have the same number of reads for each sample
• The reads are typically 251 bp long (the average length is 250-251 bp), but each file also contains some shorter reads (shortest read: 35-239 bp). This suggests that there was some pre-trimming done by the sequencing center
• The R1 reads have a higher quality than the R2 reads: 94-95% of the bases in R1 have a quality score of at least 30 (Q30(%)), compared to 89-90% in R2
• The GC content is 53-57%

Quality filtering with fastp

Run fastp

mkdir -p results/fastp

# Submitted batch job 380626
sbatch scripts/fastp.sh

# We expect 4 samples x 4 files (2 trimmed fastq.gz, 1 html, 1 json) = 16 files
ls results/fastp | wc -l

Run seqkit stats on the trimmed reads

# Submitted batch job 380627
sbatch scripts/seqkit_trimmed.sh

# We expect one table with a header line and one line per trimmed fastq file (9 lines)
ls -l results/seqkit
wc -l results/seqkit/stats_trimmed.txt

# Same summary statistics as for the raw reads
awk -F'\t' 'NR>1 {reads+=$4; bp+=$5; n++}
    END {print "files:", n; print "reads:", reads; print "bp:", bp; printf "mean reads per file: %.1f\n", reads/n; printf "mean bp per file: %.1f\n", bp/n}' results/seqkit/stats_trimmed.txt

Comments on the seqkit output post trimming:

• After cleaning we work with
  • 214,906 out of 217,486 reads (98.8% maintained)
  • 53,826,288 out of 54,493,987 bp (98.8% maintained)
• Information from the fastp html:
  • Little adapter sequence signal (< 0.05%)
  • Read quality drops at the beginning and end of the read –> reconsider re-running fastp and use --cut_front, --cut_front_window_size 4 (default setting) and --cut_front_mean_quality 20 (default setting). Note that --cut_front is disabled by default

Citation

If your work gets published and you provide the code for example on Github it’s useful to add information about how to cite your work

Contact

For questions, please contact nd@sds.net