Downsample fastqs

Question

I'm working on a project to downsample some fastqs (files that contain sequences).

Each line of the fastq bioinformatics format comprises 4 lines chunks (id, dna sequence, "+", quality score).

Downsampling a fastq is going to select n number of chunks or select x% of chunks. But I don't handle the downsampling: I just pass parameters to a tool called seqtk to handle that.

I'd like some insights on the bash I use to process the options. I especially feel like there's something to have a store_true action for the -t flag.

set -e

usage(){
    echo "Downsample fastqs by calling the seqtk tool
Usage example: sh seqtk_downsampling.sh -n {sample_name} -s {seed_for_random: int} -p {percentage_to_downsample_to: [0-100; 0-1]} -t {testing_flag}"
}


while getopts n:s:p::t:h option
do
    case "${option}"
    in
    n) sample_name=${OPTARG};;
    s) seed=${OPTARG};;
    p) perc_down=${OPTARG};;
    t) testing=${OPTARG};;
    h) usage
       exit
       ;;
    esac
done

# nothing given
if [ -z "$*" ]; then usage; exit 1; fi

# get fqs from sample name, check if they exist
if ! [ -z ${sample_name+x} ]; then
    fqs=$(find . -name "${sample_name}*fastq.gz")

    if [[ -z ${fqs} ]]; then
        current_dir=$(pwd)
        echo "Couldn't find the fqs for ${sample_name} in ${current_dir}"
        exit 1
    fi
else
    echo "No sample name given using -n flag"
    exit 1
fi

# testing mode
if [[ $testing == "true" ]]; then
    percs_down="50 60 70 80 90"
    
    for perc in $percs_down
    do
        for file in $fqs
        do
            echo $file $perc $seed
            # seqtk sample -s ${seed} ${file} ${perc}
        done
    done
    exit 0
else
    echo "Wrong argument given to -t flag"
    exit 1
fi

# given percentage and seed 
if ! [ -z ${perc_down+x} ] && ! [ -z ${seed+x} ]; then
    # format the perc
    if [[ ${perc_down} =~ ^[0][,.][0-9]*$ ]]; then
        perc=${perc_down/,/.}

    elif [[ ${perc_down} -le 100 ]] && [[ ${perc_down} -gt 0 ]]; then
        perc=$(awk -v per=$perc_down 'BEGIN {print per/100}')

    else
        echo "I wantz int between 0 and 100 (or float between 0 and 1)"
        exit 1
    fi

    for file in $fqs
    do
        echo $file $perc $seed
        # seqtk sample -s ${seed} ${file} ${perc}
    done

else
    echo "Missing required arguments -s {seed} or -p {perc_down}"
    exit 1
fi

Answer 1

If you intend to make the script executable, it needs a shebang to select the desired interpreter:

#!/bin/bash

---

I like the use of set -e; you might consider adding -u and -o pipefail.

---

We missed the error when an invalid option is specified:

'?') usage >&2; exit 1 ;;

---

Error messages should go to stream 2 (the standard error stream). Like this:

if [ -z "$*" ]; then usage >&2; exit 1; fi

    echo >&2 "Couldn't find the fqs for ${sample_name} in ${current_dir}"
    exit 1

echo >&2 "Missing required arguments -s {seed} or -p {perc_down}"
exit 1

---

Since we're using Bash, then we can use arrays:

if [[ $testing == "true" ]]; then
    percs_down=(50 60 70 80 90)

    for perc in "${percs_down[@]}"

---

We can test that both $perc_down and $seed are set with a single, positive test:

if [ "${perc_down+x}${seed+x}" = "xx" ]

---

If there's only one acceptable value for -t argument, why isn't it a simple flag argument? OTOH, if we will be adding new values, then case is going to be more suitable than if for selecting the desired code path.