INTRO

This notebook is part of the coral E5 timeseries_molecular project (GitHub repo). Whole genome bisulfite sequencing (WGBS) data was received on 20241230 (Notebook entry) and then QC’d on 20241231. Based on initial QC, I decided to trim 25bp from the 5’ end of each read. Trimming was done with fastp, followed by FastQC and MultiQC.

Note

The contents below are from markdown knitted from 01.20-D-Apul-WGBS-trimming-fastp-FastQC-MultiQC.Rmd (commit 2072811).

1 Background

This notebook will trim raw WGBS FastQs using fastp. Samples which generate an error during trimming will attempt to be repaired using BBTools repairl.sh (BBMap – Bushnell B. – sourceforge.net/projects/bbmap/). Trimming results will be analyzed with FastQC and summarized by MultiQC (Ewels et al. 2016).

Based off of the initial FastQC/MultiQC, we trimmed 15bp from each read.

1.1 Inputs

Raw WGBS FastQ files with the following pattern:

*.fastq.gz

Note

If one needs to download the raw FastQs, please see 00.20-D-Apul-WGBS-reads-FastQC-MultiQC.Rmd

1.2 Outputs

The expected outputs will be:

*_fastqc.html: Individual FastQC reports.
*_fastqc.zip: Individual FastQC ZIP files.
*fastp-trim*.fq.gz: Trimmed FastQ files.
*.md5: Individual MD5 checksums for trimmed FastQs.
*.fastp-trim.report.html: Individual fastp trimming reports. HTML format.
*.fastp-trim.report.json: Individual fastp trimming reports. JSON format.
multiqc_report.html: A summary report of the alignment results generated by MultiQC, in HTML format.
- Due to the large file sizes of FastQs, these cannot be hosted in the timeseries_molecular GitHub repo. As such these files are available for download here:
  - https://gannet.fish.washington.edu/gitrepos/urol-e5/timeseries_molecular/D-Apul/output/01.20-D-Apul-WGBS-trimming-fastp-FastQC-MultiQC/

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 repo_dir=/home/shared/8TB_HDD_01/sam/gitrepos/urol-e5/timeseries_molecular'
echo 'export output_dir_top=${repo_dir}/D-Apul/output/01.20-D-Apul-WGBS-trimming-fastp-FastQC-MultiQC'
echo 'export raw_reads_dir="${repo_dir}/D-Apul/data/wgbs-raw-fastqs"'
echo 'export trimmed_fastqs_dir="${output_dir_top}/trimmed-fastqs"'
echo ""

echo "# Paths to programs"
echo 'export programs_dir="/home/shared"'
echo 'export bbmap_repair="${programs_dir}/bbmap-39.12/repair.sh"'
echo 'export bismark_dir="${programs_dir}/Bismark-0.24.0"'
echo 'export bowtie2_dir="${programs_dir}/bowtie2-2.4.4-linux-x86_64"'
echo 'export fastp="${programs_dir}/fastp-v0.24.0/fastp"'
echo 'export fastqc="${programs_dir}/FastQC-0.12.1/fastqc"'
echo 'export multiqc="/home/sam/programs/mambaforge/bin/multiqc"'
echo 'export samtools_dir="${programs_dir}/samtools-1.12"'
echo ""


echo "# Set FastQ filename patterns"
echo "export fastq_pattern='*.fastq.gz'"
echo "export R1_fastq_pattern='*_R1_*.fastq.gz'"
echo "export R2_fastq_pattern='*_R2_*.fastq.gz'"
echo "export trimmed_fastq_pattern='*fastp-trim*.fq.gz'"
echo ""

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


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

echo "# Print formatting"
echo 'export line="--------------------------------------------------------"'
echo ""
} > .bashvars

cat .bashvars

#### Assign Variables ####

# Data directories
export repo_dir=/home/shared/8TB_HDD_01/sam/gitrepos/urol-e5/timeseries_molecular
export output_dir_top=${repo_dir}/D-Apul/output/01.20-D-Apul-WGBS-trimming-fastp-FastQC-MultiQC
export raw_reads_dir="${repo_dir}/D-Apul/data/wgbs-raw-fastqs"
export trimmed_fastqs_dir="${output_dir_top}/trimmed-fastqs"

# Paths to programs
export programs_dir="/home/shared"
export bbmap_repair="${programs_dir}/bbmap-39.12/repair.sh"
export bismark_dir="${programs_dir}/Bismark-0.24.0"
export bowtie2_dir="${programs_dir}/bowtie2-2.4.4-linux-x86_64"
export fastp="${programs_dir}/fastp-v0.24.0/fastp"
export fastqc="${programs_dir}/FastQC-0.12.1/fastqc"
export multiqc="/home/sam/programs/mambaforge/bin/multiqc"
export samtools_dir="${programs_dir}/samtools-1.12"

# Set FastQ filename patterns
export fastq_pattern='*.fastq.gz'
export R1_fastq_pattern='*_R1_*.fastq.gz'
export R2_fastq_pattern='*_R2_*.fastq.gz'
export trimmed_fastq_pattern='*fastp-trim*.fq.gz'

# Set number of CPUs to use
export threads=40

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

# Print formatting
export line="--------------------------------------------------------"

3 Trimming and Error Repair

Trimming will remove any Illumina sequencing adapters, as well as polyG and polyX (primarily polyA) sequences. Trimming will also remove the first 25bp from the 5’ end of each read.

Samples generating an error during trimming will attempt to be repaired with BBTools’ repair.sh script.

# Load bash variables into memory
source .bashvars

