Merge pull request #1247 from cfe-lab/fasta-to-fastq

FASTQ generation tools

Author: Donaim

micall/main.py

@@ -88,6 +88,9 @@
     "micall/monitor/update_qai.py",
     "micall/monitor/micall_watcher.py",
     "micall/tcr/igblast.py",
+    "micall/utils/fasta_to_fastq.py",
+    "micall/utils/append_primers.py",
+    "micall/utils/randomize_fastq.py",
 ]
 
 

micall/tests/test_append_primers.py

@@ -0,0 +1,162 @@
+
+import pytest
+from pathlib import Path
+from Bio import SeqIO
+from Bio.Seq import Seq
+from Bio.SeqRecord import SeqRecord
+
+from micall.utils.append_primers import (
+    append_primers_to_record,
+    entry,
+    parse_arguments,
+    append_primers,
+    UserError,
+    DEFAULT_FORWARD_PRIMER,
+    DEFAULT_REVERSE_PRIMER,
+)
+
+
+def test_append_primers_to_record():
+    original_seq = "ATGCGT"
+    record = SeqRecord(Seq(original_seq), id="test", description="dummy")
+    fwd_primer = "AAA"
+    rev_primer = "TTT"
+    new_record = append_primers_to_record(record, fwd_primer, rev_primer)
+    expected_seq = fwd_primer + original_seq + rev_primer
+    assert str(new_record.seq) == expected_seq
+    assert new_record.id == record.id
+    assert new_record.description == record.description
+
+
+def test_parse_arguments_defaults():
+    # simulate command-line arguments,
+    # note that the first two arguments are input and output files.
+    test_args = [
+        "input.fa",
+        "output.fa",
+    ]
+    args = parse_arguments(test_args)
+    # Check that default primers are set.
+    assert args.input_fasta == Path("input.fa")
+    assert args.output_fasta == Path("output.fa")
+    assert args.forward_primer == DEFAULT_FORWARD_PRIMER
+    assert args.reverse_primer == DEFAULT_REVERSE_PRIMER
+
+
+def test_parse_arguments_custom_primers():
+    test_args = [
+        "in.fa",
+        "out.fa",
+        "--forward-primer", "CUSTOMFWD",
+        "--reverse-primer", "CUSTOMREV"
+    ]
+    args = parse_arguments(test_args)
+    assert args.forward_primer == "CUSTOMFWD"
+    assert args.reverse_primer == "CUSTOMREV"
+
+
+@pytest.fixture
+def temp_fasta_file(tmp_path: Path) -> Path:
+    # Create a temporary FASTA file with two records.
+    fasta_path = tmp_path / "input.fa"
+    records = [
+        SeqRecord(Seq("ATGC"), id="rec1", description="first"),
+        SeqRecord(Seq("GATTACA"), id="rec2", description="second"),
+    ]
+    with fasta_path.open("w") as f:
+        SeqIO.write(records, f, "fasta")
+    return fasta_path
+
+
+def test_append_primers_default(temp_fasta_file: Path, tmp_path: Path):
+    # Use default primers.
+    output_path = tmp_path / "output.fa"
+    append_primers(
+        input=temp_fasta_file,
+        output=output_path,
+        forward_primer=DEFAULT_FORWARD_PRIMER,
+        reverse_primer=DEFAULT_REVERSE_PRIMER,
+    )
+    # Now read back the output.
+    records = list(SeqIO.parse(str(output_path), "fasta"))
+    # Check that the primers were appended.
+    for rec in records:
+        # The record sequence should start with DEFAULT_FORWARD_PRIMER
+        # and end with DEFAULT_REVERSE_PRIMER.
+        seq_str = str(rec.seq)
+        assert seq_str.startswith(DEFAULT_FORWARD_PRIMER)
+        assert seq_str.endswith(DEFAULT_REVERSE_PRIMER)
+        # The middle part should be the original sequence. For instance,
+        # for rec1, the original sequence was "ATGC".
+        # We can extract the original: sequence[len(forward): -len(reverse)]
+        original_seq = seq_str[len(DEFAULT_FORWARD_PRIMER): -len(DEFAULT_REVERSE_PRIMER)]
+        # Check that the original sequence exists in the input file.
+        # Without building a dict from the input file, we can simply check that
+        # the extracted sequence is one of the known ones.
+        assert original_seq in {"ATGC", "GATTACA"}
+
+
+def test_append_primers_default_main(temp_fasta_file: Path,
+                                     tmp_path: Path,
+                                     monkeypatch):
+    # Use default primers.
+    output_path = tmp_path / "output.fa"
+
+    monkeypatch.setattr("sys.argv", ["append_primers.py",
+                                     str(temp_fasta_file),
+                                     str(output_path),
+                                     ])
+
+    with pytest.raises(SystemExit) as exinfo:
+        entry()
+
+    ex: SystemExit = exinfo.value
+    assert ex.code == 0
+
+    # Now read back the output.
+    records = list(SeqIO.parse(str(output_path), "fasta"))
+    # Check that the primers were appended.
+    for rec in records:
+        # The record sequence should start with DEFAULT_FORWARD_PRIMER
+        # and end with DEFAULT_REVERSE_PRIMER.
+        seq_str = str(rec.seq)
+        assert seq_str.startswith(DEFAULT_FORWARD_PRIMER)
+        assert seq_str.endswith(DEFAULT_REVERSE_PRIMER)
+        # The middle part should be the original sequence. For instance,
+        # for rec1, the original sequence was "ATGC".
+        # We can extract the original: sequence[len(forward): -len(reverse)]
+        original_seq = seq_str[len(DEFAULT_FORWARD_PRIMER): -len(DEFAULT_REVERSE_PRIMER)]
+        # Check that the original sequence exists in the input file.
+        # Without building a dict from the input file, we can simply check that
+        # the extracted sequence is one of the known ones.
+        assert original_seq in {"ATGC", "GATTACA"}
+
+
+def test_append_primers_custom(temp_fasta_file: Path, tmp_path: Path):
+    output_path = tmp_path / "output_custom.fa"
+    custom_fwd = "CUSTOMFWD"
+    custom_rev = "CUSTOMREV"
+    append_primers(
+        input=temp_fasta_file,
+        output=output_path,
+        forward_primer=custom_fwd,
+        reverse_primer=custom_rev,
+    )
+    records = list(SeqIO.parse(str(output_path), "fasta"))
+    for rec in records:
+        seq_str = str(rec.seq)
+        assert seq_str.startswith(custom_fwd)
+        assert seq_str.endswith(custom_rev)
+
+
+def test_append_primers_missing_input(tmp_path: Path):
+    # Create a non-existent input file path.
+    missing_file = tmp_path / "non_existent.fa"
+    output_file = tmp_path / "output.fa"
+    with pytest.raises(UserError):
+        append_primers(
+            input=missing_file,
+            output=output_file,
+            forward_primer=DEFAULT_FORWARD_PRIMER,
+            reverse_primer=DEFAULT_REVERSE_PRIMER,
+        )

