nx_fasta.nx source
↩ module page · 164 lines · 6934 B
1// nx_fasta.nx -- FASTA single-record stream parser + strip-whitespace.
2//
3// license_tier: INDEPENDENT_REDERIVE
4// genealogy_id: international-research-sources/lipman-1985-fasta-format
5//
6// G2.0c of NISHI_GENOMICS_SUBSTRATE_ROADMAP.md. Reference-genome
7// ingestion bridge: most reference assemblies (GRCh38, T2T-CHM13,
8// every UCSC + Ensembl + NCBI download) ship as gzipped FASTA where
9// each chromosome / contig is a record with multi-line wrapped
10// sequence (canonical wrap = 60 or 80 chars per line).
11//
12// FASTA canonical format (Lipman + Pearson 1985):
13// line 1: '>' + identifier [+ optional description]
14// line 2..N: nucleotide (or amino-acid) sequence, line-wrapped
15// record ends at next '>' or EOF
16//
17// Differs from FASTQ:
18// - '>' marker (not '@')
19// - Sequence WRAPS across multiple lines -- any line not starting
20// with '>' is sequence content (no '+' / qual lines)
21// - No quality scores
22//
23// API:
24// fasta_parse_one_record(buf, buf_len, offset, out_meta) -> i64
25// out_meta layout (5 i64s):
26// [0] header_start (offset after the '>', start of identifier)
27// [1] header_len (bytes from header_start to end of line 1)
28// [2] seq_start (first byte of seq region)
29// [3] seq_end_excl (one past last byte of seq region;
30// == offset of next '>' or buf_len at EOF)
31// [4] next_offset (where the next record begins; same as seq_end_excl)
32// Return: 0 success / 1 EOF / -1 malformed
33//
34// fasta_strip_whitespace(src, src_start, src_end_excl,
35// dst, max_dst) -> i64
36// Copies src[src_start..src_end_excl) into dst, skipping
37// whitespace bytes (LF 0x0A, CR 0x0D, space 0x20, tab 0x09).
38// Returns count written or -1 if max_dst exceeded.
39//
40// Typical caller flow:
41// fasta_parse_one_record(buf, len, off, meta)
42// n = fasta_strip_whitespace(buf, meta[2], meta[3], clean_buf, cap)
43// dna_pack_ambig(clean_buf, n, out_bases, out_nbits)
44//
45// What G2.0c does NOT do (deferred):
46// - gzip / bgzip decompression -- callers pre-decompress (G2.0c.2)
47// - .fai (FASTA index) random-access lookup -- G2.0c.3
48// - Multi-FASTA mmap-streaming for chromosome-scale references -- G2.0c.4
49// - Lowercase-as-soft-masked detection -- G0.5c
50// - Protein-FASTA support (different alphabet) -- nx_protein (G7)
51//
52// nx_safety_envelope: (schema: nishi-library/seeds/safety-critical-standards.toml)
53// intended_use: "FASTA reference-genome + multi-contig
54// sequence ingestion; the canonical entry
55// point for reference data into the
56// alignment + variant-calling pipeline"
57// sil_target: SIL2
58// asil_target: QM
59// dal_target: DAL C
60// iec_62304_class: B
61// evidence: [no_floating_point, deterministic,
62// bit_equal_reproducible,
63// lipman_pearson_1985_canonical_format,
64// composes_nx_const_ASCII_markers,
65// multi_line_seq_KAT,
66// strip_whitespace_KAT_with_lf_cr_space_tab,
67// license_tier_INDEPENDENT_REDERIVE]
68// hazard_register: [bug-tape-fasta-greater-than-sign-inside-header-misread,
69// bug-tape-fasta-line-ending-cr-only-old-mac,
70// bug-tape-fasta-empty-record-zero-len-seq,
71// bug-tape-fasta-strip-skips-real-base-char,
72// bug-tape-fasta-seq-end-includes-next-record-bytes]
73// residual_risk: "G2.0c assumes Unix LF line endings; CRLF
74// is handled because CR + LF both get
75// stripped by fasta_strip_whitespace, but
76// old-Mac CR-only is not officially supported.
77// Headers containing literal '>' inside the
78// description (rare but legal in some
79// BLAST conventions) read correctly because
80// only line-start '>' terminates a record."
81// verdict: NOT_YET_EVALUATED
82
83import "nx_syscalls.nx"
84import "nx_const.nx"
85import "nx_fastq.nx" // for nx_find_lf
86
87// Whitespace test: returns 1 if byte is LF / CR / space / tab; else 0.
88func nx_is_whitespace(b: i64) -> i64 {
89 let c: i64 = b & 0xff
90 if c == NX_ASCII_LF { return 1 }
91 if c == NX_ASCII_CR { return 1 }
92 if c == NX_ASCII_SPACE { return 1 }
93 if c == NX_ASCII_TAB { return 1 }
94 return 0
95}
96
97// Find the next line that BEGINS with '>' in [start, end).
98// Returns offset of the '>' or -1 if none. Lines are LF-separated.
99// Used to find the END of a FASTA record's seq region.
100func fasta_find_next_record_start(buf: *u8, start: i64, end: i64) -> i64 {
101 // Walk LF-by-LF; after each LF check if the NEXT byte is '>'.
102 var i: i64 = start
103 while i < end {
104 if (buf[i] & 0xff) == NX_ASCII_LF {
105 let next: i64 = i + 1
106 if next < end {
107 if (buf[next] & 0xff) == NX_ASCII_GT { return next }
108 }
109 }
110 i = i + 1
111 }
112 return -1
113}
114
115// Copy src[src_start..src_end_excl) into dst, skipping whitespace.
116// Returns bytes written or -1 if dst capacity exceeded.
117func fasta_strip_whitespace(src: *u8, src_start: i64, src_end_excl: i64,
118 dst: *u8, max_dst: i64) -> i64 {
119 if src_start < 0 { return -1 }
120 if src_end_excl < src_start { return -1 }
121 if max_dst < 0 { return -1 }
122
123 var written: i64 = 0
124 var i: i64 = src_start
125 while i < src_end_excl {
126 let b: i64 = src[i] & 0xff
127 if nx_is_whitespace(b) == 0 {
128 if written >= max_dst { return -1 }
129 dst[written] = b & 0xff
130 written = written + 1
131 }
132 i = i + 1
133 }
134 return written
135}
136
137// Parse one FASTA record beginning at `offset`. See header for
138// metadata layout and return codes.
139func fasta_parse_one_record(buf: *u8, buf_len: i64, offset: i64,
140 out_meta: *i64) -> i64 {
141 if buf_len < 0 { return -1 }
142 if offset < 0 { return -1 }
143 if offset >= buf_len { return 1 }
144
145 // ---- Header line: '>' + identifier ----
146 if (buf[offset] & 0xff) != NX_ASCII_GT { return -1 }
147 let header_start: i64 = offset + 1
148 let lf1: i64 = nx_find_lf(buf, header_start, buf_len)
149 if lf1 < 0 { return -1 }
150 let header_len: i64 = lf1 - header_start
151
152 // ---- Sequence region: from after lf1 to next '>' or EOF ----
153 let seq_start: i64 = lf1 + 1
154 var seq_end_excl: i64 = buf_len
155 let next_rec: i64 = fasta_find_next_record_start(buf, seq_start, buf_len)
156 if next_rec >= 0 { seq_end_excl = next_rec }
157
158 out_meta[0] = header_start
159 out_meta[1] = header_len
160 out_meta[2] = seq_start
161 out_meta[3] = seq_end_excl
162 out_meta[4] = seq_end_excl
163 return 0
164}