Skip to content

sliding, merge-slides, mutate, translate, replace: panics, lost windows and wrong output with exit 0 #619

Open

Description

@splaisan

Version: seqkit v2.14.0 (bioconda), Linux x86_64

All examples below use small input files. Each command was run on v2.14.0.
od -c shows bytes that a terminal does not show.

sliding -c: a window longer than the sequence contains memory bytes

What happens

With -c (circular), a window that goes around the sequence more than once
is not filled with sequence letters. It contains NUL bytes or letters of an
earlier record. The name also shows an end past the sequence length
(1-15 on a 10 bp sequence). The command exits 0.

$ printf '>s\nACGTACGTAC\n' > s10.fa
$ seqkit sliding -c -W 25 -s 5 s10.fa | od -c | head -4
0000000   >   s   _   s   l   i   d   i   n   g   :   1   -   1   5  \n
0000020   A   C   G   T   A   C   G   T   A   C   A   C   G   T   A   C
0000040   G   T   A   C  \n  \0  \0  \0  \0  \n   >   s   _   s   l   i

With a first record >p AAAAAAAAAA before s, the same window of s
contains AAAA from record p.

Expected

The window takes its letters around the sequence as often as necessary
(25 letters), and the name shows an end in the range 1 to 10.

sliding -c: an empty sequence gives one record of NUL bytes

$ printf '>e\n\n' > empty.fa
$ seqkit sliding -c -W 5 -s 5 empty.fa | od -c
0000000   >   e   _   s   l   i   d   i   n   g   :   1   -   5  \n  \n
0000020  \0  \0  \0  \0  \n

Exit 0. Without -c, the same file gives no output.

Expected

No output, as without -c.

sliding -c: a very large window panics

$ seqkit sliding -c -W 99999999999 -s 1 s10.fa
panic: runtime error: slice bounds out of range [:99999999989] with capacity 10485758

Exit 2.

Expected

No panic. The window is filled as in the first section, or an error message
says that the window is too large.

merge-slides -r: an expression with fewer than 3 groups panics

$ printf 'ref_sliding:1-10\nref_sliding:2-5\nref_sliding:8-20\n' > m1.txt
$ seqkit merge-slides -r '^(.+)_sliding:(\d+)-' m1.txt
panic: runtime error: index out of range [3] with length 3

Exit 2.

Expected

An error message that says -r must contain 3 capture groups, before the
input is read.

merge-slides: a merged region can become shorter

The end of a merged region is the end of the last window that was merged,
not the largest end. A window inside the region makes the region shorter.

$ printf 'ref_sliding:1-10\nref_sliding:2-5\n' > m6.txt
$ seqkit merge-slides m6.txt
ref	0	5

With an unsorted file, 11-20, 1-10, 6-15 gives ref 10 15. Exit 0.

Expected

ref 0 10: a merged region keeps the smallest start and the largest end.

merge-slides: a window with start 0 is lost

$ printf 'ref_sliding:0-10\nref_sliding:20-30\n' > m3.txt
$ seqkit merge-slides m3.txt
ref	19	30

The window 0-10 is not in the output and no message is written. Exit 0.
A start that is not a number has the same result.

Expected

The window is kept, or an error message says that the start is not valid.

merge-slides: a file with no window writes a line

$ printf '# only comment\n\n' > m4.txt
$ seqkit merge-slides m4.txt | cat -A
^I-1^I0$

Exit 0.

Expected

No output.

merge-slides -b: a line longer than the buffer stops the reading with no message

The help of -b says: "You need increase the value when "bufio.Scanner:
token too long" error reported". This error is not reported. Reading stops
at the long line, the windows after it are lost, and the exit code is 0.

$ awk 'BEGIN{print "ref_sliding:1-10"; s=""; for(i=0;i<300;i++) s=s"x"; print s; print "ref_sliding:50-60"}' > m5.txt
$ seqkit merge-slides -b 100 m5.txt
ref	0	10

Expected

[ERRO] bufio.Scanner: token too long and a non-zero exit code.

mutate: a very large negative position panics

$ seqkit mutate -p -99999999999999999999:x s10.fa
panic: runtime error: index out of range [-9223372036854775798]

Exit 2. -i and -d have the same panic.

Expected

The "out of sequence length" warning, as for other positions out of range.

mutate: position -(L+1) edits base 1

On a 10 bp sequence, -11 is out of range, but it is used as base 1.
Position -12 gives the warning correctly.

$ seqkit mutate -p -11:x s10.fa
>s
xCGTACGTAC
$ seqkit mutate -i -11:GG s10.fa
>s
AGGCGTACGTAC
$ seqkit mutate -p -12:x s10.fa
[WARN] [s]: point mutation: position (-12) out of sequence length (10)