micall/tests/test_randomize_fastq.py

@@ -0,0 +1,199 @@
+import random
+import pytest
+from pathlib import Path
+from io import StringIO
+
+from Bio import SeqIO
+from Bio.Seq import Seq
+from Bio.SeqRecord import SeqRecord
+
+from micall.utils.randomize_fastq import (
+    introduce_errors,
+    process_records,
+    process_fastq,
+    get_parser,
+    main,
+    entry,
+    NUCLEOTIDES,
+)
+
+
+# --- Tests for introduce_errors ---
+def test_introduce_errors_no_error():
+    seq = "ACGTACGT"
+    qualities = [30] * len(seq)
+    # no errors: all error rates set to 0.
+    new_seq, new_quals = introduce_errors(seq, qualities,
+                                          subst_rate=0.0,
+                                          ins_rate=0.0,
+                                          del_rate=0.0,
+                                          ins_quality=20)
+    assert new_seq == seq
+    assert new_quals == qualities
+
+
+def test_introduce_errors_full_deletion():
+    seq = "ACGTACGT"
+    qualities = [30] * len(seq)
+    new_seq, new_quals = introduce_errors(seq, qualities,
+                                          subst_rate=0.0,
+                                          ins_rate=0.0,
+                                          del_rate=1.0,
+                                          ins_quality=20)
+    # With full deletion, expect empty sequence and empty quality list.
+    assert new_seq == ""
+    assert new_quals == []
+
+
+def test_introduce_errors_full_substitution():
+    seq = "ACGT"
+    qualities = [30, 31, 32, 33]
+    # Set subst_rate=1 so every base is substituted.
+    new_seq, new_quals = introduce_errors(seq, qualities,
+                                          subst_rate=1.0,
+                                          ins_rate=0.0,
+                                          del_rate=0.0,
+                                          ins_quality=20)
+    assert len(new_seq) == len(seq)
+    # for each base, new base should be different than original.
+    for original, new in zip(seq, new_seq):
+        assert new.upper() in NUCLEOTIDES
+        assert new.upper() != original.upper()
+    assert new_quals == qualities
+
+
+def test_introduce_errors_full_insertion():
+    seq = "ACGT"
+    qualities = [30] * len(seq)
+    new_seq, new_quals = introduce_errors(seq, qualities,
+                                          subst_rate=0.0,
+                                          ins_rate=1.0,
+                                          del_rate=0.0,
+                                          ins_quality=15)
+    # Every base should be followed by an inserted base: length 2*len(seq)
+    assert len(new_seq) == 2 * len(seq)
+    # And the bases at even indices (0-indexed) should be the original bases.
+    for i, base in enumerate(seq):
+        assert new_seq[2*i] == base
+    # Check inserted qualities:
+    # The inserted qualities (at odd indices) should equal ins_quality.
+    for i in range(len(seq)):
+        assert new_quals[2*i] == qualities[i]
+        assert new_quals[2*i+1] == 15
+
+
+# --- Test for process_records using an in-memory FASTQ string ---
+def test_process_records_no_error():
+    fastq_str = (
+        "@read1\n"
+        "ACGT\n"
+        "+\n"
+        "!!!!\n"
+    )
+
+    # in FASTQ, qualities are usually encoded by ascii, but BioPython converts
+    # them to a list of phred quality scores.
+    handle = StringIO(fastq_str)
+    # All error rates set to 0, so no changes.
+    records = list(process_records(handle,
+                                   subst_rate=0.0,
+                                   ins_rate=0.0,
+                                   del_rate=0.0,
+                                   ins_quality=20,
+                                   rng=random.Random(42)))
+    assert len(records) == 1
+    rec = records[0]
+    assert str(rec.seq) == "ACGT"
+
+    # Check phred qualities:
+    # '!' converts to 0 in Phred.
+    assert rec.letter_annotations["phred_quality"] == [0]*4
+
+
+# --- Tests for argument parser ---
+def test_get_parser_defaults():
+    parser = get_parser()
+    args = parser.parse_args(["in.fastq", "out.fastq"])
+    assert args.in_fastq == Path("in.fastq")
+    assert args.out_fastq == Path("out.fastq")
+    # default error rates from the parser.
+    assert args.subst_rate == 0.005
+    assert args.ins_rate == 0.0
+    assert args.del_rate == 0.0
+    assert args.ins_quality == 20
+    assert args.seed == 42
+
+
+@pytest.fixture
+def tmp_fastq_file(tmp_path: Path) -> Path:
+    fastq_path = tmp_path / "input.fastq"
+    records = [
+        SeqRecord(Seq("ACGTACGT"),
+                  id="read1",
+                  description="",
+                  letter_annotations={"phred_quality": [30]*8}),
+    ]
+    with fastq_path.open("w") as fh:
+        SeqIO.write(records, fh, "fastq")
+    return fastq_path
+
+
+def test_process_fastq_no_error(tmp_fastq_file: Path, tmp_path: Path):
+    output_path = tmp_path / "output.fastq"
+    # Run process_fastq with error rates set to 0.
+    process_fastq(in_fastq=tmp_fastq_file,
+                  out_fastq=output_path,
+                  subst_rate=0.0,
+                  ins_rate=0.0,
+                  del_rate=0.0,
+                  ins_quality=20,
+                  rng=random.Random(42))
+    # Read back the output file.
+    recs = list(SeqIO.parse(str(output_path), "fastq"))
+    assert len(recs) == 1
+    # Without errors, the sequence should be identical.
+    assert str(recs[0].seq) == "ACGTACGT"
+
+
+def test_main_integration(tmp_fastq_file: Path, tmp_path: Path, monkeypatch):
+    # Build a fake command line to run main.
+    out_path = tmp_path / "output2.fastq"
+    args = [
+        str(tmp_fastq_file),
+        str(out_path),
+        "--subst_rate", "0.0",
+        "--ins_rate", "0.0",
+        "--del_rate", "0.0",
+        "--ins_quality", "20",
+        "--seed", "42"
+    ]
+    # Call main directly.
+    exit_code = main(args)
+    assert exit_code == 0
+    recs = list(SeqIO.parse(str(out_path), "fastq"))
+    assert len(recs) == 1
+    assert str(recs[0].seq) == "ACGTACGT"
+
+
+def test_entry_integration(tmp_fastq_file: Path, tmp_path: Path, monkeypatch):
+    out_path = tmp_path / "output3.fastq"
+    monkeypatch.setattr("sys.argv", [
+        "error_introducer.py",
+        str(tmp_fastq_file),
+        str(out_path),
+        "--subst_rate", "0.0",
+        "--ins_rate", "0.0",
+        "--del_rate", "0.0",
+        "--ins_quality", "20",
+        "--seed", "42"
+    ])
+
+    # entry() calls sys.exit; we can catch the SystemExit exception.
+    with pytest.raises(SystemExit) as e:
+        entry()
+
+    assert e.value.code == 0
+
+    recs = list(SeqIO.parse(str(out_path), "fastq"))
+    assert len(recs) == 1
+    assert str(recs[0].seq) == "ACGTACGT"