06.1-Peve-sRNAseq-trimming-R1-only

1 Introduction

FastQC/MultiQC (ewels2016?; Andrews_undated-nf?) assessment of raw and fastp-trimmed sequences of E5 P.evermanni sRNAseq data from 20230515.

This notebook will utilize R1 reads only, per Azenta’s sRNA-seq standard pipeline. The max length of 31bp is based on the fastp insert peak size from previous trimming tests based on the the adapter and polyG trimming results, and previous evaluation of mean read lengths via FastQC and MultiQC.

Inputs:

• sRNAseq R1 reads only gzipped FastQs (e.g. *R1_001.fastq.gz)

Outputs:

• FastQC HTML reports for raw and trimmed reads.

• MultiQC HTML summaries of FastQC for raw and trimmed reads.

• Trimmed reads with final length of 31bp: *.R1.fastp_trim.31bp.fastq.gz

Libraries were prepared and sequenced by Azenta:

• Library prep: NEB nebnext-small-rna-library-prep-set-for-illumina kit (PDF)

• Sequencing: Illumina HiSeq 4000, 150bp PE

---

2 Create a Bash variables file

This allows usage of Bash variables across R Markdown chunks.

{
echo "#### Assign Variables ####"
echo ""

echo "# Data directories"
echo 'export deep_dive_dir=/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive'
echo 'export output_dir_top=${deep_dive_dir}/E-Peve/output/06.1-Peve-sRNAseq-trimming-R1-only'
echo 'export raw_fastqc_dir=${deep_dive_dir}/E-Peve/output/06-Peve-sRNAseq-trimming/raw-fastqc'
echo 'export raw_reads_dir=${deep_dive_dir}/E-Peve/data/06-Peve-sRNAseq-trimming/raw-reads'
echo 'export raw_reads_url="https://owl.fish.washington.edu/nightingales/P_evermanni/30-852430235/"'
echo 'export trimmed_fastqc_dir=${output_dir_top}/trimmed-fastqc'
echo 'export trimmed_reads_dir=${output_dir_top}/trimmed-reads'
echo ""

echo "# Paths to programs"
echo 'export fastqc=/home/shared/FastQC-0.12.1/fastqc'
echo 'export fastp=/home/shared/fastp'
echo 'export multiqc=/home/sam/programs/mambaforge/bin/multiqc'
echo ""



echo "# Set FastQ filename patterns"
echo "export fastq_pattern='*.fq.gz'"
echo "export R1_fastq_pattern='sRNA*_R1_*.fastq.gz'"
echo ""

echo "# Set number of CPUs to use"
echo 'export threads=40'
echo ""

echo "# Set maximum read length"
echo 'export max_read_length=31'
echo ""

echo "# Input/output files"
echo 'export fastq_checksums=input_fastq_checksums.md5'
echo 'export trimmed_checksums=trimmed_fastq_checksums.md5'
echo 'export NEB_adapters_fasta=NEB-adapters.fasta'
echo ""

echo "## NEB nebnext-small-rna-library-prep-set-for-illumina adapters"
echo 'export first_adapter="AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC"'
echo 'export second_adapter="GATCGTCGGACTGTAGAACTCTGAACGTGTAGATCTCGGTGGTCGCCGTATCATT"'
echo ""

echo "## Inititalize arrays"
echo 'export fastq_array_R1=()'
echo 'export fastq_array_R2=()'
echo 'export raw_fastqs_array=()'
echo 'export R1_names_array=()'
echo 'export R2_names_array=()'
echo 'export trimmed_fastqs_array=()'
echo ""

echo "# Programs associative array"
echo "declare -A programs_array"
echo "programs_array=("
echo '[fastqc]="${fastqc}" \'
echo '[multiqc]="${multiqc}" \'
echo '[fastp]="${fastp}"'
echo ")"
} > .bashvars

cat .bashvars

#### Assign Variables ####

# Data directories
export deep_dive_dir=/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive
export output_dir_top=${deep_dive_dir}/E-Peve/output/06.1-Peve-sRNAseq-trimming-R1-only
export raw_fastqc_dir=${deep_dive_dir}/E-Peve/output/06-Peve-sRNAseq-trimming/raw-fastqc
export raw_reads_dir=${deep_dive_dir}/E-Peve/data/06-Peve-sRNAseq-trimming/raw-reads
export raw_reads_url="https://owl.fish.washington.edu/nightingales/P_evermanni/30-852430235/"
export trimmed_fastqc_dir=${output_dir_top}/trimmed-fastqc
export trimmed_reads_dir=${output_dir_top}/trimmed-reads

# Paths to programs
export fastqc=/home/shared/FastQC-0.12.1/fastqc
export fastp=/home/shared/fastp
export multiqc=/home/sam/programs/mambaforge/bin/multiqc