-d 1:-11 deletes nothing and writes no warning. Exit 0 in each case.

Expected

-11 gives the same warning as -12 and 11.

translate --trim: removes everything after the first X or *

The help says: "remove all 'X' and '*' characters from the right end of the
translation".

$ printf '>t\nATGNNNAAA\n' > t1.fa
$ seqkit translate t1.fa
>t
MXK
$ seqkit translate --trim t1.fa
>t
M

Expected

MXK: only X and * at the right end are removed.

translate -s: begin is one too small on frames -1, -2, -3

On the forward strand a piece of 2 amino acids covers 6 bases. On the
reverse strand it covers 7 bases.

$ printf '>u\nATGAAATAGATGCCCTAA\n' > t2.fa
$ seqkit translate -f 1 -s t2.fa
>u_frame=1_begin=1_end=6
MK
>u_frame=1_begin=10_end=15
MP
$ seqkit seq -r -p -t dna t2.fa > t2rc.fa
$ seqkit translate -f -1 -s t2rc.fa
>u_frame=-1_begin=12_end=18
MK
>u_frame=-1_begin=3_end=9
MP

Expected

begin=13_end=18 and begin=4_end=9.

translate -m: no effect without -s

$ seqkit translate -m 10 t2.fa
>u
MK*MP*

The protein has 6 amino acids, but it is written with -m 10. With -s,
-m works.

Expected

-m also removes short proteins without -s, or the help says that -m
needs -s.

translate -e: the skipped record loses the frame in its name

$ printf '>a first\nATGAA\n>b second\nAT\n>c third\nATGAAA\n' > t5.fa
$ seqkit translate -e -F -f 1,3 t5.fa
>a_frame=1 first
M
>a_frame=3 first
E
>b second

>b second

>c_frame=1 third
MK
>c_frame=3 third
E

Record b is written twice with the same name. The _frame= part of -F
is missing. With a record that has no description, the name ends with a
space (>b ).

Expected

>b_frame=1 second and >b_frame=3 second, each with an empty sequence.

translate -l, -L, -f: a value that is not valid writes nothing and exits 0

$ seqkit translate -l 99 t2.fa; echo $?
0
$ seqkit translate -L 99 t2.fa; echo $?
0
$ seqkit translate -f '' t2.fa; echo $?
0
$ seqkit translate -T 99 t2.fa; echo $?
[ERRO] invalid translate table: 99
255

Expected

The same error as -T 99 for -l 99 and -L 99, and an error for an empty
-f. Each with a non-zero exit code.

replace -s on FASTQ: a same-length replacement with a group is refused

$ printf '@r\nACGT\n+\nABCD\n' > q.fq
$ seqkit replace -s -p '(A)C' -r '${1}G' q.fq
[ERRO] quality length mismatch after replacement (seq: 4 bp, qual: 3). Pattern may be too complex for quality tracking.
$ seqkit replace -s -p '(A)(C)' -r '$1$1' q.fq
[ERRO] quality length mismatch after replacement (seq: 4 bp, qual: 3). Pattern may be too complex for quality tracking.

Exit 255. The new sequence has the same length as the old one.

Expected

The replacement is done and the quality line keeps 4 letters.

replace -s on FASTQ: the quality does not follow the bases

Numbered groups move the quality with the bases. Named groups do not:

$ seqkit replace -s -p '(A)(C)' -r '$2$1' q.fq
@r
CAGT
+
BACD
$ seqkit replace -s -p '(?P<a>A)(?P<c>C)' -r '${c}${a}' q.fq
@r
CAGT
+
ABCD

A replacement that is longer than its match takes the quality of the next
bases:

$ printf '@r\nACGTTT\n+\nABCDEF\n' > q2.fq
$ seqkit replace -s -p 'T+|A' -r 'GG' q2.fq
@r
GGCGGG
+
ABBCDE

The first GG gets AB, but B is the quality of C. Exit 0.

Expected

Named groups move the quality as numbered groups do. A new base gets a
quality from its own match, not from the next base.

replace {nr}: the number starts again in each input file

$ printf '>x\nA\n>y\nC\n' > f1.fa
$ printf '>z\nG\n>w\nT\n>v\nA\n' > f2.fa
$ seqkit replace -p '.+' -r 'r{nr}' f1.fa f2.fa | grep '>'
>r1
>r2
>r1
>r2
>r3

The names repeat. The help says "Record number, starting from 1".

Expected

One counter over all input files: r1 to r5. If the restart is intended,
the help could say "record number in each file".

Activity

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

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions