lh3/seqtk

Toolkit for processing sequences in FASTA/Q formats

Introduction

Seqtk is a fast and lightweight tool for processing sequences in the FASTA or FASTQ format. It seamlessly parses both FASTA and FASTQ files which can also be optionally compressed by gzip. To install seqtk,

git clone https://github.com/lh3/seqtk.git;
cd seqtk; make

The only library dependency is zlib.

Seqtk Examples

• Convert FASTQ to FASTA:

    seqtk seq -a in.fq.gz > out.fa
  

• Convert ILLUMINA 1.3+ FASTQ to FASTA and mask bases with quality lower than 20 to lowercases (the 1st command line) or to N (the 2nd):

    seqtk seq -aQ64 -q20 in.fq > out.fa
    seqtk seq -aQ64 -q20 -n N in.fq > out.fa
  

• Fold long FASTA/Q lines and remove FASTA/Q comments:

    seqtk seq -Cl60 in.fa > out.fa
  

• Convert multi-line FASTQ to 4-line FASTQ:

    seqtk seq -l0 in.fq > out.fq
  

• Reverse complement FASTA/Q:

    seqtk seq -r in.fq > out.fq
  

• Extract sequences with names in file name.lst, one sequence name per line:

    seqtk subseq in.fq name.lst > out.fq
  

• Extract sequences in regions contained in file reg.bed:

    seqtk subseq in.fa reg.bed > out.fa
  

• Mask regions in reg.bed to lowercases:

    seqtk seq -M reg.bed in.fa > out.fa
  

• Subsample 10000 read pairs from two large paired FASTQ files (remember to use the same random seed to keep pairing):

    seqtk sample -s100 read1.fq 10000 > sub1.fq
    seqtk sample -s100 read2.fq 10000 > sub2.fq
  

• Trim low-quality bases from both ends using the Phred algorithm:

    seqtk trimfq in.fq > out.fq
  

• Trim 5bp from the left end of each read and 10bp from the right end:

    seqtk trimfq -b 5 -e 10 in.fa > out.fa
  

• Find telomere (TTAGGG)n repeats:

    seqtk telo seq.fa > telo.bed 2> telo.count
  

seqtk's Issues

Cannot use a fraction < 1 with seqtk sample

seqtk sample input.fastq 0.79 > /dev/null

yields:
Could not allocate enough memory for 0 sequences. Exiting...

fastq to fasta losing records

I'm trying to convert a large fastq file to fasta and the process is terminating partially through the file. I'm not running out of swap/memory and the file is a simple cat of a multi compressed fastq.gz. simple shell scripts can convert the file wo loss of records but it seems like seqtk seq -A file.fasta > file.fastq clips after a certain number. anyone else seen this?

Please tag a stable release

seqtk seems pretty stable. Could you please tag a stable release, 1.0 perhaps? Thanks.

Paired end information removed after trimming

Hi, I am using the trimfq function of seqtk on a set for paired end fastq data, but I notice that the paired end information is removed from the header after trimming. Is is possible to keep the header intact during trimming? The header is composed of read information, then a space, and then some more information like the paired end data: Example:
@TJH-ST1807:96:D0YGJAKXX:7:1101:2703:1986 2:N:0:CAAGTA

End after trimming:
@TJH-ST1807:96:D0YGJAKXX:7:1101:2703:1986

Thanks!

seqtk seq - need way split .fq into R1, R2 and orphan SE ?

Currently you can join R1 and R2 with "seqtk mergepe"
And you can get R1 and R2 separately with "seqtk seq -1" and "seqtk seq -2"
It would be more efficient to be able to get R1, R2 and SE in one pass.

Something like:

% seqtk seq -p out file.fq
% ls out.*
out_R1.fq out_R2.fq out_SE.fq

But left to your design discretion.

Trim to a specified length

I'd like to trim a bunch of variable length reads to a fixed length. I read the help page for seqtk trimfq, but didn't see what I was looking for. Could this feature please be added?

demultiplexing

I used this to demultiplex and it worked wonderfully. Given that the release is labeled 1.0-r68-dirty, I am wondering if the demultiplexing command has any known issues or needs further testing.