# Set FastQ filename patterns
export fastq_pattern='*.fq.gz'
export R1_fastq_pattern='sRNA*_R1_*.fastq.gz'

# Set number of CPUs to use
export threads=40

# Set maximum read length
export max_read_length=31

# Input/output files
export fastq_checksums=input_fastq_checksums.md5
export trimmed_checksums=trimmed_fastq_checksums.md5
export NEB_adapters_fasta=NEB-adapters.fasta

## NEB nebnext-small-rna-library-prep-set-for-illumina adapters
export first_adapter="AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC"
export second_adapter="GATCGTCGGACTGTAGAACTCTGAACGTGTAGATCTCGGTGGTCGCCGTATCATT"

## Inititalize arrays
export fastq_array_R1=()
export fastq_array_R2=()
export raw_fastqs_array=()
export R1_names_array=()
export R2_names_array=()
export trimmed_fastqs_array=()

# Programs associative array
declare -A programs_array
programs_array=(
[fastqc]="${fastqc}" \
[multiqc]="${multiqc}" \
[fastp]="${fastp}"
)

3 Download raw sRNAseq reads

Reads are downloaded from https://owl.fish.washington.edu/nightingales/P_evermanni/30-852430235/

The --cut-dirs 3 command cuts the preceding directory structure (i.e. nightingales/P_evermanni/30-852430235/) so that we just end up with the reads.

# Load bash variables into memory
source .bashvars

wget \
--directory-prefix ${raw_reads_dir} \
--recursive \
--no-check-certificate \
--continue \
--cut-dirs 3 \
--no-host-directories \
--no-parent \
--quiet \
--accept "${R1_fastq_pattern},checksums.md5" ${raw_reads_url}

ls -lh "${raw_reads_dir}"/${R1_fastq_pattern}

-rw-r--r-- 1 sam sam 953M May 17  2023 /home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/E-Peve/data/06-Peve-sRNAseq-trimming/raw-reads/sRNA-POR-73-S1-TP2_R1_001.fastq.gz
-rw-r--r-- 1 sam sam 899M May 17  2023 /home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/E-Peve/data/06-Peve-sRNAseq-trimming/raw-reads/sRNA-POR-79-S1-TP2_R1_001.fastq.gz
-rw-r--r-- 1 sam sam 934M May 17  2023 /home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/E-Peve/data/06-Peve-sRNAseq-trimming/raw-reads/sRNA-POR-82-S1-TP2_R1_001.fastq.gz

3.1 Verify raw read checksums

# Load bash variables into memory
source .bashvars

cd "${raw_reads_dir}"

# Checksums file contains other files, so this just looks for the sRNAseq files.
grep "sRNA" checksums.md5 | md5sum --check

sRNA-POR-73-S1-TP2_R1_001.fastq.gz: OK
sRNA-POR-73-S1-TP2_R2_001.fastq.gz: OK
sRNA-POR-79-S1-TP2_R1_001.fastq.gz: OK
sRNA-POR-79-S1-TP2_R2_001.fastq.gz: OK
sRNA-POR-82-S1-TP2_R1_001.fastq.gz: OK
sRNA-POR-82-S1-TP2_R2_001.fastq.gz: OK

4 Create adapters FastA for use with fastp trimming

# Load bash variables into memory
source .bashvars

echo "Creating adapters FastA."
echo ""
adapter_count=0

# Check for adapters file first
# Then create adapters file if doesn't exist
if [ -f "${output_dir_top}/${NEB_adapters_fasta}" ]; then
  echo "${output_dir_top}/${NEB_adapters_fasta} already exists. Nothing to do."
else
  for adapter in "${first_adapter}" "${second_adapter}"
  do
    adapter_count=$((adapter_count + 1))
    printf ">%s\n%s\n" "adapter_${adapter_count}" "${adapter}"
  done >> "${output_dir_top}/${NEB_adapters_fasta}"
fi

echo ""
echo "Adapters FastA:"
echo ""
cat "${output_dir_top}/${NEB_adapters_fasta}"
echo ""

Creating adapters FastA.

/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/E-Peve/output/06.1-Peve-sRNAseq-trimming-R1-only/NEB-adapters.fasta already exists. Nothing to do.

Adapters FastA:

>adapter_1
AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC
>adapter_2
GATCGTCGGACTGTAGAACTCTGAACGTGTAGATCTCGGTGGTCGCCGTATCATT

5 FastQC/MultiQC on raw reads

# Load bash variables into memory
source .bashvars

############ RUN FASTQC ############


# Create array of trimmed FastQs
raw_fastqs_array=(${raw_reads_dir}/${fastq_pattern})

