A Small Annotation Audit

A collaborator has sent us this project:

rnaseq_project/

├── annotation/
│   ├── gencode_subset.gtf
│   ├── gencode_subset.gtf.backup
│   └── notes.txt
│
├── metadata/
│   └── samples.tsv
│
└── results/

They ask:

Before we use the annotation, can you give me a quick QC report?

Questions for This Audit

1. Which annotation file do we have?

2. Which records are protein-coding genes?

3. Which of those genes are on the autosomes?

4. Which gene records are duplicated?

5. What are the chromosome, coordinates, and genomic span of each gene?

6. Can we save the results as a clean report?

What Would We Need?

→ save the result

The Goal

By the end of this lecture, we should be able to turn a messy annotation directory into a reproducible QC report using ordinary command-line tools.

Why Bioinformatics Data Works Well With CLI

The majority of genomic data is represented as plain text.

The presentation often follows:

characters → lines → fields → patterns

Example: A GTF Record

chr17   HAVANA  gene        7661779  7687550  .  +  .  gene_id "TP53";
chr17   HAVANA  transcript  7661779  7687550  .  +  .  gene_id "TP53"; transcript_id "ENST...";
chr17   HAVANA  exon        7661779  7662590  .  +  .  gene_id "TP53"; transcript_id "ENST...";

One line is one record.

Each record contains multiple fields.

The fields are separated by a delimiter.

Here the main delimiter is a tab:

        ↓      ↓      ↓
chr17<TAB>HAVANA<TAB>gene<TAB>...

Four Concepts to Remember

chr17   HAVANA   gene   ...

chr17
HAVANA
gene
7661779

TP53
gene
chr17

TAB
,
space
:
|

Mini Exercise – Read the Record

Given:

chr17   HAVANA   gene   7661779   7687550   .   +   .   gene_id "TP53";

Answer:

1. How many tab-separated fields are visible?
2. Which field tells us the chromosome?
3. Which field tells us the feature type?
4. Where would you look for the gene ID?
5. What is the delimiter?

The Five Commands

Today we will build a mental model around five commands:

A useful workflow is:

find
  ↓
grep
  ↓
cut
  ↓
sort
  ↓
uniq

Introducing a Lifesaver: tldr

Navigating the Linux filesystem was simple enough to memorize. But what happens when you start using more and more tools?

Installing tldr and tealdeer

sudo apt install tldr

sudo apt install tealdeer

tldr --update

Searching the filesystem with find

find path/to/directory -name '*.ext'

Search the current directory:

find .

Find all GTF files:

find . -name "*.gtf"

Find regular files matching a pattern.

find . -type f -name "*.fastq.gz"

Do Not Confuse find and grep!

find path/to/directory -name '*.ext'

find searches the filesystem, not the contents of files.

This distinction is important because many people, including myself, frequently confuse find and grep when they first start out with Bash.

find — Useful Options

Search recursively from PATH:

find PATH

Search only regular files:

find . -type f

Search only directories.

find . -type d

Search for name pattern:

find . -name "*.gtf"

Search using case-insensitive name matching:

find . -iname "*.GTF"

A realistic bioinformatics search:

find . -type f -name "*.fastq.gz"

Searching File Contents With grep

grep answers:

Which lines contain this pattern?

Example:

grep "TP53" genes.gtf

Only lines containing TP53 are printed.

grep "gene" genes.gtf

Find lines containing the text gene.

grep -n "TP53" genes.gtf

Also show line numbers.

grep as a Filter

grep "TP53" genes.gtf

chr17 ... BRCA1 ...
chr17 ... TP53 ...
chr17 ... EGFR ...
chr17 ... TP53 ...

chr17 ... TP53 ...
chr17 ... TP53 ...

grep – Useful Options

Show line numbers:

grep -n "TP53" genes.gtf

Use case-insensitive matching:

grep -i "tp53" genes.gtf

Invert the match:

grep -v "exon" genes.gtf

Extended regular expressions (ERE):

grep -E "TP53|BRCA1" genes.gtf

Count matching lines:

grep -c "TP53" genes.gtf

Regex: A Force Multiplier

The real power of CLI data processing tools comes from regular expressions, or regex.

A regex is a way of describing a pattern of text.

Regex Vs. Regular Strings

The second describes a much larger set of possible strings:

chr1
chr2
chr17
chr22

but not:

chrX
chromosome17
17
chr17_extra

. – Match any character

Example:

a.c matches

abc
a1c
a-c

because . means:

any single character

^ — Beginning of a line

