BCH709 Introduction to Bioinformatics: Resequencing on Pronghorn HPC (Parallelized with Slurm)

ssh <YOUR_NET_ID>@pronghorn.rc.unr.edu

# Confirm your Slurm account/partition (required for sbatch)
sacctmgr show user $USER format=user,account,partition,qos%-30 \
  withassoc where user=$USER

micromamba create -n reseq_bch709 -c conda-forge -c bioconda python=3.11 -y
micromamba activate reseq_bch709

micromamba install -c conda-forge -c bioconda \
    fastqc 'fastp>=0.24' bwa-mem2 \
    'samtools>=1.20' 'bcftools>=1.20' 'tabix>=1.11' \
    openjdk=21 'picard>=3' gatk4 snpeff plink -y

# Note: we deliberately do NOT install sra-tools. Bioconda's sra-tools 3.x is
# built against GLIBC 2.27+, which is newer than Pronghorn's system libc — the
# binary fails with "GLIBC_2.27 not found" on the compute nodes. The download
# script in Step 1 pulls FASTQ directly from ENA over HTTPS, so sra-tools is
# not needed.

# Upgrade pip first — older pip can't find the prebuilt `tiktoken`
# manylinux wheel (a transitive multiqc dep), tries to build it from
# Rust source, fails on Pronghorn (no Rust compiler).
pip install --upgrade pip

# MultiQC + pinned deps. `tiktoken<0.8` is the safety pin — older
# tiktoken has stable cp311 linux wheels.
# IMPORTANT: don't skip this — step 8 needs multiqc.
pip install --prefer-binary \
    'numpy<2.0' 'pyarrow<17' 'tiktoken<0.8' 'multiqc<1.34'

# Must be run with the env ACTIVE — $CONDA_PREFIX points at the env folder
cd "$CONDA_PREFIX/lib"
if   [ -f libcrypto.so.1.1 ]; then ln -sf libcrypto.so.1.1 libcrypto.so.1.0.0
elif [ -f libcrypto.so.3   ]; then ln -sf libcrypto.so.3   libcrypto.so.1.0.0
fi
cd - > /dev/null

# Verify the three most-used HTS tools all run
samtools --version  | head -1    # → samtools 1.xx
bcftools --version  | head -1    # → bcftools 1.xx
tabix --version     | head -1    # → tabix (htslib) 1.xx

mkdir -p ~/scratch/reseq/{raw,trim,bam,vcf,logs,scripts}
cd ~/scratch/reseq

printf 'sample1\tSRR519585\nsample2\tSRR519586\n' > ~/scratch/reseq/samples.tsv

awk -F'\t' '{print NF}' ~/scratch/reseq/samples.tsv | sort -u
# Expected output:  2
# If you see 1, you have spaces instead of tabs — re-run the printf above.

cat ~/scratch/reseq/samples.tsv          # quick visual check

#!/bin/bash
#SBATCH --job-name=download
#SBATCH --account=cpu-s5-bch709-6
#SBATCH --partition=cpu-core-0
#SBATCH --cpus-per-task=2
#SBATCH --mem=4g
#SBATCH --time=12:00:00
#SBATCH --array=1-2                          # one task per sample (rows in samples.tsv)
#SBATCH --mail-type=FAIL,END
#SBATCH --mail-user=<YOUR_EMAIL>
#SBATCH -o logs/01_download_%A_%a.out        # %A=array job id, %a=task index

set -euo pipefail

cd ~/scratch/reseq

# Pick the N-th line from samples.tsv using the array task ID
LINE=$(sed -n "${SLURM_ARRAY_TASK_ID}p" samples.tsv)
SAMPLE=$(echo "$LINE" | cut -f1)
SRR=$(echo "$LINE" | cut -f2)

echo "[task ${SLURM_ARRAY_TASK_ID}] Downloading ${SAMPLE} (${SRR}) from ENA"

# Download paired-end FASTQ from ENA (mirrors NCBI SRA, gzipped FASTQ ready-to-use).
# We deliberately avoid `prefetch` / `fasterq-dump` — bioconda's sra-tools 3.x is
# built against GLIBC 2.27+, which is newer than Pronghorn's system libc, so the
# binary fails with "GLIBC_2.27 not found" on the compute nodes.
mkdir -p raw

# Ask ENA for the exact fastq URLs (handles paired-end / single-end / multi-file).
META_URL="https://www.ebi.ac.uk/ena/portal/api/filereport?accession=${SRR}&result=read_run&fields=fastq_ftp&format=tsv"
URLS=$(curl -fsSL --retry 3 --max-time 60 "${META_URL}" \
        | tail -n +2 \
        | awk -F'\t' '{print $NF}' \
        | tr ';' '\n' \
        | sed '/^$/d')
