Feature Request: Keeping the FASTQ quality scores

[JohnMMa]

I understand splitcode has many usages, and some of them involves substituting one sequence with another, such as --assign or sub, will, in a manner of speaking, make the corresponding quality scores irrelevant, since those sequences do not come out of the sequencer.

However, I have a feeling that splitcode changes the quality code of all its output bases to K, including in places where there's no substitution. For example, in the following config file:

@extract <barcode{{a14}}>,<barcode{{bead1}}>,<barcode{{bead2}}>,<barcode{{bead3}}>

tags    ids     groups  locations       distances       previous        next    minFindsG       maxFindsG       exclude
a14_plain.txt$  a14     a14     0:0:8   0       -       {{bead1}}       1       1       0
struct/newBeads/bc1.txt$        bead1   bead1   0:8:18  0       {{a14}} {{bead2}}       1       1       0
struct/newBeads/bc2.txt$        bead2   bead2   0:18:28 0       {{bead1}}       {{bead3}}       1       1       0
struct/newBeads/bc3.txt$        bead3   bead3   0:28:38 0       {{bead2}}       -       1       1       0

There's no substitution involved, just extraction. Yet, splitcode will recode all phred scores in the output to K, i.e from

@A00563:449:HW2C7DMXY:1:1101:17345:1000:ATACATGA
AACAGACAGACGCGATATGAGAGTTCCTAATGTGAGCAATACATGA
+
??????????????????????????????????????????????

to

@A00563:449:HW2C7DMXY:1:1101:17345:1000:ATACATGA
AACAGACAGACGCGATATGAGAGTTCCTAATGTGAGCA
+
KKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK

This makes it hard for us to calculate certain QC metrics such as Q30% in barcode. Is it possible to keep the original file's phred scores if substitution is not involved?

[Yenaled]

Ooh yeah, extraction automatically resets the quality scores. I'll try to have it keep the quality scores in the next release if possible.

[clintval]

I would like this too! All my qualities are K. :)

[clintval]

Here is a bootstrappable script for updating the qualities in the extracted FASTQ from the source FASTQ:

#!/usr/bin/env -S uv run --script
#
# /// script
# requires-python = ">=3.12"
# dependencies = ["pysam==0.23.3", "syncup==0.1.0"]
# [tool.uv]
# exclude-newer = "2025-09-23T00:00:00Z"
# ///
"""Update the FASTQ quality scores from a subsequence FASTQ.

This tool exists because splitcode is unable to retain quality scores when
extracting subsequences from FASTQ files. This script takes a source FASTQ
file and a subsequence FASTQ file (produced by splitcode) and updates the
quality scores in the subsequence file based on the source file.

Upstream issue:

    - https://github.com/pachterlab/splitcode/issues/18

───────
"""

import argparse
import sys
import logging

from logging import Logger
from os import linesep
from pathlib import Path
from typing import cast
from gzip import open as gzopen

from pysam import FastxFile, FastxRecord
from syncup import sync

logging.basicConfig(
    level=logging.INFO,
    datefmt="%Y-%m-%dT%H:%M:%S%z",
    format="%(asctime)s - %(levelname)s - %(name)s - %(message)s",
)
LOGGER: Logger = logging.getLogger(Path(__file__).name)


def update_qualities(in_source: Path, in_subseq: Path, output: Path) -> None:
    """Update the qualities in `in_subseq` based on `in_source` and write to `output`."""
    with (
        FastxFile(str(in_source), persist=False) as source,
        FastxFile(str(in_subseq), persist=False) as subseq,
        gzopen(output, "wt") if output.suffix == ".gz" else open(output, "w") as writer,
    ):
        for count, (fq, fq_subseq) in enumerate(
            sync(
                iter1=source,
                iter2=subseq,
                key1=lambda rec: cast(str, rec.name),
                key2=lambda rec: cast(str, rec.name),
                cmp_func=lambda x, y: (x == y) - 1,  # only advance iter1 when non-equal
            ),
            start=1,
        ):
            assert fq.name == fq_subseq.name, "Names should be equal!"
            assert fq.quality is not None and fq_subseq.quality is not None
            assert fq.sequence is not None and fq_subseq.sequence is not None

            start = fq.sequence.find(fq_subseq.sequence)
            assert start >= 0, "Subsequence not found in source sequence!"
            quals = fq.quality[start : start + len(fq_subseq.sequence)]

            fastx = FastxRecord(
                name=fq_subseq.name,
                comment=fq_subseq.comment,
                sequence=fq_subseq.sequence,
                quality=quals,
            )

            writer.write(f"{fastx}{linesep}")

            if count % 1_000_000 == 0:
                LOGGER.info(f"Completed writing {count} FASTQ records")

        LOGGER.info(f"Completed writing all {count} FASTQ records")


def cli(cli_args: list[str] | None = None) -> int:
    """Update the FASTQ quality scores from a subsequence FASTQ."""
    cli_args = sys.argv[1:] if cli_args is None else cli_args

    parser = argparse.ArgumentParser(
        description=__doc__,
        add_help=True,
        formatter_class=argparse.RawTextHelpFormatter,
        epilog=r"Copyright © Clint Valentine 2025",
    )
    parser.add_argument(
        "--in-source",
        required=True,
        type=Path,
        help="Source FASTQ file.",
    )
    parser.add_argument(
        "--in-subseq",
        required=True,
        type=Path,
        help="Subsequence FASTQ file.",
    )
    parser.add_argument(
        "--output",
        required=True,
        type=Path,
        help="Output repaired FASTQ file.",
    )

    args = parser.parse_args()

    try:
        update_qualities(
            in_source=args.in_source,
            in_subseq=args.in_subseq,
            output=args.output,
        )
    except Exception as exception:
        LOGGER.error(f"Error: {exception}")
        return 1

    return 0


if __name__ == "__main__":
    exit(cli())

Have uv installed and run like:

./update-qualities

In my case, my extracted sequence is larger than half of the original record so Python .find()'s method works great! However, you might need custom logic if you think your extracted sequence may show up more than once in the source FASTQ record.

[Yenaled]

FastQ quality scores are now preserved when doing 'extract' in version 0.31.6