This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Getting coordinates for a given sequence motif

I would like to get all coordinates for a given series of short sequence motifs, e.g. GATTACA (and its reverse complementary TGTAATC) from a reference genome, but want to do it in a way that allows me to get the coordinates of the sampled motif (because I need to intersect these coordinates with other annotation).

My current approach is based on a strict blastn-short search of that motif. However, this is extremely inefficient, because the minimum length for the algorithm to work requires sequences of a minimum length of 15 bp, so I need to append all possible combinations of nucleotides to the flanks of the motif to reach that minimum length. As consequence, the blast search is not just on a couple or a few motifs, but on thousands of them. For shorter targets, more than inefficient is is plain impractical.

There must be a way to do this more straightforwardly, but have doubts for longer ones, and I am wondering if anybody here has a better idea on how to address this problem.

Thank you,

EDITED: what I want really is not sampling, but getting all genomic coordinates for these motifs.

genome assembly

2 answers

In addition to the options presented in ATpoint's link, seqkit locate will work also. It's become my go-to software for many common tasks for fasta and fastq files.

Yeah, this is probably the best solution as it is an efficient command line tool and the solutions in the thread I linked are rather custom code. I always forget about seqkit but it is a super handy and efficient tool, swiss army knife for fasta/q files.

From the seqkit link:

# max mismatch: 1
$ cat t.fa \
  | seqkit locate -p agc -m 1 \
  | csvtk pretty -t
seqID   patternName     pattern strand  start   end    matched
seq     agc           agc       +        1       3     agc
seq     agc           agc       +        7       9     agc
seq     agc           agc       +        11      13    acc
seq     agc           agc       -        8       10    agc
seq     agc           agc       -        2       4     agc

find positions of a short sequence in a genome

This looks promising, thanks!

Log in to answer this question.