Repository navigation
common, pair, split, split2: lost or changed data with exit 0, and inputs deleted by -O / -f #620
Description
Version: seqkit v2.14.0 (bioconda), Linux x86_64
All examples below were run with seqkit v2.14.0. Each one exits 0 unless the
text says otherwise. The worst items come first: pair -O and -f can delete
the input files. The [INFO] lines are removed where they do not help.
Test files used in more than one section. Each section starts in a new, empty
folder that holds a copy of these files.
$ awk 'BEGIN{split("ACGTACGTAA GGGGCCCCAT TTTTAAAACC ACACACACAC GTGTGTGTGT CCCCAAAATT AATTCCGGAA TTAACCGGTT CAGTCAGTCA GACTGACTGA", q, " "); for (i = 1; i <= 10; i++) printf ">s%d desc %d\n%s\n", i, i, q[i]}' > a.fa
$ gzip -k a.fa
$ awk -v s=ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGT 'BEGIN{q=s; gsub(/./,"I",q); for(i=1;i<=60000;i++) printf "@r%d\n%s\n+\n%s\n",i,s,q}' > big1.fq
$ awk -v s=TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT 'BEGIN{q=s; gsub(/./,"I",q); for(i=1;i<=60000;i++) printf "@r%d\n%s\n+\n%s\n",i,s,q}' > big2.fq
a.fa has 10 records, s1 to s10. big1.fq and big2.fq have 60 000
reads each, with the same IDs and different sequences.
pair: -O that names the input folder in another way overwrites the inputs
When -O is the input folder written as an absolute path (or a symbolic link
to it, or read2's folder), pair does not see that it is the input folder.
The outputs then open the input files for writing.
$ wc -c big1.fq big2.fq
7908894 big1.fq
7908894 big2.fq
$ seqkit pair -1 big1.fq -2 big2.fq -O "$PWD"
[WARN] outdir not empty: /.../p1, you can use --force to overwrite
[ERRO] big1.fq: fastx: unequal sequence and quality
$ wc -c big1.fq big2.fq
0 big1.fq
0 big2.fq
Exit code 255, but both inputs are now empty.
Expected: an error before anything is written when an output path is an
input file.
pair, split, split2: -f deletes inputs and other files
-f deletes the whole output folder before the input is read.
$ mkdir od && cp a.fa od/ && seqkit split od/a.fa -s 4 -O od -f
[ERRO] fastx: stat od/a.fa: no such file or directory # input deleted, exit 255
$ mkdir -p x/y && touch x/keep.txt && cd x/y
$ seqkit split2 -s 4 -O .. -f ../../a.fa
[ERRO] unlinkat ./..: directory not empty
$ cd ../.. && ls -A x # empty: keep.txt and y are gone
$ mkdir o && touch o/keep.txt && : > empty.fa
$ seqkit split2 -s 4 -O o -f empty.fa; echo $?
0
$ ls o # empty: keep.txt is gone
Expected: refuse a folder that holds an input or the current folder, and
delete nothing when the input is empty or cannot be read.
pair, split2: two outputs of one run get the same path
$ mkdir A B && cp big1.fq A/reads.fq && cp big2.fq B/reads.fq
$ seqkit pair -1 A/reads.fq -2 B/reads.fq -O out
[INFO] 60000 paired-end reads saved to out/reads.fq and out/reads.fq
$ ls out
reads.fq
$ grep -c '^ACGT' out/reads.fq; grep -c '^TTTT' out/reads.fq
337
59663 # one file, a mix of both inputs
$ printf '@a\nACGT\n+\nIIII\n@b\nGGCC\n+\nIIII\n' > pr1.fq
$ printf '@a\nTTTT\n+\nIIII\n@b\nCCCC\n+\nIIII\n' > pr2.fq
$ seqkit split2 -1 pr1.fq -2 pr2.fq -s 2 -O o -o xxx
[INFO] write 2 sequences to file: o/xxx.part_001.fq
[INFO] write 2 sequences to file: o/xxx.part_001.fq
$ ls o
xxx.part_001.fq # holds only the reads of pr2.fq
$ printf '>dup AAAA\nAAAA\n>other\nCCCC\n>dup BBBB\nGGGG\n' > ids.fa
$ seqkit split2 -s 1 -N -O n ids.fa
[INFO] write 1 sequences to file: n/dup.fa
[INFO] write 1 sequences to file: n/other.fa
[INFO] write 1 sequences to file: n/dup.fa
$ cat n/dup.fa
>dup BBBB
GGGG # record "dup AAAA" is lost
Expected: an error when two outputs of one run have the same path.
pair, split, split2: "you can use --force to overwrite", but the run overwrites
Without -f, the warning says the folder is not empty, then the run replaces
files with the same name. Old parts of an earlier run stay next to the new
ones.
$ seqkit split2 -s 2 -O o a.fa # 10 records: 5 parts
$ seqkit split2 -s 5 -O o a.fa
[WARN] outdir not empty: o, you can use --force to overwrite
[INFO] write 5 sequences to file: o/a.part_001.fa
[INFO] write 5 sequences to file: o/a.part_002.fa
$ cat o/*.fa | grep -c '>'
16 # parts 3-5 of the first run are still there
Expected: without -f, do not replace a file that existed before the run.
pair: FASTA input gives a file that seqkit cannot read
$ printf '>a\nACGT\n' > r1.fa
$ printf '>a\nTTTT\n' > r2.fa
$ seqkit pair -1 r1.fa -2 r2.fa -O po
$ od -c po/r1.fa
0000000 @ a \n A C G T \n + \n \n
$ seqkit seq po/r1.fa
[ERRO] fastx: unequal sequence and quality
Expected: FASTA input gives FASTA output.
pair: a second record with the same ID replaces the first
$ printf '@a\nAAAA\n+\nIIII\n@a\nCCCC\n+\nIIII\n' > d1.fq
$ printf '@z\nGGGG\n+\nIIII\n' > d2.fq
$ seqkit pair -1 d1.fq -2 d2.fq -O u -u
[INFO] 0 paired-end reads saved to u/d1.fq and u/d2.fq
[INFO] 1 unpaired reads saved to u/d1.unpaired.fq
[INFO] 1 unpaired reads saved to u/d2.unpaired.fq
$ cat u/d1.unpaired.fq
@a
CCCC
+
IIII # record "a AAAA" is lost, no message
Expected: no record is lost.
common -s -e: the result changes from run to run
$ printf '>q1 d\nTTGGGGTT\n>q2\nACGT\n' > e1.fa
$ printf '>v1\nCCCC\n>v2\nTTGG\n' > e2.fa
$ for i in $(seq 12); do seqkit common -s -e e1.fa e2.fa 2>/dev/null | tr '\n' ' '; echo; done
>q1:3-6 d GGGG >q1:1-4 d TTGG
>q1:1-4 d TTGG >q1:3-6 d GGGG
>q1:1-4 d TTGG >q1:3-6 d GGGG
>q1:3-6 d GGGG >q1:1-4 d TTGG
>q1:1-4 d TTGG >q1:3-6 d GGGG
In this example, 5 of 12 runs write two pieces, and 7 runs write nothing.
The order of the pieces also changes from run to run.
Expected: the same output in each run.
common -s -e: FASTQ pieces keep the full sequence and quality
$ printf '@q1 d\nTTGGGGTT\n+\nABCDEFGH\n' > q.fq; printf '>v1\nGGGG\n' > g.fa
$ seqkit common -s -e q.fq g.fa 2>/dev/null
@q1:3-6 d
TTGGGGTT
+
ABCDEFGH
Expected: GGGG and CDEF, as with FASTA input.
common -s -e: empty sequences
An empty sequence in the first file panics:
$ printf '>x\n\n>s1\nACGTACGTAA\n' > es1.fa; printf '>y\n\n>s1\nAAA\n' > es2.fa
$ seqkit common -s -e es1.fa es2.fa
panic: runtime error: slice bounds out of range [-1:]
$ echo $?
2
An empty sequence in a later file matches every record:
$ printf '>s1\nACGTACGTAA\n>s2\nGGGGCCCCAT\n' > two.fa
$ printf '>e\n\n' > es.fa
$ seqkit common -s -e two.fa es.fa 2>/dev/null
>s1:1-0
ACGTACGTAA
>s2:1-0
GGGGCCCCAT
Expected: an empty sequence is skipped.
common -s -e: the counts do not match the output
$ printf '@s1 first desc\nACGTACGTAA\n+\nIIIIIIIIII\n@s3\nTTTTAAAACC\n+\n##########\n@s7\nCCCC\n+\nIIII\n' > c1.fq
$ printf '@s3 z\nTTTTAAAACC\n+\nJJJJJJJJJJ\n@s1\nACGTACGTAA\n+\nKKKKKKKKKK\n' > c2.fq
$ seqkit common -s -e c1.fq c2.fq > /dev/null
[INFO] read file 1/2: c1.fq
[INFO] 3 seqs loaded
[INFO] read file 2/2: c2.fq
[INFO] 1 seqs left
[INFO] 2 unique sequences found in 2 files, which belong to 2 records in the first file: c1.fq
[INFO] 2 common/shared sequences saved to: -
"1 seqs left", then 2 unique sequences and 2 records written.
Expected: counts that match the records written.
common: "-" as the first file writes records from the second file
$ printf '>s1 first desc\nACGTACGTAA\n>s2 second\nGGGGCCCCAT\n>s3\nTTTTAAAACC\n' > sa.fa
$ printf '>s3 other\nTTTTAAAACC\n>s1 first desc\nAAAAAAAA\n>s2 second\nGGGGCCCCAT\n' > sb.fa
$ seqkit common sa.fa sb.fa 2>/dev/null | grep -c '>'
3
$ cat sa.fa | seqkit common - sb.fa
[INFO] read file 1/2: -
[INFO] read file 2/2: sb.fa
[INFO] find common seqs ...
[INFO] 3 unique sequence IDs found in 2 files, which belong to 3 records in the first file: sb.fa
[INFO] 2 common/shared sequences saved to: -
>s1 first desc
AAAAAAAA
>s2 second
GGGGCCCCAT
s1 gets the sequence of sb.fa, s3 is lost (its header differs), and the
message names sb.fa as the first file.
Expected: the same output as seqkit common sa.fa sb.fa.
split -2: an index made by seqkit faidx drops the descriptions
$ seqkit split a.fa -s 4 -2 -O o1
$ head -1 o1/a.part_001.fasta
>s1 desc 1
$ seqkit faidx a.fa && mv a.fa.fai a.fa.seqkit.fai
$ seqkit split a.fa -s 4 -2 -O o2
[INFO] read FASTA index from a.fa.seqkit.fai
$ head -1 o2/a.part_001.fasta
>s1 # input header is ">s1 desc 1"
Expected: the descriptions are kept, as with the index that split -2 makes.
split -2: two records with the same header line
$ printf '>a x\nAC\n>b\nAG\n>a x\nTT\n' > d2.fa
$ seqkit split d2.fa -s 1 -2 -O s
[fai warning] ignoring duplicate sequence "a x" at byte offset 19
$ ls s
d2.part_001.fasta d2.part_002.fasta # record 3 is lost
$ seqkit split d2.fa -r 1:1 -2 -O r
$ cat r/d2.part_1:1_T.fasta
>a x
AC # should be TT
Expected: the same output as without -2.
split -2: names, compression, messages and temporary files
$ seqkit split a.fa.gz -s 3 -O one # one/a.part_001.fa.gz ...
$ seqkit split a.fa.gz -s 3 -2 -O two # two/a.part_001.fasta ... (not gzip)
$ seqkit split a.fa.gz -r 1:1 -2 -O tr
[INFO] write 3 sequences to file: tr/a.part_1:1_T.fasta
[INFO] write 2 sequences to file: tr/a.part_1:1_T.fasta
[INFO] write 2 sequences to file: tr/a.part_1:1_G.fasta
[INFO] write 3 sequences to file: tr/a.part_1:1_G.fasta
$ ls tr
a.part_1:1_A.fasta a.part_1:1_C.fasta a.part_1:1_G.fasta a.part_1:1_T.fasta
$ ls
a.fa.gz a.fa.gz.fastx a.fa.gz.fastx.seqkit.fai ... # temporary files stay
With -2, part names get .fasta and lose the compression. With -r, the
messages name the wrong files (T and G twice each, A and C never), and the temporary files stay (without -k). With -s, the
split into line is missing.
Expected: the same names, compression and messages as without -2, and no
temporary files left.
split -r: an end before base 1 gives the whole sequence
$ printf '>x\nACGTA\n' > five.fa
$ seqkit split five.fa -r -7:-6 # writes five.part_-7:-6_ACGTA.fa, exit 0
$ seqkit split five.fa -r 1:-6 # writes five.part_1:-6_ACGTA.fa, exit 0
$ seqkit split five.fa -r 2:-5 # writes five.part_2:-5_.fa, exit 0
$ seqkit split five.fa -r 4:2
[ERRO] region (4:2) not match sequence (x) with length of 5
Expected: the same region error for the first three.
split -p: fewer parts than asked
$ seqkit split a.fa -p 6 -O p6 # 10 records
$ ls p6
a.part_001.fa a.part_002.fa a.part_003.fa a.part_004.fa a.part_005.fa
Expected: 6 parts.
split2 -p with -s or -l: the last part is empty and its records are lost
$ seqkit split2 -s 5 -p 2 -O sp a.fa # 10 records
[INFO] split into 5 seqs per file
[INFO] write 5 sequences to file: sp/a.part_001.fa
$ grep -c '>' sp/*
sp/a.part_001.fa:5
sp/a.part_002.fa:0 # s6-s10 are lost
$ seqkit split2 -p 2 -l 1000 -O pl a.fa
[INFO] split into 2 parts
$ grep -c '>' pl/*
0 # all records are lost
Expected: all 10 records are written.
split2 -1 -2 -l: mates go to different parts
Each file is cut by its own lengths.
$ printf '@a\nACGTACGTACGT\n+\nIIIIIIIIIIII\n@b\nAC\n+\nII\n@c\nAC\n+\nII\n' > l1.fq
$ printf '@a\nAC\n+\nII\n@b\nAC\n+\nII\n@c\nACGTACGTACGT\n+\nIIIIIIIIIIII\n' > l2.fq
$ seqkit split2 -1 l1.fq -2 l2.fq -l 12 -O o
o/l1.part_001.fq: @a
o/l1.part_002.fq: @b @c
o/l2.part_001.fq: @a @b @c
Expected: read2 follows its mate in read1.
split2 -1 -2: no check of the record counts
$ printf '@a\nAC\n+\nII\n@b\nAC\n+\nII\n@c\nAC\n+\nII\n' > m1.fq
$ printf '@b\nAC\n+\nII\n@c\nAC\n+\nII\n' > m2.fq
$ seqkit split2 -1 m1.fq -2 m2.fq -s 2 -O pm; echo $?
0
pm/m1.part_001.fq: @a @b
pm/m1.part_002.fq: @c
pm/m2.part_001.fq: @b @c
Expected: an error when read1 and read2 have different numbers of records.
split2 -l: wrong message for a bad value
$ seqkit split2 -l -5 a.fa
[ERRO] one of flags should be given: -s/-p/-l. type "seqkit split2 -h" for help
$ seqkit split2 -l 0.5 a.fa # same message
$ seqkit split2 -l 99999999999999999999 -O o a.fa
[INFO] split by sequence length: 99999999999999999999
[ERRO] one of flags should be given: -s/-p/-l. ...
$ ls -d o # an empty folder stays
o
Expected: a message that names the bad value of -l.
split2 -P: the prefix can write outside the output folder
$ mkdir w && cd w
$ seqkit split2 -s 5 -O out -P ../x_ ../a.fa
[INFO] write 5 sequences to file: x_001.fa
[INFO] write 5 sequences to file: x_002.fa
$ ls
out x_001.fa x_002.fa # out/ is empty
Expected: an error, or the parts inside out/.
Activity