[ -n "${URLS}" ] || { echo "ERROR: ENA returned no fastq URLs for ${SRR}"; exit 1; }

i=1
# --retry-all-errors + -C -: ENA's CDN occasionally drops mid-transfer with an
# SSL eof; --retry-all-errors retries on transport errors and -C - resumes from
# the partial file instead of restarting at byte 0 (see issue #64).
for U in ${URLS}; do
    OUT="raw/${SAMPLE}_R${i}.fastq.gz"
    echo "[task ${SLURM_ARRAY_TASK_ID}] -> https://${U}"
    curl -fsSL --retry 5 --retry-all-errors --retry-delay 30 --max-time 3600 -C - -o "${OUT}" "https://${U}"
    # Don't use a `[ -s "${OUT}" ]` partial-file skip-check here: with -C - a
    # short file from a prior failed attempt would falsely pass. curl -f exits
    # non-zero on real failures, which set -e turns into a hard stop.
    i=$((i+1))
done

ls -lh raw/${SAMPLE}_R*.fastq.gz

# Make sure the env is active in THIS shell (once per login session)
micromamba activate reseq_bch709
cd ~/scratch/reseq

DL_JID=$(sbatch --parsable scripts/01_download.sh)
echo "Download job: ${DL_JID}"
squeue -u $USER   # you should see 2 tasks (download_1, download_2) running in parallel

#!/bin/bash
#SBATCH --job-name=reference
#SBATCH --account=cpu-s5-bch709-6
#SBATCH --partition=cpu-core-0
#SBATCH --cpus-per-task=4
#SBATCH --mem=16g
#SBATCH --time=06:00:00
#SBATCH --mail-type=FAIL,END
#SBATCH --mail-user=<YOUR_EMAIL>
#SBATCH -o logs/02_reference_%j.out

set -euo pipefail
# Environment is activated in the login shell before `sbatch` is called —
# the PATH (with bwa-mem2, samtools, gatk, …) is inherited automatically.

cd ~/scratch/reseq

# ---- Download TAIR10 from TAIR (www.arabidopsis.org) ----
# Primary: TAIR (canonical source) — uses self-signed cert, so curl needs -k
# Fallback: EBI Ensembl Plants mirror (most stable mirror; NCBI is last resort)
TAIR_URL="https://www.arabidopsis.org/api/download-files/download?filePath=Genes/TAIR10_genome_release/TAIR10_chromosome_files/TAIR10_chr_all.fas.gz"
EBI_URL="https://ftp.ebi.ac.uk/ensemblgenomes/pub/plants/release-60/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz"
NCBI_URL="https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/735/GCF_000001735.4_TAIR10.1/GCF_000001735.4_TAIR10.1_genomic.fna.gz"

OUT=reference.fasta.gz
rm -f "$OUT"
for URL in "$TAIR_URL" "$EBI_URL" "$NCBI_URL"; do
    echo "[ref] trying $URL"
    if curl -kfsSL --retry 3 --max-time 900 -o "$OUT" "$URL" && [ -s "$OUT" ]; then
        echo "[ref] downloaded from $URL"
        break
    fi
    rm -f "$OUT"
done
[ -s "$OUT" ] || { echo "ERROR: TAIR10 download failed from all sources"; exit 1; }
gunzip -f "$OUT"

# Normalize chromosome names to match the 1001genomes VCF (1..5, Mt, Pt):
#   - TAIR fasta uses Chr1..Chr5, ChrM, ChrC
#   - NCBI RefSeq uses NC_003070.9 …
#   - EBI Ensembl already uses 1..5, Mt, Pt (no rename needed)
if grep -q '^>Chr' reference.fasta; then
    sed -i -E 's/^>Chr([0-9]+).*/>\1/; s/^>ChrM.*/>Mt/; s/^>ChrC.*/>Pt/' reference.fasta
elif grep -q '^>NC_' reference.fasta; then
    awk 'BEGIN{
        m["NC_003070.9"]="1"; m["NC_003071.7"]="2"; m["NC_003074.8"]="3"
        m["NC_003075.7"]="4"; m["NC_003076.8"]="5"
        m["NC_037304.1"]="Mt"; m["NC_000932.1"]="Pt"
    }
    /^>/{
        name=substr($1,2); print ">" (name in m ? m[name] : name); next
    }
    {print}' reference.fasta > reference.renamed.fa && mv -f reference.renamed.fa reference.fasta