^chr

means:

chr must appear at the beginning of the line.

For example:

chr17
chrX
chrM

match.

But my_chr17 does not.

This is useful for identifying structured lines such as FASTA headers or chromosome names.

$ — End of a line

_old$

means:

_old must appear at the end of the line.

So:

gene_old
sample_old

match, but:

gene_old_version2

does not.

Together, ^ and $ are called anchors.

They describe a position in the line rather than a character.

[] – Character Classes

For example:

chr[XY]

matches

chrX
chrY

but chr17.

Using Ranges in Character Classes

Ranges can also be used: [0-9] means

any single digit from 0 to 9

[A-Z] means

any uppercase letter

Quantifiers — Matching Repeated Characters

Regex can also describe how many times something should appear.

| — OR

The | operator means:

Match either the pattern on the left or the pattern on the right.

For example:

chrX|chrY

matches:

chrX
chrY

but not:

chr17
chrM

It is useful when a pattern can have multiple alternatives.

Combining Regex Components

Now we can build more useful patterns.

For example:

^chr[0-9]+$

means:

Extracting Fields With cut

Suppose a TSV file contains:

gene_id gene_name   length
ENSG001 TP53    19149
ENSG002 BRCA1   134289
ENSG003 EGFR    121110

To extract the first field:

cut -f1 genes.tsv

Output:

gene_id
ENSG001
ENSG002
ENSG003

Extract field 2:

cut -f2 genes.tsv

cut and Delimiters

For tab-separated data:

cut -f1 file.tsv

For comma-separated data:

cut -d',' -f2 file.csv

For colon-separated data:

cut -d':' -f1 file.txt

General form:

cut -d'DELIMITER' -fFIELD_NUMBER file

For tab-separated files, -f implies tab in common usage.

cut Can Select Multiple Fields

cut -f1,3 genes.tsv

Fields 1 and 3.

cut -f1-3 genes.tsv

Fields 1 through 3.

cut -f2- genes.tsv

Field 2 through the end.

The key question is always:

What is the delimiter?

Ordering Text With sort

Order in alphabetical order:

sort genes.txt

Order in reverse order:

sort -r genes.txt

sort -n lengths.txt

Numeric Sorting With sort

sort -n lengths.txt

Example:

10
2
30
4

Alphabetical sort can produce:

10
2
30
4

Numerical sort:

2
4
10
30

Sorting by a Field

For TSV-like data:

sort -t$'\t' -k3,3 genes.tsv

Breakdown:

-t$'\t'   delimiter = TAB
-k3,3      sort using field 3

Numerical field sorting:

sort -t$'\t' -k3,3n genes.tsv

Reverse numerical sorting:

sort -t$'\t' -k3,3nr genes.tsv

Removing Adjacent Duplicates With uniq

Given:

BRCA1
BRCA1
EGFR
EGFR
TP53
TP53
TP53

uniq genes.txt

Output:

BRCA1
EGFR
TP53

But there is an important catch:

uniq only removes adjacent duplicates.

What if Duplicates Aren’t Adjacent?

Input:

TP53
BRCA1
TP53
EGFR
BRCA1

Running:

uniq genes.txt

does not remove all duplicate values.

Instead:

TP53
BRCA1
TP53
EGFR
BRCA1

Nothing is adjacent.

uniq Often Needs sort

Input:

TP53
BRCA1
TP53
EGFR
BRCA1

Solution? Use sort first:

sort genes.txt | uniq

first groups duplicates, then removes them.

Counting Duplicates

uniq -c counts adjacent identical lines:

sort genes.txt | uniq -c

Example:

      2 BRCA1
      2 EGFR
      3 TP53

This is already useful bioinformatics information.

For example:

Which gene IDs occur most often in this annotation subset?

The Pipe |

The pipe sends the output of one command to the input of another.

command_1 | command_2

Example:

grep "TP53" genes.gtf | cut -f1

First:

grep → keep matching records

Then:

cut → extract field 1

Reading a Pipeline Left → Right

Consider:

grep "gene" genes.gtf | cut -f1 | sort | uniq

Read it as:

Keep gene records

↓

Extract chromosome

↓

Sort chromosomes

↓

Remove duplicates

A More Useful Pipeline

cut -f3 genes.gtf | sort | uniq -c

    120 exon
     25 gene
     40 transcript

Pipeline Debugging

When a long pipeline fails, split it.

Instead of:

grep "gene_id" genes.gtf | cut -f9 | sort | uniq -c

test:

grep "gene_id" genes.gtf

