Essential Shell Scripts for Bioinformatics & HPC Automation

#!/usr/bin/env bash
#SBATCH --job-name=omics_array
#SBATCH --output=logs/array_%A_%a.out
#SBATCH --error=logs/array_%A_%a.err
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=8
#SBATCH --mem=32G
#SBATCH --time=12:00:00
#SBATCH --partition=standard
#SBATCH --array=1-50%10

set -euo pipefail

MANIFEST="${1:-samples.tsv}"
if [[ ! -f "$MANIFEST" ]]; then
    echo "ERROR: Manifest file '$MANIFEST' not found!" >&2
    exit 1
fi

LINE=$(sed -n "${SLURM_ARRAY_TASK_ID}p" "$MANIFEST")
SAMPLE=$(echo "$LINE" | cut -f1)
R1=$(echo "$LINE" | cut -f2)
R2=$(echo "$LINE" | cut -f3)

SCRATCH="${TMPDIR:-/tmp}/${SAMPLE}_${SLURM_JOB_ID}"
mkdir -p "$SCRATCH" logs results/"$SAMPLE"

cleanup() {
    echo "Cleaning up node scratch space: $SCRATCH"
    rm -rf "$SCRATCH"
}
trap cleanup EXIT INT TERM

echo "=========================================="
echo "Job ID: $SLURM_JOB_ID | Task ID: $SLURM_ARRAY_TASK_ID"
echo "Processing Sample: $SAMPLE on host: $(hostname)"
echo "R1: $R1"
echo "R2: $R2"
echo "=========================================="

cp "$R1" "$SCRATCH/in_R1.fq.gz"
cp "$R2" "$SCRATCH/in_R2.fq.gz"

fastp -i "$SCRATCH/in_R1.fq.gz" -I "$SCRATCH/in_R2.fq.gz" \
      -o "$SCRATCH/${SAMPLE}_clean_R1.fq.gz" -O "$SCRATCH/${SAMPLE}_clean_R2.fq.gz" \
      --thread "$SLURM_CPUS_PER_TASK" \
      --html "results/${SAMPLE}/${SAMPLE}_fastp.html" \
      --json "results/${SAMPLE}/${SAMPLE}_fastp.json"

cp "$SCRATCH/${SAMPLE}_clean_R1.fq.gz" results/"$SAMPLE"/
cp "$SCRATCH/${SAMPLE}_clean_R2.fq.gz" results/"$SAMPLE"/

echo "Sample $SAMPLE processing successfully completed!"

#!/usr/bin/env bash
#SBATCH --job-name=deepvariant_gpu
#SBATCH --output=logs/dv_gpu_%j.out
#SBATCH --error=logs/dv_gpu_%j.err
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=16
#SBATCH --mem=64G
#SBATCH --gres=gpu:1
#SBATCH --time=04:00:00
#SBATCH --partition=gpu

set -euo pipefail

BAM="${1:-input.bam}"
REF="${2:-reference.fa}"
OUTDIR="${3:-output_dv}"

mkdir -p "$OUTDIR" logs

SIF_IMAGE="/shared/containers/deepvariant_gpu_1.6.0.sif"
if [[ ! -f "$SIF_IMAGE" ]]; then
    echo "Pulling DeepVariant GPU singularity image..."
    apptainer pull "$SIF_IMAGE" docker://google/deepvariant:1.6.0-gpu
fi

echo "CUDA Device: $CUDA_VISIBLE_DEVICES"
nvidia-smi

apptainer exec --nv \
    -B /usr/lib/locale/:/usr/lib/locale/ \
    -B "$(pwd)":/data \
    "$SIF_IMAGE" \
    /opt/deepvariant/bin/run_deepvariant \
    --model_type=WGS \
    --ref="/data/$REF" \
    --reads="/data/$BAM" \
    --output_vcf="/data/$OUTDIR/output.vcf.gz" \
    --output_gvcf="/data/$OUTDIR/output.g.vcf.gz" \
    --num_shards="$SLURM_CPUS_PER_TASK"

echo "DeepVariant GPU calling finished successfully!"

#!/usr/bin/env bash
# Real-time process memory & CPU profiler for HPC tasks
set -euo pipefail

TARGET_PID="${1:-}"
LOG_FILE="${2:-resource_profile.tsv}"
INTERVAL="${3:-10}"

if [[ -z "$TARGET_PID" ]]; then
    echo "Usage: $0 <PID> [log_file.tsv] [interval_sec]" >&2
    exit 1