fi

# Build all indices.  bwa-mem2 index / samtools faidx overwrite their outputs by
# default, but Picard CreateSequenceDictionary refuses to overwrite an existing
# .dict and aborts with `file:///.../reference.dict already exists.  Delete this
# file and try again`.  rm -f makes the whole step idempotent so a re-run of
# 02_reference.sh after a partial failure just works.
bwa-mem2 index reference.fasta
samtools faidx reference.fasta
rm -f reference.dict
picard CreateSequenceDictionary R=reference.fasta O=reference.dict

# Download and prepare known variants (sites-only to avoid malformed VCF)
# 1001genomes.org is sometimes slow — give curl up to 30 min and retry 3×
KSITES_URL="https://1001genomes.org/data/GMI-MPI/releases/v3.1/1001genomes_snp-short-indel_only_ACGTN.vcf.gz"
# The 1001genomes VCF is ~19 GB; on Pronghorn outbound (~3 MB/s) the download
# takes ~3 hours. We use `-C -` to resume on retry instead of restarting at
# byte 0, and drop --max-time so a single attempt isn't capped (Slurm --time
# already bounds the whole job).
# --retry-all-errors + -C -: 1001genomes occasionally drops the connection
# with an SSL eof on a 19 GB download; --retry-all-errors retries on transport
# errors and -C - resumes from the partial file (see issue #64).
curl -fsSL --retry 5 --retry-all-errors --retry-delay 30 -C - \
    -o 1001genomes_snp-short-indel_only_ACGTN.vcf.gz "${KSITES_URL}"
# Atomic rename: write to a .tmp file, only promote to known_sites.vcf.gz once
# bgzip + tabix BOTH succeed and outputs are non-empty. If anything in this
# block fails, the original 1001genomes file is preserved so the student can
# re-run step 02 without re-downloading 19 GB.
zcat 1001genomes_snp-short-indel_only_ACGTN.vcf.gz \
    | awk 'BEGIN{OFS="\t"} /^##/{print; next} /^#CHROM/{print $1,$2,$3,$4,$5,$6,$7,$8; next} {print $1,$2,$3,$4,$5,$6,$7,$8}' \
    | bgzip > known_sites.vcf.gz.tmp
tabix -p vcf known_sites.vcf.gz.tmp
[ -s known_sites.vcf.gz.tmp ] && [ -s known_sites.vcf.gz.tmp.tbi ] || {
    echo "ERROR: known_sites.vcf.gz build failed (empty output). Re-run step 02." >&2
    rm -f known_sites.vcf.gz.tmp known_sites.vcf.gz.tmp.tbi
    exit 1
}
mv -f known_sites.vcf.gz.tmp known_sites.vcf.gz
mv -f known_sites.vcf.gz.tmp.tbi known_sites.vcf.gz.tbi
rm -f 1001genomes_snp-short-indel_only_ACGTN.vcf.gz

echo "Reference prep done."

REF_JID=$(sbatch --parsable scripts/02_reference.sh)
echo "Reference job: ${REF_JID}"

ls -lh ~/scratch/reseq/reference.*

ls -lh ~/scratch/reseq/known_sites.vcf.gz*

#!/bin/bash
#SBATCH --job-name=align
#SBATCH --account=cpu-s5-bch709-6
#SBATCH --partition=cpu-core-0
#SBATCH --cpus-per-task=8
#SBATCH --mem=32g
#SBATCH --time=12:00:00
#SBATCH --array=1-2
#SBATCH --mail-type=FAIL,END
#SBATCH --mail-user=<YOUR_EMAIL>
#SBATCH -o logs/03_align_%A_%a.out

set -euo pipefail
# Environment is activated in the login shell before `sbatch` is called —
# the PATH (with bwa-mem2, samtools, gatk, …) is inherited automatically.

cd ~/scratch/reseq

LINE=$(sed -n "${SLURM_ARRAY_TASK_ID}p" samples.tsv)
SAMPLE=$(echo "$LINE" | cut -f1)

# ---- Trim ----
fastp \
    --in1  raw/${SAMPLE}_R1.fastq.gz \
    --in2  raw/${SAMPLE}_R2.fastq.gz \
    --out1 trim/${SAMPLE}_R1.trimmed.fq.gz \
    --out2 trim/${SAMPLE}_R2.trimmed.fq.gz \
    --detect_adapter_for_pe \
    --qualified_quality_phred 20 \
    --length_required 50 \
    --thread ${SLURM_CPUS_PER_TASK} \
    --html trim/${SAMPLE}_fastp.html \
    --json trim/${SAMPLE}_fastp.json