Add a LICENSE to the code

Please choose a license for the code. I would like to build tools that use seqtk, and would like to know what permissions are given.

If you have no preference, may I suggest the MIT license? :)

Ability to exclude reads

It would be useful to have a function to exclude sequences from a file. Basically it would be the exact opposite of subseq.

seqtk remove in.fq name.lst > out.fq

Only sequences with name NOT in name.lst are extracted.

seqtk

I did the following:

git clone https://github.com/lh3/misc.git
cd misc/seq/seqtk/

gcc -O2 seqtk.c -o seqtk -lz

But then got:

/tmp/ccErWOiV.o: In function stk_trimfq': seqtk.c:(.text+0x16ca): undefined reference topow'
collect2: ld returned 1 exit status

Any clue?

seqtk subseq and indexed FASTA files

I was surprised to find that seqtk subseq doesn't makes use of .fai indexed FASTA files. It warns that the user should use samtools faidx for a few sequences. Is there an intentional reason why indexed FASTA files aren't used? If not, I could work to implement this for seqk.

Add installation instructions to the README

My C experience is very limited. I did find install instructions that seems to work in #27 but adding them to the README would be nice.

$ git clone https://github.com/lh3/seqtk.git
$ cd seqtk
$ make

trimfq option not visible

I have cloned seqtk by

git clone https://github.com/lh3/seqtk.git

Then compiled the source using make on Linux 2.6.32-46-generic #107-Ubuntu.

I do not see the trimfq and seq options when running seqtk

Usage: seqtk

Command: comp get the nucleotide composite of FASTA/Q
hety regional heterozygosity
fq2fa convert FASTQ to FASTA
subseq extract subsequences from FASTA/Q
maskseq mask sequences
mutfa point mutate FASTA at specified positions
mergefa merge two FASTA/Q files
randbase choose a random base from hets
cutN cut sequence at long N
listhet extract the position of each het

suggestion for trimfq q=0

Would be possible that when running seqtk with "-q 0" (that is "seqtk trimfq -q0 -l N input.fq") the input FASTQ reads are trimmed to the first N bp automatically (and no trimming based quality score is done)?

This would be useful for input FASTQ files which have reads with different sizes.

Daniel

seqtk subseq -l option

UPDATE:
Forget it. I was calling a bad version.. !!! Sorry about that.

Hi there, thanks for that tool.

I'm trying to use the "-l" option to specify a length for the fastq sequence and qualities lines.

In the code, this option seems to be there, but when we execute it, it doesn't like it.
Any idea?

Thanks

trimming too much -> fail

Hello lh3,

great job on seqtk - thanks! But I think I just found a bug.

A collaborator handed me the fastq file - I guess some trimming/cleaning/removal has already been done. But I need to do more.

CCTCTTCTACTCACAGTCTGCCAGTGCCCGGAAAACTCGATGTTTCGTTCCAAGTCAAATACAGGACAAGATGTATTTTGTATAGATACAGTCATA
+
=DFFFFHHHHHJJJJJIIJJJJJJJJJJJJJJJJIJJJJJJJGIJJJIJJJJJJHIJJJIIJJJGHHHFFFFFECEEFEEEEEFEDEEDDDFDDEF
@HWI-ST885:90:D1B33ACXX:1:1101:2155:1995
AGGACC
+
=DDFFF
@HWI-ST885:90:D1B33ACXX:1:1101:2480:1958
AGAATCCCCCTTTACAATTGCTGACTCAATCTATGGATCTACCTTCTTTATTGCTACAGGCTTTCATGGCCTACATGTATTAATTGGAACAATTTT
+
=DDFFFHHHHHJJJJJJJJJJJJJJJJIJJJJJJJJJJJJJJJIJJJJJJJJJJJJJJJJIJJIJJJJJJJJIEHHHHHFFFFFFFEEEEEEDDDD

launching seqtk trimfq -b 2 -e 2 input.fastq > output.fastq leads to something horrible for the read that is too short.
Because this read isn't being eliminated (and I dont want it to be because otherwise the paired read will be an orphan), it should instead become an empty read (length 0)

