You are ready for this if...
You have finished Command-Line Basics and you can move around folders, look inside files, and use grep and the pipe |. If that still feels new, go back and practise first. There is no rush here.
1. Variables and quoting
A variable is just a name for a value. Set it with = (no spaces), and read it back with $.
name="genes.fasta" echo "$name" # prints: genes.fasta echo "I am in $PWD" # $PWD is the current folder
Always quote your variables
Write "$name", not $name. The quotes keep filenames with spaces in one piece. Forgetting them is the single most common beginner bug in shell scripts.
2. Capture a command's output
Command substitution with $( ) runs a command and hands you its result to store or print.
count=$(grep -c ">" genes.fasta) echo "The file has $count sequences"
3. Loops: the same job on many files
This is where the terminal starts to save real time. A for loop repeats a command for every item you give it.
for f in *.fasta; do echo "$f has $(grep -c '>' "$f") sequences" done
Read it as a sentence
"For each file f ending in .fasta, do this, then move to the next." One instruction, applied to every matching file, however many there are.
4. Making decisions: if and test
if [ -f genes.fasta ]; then echo "the file is here" else echo "missing" fi
Handy tests
[ -f file ] is true if a file exists, [ -d folder ] for a directory, and [ -z "$x" ] if a string is empty. Mind the spaces inside the brackets: they are required.
5. Writing a script you can run
Put commands in a .sh file and you have a reusable tool. The first line (the "shebang") tells the system to run it with bash.
# file: count_seqs.sh #!/bin/bash # count sequences in the FASTA file given as the first argument grep -c ">" "$1"
Then make it executable and run it:
chmod +x count_seqs.sh
./count_seqs.sh genes.fasta
What "permissions" mean
chmod +x marks a file as runnable. Those number codes you see (like chmod 755) are just read = 4, write = 2, execute = 1 added up for the owner, group, and everyone else. 755 means the owner can do everything and others can read and run.
6. sed and awk: editing text and columns
Two classic tools that handle most text wrangling. sed does find-and-replace; awk works with columns, which is perfect for table formats like BED, VCF, and GTF.
# sed: strip the leading ">" from FASTA headers grep ">" genes.fasta | sed 's/>//' # awk: print columns 1 and 3 of a table (e.g. a BED file) awk '{print $1, $3}' regions.bed # awk: keep only rows where column 3 is above 1000 awk '$3 > 1000' regions.bed
Why awk is everywhere in bioinformatics
So many formats are just columns of tab-separated data. awk lets you pick columns, filter rows, and do quick maths on them without opening anything. A little awk replaces a lot of manual work.
7. Finding files: find and xargs optional on first pass
find . -name "*.fastq" # every FASTQ in this folder and below find . -name "*.sam" | xargs rm # feed the results to another command (careful: rm!)
Look before you delete
A find ... | xargs rm deletes every match with no undo. Run the find on its own first to see exactly what it lists, then add the rm.
8. Getting data and unpacking it optional on first pass
wget https://example.org/data.tar.gz # download a file (or: curl -O ...) tar -xzf data.tar.gz # unpack a .tar.gz archive gunzip reads.fastq.gz # unzip a single .gz file
When something runs too long
Check what is running with top (press q to quit) or ps, and stop a stuck command with Ctrl-C. Big sequencing files can take a while; that is normal.
Practice exercises
Try each at your own terminal before opening the solution.
Exercise 1: loop over files. Print the name of every .fasta file in the current folder.
Show a worked solution
for f in *.fasta; do echo "$f" done
Exercise 2: count reads in a FASTQ. Store the line count of reads.fastq in a variable, then print the number of reads (lines divided by 4).
Show a worked solution
n=$(wc -l < reads.fastq) echo $(( n / 4 ))
Exercise 3: pick a column. Use awk to print just the second column of a table called table.txt.
Show a worked solution
awk '{print $2}' table.txt
Check yourself
R and the terminal do not run in this page, so here is what correct output looks like. Type each one and compare.
| You type | You should see |
|---|---|
| count=$(grep -c ">" genes.fasta); echo $count | the sequence count |
| echo $(( $(wc -l < reads.fastq) / 4 )) | the number of reads |
| awk '{print $2}' table.txt | the second column |
Check your understanding
awk '{print $2}' table.txt do?$2 is the second field (column) of each line.#!/bin/bash (the "shebang"), do?./script.sh knows to run with bash..fastq file in this folder and all folders below it. Which command?find searches recursively down the folder tree; ls *.fastq only looks in the current folder.for f in *.fasta; do mv $f data/; done fails on a file named my reads.fasta. What is the fix?my reads.fasta into two words. "$f" keeps the filename whole.count=$(grep -c ">" genes.fasta) do?$( ) runs the command inside and hands back its output, which you then store in the variable."$f" with quotes inside a loop?my reads.fastq would be split into two pieces. Quoting keeps it whole.chmod +x script.sh do?+x adds the "execute" permission, which a file needs before the shell will run it as a program.Project: wrangle a real data file from the command line →
Level 3 (Advanced): script functions, safe error handling, parallel jobs, and working on a remote server is coming next.