Raw Sequencing Data

You already know the molecular biology — this is a refresher focused on the steps that matter once you're sitting in front of a terminal rather than a bench:

Before continuing, note in your own words: which of these steps do you think determines how short and how damaged the DNA fragments end up being — and why does that matter for everything downstream, from this module through mapping?

At which step does "the specimen" stop being a unique physical object and become digital data you can copy, share, and re-analyse endlessly without touching the original again?

Just enough to make a FASTQ file meaningful, not a full chemistry lecture:

• DNA is fragmented and ligated to adapters on both ends
• Fragments bind to a flow cell and are locally amplified into clusters — many identical copies, needed so the optical signal is strong enough to detect
• Sequencing by synthesis: fluorescently labelled bases are added cycle by cycle and imaged; each cycle reveals one position
• For paired-end runs, after reading from one end the process re-primes and reads from the other end too (more on this in Exercise 6)
• Each cycle's signal becomes a base call plus a confidence value — the Phred quality score

In one sentence: why is a cluster of many identical copies needed, rather than sequencing one DNA molecule directly?

If the fluorescence signal at a given cycle is weak or ambiguous between two bases, would you expect the resulting Phred quality score to be high or low? Why?

A library is a prepared, adapter-ligated, sequencer-ready collection of DNA molecules derived from a sample. A single specimen can genuinely yield several distinct libraries, for reasons including:

• Multiple independent extractions or preparations
• Different library protocols (e.g. single-stranded vs. double-stranded library chemistry, common in ancient DNA work)
• Independent sequencing batches run at different times
• Deliberate attempts to increase total sequencing depth
• Resequencing

Don't assume which of these applies to any particular real dataset unless the metadata or publication actually says so — that's exactly what Exercise 8 will have you work out.

The real Great Auk genome project (ENA BioProject PRJNA1228327) consists of 13 separately named libraries — G1–G4, GH1–GH4, lib_12, lib_15, lib_20, lib2_13, lib2_17 — all from one specimen (BioSample SAMN47047800).

Before Exercise 8 tries to pin this down carefully, write your own first-pass hypotheses: what are 2–3 plausible, genuine reasons a single historical specimen might have been prepared into 13 separate libraries rather than one?

Does a larger number of libraries automatically mean a larger number of individuals were sampled? Why is that assumption risky?

This distinction matters enough to be explicit, using ENA's own terms:

• Specimen — the physical biological source (here: one museum-vouchered Great Auk)
• Sample (ENA/NCBI BioSample, accession prefix SAMN/SAMEA) — the database record describing that specimen/extraction
• Library (ENA experiment, prefix SRX/ERX) — one prepared, sequenceable DNA pool
• Run (prefix SRR/ERR) — one actual execution of sequencing on a library, producing FASTQ file(s)

In general, this hierarchy can branch at more than one level — illustrated generically below (not real data):

Specimen
└── Sample / extraction
    ├── Library A
    │   ├── Run 1 → FASTQ
    │   └── Run 2 → FASTQ
    └── Library B
        └── Run 3 → FASTQ

Now compare that to the real Great Auk project, which turns out to be flatter than the generic picture above — every library here has exactly one run:

Specimen (NHMD153070)
└── Sample SAMN47047800
    ├── Library G1  → Run SRR32500060 → FASTQ
    ├── Library G2  → Run SRR32500059 → FASTQ
    ├── Library GH1 → Run SRR32500064 → FASTQ
    ├── …
    └── Library lib2_17 → Run SRR32500066 → FASTQ
        (13 libraries total, each with exactly one run)

Match each ENA/NCBI accession prefix to the right level of the hierarchy: PRJNA, SAMN, SRX, SRR. Then, before checking in Exercise 7, say how many of each you'd expect to see for the Great Auk project.

"One specimen is not necessarily one library" and "one library is not necessarily one run" — which of these two statements does the real Great Auk data confirm, and which does it leave untested (because the data happens not to show that particular case)?

FASTQ represents each read as exactly four lines. Here's an illustrative example (not a real Great Auk read — you'll inspect real ones in Exercise 10):

Line 1 starts with @ and identifies the read. Line 2 is the called sequence. Line 3 is almost always a bare + (historically it could repeat the identifier). Line 4 is the same length as line 2 — one character per base — each character encoding a Phred quality score as an ASCII character (higher letters later in the ASCII range generally mean higher quality under the near-universal Phred+33 encoding). Real FASTQ files are almost always distributed as .fastq.gz.

Using the example above, note in your own words which line is which. Then look at the quality line: does it stay constant across the read, or change? Where would you expect quality to be worst, and why?

If you opened a FASTQ file and every single quality character was identical, what would that suggest to you — ordinary real sequencing, or something worth double-checking (e.g. a placeholder or simulated file)?

In paired-end sequencing, a fragment is read from both ends, producing two files per run — R1 and R2. R1 and R2 are mates: the same physical fragment, observed from each end, not two independent reads. Read n in R1 must correspond to read n in R2 — pairing must stay synchronized. This is the common Illumina design, and you will likely meet paired-end data elsewhere in this project. But don't assume it applies here — check.

