Day 15: Whole Genome Sequencing (WGS) Principles & Microbial Genome Analysis

1. Whole Genome Sequencing (WGS) Paradigms

Whole Genome Sequencing (WGS) determines the complete nucleotide sequence of an organism's chromosomal and extrachromosomal (plasmid, mitochondrial) DNA. In molecular microbiology and biotechnology, microbial WGS underpins pathogen surveillance, antimicrobial resistance (AMR) tracking, synthetic biology strain verification, and industrial enzyme discovery.

2. De Novo Assembly Mathematics: De Bruijn Graphs vs. Overlap-Layout-Consensus

Because short-read sequencing cannot reconstruct genomes by simple visual alignment, computational assemblers decompose reads into overlapping mathematical graphs:

De Bruijn Graph (DBG) Assembly

• Used by short-read assemblers (SPAdes, MEGAHIT, Velvet).
• Decomposes each read of length $L$ into overlapping substrings of fixed length $k$ (termed $k$-mers).
• Nodes in the graph represent $(k-1)$-mers; directed edges represent observed $k$-mers in the reads.
• Assembling the genome corresponds to finding an Eulerian path (visiting every edge exactly once) through the graph.

Assembly Quality Assessment: $N_{50}$ and $L_{50}$

A de novo assembly yields a set of contiguous genomic fragments termed contigs. Assembly contiguity is quantified using $N_{50}$:

$$\sum_{i=1}^{M} \text{Length}(\text{Contig}_i) \ge 0.50 \times \text{Total Assembly Length}$$

3. Microbial Genome Annotation Pipeline

Once contigs are assembled, structural and functional annotation classifies every base pair:

1. Ribosomal RNA (rRNA) Detection: Identified via covariance models using barrnap (identifies 5S, 16S, and 23S rRNAs).
2. Transfer RNA (tRNA) Identification: Predicted via tRNAscan-SE based on secondary cloverleaf structure models.
3. Protein-Coding Genes (CDS): Predicted using ab initio prokaryotic gene finders (Prodigal) based on hexamer frequencies and Shine-Dalgarno motifs.
4. Functional Assignment: Assigned by aligning predicted CDS against curated databases (UniProt, COG, Pfam) using Prokka or Bakta.

4. Production Python Script: Contig $N_{50}$, $L_{50}$ & Assembly QC Calculator

#!/usr/bin/env python3
"""
Day 15 Practical: De Novo Genome Assembly Quality Profiler
Computes N50, L50, GC content, and assembly size from multi-FASTA contigs.
"""

import sys
from Bio import SeqIO

def calculate_assembly_metrics(contig_fasta: str):
    """Calculate total assembly length, contig count, N50, and L50."""
    lengths = []
    total_bases = 0
    total_gc = 0
    
    for record in SeqIO.parse(contig_fasta, "fasta"):
        seq_str = str(record.seq).upper()
        l = len(seq_str)
        lengths.append(l)
        total_bases += l
        total_gc += seq_str.count("G") + seq_str.count("C")
        
    if not lengths:
        print("No sequences found in FASTA.", file=sys.stderr)
        return
        
    # Sort contigs descending by length
    lengths.sort(reverse=True)
    num_contigs = len(lengths)
    gc_percent = (total_gc / total_bases) * 100
    
    # Calculate N50 and L50
    half_length = total_bases / 2.0
    cumulative_sum = 0
    n50 = 0
    l50 = 0
    
    for rank, length in enumerate(lengths, 1):
        cumulative_sum += length
        if cumulative_sum >= half_length:
            n50 = length
            l50 = rank
            break
            
    print("=== WHOLE GENOME ASSEMBLY QUALITY REPORT ===")
    print(f"Total Assembly Size: {total_bases:,} bp ({round(total_bases / 1e6, 2)} Mb)")
    print(f"Total Number of Contigs: {num_contigs}")
    print(f"Largest Contig: {lengths[0]:,} bp")
    print(f"Shortest Contig: {lengths[-1]:,} bp")
    print(f"Genome GC Content: {round(gc_percent, 2)}%")
    print(f"Assembly N50: {n50:,} bp")
    print(f"Assembly L50: {l50} contig(s)")
    
    return {
        "size": total_bases,
        "contigs": num_contigs,
        "n50": n50,
        "l50": l50,
        "gc": round(gc_percent, 2)
    }

if __name__ == "__main__":
    # Test on synthetic contig distribution
    sample_assembly = """\
>contig_1 [length=250000]
""" + "GATC" * 62500 + """
>contig_2 [length=150000]
""" + "ATCG" * 37500 + """
>contig_3 [length=80000]
""" + "CGAT" * 20000 + """
>contig_4 [length=20000]
""" + "TAGC" * 5000
    
    with open("sample_microbial_contigs.fasta", "w") as f:
        f.write(sample_assembly)
        
    calculate_assembly_metrics("sample_microbial_contigs.fasta")