7: Columnar Data Wrangling & Genomic Arithmetic with awk

Chapter 7: Columnar Data Wrangling & Genomic Arithmetic with awk

Almost all biological annotation and variation files—including BED, VCF, GFF3, GTF, and SAM—are structured as tab-separated columnar tables. awk is a complete Turing-complete programming language designed specifically for columnar data processing and report generation. With awk, you can perform coordinate arithmetic, calculate mean read depths, and filter multi-million-row genomic callsets in a single command line.

1. The awk Program Structure

An awk script consists of three coordinated execution blocks:

awk 'BEGIN { print "Header" }  { pattern / action }  END { print "Summary" }' file.tsv

• BEGIN { ... } : Executes once before any input lines are read. Used to set field separators and print headers.
• { ... } : Executes once for every input line matching the optional condition.
• END { ... } : Executes once after the entire file has been processed. Used to print totals, means, and summaries.

2. Special Variables & Field Splitting

• $0 : The entire current line.
• $1, $2, $3 ... $NF : Specific fields/columns across the line.
• NF : Number of Fields in the current line (total columns).
• NR : Current Record Number (line counter across input).
• FS : Field Separator (e.g. -F' ' for tab-delimited files).

3. BED Interval Arithmetic (Calculating Region Spans)

A standard 3-column BED file defines chrom ($1), chromStart ($2), and chromEnd ($3). The length of any interval is calculated as $3 - $2:

# Print chromosome, start, end, and calculated interval length
awk -F'	' '{ print $1, $2, $3, $3 - $2 }' intervals.bed

# Filter for large structural intervals greater than 10,000 base pairs
awk -F'	' '($3 - $2) > 10000 { print $0 }' intervals.bed

# Calculate total cumulative genomic footprint (sum of all intervals in bp):
awk -F'	' '{ total += ($3 - $2) } END { print "Total Genomic Span:", total, "bp" }' intervals.bed

4. Filtering VCF Records with Mathematical Precision

VCF columns: $1=CHROM, $2=POS, $3=ID, $4=REF, $5=ALT, $6=QUAL, $7=FILTER:

# Skip comment header lines (^#), ensure QUAL >= 30, and require FILTER == "PASS"
awk -F'	' '!/^#/ && $6 >= 30 && $7 == "PASS" { print $1, $2, $4, $5, $6 }' cohort.vcf

# Count total high-quality SNVs vs Indels:
awk -F'	' '!/^#/ && $7 == "PASS" {
    if (length($4) == 1 && length($5) == 1) snv++;
    else indel++;
} END {
    print "High-Quality SNVs:", snv;
    print "High-Quality Indels:", indel;
}' cohort.vcf

5. Computing Mean Sequencing Depth & Custom Formatting

# Given a 3-column depth file from "samtools depth" ($1=chr, $2=pos, $3=depth):
# Calculate average genome-wide sequencing coverage depth:
awk '{ sum += $3; count++ } END { printf "Mean Coverage Depth: %.2fX across %d bases
", sum/count, count }' coverage.depth