# Pass array contents to new variable as space-delimited list
raw_fastqc_list=$(echo "${raw_fastqs_array[*]}")

echo "Beginning FastQC on raw reads..."
echo ""

# Run FastQC
### NOTE: Do NOT quote raw_fastqc_list
${programs_array[fastqc]} \
--threads ${threads} \
--outdir ${raw_fastqc_dir} \
--quiet \
${raw_fastqc_list}

echo "FastQC on raw reads complete!"
echo ""

############ END FASTQC ############

############ RUN MULTIQC ############
echo "Beginning MultiQC on raw FastQC..."
echo ""

${programs_array[multiqc]} ${raw_fastqc_dir} -o ${raw_fastqc_dir}

echo ""
echo "MultiQC on raw FastQs complete."
echo ""

############ END MULTIQC ############

echo "Removing FastQC zip files."
echo ""
rm ${raw_fastqc_dir}/*.zip
echo "FastQC zip files removed."
echo ""

# View directory contents
ls -lh ${raw_fastqc_dir}

6 Trimming with fastp

# Load bash variables into memory
source .bashvars

# Create output directory, if it doesn't exist.
mkdir --parents "${trimmed_reads_dir}"

# Change to directory with raw reads
cd "${raw_reads_dir}"

# Create arrays of FastQ R1 files and sample names
# Do NOT quote R1_fastq_pattern variable
for fastq in ${R1_fastq_pattern}
do
  fastq_array_R1+=("${fastq}")

  # Use parameter substitution to remove all text up to and including last "." from
  # right side of string.
  R1_names_array+=("${fastq%%.*}")
done



############ RUN FASTP ############
# Uses parameter substitution (e.g. ${R1_sample_name%%_*})to rm the _R[12]

# Run fastp on files
echo "Beginning fastp trimming."
echo ""

time \
for index in "${!fastq_array_R1[@]}"
do
  # Get sample name
  R1_sample_name="${R1_names_array[index]}"
  
  trimmed_fastp_name="${R1_sample_name}.fastp-R1-${max_read_length}bp-auto_adapters-polyG"


  # Begin fastp trimming
  ${programs_array[fastp]} \
  "${fastp}" \
  --in1 ${fastq_array_R1[index]} \
  --out1 ${trimmed_reads_dir}/${trimmed_fastp_name}.fq.gz \
  --trim_poly_g \
  --length_limit ${max_read_length} \
  --thread ${threads} \
  --html "${trimmed_fastp_name}.html" \
  --json "${trimmed_fastp_name}.json" \
  --report_title "${trimmed_fastp_name}"
    
  # Move to trimmed directory
  # This is done so checksums file doesn't include excess path
  cd ${trimmed_reads_dir}

  # Generate md5 checksums for newly trimmed files
  {
    md5sum "${trimmed_fastp_name}.fq.gz"
  } >> "${trimmed_checksums}"
    
  # Change back to to raw reads directory
  cd "${raw_reads_dir}"

done

echo ""
echo "fastp trimming complete."
echo ""

echo "Trimmed FastQs MD5 checksums:"
echo ""

cat "${trimmed_reads_dir}/${trimmed_checksums}"

############ END fastp ############

7 FastQC/MultiQC on trimmed reads

# Load bash variables into memory
source .bashvars

# Create output directory, if it doesn't exist.
mkdir --parents "${trimmed_fastqc_dir}"

############ RUN FASTQC ############

### NOTE: Do NOT quote raw_fastqc_list
# Create array of trimmed FastQs
trimmed_fastqs_array=(${trimmed_reads_dir}/${fastq_pattern})

# Pass array contents to new variable as space-delimited list
trimmed_fastqc_list=$(echo "${trimmed_fastqs_array[*]}")

echo "Beginning FastQC on raw reads..."
echo ""

# Run FastQC
${programs_array[fastqc]} \
--threads ${threads} \
--outdir ${trimmed_fastqc_dir} \
--quiet \
${trimmed_fastqc_list}

echo "FastQC on trimmed reads complete!"
echo ""

############ END FASTQC ############

############ RUN MULTIQC ############
echo "Beginning MultiQC on raw FastQC..."
echo ""

${programs_array[multiqc]} ${trimmed_fastqc_dir} -o ${trimmed_fastqc_dir}

echo ""
echo "MultiQC on trimmed FastQs complete."
echo ""

############ END MULTIQC ############

echo "Removing FastQC zip files."
echo ""
rm ${trimmed_fastqc_dir}/*.zip
echo "FastQC zip files removed."
echo ""

# View directory contents
ls -lh ${trimmed_fastqc_dir}

8 Summary

A quick comparison of raw and trimmed reads to show trimming worked:

• quality is improved
• length is 31bp
• adapters removed

RAW	TRIMMED

9 Citations