fi

echo -e "Timestamp\tElapsed_Sec\tPID\tCPU_Pct\tRSS_MB\tVMS_MB" > "$LOG_FILE"
START_TIME=$(date +%s)

echo "Monitoring PID $TARGET_PID every ${INTERVAL}s -> $LOG_FILE"

while kill -0 "$TARGET_PID" 2>/dev/null; do
    NOW=$(date +%s)
    ELAPSED=$(( NOW - START_TIME ))
    TIMESTAMP=$(date +"%Y-%m-%d %H:%M:%S")
    
    METRICS=$(ps -p "$TARGET_PID" -o %cpu,rss,vsz --no-headers 2>/dev/null || echo "0 0 0")
    CPU=$(echo "$METRICS" | awk '{print $1}')
    RSS_KB=$(echo "$METRICS" | awk '{print $2}')
    VSZ_KB=$(echo "$METRICS" | awk '{print $3}')
    
    RSS_MB=$(awk "BEGIN {printf \"%.2f\", $RSS_KB/1024}")
    VSZ_MB=$(awk "BEGIN {printf \"%.2f\", $VSZ_KB/1024}")
    
    echo -e "${TIMESTAMP}\t${ELAPSED}\t${TARGET_PID}\t${CPU}\t${RSS_MB}\t${VSZ_MB}" >> "$LOG_FILE"
    sleep "$INTERVAL"
done

echo "Process $TARGET_PID terminated. Monitoring complete."

#!/usr/bin/env bash
# Provision Ubuntu / Debian workstation for computational biology
set -euo pipefail

if [[ $EUID -ne 0 ]]; then
   echo "ERROR: This script must be run as root (e.g. sudo bash $0)" >&2
   exit 1
fi

echo "Updating APT package indices..."
apt-get update -y && apt-get upgrade -y

echo "Installing essential build toolchains & core libraries..."
apt-get install -y --no-install-recommends \
    build-essential \
    cmake \
    git \
    curl \
    wget \
    pigz \
    pbzip2 \
    htop \
    tree \
    jq \
    aria2 \
    pkg-config \
    zlib1g \
    zlib1g-dev \
    libbz2-dev \
    liblzma-dev \
    libcurl4-openssl-dev \
    libssl-dev \
    libncurses5-dev \
    libncursesw5-dev \
    libdeflate-dev

echo "Installing Python 3 pip & virtual environment modules..."
apt-get install -y python3-pip python3-dev python3-venv

echo "Cleaning up APT caches..."
apt-get autoremove -y && apt-get clean

echo "======================================================"
echo "Workstation core libraries provisioned successfully!"
echo "GCC Version:    $(gcc --version | head -n1)"
echo "CMake Version:  $(cmake --version | head -n1)"
echo "Python Version: $(python3 --version)"
echo "======================================================"

#!/usr/bin/env bash
# Non-interactive Miniforge3 & Bioconda environment installer
set -euo pipefail

INSTALL_DIR="${HOME}/miniforge3"
MINIFORGE_URL="https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-Linux-x86_64.sh"

if [[ -d "$INSTALL_DIR" ]]; then
    echo "Miniforge is already installed at $INSTALL_DIR"
else
    echo "Downloading Miniforge installer..."
    curl -fsSL "$MINIFORGE_URL" -o /tmp/miniforge.sh
    
    echo "Installing Miniforge non-interactively to $INSTALL_DIR..."
    bash /tmp/miniforge.sh -b -p "$INSTALL_DIR"
    rm -f /tmp/miniforge.sh
    
    echo "Initializing shell profiles..."
    "$INSTALL_DIR/bin/conda" init bash
fi

export PATH="$INSTALL_DIR/bin:$PATH"

echo "Configuring Bioconda channels with strict priority..."
conda config --add channels defaults
conda config --add channels bioconda
conda config --add channels conda-forge
conda config --set channel_priority strict

echo "Creating core genomics environment..."
mamba create -y -n omics-core \
    samtools \
    bcftools \
    bedtools \
    seqkit \
    fastp \
    multiqc

echo "=========================================================="
echo "Installation complete! To activate your environment:"
echo "source ~/.bashrc && conda activate omics-core"
echo "=========================================================="

#!/usr/bin/env bash
# Generate publication-grade reproducible Conda environment lockfiles
set -euo pipefail

ENV_NAME="${1:-}"
OUTDIR="${2:-.}"

