lh3 / seqtk Goto Github PK
View Code? Open in Web Editor NEWToolkit for processing sequences in FASTA/Q formats
License: MIT License
Toolkit for processing sequences in FASTA/Q formats
License: MIT License
seqtk sample adds ' correct' suffix to read names at random.
Command used:
seqtk sample -s 100 fwd.fq 350000
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
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.
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
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.
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.
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
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
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
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
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.
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
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 !
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!
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.
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
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!
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?
$ 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.
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.
Thanks!
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.
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
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 !
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
seqtk seems pretty stable. Could you please tag a stable release, 1.0 perhaps? Thanks.
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
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?
Release 1.0 corresponds with r31, and the current commit of seqtk is r68. Please tag a recent stable release.
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
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.
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
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?
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.
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.
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
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.
Can we strip whitespace from sequences?
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
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 to
pow'
collect2: ld returned 1 exit status
Any clue?
I think the (experimental) bwa "mergepe" module would be more logically placed in seqtk.
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
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!!
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!
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? :)
seqtk sample input.fastq 0.79 > /dev/null
yields:
Could not allocate enough memory for 0 sequences. Exiting...
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.
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.
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?
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.