Cheers
yannick

A question ： about fastq quality encode

Hi ,teacher ,
when I use GATK to treat my fastq reads , I meet such this error : SAM/BAM/CRAM file htsjdk.samtools.SamReader$PrimitiveSamReaderToSamReaderAdapter@60b2082a appears to be using the wrong encoding for quality scores: we encountered an extremely high quality score of 62

I search on the website ,I found somebody said I can use Seqtk to fix such this misencoded quals , But I
do not know how to do that ????
And why will appear such these misencoded quals in fastq reads ????

Expecting your answers !

Ability to process quality (.qual) files

It would be nice to have seqtk to be able to process quality files. Quality files are in the same format as FASTA files, except they have a space separating their quality scores. If seqtk is used to process any quality files, it will wrap the lines and merge two quality score values.

[deprekate@anthill ~]$ cat seq.fna 
>seq1
CCGAATGGATCATCCCGACTTTCAGGCCGGGATGGCCGGCCTGAAAGGGGACTGGGAACT
CCTCTGCCGCCCCTTGTGCGACCCGGATGCCCCGCGCGGCTGGCTGGGGGTCTGGGCGCT
[deprekate@anthill ~]$ cat seq.qual
>seq1
35 35 50 50 44 44 44 43 44 43 55 55 55 44 50 42 42 42 42 43 52 52 52 52 52
52 42 52 52 44 52 39 39 40 43 43 55 44 52 42 42 42 42 42 52 55 55 55 55 55

subseq shows the bug, The sampled nucleotide sequence has 120 bases, while the quality score has 119 qualities. You can see the 5252 in the middle of the quality that got merged.

[deprekate@anthill ~]$ seqtk sample seq.fna 1
>seq1
CCGAATGGATCATCCCGACTTTCAGGCCGGGATGGCCGGCCTGAAAGGGGACTGGGAACTCCTCTGCCGCCCCTTGTGCGACCCGGATGCCCCGCGCGGCTGGCTGGGGGTCTGGGCGCT
[deprekate@anthill ~]$ seqtk sample seq.qual 1
>seq1
35 35 50 50 44 44 44 43 44 43 55 55 55 44 50 42 42 42 42 43 52 52 52 52 5252 42 52 52 44 52 39 39 40 43 43 55 44 52 42 42 42 42 42 52 55 55 55 55 55

An easy fix would be to have the option to replace newlines with spaces instead of just removing them entirely?

Differences between setk and Biopython fastq-solexa conversion

Hi Heng,

I am using this toolkit to convert solexa (Illumina <1.3) to the newer format. I compared the output to the output generated by the Bio.SeqIO.convert function (Biopython toolkit) and I am seeing consistent differences of between 1-2 at low quality. Do you know why this might be happening?

Thanks,
Taylor

How to remove ambiguities?

I know you can do this with seqtk, but I cannot find the command in the readme file? I want to change all the ambiguious codes to "N" in my fasta file?

Maybe
seqtk seq -n N in.fasta > out.fasta??

Thanks!
B

fastq file with empty sequences

I am doing trimming/filtering using merTrim from Celera. This program does not remove sequences that it decides to delete from the input. Instead, it leaves an entry of length 0 in the fastq file. 'seqtk sample' stumbles upon these and reports a fasta header and one empty line for these in the fastq output. I can provide example data when needed.

unsorted output when using seqtk subseq

Is there a way to keep the subsetted sequences of say an extracted fasta file in the same order as the .lst file rather than sorting happen? For example I have contigs00xx and scaffolds00xx, and so even though contig and scaffold identifiers are interspersed throughout the .lst file, the final resulting fasta always has contigs first (essentially alpha numeric sorting). The issue arises as I'm trying to concatenate chunks of the draft assembly into linkage groups for downstream pop gen analysis. I've found alternatives using perl but generally find seqtk much much faster.

Explain seqtk comp columns

The usage says:

seqtk comp
Usage:  seqtk comp [-u] [-r in.bed] <in.fa>

Output format: chr, length, #A, #C, #G, #T, #2, #3, #4, #CpG, #tv, #ts, #CpG-ts

