Comprehensive Notes on Microbial Genome Assembly, Quality Assessment, Annotation, Resequencing, and Pangenomics

Principles of Genome Assembly, Scaffolding, and Contig Ordering

• Paired-End Sequencing and Gap Representation:

  • Paired-end sequencing produces a forward read and a reverse read from the opposing ends of a single DNA fragment.

  • When forward and reverse reads map to distinct contigs, their pairing demonstrates physical linkage and defines the relative orientation and proximity of the contigs.

  • If the middle span of a sequenced fragment fails to sequence through due to high sequence complexity, low quality, or technical limitations, the intervening sequence remains unresolved.

  • Unresolved spanning regions within a scaffold are represented by ambiguous base characters, designated as $N$s (e.g., $NNN$).

• Contigs Versus Scaffolds:

  • Contig: A continuous, unbroken consensus sequence assembled from overlapping reads without intervening gaps.

  • Scaffold: A series of oriented and ordered contigs connected across unsequenced gaps of estimated lengths represented by runs of $N$ bases.

  • Impact on Annotation: Contigs and scaffolds yield identical gene annotation outputs because annotation pipelines only process discrete, well-defined nucleotide regions and skip non-informative ambiguous ($N$) stretches.

  • Impact on Phylogenetics: Scaffolds provide positional and syntenic context. Even when problematic regions persist as ambiguous stretches across related organisms, knowing relative contig orientation enables more reliable comparative alignments and superior phylogenetic statistics.

• Genome Completeness Criteria:

  • Complete (Closed) Genome: Fully resolved, circularized, or end-to-end chromosomal sequence with zero undefined gaps, requiring no scaffolding prediction ($0$ ambiguous $N$ residues).

  • Incomplete (Draft) Genome: An assembly composed of multiple contigs placed into predictive scaffolds, containing unresolved crossover, repeat, or low-coverage gaps that cannot be resolved using short-read data alone.

  • Resolution with Long-Read Sequencing: Resolving repeat and crossover structures requires adding long-read platforms (e.g., Oxford Nanopore Technologies or Pacific Biosciences [PacBio]) to bridge repetitive genomic regions.

Assembly Quality Evaluation and Metrics

• De Novo Assembler Benchmarking:

  • Typical read inputs must first be preprocessed using cutadapt to discard low-quality bases and trim sequencing adapters.

  • SPAdes: Consistently yields the highest-quality assemblies with the highest $N_{50}$ values, making it the preferred standard assembler.

  • Velvet and Sischia: Often generate a higher overall number of contigs, but produce lower $N_{50}$ values, indicating greater fragmentation.

• Evaluation via QUAST (Quality Assessment Tool for Genome Assemblies):

  • Assembler output logs and raw files (e.g., assembly graphs, contig paths) are difficult to parse directly; QUAST standardizes the evaluation.

  • QUAST ingests contig or scaffold FASTA files to compute assembly length, contig counts above size thresholds (e.g., contigs $> 5000\,\text{bp}$), and $N_{50}$ metrics.

  • QUAST supports simultaneous multi-sample benchmarking by presenting comparative metrics in side-by-side tabular reports.

• Definition and Mathematical Derivation of $N_{50}$:

  • The $N_{50}$ metric is a standard statistical measure of assembly contiguity, reflecting sequence continuity relative to total assembled mass and indicating sequencing quality and coverage depth.

  • In a fully closed genome, the $N_{50}$ is exactly equal to the total length of the genome.

  • A high $N_{50}$ indicates an assembly dominated by long contigs approaching completion, whereas a low $N_{50}$ indicates high fragmentation into many small pieces.

  • Calculation Algorithm:

    1. Sort all generated contigs in descending order of length (from largest to smallest).

    2. Sum the lengths of all contigs to determine the total assembly size in kilobase pairs ($\text{kb}$).

    3. Calculate the $50\%$ threshold of the total assembly size: $\text{Threshold} = \frac{\text{Total Genome Length}}{2}$.

    4. Sum the lengths of the sorted contigs cumulatively, beginning with the largest.

    5. Identify the specific contig whose addition causes the cumulative sum to equal or exceed the $50\%$ threshold.

    6. The length of that individual contig is defined as the $N_{50}$.

  • Step-by-Step Calculation Example:

    • Total assembled genome size: $600\,\text{kb}$.

    • $50\%$ threshold: $300\,\text{kb}$.

    • Contig 1: $150\,\text{kb}$ (Cumulative: $150\,\text{kb}$; $< 300\,\text{kb}$).

    • Contig 2: $80\,\text{kb}$ (Cumulative: $230\,\text{kb}$; $< 300\,\text{kb}$).

    • Contig 3: $50\,\text{kb}$ (Cumulative: $280\,\text{kb}$; $< 300\,\text{kb}$).

    • Contig 4: $35\,\text{kb}$ (Cumulative: $315\,\text{kb}$; $\ge 300\,\text{kb}$).

    • Result: Contig 4 crosses the threshold, so $N_{50} = 35\,\text{kb}$.