Then:

grep "gene_id" genes.gtf | cut -f9

Then:

grep "gene_id" genes.gtf | cut -f9 | sort

Then:

grep "gene_id" genes.gtf | cut -f9 | sort | uniq -c

Build complicated commands from verified small steps.

Redirecting Output with >

> operator means Write to a file.

grep "TP53" genes.gtf > tp53_records.txt

Meaning:

Run command

↓

Send stdout to file

Appending to a File With >>

Suppose selected_genes.txt already contains:

TP53

Run:

grep "BRCA1" genes.gtf >> selected_genes.txt

Now the file contains:

TP53
BRCA1

>> appends (= adds to the end) instead of replacing the existing content.

Why sed?

grep is mainly for selecting lines.

cut is mainly for selecting fields.

Simple sed Example

Replace text:

sed 's/TP53/TUMOR_SUPPRESSOR/g' file.txt

General idea:

s/old/new/g

sed Without g

Without g:

sed 's/TP53/TUMOR_SUPPRESSOR/' genes.txt

only the first occurrence on each line is replaced.

sed reads the input, applies a transformation, and prints the result.

By default, it does not modify the original file.

sed Example: Remove a Chromosome Prefix

chr17
chr18
chrX

17
18
X

sed 's/^chr//'

So:

chr17

becomes:

17

The command does not remove every occurrence of chr.

It removes only a chr prefix.

sed Example: Remove Unwanted Lines

GTF files often contain comment lines beginning with #:

##description: evidence-based annotation
# source: HAVANA
chr17 HAVANA gene ...
chr17 HAVANA transcript ...

To remove all comment lines:

sed '/^#/d' annotations.gtf

The important pieces are:

• ^#: line starts with #
• d: delete the matching line

So the command means:

Delete every line that begins with #.

sed and Regex Capture Groups

Regex can also divide a match into groups.

For example:

chr17
chr18
chr22

Suppose we want to transform them into:

chromosome_17
chromosome_18
chromosome_22

sed -E 's/^chr([0-9]+)$/chromosome_\1/'

Inserting Captured Text

sed -E 's/^chr([0-9]+)$/chromosome_\1/'

After capture, the remembered text can then be inserted into the replacement: \1

So:

chr17

is interpreted as:

chr + 17

and becomes:

chromosome_17

This is one of the most useful regex techniques for reformatting structured biological data.

More sed Examples

Remove a prefix:

sed 's/^chr//'

Convert:

chr17
chr18
chrX

into:

17
18
X

Remove a suffix:

sed 's/_old$//'

What is awk for?

awk is useful for structured text.

It can:

• Access multiple fields
• Apply conditions
• Perform calculations
• Create new output

Think of it as:

filter + calculate + format

Basic awk

Consider a GTF record:

chr17   HAVANA   gene   7661779   7687550

awk automatically treats the columns as fields:

$1   chr17
$2   HAVANA
$3   gene
$4   7661779
$5   7687550

awk '{print $1, $3}' genes.gtf

Basic awk (Cont.)

awk '{print $1, $3}' genes.gtf

Example:

awk '{print $1, $3}' genes.gtf

Print field 1 and field 3. Output:

chr17 gene

awk + Conditions

Only print gene records:

awk '$3 == "gene" {print $1, $4, $5}' genes.gtf

Interpretation:

IF field 3 equals "gene"
THEN print fields 1, 4, and 5

This is much more expressive than a simple field extraction.

awk + Calculations

Suppose:

chr17   gene   7661779   7687550

The fields are:

$1      $2     $3        $4
chr17   gene   start     end

For inclusive genomic coordinates:

length = end - start + 1

With awk:

awk '{print $1, $3, $4, $4-$3+1}' file.tsv

Output:

chr17 7661779 7687550 25772

Variables: Remembering Values

awk can store values in variables.

For example, suppose we want the total length of all regions:

awk '{total += $4-$3+1} END {print total}' file.tsv

Here:

total += $4-$3+1

adds the length of each record to total.

Useful Patterns for awk

This gives awk a very useful pattern:

awk '{total += $5-$4+1} END {print total}' file.tsv

• FOR every record:
  • calculate its length
  • add it to total
• AFTER all records:
  • print total

This is something tools such as cut cannot do by themselves.

Formatting the Output

awk can also combine fields into a new format.

Suppose:

chr17   gene   7661779   7687550
chr18   gene   46841871  46854673

We can select specific fields to produce a BED-like representation:

awk '{print $1, $3, $4}' file.tsv