Can you please explain the meaning of #2 #3 #4 #CpG #tv #ts #CpG-ts? I read your code, but it is hard to digest for me.

The subseq function adds tabs to the sequence headers

The current version inserts a tab between the header ID and the comment. (I use the @ to show where tabs are)

$ seqtk 
Version: 1.0-r82-dirty

$ cat test.fna | tr "\t" "@"
>one comment
ACTG

$ cat test.lst | tr "\t" "@"
two

$ seqtk seq -A test.fna | tr "\t" "@"
>one comment
ACTG

$ seqtk subseq test.fna test.lst | tr "\t" "@"
>one@comment
ACTG

Line 569 of seqtk.c uses a tab instead of a space.

tabtk cut: loss of empty columns

When using tabtk cut for cutting and reordering columns, empty fields are lost and columns are shifted to the left in those rows that contain the empty fields.

the file I am using is an expression matrix, with in the first 2 rows sample identifiers, therefore no entry is present in the first field of those rows.
As simplified example I used following command:
startfile:

tabtk cut -r -f 1,5,2 startfile.txt > resultfile.txt

resultfile:

As can be seen, the tab in the first column is lost. Adding TAB as delimiter option did not solve this problem.
This is the wanted result. For this case it is easy enough to manually edit the file, but in more complex files this is not possible.

Fastq comment/pair info removed after subsampling

Hey guys,
subsampling fastq files to a specific number of reads removes the comment in the fastq header. I see that was also the case previously with trimming in issue #1.
$ zcat read1.fastq.gz | head
@H108:305:C0HA1ACXX:1:1101:1056:2056 1:N:0:TAAGGCGACTCTCTAT

$ seqtk sample -s100 read1.fastq.gz 10000 > test.fastq

$ head test.fastq
@H108:305:C0HA1ACXX:1:1202:3743:13412

All other subsampling functions work as expected, so this bug was probably missed. Thanks a lot.
Cheers

Move "bwa mergepe" to seqtk

I think the (experimental) bwa "mergepe" module would be more logically placed in seqtk.

subseq doesn't work

Hi, the following command doesn't give any output. Why?

seqtk subseq seq.fastq name.lst

The contents of the files look like this:

cat name.lst
@D2FC08P1:184:C0W67ACXX:1:1101:1173:2157
@D2FC08P1:184:C0W67ACXX:1:1101:1227:2179
cat seq.fastq
@D2FC08P1:184:C0W67ACXX:1:1101:1173:2157 1:N:0:
CTTCCTTCTGTCATTCCTTCTGGGCAAATTTTGGTTACACCCTTAGTAGAGCGTTTGGTGTGACCGTAGAGCTTCAAGACTTTGAATTTAATTGTAGATGTGATGAGAATGAGCTTTGTTATACGAAGAATGTCTTCAGTTTGCGTATCATGAGACAGTCTATGTAGGTAATATCTGTAGGGAATGTAC
+
BCCFFFFFHGHHHHIIIIIIIIIIIGIIIIIIIIHHHIIIIIIICHHIIHIIIIIIIIIIIIGGIHFEHFFDFFDC@>>CDDDD;;CDIIIIIIIIIIIIEDEFFFFFFFHHGHGIIIHEHHIIGIIIIGCIIIHHIHIIICIIIIIHIIIIIIIIIIIIIIIIIIIIIIIIIIIIGHHHHFFDFF@CB
@D2FC08P1:184:C0W67ACXX:1:1101:1227:2179 1:N:0:
TTTTCTGATGGTTCAACAAATAGGAAAGGCTATTCTAAGAAAGCAAAAATAATTGTCAAACATTTGACTCATGCATGCGATGAATCGCTCGCATCATCAGCTCCATGCGTGTGCCAATGCGTGCCCCTGTGATCGCACGCATTAAAAGACGAAGACTGAAAATATAATGGTTACATGTTTTAT
+
B@@ffff:CFHHHIIIIIIIGHIIGIIIIIIIIIIIIIIIIIIIIIFIIIIIIIIIIGGIIIIIIIGGHIIIIIIIIGFHHFIIIIIIIIIIIIIIIIIIICDDBBDDDDCCEEEFFDHHHIIIIIIIIHIHGGGG@IIGGGIIIIIIIHFIIIIIIIIIIIIIIIIIIIHHHHHFFDFFCCB

