faidx, locate, amplicon, fish: wrong output, panics and lost data, often with exit 0
Nobody has claimed this yet.
Assessment
- Difficulty
- 5/5
- Estimated time
- Over a week
- Newbie friendliness
- 20/100
- Issue type
- Bug
- Clarity
- Clearly specified
- Activity status
- Active
- Tech stack
- go
- Domain
- bioinformatics
Research direction
The report covers faidx, locate, amplicon, and fish across many independent failures. Start by reading the implementations and their tests for each command, then reproduce one listed example with seqkit v2.14.0 on Linux. Done means each reported case matches its expected behavior, including correct errors and exit statuses.
Written by the indexing model from the issue text.
Description
Version: seqkit v2.14.0 (bioconda), Linux x86_64
Each section below is one problem. Each example uses a small input file that
the section shows.
faidx: a negative begin before the sequence start gives an empty record
The help shows -12:-1 on a 10-base sequence as "the whole sequence". The
command writes an empty record.
$ printf '>s\nACGTACGTAC\n' > s.fa
$ seqkit faidx s.fa s:-12--1 s:-10--1
>s:-12--1
>s:-10--1
ACGTACGTAC
Expected: the begin is clipped to base 1, as the help example shows.
faidx: a failed index build leaves an empty index, and the next run exits 0
$ printf '>a\nACGT\nAC\nACGT\n' > rag.fa
$ seqkit faidx rag.fa a
[ERRO] different line length in sequence: a. Please format the file with 'seqkit seq'
$ seqkit faidx rag.fa a; echo $?
[INFO] read FASTA index from rag.fa.fai
[WARN] 0 records loaded from rag.fa.fai, please check if it matches the FASTA file, or switch on the flag -U/--update-faidx
0
The first run leaves a 0-byte rag.fa.fai. The second run writes nothing and
exits 0. An empty input file gives the same result.
Expected: no index file stays after a failed build.
faidx: a record whose name holds a colon cannot be got as a whole
$ printf '>chr1:100-200\nACGTACGT\n' > col.fa
$ seqkit faidx col.fa 'chr1:100-200'; echo $?
[WARN] sequence not found: chr1
0
Expected: the full text is tried as a record name first, as samtools faidx does.
faidx: the forms b--e and b- write a header with no region
$ seqkit faidx s.fa s:2--2 s:3-
>s
CGTACGTA
>s
GTACGTAC
Two different regions get the same name.
Expected: >s:2--2 and >s:3-, as for the other forms.
faidx: a region that is not found gives exit 0
$ seqkit faidx s.fa nope; echo $?
[WARN] sequence not found: nope
0
A script that checks only the exit code does not see the missing record.
Expected: a non-zero exit code after the other regions are written.
locate -c: a minus-strand hit across the end gets a negative start
$ printf '>seq\nACGTTGCA\n' > circ.fa
$ seqkit locate -p TTGCAACG -c circ.fa | tail -1
seq TTGCAACG TTGCAACG - 2 9 TTGCAACG
$ seqkit locate -p TTGCAACG -c -m 1 circ.fa | tail -1
seq TTGCAACG TTGCAACG - -6 1 TTGCAACG
$ seqkit locate -F -p TTGCAACG -c circ.fa | tail -1
seq TTGCAACG TTGCAACG - -6 1 TTGCAACG
$ seqkit locate -p TTGCAACG -c -m 1 --bed circ.fa | tail -1
seq -7 1 TTGCAACG 0 -
The three modes give two different positions for one hit.
Expected: 2 9 in every mode, as in the plain mode.
locate -f: two patterns with the same name, only the last is searched
$ printf '>m1\nGAATTC\n>m1\nAAGA\n' > patdup.fa
$ printf '>s1\nAAGAATTCTT\n' > s1.fa
$ seqkit locate -f patdup.fa s1.fa
seqID patternName pattern strand start end matched
s1 m1 AAGA + 1 4 AAGA
s1 m1 AAGA - 7 10 AAGA
GAATTC is in s1 but is not reported.
Expected: each pattern record is searched.
locate: -p is ignored when -f is given
$ printf '>m1\nGAATTC\n' > pat.fa
$ seqkit locate -f pat.fa -p AAGA s1.fa
The output has only GAATTC hits and no message.
Expected: an error that -p and -f cannot be given together, or a search for both.
locate: an empty pattern or an empty regex match gives end < start
$ printf '>e\n\n' > patempty.fa
$ seqkit locate -f patempty.fa circ.fa | head -3
seqID patternName pattern strand start end matched
seq e + 1 0
seq e + 2 1
$ seqkit locate -p 'C*' -r circ.fa | sed -n 2p
seq C* C* + 1 0
The minus strand also gets start 9 on a sequence of 8 bases.
Expected: an empty pattern is skipped with a warning; empty matches are not reported.
amplicon: no primer gives a panic
$ printf '>seq\nacgcccactgaaatga\n' > ex.fa
$ seqkit amplicon ex.fa; echo $?
[INFO] 1 primer pair loaded
panic: runtime error: invalid memory address or nil pointer dereference
2
Expected: an error that -F, -R or -p is needed.
amplicon: N in a primer matches only the letter N
The help says "degenerate bases allowed" for -F and -R.
$ seqkit amplicon -F ccc --bed ex.fa
seq 3 6 . 0 + ccc
$ seqkit amplicon -F ccN --bed ex.fa; echo $?
0
ccN must match ccc at the same place, but gives no result.
Expected: N matches any base, as the other degenerate codes do.
amplicon -m: a reverse site before the forward site gives an empty record
$ printf '>rev\naaagggccc\n' > rev.fa
$ seqkit amplicon -F ccc -R ttt -m 0 -M --bed rev.fa
$ seqkit amplicon -F ccc -R ttt -m 1 -M --bed rev.fa
rev 5 4 . 0 + 2 1 1
With -m 1 the record has start > end and no sequence. With -m 0 there is no result.
Expected: no result, as with -m 0.
amplicon -r: a region outside the amplicon gives one base
$ seqkit amplicon -F ccc -R ttt --bed ex.fa
seq 3 13 . 0 + cccactgaaa
$ seqkit amplicon -F ccc -R ttt -r 14:20 --bed ex.fa
seq 15 16 . 0 + a
$ seqkit amplicon -F ccc -R ttt -r -30:-20 --bed ex.fa
seq 0 1 . 0 + a
Expected: no output, as the flanking mode -f gives.
amplicon -p: bad primer-file lines are skipped with no message
$ printf 'p1\n' > p1.tsv
$ printf 'p1\tccc\tttt\textra\n' > p4.tsv
$ printf 'p1 ccc ttt\n' > psp.tsv
$ seqkit amplicon -p p1.tsv ex.fa; echo $?
[INFO] 0 primer pair loaded
0
The same for p4.tsv (4 columns) and psp.tsv (spaces, not tabs).
Expected: an error that names the file and the line.
fish: lower-case input gives Acc 0.00 for a perfect hit
$ printf '>lc\nttgcatgcaagtccgataggcttacagtcggatccaagtt\n' > lc.fa
$ seqkit fish -F GTCCGATAGGCTTAC lc.fa
Ref RefStart RefEnd Query QueryStart QueryEnd Strand MapQual RawScore Acc ClipAcc QueryCov
lc 10 25 q0 0 15 + 60.00 60 0.00 0.00 100.00
The upper-case reference gives 100.00 100.00. A lower-case query gives the same 0.00.
Expected: Acc and ClipAcc ignore case, as the score does.
fish -p, -r: bad values give a panic
$ printf '>ref1\nTTGCATGCAAGTCCGATAGGCTTACAGTCGGATCCAAGTT\n' > ref.fa
$ seqkit fish -p 4,-4,-2 -F GTCCGATAGGCTTAC ref.fa
panic: Invalid alignment parameters: 4,-4,-2
$ seqkit fish -r 5-30 -F GTCCGATAGGCTTAC ref.fa
panic: invalid range: 5-30
$ seqkit fish -r 30:5 -F GTCCGATAGGCTTAC ref.fa
panic: invalid range: 30:5
Each exits 2 with a Go stack trace.
Expected: the same message as an [ERRO] line, with no stack trace.
fish -p: positive mismatch and gap values give MapQual NaN
$ seqkit fish -p 4,4,2,1 -q 50 -F GTCCGATTGGCTTAC ref.fa
Ref RefStart RefEnd Query QueryStart QueryEnd Strand MapQual RawScore Acc ClipAcc QueryCov
ref1 0 28 q0 0 15 + NaN 99 2.44 2.44 100.00
ref1 0 28 q0 0 15 - NaN 99 2.44 2.44 100.00
NaN passes -q 50. With -b, the BAM record gets MAPQ 60.
Expected: an error for a mismatch or gap value above 0.
fish -b: the BAM has a PG:i:0 tag and PP:-
$ seqkit fish -b f.bam -F GTCCGATAGGCTTAC ref.fa
$ samtools view -h f.bam
[W::sam_hdr_link_pg] PG line with ID:seqkit has a PP link to missing program '-'
@SQ SN:ref1 LN:40
@PG ID:seqkit PN:seqkit CL:seqkit fish PP:- VN:1.0
q0 0 ref1 11 60 15M * 0 0 GTCCGATAGGCTTAC * PG:i:0 NM:i:0
The SAM specification gives PG the type Z (the program ID). No program - exists.
Expected: PG:Z:seqkit, and no PP field.
fish -b -D: @SQ SN holds the full header with spaces
$ printf '>ref1 first ref\nTTGCATGCAAGTCCGATAGGCTTACAGTCGGATCCAAGTT\n' > refd.fa
$ seqkit fish -D -b d.bam -F GTCCGATAGGCTTAC refd.fa
$ samtools view -H d.bam | head -1
@SQ SN:ref1 first ref LN:40
Expected: SN:ref1, the sequence ID.
- Dominant language
- Go
- Stars
- 1.6k
- Forks
- 189
- PR merge metrics
- No merged PRs in 30d
Getting set up
- Ships a Dockerfile or Docker Compose file
- No pull request template
- No contributing guide
First steps
- Read the whole issue, then the project's contributing guide.
- Comment on the issue to say you are picking it up — it saves two people doing the same work.
- Fork the repository and make your change on a branch.
- Open a pull request that references the issue number.
More from shenwei356/seqkit
-
convert: no range check between Phred encodings, so a wrong --from writes a broken FASTQ (exit 0)Open
Difficulty 2/5 1-3 hours Newbie friendliness 82/100
shenwei356/seqkit#612 · 1 comment ·
-
Difficulty 5/5 Over a week Newbie friendliness 35/100
shenwei356/seqkit#620 ·
-
sliding, merge-slides, mutate, translate, replace: panics, lost windows and wrong output with exit 0Open
Difficulty 5/5 Over a week Newbie friendliness 35/100
shenwei356/seqkit#619 ·
-
Difficulty 5/5 Over a week Newbie friendliness 25/100
shenwei356/seqkit#617 ·
-
Difficulty 5/5 Over a week Newbie friendliness 35/100
shenwei356/seqkit#616 ·
All issues in shenwei356/seqkit
Similar issues
-
cleanup hacktoberfest priority: low skill: good first issue testing
Difficulty 2/5 1-3 hours Newbie friendliness 78/100
hyperledger-labs/fabric-smart-client#2197 ·
Maintainers usually reply within 1 day
-
New term request
Difficulty 2/5 1-3 hours Newbie friendliness 62/100
geneontology/go-ontology#32791 ·
Maintainers usually reply within 1 day
-
Difficulty 2/5 1-3 hours Newbie friendliness 69/100
aws/amazon-ecs-agent#5177 ·
Maintainers usually reply within 2 days
-
release-followup
Difficulty 2/5 1-3 hours Newbie friendliness 82/100
Maintainers usually reply within 1 day
-
kind/bug
Difficulty 2/5 1-3 hours Newbie friendliness 78/100
Maintainers usually reply within 1 day