Working with Sequencing Data in Bash

After Completing This Topic

You'll be able to open FASTA/FASTQ files in the terminal, count sequences, and search for specific sequences. This is your first experience working with data using only the keyboard — no mouse required.

What Is a Terminal?

Just as you use a pipette to precisely handle reagents in the lab, the terminal is a tool for precisely handling data on your computer. Instead of clicking, you type commands.

In Google Colab, you can run Bash commands by prefixing them with !.

Creating and Examining Files

Let's start by creating a practice FASTA file.

# Create a practice FASTA file
echo ">BRCA1_human
ATGGATTTATCTGCTCTTCGCGTTGAAGAAGTACAAAATGTC
>TP53_human
ATGGAGGAGCCGCAGTCAGATCCTAGCGTGAGTTTGCTGTGA
>EGFR_human
ATGCGACCCTCCGGGACGGCCGGGGCAGCGCTCCTGGCGCTG" > genes.fasta

# View the file contents
cat genes.fasta

# Count the number of lines
wc -l genes.fasta

wc stands for Word Count. The -l option counts only lines.

grep: Searching for Sequences

grep is a command that finds specific patterns in a file. It's like searching for keywords in a paper.

# Print only FASTA header (>) lines
echo ">BRCA1_human
ATGGATTTATCTGCTCTTCGCGTTGAAGAAGTACAAAATGTC
>TP53_human
ATGGAGGAGCCGCAGTCAGATCCTAGCGTGAGTTTGCTGTGA
>EGFR_human
ATGCGACCCTCCGGGACGGCCGGGGCAGCGCTCCTGGCGCTG" > genes.fasta

grep ">" genes.fasta

# Count the number of sequences (header lines = number of sequences)
grep -c ">" genes.fasta

The -c option outputs only the count of matched lines.

# Find a specific gene
grep "BRCA1" genes.fasta

Understanding the FASTQ Format

The data from actual sequencing instruments comes in FASTQ format. Unlike FASTA, it includes quality scores.

# Create a FASTQ file (4 lines per set)
echo "@read_001
ATGCGATCGATCGATCGATCG
+
IIIIIIIIIIIIIIIIIIIII
@read_002
GCTAGCTAGCTAGCTAGCTAG
+
IIIIIIIIIIIIIIIIIIIII
@read_003
NNNNNATCGATCGATCGATCG
+
!!!!!IIIIIIIIIIIIIIIII" > reads.fastq

cat reads.fastq

The 4-line structure of FASTQ:

1. Read ID starting with @
2. Nucleotide sequence
3. + (separator)
4. Quality scores (higher is better; I = best, ! = worst)

awk: Processing Data

awk is a tool that processes text data column by column. It's similar to selecting specific columns in a spreadsheet.

# Create tab-delimited data
echo "BRCA1	chr17	43044295	43170245
TP53	chr17	7661779	7687538
EGFR	chr7	55019017	55211628" > gene_locations.tsv

# Print only gene names (column 1)
awk '{print $1}' gene_locations.tsv

# Calculate gene length (column 4 - column 3)
awk '{print $1, $4 - $3, "bp"}' gene_locations.tsv

Pipes (|): Chaining Commands

A pipe (|) connects the output of one command to the input of the next. It's like feeding the product of Reagent A into Reagent B in an experiment.

# Extract headers from FASTA and count them
echo ">BRCA1_human
ATGGATTTATCTGCTCTTCGCGTTGAAGAAGTACAAAATGTC
>TP53_human
ATGGAGGAGCCGCAGTCAGATCCTAGCGTGAGTTTGCTGTGA
>EGFR_human
ATGCGACCCTCCGGGACGGCCGGGGCAGCGCTCCTGGCGCTG" > genes.fasta

grep ">" genes.fasta | wc -l

Try It Yourself (Faded Example)

Fill in the blanks to complete a command that counts the number of reads in a FASTQ file.

In FASTQ, read IDs start with @.

# How do you count lines starting with @?
grep  "" test.fastq

Common Errors & Solutions

Q: I get command not found

In Colab, you need the ! prefix: !grep ">" genes.fasta

Q: Non-ASCII characters are garbled

This is a file encoding issue. Check the encoding with file genes.fasta. If it's not UTF-8, convert it with iconv -f EUC-KR -t UTF-8 input.txt > output.txt.

Q: I get Permission denied

This happens when the file lacks execution permission. Grant it with chmod +x script.sh.

In the next article, we'll learn how to analyze this data more precisely with Python.