Long Read QC and mapping with slurm! - mestato/EPP622_2026 GitHub Wiki

Background

Two dogwood (Cornus florida) cultivars from UTK, 'Cherokee Brave' and 'Appalachian Spring' were sequenced with PacBio Hifi with the goal of assembling high quality, haplotype resolved genome assemblies and mapping the locus for leaf and bract color. See Hamm et al 2025.

Project Set up

We need software. Lets start a new conda environment

conda create -n svs fastplong minimap2 cutesv samtools
conda activate svs

All the files for this lab are in /lustre/isaac24/proj/UTK0505/module3_data .

But if you ever need to download ncbi files, here's how to do it with sra-toolkit: You do not need to do this

conda install sra-tools

# HiFi sequencing of Appalachian Spring Run 2 (SRR34660596)
prefetch SRR34660596
fasterq-dump SRR34660596

# Cherokee Brave HiFi Run 3 (SRR34635021)
prefetch SRR34635021
fasterq-dump SRR34635021

# Compress
pigz *.fastq 

# Clean up the large temporary directories left behind
rm -rf SRR34635021 SRR34660596

# Rename read files to have meaningful names
mv SRR34660596.fastq.gz AppSpring2.fq.gz
mv SRR34635021.fastq.gz CherBrave3.fq.gz

Directory Set up

Lets start a new space for this lab session.

cd /lustre/isaac24/proj/UTK0505/<yourusername>
mkdir module3
cd module3
mkdir 01_fastplong

Fastplong - QC and trimming in one!

Lets use fastplong to quality trim, adapter trim, and output qc stats in one go. The default number of threads is 3, we'll increase that to 10. There are many options, but we will stick with the defaults for the rest.

We'll create a slurm script to run fastplong for both files sequentially. Name it fastplong.qsh. Note that since i am running it from the directory /lustre/isaac24/proj/UTK0505/mstaton1/01_fastplong, the output files will land in that directory.

Don't forget to change yourusername to your user name.

#!/bin/bash
#SBATCH --job-name fastplong
#SBATCH --nodes=1
#SBATCH --cpus-per-task=10
#SBATCH -A ISAAC-UTK0505
#SBATCH --partition=short
#SBATCH --qos=short
#SBATCH --time=02:00:00
#SBATCH --mail-type=BEGIN,END,FAIL
#SBATCH [email protected]

# Activate the conda environment
eval "$(conda shell.bash hook)"
conda activate svs

# Run fastplong on each of two sequence sets
fastplong \
  --in /lustre/isaac24/proj/UTK0505/module3_data/AppSpring2.fq.gz \
  --out AppSpring2.trimmed.fq.gz \
  --html AppSpring2.trimmed.html \
  --thread 10
fastplong \
  --in /lustre/isaac24/proj/UTK0505/module3_data/CherBrave3.fq.gz \
  --out CherBrave3.trimmed.fq.gz \
  --html CherBrave3.trimmed.html \
  --thread 10

Run it

sbatch fastplong.qah

Check on it

squeue -u <yourusername>

Explore your output/error files.

Minimap2

We will map reads with minimap2. I already made the reference index in /lustre/isaac24/proj/UTK0505/module3_data so you don't have to do this. But here are the commands to index with minimap2 and samtools.

samtools faidx AppSpring_Hap1_Chr9.fasta
minimap2 -d AppSpring_Hap1_Chr9.mmi  AppSpring_Hap1_Chr9.fasta

We'll use a for loop even though there are only two samples - just to show how you would scale up. We'll create a slurm script minimap2.qsh. Run from new directory, /lustre/isaac24/proj/UTK0505/mstaton1/02_minimap2, the output files will land in that directory.

Don't forget to change yourusername to your user name.

#!/bin/bash
#SBATCH --job-name minimap2
#SBATCH --nodes=1
#SBATCH --cpus-per-task=8
#SBATCH -A ISAAC-UTK0505
#SBATCH --partition=campus
#SBATCH --qos=campus
#SBATCH --time=08:00:00
#SBATCH --mem=60G
#SBATCH --mail-type=BEGIN,END,FAIL
#SBATCH [email protected]

# Activate the conda environment
eval "$(conda shell.bash hook)"
conda activate svs

# Run minimap2 on two sequence sets
for FASTQ in /lustre/isaac24/proj/UTK0505/youruserid/module3/01_fastplong/*trimmed.fq.gz
do

  # The */ matches any text ending with a forward slash.
  SAMPLE="${FASTQ##*/}"
  # second, remove trimmed.fq.gz from the file name
  SAMPLE="${SAMPLE%.trimmed.fq.gz}"

  # minimap2 memory is based on index size; samtools is memory per thread

  minimap2 \
  -a -x map-hifi \
  -t 7 \
  /lustre/isaac24/proj/UTK0505/module3_data/AppSpring_Hap1_Chr9.fasta \
  $FASTQ \
  | samtools sort -@ 7 -m 6G -O BAM \
  -o $SAMPLE.bam

  samtools index $SAMPLE.bam
  samtools stats $SAMPLE.bam > $SAMPLE.stats.txt
  samtools coverage $SAMPLE.bam > $SAMPLE.coverage.txt
done