Galaxy Assembly and Hybrid Workflow Execution

• Preprocessing with Cutadapt:

  • Input: Paired-end fastq collection (e.g., SRA accession DRX4405...).

  • Quality parameter thresholds: 30 30 22 22 to ensure quality trimming on short-read datasets.

  • Output: Filtered paired-end read collections and statistical trimming reports.

• Assembly Execution in SPAdes:

  • Short-read paired collections are loaded into the standard paired library input.

  • For Hybrid Assemblies, long reads are incorporated via the "additional read files" interface, designating whether the platform is Oxford Nanopore or PacBio.

  • SPAdes generates multiple output entities: corrected reads, assembly graph, contigs, contig paths, and scaffolds.

• Assembly Benchmarking Results:

  • Contig count for contigs $> 0\,\text{bp}$: Yields $455$ contigs in standard runs ($71\%$ of runs); select runs generate $415$ due to scaffold-versus-contig boundary differences.

  • Total assembled genome size: Approximately $697{,}452\,\text{bp}$ ($\approx 697\,\text{kb}$, reported by $67\%$ of standardized runs).

  • Minor base-pair fluctuations occur across runs due to stochastic sequence read trimming at terminal contig ends.

Prokaryotic Genome Annotation with Prokka

• Input Specifications and Reference Identification:

  • Annotation tools accept assembled FASTA files containing contigs or scaffolds.

  • Providing specific genus and species designations refines functional gene descriptions (e.g., designating organism-specific oxidoreductases) without altering total gene counts.

  • For unclassified organisms, run contigs or scaffolds through RefSeq Masher Contains to match against RefSeq reference genomes.

• Formula for Calculating Total Gene Content:

  • The total number of functional genes in a prokaryotic annotation is calculated by summing all discrete RNA and protein-coding features:         $\text{Total Genes} = \text{CDS} + \text{rRNA} + \text{tRNA} + \text{tmRNA}$

  • CDS: Coding sequence (protein-encoding open reading frames).

  • rRNA: Ribosomal RNA features.

  • tRNA: Transfer RNA features.

  • tmRNA: Transfer-messenger RNA features.

  • Non-Gene Features: Repeat regions are structurally annotated elements but are never enumerated as functional genes.

• Output Parsing (GFF and TSV):

  • Prokka generates GFF (General Feature Format) and TSV (Tab-Separated Values) summary files.

  • Parsing the GFF requires differentiating distinct functional annotations (gene, CDS) from the same physical locus to avoid double-counting overlapping features.

Sequence Resequencing and Variant Calling Pipelines

• Clinical Resequencing and Pathogen Microevolution:

  • Analysis of longitudinal bacterial isolates (e.g., Mycobacterium abscessus isolated from cystic fibrosis patients via the Emory Cystic Fibrosis Center) reveals in-host evolutionary dynamics.

  • Comparing an ancestral isolate against subsequent persistent isolates identifies acquired single nucleotide polymorphisms (SNPs), insertions, and deletions.

  • Functional mutations correlate directly with clinical declines, hyper-virulence, and the emergence of antimicrobial resistance (AMR), analogous to the surveillance of single mutations in SARS-CoV-2 (COVID-19) variants driving increased transmissibility and mortality.

• Reference Genome Acquisition and Mapping:

  • Reference genomes must be retrieved as FASTA genome sequence files using tools such as NCBI Datasets Genome rather than GFF annotation files.

  • Reads are aligned to reference sequences using Bowtie/Bowtie2, generating alignment statistics that specify the exact percentage of mapped reads.