grep does find a match.

grep $(head -1 name.lst) seq.fastq
@D2FC08P1:184:C0W67ACXX:1:1101:1173:2157 1:N:0:

The example mentioned in #5 does work in my environment though.

Thanks in advance.

Takeo

Not an issue... Permission to add seqtk to Archlinux User Repositories?

Hi,
with your permission I would like to add this software to the Archlinux distribution user repository (https://aur.archlinux.org/packages/seqtk-git). It means Archlinux users will be able to install your software seamlessly and have it updated when you make commits.
I assumed that this software is distributed under GPL license (from your other Github contributions), but I have not found any paper presenting it- if you want to be cited, please tell me how. I have made sure to mention Heng Li as the author.
If instead you would like me to take this Archlinux package down, I will do it at once (answer here or email me at [email protected]).
Thanks again for your great work!
Cyril

Read input from stdin

I was toying with seqtk as an efficient and simple fastq-to-fasta converter and think it would be really nice if it could read from stdin, like many GNU utilities (and also e.g. EMBOSS seqret).
Now I resorted to creating a FIFO, piping data through that and that works perfectly well, but it would be nice not to have to pipe data through a named pipe like that.

subseq not using the strand information in bed

Hi,
seqtk is so powerful and helpful to my regular work.
One problem for subseq tool is that it does not using the strand information in the bed file. I have tested using the same position with different strand information with same returned results.

Best,
Pengcheng

writing XA tags for primary alignments in bwa-mem bams

Hi all,

As I am trying to switch from bwa to bwa-mem I realize that bwa-mem does not populate the XA tag in the primary alignment to report alternative mapping positions.
I understand that this information is implicit if a read has many entries, as each corresponds to an alternative mapping position, however I think that it would be useful to preserve this format since some downstream analyses rely on it.
Perhaps I am missing a parameter that allows to turn this behavior on, but if that is not implemented I suggest it as a feature :)

Cheers,

Elizabeth

A bug about seqtk

Hi
I use seqtk to convert my fastq : seqtk seq -V -Q 64 $1.fq > $1.s.fq
and Use fastqQc , but I meet such this error :
"FATAL ERROR (file fastq_qc_stat.c line 80): quality (-1) not in [0, 45]"

What happend about seqtk ? Thank you !

How to interlace pair-end reads via seqtk and make the headers end in /1 and /2 corresponding to its read orientation (R1 and R2, respectively)?

Hello, sorry for the empty previous post cuz I tapped Enter key so fast...
I am preparing datasets for kmernorm which requires the files look like
@seq_1/1
AAAAAAAACCCCCCCTTTTTTTTTGGGGGGGG
+
&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
@seq_1/2
AAAAAAAACCCCCCCTTTTTTTTTGGGGGGGG
+
&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
@seq_2/1
AAAAAAAACCCCCCCTTTTTTTTTGGGGGGGG
+
&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
@seq_2/2
AAAAAAAACCCCCCCTTTTTTTTTGGGGGGGG
+
&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

And a bioinformatician told me seqtk could archive this. But after screening the examples of seqtk, I didn't find the command. Could someone tell me how to do this in seqtk? Many thanks!

problem installing/using seqtk

I am new to seqtk, and installed it following the directions of a previous thread. I am using a MacOSX, and in my terminal typed the following:

$ git clone https://github.com/lh3/seqtk.git
$ cd seqtk
$ make

however, I get the following error message: "-bash: make: command not found"

I would appreciate some advice on what I may be doing wrong or what I can do to use seqtk. Thank you

seqtk trimfq -l 50 in.fq > out.fq have much reads longer than 50bp.

subseq => core dumped

It looks like when there is not enough free memory available for "seqtk subseq", then a crash happens with this error message: "Segmentation fault (core dumped)".

For example running this:

seqtk subseq reads-of-size-35-GB.fq list-of-size-20-GB.txt > output.fq

on a computer with 16GB RAM will produce this error messge "Segmentation fault (core dumped)".

The expected behavior would be to end graciously with a message like "not enough memory" and so on!

seqtk truncating file

I am having issues with seqtk truncating a file while changing the quality encoding.

I have describibed this on repo of the pipeline that is using seqtk, bcbio/bcbio-nextgen#705

Given that you develop seqtk, your help on this would be great!!

Removing whitespace from sequence

Can we strip whitespace from sequences?

mutfa: Add option to allow mutating a nucleotide with the same case as the reference nucleotide

Can an option (e.g. --keepcase) being added to mutfa to allow keeping the case of the reference nucleotide for the introduced mutation?

FASTA file:

$ cat test.fa 
>seq1
AAAACCCCGGGGTTTT
>seq2
aaaaccccggggtttt

File with mutations:

$ cat in.snp
seq1    2   bla T
seq1    8   bla A
seq1    15  bla a
seq2    2   bla T
seq2    8   bla A
seq2    15  bla a

Normal mutfa command:

$ seqtk mutfa test.fa in.snp 
>seq1
ATAACCCAGGGGTTaT
>seq2
aTaacccAggggttat

Change reference nucleotide, but keep case (uppercase/lowercase) of reference nucleotide:

$ seqtk mutfa --keepcase test.fa in.snp 
>seq1
ATAACCCAGGGGTTAT
>seq2
ataacccaggggttat

Force reduced memory usage in random sampling

According to https://www.biostars.org/p/6544/#42695, seqtk will only use the memory footprint of one record when using a fractional sample size. Is there a way to force this behaviour but with giving an integer part instead of a fractional part?

Allow for appending to sequence id with RENAME

With the RENAME function you replace the sequence id with an incrementing integer (and optional prefix).

It would be nice to be able to just append the incrementing integer to the sequence id. Other wise a separate file is needed to map an original read name to the new read name.

Please tag a release that includes tabtk

Thanks!

seqtk doesn't support the full range of ambiguity codes.

E.g. for a test sequence:

>ambig
ACGTMRWSYKVHDBN

We'll see the following bases identified as ambiguous:

#chr    position  base
ambig   5         M
ambig   6         R
ambig   7         W
ambig   8         S
ambig   9         Y
ambig   10        K

It would be great if seqtk could support the following IUPAC ambiguities:

B   C or G or T
D   A or G or T
H   A or C or T
V   A or C or G
N   any base

If I wrote any C at all, I'd make a PR, but as I don't I'll just leave an issue. It looks like the offending line is code like L812 and L438 which checks if the value obtained from bitcnt_table == 2, whereas there are ambiguities with values 3 and 4.

Tag a recent stable release

Release 1.0 corresponds with r31, and the current commit of seqtk is r68. Please tag a recent stable release.

How to interlace pair-end reads via seqtk?

sample adds ' correct' suffix to read name

seqtk sample adds ' correct' suffix to read names at random.

Command used:

seqtk sample -s 100 fwd.fq 350000

problem of using subseq

Hi,
I try to extract sequences from a fastq file using the "subseq" in seqtk.
But the extract file contains only the 1st sequence but no others.
I am wondering whether my name.lst file does not fit with what seqtk needs.
I have names of each sequence without other symbols each line in the name.lst. But the fastq file starts each sequence name with a @. Should I add @ in front of each sequence name?
Or what other problem it can be?

Any suggestion is welcome.
Thanks,

Chih-Ming

ksort.h is missing

$ make
gcc -g -Wall -O2 tabtk.c -o tabtk -lz -lm
tabtk.c:14:19: fatal error: ksort.h: No such file or directory
compilation terminated.
make: *** [tabtk] Error 1

I copied ksort.h from fermi2 repository, and compiled without error. The resulting tabtk works too.

Add in support for Mersenne Prime Twister PRNG

The current random number generator used is a linear congruential generator and is not a high quality random number generator. A brief discussion on random number generators can be found here with references to additional sources. If possible, it would be fantastic if seqtk relied on the Mersenne Twister, the source of which can be found here or alternatively, the algorithm can be used from the GNU Scientific Library.