diff --git a/pawsey_shortread/megahit_map_reads.slurm b/pawsey_shortread/megahit_map_reads.slurm new file mode 100644 index 0000000..0f34bf8 --- /dev/null +++ b/pawsey_shortread/megahit_map_reads.slurm @@ -0,0 +1,102 @@ +#!/bin/bash +#SBATCH --job-name=MegaHitMap +#SBATCH --time=0-8 +#SBATCH --ntasks=1 +#SBATCH --cpus-per-task=16 +#SBATCH --mem=32G +#SBATCH -o slurm_output/megahit_slurm/megahit-map-%A_%a.out +#SBATCH -e slurm_output/megahit_slurm/megahit-map-%A_%a.err + +set -euo pipefail +trap 's=$?; echo "$0: Error on line $LINENO: $BASH_COMMAND" >&2; exit $s' ERR + +# Test we have atavide_lite conda installed, and then activate it. +eval "$(conda shell.bash hook)" +"$HOME/GitHubs/atavide_lite/pawsey_lib/check_atavide_lite_env.sh" +ATAVIDE_CONDA="/scratch/$PAWSEY_PROJECT/$USER/software/miniconda3/atavide_lite" +conda activate "$ATAVIDE_CONDA" + +if [[ ! -f R1_reads.txt ]]; then + echo "Please make a file with the R1 reads using this command:" >&2 + echo 'find fastq -name \*R1\* -printf "%f\n" > R1_reads.txt' >&2 + exit 2 +fi + +if [[ ! -f DEFINITIONS.sh ]]; then + echo "Please create a DEFINITIONS.sh file with FILEEND and HOSTREMOVED" >&2 + exit 2 +fi + +source DEFINITIONS.sh +: "${FILEEND:?FILEEND must be defined in DEFINITIONS.sh}" +: "${HOSTREMOVED:?HOSTREMOVED must be defined in DEFINITIONS.sh}" +: "${SLURM_ARRAY_TASK_ID:?Submit this script as a SLURM array job}" +: "${SLURM_CPUS_PER_TASK:?SLURM_CPUS_PER_TASK is not set}" +if (( SLURM_CPUS_PER_TASK < 5 )); then + echo "This job requires at least 5 CPUs, but SLURM_CPUS_PER_TASK=$SLURM_CPUS_PER_TASK" >&2 + exit 2 +fi + +# The mapper, view, and sort stages run concurrently. Account for each +# samtools process's main thread as well as its additional -@ threads. +MAP_THREADS=$((SLURM_CPUS_PER_TASK - 4)) +VIEW_THREADS=1 +SORT_THREADS=1 +INDEX_THREADS=$((SLURM_CPUS_PER_TASK - 1)) + +R1=$(sed -n "${SLURM_ARRAY_TASK_ID}p" R1_reads.txt) +if [[ -z "$R1" ]]; then + echo "No sample was found on line $SLURM_ARRAY_TASK_ID of R1_reads.txt" >&2 + exit 2 +fi +R2=${R1/_R1/_R2} +SAMPLE=${R1%"$FILEEND"} +if [[ "$SAMPLE" == "$R1" ]]; then + echo "ERROR: $R1 does not end with FILEEND ($FILEEND)" >&2 + exit 2 +fi + +ASSEMBLY_DIR="megahit/$SAMPLE" +REFERENCE="$ASSEMBLY_DIR/final.contigs.fa" +BAM="$ASSEMBLY_DIR/mapped_reads.bam" +TMP_BAM="$ASSEMBLY_DIR/.mapped_reads.${SLURM_JOB_ID:-$$}.bam" +TMP_PREFIX="${TMPDIR:-/tmp}/megahit-map-${SLURM_JOB_ID:-$$}-${SLURM_ARRAY_TASK_ID}" + +for INPUT in "$HOSTREMOVED/$R1" "$HOSTREMOVED/$R2" "$REFERENCE"; do + if [[ ! -s "$INPUT" ]]; then + echo "Missing or empty input: $INPUT" >&2 + exit 2 + fi +done + +if [[ -s "$BAM" && -s "$BAM.bai" ]] \ + && samtools quickcheck "$BAM" \ + && samtools idxstats "$BAM" >/dev/null; then + echo "$BAM and its index already exist and passed validation. Nothing to do." >&2 + exit 0 +fi +if [[ -e "$BAM" || -e "$BAM.bai" ]]; then + echo "Incomplete output exists at $BAM or $BAM.bai; refusing to overwrite it" >&2 + exit 2 +fi + +cleanup() { + status=$? + rm -f -- "$TMP_BAM" "$TMP_BAM.bai" + exit "$status" +} +trap cleanup EXIT + +echo "Mapping $HOSTREMOVED/$R1 and $HOSTREMOVED/$R2 to $REFERENCE" >&2 +minimap2 -t "$MAP_THREADS" --secondary=no -ax sr \ + --split-prefix "$TMP_PREFIX" "$REFERENCE" "$HOSTREMOVED/$R1" "$HOSTREMOVED/$R2" \ + | samtools view -@ "$VIEW_THREADS" -b -F 4 - \ + | samtools sort -@ "$SORT_THREADS" -T "${TMP_PREFIX}.sort" -o "$TMP_BAM" - + +samtools index -@ "$INDEX_THREADS" "$TMP_BAM" +samtools quickcheck "$TMP_BAM" +samtools idxstats "$TMP_BAM" >/dev/null +mv -- "$TMP_BAM" "$BAM" +mv -- "$TMP_BAM.bai" "$BAM.bai" + +echo "Created $BAM and $BAM.bai" >&2