if [[ -z "$ENV_NAME" ]]; then
    echo "Usage: $0 <environment_name> [output_directory]" >&2
    exit 1
fi

mkdir -p "$OUTDIR"
DATE=$(date +%Y%m%d)

echo "Exporting portable cross-platform YAML (no builds)..."
conda env export -n "$ENV_NAME" --no-builds | grep -v "prefix:" > "${OUTDIR}/${ENV_NAME}_portable_${DATE}.yml"

echo "Exporting explicit URL & hash lockfile..."
conda list -n "$ENV_NAME" --explicit > "${OUTDIR}/${ENV_NAME}_explicit_${DATE}.txt"

echo "Lockfiles generated:"
ls -lh "${OUTDIR}/${ENV_NAME}"*
echo "To recreate on another system:"
echo "conda env create -f ${OUTDIR}/${ENV_NAME}_portable_${DATE}.yml"

#!/usr/bin/env bash
# Batch paired-end FASTQ trimming & MultiQC aggregator
set -euo pipefail

RAW_DIR="${1:-raw_data}"
OUT_DIR="${2:-clean_data}"
THREADS="${3:-8}"

mkdir -p "$OUT_DIR" "$OUT_DIR/reports"

echo "Scanning $RAW_DIR for paired-end FASTQs..."

for R1 in "$RAW_DIR"/*_R1*.fastq.gz "$RAW_DIR"/*_1.fastq.gz; do
    [[ -e "$R1" ]] || continue
    
    # Derive corresponding R2 path
    if [[ "$R1" == *"_R1"* ]]; then
        R2="${R1/_R1/_R2}"
    else
        R2="${R1/_1.fastq.gz/_2.fastq.gz}"
    fi
    
    if [[ ! -f "$R2" ]]; then
        echo "WARNING: Pair for $R1 not found, skipping." >&2
        continue
    fi
    
    SAMPLE=$(basename "$R1" | sed -E 's/(_R1|_1).*//')
    echo "Processing Sample: $SAMPLE..."
    
    fastp \
        --in1 "$R1" \
        --in2 "$R2" \
        --out1 "$OUT_DIR/${SAMPLE}_clean_R1.fq.gz" \
        --out2 "$OUT_DIR/${SAMPLE}_clean_R2.fq.gz" \
        --thread "$THREADS" \
        --detect_adapter_for_pe \
        --qualified_quality_phred 20 \
        --length_required 35 \
        --html "$OUT_DIR/reports/${SAMPLE}_fastp.html" \
        --json "$OUT_DIR/reports/${SAMPLE}_fastp.json"
done

echo "Compiling MultiQC aggregated quality summary..."
multiqc "$OUT_DIR/reports" -o "$OUT_DIR/multiqc_report"

echo "QC and Trimming Pipeline Complete! View report at: $OUT_DIR/multiqc_report/multiqc_report.html"

#!/usr/bin/env bash
# Parallel NCBI SRA fasterq-dump batch downloader
set -euo pipefail

ACC_LIST="${1:-srr_list.txt}"
DEST_DIR="${2:-raw_fastq}"
THREADS="${3:-8}"

if [[ ! -f "$ACC_LIST" ]]; then
    echo "ERROR: Accession list '$ACC_LIST' not found!" >&2
    exit 1
fi

mkdir -p "$DEST_DIR"

