Create convert_fast5_to_blow5.py

Add fast5 convert script (for signal only fast5)

Author: Runsheng

convert_fast5_to_blow5.py

@@ -0,0 +1,287 @@
+#!/usr/bin/env python
+"""Convert Oxford Nanopore FAST5 files to BLOW5 using pyslow5.
+
+This script is intended for FAST5 files that only contain raw signal data
+(current signal) and minimal metadata. Missing header attributes are
+populated with a user configurable default string ("CycloneSeq" by default)
+so the resulting BLOW5 can be processed by slow5tools and downstream tools.
+
+Example usage:
+
+    python convert_fast5_to_blow5.py \
+        --input fast5_WGA_p \
+        --output blow5_out
+
+Requirements: pyslow5, h5py
+"""
+
+from __future__ import annotations
+
+import argparse
+import logging
+import sys
+from pathlib import Path
+from typing import Dict, Iterable, Optional, Tuple
+
+import h5py
+import numpy as np
+import pyslow5
+
+
+LOGGER = logging.getLogger(__name__)
+
+
+def parse_args(argv: Optional[Iterable[str]] = None) -> argparse.Namespace:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument(
+        "--input",
+        required=True,
+        type=Path,
+        help="Path to a FAST5 file or directory containing FAST5 files.",
+    )
+    parser.add_argument(
+        "--output",
+        required=True,
+        type=Path,
+        help="Output directory for generated .blow5 files.",
+    )
+    parser.add_argument(
+        "--default-meta",
+        default="CycloneSeq",
+        help="Fallback string used when a header attribute is missing.",
+    )
+    parser.add_argument(
+        "--record-compression",
+        default="zlib",
+        choices=["none", "zlib", "zstd"],
+        help="Record compression codec for the BLOW5 container.",
+    )
+    parser.add_argument(
+        "--signal-compression",
+        default="svb_zd",
+        choices=["none", "svb_zd", "svb-zd", "ex_zd", "ex-zd"],
+        help="Signal compression codec for the raw signal array.",
+    )
+    parser.add_argument(
+        "--verbose",
+        action="store_true",
+        help="Enable verbose logging for debugging conversions.",
+    )
+    return parser.parse_args(argv)
+
+
+def discover_fast5_files(path: Path) -> Iterable[Path]:
+    if path.is_file():
+        if path.suffix != ".fast5":
+            raise ValueError(f"Input file must have .fast5 extension: {path}")
+        yield path
+        return
+
+    if not path.is_dir():
+        raise FileNotFoundError(f"Input path does not exist: {path}")
+
+    for fast5_path in sorted(path.rglob("*.fast5")):
+        if fast5_path.is_file():
+            yield fast5_path
+
+
+def decode_attr(value) -> Optional[str]:
+    if value is None:
+        return None
+    if isinstance(value, (bytes, np.bytes_)):
+        return value.decode("utf-8", errors="ignore")
+    if isinstance(value, np.generic):
+        return str(value.item())
+    return str(value)
+
+
+def gather_header_metadata(read_group) -> Dict[str, str]:
+    metadata: Dict[str, str] = {}
+
+    tracking = read_group.get("tracking_id")
+    if tracking is not None:
+        for key, value in tracking.attrs.items():
+            metadata[key] = decode_attr(value) or ""
+
+    context = read_group.get("context_tags")
+    if context is not None:
+        for key, value in context.attrs.items():
+            metadata[key] = decode_attr(value) or ""
+
+    channel = read_group.get("channel_id")
+    if channel is not None:
+        for key, value in channel.attrs.items():
+            metadata[key] = decode_attr(value) or ""
+        if "sampling_rate" in channel.attrs:
+            sample_rate = decode_attr(channel.attrs.get("sampling_rate"))
+            metadata.setdefault("sample_frequency", sample_rate)
+            metadata.setdefault("sampling_rate", sample_rate)
+        if "digitisation" in channel.attrs:
+            metadata.setdefault("digitisation", decode_attr(channel.attrs.get("digitisation")))
+        if "offset" in channel.attrs:
+            metadata.setdefault("offset", decode_attr(channel.attrs.get("offset")))
+        if "range" in channel.attrs:
+            metadata.setdefault("range", decode_attr(channel.attrs.get("range")))
+
+    return metadata
+
+
+def build_header(writer: pyslow5.Open, metadata: Dict[str, str], default_value: str) -> Dict[str, str]:
+    header = writer.get_empty_header()
+    for key in header:
+        value = metadata.get(key)
+        if value is None or value == "":
+            if key in {"exp_start_time", "protocol_start_time"}:
+                value = "1970-01-01T00:00:00Z"
+            elif key in {"sample_frequency", "digitisation", "offset", "range", "sampling_rate"}:
+                value = metadata.get(key, metadata.get("sampling_rate", "0"))
+            else:
+                value = default_value
+        header[key] = str(value)
+    header["slow5_version"] = "1.0.0"
+    header["num_read_groups"] = "1"
+    header.setdefault("run_id", default_value)
+    header.setdefault("exp_start_time", "1970-01-01T00:00:00Z")
+    header.setdefault("flow_cell_id", default_value)
+    header.setdefault("sample_id", default_value)
+    return header
+
+
+def read_signal_group(group) -> Dict[str, np.ndarray]:
+    raw = group.get("Raw")
+    if raw is None or "Signal" not in raw:
+        raise KeyError("Raw/Signal dataset is missing in FAST5 read group")
+    signal = raw["Signal"][:]
+    if signal.dtype != np.int16:
+        signal = signal.astype(np.int16)
+    return {
+        "signal": signal,
+        "len_raw_signal": int(signal.shape[0]),
+    }
+
+
+def _safe_int(value, default: int = 0) -> int:
+    try:
+        return int(value)
+    except Exception:
+        return default
+
+
+def build_record(
+    writer: pyslow5.Open,
+    read_id: str,
+    group,
+    default_group: int = 0,
+) -> Tuple[Dict[str, object], Dict[str, object]]:
+    record = writer.get_empty_record()
+    channel_attrs = group["channel_id"].attrs if "channel_id" in group else {}
+    raw_info = read_signal_group(group)
+    record["read_id"] = read_id
+    record["read_group"] = default_group
+    record["digitisation"] = float(channel_attrs.get("digitisation", 0.0))
+    record["offset"] = float(channel_attrs.get("offset", 0.0))
+    record["range"] = float(channel_attrs.get("range", 0.0))
+    record["sampling_rate"] = float(channel_attrs.get("sampling_rate", 0.0))
+    record["len_raw_signal"] = int(raw_info["len_raw_signal"])
+    record["signal"] = np.ascontiguousarray(raw_info["signal"], dtype=np.int16)
+
+    aux: Dict[str, object] = {}
+
+    if "channel_id" in group:
+        channel_group = group["channel_id"].attrs
+        if "channel_number" in channel_group:
+            channel_number = channel_group.get("channel_number")
+            if isinstance(channel_number, (bytes, np.bytes_)):
+                aux["channel_number"] = channel_number.decode("utf-8", errors="ignore")
+            else:
+                aux["channel_number"] = str(channel_number)
+
+    raw_attrs = group["Raw"].attrs if "Raw" in group else {}
+    if "median_before" in raw_attrs:
+        aux["median_before"] = float(raw_attrs.get("median_before"))
+    if "read_number" in raw_attrs:
+        aux["read_number"] = _safe_int(raw_attrs.get("read_number"))
+    if "start_mux" in raw_attrs:
+        aux["start_mux"] = max(0, min(255, _safe_int(raw_attrs.get("start_mux"))))
+    if "start_time" in raw_attrs:
+        aux["start_time"] = max(0, _safe_int(raw_attrs.get("start_time")))
+
+    return record, aux
+
+
+def convert_fast5_to_blow5(
+    fast5_path: Path,
+    output_blow5: Path,
+    default_meta: str,
+    record_compression: str,
+    signal_compression: str,
+) -> None:
+    LOGGER.info("Converting %s -> %s", fast5_path, output_blow5)
+    output_blow5.parent.mkdir(parents=True, exist_ok=True)
+
+    writer: Optional[pyslow5.Open] = None
+
+    with h5py.File(fast5_path, "r") as fast5_file:
+        read_ids = list(fast5_file.keys())
+        if not read_ids:
+            LOGGER.warning("No read groups found in %s; skipping", fast5_path)
+            return
+
+        first_group = fast5_file[read_ids[0]]
+        metadata = gather_header_metadata(first_group)
+
+        writer = pyslow5.Open(str(output_blow5), "w", record_compression, signal_compression)
+        header = build_header(writer, metadata, default_meta)
+        writer.write_header(header)
+
+        for read_id in read_ids:
+            group = fast5_file[read_id]
+            record, aux = build_record(writer, read_id, group)
+            if aux:
+                writer.write_record(record, aux)
+            else:
+                writer.write_record(record)
+
+    if writer is not None:
+        writer.close()
+
+
+def main(argv: Optional[Iterable[str]] = None) -> int:
+    args = parse_args(argv)
+
+    logging.basicConfig(
+        level=logging.DEBUG if args.verbose else logging.INFO,
+        format="%(asctime)s - %(levelname)s - %(message)s",
+    )
+
+    try:
+        fast5_files = list(discover_fast5_files(args.input))
+    except Exception as exc:  # pragma: no cover - surfaced to user
+        LOGGER.error("Failed to discover FAST5 files: %s", exc)
+        return 1
+
+    if not fast5_files:
+        LOGGER.error("No FAST5 files found under %s", args.input)
+        return 1
+
+    for fast5_path in fast5_files:
+        output_blow5 = args.output / (fast5_path.stem + ".blow5")
+        try:
+            convert_fast5_to_blow5(
+                fast5_path,
+                output_blow5,
+                args.default_meta,
+                args.record_compression,
+                args.signal_compression,
+            )
+        except Exception as exc:
+            LOGGER.error("Failed to convert %s: %s", fast5_path, exc)
+            return 1
+
+    LOGGER.info("Conversion completed for %d FAST5 file(s)", len(fast5_files))
+    return 0
+
+
+if __name__ == "__main__":  # pragma: no cover - CLI entry point
+    sys.exit(main())
+