seqkit locate misses embedded patterns

[ZacharyArdern]

Describe your issue

Seqkit locate does not seem to find fully overlapping patterns on the same strand (whether in greedy or non-greedy mode)

For instance, the regex below for searching for ORFs, adapted from that provided, does not find fully overlapping (embedded) ORFs on the same strand.

#with genome file "test.txt":
cat test.txt | seqkit locate -i -p "[AT]TG(?:.{3})+?T" -r | awk '{if (($6-$5)>=60 && NR>=2 ) print }' > ORFs.txt ;

e.g. this predicts an ORF from positions 123-260 on the plus strand, but misses the embedded ORF from 128-247 (found e.g. with prodigal, in test_all.txt)
prodigal -i test.txt -g 1 -f gff -s test_all.txt > test.gff ;

faidx test.txt LVMU01000403.1:123-260 ;
faidx test.txt LVMU01000403.1:128-247 ;

test.txt

[shenwei356]

For some reasons I can't recall, I chose to remove these "embedded" regions when searching with regular expressions before. It looks unnecessary now, cause we already have the -G/--non-greedy flag.

Thanks for reporting :)

• seqkit_darwin_arm64.tar.gz
• seqkit_linux_amd64.tar.gz
• seqkit_darwin_amd64.tar.gz

$ seqkit faidx test.txt LVMU01000403.1:123-260 \
    | seqkit locate -r -i -p "A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)" -M \
    | csvtk cut -t -f -patternName

seqID   pattern strand  start   end
LVMU01000403.1:123-260  A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)        +       1       138
LVMU01000403.1:123-260  A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)        +       6       125   * i.e. region 128-247
LVMU01000403.1:123-260  A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)        +       82      138
LVMU01000403.1:123-260  A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)        +       106     138
LVMU01000403.1:123-260  A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)        -       36      47


$ seqkit faidx test.txt LVMU01000403.1:123-260 \
    | seqkit locate -r -i -p "A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)" -M -G \
    | csvtk cut -t -f -patternName

seqID   pattern strand  start   end
LVMU01000403.1:123-260  A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)        +       1       138
LVMU01000403.1:123-260  A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)        -       36      47

[ZacharyArdern]

Thank you for the quick response again @shenwei356 !

I have tested it, and see just one issue, with the example you give for finding ORFs, rather than seqkit locate

• the regex you give also finds sequences where a stop immediately follows the start codon.

e.g.

LVMU01000403.1:12414-12659(+)
ATGTAAACCGACACCCGTCGCCTGAACGATACGTTTAAGTGTCCTTTGTTTGATCATCGTATTATCTCGCCAAATTACCTATCCAACCGAAGTGTACTATACATTCGGCGGGCCAGTTTAGCACAAAGAGCCTCGAAACCCAAATTCCAGTCAATTCTTAATCAGCTTGCTTACGCAGGAATGCTGGGATATCCAGATAATCCGGCTCTTTCGCAGTTTGCGGCGCATTGTCATTCACGACTTTAG

These few cases can be filtered out easily afterwards, so this is just a minor comment in case you know how to improve the regex so it excludes these non-ORF cases.

[shenwei356]

the regex you give also finds sequences where a stop immediately follows the start codon

I'm afraid it's not.

A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)
      ----------
      It means at least one triple-base codon.
      +  means one or more

$ echo -ne ">s\nATGTAA\n" | seqkit locate -r -i -p "A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)"
seqID   patternName     pattern strand  start   end     matched

$ echo -ne ">s\nATGcccTAA\n" | seqkit locate -r -i -p "A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)"
seqID   patternName     pattern strand  start   end     matched
s       A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)        A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)        +       1       9       ATGcccTAA

[ZacharyArdern]

@shenwei356 my statement wasn't fully clear, but this regex includes a stop within an "ORF" if it is immediately after the start codon:

echo -ne ">s\nATGTAAGGGTGA\n" | seqkit locate -r -i -p "A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)"
seqID	patternName	pattern	strand	start	end	matched
s	A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)	A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)	+	1	12	ATGTAAGGGTGA

[shenwei356]

I see. The regular expression is not perfect for finding ORF.