Skip to content

find

Jump to bottom
markschl edited this page Aug 15, 2018 · 8 revisions

Fast searching for one or more patterns in sequences or ids/descriptions, with optional multithreading.

Usage:
  st find [options] [-a <attr>...][-l <list>...] <pattern> [<input>...]
  st find (-h | --help)
  st find --help-vars

Search Options:
    <pattern>           Pattern string or 'file:<patterns.fasta>'
    -r, --regex         Treat the pattern(s) as regular expressions.
    -d, --dist <dist>   Approximative string matching with maximum edit distance
                        of <dist> [default: 0]
    --in-order          Report hits in the order of their occurrence instead
                        of sorting by distance (with -d > 0)
    --seqtype <type>    Sequence type {dna/rna/protein/other}
    -t, --threads <N>   Number of threads to use [default: 1]
    --ambig <yn>        Override choice of whether DNA ambiguity codes (IUPAC)
                        are recognized or not {yes/no}.
    --algo <algorithm>  Override decision of algorithm for testing
                        (regex/exact/myers/auto) [default: auto]

Search range:
    --rng <range>       Search within the given range ('start..end', 'start..'
                        or '..end'). Using variables is not possible.
    --max-shift-l <n>   Consider only matches with a maximum distance of <n> from
                        the search start (eventually > 1 if using --rng)
    --max-shift-r <n>   Consider only matches with a maximum distance from the
                        end of the search range

Attributes:
    -i, --id            Search / replace in IDs instead of sequences
    --desc              Search / replace in descriptions

Actions:
    -f, --filter        Keep only matching sequences
    -e, --exclude       Exclude sequences that matched
    --dropped <file>    Output file for sequences that were removed by filtering.
                        The extension is autorecognized if possible, fallback
                        is the input format.
    --rep <with>        Replace by a composable string

See this page for the options common to all commands.

Contents

Description

The find command allows searching for patterns in sequences or sequence headers. It can be used for simple filtering, but results can also be returned by using the variables. Occurrences can be replaced.

Exact searching

This is the simplest and also the fastest way of searching. In the following example, a pattern is searched and only hits matching sequences are returned (-f). Unmatched sequences are written to a different file.

st find -s AATG -f --dropped not_matched.fa seqs.fa

Regular expressions

Especially useful when searching in sequence headers. Different regex groups can be accessed as output variables.

The following example extracts Genbank accessions from sequences in an old style Genbank format by using a regex group named 'acc' and replaces the ID with this group:

st find -ir "gi\|\d+\|[a-z]+\|(?P<acc>.+?)\|.*" seqs.fa \
   --rep "{f:match::acc}" > seqs_accession.fa

Headers like this one: >gi|1031916024|gb|KU317675.1| are transformed to >KU317675.

Note: Only ASCII is currently supported, unicode is not recognized in regular expressions.

Approximative searching

Approximative matching is particularly useful for finding primers and adapters or other short patterns. The tool is not specialized for removing adapters like other software, but still aims at being very fast and useful for many purposes. Note: The length of search patterns for approximative matching is limited to 64 characters, currently.

It is possible to search for all matches up to a maximum edit distance (-d/--dist argument). By default, the best hit is reported first. Example:

st find -d 4 ATTAGCG seqs.fa \
     -a hit="{f:range}_(dist:{f:dist})" -a range="{m:range}"

Possible output:

>seq464 hit=2-9_(dist:1) matched=ATCAGCG
GGATCAGCGATCC
(...)

The second best hit can be returned by using {f:range:2} or {f:match:2}, etc. Use --in-order to report hits in order from left to right instead.

Substantial speedups can be achieved by restricting the search range and by multithreading (-t).

Search performance

Approximative matching uses the fast bit-parallel algorithm by Myers. The runtimes for searching two forward primers (22 bp) with up to 4 mismatches (-d 4) in a 1.2 GB file vary depending on the options used. In the following table, different settings as well as other tools are compared, using 1 or 4 threads or processes ('cores').

