Reference Genome

Comparative genomics starts with knowing your comparison group. The Great Auk (Pinguinus impennis, extinct since the mid-19th century) belonged to family Alcidae — the auks — within the order Charadriiformes. Before picking any genome, you need to know which living species are close enough relatives to be scientifically meaningful comparators, not just "any bird genome." Taxonomic distance affects everything downstream: how useful a species is as a mapping reference, how informative it is phylogenetically, and how interpretable any comparison ends up being.

Pinguinus was a single-species genus — the Great Auk had no other living members of its own genus. Phylogenetic work places Alca torda (the razorbill) as its closest living relative, which is why the razorbill will matter a great deal later in this module.

Both diagrams above are original illustrations redrawn from the published topologies, not reproductions of the original figures (see sources below) — the original figures weren't clearly licensed for reuse here, and a simplified version better suits a two-species zoom-in anyway.

Using the NCBI Taxonomy browser, look up family Alcidae and compile your own list of extant genera and representative species in your notes. You don't need to work out the exact phylogeny — focus on recognising who the family's living members actually are, beyond the razorbill.

Search "Alcidae" directly in the NCBI Taxonomy browser and look at its child taxa — this is faster and more reliable than trying to recall species names from memory.

Why isn't "any bird genome" an adequate comparator here? What specifically would you lose — for mapping, and for a phylogenetic tree — if you picked a distantly related seabird instead of a close Alcidae relative?

Genome assemblies for most species live in NCBI Assembly (searchable via NCBI Datasets), with GenBank and RefSeq as the two underlying databases. An accession starting GCA_ is a GenBank assembly — submitted by any research group, not independently curated by NCBI. One starting GCF_ is a RefSeq assembly — a GenBank assembly that NCBI has additionally processed and annotated. A species can have several assemblies at once, from different groups, of very different quality — "a genome exists for this species" does not mean every assembly of it is equally usable. ENA (European Nucleotide Archive) mirrors much of the same data if you ever need an alternative source, but we'll centre this workflow on NCBI.

Using NCBI Datasets (Genome), search for 2–3 species from your Exercise 1 list, including the razorbill. For each, record in your notes: how many assemblies exist, their accessions, and which (if any) is flagged as RefSeq.

The NCBI Datasets results table has sortable columns for assembly level and RefSeq category — use them to compare candidates at a glance rather than opening each record individually.

For the razorbill, did you find more than one assembly? Independent of quality, what does the GCA_ vs GCF_ prefix tell you about an assembly's status — and what does it not tell you?

Reference selection is a biological and analytical decision, not "pick the assembly with the biggest N50." Metrics worth weighing together include:

• Assembly level — contig, scaffold, or chromosome-level (how far sequence has been ordered and joined)
• Contig N50 / Scaffold N50 — the length at which half the assembly's total size is contained in pieces that size or longer; higher usually means more contiguous, but contig N50 and scaffold N50 measure different things
• Number of contigs/scaffolds — fewer, larger pieces are generally easier to work with
• BUSCO completeness — the proportion of expected single-copy genes found intact, a proxy for how much real biology is actually present and unbroken
• Gaps — regions of unknown sequence (often shown as runs of N) within an otherwise joined assembly
• RefSeq vs GenBank status, and provenance — who produced it, whether it's linked to a peer-reviewed paper, and how it was generated (short-read vs long-read, Hi-C scaffolded, etc.)
• Primary vs alternate/haplotype-resolved — some modern assemblies split a diploid genome into two haplotypes rather than collapsing it into one representative sequence

Two fictional candidate assemblies for the same species, to illustrate that these metrics can pull in different directions:

Now do this for real. Using the candidates you found in Exercise 2 for the razorbill (or another species with more than one assembly), build the same comparison for your actual options using NCBI Datasets' reported metrics, then write a justified decision: which would you choose, and why?

NCBI Datasets' genome record page reports most of these metrics directly (assembly stats, BUSCO where available) — you shouldn't need to calculate any of them yourself.

Which single metric, if you'd looked at only that one, would have given you a misleading answer? And separately: why shouldn't every Alcidae species automatically be mapped to the razorbill genome, even though it's well-studied? What's the general strategy for choosing a reference for an extant species, versus for the Great Auk specifically?

The Great Auk has no living member of its own genus, so there is no conspecific reference genome to map its reads against. Instead, we use the razorbill — a closely related but distinct species. Mapping reads to a different species' genome is called heterospecific mapping, and it has real, predictable consequences:

Ancient DNA compounds this: Great Auk reads will typically be short and chemically damaged (characteristic cytosine-to-thymine misincorporation near fragment ends), which makes divergent regions even harder to map correctly. We won't run any of this yet — the quantitative test comes once you reach the mapping module — but the causal logic above is worth understanding now, since it's the reason reference choice matters at all.

In your own words, walk through what happens to a single Great Auk sequencing read that falls in a genomic region where the Great Auk and razorbill sequences differ substantially — follow it through each step of the chain above.

Think about what an aligner does when a read has too many mismatches to confidently place — does it map with low confidence, map incorrectly, or simply fail to map at all?

