Skip to content

FASTA reads, blank lines and mate names as STAR; spill folder; --outBAMcompression -1 - #358

Open
BenjaminDEMAILLE wants to merge 2 commits into
claude/brave-hamilton-c5d377from
fix/input-output-small
Open

BenjaminDEMAILLE wants to merge 2 commits into
claude/brave-hamilton-c5d377from
fix/input-output-small

Conversation

@BenjaminDEMAILLE

@BenjaminDEMAILLE BenjaminDEMAILLE commented Nov 2, 2026 •

Copy link
Copy Markdown
Contributor

Part of #354 (rest of finding 17, finding 18, items 28 and 29).

Stacked on #353, which is stacked on #352, because all of this touches next_in_current_file and #352 already adds the quality-length check there. Only the last commit is this PR. I will rebase it on main once those two are merged.

1. FASTA reads

Problem

--readFilesIn r.fa stopped with invalid name prefix, and no file or read was named, while STAR accepts FASTA.

What STAR does

Observed with Homebrew STAR 2.7.11b:

  • A sequence may span several lines, and the name is cut at the first space and at /.
  • FASTA and FASTQ records may mix in one file.
  • The record's QUAL is *.
  • --outReadsUnmapped Fastx writes >name and the sequence.

Change

  • next_in_current_file peeks at the first byte, and > goes to the new next_fasta. Quality stays empty.
  • Four places sliced quality by the clip lengths, which panics on an empty slice: clip_read and three in src/io/sam.rs, which now use aligned_quality.
  • UnmappedFastqWriter::write_record writes FASTA when there is no quality.

Measured against STAR

Same index contents, 202 FASTA reads (forward, reverse, one random, one 3 bases long) and 101 FASTA pairs. SAM fields 1-11, sorted:

Run SAM Exit codes
SE, --outSAMunmapped Within identical 0 / 0
PE, --outSAMunmapped Within identical 0 / 0
--clip5pNbases 3 --clip3pNbases 4 identical 0 / 0
--clip3pAdapterSeq ACGTAC identical 0 / 0
--outSAMmode NoQS identical 0 / 0
--outReadsUnmapped Fastx, SE identical 0 / 0
--outReadsUnmapped Fastx, PE differs by one record, see below 0 / 0

BAM Unsorted SortedByCoordinate and --chimOutType WithinBAM also ran; samtools quickcheck passes on those files and the unsorted BAM has STAR's 200 records. I did not compare their contents field by field.

The one PE difference is a half-mapped pair's extra unmapped-mate record, which has nothing to do with FASTA. It is finding 2 of #354 and happens with FASTQ too.

Not the same as STAR

  • STAR's Unmapped.out.mate1 names are >junk 0:N: ; here they are >junk. This also happens with FASTQ on main (@q9 0:N: against @q9), so I left it out of this PR and added it to Tracking: fixes from the whole-repo audit (main at 32fbbc1) #354.
  • No --soloType run was made with FASTA. A barcode read with no quality gets empty Q30 counts.

2. Blank lines

Problem

One blank line after the last record made the whole run exit 1 with invalid name prefix. It only surfaced at the end of the input, so a long run died at the very end.

What STAR does

It stops at the first line starting with a newline or a space, and exits 0. So a trailing blank line is harmless. But @q1 ... \n\n@q2 ... maps q1, ignores q2, and exits 0.

Change

Blank lines at the end of a file end it, including \r\n and spaces. Anything after a blank line is an error that names the file.

Deliberate divergence: DIVERGENCE.md §4.8

Reads after a blank line would be lost with exit code 0, so this stops the run where STAR carries on. This needs maintainer sign-off. Following STAR is a two-line change.

3. Mate names and QNAME

Problem

@SRR1.5.1 and @SRR1.5.2 were refused with read names do not match, while STAR accepts them.

While testing that, I found a second difference. PairedRead.name was the stripped base name, so @p_1 / @p_2 got QNAME p. STAR writes mate 1's name.

Input STAR main This branch
p_1 / p_2 QNAME p_1 p p_1
S.0.1 / S.0.2 QNAME S.0.1 exit 1 S.0.1

Change

  • strip_mate_suffix strips .1 / .2 only when another . comes before it. SRR1.1 / SRR1.2 are still refused: there the number is the read's, and stripping it would make every out-of-step pair pass the check.
  • PairedRead.name is mate 1's name.

STAR has no such check at all: it exits 0 on mates that are out of step. The check here is older than this PR and is useful, so I kept it. It should have a DIVERGENCE.md entry; I did not add one, since this PR only narrows it.

4. Sort spill folder

Path::new("out/sub/").parent() is out, so --outFileNamePrefix /data/run1/ spilled into /data, against DIVERGENCE.md §4.3 ("beside the output"). A small or read-only parent would fail the run at its first spill. sort_temp_dir now returns the prefix itself when it ends in a separator.

5. --outBAMcompression -1

bgzf_compression mapped -1 to no compression, and its comment said that is what STAR does. Measured on 3000 reads:

Level STAR main This branch
-1 59,925 B uncompressed 45,935 B
0 635,969 B 636,292 B 636,292 B
1 59,924 B 45,928 B 45,928 B

So STAR's -1 is its level 1. Our level 1 is smaller than STAR's because of the deflate library; that is older than this PR. STAR's manual calls -1 "default compression", which is not what the binary does, so I followed the binary and said so in the comment.

Tests

  • src/io/fastq.rs: fasta_reads_have_empty_quality_and_may_span_lines, blank_lines_end_a_file_but_not_before_more_reads, clip_read_of_a_fasta_read_has_no_quality, dotted_mate_suffix_needs_a_read_number_before_it, pair_is_named_after_mate1_as_star_does
  • src/io/bam.rs: minus_one_compresses_like_level_one_and_zero_does_not, spill_files_stay_inside_a_directory_prefix

cargo test --release (851 passed, with #352 and #353's tests), cargo clippy --all-targets and cargo fmt --check are clean on macOS. I did not run test/parity_suite.sh. FASTQ input without a blank line takes the old path after one peeked byte; I did not measure that byte's cost, and I expect none, since fill_buf is what the parser calls next anyway.

Second commit, after review

Review found that a FASTA record with no sequence (>e0, ACGT…, >e1, >e2, …) was read as an empty read: exit 0, e1 written with SEQ *. STAR stops with short read sequence line: 0 (readLoad.cpp:35-39, exit 104). cd5e39f stops both readers with the same message, an empty FASTQ sequence included, locked by a_read_with_no_bases_is_an_error. The exit code is 1, as for every error in rustar-aligner. The review confirmed the rest against STAR: mixed FASTA/FASTQ, CRLF FASTA and FASTQ, gzip FASTA with trailing blank lines, and S.i.1/S.i.2 names.

🤖 Generated with Claude Code

…ssion

FASTA input is accepted and written with `*` for quality. Trailing blank
lines end a file; data after one is an error. `.1`/`.2` after a read number
is a mate suffix, and a pair's QNAME is mate 1's name. Sort spills stay in a
prefix ending in `/`. `--outBAMcompression -1` compresses like level 1.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
A FASTA record with no sequence line (`>e1` followed by `>e2`) was read as
an empty read, written with SEQ `*`, and the run exited 0. STAR stops with
"short read sequence line: 0" (readLoad.cpp, exit 104). Both readers now
stop with the same message; an empty FASTQ sequence too.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant