AurinVos · Jun 8, 2026 · Jun 8, 2026 · Jun 9, 2026 · Jul 2, 2026
diff --git a/PerlLib/GFF3_utils2.pm b/PerlLib/GFF3_utils2.pm
diff --git a/testing/compare_gff3_file_to_proteins_memory.sh b/testing/compare_gff3_file_to_proteins_memory.sh
@@ -0,0 +1,66 @@
+#!/usr/bin/env bash
+set -euo pipefail
+
+if [[ $# -ne 2 ]]; then
+    cat >&2 <<'USAGE'
+Usage: testing/compare_gff3_file_to_proteins_memory.sh <num_records> <seq_len>
+
+Generates synthetic transcript FASTA/GFF3 data and compares maximum resident
+set size for two streaming batch sizes. Requires /usr/bin/time with the -v
+option.
+USAGE
+    exit 2
+fi
+
+num_records="$1"
+seq_len="$2"
+
+if [[ ! -x /usr/bin/time ]]; then
+    echo "Error: /usr/bin/time is required for memory comparison" >&2
+    exit 1
+fi
+
+repo_dir="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)"
+workdir="$(mktemp -d "${TMPDIR:-/tmp}/gff3_to_proteins_memory.XXXXXX")"
+trap 'rm -rf "$workdir"' EXIT
+
+python3 - "$num_records" "$seq_len" "$workdir/test.fa" "$workdir/test.gff3" <<'PY'
+import sys
+num_records = int(sys.argv[1])
+seq_len = int(sys.argv[2])
+fasta = sys.argv[3]
+gff3 = sys.argv[4]
+if seq_len < 6:
+    raise SystemExit("seq_len must be at least 6")
+seq = "ATG" + ("A" * (seq_len - 6)) + "TAA"
+with open(fasta, "w") as fa, open(gff3, "w") as gf:
+    gf.write("##gff-version 3\n")
+    for i in range(1, num_records + 1):
+        tid = f"tx{i}"
+        gid = f"gene{i}"
+        mid = f"mrna{i}"
+        fa.write(f">{tid}\n{seq}\n")
+        gf.write(f"{tid}\tTD\tgene\t1\t{seq_len}\t.\t+\t.\tID={gid}\n")
+        gf.write(f"{tid}\tTD\tmRNA\t1\t{seq_len}\t.\t+\t.\tID={mid};Parent={gid}\n")
+        gf.write(f"{tid}\tTD\texon\t1\t{seq_len}\t.\t+\t.\tID=exon{i};Parent={mid}\n")
+        gf.write(f"{tid}\tTD\tCDS\t1\t{seq_len}\t.\t+\t0\tID=cds{i};Parent={mid}\n")
+PY
+
+/usr/bin/time -v "$repo_dir/util/gff3_file_to_proteins.pl" \
+    --batch_size 1 \
+    --gff3 "$workdir/test.gff3" \
+    --fasta "$workdir/test.fa" \
+    > "$workdir/batch_1.pep" 2> "$workdir/batch_1.time"
+
+/usr/bin/time -v "$repo_dir/util/gff3_file_to_proteins.pl" \
+    --batch_size 1000 \
+    --gff3 "$workdir/test.gff3" \
+    --fasta "$workdir/test.fa" \
+    > "$workdir/batch_1000.pep" 2> "$workdir/batch_1000.time"
+
+diff -u "$workdir/batch_1.pep" "$workdir/batch_1000.pep"
+
+echo "Streaming mode memory (--batch_size 1):"
+grep "Maximum resident set size" "$workdir/batch_1.time"
+echo "Streaming mode memory (--batch_size 1000):"
+grep "Maximum resident set size" "$workdir/batch_1000.time"
diff --git a/testing/test_gff3_file_to_proteins_streaming.sh b/testing/test_gff3_file_to_proteins_streaming.sh
@@ -0,0 +1,86 @@
+#!/usr/bin/env bash
+set -euo pipefail
+
+repo_dir="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)"
+workdir="$(mktemp -d "${TMPDIR:-/tmp}/gff3_to_proteins_streaming.XXXXXX")"
+trap 'rm -rf "$workdir"' EXIT
+
+cat > "$workdir/test.fa" <<'FASTA'
+>tx1 first transcript
+ATGGCCGCTTAA
+>tx2 second transcript
+ATGAAACCCGGGTAA
+>tx3 third transcript
+ATGTTTCCCTAG
+FASTA
+
+cat > "$workdir/test.gff3" <<'GFF3'
+##gff-version 3
+tx1	TD	gene	1	12	.	+	.	ID=gene1;Name=gene one
+tx1	TD	mRNA	1	12	.	+	.	ID=mrna1;Parent=gene1
+tx1	TD	exon	1	12	.	+	.	ID=exon1;Parent=mrna1
+tx1	TD	CDS	1	12	.	+	0	ID=cds1;Parent=mrna1
+tx2	TD	gene	1	15	.	+	.	ID=gene2;Name=gene two
+tx2	TD	mRNA	1	15	.	+	.	ID=mrna2;Parent=gene2
+tx2	TD	exon	1	15	.	+	.	ID=exon2;Parent=mrna2
+tx2	TD	CDS	1	15	.	+	0	ID=cds2;Parent=mrna2
+tx3	TD	gene	1	12	.	+	.	ID=gene3;Name=gene three
+tx3	TD	mRNA	1	12	.	+	.	ID=mrna3;Parent=gene3
+tx3	TD	exon	1	12	.	+	.	ID=exon3;Parent=mrna3
+tx3	TD	CDS	1	12	.	+	0	ID=cds3;Parent=mrna3
+GFF3
+
+for seq_type in prot CDS cDNA gene; do
+    "$repo_dir/util/gff3_file_to_proteins.pl" \
+        --gff3 "$workdir/test.gff3" \
+        --fasta "$workdir/test.fa" \
+        --seqType "$seq_type" \
+        > "$workdir/default.${seq_type}"
+
+    for batch_size in 1 2 5000; do
+        "$repo_dir/util/gff3_file_to_proteins.pl" \
+            --batch_size "$batch_size" \
+            --gff3 "$workdir/test.gff3" \
+            --fasta "$workdir/test.fa" \
+            --seqType "$seq_type" \
+            > "$workdir/streaming.${seq_type}.${batch_size}"
+
+        diff -u "$workdir/default.${seq_type}" "$workdir/streaming.${seq_type}.${batch_size}"
+    done
+done
+
+if "$repo_dir/util/gff3_file_to_proteins.pl" \
+    --legacy_memory_mode \
+    --gff3 "$workdir/test.gff3" \
+    --fasta "$workdir/test.fa" \
+    > "$workdir/legacy_mode.out" 2> "$workdir/legacy_mode.err"; then
+    echo "Expected removed --legacy_memory_mode option to fail" >&2
+    exit 1
+fi
+
+if ! grep -q "Unknown option: legacy_memory_mode" "$workdir/legacy_mode.err"; then
+    echo "Removed --legacy_memory_mode option was not rejected as expected" >&2
+    cat "$workdir/legacy_mode.err" >&2
+    exit 1
+fi
+
+cat > "$workdir/missing.fa" <<'FASTA'
+>tx1 first transcript
+ATGGCCGCTTAA
+FASTA
+
+if "$repo_dir/util/gff3_file_to_proteins.pl" \
+    --gff3 "$workdir/test.gff3" \
+    --fasta "$workdir/missing.fa" \
+    > "$workdir/missing.out" 2> "$workdir/missing.err"; then
+    echo "Expected missing FASTA entry check to fail" >&2
+    exit 1
+fi
+
+if ! grep -q "no FASTA entry found for GFF3 sequence id: tx2" "$workdir/missing.err"; then
+    echo "Missing FASTA entry error was not reported as expected" >&2
+    cat "$workdir/missing.err" >&2
+    exit 1
+fi
+
+echo "gff3_file_to_proteins streaming tests passed"
diff --git a/util/PWM/build_atgPWM_+-.pl b/util/PWM/build_atgPWM_+-.pl
@@ -6,7 +6,7 @@
 use Getopt::Long qw(:config posix_default no_ignore_case bundling pass_through);
 use FindBin;
 use lib ("$FindBin::Bin/../../PerlLib");
-use Fasta_reader;
+use Fasta_retriever;
 use Nuc_translator;
 use PWM;
 
@@ -66,12 +66,9 @@
 
 main: {
 
-    my $fasta_reader = new Fasta_reader($transcripts_file);
-
-    my %seqs = $fasta_reader->retrieve_all_seqs_hash();
-
     my @starts = &parse_starts($selected_orfs_file);
 
+    my $fasta_retriever = new Fasta_retriever($transcripts_file);
 
 
     my $pwm_plus = new PWM();
@@ -88,7 +85,7 @@
     foreach my $start_info (@starts) {
         my ($transcript_acc, $start_coord, $orient) = @$start_info;
 
-        my $transcript_seq = uc $seqs{$transcript_acc};
+        my $transcript_seq = uc $fasta_retriever->get_seq($transcript_acc);
 
         if ($orient eq '-') {
             # convert info to (+) strand

diff --git a/util/TransDecoder.Predict b/util/TransDecoder.Predict
@@ -223,6 +223,8 @@ main: {
     unless (-d $workdir) {
         die "Error, cannot find directory: $workdir,  be sure to first run TransDecoder.LongOrfs before TransDecoder.Predict\n\n";
     }
+
+    &ensure_fasta_fai_if_possible($transcripts_file);
 
     my $checkpoints_dir = "$workdir/__checkpoints_TDpredict";
     if (! -d $checkpoints_dir) {
@@ -366,6 +368,58 @@ main: {
 }
 
 
+####
+sub ensure_fasta_fai_if_possible {
+    my ($fasta_file) = @_;
+
+    my $fai_file = "$fasta_file.fai";
+    if (-s $fai_file) {
+        return;
+    }
+
+    my $samtools = &find_program_on_path("samtools");
+    unless ($samtools) {
+        warn "WARNING: samtools not found on PATH; cannot create FASTA index $fai_file. Continuing without it.\n";
+        return;
+    }
+
+    print STDERR "CMD: $samtools faidx $fasta_file\n";
+    my $ret = system($samtools, "faidx", $fasta_file);
+
+    if ($ret) {
+        my $failure = ($ret == -1) ? "failed to execute: $!" : "exit status " . ($ret >> 8);
+        warn "WARNING: samtools faidx failed for $fasta_file ($failure). Continuing without FASTA index $fai_file.\n";
+    }
+
+    return;
+}
+
+####
+sub find_program_on_path {
+    my ($prog) = @_;
+
+    my @extensions = ("");
+    if ($^O =~ /^(MSWin32|cygwin|msys)$/i) {
+        push(@extensions, split(/;/, $ENV{PATHEXT} || ".COM;.EXE;.BAT;.CMD"));
+    }
+
+    my %seen;
+    @extensions = grep { ! $seen{lc($_)}++ } @extensions;
+
+    foreach my $dir (File::Spec->path()) {
+        next unless defined $dir && length $dir;
+
+        foreach my $extension (@extensions) {
+            my $path = File::Spec->catfile($dir, $prog . $extension);
+            if (-f $path && -x _) {
+                return($path);
+            }
+        }
+    }
+
+    return;
+}
+
 ####
 sub process_cmd {
 	my ($cmd) = @_;