# Make output directory, if it doesn't exist
mkdir --parents ${output_dir_top}

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

# Create array of fastq R2 files
# Do NOT quote R2_fastq_pattern variable
for fastq in ${R2_fastq_pattern}
do
  fastq_array_R2+=("${fastq}")

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



# Run fastp on files
# Adds JSON report output for downstream usage by MultiQC

echo "Beginning fastp trimming."
echo ""

for index in "${!fastq_array_R1[@]}"
do
  R1_sample_name="${R1_names_array[index]}"
  R2_sample_name="${R2_names_array[index]}"
  
  stderr_PE_name=$(echo ${R1_sample_name} | awk -F"_" '{print $1}')
  
  ${fastp} \
  --in1 ${fastq_array_R1[index]} \
  --in2 ${fastq_array_R2[index]} \
  --detect_adapter_for_pe \
  --trim_poly_g \
  --trim_poly_x \
  --thread 16 \
  --trim_front1 25 \
  --trim_front2 25 \
  --html ${output_dir_top}/"${R1_sample_name}".fastp-trim.report.html \
  --json ${output_dir_top}/"${R1_sample_name}".fastp-trim.report.json \
  --out1 ${output_dir_top}/"${R1_sample_name}".fastp-trim.fq.gz \
  --out2 ${output_dir_top}/"${R2_sample_name}".fastp-trim.fq.gz \
  2> ${output_dir_top}/"${stderr_PE_name}".fastp-trim.stderr
  
  grep --before-context 5 "ERROR" ${output_dir_top}/"${stderr_PE_name}".fastp-trim.stderr
  
  # Check for fastp errors and then repair
  if grep --quiet "ERROR" ${output_dir_top}/"${stderr_PE_name}".fastp-trim.stderr; then
      
    rm ${output_dir_top}/"${R1_sample_name}".fastp-trim.fq.gz
    rm ${output_dir_top}/"${R2_sample_name}".fastp-trim.fq.gz
    
    
    ${bbmap_repair} \
    in1=${fastq_array_R1[index]} \
    in2=${fastq_array_R2[index]} \
    out1="${R1_sample_name}".REPAIRED.fastq.gz \
    out2="${R2_sample_name}".REPAIRED.fastq.gz \
    outs=/dev/null \
    2> "${R1_sample_name}".REPAIRED.stderr
    
    ${fastp} \
    --in1 "${R1_sample_name}".REPAIRED.fastq.gz \
    --in2 "${R2_sample_name}".REPAIRED.fastq.gz \
    --detect_adapter_for_pe \
    --trim_poly_g \
    --trim_poly_x \
    --thread ${threads} \
    --trim_front1 25 \
    --trim_front2 25 \
    --html ${output_dir_top}/"${R1_sample_name}".fastp-trim.REPAIRED.report.html \
    --json ${output_dir_top}/"${R1_sample_name}".fastp-trim.REPAIRED.report.json \
    --out1 ${output_dir_top}/"${R1_sample_name}".fastp-trim.REPAIRED.fq.gz \
    --out2 ${output_dir_top}/"${R2_sample_name}".fastp-trim.REPAIRED.fq.gz \
    2> ${output_dir_top}/"${stderr_PE_name}".fastp-trim.REPAIRED.stderr
 
    
    if grep --quiet "ERROR" ${output_dir_top}/"${stderr_PE_name}".fastp-trim.REPAIRED.stderr; then
      echo "These ${stderr_PE_name} samples are broken."
      echo "Just give up.  :("
      echo ""
    fi
  fi

  echo "Finished trimming:"
  echo "${fastq_array_R1[index]}"
  echo "${fastq_array_R1[index]}"
  echo ""
  
  # Generate md5 checksums for newly trimmed files
  cd "${output_dir_top}"
  md5sum "${R1_sample_name}".fastp-trim.fq.gz > "${R1_sample_name}".fastp-trim.fq.gz.md5
  md5sum "${R2_sample_name}".fastp-trim.fq.gz > "${R2_sample_name}".fastp-trim.fq.gz.md5
  
  cd "${raw_reads_dir}"
done

4 FastQC

# Load bash variables into memory
source .bashvars

cd "${output_dir_top}"

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


# Create array of trimmed FastQs
trimmed_fastqs_array=(${trimmed_fastq_pattern})

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

echo "Beginning FastQC on trimmed reads..."
echo ""

# Run FastQC
### NOTE: Do NOT quote raw_fastqc_list
${fastqc} \
--threads ${threads} \
--outdir "${output_dir_top}" \
--quiet \
${trimmed_fastqc_list}

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

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

5 MultiQC

# Load bash variables into memory
source .bashvars

cd "${output_dir_top}"

${multiqc} \
--interactive \
.

RESULTS

Trimming worked, and all samples show < 5% adapter contamination (which is considered acceptable by FastQC). However, it is still annoying that any samples show adapter contamination. Additionally, there are a number of samples with relatively high levels of polyA contamination (approaching 5% of reads), despite specifying that fastp should trim polyA sequences… Regardless, the results are acceptable.

Output files

Output files are here:

https://gannet.fish.washington.edu/gitrepos/urol-e5/timeseries_molecular/D-Apul/output/01.20-D-Apul-WGBS-trimming-fastp-FastQC-MultiQC/
- MultiQC report (HTML)
  - multiqc_report.html

Ewels, Philip, Måns Magnusson, Sverker Lundin, and Max Käller. 2016. “MultiQC: Summarize Analysis Results for Multiple Tools and Samples in a Single Report.” Bioinformatics 32 (19): 3047–48. https://doi.org/10.1093/bioinformatics/btw354.