1 core 4 cores
Find the position of best hit 31.3s 7.73s
Report first hit from left1 31.1s 7.72s
Best hit in range where the primer should occur2 8.02s 2.34s
Filter by occurrence (-f), no position determined 14.0s 3.24s
Exact search (no mismatches/ambiguities)3 11.1s 2.57s
cutadapt 1min 16s 24.2s
AdapterRemoval4 35.0s 8.75s

1 --in-order option. Normally, hits are sorted by decreasing distance.

2 --rng ..25

3 Using the Two Way algorithm for comparison, not Myers

4 Single-end mode

Ambiguities

DNA (RNA) ambiguity codes according to the IUPAC nomenclature are accepted and automatically recognised in search patterns. For Proteins, X is recognised as wildcard for all amino acids. The molecule type is automatically determined from the pattern. Use -v to show which search settings are being used. Example:

st find -v file:primers.fasta -a primer={f:name} -a rng={f:range} input.fasta > output.fasta
primer1: DNA, search algorithm: Exact
primer2: DNA with ambiguities, search algorithm: Myers
Sorting by distance: true, searching start position: true

In case of wrongly recognised patterns, specify --seqtype. The tool is very cautious and will warn about any inconsistencies between multiple patterns.

Note: Matching is asymmetric: R in a search pattern matches [A, G, R] in sequences, but R in a sequence will only match ambiguities sharing the same set of bases (R, V, D, N) in the pattern. This should prevent false positive matches in sequences with many ambiguous characters.

Multiple patterns

Several patterns can be searched in one command. They have to be supplied in a separate FASTA file. The best matching pattern with the smallest edit distance is always reported first. Other patterns are accessed using <variable>.<pattern_num>:

st find -d6 file:f_primers.fa seqs.fa \
    -p f_primer={f:name} -p f_dist={f:dist} \
    -p second_best={f:name.2}_({f:dist.2}) > primer_search.fa

Example output:

>seq1 f_primer=primer2 f_dist=1 second_best=primer1_(5)
SEQUENCE

Restrict search range

It is possible to search only part of the target string by using --rng. This can substantially speed up the search.

st find -d6 --rng ..23 file:f_primers.fa seqs.fa \
    -p f_primer={f:name} > primer_search.fa

If the hit is known to occur at the start or end of the search range, --max-shift-l and --max-shift-r can be used to report only those hits.

Replace matches

Hits can be replaced by other text (--repl). Variables are allowed as well (in contrast to the replace command). Backreferences to regex groups (e.g. $1) are not supported like the replace command does. Instead, they can be accessed using variables (<variable>::<group>)

Variables

Variables for the matched sequence (f:match), coordinates (f:start, f:end, f:range, etc.) are available. Selecting another match than the best/first one is possible by adding the match number after a colon (:): f:range:2 will return the second match. Match groups of regular expressions can be specified as well by the addition of a second colon: f:range::2 will select the second match group of the first match. Even named groups are possible. In that case, specify the group name instead of the index: f:range::<group_name>.

A more generalized scheme:

m:<variable>.<pattern_rank>:<hit_num>:<match_group>

This is admittedly quite complicated, but adds a lot of flexibility.

It is also possible to return all hits instead of a specific one by using f:<variable>:all. This will return a comma delimited list. The following command searches for all occurrences of a pattern and converts them to lowercase:

st find -r -p rng={f:drange:all} [AG]GA seqs.fa \
  | st mask {a:rng}

Exmaple output:

>seq464 rng=6..8,14..16
AGTTAagaCTTAAggaT

Provided variables

Pattern finding variables. Usage: f:<variable>[.pattern_rank][:match_num][:group]
---------------------------------------------------------------------------------
f:match      The matched sequence/pattern
f:start      Start of the match.
f:end        End of the match.
f:dist       Distance of the matched sequence compared to the pattern. Normally,
             this is the edit distance, unless --gapw is used
f:neg_start  Start of the match relative to sequence end (negative number)
f:neg_end    End of the match relative to sequence end (negative number)
f:range      Range of the match in the form start-end
f:drange     Range of the match with two dots as delimiter (start..end). Useful
             with 'trim'and 'mask'
f:neg_drange Range of the match (dot delimiter) relative to the sequence end
             (-<start>..-<end>)
f:name       Name of the best matching pattern if there are multiple (read from
             pattern file)

Example:

Clone this wiki locally