Output:

chr17 7661779 7687550
chr18 46841871 46854673

awk '{print "Region:", $1 ":" $3 "-" $4}' file.tsv

Region: chr17:7661779-7687550
Region: chr18:46841871-46854673

Regex in awk

Just like grep and sed, awk can use regex.

For example:

awk '$1 ~ /^chr[0-9]+$/ {print $1}' file.tsv

The operator:

~

means:

matches this regular expression

So:

$1 ~ /^chr[0-9]+$/

Regex in awk (Cont.)

$1 ~ /^chr[0-9]+$/

This could match:

chr1
chr17
chr22

but not:

chrX
chromosome17

The same regex concepts used with grep and sed therefore carry over to awk.

awk as a Mini Data-Processing Language

awk '$3 == "gene" {print $1, $4, $5, $5-$4+1}' annotations.gtf

Tool Selection Mental Model

Use:

• find – locate files
• grep – filter lines
• cut – select fields
• sort – order lines
• uniq – collapse adjacent duplicates

Then:

• sed – transform text
• awk – process structured text with logic

Choose the simplest tool that matches the problem.

What Is a Variable?

A variable stores a value.

project_dir="$HOME/bioinformatics"

Then:

echo "$project_dir"

Output:

/home/student/bioinformatics

Use the variable elsewhere:

find "$project_dir" -type f -name "*.gtf"

Now one value controls multiple commands.

Parameterizing Workflows With Variables

Imagine a workflow using:

project directory
reference genome
annotation file
sample name
output directory

Hard-coding every path makes scripts difficult to reuse.

Instead:

project_dir="/data/project"
sample_id="sample_A"
output_dir="$project_dir/results"

Then:

grep "$sample_id" "$project_dir/metadata.tsv"

This makes the workflow parameterized (= tham số hóa).

Persistent Variables with ~/.bashrc

~/.bashrc is commonly used for interactive shell configuration.

Example:

export BIOINFO_HOME="$HOME/bioinformatics"

After reloading:

source ~/.bashrc

The variable will be available in future shell sessions.

Manual Workflow vs. Script

Suppose every sample requires:

grep "gene" annotation.gtf > genes.txt
cut -f1 genes.txt | sort | uniq > chromosomes.txt

Doing this manually for:

sample_A
sample_B
sample_C
sample_D
...
sample_Z

is slow and error-prone.

A script turns the procedure into a reusable program.

A Minimal Bash Script

Create:

nano summarize_genes.sh

Contents:

#!/usr/bin/env bash

grep "gene" annotation.gtf \
  | cut -f1 \
  | sort \
  | uniq \
  > chromosomes.txt

The first line is the shebang.

It tells the operating system how to interpret the script.

Bash Strict Mode

A common safety-oriented starting point is:

#!/usr/bin/env bash
set -euo pipefail

Meaning:

-e  stop when a command fails
-u  treat unset variables as errors
-o pipefail  detect failures inside pipelines

This doesn’t magically make a script safe.

It makes some classes of mistakes easier to detect.

Variables Inside Scripts

#!/usr/bin/env bash
set -euo pipefail

input_file="annotation.gtf"
output_file="chromosomes.txt"

grep "gene" "$input_file" \
  | cut -f1 \
  | sort \
  | uniq \
  > "$output_file"

Now the filenames are easy to change.

Making the Script Executable

Use:

chmod +x summarize_genes.sh

Run it with:

./summarize_genes.sh

Alternative:

bash summarize_genes.sh

Turn a Repeated Task Into a Parameterized Script

Suppose we have:

annotation/
├── sample_A.gtf
├── sample_B.gtf
└── sample_C.gtf

A script can take the input filename as an argument.

#!/usr/bin/env bash
set -euo pipefail

input_file="$1"
output_file="$2"

grep "gene" "$input_file" \
  | cut -f1 \
  | sort \
  | uniq \
  > "$output_file"

Run:

./summarize_genes.sh annotation/sample_A.gtf results/sample_A_chromosomes.txt

What We Have Built

Back to Our Opening Problem

At the beginning of the lecture, we were given:

rnaseq_project/

├── annotation/
│   ├── gencode_subset.gtf
│   ├── gencode_subset.gtf.backup
│   └── notes.txt
│
├── metadata/
│   └── samples.tsv
│
└── results/

Our collaborator asked for a quick QC report.

Let’s solve the problem one question at a time.

Step 1 — Which Annotation File Do We Have?

Start by finding GTF files:

find annotation -type f -name "*.gtf"