Check the library_layout field for the Great Auk runs (ENA browser or the filereport API — you'll use both properly in Exercise 7). Is it PAIRED or SINGLE? What does that tell you about how many FASTQ files to expect per run?

A single field in the run metadata answers this directly — you don't need to infer it from file counts.

If this dataset were paired-end instead, what failure mode would you need to watch for that simply doesn't exist for single-end data?

You can browse ENA records through its web interface, but for anything you intend to record in a manifest, its filereport API is more reproducible than copying numbers off a webpage by hand — the same philosophy as downloading by accession in Module 01.

Using the ENA browser page for PRJNA1228327 and the filereport API, build your own table recording: the BioProject and BioSample accessions, the specimen's voucher/locality/collection date, and for each of the 13 runs — run accession, experiment accession, library name, platform, layout, read count, base count, FASTQ file size, and MD5 checksum.

The filereport endpoint takes an accession and a comma-separated list of fields. Construct the URL yourself rather than guessing field names from memory — fields=run_accession,experiment_accession,library_name,library_layout,instrument_platform,read_count,base_count,fastq_ftp,fastq_md5,fastq_bytes is a good starting set.

Does the BioSample accession differ across any of the 13 runs, or is it identical for all of them? What does that tell you, on its own, about how many individual birds this project sequenced?

Now test your Exercise 3 hypotheses against what the metadata you just built in Exercise 7 actually shows. No publication is linked to this project, so your answers have to come entirely from the metadata — say clearly when something can't be determined, rather than guessing.

Using your Exercise 7 table, answer each of these — and for each answer, note whether it's a fact the metadata directly states, or your own reasonable inference:

• Are all 13 runs from the same specimen?
• Are they paired-end or single-end?
• Were they generated on the same sequencing platform?
• Are read lengths identical across runs?
• Are the libraries named in a way that suggests distinct batches?
• Is there a sequencing date recorded per run (as opposed to the specimen's collection date)?
• Is there any circumstantial evidence of resequencing to increase depth?

Read counts and file sizes vary enormously across the 13 runs — that variation itself is a clue, but it's circumstantial, not proof of any single explanation.

Write one short paragraph separating what you've confirmed from the metadata, what you're inferring reasonably but can't prove, and what remains genuinely unknown. You'll want this distinction later if anyone asks why you treated these 13 libraries separately rather than pooling them.

Unlike NCBI's datasets tool in Module 01, ENA needs no special client at all — the filereport gives you direct HTTPS URLs and MD5 checksums, so wget plus md5sum is the whole toolchain. All 13 runs together total roughly 96.5 GB — too much to treat as a single step, and too much to just run on the login node as an ordinary transfer. Mjolnir's storage documentation draws exactly this distinction: everyday transfers ("a few files, actively working") belong on the login node, but a large, long-running transfer has its own dedicated filetransfer partition and QoS — the same one Module 00 mentioned existed without going into detail. It allows up to 4 days of runtime with 1 CPU and 20 GB memory, and is restricted to transfer commands specifically — "running compute work through it is a misuse."

Mjolnir explicitly documents this partition/QoS for transfer commands such as rsync, scp, and sftp — it does not list wget by name. Using wget here is our reasonable extension of that documented transfer pattern — it is a plain network-transfer command, not compute — rather than an example shown verbatim in the documentation. Worth knowing if you ever need to justify it.

Stage 1 — on the login node, download the smallest run (SRR32500061, library GH4, ~1.6 GB) into its own subdirectory, and verify it against its MD5. Stage 2 — use your saved filereport TSV from Exercise 7 to build a script that loops over all 13 runs, giving each its own subdirectory, without concatenating anything. Stage 3 — wrap that script as a Slurm batch job using the filetransfer partition/QoS. Stage 4 — submit it with sbatch so the full download runs independently while you continue with later exercises.

The fastq_ftp field from ENA is given without a URL scheme — you need to prepend https:// yourself. The filetransfer QoS needs its own #SBATCH directives, the same pattern you wrote in Module 00 — just with --partition=filetransfer --qos=filetransfer in place of cpuqueue/normal.

Stage 1's single small file ran directly on the login node; the full 96.5 GB download goes through Slurm instead. What actually changed between the two — is it the duration, the data volume, something Mjolnir defines explicitly, or some combination? Why does Mjolnir bother with a separate QoS for this rather than just telling you to use cpuqueue like any other job?

Module 01 already created an empty raw_reads/ directory alongside reference_genomes/ — this is where it finally gets used. Verification should be streaming and safe: checking a 27 GB file shouldn't require fully decompressing it onto disk.

Confirm all 13 runs are present as their own subdirectories under raw_reads/great_auk/, that each file's size matches what the manifest/filereport reported, that each file passes a gzip integrity check, and — connecting back to Exercise 5 — peek at the first couple of real reads from one file.

gzip -t tests integrity without writing any decompressed output. zcat file.gz | head -n 8 streams and decompresses on the fly, safely, even for the largest file here.

Does the real read identifier line match the anatomy from Exercise 5? If one file failed its gzip integrity check, what would that suggest — a still-incomplete download, a truly corrupt file, or something else — and how would you tell the difference?

Same approach as Module 01's genome manifest, applied to raw reads this time. Here the data are fixed, verified facts you already looked up in Exercise 7 — this exercise is about recording them correctly, not discovering them again.

Create metadata/raw_reads_manifest.tsv with one row per run, covering: species, specimen_id, sample_accession, experiment_accession, run_accession, library_name, platform, library_layout, read_length, fastq_1, file_size_1, md5_1, source, download_date, notes.

If someone asked in six months "did you use all 13 Great Auk runs, or only some of them, and how do you know?" — could this manifest answer that on its own? Separately: why record library_name as well as run_accession, when the run accession alone is the "real" unique identifier?