02.00-bismark-bowtie2-alignment

1 Intro

This notebook will align trimmed/repaired FastQs generated by 00.00-trimming-fastp.Rmd to the C.viginica genome using Bismark (Krueger and Andrews 2011) and Bowtie2 (Langmead et al. 2018; Langmead and Salzberg 2012). Alignment results will be summarized by MultiQC (Ewels et al. 2016).

Alignments were run with the Bowtie score setting of L,0,-0.6.

NOTE: This process will take ~2 weeks to complete, when run on a computer with 48 CPUs!

NOTE: This notebook was not used for this project. It was superseded by 02.01-bismark-bowtie2-alignment-SLURM-array.sh, which was executed as part of an array of nodes on the Univ. of Washington HPC. However, the code has been left here to use for those who do not have access to such resources.

1.1 Inputs

Expects trimmed/repaired FastQs generated by 00.00-trimming-fastp.Rmd, with the following format:

• *fastp-trim*.fq.gz

1.2 Outputs

The expected outputs will be:

• *_R1_001.fastp-trim_bismark_bt2_PE_report.txt: A text file summarizing the alignment input and results. Despite the R1 naming, these reports are based on paired reads; the R1 naming is a quirk of Bismark.

• *_R1_001.fastp-trim_bismark_bt2_pe.bam:A BAM alignment. Despite the R1 naming, these BAMs are paired reads; the R1 naming is a quirk of Bismark.

• bismark_summary.txt: An overall summary report of the alignment process. Essentially, this is all of the individual *report.txt files combined into a single file.

• multiqc_report.html: A summary report of the alignment results generated by MultiQC, in HTML format.

Due to the large file sizes of BAMS, these cannot be hosted in the ceasmallr GitHub repo. As such these files are available for download here:

• https://gannet.fish.washington.edu/gitrepos/ceasmallr/

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/ceasmallr'
echo 'export output_dir_top=${repo_dir}/output/02.00-bismark-bowtie2-alignment'
echo 'export trimmed_reads_url="https://gannet.fish.washington.edu/Atumefaciens/20220826-cvir-larvae_zygote-BSseq-fastp_trimming/"'
echo 'export trimmed_fastqs_dir="${repo_dir}/output/00.00-trimming-fastp"'
echo ""

echo "# Input files"
echo 'export bisulfite_genome_url="http://owl.fish.washington.edu/halfshell/genomic-databank/Cvirginica_v300_bisulfite.tar.gz"'
echo 'export bisulfite_genome_gz="Cvirginica_v300_bisulfite.tar.gz"'
echo 'export bisulfite_genome_dir="${repo_dir}/data/Cvirginica_v300"'
echo ""

echo "# Paths to programs"
echo 'export programs_dir="/home/shared"'
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 multiqc="/home/sam/programs/mambaforge/bin/multiqc"'
echo 'export samtools_dir="${programs_dir}/samtools-1.12"'
echo ""

echo "# Program settings"
echo 'export bowtie2_min_score="L,0,-0.6"'
echo ""

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

echo "export trimmed_fastq_pattern='*fastp-trim*.fq.gz'"
echo ""

echo "# Set number of CPUs to use"
echo 'export threads=10'
echo "# Bismark already spawns multiple instances and additional threads are multiplicative."
echo 'export bismark_threads=8' 
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/ceasmallr
export output_dir_top=${repo_dir}/output/02.00-bismark-bowtie2-alignment
export trimmed_reads_url="https://gannet.fish.washington.edu/Atumefaciens/20220826-cvir-larvae_zygote-BSseq-fastp_trimming/"
export trimmed_fastqs_dir="${repo_dir}/output/00.00-trimming-fastp"

# Input files
export bisulfite_genome_url="http://owl.fish.washington.edu/halfshell/genomic-databank/Cvirginica_v300_bisulfite.tar.gz"
export bisulfite_genome_gz="Cvirginica_v300_bisulfite.tar.gz"
export bisulfite_genome_dir="${repo_dir}/data/Cvirginica_v300"

# Paths to programs
export programs_dir="/home/shared"
export bismark_dir="${programs_dir}/Bismark-0.24.0"
export bowtie2_dir="${programs_dir}/bowtie2-2.4.4-linux-x86_64"
export multiqc="/home/sam/programs/mambaforge/bin/multiqc"
export samtools_dir="${programs_dir}/samtools-1.12"

