diff --git a/Makefile b/Makefile index 467aec9..2a7817b 100644 --- a/Makefile +++ b/Makefile @@ -156,8 +156,8 @@ bump-version: jjnew: @echo "$(YELLOW)→ Creating a new commit...$(NC)" - @echo "$(BLUE)→ Documenting current commit...$(NC)" - @jj auto-describe + @echo "$(BLUE)→ Documenting undocumented commits...$(NC)" + @jj auto-doc @echo "$(BLUE)→ Done.$(NC)" @jj new @echo "$(GREEN)✓ New commit created$(NC)" @@ -171,8 +171,8 @@ jjpush: @echo "$(GREEN)✓ Release complete$(NC)" jjpush-describe: - @echo "$(BLUE)→ Documenting current commit...$(NC)" - @jj auto-describe + @echo "$(BLUE)→ Documenting undocumented commits...$(NC)" + @jj auto-doc jjpush-bump: @echo "$(BLUE)→ Creating new commit for version bump...$(NC)" diff --git a/pkg/obialign/locatepattern.go b/pkg/obialign/locatepattern.go index d5170e6..304cfbe 100644 --- a/pkg/obialign/locatepattern.go +++ b/pkg/obialign/locatepattern.go @@ -28,10 +28,18 @@ func buffIndex(i, j, width int) int { // // The function returns the start and end positions of the best // match, as well as the number of errors in the best match. +// +// When the sequence is too short relative to the pattern for the +// backtracking to reconstruct a valid alignment (e.g. the pattern +// is longer than the sequence, or the match sits too close to a +// sequence boundary), no reliable position can be computed. In that +// case the function returns the sentinel (-1, -1, -1) instead of a +// guessed, potentially wrong, position: callers must treat this as +// "no match" rather than use the returned coordinates. func LocatePattern(id string, pattern, sequence []byte) (int, int, int) { - if len(pattern) >= len(sequence) { - log.Panicf("Sequence %s:Pattern %s must be shorter than sequence %s", id, pattern, sequence) + if len(sequence) == 0 { + log.Panicf("Sequence %s:Pattern %s must not be empty", id, pattern) } // Pattern spreads over the columns @@ -158,5 +166,18 @@ func LocatePattern(id string, pattern, sequence []byte) (int, int, int) { // obilog.Warnf("from : %d to: %d error: %d match: %v", // i, end+1, -buffer[buffIndex(len(sequence)-1, len(pattern)-1, width)], // string(sequence[i:(end+1)])) + + if i < 0 || end == -1 { + // i < 0: the backtracking ran off the start of the sequence + // without fully consuming the pattern. + // end == -1: the backtracking loop never ran at all (e.g. a + // single-base pattern, jmax == 0), so no alignment boundary + // was ever established. + // Either way, no valid alignment exists for this (pattern, + // sequence) pair: signal it explicitly instead of returning + // an out-of-bounds or uncomputed position. + return -1, -1, -1 + } + return i, end + 1, -buffer[buffIndex(len(sequence)-1, len(pattern)-1, width)] } diff --git a/pkg/obialign/locatepattern_test.go b/pkg/obialign/locatepattern_test.go new file mode 100644 index 0000000..70fcbfe --- /dev/null +++ b/pkg/obialign/locatepattern_test.go @@ -0,0 +1,123 @@ +package obialign + +import ( + "math/rand" + "testing" +) + +func TestLocatePatternNormal(t *testing.T) { + // Pattern fully and exactly present in the middle of a longer sequence. + start, end, nerr := LocatePattern("id", []byte("ACGT"), []byte("TTTTACGTTTTT")) + if start != 4 || end != 8 || nerr != 0 { + t.Errorf("got start=%d end=%d nerr=%d, want start=4 end=8 nerr=0", start, end, nerr) + } +} + +func TestLocatePatternOneMismatch(t *testing.T) { + start, end, nerr := LocatePattern("id", []byte("ACGT"), []byte("TTTTACTTTTTT")) + if nerr != 1 { + t.Errorf("got nerr=%d, want 1 (start=%d end=%d)", nerr, start, end) + } +} + +// The real-world case that used to panic: pattern longer than the sequence +// fragment extracted for indel relocation. +func TestLocatePatternPatternLongerThanSequence(t *testing.T) { + start, end, nerr := LocatePattern("id", + []byte("GGGCAATCCTGAGCCAAATC"), + []byte("tcctgagccaaatcacgtt")) + + if start != -1 || end != -1 || nerr != -1 { + t.Errorf("got start=%d end=%d nerr=%d, want the (-1,-1,-1) sentinel", start, end, nerr) + } +} + +func TestLocatePatternSequenceLengthOne(t *testing.T) { + start, end, nerr := LocatePattern("id", []byte("AB"), []byte("A")) + + if start < 0 || end < 0 || nerr < 0 { + t.Fatalf("got start=%d end=%d nerr=%d, expected a valid (non-sentinel) result", start, end, nerr) + } + if start != 0 || end != 1 || nerr != 1 { + t.Errorf("got start=%d end=%d nerr=%d, want start=0 end=1 nerr=1", start, end, nerr) + } +} + +// A pattern much longer than the sequence can still yield a mathematically +// valid (in-bounds) alignment: the extra pattern length is absorbed as gaps, +// driving the error count high enough that the caller's maxerr threshold +// rejects it. The function itself must still return consistent bounds. +func TestLocatePatternPatternMuchLongerThanSequence(t *testing.T) { + start, end, nerr := LocatePattern("id", []byte("ACGTACGTACGTACGTACGT"), []byte("ACG")) + + isSentinel := start == -1 && end == -1 && nerr == -1 + isValid := start >= 0 && end > start && end <= 3 && nerr >= 0 + + if !isSentinel && !isValid { + t.Errorf("got start=%d end=%d nerr=%d, want either the sentinel or consistent in-bounds values", start, end, nerr) + } +} + +func TestLocatePatternNeverReturnsOutOfBounds(t *testing.T) { + patterns := []string{"A", "AC", "ACG", "ACGT", "ACGTA", "ACGTAC", "ACGTACG", "ACGTACGT"} + sequences := []string{"A", "AC", "ACG", "ACGT", "ACGTA", "ACGTAC", "ACGTACG", "ACGTACGT"} + + for _, p := range patterns { + for _, s := range sequences { + start, end, nerr := LocatePattern("id", []byte(p), []byte(s)) + + if start == -1 && end == -1 && nerr == -1 { + // Explicit "no reliable match" sentinel: always acceptable. + continue + } + + if start < 0 || end < 0 || start >= end || end > len(s) || nerr < 0 { + t.Errorf("pattern=%q sequence=%q -> start=%d end=%d nerr=%d is out of bounds / inconsistent", + p, s, start, end, nerr) + } + } + } +} + +// Randomized property test over a wide range of pattern/sequence length +// combinations, including pattern >= sequence, to make sure the function +// never panics and never returns anything but the sentinel or fully +// consistent, in-bounds coordinates. +func TestLocatePatternRandomizedNeverInvalid(t *testing.T) { + const bases = "ACGT" + rng := rand.New(rand.NewSource(42)) + + randSeq := func(n int) []byte { + b := make([]byte, n) + for i := range b { + b[i] = bases[rng.Intn(len(bases))] + } + return b + } + + for trial := 0; trial < 5000; trial++ { + patLen := 1 + rng.Intn(15) + seqLen := 1 + rng.Intn(15) + + pattern := randSeq(patLen) + sequence := randSeq(seqLen) + + func() { + defer func() { + if r := recover(); r != nil { + t.Fatalf("panic for pattern=%q sequence=%q: %v", pattern, sequence, r) + } + }() + + start, end, nerr := LocatePattern("id", pattern, sequence) + + isSentinel := start == -1 && end == -1 && nerr == -1 + isValid := start >= 0 && end > start && end <= len(sequence) && nerr >= 0 + + if !isSentinel && !isValid { + t.Errorf("pattern=%q sequence=%q -> start=%d end=%d nerr=%d is neither the sentinel nor consistent", + pattern, sequence, start, end, nerr) + } + }() + } +} diff --git a/pkg/obiapat/pattern.go b/pkg/obiapat/pattern.go index 714c978..b58473f 100755 --- a/pkg/obiapat/pattern.go +++ b/pkg/obiapat/pattern.go @@ -373,6 +373,7 @@ func (pattern ApatPattern) BestMatch(sequence ApatSequence, begin, length int) ( cpattern := (*[1 << 30]byte)(unsafe.Pointer(pattern.pointer.pointer.cpat)) frg := sequence.pointer.reference.Sequence()[start:end] + fragStart := start log.Debugln( string(frg), @@ -384,11 +385,20 @@ func (pattern ApatPattern) BestMatch(sequence ApatSequence, begin, length int) ( (*cpattern)[0:int(pattern.pointer.pointer.patlen)], frg) - // olderr := m[2] + if from < 0 { + // obialign.LocatePattern could not reconstruct a reliable + // alignment (e.g. the fragment is too short relative to the + // pattern). Reporting a guessed position would risk placing + // the primer boundary incorrectly, so treat it as no match + // at all rather than falling back to an unrefined position. + matched = false + log.Debugln("No reliable indel relocation, discarding match", sequence.pointer.reference.Id()) + return + } nerr = score - start = start + from - end = start + to + start = fragStart + from + end = fragStart + to log.Debugf("BestMatch on %s : score=%d [%d..%d]", sequence.pointer.reference.Id(), score, start, nerr) return } @@ -467,6 +477,7 @@ func (pattern ApatPattern) AllMatches(sequence ApatSequence, begin, length int) for _, m := range res { // Recompute the start and end position of the match // when the pattern allows for indels + valid := true if m[2] > 0 && pattern.pointer.pointer.hasIndel { // obilog.Warnf("Locating indel on sequence %s[%s]", sequence.pointer.reference.Id(), pattern.String()) start := m[0] - m[2]*2 @@ -485,16 +496,26 @@ func (pattern ApatPattern) AllMatches(sequence ApatSequence, begin, length int) (*cpattern)[0:int(pattern.pointer.pointer.patlen)], frg) - // olderr := m[2] - m[2] = score - m[0] = start + pb - m[1] = start + pe + if pb < 0 { + // obialign.LocatePattern could not reconstruct a + // reliable alignment (e.g. the match sits too close + // to a sequence end for the fragment to be usable). + // Reporting a guessed position risks placing the + // primer boundary incorrectly, so drop the match + // entirely instead of keeping an unrefined guess. + valid = false + } else { + // olderr := m[2] + m[2] = score + m[0] = start + pb + m[1] = start + pe - // obilog.Warnf("seq[%d@%d:%d] %d: %s %d - %s:%s:%s", i, m[0], m[1], olderr, sequence.pointer.reference.Id(), score, - // frg, (*cpattern)[0:int(pattern.pointer.pointer.patlen)], sequence.pointer.reference.Sequence()[m[0]:m[1]]) + // obilog.Warnf("seq[%d@%d:%d] %d: %s %d - %s:%s:%s", i, m[0], m[1], olderr, sequence.pointer.reference.Id(), score, + // frg, (*cpattern)[0:int(pattern.pointer.pointer.patlen)], sequence.pointer.reference.Sequence()[m[0]:m[1]]) + } } - if int(pattern.pointer.pointer.maxerr) >= m[2] { + if valid && int(pattern.pointer.pointer.maxerr) >= m[2] { res[j] = m j++ } diff --git a/pkg/obiformats/fastqseq_read.go b/pkg/obiformats/fastqseq_read.go index 861705f..125c838 100644 --- a/pkg/obiformats/fastqseq_read.go +++ b/pkg/obiformats/fastqseq_read.go @@ -131,7 +131,7 @@ func _storeSequenceQuality(bytes *bytes.Buffer, out *obiseq.BioSequence, quality out.SetQualities(q) } -func FastqChunkParser(quality_shift byte, with_quality bool, UtoT bool) func(string, io.Reader) (obiseq.BioSequenceSlice, error) { +func FastqChunkParser(quality_shift byte, with_quality bool, UtoT bool, fileName string) func(string, io.Reader) (obiseq.BioSequenceSlice, error) { parser := func(source string, input io.Reader) (obiseq.BioSequenceSlice, error) { var identifier string @@ -160,12 +160,12 @@ func FastqChunkParser(quality_shift byte, with_quality bool, UtoT bool) func(str // Beginning of sequence state = 1 } else { - log.Fatalf("%s : sequence entry is not starting with @", source) + log.Fatalf("file %s: sequence entry is not starting with @", fileName) } case 1: // Beginning of identifier (Mandatory) if is_sep { // No identifier -> ERROR - log.Fatalf("%s : sequence identifier is empty", source) + log.Fatalf("file %s: sequence identifier is empty", fileName) } else { // Beginning of identifier state = 2 @@ -221,7 +221,7 @@ func FastqChunkParser(quality_shift byte, with_quality bool, UtoT bool) func(str // End of sequence rawseq := seqBytes.Bytes() if len(rawseq) == 0 { - log.Fatalf("@%s[%s] : sequence is empty", identifier, source) + log.Fatalf("file %s: record @%s has an empty sequence line", fileName, identifier) } s := obiseq.NewBioSequence(identifier, rawseq, definition) s.SetSource(source) @@ -241,8 +241,8 @@ func FastqChunkParser(quality_shift byte, with_quality bool, UtoT bool) func(str context = append( append([]byte{previous}, C), context...) - log.Fatalf("%s [%s]: sequence contains invalid character %c (%s)", - source, identifier, C, string(context)) + log.Fatalf("file %s: record @%s contains invalid character %c (%s)", + fileName, identifier, C, string(context)) } } case 7: @@ -251,7 +251,7 @@ func FastqChunkParser(quality_shift byte, with_quality bool, UtoT bool) func(str } else if C == '+' { state = 8 } else { - log.Fatalf("@%s[%s] : sequence data not followed by a line starting with + but a %c", identifier, source, C) + log.Fatalf("file %s: record @%s: sequence data not followed by a line starting with + but a %c", fileName, identifier, C) } case 8: // State consuming the + internal header line @@ -282,7 +282,7 @@ func FastqChunkParser(quality_shift byte, with_quality bool, UtoT bool) func(str } else if C == '@' { state = 1 } else { - log.Fatalf("%s[%s] : sequence record not followed by a line starting with @", identifier, source) + log.Fatalf("file %s: record @%s not followed by a line starting with @", fileName, identifier) } } @@ -304,7 +304,7 @@ func FastqChunkParser(quality_shift byte, with_quality bool, UtoT bool) func(str } // FastqChunkParserRope parses a FASTQ chunk directly from a rope without Pack(). -func FastqChunkParserRope(source string, rope *PieceOfChunk, quality_shift byte, with_quality, UtoT bool) (obiseq.BioSequenceSlice, error) { +func FastqChunkParserRope(source string, rope *PieceOfChunk, quality_shift byte, with_quality, UtoT bool, fileName string) (obiseq.BioSequenceSlice, error) { scanner := newRopeScanner(rope) sequences := obiseq.MakeBioSequenceSlice(100)[:0] @@ -334,7 +334,7 @@ func FastqChunkParserRope(source string, rope *PieceOfChunk, quality_shift byte, // Line 2: sequence sline := scanner.ReadLine() if sline == nil { - log.Fatalf("@%s[%s]: unexpected EOF after header", id, source) + log.Fatalf("file %s: record @%s is truncated (header line with no sequence line following) — the FASTQ file appears incomplete", fileName, id) } seqDest := make([]byte, len(sline)) w := 0 @@ -350,7 +350,7 @@ func FastqChunkParserRope(source string, rope *PieceOfChunk, quality_shift byte, } seqDest = seqDest[:w] if len(seqDest) == 0 { - log.Fatalf("@%s[%s]: sequence is empty", id, source) + log.Fatalf("file %s: record @%s has an empty sequence line", fileName, id) } // Line 3: + (skip) @@ -382,16 +382,17 @@ func _ParseFastqFile( out obiiter.IBioSequence, quality_shift byte, with_quality, UtoT bool, + fileName string, ) { - parser := FastqChunkParser(quality_shift, with_quality, UtoT) + parser := FastqChunkParser(quality_shift, with_quality, UtoT, fileName) for chunks := range input { var sequences obiseq.BioSequenceSlice var err error if chunks.Rope != nil { - sequences, err = FastqChunkParserRope(chunks.Source, chunks.Rope, quality_shift, with_quality, UtoT) + sequences, err = FastqChunkParserRope(chunks.Source, chunks.Rope, quality_shift, with_quality, UtoT, fileName) } else { sequences, err = parser(chunks.Source, chunks.Raw) } @@ -423,6 +424,8 @@ func ReadFastq(reader io.Reader, options ...WithOption) (obiiter.IBioSequence, e false, ) + fileName := opt.FileName() + for i := 0; i < nworker; i++ { out.Add(1) go _ParseFastqFile( @@ -431,6 +434,7 @@ func ReadFastq(reader io.Reader, options ...WithOption) (obiiter.IBioSequence, e obidefault.ReadQualitiesShift(), opt.ReadQualities(), opt.UtoT(), + fileName, ) } @@ -456,7 +460,9 @@ func ReadFastq(reader io.Reader, options ...WithOption) (obiiter.IBioSequence, e } func ReadFastqFromFile(filename string, options ...WithOption) (obiiter.IBioSequence, error) { - options = append(options, OptionsSource(obiutils.RemoveAllExt((path.Base(filename))))) + options = append(options, + OptionsSource(obiutils.RemoveAllExt((path.Base(filename)))), + OptionsFileName(filename)) file, err := obiutils.Ropen(filename) diff --git a/pkg/obiformats/fastseq_json_header.go b/pkg/obiformats/fastseq_json_header.go index af8de1c..3928613 100644 --- a/pkg/obiformats/fastseq_json_header.go +++ b/pkg/obiformats/fastseq_json_header.go @@ -199,8 +199,111 @@ func _parse_json_array_interface(str []byte) ([]interface{}, error) { return values, nil } -func _parse_json_header_(header string, sequence *obiseq.BioSequence) string { +// _parse_json_annotation_field parses a single key/value pair coming from a +// JSON object (either a FASTA/FASTQ inline JSON header, or the "annotations" +// field of a JSON sequence record) and applies it to the sequence, special +// casing the well-known OBITools attributes (id, definition, count, taxid, +// obiclean_*, merged_*). +func _parse_json_annotation_field(key []byte, value []byte, dataType jsonparser.ValueType, sequence *obiseq.BioSequence) error { annotations := sequence.Annotations() + var err error + + skey := obiutils.UnsafeString(key) + + switch { + case skey == "id": + sequence.SetId(string(value)) + case skey == "definition": + sequence.SetDefinition(string(value)) + + case skey == "count": + if dataType != jsonparser.Number { + log.Fatalf("%s: Count attribut must be numeric: %s", sequence.Id(), string(value)) + } + count, err := jsonparser.ParseInt(value) + if err != nil { + log.Fatalf("%s: Cannot parse count %s", sequence.Id(), string(value)) + } + sequence.SetCount(int(count)) + + case skey == "obiclean_weight": + weight, err := _parse_json_map_int(value) + if err != nil { + log.Fatalf("%s: Cannot parse obiclean weight %s", sequence.Id(), string(value)) + } + annotations[skey] = weight + + case skey == "obiclean_status": + status, err := _parse_json_map_string(value) + if err != nil { + log.Fatalf("%s: Cannot parse obiclean status %s", sequence.Id(), string(value)) + } + annotations[skey] = status + + case strings.HasPrefix(skey, "merged_"): + if dataType == jsonparser.Object { + data, err := _parse_json_map_int(value) + if err != nil { + log.Fatalf("%s: Cannot parse merged slot %s: %v", sequence.Id(), skey, err) + } else { + annotations[skey] = obiseq.MapAsStatsOnValues(data) + } + } else { + log.Fatalf("%s: Cannot parse merged slot %s", sequence.Id(), skey) + } + + case skey == "taxid": + if dataType == jsonparser.Number || dataType == jsonparser.String { + taxid := string(value) + sequence.SetTaxid(taxid) + } else { + log.Fatalf("%s: Cannot parse taxid %s", sequence.Id(), string(value)) + } + + case strings.HasSuffix(skey, "_taxid"): + if dataType == jsonparser.Number || dataType == jsonparser.String { + rank := skey[:len(skey)-len("_taxid")] + + taxid := string(value) + sequence.SetTaxid(taxid, rank) + } else { + log.Fatalf("%s: Cannot parse taxid %s", sequence.Id(), string(value)) + } + + default: + skey = strings.Clone(skey) + switch dataType { + case jsonparser.String: + annotations[skey] = string(value) + case jsonparser.Number: + // Try to parse the number as an int at first then as float if that fails. + annotations[skey], err = jsonparser.ParseInt(value) + if err != nil { + annotations[skey], err = strconv.ParseFloat(obiutils.UnsafeString(value), 64) + } + case jsonparser.Array: + annotations[skey], err = _parse_json_array_interface(value) + case jsonparser.Object: + annotations[skey], err = _parse_json_map_interface(value) + case jsonparser.Boolean: + annotations[skey], err = jsonparser.ParseBoolean(value) + case jsonparser.Null: + annotations[skey] = nil + default: + log.Fatalf("Unknown data type %v", dataType) + } + } + + if err != nil { + annotations[skey] = "NaN" + log.Fatalf("%s: Cannot parse value %s assicated to key %s into a %s value", + sequence.Id(), string(value), skey, dataType.String()) + } + + return err +} + +func _parse_json_header_(header string, sequence *obiseq.BioSequence) string { start := -1 stop := -1 level := 0 @@ -240,101 +343,7 @@ func _parse_json_header_(header string, sequence *obiseq.BioSequence) string { jsonparser.ObjectEach(obiutils.UnsafeBytes(header[start:stop]), func(key []byte, value []byte, dataType jsonparser.ValueType, offset int) error { - var err error - - skey := obiutils.UnsafeString(key) - - switch { - case skey == "id": - sequence.SetId(string(value)) - case skey == "definition": - sequence.SetDefinition(string(value)) - - case skey == "count": - if dataType != jsonparser.Number { - log.Fatalf("%s: Count attribut must be numeric: %s", sequence.Id(), string(value)) - } - count, err := jsonparser.ParseInt(value) - if err != nil { - log.Fatalf("%s: Cannot parse count %s", sequence.Id(), string(value)) - } - sequence.SetCount(int(count)) - - case skey == "obiclean_weight": - weight, err := _parse_json_map_int(value) - if err != nil { - log.Fatalf("%s: Cannot parse obiclean weight %s", sequence.Id(), string(value)) - } - annotations[skey] = weight - - case skey == "obiclean_status": - status, err := _parse_json_map_string(value) - if err != nil { - log.Fatalf("%s: Cannot parse obiclean status %s", sequence.Id(), string(value)) - } - annotations[skey] = status - - case strings.HasPrefix(skey, "merged_"): - if dataType == jsonparser.Object { - data, err := _parse_json_map_int(value) - if err != nil { - log.Fatalf("%s: Cannot parse merged slot %s: %v", sequence.Id(), skey, err) - } else { - annotations[skey] = obiseq.MapAsStatsOnValues(data) - } - } else { - log.Fatalf("%s: Cannot parse merged slot %s", sequence.Id(), skey) - } - - case skey == "taxid": - if dataType == jsonparser.Number || dataType == jsonparser.String { - taxid := string(value) - sequence.SetTaxid(taxid) - } else { - log.Fatalf("%s: Cannot parse taxid %s", sequence.Id(), string(value)) - } - - case strings.HasSuffix(skey, "_taxid"): - if dataType == jsonparser.Number || dataType == jsonparser.String { - rank := skey[:len(skey)-len("_taxid")] - - taxid := string(value) - sequence.SetTaxid(taxid, rank) - } else { - log.Fatalf("%s: Cannot parse taxid %s", sequence.Id(), string(value)) - } - - default: - skey = strings.Clone(skey) - switch dataType { - case jsonparser.String: - annotations[skey] = string(value) - case jsonparser.Number: - // Try to parse the number as an int at first then as float if that fails. - annotations[skey], err = jsonparser.ParseInt(value) - if err != nil { - annotations[skey], err = strconv.ParseFloat(obiutils.UnsafeString(value), 64) - } - case jsonparser.Array: - annotations[skey], err = _parse_json_array_interface(value) - case jsonparser.Object: - annotations[skey], err = _parse_json_map_interface(value) - case jsonparser.Boolean: - annotations[skey], err = jsonparser.ParseBoolean(value) - case jsonparser.Null: - annotations[skey] = nil - default: - log.Fatalf("Unknown data type %v", dataType) - } - } - - if err != nil { - annotations[skey] = "NaN" - log.Fatalf("%s: Cannot parse value %s assicated to key %s into a %s value", - sequence.Id(), string(value), skey, dataType.String()) - } - - return err + return _parse_json_annotation_field(key, value, dataType, sequence) }, ) diff --git a/pkg/obiformats/json_read.go b/pkg/obiformats/json_read.go new file mode 100644 index 0000000..7658ff1 --- /dev/null +++ b/pkg/obiformats/json_read.go @@ -0,0 +1,147 @@ +package obiformats + +import ( + "io" + "os" + "path" + + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obidefault" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiiter" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq" + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiutils" + "github.com/buger/jsonparser" + "github.com/goccy/go-json" + log "github.com/sirupsen/logrus" +) + +// _parse_json_record parses a single JSON object describing a sequence +// (as produced by JSONRecord in json_writer.go) into a *obiseq.BioSequence. +func _parse_json_record(raw []byte, shift byte) *obiseq.BioSequence { + sequence := obiseq.NewEmptyBioSequence(0) + + if id, err := jsonparser.GetString(raw, "id"); err == nil { + sequence.SetId(id) + } + + if seq, err := jsonparser.GetString(raw, "sequence"); err == nil { + sequence.SetSequence([]byte(seq)) + } + + if qual, err := jsonparser.GetString(raw, "qualities"); err == nil { + q := []byte(qual) + for i := 0; i < len(q); i++ { + q[i] -= shift + } + sequence.SetQualities(q) + } + + if annot, dataType, _, err := jsonparser.Get(raw, "annotations"); err == nil && dataType == jsonparser.Object { + jsonparser.ObjectEach(annot, + func(key []byte, value []byte, valType jsonparser.ValueType, offset int) error { + return _parse_json_annotation_field(key, value, valType, sequence) + }, + ) + } + + return sequence +} + +// _ParseJsonFile streams the top-level JSON array, decoding and pushing one +// batch of sequences at a time, without ever loading the whole document in +// memory. Only one raw record at a time is buffered by the decoder. +func _ParseJsonFile(source string, + reader io.Reader, + out obiiter.IBioSequence, + shift byte, + batchSize int) { + + dec := json.NewDecoder(reader) + + if _, err := dec.Token(); err != nil { + if err == io.EOF { + out.Done() + return + } + log.Fatalf("cannot parse JSON data: %v", err) + } + + slice := obiseq.MakeBioSequenceSlice() + o := 0 + + for dec.More() { + var raw json.RawMessage + + if err := dec.Decode(&raw); err != nil { + log.Fatalf("cannot parse JSON data: %v", err) + } + + sequence := _parse_json_record(raw, shift) + + slice = append(slice, sequence) + if len(slice) >= batchSize { + out.Push(obiiter.MakeBioSequenceBatch(source, o, slice)) + o++ + slice = obiseq.MakeBioSequenceSlice() + } + } + + if len(slice) > 0 { + out.Push(obiiter.MakeBioSequenceBatch(source, o, slice)) + } + + out.Done() +} + +func ReadJSON(reader io.Reader, options ...WithOption) (obiiter.IBioSequence, error) { + + opt := MakeOptions(options) + out := obiiter.MakeIBioSequence() + + out.Add(1) + go _ParseJsonFile(opt.Source(), + reader, + out, + obidefault.ReadQualitiesShift(), + opt.BatchSize()) + + go func() { + out.WaitAndClose() + }() + + return out, nil + +} + +func ReadJSONFromFile(filename string, options ...WithOption) (obiiter.IBioSequence, error) { + + options = append(options, OptionsSource(obiutils.RemoveAllExt((path.Base(filename))))) + file, err := obiutils.Ropen(filename) + + if err == obiutils.ErrNoContent { + log.Infof("file %s is empty", filename) + return ReadEmptyFile(options...) + } + + if err != nil { + return obiiter.NilIBioSequence, err + } + + return ReadJSON(file, options...) +} + +func ReadJSONFromStdin(reader io.Reader, options ...WithOption) (obiiter.IBioSequence, error) { + options = append(options, OptionsSource(obiutils.RemoveAllExt("stdin"))) + input, err := obiutils.Buf(os.Stdin) + + if err == obiutils.ErrNoContent { + log.Infof("stdin is empty") + return ReadEmptyFile(options...) + } + + if err != nil { + log.Fatalf("open file error: %v", err) + return obiiter.NilIBioSequence, err + } + + return ReadJSON(input, options...) +} diff --git a/pkg/obiformats/options.go b/pkg/obiformats/options.go index 8190801..e0c4525 100644 --- a/pkg/obiformats/options.go +++ b/pkg/obiformats/options.go @@ -35,6 +35,7 @@ type __options__ struct { csv_auto bool paired_filename string source string + filename string with_feature_table bool with_pattern bool with_parent bool @@ -216,6 +217,16 @@ func (opt Options) Source() string { return opt.pointer.source } +// FileName returns the full path of the file being read, for use in +// diagnostic messages. It falls back to Source() when no explicit +// file name has been set (e.g. reading from stdin or a raw reader). +func (opt Options) FileName() string { + if opt.pointer.filename == "" { + return opt.pointer.source + } + return opt.pointer.filename +} + func (opt Options) WithFeatureTable() bool { return opt.pointer.with_feature_table } @@ -421,6 +432,14 @@ func OptionsSource(source string) WithOption { return f } +func OptionsFileName(filename string) WithOption { + f := WithOption(func(opt Options) { + opt.pointer.filename = filename + }) + + return f +} + func OptionsWithProgressBar() WithOption { f := WithOption(func(opt Options) { opt.pointer.with_progress_bar = true diff --git a/pkg/obiformats/universal_read.go b/pkg/obiformats/universal_read.go index 62e967e..7df141c 100644 --- a/pkg/obiformats/universal_read.go +++ b/pkg/obiformats/universal_read.go @@ -145,6 +145,8 @@ func ReadSequencesFromFile(filename string, return ReadGenbank(reader, options...) case "text/csv": return ReadCSV(reader, options...) + case "application/json": + return ReadJSON(reader, options...) default: log.Fatalf("File %s has guessed format %s which is not yet implemented", filename, mime.String()) diff --git a/pkg/obifp/uint128_test.go b/pkg/obifp/uint128_test.go index bc834ee..411416d 100644 --- a/pkg/obifp/uint128_test.go +++ b/pkg/obifp/uint128_test.go @@ -134,7 +134,7 @@ func TestUint128_QuoRem(t *testing.T) { u := Uint128{w1: 3, w0: 8} v := Uint128{w1: 0, w0: 4} q, r := u.QuoRem(v) - assert.Equal(t, Uint128{w1: 0, w0: 2}, q) + assert.Equal(t, Uint128{w1: 0, w0: 13835058055282163714}, q) assert.Equal(t, Uint128{w1: 0, w0: 0}, r) } @@ -150,7 +150,7 @@ func TestUint128_Div(t *testing.T) { u := Uint128{w1: 3, w0: 8} v := Uint128{w1: 0, w0: 4} q := u.Div(v) - assert.Equal(t, Uint128{w1: 0, w0: 2}, q) + assert.Equal(t, Uint128{w1: 0, w0: 13835058055282163714}, q) } func TestUint128_Div64(t *testing.T) { @@ -183,7 +183,7 @@ func TestUint128_Cmp(t *testing.T) { func TestUint128_Cmp64(t *testing.T) { u := Uint128{w1: 1, w0: 2} v := uint64(3) - assert.Equal(t, -1, u.Cmp64(v)) + assert.Equal(t, 1, u.Cmp64(v)) } func TestUint128_Equals(t *testing.T) { diff --git a/pkg/obingslibrary/multimatch.go b/pkg/obingslibrary/multimatch.go index b28befa..3d32942 100644 --- a/pkg/obingslibrary/multimatch.go +++ b/pkg/obingslibrary/multimatch.go @@ -777,7 +777,7 @@ func (library *NGSLibrary) ExtractMultiBarcodeSliceWorker(options ...WithOption) library.SetAllowsIndels(true) } - if opt.AllowedMismatches() > 0 { + if opt.AllowedMismatchesIsSet() { library.SetAllowedMismatches(opt.AllowedMismatches()) } diff --git a/pkg/obingslibrary/worker.go b/pkg/obingslibrary/worker.go index 2cfe391..75347f6 100644 --- a/pkg/obingslibrary/worker.go +++ b/pkg/obingslibrary/worker.go @@ -6,13 +6,14 @@ import ( ) type _Options struct { - discardErrors bool - unidentified string - allowedMismatch int - allowsIndel bool - withProgressBar bool - parallelWorkers int - batchSize int + discardErrors bool + unidentified string + allowedMismatch int + allowedMismatchSet bool + allowsIndel bool + withProgressBar bool + parallelWorkers int + batchSize int } // Options stores a set of option usable by the @@ -52,6 +53,7 @@ func OptionWithProgressBar(yes bool) WithOption { func OptionAllowedMismatches(count int) WithOption { f := WithOption(func(opt Options) { opt.pointer.allowedMismatch = count + opt.pointer.allowedMismatchSet = true }) return f @@ -97,6 +99,12 @@ func (options Options) AllowedMismatches() int { return options.pointer.allowedMismatch } +// AllowedMismatchesIsSet returns true if OptionAllowedMismatches +// was explicitly applied to these options. +func (options Options) AllowedMismatchesIsSet() bool { + return options.pointer.allowedMismatchSet +} + func (options Options) AllowsIndels() bool { return options.pointer.allowsIndel } diff --git a/pkg/obioptions/version.go b/pkg/obioptions/version.go index 3229975..ac10bdb 100644 --- a/pkg/obioptions/version.go +++ b/pkg/obioptions/version.go @@ -3,7 +3,7 @@ package obioptions // Version is automatically updated by the Makefile from version.txt // The patch number (third digit) is incremented on each push to the repository -var _Version = "Release 4.4.46" +var _Version = "Release 4.5.0" // Version returns the version of the obitools package. // diff --git a/pkg/obitools/obiconvert/options.go b/pkg/obitools/obiconvert/options.go index a860194..a7fbe02 100644 --- a/pkg/obitools/obiconvert/options.go +++ b/pkg/obitools/obiconvert/options.go @@ -24,6 +24,7 @@ var __input_genbank_format__ = false var __input_fastq_format__ = false var __input_fasta_format__ = false var __input_csv_format__ = false +var __input_json_format__ = false var __output_in_fasta__ = false var __output_in_fastq__ = false @@ -71,6 +72,9 @@ func InputOptionSet(options *getoptions.GetOpt) { options.BoolVar(&__input_csv_format__, "csv", __input_csv_format__, options.Description("Read data following the CSV format.")) + options.BoolVar(&__input_json_format__, "json", __input_json_format__, + options.Description("Read data following the JSON format.")) + options.BoolVar(&__no_ordered_input__, "no-order", __no_ordered_input__, options.Description("When several input files are provided, "+ "indicates that there is no order among them.")) @@ -158,6 +162,8 @@ func CLIInputFormat() string { return "genbank" case __input_csv_format__: return "csv" + case __input_json_format__: + return "json" default: return "guessed" } diff --git a/pkg/obitools/obiconvert/sequence_reader.go b/pkg/obitools/obiconvert/sequence_reader.go index cbda802..d16f39e 100644 --- a/pkg/obitools/obiconvert/sequence_reader.go +++ b/pkg/obitools/obiconvert/sequence_reader.go @@ -73,7 +73,9 @@ func ExpandListOfFiles(check_ext bool, filenames ...string) ([]string, error) { strings.HasSuffix(path, "dat") || strings.HasSuffix(path, "dat.gz") || strings.HasSuffix(path, "ecopcr") || - strings.HasSuffix(path, "ecopcr.gz") { + strings.HasSuffix(path, "ecopcr.gz") || + strings.HasSuffix(path, "json") || + strings.HasSuffix(path, "json.gz") { log.Debugf("Appending %s file\n", path) list_of_files.Add(path) } @@ -142,6 +144,8 @@ func CLIReadBioSequences(filenames ...string) (obiiter.IBioSequence, error) { iterator, err = obiformats.ReadFastq(os.Stdin, opts...) case "csv": iterator, err = obiformats.ReadCSV(os.Stdin, opts...) + case "json": + iterator, err = obiformats.ReadJSON(os.Stdin, opts...) default: iterator, err = obiformats.ReadSequencesFromStdin(opts...) } @@ -163,6 +167,8 @@ func CLIReadBioSequences(filenames ...string) (obiiter.IBioSequence, error) { reader = obiformats.ReadFastaFromFile case "csv": reader = obiformats.ReadCSVFromFile + case "json": + reader = obiformats.ReadJSONFromFile case "ecopcr": reader = obiformats.ReadEcoPCRFromFile case "embl": diff --git a/pkg/obitools/obilandmark/obilandmark.go b/pkg/obitools/obilandmark/obilandmark.go index 175ee0e..20e6db9 100644 --- a/pkg/obitools/obilandmark/obilandmark.go +++ b/pkg/obitools/obilandmark/obilandmark.go @@ -170,7 +170,7 @@ func CLISelectLandmarkSequences(iterator obiiter.IBioSequence) obiiter.IBioSeque for i, seq := range library { taxon := seq.Taxon(taxo) if taxon == nil { - log.Fatal("%s: Cannot identify taxid %s in %s", seq.Id(), seq.Taxid(), taxo.Name()) + log.Fatalf("%s: Cannot identify taxid %s in %s", seq.Id(), seq.Taxid(), taxo.Name()) } taxa.Set(i, taxon) } diff --git a/pkg/obitools/obimultiplex/demultiplex.go b/pkg/obitools/obimultiplex/demultiplex.go index 064e4d5..996843b 100644 --- a/pkg/obitools/obimultiplex/demultiplex.go +++ b/pkg/obitools/obimultiplex/demultiplex.go @@ -15,7 +15,6 @@ func IExtractBarcode(iterator obiiter.IBioSequence) (obiiter.IBioSequence, error opts := make([]obingslibrary.WithOption, 0, 10) opts = append(opts, - obingslibrary.OptionAllowedMismatches(CLIAllowedMismatch()), obingslibrary.OptionAllowedIndel(CLIAllowsIndel()), obingslibrary.OptionUnidentified(CLIUnidentifiedFileName()), obingslibrary.OptionDiscardErrors(!CLIConservedErrors()), @@ -23,6 +22,14 @@ func IExtractBarcode(iterator obiiter.IBioSequence) (obiiter.IBioSequence, error obingslibrary.OptionBatchSize(obidefault.BatchSize()), ) + // Only propagate the CLI --allowed-mismatches value if the user + // explicitly set it: otherwise the per-primer values defined in + // the NGSFilter config file (@primer_mismatches, @forward_mismatches, + // @reverse_mismatches) must be preserved. + if CLIAllowedMismatchIsSet() { + opts = append(opts, obingslibrary.OptionAllowedMismatches(CLIAllowedMismatch())) + } + ngsfilter, err := CLINGSFIlter() if err != nil { log.Fatalf("%v", err) diff --git a/pkg/obitools/obimultiplex/options.go b/pkg/obitools/obimultiplex/options.go index 3809d76..7fd4897 100644 --- a/pkg/obitools/obimultiplex/options.go +++ b/pkg/obitools/obimultiplex/options.go @@ -18,6 +18,7 @@ var _UnidentifiedFile = "" var _AllowedMismatch = 2 var _AllowsIndel = false var _ConservedError = false +var _optionsParser *getoptions.GetOpt // PCROptionSet defines every options related to a simulated PCR. // @@ -29,6 +30,8 @@ var _ConservedError = false // - option : is a pointer to a getoptions.GetOpt instance normaly // produced by the func MultiplexOptionSet(options *getoptions.GetOpt) { + _optionsParser = options + options.StringVar(&_NGSFilterFile, "tag-list", _NGSFilterFile, options.Alias("s"), options.Description("File name of the NGSFilter file describing PCRs.")) @@ -62,6 +65,15 @@ func CLIAllowedMismatch() int { return _AllowedMismatch } +// CLIAllowedMismatchIsSet returns true if the user explicitly +// specified --allowed-mismatches on the command line, as opposed +// to relying on its default value. This allows per-primer mismatch +// settings from the NGSFilter config file to take precedence unless +// the user explicitly overrides them from the CLI. +func CLIAllowedMismatchIsSet() bool { + return _optionsParser != nil && _optionsParser.Called("allowed-mismatches") +} + func CLIAllowsIndel() bool { return _AllowsIndel } diff --git a/pkg/obitools/obitagpcr/pcrtag.go b/pkg/obitools/obitagpcr/pcrtag.go index bd32402..a6c7f31 100644 --- a/pkg/obitools/obitagpcr/pcrtag.go +++ b/pkg/obitools/obitagpcr/pcrtag.go @@ -55,7 +55,7 @@ func IPCRTagPESequencesBatch(iterator obiiter.IBioSequence, ngsfilter.SetAllowsIndels(true) } - if obimultiplex.CLIAllowedMismatch() > 0 { + if obimultiplex.CLIAllowedMismatchIsSet() { ngsfilter.SetAllowedMismatches(obimultiplex.CLIAllowedMismatch()) } diff --git a/pkg/obiutils/mimetypes.go b/pkg/obiutils/mimetypes.go index 52cf083..901dea2 100644 --- a/pkg/obiutils/mimetypes.go +++ b/pkg/obiutils/mimetypes.go @@ -102,6 +102,11 @@ func RegisterOBIMimeType() { return ok } + jsonDetector := func(raw []byte, limit uint32) bool { + raw = bytes.TrimLeft(raw, " \t\r\n") + return len(raw) > 0 && (raw[0] == '[' || raw[0] == '{') + } + mimetype.Lookup("text/plain").Extend(fastaDetector, "text/fasta", ".fasta") mimetype.Lookup("text/plain").Extend(fastqDetector, "text/fastq", ".fastq") mimetype.Lookup("text/plain").Extend(ecoPCR2Detector, "text/ecopcr2", ".ecopcr") @@ -115,6 +120,7 @@ func RegisterOBIMimeType() { mimetype.Lookup("application/octet-stream").Extend(genbankDetector, "text/genbank", ".seq") mimetype.Lookup("application/octet-stream").Extend(emblDetector, "text/embl", ".dat") mimetype.Lookup("application/octet-stream").Extend(csv, "text/csv", ".csv") + mimetype.Lookup("application/octet-stream").Extend(jsonDetector, "application/json", ".json") } __obimimetype_registred__ = true } diff --git a/version.txt b/version.txt index 9e1350a..a84947d 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -4.4.46 +4.5.0