while IFS= read -r ACC || [[ -n "$ACC" ]]; do
    # Skip comments and empty lines
    [[ "$ACC" =~ ^#.*$ || -z "$ACC" ]] && continue
    
    echo "=========================================="
    echo "Downloading SRA Run: $ACC..."
    echo "=========================================="
    
    fasterq-dump "$ACC" \
        --outdir "$DEST_DIR" \
        --split-files \
        --threads "$THREADS" \
        --temp "/tmp" \
        --progress
        
    echo "Compressing $ACC FASTQs with pigz..."
    pigz -p "$THREADS" "$DEST_DIR/${ACC}"*.fastq
    
    echo "$ACC successfully downloaded & compressed."
done < "$ACC_LIST"

echo "All accessions downloaded successfully!"

#!/usr/bin/env bash
# Deterministic paired-end FASTQ subsampling script
set -euo pipefail

R1="${1:-}"
R2="${2:-}"
NUM_READS="${3:-100000}"
OUT_R1="${4:-subsampled_R1.fq.gz}"
OUT_R2="${5:-subsampled_R2.fq.gz}"
SEED="${6:-100}"

if [[ -z "$R1" || -z "$R2" ]]; then
    echo "Usage: $0 <in_R1.fq.gz> <in_R2.fq.gz> [num_reads] [out_R1.fq.gz] [out_R2.fq.gz] [seed]" >&2
    exit 1
fi

echo "Subsampling $NUM_READS pairs with random seed $SEED..."

seqtk sample -s "$SEED" "$R1" "$NUM_READS" | gzip > "$OUT_R1"
seqtk sample -s "$SEED" "$R2" "$NUM_READS" | gzip > "$OUT_R2"

echo "Verifying output read counts..."
C1=$(zcat "$OUT_R1" | awk '{s++}END{print s/4}')
C2=$(zcat "$OUT_R2" | awk '{s++}END{print s/4}')

echo "R1 reads: $C1 | R2 reads: $C2"
if [[ "$C1" -eq "$C2" ]]; then
    echo "SUCCESS: Read synchronization confirmed ($C1 pairs)!"
else
    echo "ERROR: Read count mismatch ($C1 vs $C2)!" >&2
    exit 1
fi

#!/usr/bin/env bash
# Automated coordinate sorting, indexing, and QC statistics for BAM files
set -euo pipefail

IN_DIR="${1:-raw_bams}"
OUT_DIR="${2:-sorted_bams}"
THREADS="${3:-8}"

mkdir -p "$OUT_DIR" "$OUT_DIR/stats"

for IN_BAM in "$IN_DIR"/*.bam "$IN_DIR"/*.sam; do
    [[ -e "$IN_BAM" ]] || continue
    
    BASE=$(basename "$IN_BAM" | sed -E 's/\.(bam|sam)$//')
    OUT_BAM="$OUT_DIR/${BASE}.sorted.bam"
    
    echo "Processing $BASE..."
    samtools sort -@ "$THREADS" -o "$OUT_BAM" "$IN_BAM"
    samtools index -@ "$THREADS" "$OUT_BAM"
    
    samtools flagstat -@ "$THREADS" "$OUT_BAM" > "$OUT_DIR/stats/${BASE}.flagstat.txt"
    samtools idxstats "$OUT_BAM" > "$OUT_DIR/stats/${BASE}.idxstats.tsv"
    
    echo "Finished $BASE -> $OUT_BAM"
done

echo "Batch BAM processing complete!"

#!/usr/bin/env bash
# Calculate target region depth & coverage fractions across BED intervals
set -euo pipefail

BAM="${1:-}"
BED="${2:-}"
OUT_FILE="${3:-coverage_summary.tsv}"

if [[ -z "$BAM" || -z "$BED" ]]; then
    echo "Usage: $0 <sample.bam> <targets.bed> [out_summary.tsv]" >&2
    exit 1
fi

echo "Computing per-base coverage on $BAM over $BED..."

# Extract depths: col1=chrom, col2=pos, col3=depth
samtools depth -a -b "$BED" "$BAM" | awk -v out="$OUT_FILE" '
BEGIN {
    total_bases=0; sum_depth=0; c1=0; c10=0; c20=0; c30=0;
}
{
    depth = $3;
    total_bases++;
    sum_depth += depth;
    if (depth >= 1)  c1++;
    if (depth >= 10) c10++;
    if (depth >= 20) c20++;
    if (depth >= 30) c30++;
}
END {
    if (total_bases > 0) {
        mean = sum_depth / total_bases;
        pct1  = (c1 / total_bases) * 100;
        pct10 = (c10 / total_bases) * 100;
        pct20 = (c20 / total_bases) * 100;
        pct30 = (c30 / total_bases) * 100;
        
        printf "Metric\tValue\n" > out;
        printf "Target_Bases\t%d\n", total_bases >> out;
        printf "Mean_Depth\t%.2fx\n", mean >> out;
        printf "Pct_Bases_ge_1x\t%.2f%%\n", pct1 >> out;
        printf "Pct_Bases_ge_10x\t%.2f%%\n", pct10 >> out;
        printf "Pct_Bases_ge_20x\t%.2f%%\n", pct20 >> out;
        printf "Pct_Bases_ge_30x\t%.2f%%\n", pct30 >> out;
    }
}'

cat "$OUT_FILE"

#!/usr/bin/env bash
# Extract unmapped read pairs from BAM into paired FASTQ
set -euo pipefail

BAM="${1:-}"
OUT_R1="${2:-unmapped_R1.fq.gz}"
OUT_R2="${3:-unmapped_R2.fq.gz}"
THREADS="${4:-4}"

if [[ -z "$BAM" ]]; then
    echo "Usage: $0 <in.bam> [out_R1.fq.gz] [out_R2.fq.gz] [threads]" >&2
    exit 1
fi

echo "Extracting unmapped paired reads (flag 12 = read & mate unmapped)..."

samtools view -u -f 12 -@ "$THREADS" "$BAM" \
    | samtools collate -u -O -@ "$THREADS" - \
    | samtools fastq -@ "$THREADS" \
        -1 "$OUT_R1" \
        -2 "$OUT_R2" \
        -0 /dev/null \
        -s /dev/null -

echo "Unmapped read pairs written to $OUT_R1 and $OUT_R2"

#!/usr/bin/env bash
# Normalize, split multiallelics, filter PASS variants, and tabix index
set -euo pipefail

IN_VCF="${1:-}"
REF_FA="${2:-}"
OUT_VCF="${3:-filtered_normalized.vcf.gz}"
MIN_QUAL="${4:-30}"
MIN_DP="${5:-10}"

if [[ -z "$IN_VCF" || -z "$REF_FA" ]]; then
    echo "Usage: $0 <raw.vcf.gz> <ref.fa> [out.vcf.gz] [min_qual] [min_dp]" >&2
    exit 1
fi

echo "Step 1: Normalizing & splitting multiallelic variants..."
bcftools norm -m -any -f "$REF_FA" "$IN_VCF" -Ou \
    | bcftools view -i "FILTER='PASS' && QUAL>=${MIN_QUAL} && INFO/DP>=${MIN_DP}" -Oz -o "$OUT_VCF"

echo "Step 2: Building tabix index..."
tabix -p vcf "$OUT_VCF"

echo "Step 3: Verification stats:"
bcftools stats "$OUT_VCF" | grep "^SN" | cut -f 3,4

echo "Standardized VCF created at: $OUT_VCF"

#!/usr/bin/env bash
# Generate comprehensive Ti/Tv and variant metric summary
set -euo pipefail

VCF="${1:-}"
REPORT="${2:-vcf_qc_summary.txt}"

if [[ -z "$VCF" ]]; then
    echo "Usage: $0 <sample.vcf.gz> [out_report.txt]" >&2
    exit 1
fi

TMP_STATS=$(mktemp)
trap 'rm -f "$TMP_STATS"' EXIT

echo "Generating stats for $VCF..."
bcftools stats "$VCF" > "$TMP_STATS"

TITV=$(grep "^TSTV" "$TMP_STATS" | cut -f 5)
SNPS=$(grep "^SN.*number of SNPs:" "$TMP_STATS" | cut -f 4)
INDELS=$(grep "^SN.*number of indels:" "$TMP_STATS" | cut -f 4)
MULTIALLELIC=$(grep "^SN.*number of multiallelic sites:" "$TMP_STATS" | cut -f 4)

cat <<EOF | tee "$REPORT"
======================================================
           VCF CALLSET QUALITY BENCHMARKS
======================================================
VCF File:               $VCF
Total SNVs:             $SNPS
Total Indels:           $INDELS
Multiallelic Sites:     $MULTIALLELIC
Overall Ti/Tv Ratio:    $TITV

Expected Ti/Tv Reference Benchmarks:
  - Whole Genome (WGS):  ~2.0 - 2.1
  - Whole Exome (WES):   ~2.8 - 3.0
======================================================
EOF

#!/usr/bin/env bash
# HPC scratch space audit and orphan file cleanup guard
set -euo pipefail

TARGET_DIR="${1:-$SCRATCH}"
DAYS_OLD="${2:-30}"

if [[ ! -d "$TARGET_DIR" ]]; then
    echo "ERROR: Directory '$TARGET_DIR' does not exist!" >&2
    exit 1
fi

echo "=========================================="
echo "HPC SCRATCH DISK AUDIT: $TARGET_DIR"
echo "=========================================="

echo "Top 10 Largest Directories:"
du -h --max-depth=2 "$TARGET_DIR" 2>/dev/null | sort -hr | head -n 10

echo -e "\nFinding intermediate temporary artifacts older than $DAYS_OLD days..."
find "$TARGET_DIR" -type f \( -name "*.tmp.bam" -o -name "*.sam" -o -name "*_temp*" \) -mtime +"$DAYS_OLD" -print -delete

echo -e "\nFinding empty directories..."
find "$TARGET_DIR" -type d -empty -delete

echo "Scratch audit and cleanup completed!"