Linux for Bioinformatics

Lesson 13 of 14 · 11 min

Writing Your First Bash Script

Throughout this course you have typed commands one at a time. A bash script is a plain text file that stores those commands so you can run the whole sequence again with a single command. In this lesson you will build a script that counts the sequencing reads in every FASTQ file in a folder, letting automation replace repetitive typing.

Build The Script

bash
#!/bin/bash
# run.sh - count the reads in every FASTQ file in a folder

# The folder that holds our FASTQ files
data_dir="fastq"

# Go through each file whose name ends in .fastq
for file in "$data_dir"/*.fastq
do
    # wc -l counts lines; < feeds in the file so we get just a number
    lines=$(wc -l < "$file")

    # Each read in a standard FASTQ file is 4 lines, so divide to get the read count
    reads=$(( lines / 4 ))

    echo "$file has $reads reads"
done

The first line, #!/bin/bash, is called the shebang, and it tells your computer to run the file with the bash program. Any line starting with # is a comment that bash ignores, so you can leave notes for yourself. The line data_dir="fastq" creates a variable named data_dir that stores the folder name, and the quotes keep the value safe if it ever contains spaces.

The line beginning with for starts a loop, where the pattern "$data_dir"/*.fastq expands to every file whose name ends in .fastq and the variable file holds one filename each time around. Everything between do and done runs once per file, so you write the steps a single time and bash repeats them for you.

Inside the loop, wc -l counts lines, the < symbol feeds the file in so you get only a number, and wrapping it in $(...) stores that number in a variable called lines. Because each read in a standard FASTQ file is written as 4 lines, $(( lines / 4 )) divides to turn the line count into a read count. Finally, echo prints a sentence showing each filename and how many reads it holds.

Make It Executable

By default a new text file is not allowed to run as a program. The command chmod +x run.sh changes its permissions to add the executable flag (+x), giving the file permission to run. Then ./run.sh runs it, where the ./ tells bash to look in the current folder for a file named run.sh.

bash
chmod +x run.sh
./run.sh
fastq/sample_A.fastq has 25000 reads
fastq/sample_B.fastq has 18200 reads
fastq/sample_C.fastq has 30450 reads

Try it yourself: Make a folder called fastq, drop two or three .fastq files into it, then save the script above as run.sh. Run chmod +x run.sh followed by ./run.sh and confirm it prints a read count for every file. Now change data_dir to point at a different folder to prove you can reuse the same script without retyping a single command.