Why might heterospecific mapping specifically cause heterozygosity to be underestimated, rather than just noisier in both directions? Why does ancient DNA make this worse rather than just adding random noise?

A filename tells a future reader nothing, and a link copied from a browser can break, redirect, or silently point at a different assembly version later. An accession is a stable, citable identifier that anyone can resolve back to the exact same assembly, indefinitely. NCBI's own datasets command-line tool is built around exactly this: datasets download genome accession <ACCESSION> fetches the official data package for that accession, with no manual URL involved at all.

Mjolnir's software documentation doesn't list datasets as a pre-installed module (it names samtools as its example tool, not this one) — so, exactly as in Module 00, the first step is to check rather than assume.

Check whether datasets is already available on Mjolnir. If it isn't, fetch NCBI's own static binary directly into your personal space (no admin rights needed) and use that instead. Either way, use it to download your chosen razorbill assembly by its accession from Exercise 3.

Two commands from Module 00 check whether a tool is available. If neither finds anything, NCBI publishes the tool itself as a plain downloadable binary — no module system or conda environment required.

Why is downloading by accession more reproducible than right-clicking "download" on a web page? Exactly what would someone else need — written down, nothing more — to redo this identical download in a year's time?

NCBI Datasets packages arrive as a .zip archive with a nested structure — typically ncbi_dataset/data/<ACCESSION>/ containing the genomic FASTA alongside metadata files (an assembly data report, sequence report, and similar). You'll need to unzip it, find the actual FASTA inside, and confirm it before trusting it — the same "inspect before you use it" habit from Module 00, just applied to a real, messier package this time. From here, each species gets its own named subdirectory inside your project's reference collection, so the structure stays legible as it grows:

great_auk/
  reference_genomes/
    Alca_torda/
    Fratercula_arctica/
    Uria_aalge/
  metadata/
  analysis/
  raw_reads/

Species subfolders are named as they're added — you'll only have Alca_torda/ after this exercise; the others appear later as the Alcidae collection grows.

Unzip the downloaded package, locate the genomic FASTA inside it, confirm it with the file tools you already know, then move it into a new Alca_torda/ subdirectory inside reference_genomes/. Also create metadata/, analysis/, and raw_reads/ as siblings of reference_genomes/ — you won't put anything in the last two yet, but the structure should exist.

unzip <file>.zip extracts an archive (new command — your first .gz files in Module 01 didn't need it, since gzip and zip are different formats). find . -name "*.fna" will locate the FASTA regardless of how deeply it's nested.

How many sequences did grep -c '^>' report, and does that roughly match the assembly level you recorded back in Exercise 2/3 — a small number close to a chromosome count, or a much larger number for a scaffold-level assembly? Separately: why keep each species in its own subfolder, rather than putting every FASTA directly inside reference_genomes/?

A folder full of FASTA files means nothing on its own — without a record of what each file is, where it came from, and why it was chosen, neither you nor anyone else can reconstruct this work later. A plain-text, tab-separated manifest is simple, diffable, and easy to extend as you add more species to the Alcidae collection.

Create metadata/genome_manifest.tsv with a header row and one row for the razorbill assembly you just downloaded. Include at minimum: species, common_name, assembly_accession, assembly_name, source, assembly_level, genome_size, contig_N50, scaffold_N50, download_date, intended_use (e.g. "Great Auk mapping reference"), and notes. You'll add a new row for each additional species as the Alcidae collection grows — the file just needs to exist and be correct for now.

Typing literal tab characters into a heredoc is error-prone. printf with explicit \t and \n is more reliable than echo for building a TSV line by line.

Six months from now, could someone redownload exactly this file using only this manifest row — without searching NCBI's website again? If not, what's missing? Also: why record both assembly_accession and assembly_name, when the accession alone is the "real," permanent identifier?

Exercise 4 argued, conceptually, that mapping to a related-but-different species' genome should cost you something. We're not going to take that on faith. Later, once you've learned mapping, you'll take the same raw sequencing dataset from one extant Alcidae species and map it twice:

• (A) against that species' own reference genome
• (B) against a related species' reference genome (the razorbill's situation, simulated on a species where you have the choice)

Then compare mapping rate, MAPQ distribution, properly-paired rate, depth, callable genome size, and heterozygosity between the two. The teaching point: same reads, same pipeline, different reference — any difference you see is reference bias, measured directly rather than assumed. Nothing about that experiment happens in this module; the reference you've chosen, downloaded, and documented here is what makes it possible later.

This also closes the loop on this module's opening point: even with a carefully chosen assembly, a heterospecific reference is still a heterospecific reference. That's precisely why the project will eventually move to raw reads processed through one shared pipeline for the diversity comparisons — it standardises the methodology, even though, as you now know, it can't fully standardise away the fact that the Great Auk's reference isn't its own species.

Write a short prediction in your notes: for the own-reference vs. related-reference experiment above, which direction do you expect mapping rate and heterozygosity to move when using the non-native reference, and why? You're not running anything yet — this is a prediction worth revisiting once the real experiment happens.

What result from that future experiment would make you more confident in the Great Auk's razorbill-mapped heterozygosity estimate? What result would make you less confident?