# Program settings
export bowtie2_min_score="L,0,-0.6"

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

# Set number of CPUs to use
export threads=10
# Bismark already spawns multiple instances and additional threads are multiplicative.
export bismark_threads=8

## 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 Download bisulfite genome

# Load bash variables into memory
source .bashvars

cd "${repo_dir}"/data

wget --quiet \
--continue \
"${bisulfite_genome_url}"

ls -ltrh

3.1 Unpack bisulfite genome

# Load bash variables into memory
source .bashvars

cd "${repo_dir}"/data

if [ ! -d "${bisulfite_genome_dir}" ]; then
  tar -xzf "${bisulfite_genome_gz}"
fi

tree "${bisulfite_genome_dir}"

4 Alignment

4.1 Find read pairs

This step is not needed and was only used as a “quick fix” measure to generate a file for use by 02.01-bismark-bowtie2-alignment-SLURM-array.sh.

# Load bash variables into memory
source .bashvars


cd "${trimmed_fastqs_dir}"

if [[ -f "${output_dir_top}"/fastq_pairs.txt ]]; then
  rm "${output_dir_top}"/fastq_pairs.txt
fi

# Find all _R1_ files and match them with their corresponding _R2_ files
for R1_file in *_R1_*.fq.gz; do
    R2_file="${R1_file/_R1_/_R2_}"
    if [[ -f "$R2_file" ]]; then
        echo "$R1_file $R2_file" >> "${output_dir_top}"/fastq_pairs.txt
    else
        echo "Warning: No matching R2 file for $R1_file"
    fi
done


cat "${output_dir_top}"/fastq_pairs.txt

CF01-CM01-Zygote_R1_001.fastp-trim.REPAIRED.fq.gz CF01-CM01-Zygote_R2_001.fastp-trim.REPAIRED.fq.gz
CF01-CM02-Larvae_R1_001.fastp-trim.fq.gz CF01-CM02-Larvae_R2_001.fastp-trim.fq.gz
CF02-CM02-Zygote_R1_001.fastp-trim.fq.gz CF02-CM02-Zygote_R2_001.fastp-trim.fq.gz
CF03-CM03-Zygote_R1_001.fastp-trim.fq.gz CF03-CM03-Zygote_R2_001.fastp-trim.fq.gz
CF03-CM04-Larvae_R1_001.fastp-trim.fq.gz CF03-CM04-Larvae_R2_001.fastp-trim.fq.gz
CF03-CM05-Larvae_R1_001.fastp-trim.fq.gz CF03-CM05-Larvae_R2_001.fastp-trim.fq.gz
CF04-CM04-Zygote_R1_001.fastp-trim.fq.gz CF04-CM04-Zygote_R2_001.fastp-trim.fq.gz
CF05-CM02-Larvae_R1_001.fastp-trim.fq.gz CF05-CM02-Larvae_R2_001.fastp-trim.fq.gz
CF05-CM05-Zygote_R1_001.fastp-trim.fq.gz CF05-CM05-Zygote_R2_001.fastp-trim.fq.gz
CF06-CM01-Zygote_R1_001.fastp-trim.fq.gz CF06-CM01-Zygote_R2_001.fastp-trim.fq.gz
CF06-CM02-Larvae_R1_001.fastp-trim.fq.gz CF06-CM02-Larvae_R2_001.fastp-trim.fq.gz
CF07-CM02-Zygote_R1_001.fastp-trim.fq.gz CF07-CM02-Zygote_R2_001.fastp-trim.fq.gz
CF08-CM03-Zygote_R1_001.fastp-trim.fq.gz CF08-CM03-Zygote_R2_001.fastp-trim.fq.gz
CF08-CM04-Larvae_R1_001.fastp-trim.REPAIRED.fq.gz CF08-CM04-Larvae_R2_001.fastp-trim.REPAIRED.fq.gz
CF08-CM05-Larvae_R1_001.fastp-trim.fq.gz CF08-CM05-Larvae_R2_001.fastp-trim.fq.gz
EF01-EM01-Zygote_R1_001.fastp-trim.fq.gz EF01-EM01-Zygote_R2_001.fastp-trim.fq.gz
EF02-EM02-Zygote_R1_001.fastp-trim.fq.gz EF02-EM02-Zygote_R2_001.fastp-trim.fq.gz
EF03-EM03-Zygote_R1_001.fastp-trim.REPAIRED.fq.gz EF03-EM03-Zygote_R2_001.fastp-trim.REPAIRED.fq.gz
EF03-EM04-Larvae_R1_001.fastp-trim.fq.gz EF03-EM04-Larvae_R2_001.fastp-trim.fq.gz
EF03-EM05-Larvae_R1_001.fastp-trim.fq.gz EF03-EM05-Larvae_R2_001.fastp-trim.fq.gz
EF04-EM04-Zygote_R1_001.fastp-trim.fq.gz EF04-EM04-Zygote_R2_001.fastp-trim.fq.gz
EF04-EM05-Larvae_R1_001.fastp-trim.fq.gz EF04-EM05-Larvae_R2_001.fastp-trim.fq.gz
EF05-EM01-Larvae_R1_001.fastp-trim.fq.gz EF05-EM01-Larvae_R2_001.fastp-trim.fq.gz
EF05-EM05-Zygote_R1_001.fastp-trim.fq.gz EF05-EM05-Zygote_R2_001.fastp-trim.fq.gz
EF05-EM06-Larvae_R1_001.fastp-trim.fq.gz EF05-EM06-Larvae_R2_001.fastp-trim.fq.gz
EF06-EM01-Larvae_R1_001.fastp-trim.fq.gz EF06-EM01-Larvae_R2_001.fastp-trim.fq.gz
EF06-EM02-Larvae_R1_001.fastp-trim.fq.gz EF06-EM02-Larvae_R2_001.fastp-trim.fq.gz
EF06-EM06-Larvae_R1_001.fastp-trim.fq.gz EF06-EM06-Larvae_R2_001.fastp-trim.fq.gz
EF07-EM01-Zygote_R1_001.fastp-trim.fq.gz EF07-EM01-Zygote_R2_001.fastp-trim.fq.gz
EF07-EM03-Larvae_R1_001.fastp-trim.fq.gz EF07-EM03-Larvae_R2_001.fastp-trim.fq.gz
EF08-EM03-Larvae_R1_001.fastp-trim.fq.gz EF08-EM03-Larvae_R2_001.fastp-trim.fq.gz
EF08-EM04-Larvae_R1_001.fastp-trim.fq.gz EF08-EM04-Larvae_R2_001.fastp-trim.fq.gz