# ---- Align + sort ----
set -o pipefail
bwa-mem2 mem \
    -t ${SLURM_CPUS_PER_TASK} \
    -R "@RG\tID:${SAMPLE}\tSM:${SAMPLE}\tPL:ILLUMINA\tLB:lib_${SAMPLE}\tPU:unit_${SAMPLE}" \
    reference.fasta \
    trim/${SAMPLE}_R1.trimmed.fq.gz \
    trim/${SAMPLE}_R2.trimmed.fq.gz \
  | samtools sort -@ ${SLURM_CPUS_PER_TASK} -T /tmp/${SAMPLE}_sort \
      -o bam/${SAMPLE}.bam

samtools index -@ ${SLURM_CPUS_PER_TASK} bam/${SAMPLE}.bam
samtools quickcheck bam/${SAMPLE}.bam && echo "BAM OK"

ls -lh ~/scratch/reseq/trim/sample1_fastp.* ~/scratch/reseq/bam/sample1.bam

samtools flagstat ~/scratch/reseq/bam/sample1.bam

ALIGN_JID=$(sbatch --parsable --dependency=afterok:${DL_JID}:${REF_JID} scripts/03_align.sh)
echo "Align: ${ALIGN_JID}"
squeue -u $USER   # ALIGN_JID should show state PD with reason (Dependency)

#!/bin/bash
#SBATCH --job-name=markdup_bqsr
#SBATCH --account=cpu-s5-bch709-6
#SBATCH --partition=cpu-core-0
#SBATCH --cpus-per-task=4
#SBATCH --mem=16g
#SBATCH --time=08:00:00
#SBATCH --array=1-2
#SBATCH --mail-type=FAIL,END
#SBATCH --mail-user=<YOUR_EMAIL>
#SBATCH -o logs/04_markdup_bqsr_%A_%a.out

set -euo pipefail
# Environment is activated in the login shell before `sbatch` is called —
# the PATH (with bwa-mem2, samtools, gatk, …) is inherited automatically.

cd ~/scratch/reseq

SAMPLE=$(sed -n "${SLURM_ARRAY_TASK_ID}p" samples.tsv | cut -f1)

# ---- MarkDuplicates ----
picard -Xmx12g MarkDuplicates \
    I=bam/${SAMPLE}.bam \
    O=bam/${SAMPLE}.markdup.bam \
    M=bam/${SAMPLE}.markdup.metrics \
    VALIDATION_STRINGENCY=SILENT
samtools index -@ ${SLURM_CPUS_PER_TASK} bam/${SAMPLE}.markdup.bam

# ---- Preflight: known_sites.vcf.gz must exist (built by step 02) ----
if [ ! -s known_sites.vcf.gz ] || [ ! -s known_sites.vcf.gz.tbi ]; then
    echo "ERROR: known_sites.vcf.gz (and .tbi) not found in $(pwd)." >&2
    echo "       Step 02 (02_reference.sh) did not finish — re-run it before step 04." >&2
    ls -lh 1001genomes_*.vcf.gz known_sites.vcf.gz* 2>/dev/null || true
    exit 1
fi

# ---- BaseRecalibrator ----
gatk --java-options "-Xmx12g" BaseRecalibrator \
    -I bam/${SAMPLE}.markdup.bam \
    -R reference.fasta \
    --known-sites known_sites.vcf.gz \
    -O bam/${SAMPLE}.recal.table

# ---- ApplyBQSR ----
gatk --java-options "-Xmx12g" ApplyBQSR \
    -I bam/${SAMPLE}.markdup.bam \
    -R reference.fasta \
    --bqsr-recal-file bam/${SAMPLE}.recal.table \
    -O bam/${SAMPLE}.recal.bam
samtools index -@ ${SLURM_CPUS_PER_TASK} bam/${SAMPLE}.recal.bam

# Clean up intermediate markdup BAM
rm -f bam/${SAMPLE}.bam bam/${SAMPLE}.bam.bai
rm -f bam/${SAMPLE}.markdup.bam bam/${SAMPLE}.markdup.bam.bai

head -10 ~/scratch/reseq/bam/sample1.markdup.metrics

ls -lh ~/scratch/reseq/bam/sample1.recal.bam*

MARKDUP_JID=$(sbatch --parsable --dependency=afterok:${ALIGN_JID} scripts/04_markdup_bqsr.sh)
echo "MarkDup+BQSR: ${MARKDUP_JID}"

