search

lh3 / seqtk

Toolkit for processing sequences in FASTA/Q formats
MIT License
1.35k stars 310 forks source link

buggy behavior with seqtk subseq command #208

Closed franciscozorrilla closed 4 months ago

franciscozorrilla commented 4 months ago

I noticed some buggy behavior with seqtk (Version: 1.4-r122) subseq command misfolding sequences for contig headers containing - in list. As you can see below, contigs/sequences with a hyphen end up with modified headers and single amino acid length:

$ head headers_clean.tsv 
MGYG000000001_00051 hypothetical protein
MGYG000000001_00092 Hypoxanthine phosphoribosyltransferase
MGYG000000001_00096 4-hydroxy-tetrahydrodipicolinate reductase
MGYG000000001_00097 4-hydroxy-tetrahydrodipicolinate synthase
MGYG000000001_00098 Aspartate-semialdehyde dehydrogenase 2
MGYG000000001_00111 D-hydantoinase
MGYG000000001_00123 Pyruvate synthase subunit PorC
MGYG000000001_00124 Pyruvate synthase subunit PorD
MGYG000000001_00125 Pyruvate synthase subunit PorA
MGYG000000001_00126 Pyruvate synthase subunit PorB

$ seqtk subseq -l 60 fasta.faa headers_clean.tsv > misfolded.faa

$ head misfolded.faa -n20
>MGYG000000001_00051 hypothetical protein
MKTKRLTFLSLMIGYSLILYIVETYIPNPLIILFPGAKLGLTNIITLISLILLGGKDTFI
IVTIRVILSSIFSGPLSYLLFSIGGAYLSLVLMYIISKIDGFSLIGVSIIGAIGHNIGQL
IVASIIVENALTIGYLPFMLTASLVTGFFVGMVTKYCIPLTGKLCKIH
>MGYG000000001_00092 Hypoxanthine phosphoribosyltransferase
MDIETKKWEVLCSKEQIDERLKELGSQISEDYKDKNLLVVSLLKGSFIFCADLVRQINVP
VRIEFMTTASYGHGEESSGTVNVVSDVKADLEGYDVLVVDDITDSALTMDFVMKHLSSKN
PASIKSCVLLDKPSRRKVDLVPDYCGFTIEDKFVVGYGLNFGDYYRNVPYVFNVTNEDR
>MGYG000000001_00096:4-4 4-hydroxy-tetrahydrodipicolinate reductase
V
>MGYG000000001_00097:4-4 4-hydroxy-tetrahydrodipicolinate synthase
K
>MGYG000000001_00098 Aspartate-semialdehyde dehydrogenase 2
MFKLSFSKYILIGVIKMKKVNVAIVGATGMVGRTFLKVLEERKFPINNLYLFSSAKSAGK
VVKFAGKDYVVEELTENSFDRDIQIALFSAGGAISEKFAPIAASKGVLVVDNSSYWRMDP
KVPLVVPEVNPEAVKEHNGIIANPNCSTIQAVVPLKPLHEKYKLKRIVYSTYQAVSGSGV
KGVKDLENGINGEANKFYPYQIAYNCLPHIDVFTDNGYTKEEMKMVNETMKIFNDYDIKV
TATTVRVPVKNCHSESINVELKNPFELNDVVEELKNSDGIVVIDNTEKLEYPTTLDSNGH
DEVFVGRIRRDFSIDNGVNFWCVADNVRKGAATNAVQIAELLIKYDLV
>MGYG000000001_00111 D-hydantoinase

Same thing happens when just using the contig headers as they are from the fasta file:

$ head headers.txt 
MGYG000000001_00051 hypothetical protein
MGYG000000001_00092 Hypoxanthine phosphoribosyltransferase
MGYG000000001_00096 4-hydroxy-tetrahydrodipicolinate reductase
MGYG000000001_00097 4-hydroxy-tetrahydrodipicolinate synthase
MGYG000000001_00098 Aspartate-semialdehyde dehydrogenase 2
MGYG000000001_00111 D-hydantoinase
MGYG000000001_00123 Pyruvate synthase subunit PorC
MGYG000000001_00124 Pyruvate synthase subunit PorD
MGYG000000001_00125 Pyruvate synthase subunit PorA
MGYG000000001_00126 Pyruvate synthase subunit PorB

Correct behavior by removing extra info from contig headers:

$ head headers_clean2.tsv 
MGYG000000001_00051
MGYG000000001_00092
MGYG000000001_00096
MGYG000000001_00097
MGYG000000001_00098
MGYG000000001_00111
MGYG000000001_00123
MGYG000000001_00124
MGYG000000001_00125
MGYG000000001_00126

$ seqtk subseq -l 60 fasta.faa headers_clean2.tsv > folded.faa

$ head gutdbx_v2_folded.faa -n 20
>MGYG000000001_00051 hypothetical protein
MKTKRLTFLSLMIGYSLILYIVETYIPNPLIILFPGAKLGLTNIITLISLILLGGKDTFI
IVTIRVILSSIFSGPLSYLLFSIGGAYLSLVLMYIISKIDGFSLIGVSIIGAIGHNIGQL
IVASIIVENALTIGYLPFMLTASLVTGFFVGMVTKYCIPLTGKLCKIH
>MGYG000000001_00092 Hypoxanthine phosphoribosyltransferase
MDIETKKWEVLCSKEQIDERLKELGSQISEDYKDKNLLVVSLLKGSFIFCADLVRQINVP
VRIEFMTTASYGHGEESSGTVNVVSDVKADLEGYDVLVVDDITDSALTMDFVMKHLSSKN
PASIKSCVLLDKPSRRKVDLVPDYCGFTIEDKFVVGYGLNFGDYYRNVPYVFNVTNEDR
>MGYG000000001_00096 4-hydroxy-tetrahydrodipicolinate reductase
MLKVIINGCSGKMGRVLARCINEESDLELLCGVSPQNTEDLNFKTYSTMSEINEKSDVVI
DFSHHSTLDDVLSYTIKTKTPLVIATTGYNDEELGKIHEASKIIPIFHSYNMSLGVNVLL
RLVKEAANVLSNFDIEVIEKHHNKKVDAPSGTAVMIANAIKEVLPSLENNYGRYGRESKR
KENEIGIHAIRGGTIVGEHDVIFAGHDEIVEIKHTAQSKDIFAKGSIAAAKFIVVQEAGY
YNMDNMF
>MGYG000000001_00097 4-hydroxy-tetrahydrodipicolinate synthase
MIFKGSAVALITPFNENNEVDFEKLGELVEYHIENNTDAIVACGTTGEASTLDDAEHLAV
IKYVIDKVNKRIPVIAGTGGNDTQHSIYMSQEAEKFGVDGLLVITPYYNKGNKSGIKKHF
EAIAASVKTPIIMYNVPGRTCVNMSPALVAELAKIDNIVAIKEASGDLAQVAEIARLVPD
DFAIYSGNDDSILPLLSLGGSGVISVLANVCPQETHDLVYKFFDGDLKGSRDLQLGMKPL
IDALFIEVNPVPVKTAANLLGFNVGGLRLPLAEMDSKNVEILKKELINWGLNVQEETC

Hope this helps prevent issues for other users. Best, Francisco

lh3 commented 4 months ago

This is expected. Either one-column list or three-column BED.