4.2 Perform alignments

# Load bash variables into memory
source .bashvars

cd "${trimmed_fastqs_dir}"

find ${trimmed_fastqs_dir}/${R1_fastq_pattern} \
| xargs basename -s "${R1_reads_basename}" \
| xargs -I {} ${bismark_dir}/bismark \
--path_to_bowtie2 ${bowtie2_dir} \
--genome ${bisulfite_genome_dir} \
--score_min "${bowtie2_min_score}" \
--parallel "${bismark_threads}" \
--non_directional \
--samtools_path "${samtools_dir}" \
--gzip \
-p "${threads}" \
-1 ${trimmed_fastqs_dir}/{}"${R1_reads_basename}" \
-2 ${trimmed_fastqs_dir}/{}"${R2_reads_basename}" \
--output_dir "${output_dir_top}" \
2> "${output_dir_top}"/bismark_summary.txt

4.3 List outputs

# Load bash variables into memory
source .bashvars

cd "${output_dir_top}"

ls -lh

5 MultiQC

# Load bash variables into memory
source .bashvars

cd "${output_dir_top}"

${multiqc} .

6 Checksums

# Load bash variables into memory
source .bashvars

cd "${output_dir_top}"

for file in *; do
  if [ "${file}" != "checksums.md5" ]; then
    md5sum "${file}" >> checksums.md5
  fi
done

Citations

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.

Krueger, Felix, and Simon R. Andrews. 2011. “Bismark: A Flexible Aligner and Methylation Caller for Bisulfite-Seq Applications.” Bioinformatics 27 (11): 1571–72. https://doi.org/10.1093/bioinformatics/btr167.

Langmead, Ben, and Steven L Salzberg. 2012. “Fast Gapped-Read Alignment with Bowtie 2.” Nature Methods 9 (4): 357–59. https://doi.org/10.1038/nmeth.1923.

Langmead, Ben, Christopher Wilks, Valentin Antonescu, and Rone Charles. 2018. “Scaling Read Aligners to Hundreds of Threads on General-Purpose Processors.” Edited by John Hancock. Bioinformatics 35 (3): 421–32. https://doi.org/10.1093/bioinformatics/bty648.