Tutorial

Using Variables

30 Minutes

Objectives

• Read a for loop that runs commands on each file
• Use variables to make a script reusable
• Use variables to improve script robustness
• Use variables to improve script readability

Prerequisite

• The script we are using is introduced in Evaluating Job Efficiency.

The script we will be working from, pipeline.sl, is below:

#!/bin/bash -e

#SBATCH --job-name      pipeline
#SBATCH --account       nesi99991
#SBATCH --time          04:00:00
#SBATCH --mem           5G
#SBATCH --cpus-per-task 4
#SBATCH --profile       task

module purge
module load FastQC/0.12.1
module load Trimmomatic/0.39-Java-1.8.0_144
module load bwa-mem2/2.3-GCC-12.3.0
module load SAMtools/1.22-GCC-12.3.0
module load BCFtools/1.22-GCC-12.3.0

# Limit Java's memory, or Trimmomatic sizes it from the whole node.
export JAVA_TOOL_OPTIONS="-Xmx1g"

# Reference genome indexing
gunzip -kf ref_genome/ecoli_rel606.fasta.gz
bwa-mem2 index ref_genome/ecoli_rel606.fasta

# Quality control
cd untrimmed_fastq
fastqc *.fastq*

mkdir -p ../trimmed
cd ../trimmed
for infile in ../untrimmed_fastq/*_1.fastq.gz
do
    base=$(basename ${infile} _1.fastq.gz)
    trimmomatic PE ${infile} ../untrimmed_fastq/${base}_2.fastq.gz \
        ${base}_1.trim.fastq.gz ${base}_1un.trim.fastq.gz \
        ${base}_2.trim.fastq.gz ${base}_2un.trim.fastq.gz \
        SLIDINGWINDOW:4:20 MINLEN:25 \
        ILLUMINACLIP:${EBROOTTRIMMOMATIC}/adapters/NexteraPE-PE.fa:2:40:15
done

cd ..
mkdir -p results/sam results/bam results/bcf results/vcf
gunzip trimmed/*.fastq.gz

# Alignment and variant calling
for infile in trimmed/*_1.trim.fastq
do
    base=$(basename ${infile} _1.trim.fastq)
    bwa-mem2 mem -t 4 ref_genome/ecoli_rel606.fasta \
        ${infile} trimmed/${base}_2.trim.fastq > results/sam/${base}.aligned.sam
    samtools view -S -b results/sam/${base}.aligned.sam \
        > results/bam/${base}.aligned.bam
    samtools sort results/bam/${base}.aligned.bam \
        -o results/bam/${base}.aligned.sorted.bam
    bcftools mpileup -O b -o results/bcf/${base}.bcf \
        -f ref_genome/ecoli_rel606.fasta results/bam/${base}.aligned.sorted.bam
    bcftools call --ploidy 1 -m -v \
        -o results/vcf/${base}.vcf results/bcf/${base}.bcf
    vcfutils.pl varFilter results/vcf/${base}.vcf > results/vcf/${base}.var.vcf
done

This script already uses variables in the for loop with $infile and $base. The for loop runs its commands once for each file that matches trimmed/*_1.trim.fastq, with $infile set to that file. The line base=$(basename ${infile} _1.trim.fastq) runs the command inside $( ) and stores its output in base. basename removes the directory and the ending, so for trimmed/SRR2584863_1.trim.fastq, $base is SRR2584863. But how can they help us improve our script further?

Using variables helps with several aspects of writing and using code. Variables can help with ensuring our script is readable, reusable and robust.

Formatting variables

You will likely see variables called several different ways: $var, ${var} and "${var}" are all common options. All of these are technically correct, and in some situations you can use whatever option you want. BUT, in some situations how you format your variable calls makes a difference.

In some cases, it can be unclear where your variable ends and where any adjacent string begins. To help with that, we can wrap the variable name in curly braces: ${var}.

ls

01sample.txt  02sample.txt

number=01
wc $numbersample.txt

wc: .txt: No such file or directory

Instead of combining our variable (01) with the rest of the filename (sample), bash has instead used what it has stored in the variable $numbersample which is currently nothing.

wc ${number}sample.txt

0 0 0 01sample.txt

#

Braces in our script

In pipeline.sl, $base is set to SRR2584863. What does this print?

echo $base_1.trim.fastq

.trim.fastq

Bash reads the variable name as base_1, which is not set, so it expands to nothing.

SRR2584863_1.trim.fastq

That needs braces: ${base}_1.trim.fastq, as the script uses.

An error, because base_1 is not set.

Bash carries on silently, unless the script uses set -u (see below).

Using quotations around a call of a variable name can help avoid issues caused by spaces or other special characters in the variable. For example, if we try to use a variable to create a directory with a space in the name, we run into issues that don't immediately raise any errors.

dir="Folder Name"
mkdir $dir
ls

Folder  Name

Instead of creating a directory called 'Folder Name', bash instead interpreted the command as mkdir Folder Name and created two directories: Folder and Name.

mkdir "$dir"
ls

 Folder  'Folder Name'   Name

#

Quotes

What happens when you run these commands?

file="my data.txt"
wc -l $file

wc looks for two files, my and data.txt.

Without quotes, bash splits the value at the space. wc -l "$file" counts the lines in my data.txt.

wc counts the lines in my data.txt.

That needs quotes: wc -l "$file".

Bash stops with an error about the space.

Variables in file paths

Let's use variables to help make the for loop more readable. We will define our filenames once, and then use variables to call the files where they are needed.

For each loop we have the following files, almost all of which are used multiple times:

• $infile: trimmed/${base}_1.trim.fastq and its partner trimmed/${base}_2.trim.fastq
• the reference genome: ref_genome/ecoli_rel606.fasta
• the aligned .sam
• the aligned .bam
• the sorted .bam
• the .bcf
• the .vcf
• the .var.vcf

Rewrite the loop

Before reading on, try rewriting the loop in pipeline.sl so that each of these files has a variable, set once at the start of the loop, and the commands use the variables.

Let's give each of these a simple but unique name, and see how the for loop looks.

for infile in trimmed/*_1.trim.fastq
do
    base=$(basename ${infile} _1.trim.fastq)

    # setting the variables
    ref_genome=ref_genome/ecoli_rel606.fasta
    fq1=trimmed/${base}_1.trim.fastq
    fq2=trimmed/${base}_2.trim.fastq
    sam=results/sam/${base}.aligned.sam
    bam=results/bam/${base}.aligned.bam
    sorted_bam=results/bam/${base}.aligned.sorted.bam
    bcf=results/bcf/${base}.bcf
    vcf=results/vcf/${base}.vcf
    final_variants=results/vcf/${base}.var.vcf

    bwa-mem2 mem -t 4 ${ref_genome} ${fq1} ${fq2} > ${sam}
    samtools view -S -b ${sam} > ${bam}
    samtools sort ${bam} -o ${sorted_bam}
    bcftools mpileup -O b -o ${bcf} -f ${ref_genome} ${sorted_bam}
    bcftools call --ploidy 1 -m -v -o ${vcf} ${bcf}
    vcfutils.pl varFilter ${vcf} > ${final_variants}
done

The quality control loop gets the same treatment. With every path written from the folder the job runs in, the cd lines are no longer needed, and FastQC's reports go to results/fastqc instead of in with the raw data.

pipeline.sl with variables

#!/bin/bash -e

#SBATCH --job-name      pipeline
#SBATCH --account       nesi99991
#SBATCH --time          04:00:00
#SBATCH --mem           5G
#SBATCH --cpus-per-task 4
#SBATCH --profile       task

module purge
module load FastQC/0.12.1
module load Trimmomatic/0.39-Java-1.8.0_144
module load bwa-mem2/2.3-GCC-12.3.0
module load SAMtools/1.22-GCC-12.3.0
module load BCFtools/1.22-GCC-12.3.0

# Limit Java's memory, or Trimmomatic sizes it from the whole node.
export JAVA_TOOL_OPTIONS="-Xmx1g"

adapters=${EBROOTTRIMMOMATIC}/adapters/NexteraPE-PE.fa

# Reference genome indexing
gunzip -kf ref_genome/ecoli_rel606.fasta.gz
bwa-mem2 index ref_genome/ecoli_rel606.fasta

# Quality control
mkdir -p results/fastqc trimmed
for infile in untrimmed_fastq/*_1.fastq.gz
do
    base=$(basename ${infile} _1.fastq.gz)
    raw_2=untrimmed_fastq/${base}_2.fastq.gz
    fastqc --outdir results/fastqc ${infile} ${raw_2}
    trimmomatic PE ${infile} ${raw_2} \
        trimmed/${base}_1.trim.fastq.gz trimmed/${base}_1un.trim.fastq.gz \
        trimmed/${base}_2.trim.fastq.gz trimmed/${base}_2un.trim.fastq.gz \
        SLIDINGWINDOW:4:20 MINLEN:25 ILLUMINACLIP:${adapters}:2:40:15
done

mkdir -p results/sam results/bam results/bcf results/vcf
gunzip -kf trimmed/*.fastq.gz

# Alignment and variant calling
for infile in trimmed/*_1.trim.fastq
do
    base=$(basename ${infile} _1.trim.fastq)

    # setting the variables
    ref_genome=ref_genome/ecoli_rel606.fasta
    fq1=trimmed/${base}_1.trim.fastq
    fq2=trimmed/${base}_2.trim.fastq
    sam=results/sam/${base}.aligned.sam
    bam=results/bam/${base}.aligned.bam
    sorted_bam=results/bam/${base}.aligned.sorted.bam
    bcf=results/bcf/${base}.bcf
    vcf=results/vcf/${base}.vcf
    final_variants=results/vcf/${base}.var.vcf

    bwa-mem2 mem -t 4 ${ref_genome} ${fq1} ${fq2} > ${sam}
    samtools view -S -b ${sam} > ${bam}
    samtools sort ${bam} -o ${sorted_bam}
    bcftools mpileup -O b -o ${bcf} -f ${ref_genome} ${sorted_bam}
    bcftools call --ploidy 1 -m -v -o ${vcf} ${bcf}
    vcfutils.pl varFilter ${vcf} > ${final_variants}
done

Our for loop is definitely longer, but the commands are easier to read because we don't need to sift through the file paths to know what file is being used where. Beyond readability, this also helps reduce errors and make any errors that do occur easier to fix.

One more improvement

One of these variables is set to the same value every time the loop runs. Which one, and where could it go instead?

Solution

ref_genome is the same for every sample, so it can be set once, before the loop starts.

Let's compare the impact of a simple typo in each version of this loop. For example, what if I miss a letter when typing in the path for the .sam file?

Without variables this returns:

[E::hts_open_format] Failed to open file "results/sam/SRR2584863.aligned.sam" : No such file or directory
samtools view: failed to open "results/sam/SRR2584863.aligned.sam" for reading: No such file or directory

Exercise

In the original script, the path for .sam is written out twice. Which call of the .sam file had the typo that created that error?

Answer

The typo was made in the bwa-mem2 call where the .sam file was created, which means the samtools view call requests a file that doesn't exist.

When we use variables, the same mistake in the .sam variable (sam=results/sam/${base}.aliged.sam) does not stop the script: every command uses the same misspelled path, so the file is named wrongly, but consistently. And instead of needing to search for which call of the .sam path is incorrect, there is only one place that the path is defined.

Unset variables

A different typo is harder to spot: misspelling the name of a variable. If a script uses ${smaple} instead of ${sample}, bash quietly replaces it with nothing, and the command uses the wrong path without any error.

Adding set -u near the top of a script makes it stop as soon as it uses a variable that has not been set, with an error such as smaple: unbound variable.

set -u

The -e in #!/bin/bash -e does the same for commands: the script stops at the first command that fails, instead of carrying on without the files that command should have made.

#

Unset variables

A script starts with set -u, then runs mkdir "${resluts}/vcf", misspelling results. What happens?

The script stops with the error resluts: unbound variable.

It tries to create /vcf, at the top of the filesystem.

That is what happens without set -u, as the misspelled variable expands to nothing.

It creates results/vcf.

Environment variables

In addition to defining variables ourselves, there are also variables that are defined by the operating system and software. We can use these to make scripts flexible for use by other users or in other contexts. Rather than having to update your script to run in a new environment etc, we can have the script adapt to the new environment automatically.

Matching threads to CPUs

In this case, we are going to make sure that the software we are running knows how many CPUs are available to it. Currently, we pass the argument -t 4 to the command bwa-mem2 mem. This argument tells bwa-mem2 how many threads it can use, and we want it to be able to use as many threads as are available to it. Currently, if we change the number of CPUs we are requesting, we need to change this argument to see any difference in the behaviour of bwa-mem2. But Slurm has an environment variable called $SLURM_CPUS_PER_TASK that can be used to pass this information programmatically. We can update our for loop to use this variable:

for infile in trimmed/*_1.trim.fastq
do
    base=$(basename ${infile} _1.trim.fastq)

    # setting the variables
    fq1=trimmed/${base}_1.trim.fastq
    fq2=trimmed/${base}_2.trim.fastq
    sam=results/sam/${base}.aligned.sam
    bam=results/bam/${base}.aligned.bam
    sorted_bam=results/bam/${base}.aligned.sorted.bam
    bcf=results/bcf/${base}.bcf
    vcf=results/vcf/${base}.vcf
    final_variants=results/vcf/${base}.var.vcf

    bwa-mem2 mem -t ${SLURM_CPUS_PER_TASK} ${ref_genome} ${fq1} ${fq2} > ${sam}
    samtools view -S -b ${sam} > ${bam}
    samtools sort ${bam} -o ${sorted_bam}
    bcftools mpileup -O b -o ${bcf} -f ${ref_genome} ${sorted_bam}
    bcftools call --ploidy 1 -m -v -o ${vcf} ${bcf}
    vcfutils.pl varFilter ${vcf} > ${final_variants}
done

pipeline.sl so far, with set -u and the reference genome set once

#!/bin/bash -e

#SBATCH --job-name      pipeline
#SBATCH --account       nesi99991
#SBATCH --time          04:00:00
#SBATCH --mem           5G
#SBATCH --cpus-per-task 4
#SBATCH --profile       task

# Stop if a variable is used before it is set.
set -u

module purge
module load FastQC/0.12.1
module load Trimmomatic/0.39-Java-1.8.0_144
module load bwa-mem2/2.3-GCC-12.3.0
module load SAMtools/1.22-GCC-12.3.0
module load BCFtools/1.22-GCC-12.3.0

# Limit Java's memory, or Trimmomatic sizes it from the whole node.
export JAVA_TOOL_OPTIONS="-Xmx1g"

adapters=${EBROOTTRIMMOMATIC}/adapters/NexteraPE-PE.fa
ref_genome=ref_genome/ecoli_rel606.fasta

# Reference genome indexing
gunzip -kf ref_genome/ecoli_rel606.fasta.gz
bwa-mem2 index ref_genome/ecoli_rel606.fasta

# Quality control
mkdir -p results/fastqc trimmed
for infile in untrimmed_fastq/*_1.fastq.gz
do
    base=$(basename ${infile} _1.fastq.gz)
    raw_2=untrimmed_fastq/${base}_2.fastq.gz
    fastqc --outdir results/fastqc ${infile} ${raw_2}
    trimmomatic PE ${infile} ${raw_2} \
        trimmed/${base}_1.trim.fastq.gz trimmed/${base}_1un.trim.fastq.gz \
        trimmed/${base}_2.trim.fastq.gz trimmed/${base}_2un.trim.fastq.gz \
        SLIDINGWINDOW:4:20 MINLEN:25 ILLUMINACLIP:${adapters}:2:40:15
done

mkdir -p results/sam results/bam results/bcf results/vcf
gunzip -kf trimmed/*.fastq.gz

# Alignment and variant calling
for infile in trimmed/*_1.trim.fastq
do
    base=$(basename ${infile} _1.trim.fastq)

    # setting the variables
    fq1=trimmed/${base}_1.trim.fastq
    fq2=trimmed/${base}_2.trim.fastq
    sam=results/sam/${base}.aligned.sam
    bam=results/bam/${base}.aligned.bam
    sorted_bam=results/bam/${base}.aligned.sorted.bam
    bcf=results/bcf/${base}.bcf
    vcf=results/vcf/${base}.vcf
    final_variants=results/vcf/${base}.var.vcf

    bwa-mem2 mem -t ${SLURM_CPUS_PER_TASK} ${ref_genome} ${fq1} ${fq2} > ${sam}
    samtools view -S -b ${sam} > ${bam}
    samtools sort ${bam} -o ${sorted_bam}
    bcftools mpileup -O b -o ${bcf} -f ${ref_genome} ${sorted_bam}
    bcftools call --ploidy 1 -m -v -o ${vcf} ${bcf}
    vcfutils.pl varFilter ${vcf} > ${final_variants}
done

Now if we request more CPUs for this job, bwa-mem2 will have access to all the CPUs allocated. Slurm only sets $SLURM_CPUS_PER_TASK when the job requests --cpus-per-task, as this script does. Slurm sets its variables only inside a job, so they are empty if you try them on the login node.

#

Where Slurm variables exist

You run echo ${SLURM_CPUS_PER_TASK} on the login node. What do you see?

An empty line.

Slurm only sets this variable inside a job, and only when the job requests --cpus-per-task.

The number of CPUs on the login node.

The number of CPUs your last job requested.

Solo learners

Replace nesi99991 in the --account line with your project code.

Check that it works

Submit the updated script with sbatch pipeline.sl, and once it has finished, check it with seff JOB_ID. Then submit it again with sbatch --cpus-per-task 2 pipeline.sl, and compare the elapsed time.

Replace JOB_ID with the job ID that sbatch printed.

Any #SBATCH option can also be given on the sbatch command line, before the script name. Options on the command line override the same option in the script, which is how sbatch --cpus-per-task 2 pipeline.sl runs the script with 2 CPUs, without editing it.

#

Command line or script

pipeline.sl has #SBATCH --mem 5G. You submit it with sbatch --mem 8G pipeline.sl. How much memory does the job get?

8G, because options on the command line override the script.

5G, because the script is read last.

13G, because the two are added together.

Which other commands could use the CPUs?

bwa-mem2 is not the only command in the loop that can use more than one CPU. Check the documentation for samtools sort: how could it use the CPUs too?

Solution

samtools sort has a -@ option for the number of threads. Each thread uses up to 768 MB of memory by default, so also limit that with -m to stay within the job's memory:

samtools sort -@ ${SLURM_CPUS_PER_TASK} -m 200M ${bam} -o ${sorted_bam}

Other useful variables

There are lots of environment variables that can be useful in scripts:

Variable name	Meaning	Potential use
$USER	The username of the current user	Paths that work for everyone, such as /nesi/nobackup/nesi99991/$USER
$HOME	The path to the user's home directory	Finding your own configuration files or scripts
$PWD	Present Working Directory	Recording where a job ran
$SHELL	The path to the current shell interpreter	Checking which shell a script runs in
$PATH	Directories searched for executable commands	Checking which programs a loaded module has made available
$HOSTNAME	The network name of the machine (within a job this will be the compute node being used)	Recording which node a job ran on

Slurm variables

Slurm has some specific environment variables that it sets when running a job as well.

Variable name	Meaning	Potential use
$SLURM_JOB_ID	ID for the current job	Useful for naming output files that won't clash.
$SLURM_JOB_NAME	Name of the job	Could be used for creating files/directories with useful names
$SLURM_ARRAY_TASK_ID	The current index of your job array	Helpful when writing job arrays
$SLURM_CPUS_PER_TASK	Number of CPUs	Useful as an input for multithreaded functions.
$SLURM_NTASKS	Number of tasks	Useful as an input for MPI functions.
$SLURM_SUBMIT_DIR	Directory where sbatch was called.	Helpful for being able to use relative paths

Variables in #SBATCH lines

#SBATCH lines are read by Slurm before the job starts, not by bash, so variables in them are not replaced. #SBATCH --output ${SLURM_JOB_ID}.out creates a file literally named ${SLURM_JOB_ID}.out. Use tokens instead, such as %j for the job ID and %x for the job name, as in Output files.

Finding out about more useful environment variables

You can't use environment variables that you don't know about! So how can we learn about new environment variables that might help us?

One option is the classic advice to RTFM (read the frickin manual). Lots of documentation includes information about environment variables that are set or used by the software/tool. For example, the samtools documentation includes a section called 'Environment variables' which includes paths to reference files and plugins that can be used.

Another option is to look at what environment variables exist in your current environment. Running the command env will print out the details of your environment, and can be a little overwhelming. You can combine env with the command grep to search and narrow things down. For example, after loading the SAMtools module, we can use env | grep SAM to only get results that include the string SAM.

Output files

pipeline.sl doesn't say where to put its output, so Slurm writes everything to one file, slurm-JOB_ID.out, in the folder the job was submitted from. With many jobs, these files soon pile up and are hard to tell apart.

Adding --output to the header puts each job's output in log/, named after the job name (%x) and job ID (%j):

#SBATCH --output        log/%x_%j.out

For pipeline.sl, that is log/pipeline_JOB_ID.out. Slurm doesn't create the log folder, so create it before submitting with mkdir -p log. From here on, pipeline.sl has this line.

pipeline.sl still runs every step, for every sample, in one job that holds 4 CPUs and 5 GB the whole time. Job Dependencies splits it into a job for each step, each with its own resources.

Keypoints

• Using variables can help make code more readable and robust to errors
• Using environment variables can help make your code more portable
• Braces and quotes make variables expand the way you expect
• set -u stops a script when it uses a variable that is not set
• A for loop runs the same commands on each file
• Options on the sbatch command line override the same option in the script
• Variables do not work in #SBATCH lines. Use tokens like %x and %j

What next

• Next lesson: Debugging Failed Jobs.
• For reference: the Batch Computing Guide and the Slurm Reference Sheet.
• For more shell practice: the Software Carpentry Unix Shell lessons on loops and shell scripts.

Quiz Progress

• Answered: 0 / 0 (0%)
• Correct: 0 / 0 (0%) Reset