Long-Read Sequencing Files

Long-read sequencing changes what the files carry. The reads are kilobases rather than hundreds of bases, the raw record is a signal rather than a sequence of base calls, and modified bases such as 5-methylcytosine are measured natively alongside the sequence. Two platforms dominate, Oxford Nanopore and PacBio, and they converge on BAM for the aligned output.

Oxford Nanopore: from signal to sequence

A Nanopore run writes its raw record as POD5, an HDF5-derived format holding the electrical current trace of every read. POD5 replaced the older FAST5 format, which was slower to write and read in parallel; both hold the same signal. The basecaller, dorado today, guppy before it, turns the signal into base calls and writes FASTQ or BAM, and the run report travels as a tab-delimited sequencing_summary file with one row per read: read ID, channel, flow cell position, mean Q score, length.

# Inspect the reads in a POD5 file
pod5 view reads.pod5 | head

# Basecall, calling 5mC in CpG context at the same time
dorado basecaller hac,5mCG_5hmCG pod5_dir/ > calls.bam

The basecaller writes BAM directly when it calls modifications, because the modification calls have nowhere to go in FASTQ.

The modified-base BAM tags

Base modifications travel inside the BAM as two per-read tags, defined in the SAMtags specification. MM says which base type is modified and where the modified bases sit, encoded as skips across the candidate bases of the read. ML gives the probability of each called modification, scaled to the integers 0 to 255. The fixture below is a small SAM written the way dorado writes one, with the tags on three of six reads. Rsamtools reads it.

library(Rsamtools)

# Convert the SAM fixture to BAM, as it would arrive from the basecaller.
dir.create("outputs", showWarnings = FALSE, recursive = TRUE)
asBam("/opt/data/methylated.sam", destination = "outputs/methylated",
      overwrite = TRUE)

# Read the alignments back with the modification tags.
param <- ScanBamParam(what = c("qname", "strand", "mapq"),
                      tag = c("MM", "ML"))
reads <- scanBam("outputs/methylated.bam", param = param)[[1]]

data.frame(
  name = reads$qname,
  strand = as.character(reads$strand),
  mapq = reads$mapq,
  mm_tag = ifelse(is.na(reads$tag$MM), "-", reads$tag$MM)
)

[1] "outputs/methylated.bam"
      name strand mapq    mm_tag
1 read_001      +   60 C+m,4,10;
2 read_002      +   60         -
3 read_003      -   45    C+m,2;
4 read_004      +   60  C+m,0,6;
5 read_005      +   25         -
6 read_006      -   60    C+m,8;

Read 1 carries C+m,4,10. The C+m says the tag is about cytosine and the modification is 5mC. The numbers count cytosines to skip on the read: skip four C’s, the next C is modified, skip ten more, the next C is modified. The ML tag then holds one probability per called site.

# The probability of each called modification on the first read.
ml <- reads$tag$ML[[1]]
cat("MM tag:", reads$tag$MM[[1]], "\n")
cat("probabilities (0-255):", ml, "\n")
cat("as fractions:", round(ml / 255, 2), "\n")

MM tag: C+m,4,10;
probabilities (0-255): 230 245
as fractions: 0.9 0.96

A value near 255 is a confident modification and a value near 0 is a confident absence. The per-base calls aggregate into per-CpG counts with modkit pileup, which is where methylation analysis actually starts.

# Aggregate the per-read calls into per-site counts
modkit pileup calls.bam cpg_counts.bed --cpg --ref reference.fa

PacBio: from subreads to HiFi

PacBio’s record of a run is also BAM. The raw reads are subreads.bam, the multiple passes of the polymerase over each circular molecule. The ccs tool collapses the passes of one molecule into a circular consensus read, and a CCS read with enough passes is a HiFi read: long and accurate at once. HiFi reads arrive as BAM or as FASTQ, and the same MM and ML tags carry their modification calls. The older HDF5 containers, .bas.h5 files, survive in archives.

Summary

• Nanopore raw signal lives in POD5 (or the older FAST5); dorado basecalls it to FASTQ or BAM.
• Base modifications have no place in FASTQ. They travel as the MM and ML tags in BAM, per the SAMtags specification.
• MM encodes which bases are modified as skip counts, and ML carries a 0 to 255 probability per called site. Rsamtools reads both.
• PacBio writes subreads.bam, and ccs collapses subreads into accurate HiFi reads, also as BAM.
• Both platforms converge on BAM for aligned output, so samtools and the rest of the alignment toolchain apply unchanged.