#!/bin/bash
#SBATCH --job-name=hc_scatter
#SBATCH --account=cpu-s5-bch709-6
#SBATCH --partition=cpu-core-0
#SBATCH --cpus-per-task=4
#SBATCH --mem=16g
#SBATCH --time=12:00:00
#SBATCH --array=1-14                  # N_SAMPLES * N_CHROMS = 2 * 7
#SBATCH --mail-type=FAIL,END
#SBATCH --mail-user=<YOUR_EMAIL>
#SBATCH -o logs/05a_hc_%A_%a.out

set -euo pipefail
# Environment is activated in the login shell before `sbatch` is called —
# the PATH (with bwa-mem2, samtools, gatk, …) is inherited automatically.

cd ~/scratch/reseq

# Map the linear task ID to (sample_index, chr_index)
CHROMS=(1 2 3 4 5 Mt Pt)
N_CHROMS=${#CHROMS[@]}

# 0-indexed within the array
IDX=$((SLURM_ARRAY_TASK_ID - 1))
SAMPLE_IDX=$((IDX / N_CHROMS + 1))     # 1-indexed line in samples.tsv
CHR_IDX=$((IDX % N_CHROMS))

SAMPLE=$(sed -n "${SAMPLE_IDX}p" samples.tsv | cut -f1)
CHR=${CHROMS[$CHR_IDX]}

echo "Task ${SLURM_ARRAY_TASK_ID}: sample=${SAMPLE}, chr=${CHR}"

mkdir -p vcf/scatter

gatk --java-options "-Xmx12g" HaplotypeCaller \
    -R reference.fasta \
    -I bam/${SAMPLE}.recal.bam \
    -O vcf/scatter/${SAMPLE}.${CHR}.g.vcf.gz \
    -L ${CHR} \
    -ERC GVCF \
    --native-pair-hmm-threads ${SLURM_CPUS_PER_TASK}

ls -lh ~/scratch/reseq/vcf/scatter/

#!/bin/bash
#SBATCH --job-name=hc_gather
#SBATCH --account=cpu-s5-bch709-6
#SBATCH --partition=cpu-core-0
#SBATCH --cpus-per-task=2
#SBATCH --mem=8g
#SBATCH --time=02:00:00
#SBATCH --array=1-2                   # one task per sample
#SBATCH --mail-type=FAIL,END
#SBATCH --mail-user=<YOUR_EMAIL>
#SBATCH -o logs/05b_gather_%A_%a.out

set -euo pipefail
# Environment is activated in the login shell before `sbatch` is called —
# the PATH (with bwa-mem2, samtools, gatk, …) is inherited automatically.

cd ~/scratch/reseq

SAMPLE=$(sed -n "${SLURM_ARRAY_TASK_ID}p" samples.tsv | cut -f1)

# Merge the 7 chromosome GVCFs for this sample into one
CHROMS=(1 2 3 4 5 Mt Pt)
INPUTS=""
for CHR in "${CHROMS[@]}"; do
    INPUTS+=" -I vcf/scatter/${SAMPLE}.${CHR}.g.vcf.gz"
done

# MergeVcfs sorts inputs against the reference sequence dictionary, so
# the order of -I arguments does not matter. (GatherVcfs is faster but
# requires inputs to be in strict dict order — TAIR10 has Pt before Mt,
# which is easy to get wrong.)
gatk --java-options "-Xmx6g" MergeVcfs \
    ${INPUTS} \
    -O vcf/${SAMPLE}.g.vcf.gz

# Verify and clean up scatter files
ls -lh vcf/${SAMPLE}.g.vcf.gz
rm -rf vcf/scatter/${SAMPLE}.*.g.vcf.gz vcf/scatter/${SAMPLE}.*.g.vcf.gz.tbi

HC_JID=$(sbatch --parsable --dependency=afterok:${MARKDUP_JID} scripts/05a_haplotypecaller.sh)
GATHER_JID=$(sbatch --parsable --dependency=afterok:${HC_JID} scripts/05b_gather_gvcf.sh)
echo "HaplotypeCaller: ${HC_JID} | Gather: ${GATHER_JID}"

#!/bin/bash
#SBATCH --job-name=joint_geno
#SBATCH --account=cpu-s5-bch709-6
#SBATCH --partition=cpu-core-0
#SBATCH --cpus-per-task=4
#SBATCH --mem=32g
#SBATCH --time=04:00:00
#SBATCH --mail-type=FAIL,END
#SBATCH --mail-user=<YOUR_EMAIL>
#SBATCH -o logs/06_joint_%j.out

set -euo pipefail
# Environment is activated in the login shell before `sbatch` is called —
# the PATH (with bwa-mem2, samtools, gatk, …) is inherited automatically.

cd ~/scratch/reseq

# Build the -V argument list from samples.tsv (auto-scales with sample count)
INPUTS=""
while read SAMPLE SRR; do
    INPUTS+=" -V vcf/${SAMPLE}.g.vcf.gz"
done < samples.tsv

gatk --java-options "-Xmx24g" CombineGVCFs \
    -R reference.fasta \
    ${INPUTS} \
    -O vcf/cohort.g.vcf.gz

gatk --java-options "-Xmx24g" GenotypeGVCFs \
    -R reference.fasta \
    -V vcf/cohort.g.vcf.gz \
    -O vcf/cohort.vcf.gz

echo "Joint genotyping done. Variants:"
bcftools stats vcf/cohort.vcf.gz | grep "^SN"

bcftools view -H ~/scratch/reseq/vcf/cohort.vcf.gz | wc -l

JOINT_JID=$(sbatch --parsable --dependency=afterok:${GATHER_JID} scripts/06_joint_genotype.sh)

#!/bin/bash
#SBATCH --job-name=filter_ann
#SBATCH --account=cpu-s5-bch709-6
#SBATCH --partition=cpu-core-0
#SBATCH --cpus-per-task=4
#SBATCH --mem=16g
#SBATCH --time=04:00:00
#SBATCH --mail-type=FAIL,END
#SBATCH --mail-user=<YOUR_EMAIL>
#SBATCH -o logs/07_filter_%j.out

set -euo pipefail
# Environment is activated in the login shell before `sbatch` is called —
# the PATH (with bwa-mem2, samtools, gatk, …) is inherited automatically.

cd ~/scratch/reseq

# Pre-download the snpEff Arabidopsis database if it's not already present.
# snpEff would auto-download on first use, but doing it explicitly here makes
# the failure mode obvious if the compute node has no outbound HTTPS, and
# prevents the download from adding minutes to every fresh student's first
# step-7 run silently.
if [ ! -d "${HOME}/snpeff_data/Arabidopsis_thaliana" ]; then
    snpEff download -dataDir "${HOME}/snpeff_data" Arabidopsis_thaliana
fi

# ---- Split SNPs / Indels ----
gatk SelectVariants -R reference.fasta -V vcf/cohort.vcf.gz \
    --select-type-to-include SNP   -O vcf/cohort.snps.vcf.gz
gatk SelectVariants -R reference.fasta -V vcf/cohort.vcf.gz \
    --select-type-to-include INDEL -O vcf/cohort.indels.vcf.gz

# ---- Hard filter ----
gatk VariantFiltration -R reference.fasta -V vcf/cohort.snps.vcf.gz \
    --filter-expression "QD < 2.0"               --filter-name "QD2" \
    --filter-expression "FS > 60.0"              --filter-name "FS60" \
    --filter-expression "MQ < 40.0"              --filter-name "MQ40" \
    --filter-expression "MQRankSum < -12.5"      --filter-name "MQRankSum-12.5" \
    --filter-expression "ReadPosRankSum < -8.0"  --filter-name "ReadPosRankSum-8" \
    -O vcf/cohort.snps.filtered.vcf.gz

gatk VariantFiltration -R reference.fasta -V vcf/cohort.indels.vcf.gz \
    --filter-expression "QD < 2.0"                --filter-name "QD2" \
    --filter-expression "FS > 200.0"              --filter-name "FS200" \
    --filter-expression "ReadPosRankSum < -20.0"  --filter-name "ReadPosRankSum-20" \
    -O vcf/cohort.indels.filtered.vcf.gz

# ---- PASS-only VCF ----
bcftools view -f PASS -Oz -o vcf/cohort.snps.pass.vcf.gz vcf/cohort.snps.filtered.vcf.gz
tabix -p vcf vcf/cohort.snps.pass.vcf.gz

# ---- SnpEff annotation ----
# -csvStats writes vcf/snpEff_summary.csv — REQUIRED for MultiQC's snpeff
# module (it parses the CSV, not the HTML summary).
snpEff -csvStats vcf/snpEff_summary.csv \
    -dataDir ${HOME}/snpeff_data Arabidopsis_thaliana \
    vcf/cohort.snps.pass.vcf.gz > vcf/cohort.snps.annotated.vcf 2> logs/snpeff.log

# ---- PLINK (requires 2+ samples) ----
cd vcf
plink --vcf cohort.snps.pass.vcf.gz --make-bed --out cohort       --allow-extra-chr
plink --bfile cohort --pca 10                    --out cohort.pca --allow-extra-chr
plink --bfile cohort --indep-pairwise 50 10 0.2  --out cohort.ld  --allow-extra-chr

echo "Pipeline complete."

head -25 ~/scratch/reseq/logs/snpeff.log

cat ~/scratch/reseq/vcf/cohort.pca.eigenvec

bcftools view -H ~/scratch/reseq/vcf/cohort.snps.pass.vcf.gz | wc -l
bcftools stats ~/scratch/reseq/vcf/cohort.snps.pass.vcf.gz | grep "^SN"

FILTER_JID=$(sbatch --parsable --dependency=afterok:${JOINT_JID} scripts/07_filter_annotate.sh)
echo "Filter+Annotate: ${FILTER_JID}"

#!/bin/bash
#SBATCH --job-name=multiqc
#SBATCH --account=cpu-s5-bch709-6
#SBATCH --partition=cpu-core-0
#SBATCH --cpus-per-task=2
#SBATCH --mem=8g
#SBATCH --time=01:00:00
#SBATCH --mail-type=FAIL,END
#SBATCH --mail-user=<YOUR_EMAIL>
#SBATCH -o logs/08_multiqc_%j.out

set -euo pipefail
cd ~/scratch/reseq

mkdir -p qc

# bcftools stats over the final cohort VCF — MultiQC parses the .vchk
mkdir -p qc/bcftools

# Step 7 emits SNPs and INDELs as two separate filtered files. Concat them
# into one cohort.filtered.vcf.gz so bcftools stats (and MultiQC's bcftools
# module) sees a single record set covering both variant classes. The
# guard makes the step idempotent on re-runs.
if [ ! -s vcf/cohort.filtered.vcf.gz ]; then
    bcftools concat -a -Oz \
        -o vcf/cohort.filtered.vcf.gz \
        vcf/cohort.snps.filtered.vcf.gz \
        vcf/cohort.indels.filtered.vcf.gz
    bcftools index -t vcf/cohort.filtered.vcf.gz
fi
bcftools stats vcf/cohort.filtered.vcf.gz > qc/bcftools/cohort.filtered.vchk

# Aggregate everything multiqc can find under the working dir.
# NOTE: --module gatk and --module rsem are intentionally OMITTED, and
# also explicitly excluded with --exclude as belt-and-suspenders.
# MultiQC 1.34's GATK base_recalibrator parser hits a pydantic
# ValidationError on the BQSR scatter (points.X.name is None) and
# crashes the whole report — the BQSR table is human-readable, inspect
# bam/*.recal.table directly (see Section 8 below).
# MultiQC 1.34's RSEM parser also misreads stray '#'-prefixed files in
# sibling dirs (e.g. ~/scratch/rnaseq/) and ValueErrors out. We don't
# run RSEM in this pipeline at all, so excluding it is free insurance.
# Re-enable both modules once the upstream MultiQC bugs are fixed.
multiqc . -o qc/ -n reseq_report --force \
    --module fastp \
    --module picard \
    --module snpeff \
    --module bcftools \
    --exclude rsem \
    --exclude gatk

echo "MultiQC report: qc/reseq_report.html"

MQC_JID=$(sbatch --parsable --dependency=afterok:${FILTER_JID} scripts/08_multiqc.sh)
echo "MultiQC: ${MQC_JID}"

scp <netid>@pronghorn.rc.unr.edu:~/scratch/reseq/qc/reseq_report.html ./
open reseq_report.html      # or: double-click in your file browser

#!/bin/bash
set -euo pipefail
cd ~/scratch/reseq

# Activate the env in THIS shell so every sbatch below inherits the PATH
# (sbatch --export=ALL is the default — the submitted jobs see the same tools)
export MAMBA_ROOT_PREFIX="${MAMBA_ROOT_PREFIX:-$HOME/micromamba}"
eval "$(micromamba shell hook --shell=bash)"
micromamba activate reseq_bch709

# --parsable returns just the job ID, perfect for chaining
DL_JID=$(sbatch --parsable                                  scripts/01_download.sh)
REF_JID=$(sbatch --parsable                                 scripts/02_reference.sh)
ALIGN_JID=$(sbatch --parsable --dependency=afterok:${DL_JID}:${REF_JID} scripts/03_align.sh)
MARK_JID=$(sbatch --parsable --dependency=afterok:${ALIGN_JID}          scripts/04_markdup_bqsr.sh)
HC_JID=$(sbatch --parsable --dependency=afterok:${MARK_JID}             scripts/05a_haplotypecaller.sh)
GATHER_JID=$(sbatch --parsable --dependency=afterok:${HC_JID}           scripts/05b_gather_gvcf.sh)
JOINT_JID=$(sbatch --parsable --dependency=afterok:${GATHER_JID}        scripts/06_joint_genotype.sh)
FILT_JID=$(sbatch --parsable --dependency=afterok:${JOINT_JID}          scripts/07_filter_annotate.sh)
MQC_JID=$(sbatch --parsable --dependency=afterany:${FILT_JID}           scripts/08_multiqc.sh)

echo "Submitted pipeline:"
printf '  %-20s %s\n' "download"       "${DL_JID}"
printf '  %-20s %s\n' "reference"      "${REF_JID}"
printf '  %-20s %s\n' "align"          "${ALIGN_JID}"
printf '  %-20s %s\n' "markdup+bqsr"   "${MARK_JID}"
printf '  %-20s %s\n' "haplotypecaller" "${HC_JID}"
printf '  %-20s %s\n' "gather"         "${GATHER_JID}"
printf '  %-20s %s\n' "joint_genotype" "${JOINT_JID}"
printf '  %-20s %s\n' "filter+annotate" "${FILT_JID}"
printf '  %-20s %s\n' "multiqc"        "${MQC_JID}"

echo ""
echo "Monitor with:  squeue -u \$USER"
echo "Cancel all:    scancel ${DL_JID} ${REF_JID} ${ALIGN_JID} ${MARK_JID} ${HC_JID} ${GATHER_JID} ${JOINT_JID} ${FILT_JID} ${MQC_JID}"
echo ""
echo "Final report (after pipeline finishes): ~/scratch/reseq/qc/reseq_report.html"

chmod +x scripts/*.sh
bash scripts/run_all.sh

micromamba activate reseq_bch709
cd ~/scratch/reseq

# --- Step 1: download FASTQs (no prerequisites) ---
DL_JID=$(sbatch --parsable scripts/01_download.sh)
echo "download     → $DL_JID"

# --- Step 2: reference prep (no prerequisites, runs in parallel with Step 1) ---
REF_JID=$(sbatch --parsable scripts/02_reference.sh)
echo "reference    → $REF_JID"

# --- Step 3: align (waits for BOTH download and reference) ---
ALIGN_JID=$(sbatch --parsable --dependency=afterok:${DL_JID}:${REF_JID} scripts/03_align.sh)
echo "align        → $ALIGN_JID"

# --- Step 4: markdup + BQSR (waits for align) ---
MARK_JID=$(sbatch --parsable --dependency=afterok:${ALIGN_JID} scripts/04_markdup_bqsr.sh)
echo "markdup+bqsr → $MARK_JID"

# --- Step 5a: HaplotypeCaller (waits for markdup) ---
HC_JID=$(sbatch --parsable --dependency=afterok:${MARK_JID} scripts/05a_haplotypecaller.sh)
echo "haplocaller  → $HC_JID"

# --- Step 5b: gather per-sample GVCFs (waits for HaplotypeCaller) ---
GATHER_JID=$(sbatch --parsable --dependency=afterok:${HC_JID} scripts/05b_gather_gvcf.sh)
echo "gather       → $GATHER_JID"

# --- Step 6: joint genotyping (waits for gather) ---
JOINT_JID=$(sbatch --parsable --dependency=afterok:${GATHER_JID} scripts/06_joint_genotype.sh)
echo "joint_geno   → $JOINT_JID"

# --- Step 7: filter + annotate (waits for joint genotyping) ---
FILT_JID=$(sbatch --parsable --dependency=afterok:${JOINT_JID} scripts/07_filter_annotate.sh)
echo "filter+ann   → $FILT_JID"

# Check that everything is queued
squeue -u $USER
# Steps 3-7 should show state PD with reason (Dependency)

squeue -u $USER                                        # see all your jobs
squeue -u $USER --start                                # see expected start times for pending jobs
sacct -u $USER --starttime=today --format=JobID,JobName,State,Elapsed,MaxRSS,ExitCode

# Latest log from any step (full path so this works from any cwd)
tail -f ~/scratch/reseq/logs/03_align_*.out | head -200

# Specifically task 2 of the align array job:
tail -f ~/scratch/reseq/logs/03_align_${ALIGN_JID}_2.out

# Cancel all your jobs at once
scancel -u $USER

# Or cancel just the downstream ones (if e.g. alignment is still fine)
scancel ${MARK_JID} ${HC_JID} ${GATHER_JID} ${JOINT_JID} ${FILT_JID}

sacct -u $USER --starttime=today \
      --format=JobID,JobName,State,Elapsed,MaxRSS,ReqMem,ReqCPUs