• Variant Calling with Snippy Versus Standard VCF Workflows:

  • Snippy Workflow: Integrates mapping and variant calling into a single automated execution, bypassing manual intermediate Bowtie runs.

  • Reference Formatting Requirement: Snippy requires uncompressed FASTA reference files; passing compressed archives (.gz) triggers pipeline errors.

  • Variant Count Metrics: In benchmark testing, the dataset contains $12$ verified SNPs, within a broader set of $18$ total sequence variations (incorporating indels).

  • Snippy Limitations vs. VCFtools: Snippy relies on hardcoded internal filtering cutoffs that can suppress legitimate low-frequency variants. Formal research publications favor comprehensive variant calling workflows evaluated with VCFtools due to transparent and flexible statistical threshold filtering.

Pangenomics and Microbial Genome Plasticity

• Core Concepts in Pangenomics:

  • Pangenomic analysis involves the alignment and comparison of all genes across multiple isolate genomes to characterize gene presence and absence patterns.

  • Annotated Prokka output files (GFF format) serve as inputs to construct orthologous gene clusters.

• Genome-Wide Association Studies (GWAS) and Trait Linkage:

  • Pangenome presence/absence matrices enable GWAS to associate specific accessory genes with clinical or phenotypic traits.

  • Pathogenicity in Escherichia coli: Commensal E. coli strains lack toxin operons and cause mild or negligible illness. In contrast, enterohemorrhagic E. coli isolates harbor the Shiga toxin operon, which causes severe hemorrhagic colitis (bloody diarrhea) and systemic complications.

  • Therapeutic Resistance in Cystic Fibrosis Isolates: Comparative pangenomics on Emory Cystic Fibrosis Center patient strains revealed an operon that metabolizes or confers resistance to nitric oxide ($\text{NO}$), an inhaled antimicrobial treatment, thereby directly causing therapeutic failure.

• Open Versus Closed Genomes:

  • Open Genome:

    • Organisms display substantial genomic plasticity, continually acquiring novel genes from the environment via horizontal gene transfer (HGT).

    • Characterized by a large, expanding accessory genome and a relatively small core genome fraction.

    • Typical of organisms that colonize multiple variable ecological niches (e.g., Pseudomonas, which persists in soil, freshwater systems, and hospital settings).

  • Closed Genome:

    • Organisms display minimal genomic variation between strains and rarely incorporate exogenous genetic elements from the environment.

    • Characterized by a conserved gene repertoire across the entire species, with almost all genes belonging to the core genome.

    • Typical of specialized, host-restricted pathogens adapted to a single niche (e.g., Mycobacterium tuberculosis, which does not survive as a free-living environmental saprophyte).

Questions & Discussion

• Canadian Galaxy Server Discrepancies:

  • Question/Concern: Can the Canadian Galaxy server mirror be used for homework and examination workflows?

  • Response: Third-party or regional server mirrors use divergent tool versions, dependencies, and analysis pipelines, which introduce discrepancies in assembly and variant calling metrics. Standardized instances must be used for course assessments, though independent platforms are acceptable for external research.

• Import Tool Identifier Failures (DRX4405...):

  • Question/Concern: Tool runs fail when attempting to import dataset accession DRX4405... using default tool names.

  • Response: The tool version failed to parse that accession string format. The solution is to use the preceding software version of the import tool to correctly pull down the raw sequence data.

• Assembly Discrepancies in Contig Count and Genome Size (Dylan's Query on Lab 6, Question 9):

  • Question/Concern: Why do distinct runs yield contig counts differing between $415$ and $455$, and slight length shifts around $697{,}452\,\text{bp}$?

  • Response: These minor variances stem from differences between contig versus scaffold boundaries and minor stochastic clipping of read ends during preprocessing. Any result falling within the approximately $697\,\text{kb}$ range is scientifically accurate.

• Data Download Performance: Faster Download Versus NCBI Datasets:

  • Question/Concern: Does using faster-download utilities introduce variations compared to standard NCBI dataset extraction tools?

  • Response: The retrieved data is fundamentally identical across downloaders. Standard NCBI downloads are more reliable, as tools marketed as "faster" frequently experience stability issues.

• Snippy Variant Output Discrepancy (Will's Observation):

  • Question/Concern: Why are some analyses showing $8$, others $12$, and some $18$ variants in Snippy?

  • Response: The value $18$ represents the total number of all sequence variations combined (including indels and complex mutations). The value $12$ represents the true, isolated single nucleotide polymorphisms (SNPs). The value $8$ was an unupdated key from an older software release.