We can now store the path in a variable:

annotation_file="annotation/gencode_subset.gtf"

Step 2 — Keep Only Gene Records

A GTF contains many types of records:

gene
transcript
exon
CDS
...

Our audit is about gene records.

The feature type is field 3:

awk -F '\t' '$3 == "gene" {print}' "$annotation_file"

Keep the record when field 3 is gene.

Step 3 — Keep Protein-Coding Genes

GTF stores the gene biotype in the attributes:

gene_type "protein_coding"

We can add this condition:

awk -F '\t' \
  '$3 == "gene" && $9 ~ /gene_type "protein_coding"/' \
  "$annotation_file"

Now we’re asking:

Is this a gene record? AND Is it protein-coding?

Step 4 — Restrict the Audit to Autosomes

Our audit defines autosomes as chromosomes:

chr1–chr22

We can describe that using regex:

^chr([1-9]|1[0-9]|2[0-2])$

Then:

awk -F '\t' \
  '$3 == "gene" &&
  $9 ~ /gene_type "protein_coding"/ &&
  $1 ~ /^chr([1-9]|1[0-9]|2[0-2])$/' \
  "$annotation_file"

Now the filter means:

• field 3 = gene AND

• gene_type = protein_coding AND

• field 1 = chr1–chr22

Notice how the regex is being used as part of the logic.

Step 5 — Find Duplicated Gene IDs

The gene ID is stored inside field 9:

gene_id "ENSG00000123456";

First, keep only the attribute field:

... | cut -f9

Then extract the gene ID:

... | sed -n 's/.*gene_id "\([^"]*\)".*/\1/p'

Now we have a simple list:

ENSG00000123456
ENSG00000198765
ENSG00000123456
ENSG00000234567

... | sort

Step 5 — Find Duplicated Gene IDs (Cont.)

... | sort

Then count identical adjacent values:

... | sort | uniq -c

For example:

    1 ENSG00000198765
    2 ENSG00000123456
    1 ENSG00000234567

The duplicated gene is therefore the ID with a count greater than 1.

We can filter for those:

... | sort | uniq -c | awk '$1 > 1'

This gives us:

Which gene IDs occur more than once?

Putting the Duplicate Check Together

The whole pipeline is now:

awk -F '\t' \
  '$3 == "gene" &&
   $9 ~ /gene_type "protein_coding"/ &&
   $1 ~ /^chr([1-9]|1[0-9]|2[0-2])$/' \
  "$annotation_file" \
  | cut -f9 \
  | sed -n 's/.*gene_id "\([^"]*\)".*/\1/p' \
  | sort \
  | uniq -c \
  | awk '$1 > 1'

Step 6 — Extract Coordinates and Calculate Genomic Span

Now we need:

• chromosome
• start
• end
• span

For a GTF record:

$1 = chromosome
$4 = start
$5 = end

For inclusive GTF coordinates:

span = end - start + 1

Calculating Genomic Span (Cont.)

awk -F '\t' '$3 == "gene" &&
$9 ~ /gene_type "protein_coding"/ &&
$1 ~ /^chr([1-9]|1[0-9]|2[0-2])$/ {
  print $1, $4, $5, $5-$4+1
}
' "$annotation_file"

The last column is calculated rather than extracted.

That is where awk becomes especially useful:

read fields → apply logic → calculate → print

Nice. Where would we go from here?

With these tools, you can perform further analysis by yourself, then redirect the output directly to a file using output redirection.

By now, you probably noticed that some commands, especially the awk pipeline, were getting harder to read. This is often the point where a script or a higher-level language such as Python becomes useful.

Example: When You Should Consider Using Another Tool

Instead of printing the results to the terminal, we want to create a TSV report.

report_file="results/annotation_audit.tsv"

Then:

awk -F '\t' '
BEGIN {
    OFS = "\t"
    print "gene_id", "gene_name", "chromosome", "start", "end", "span_bp"
}
$3 == "gene" && $9 ~ /gene_type "protein_coding"/ && $1 ~ /^chr([1-9]|1[0-9]|2[0-2])$/ {
    gene_id = $9
    sub(/.*gene_id "/, "", gene_id)
    sub(/".*/, "", gene_id)
    gene_name = $9
    sub(/.*gene_name "/, "", gene_name)
    sub(/".*/, "", gene_name)
    span = $5 - $4 + 1
    print gene_id, gene_name, $1, $4, $5, span
}
' "$annotation_file" > "$report_file"

The results are good, but reading and maintaining this kind of code would be a pain.