From 8ba6abe9e1f6d8c2a0d40ab165880adda5e27241 Mon Sep 17 00:00:00 2001 From: Jakob Nybo Andersen Date: Sat, 29 Aug 2026 15:07:43 +0200 Subject: [PATCH 1/9] Documentation: Language and grammar fixes Had an AI agent fix up grammar, typos and formatting. --- docs/src/construction.md | 56 +++++++++++++++------------ docs/src/counting.md | 14 +++---- docs/src/index.md | 8 ++-- docs/src/interfaces.md | 22 +++++------ docs/src/io.md | 6 +-- docs/src/random.md | 6 +-- docs/src/recipes.md | 1 - docs/src/sequence_search.md | 57 ++++++++++++++-------------- docs/src/symbols.md | 6 +-- docs/src/transforms.md | 23 ++++++----- docs/src/types.md | 39 ++++++++++--------- src/BioSequences.jl | 8 ++-- src/alphabet.jl | 52 ++++++++++++------------- src/biosequence/biosequence.jl | 28 +++++++------- src/biosequence/copying.jl | 10 ++--- src/biosequence/predicates.jl | 8 ++-- src/biosequence/transformations.jl | 34 ++++++++--------- src/bit-manipulation/bitindex.jl | 2 +- src/counting.jl | 25 ++++++------ src/longsequences/copying.jl | 10 ++--- src/longsequences/longsequence.jl | 12 +++--- src/longsequences/operators.jl | 32 ++++++++-------- src/longsequences/randseq.jl | 26 ++++++------- src/longsequences/stringliterals.jl | 12 +++--- src/longsequences/transformations.jl | 10 ++--- src/search/ApproxSearchQuery.jl | 8 ++-- src/search/ExactSearchQuery.jl | 14 +++---- src/search/pwm.jl | 4 +- src/search/re.jl | 30 ++++++--------- test/search/regex.jl | 4 ++ 30 files changed, 283 insertions(+), 284 deletions(-) diff --git a/docs/src/construction.md b/docs/src/construction.md index d3bae0058..966229d3f 100644 --- a/docs/src/construction.md +++ b/docs/src/construction.md @@ -7,8 +7,8 @@ end # Construction & conversion -Here we will showcase the various ways you can construct the various sequence -types in BioSequences. +Here we will showcase the various ways you can construct the sequence types in +BioSequences. ## Constructing sequences @@ -35,7 +35,7 @@ UUAGC ``` -Type alias' can also be used for brevity. +Type aliases can also be used for brevity. ```jldoctest julia> LongDNA{4}("TTANC") @@ -72,7 +72,8 @@ TTAGC ### Constructing sequences from other sequences -You can create sequences, by concatenating other sequences together: +You can create sequences by concatenating other sequences together: + ```jldoctest julia> LongDNA{2}("ACGT") * LongDNA{2}("TGCA") 8nt DNA Sequence: @@ -88,7 +89,7 @@ TATATATATATATATATATA ``` -Sequence views (`LongSubSeq`s) are special, in that they do not own their own data, +Sequence views (`LongSubSeq`s) are special in that they do not own their own data, and must be constructed from a `LongSequence` or another `LongSubSeq`: ```jldoctest @@ -112,7 +113,9 @@ true ## Conversion of sequence types -You can convert between sequence types, if the sequences are compatible - that is, if the source sequence does not contain symbols that are un-encodable by the destination type. +You can convert between sequence types if the sequences are compatible—that is, +if the source sequence does not contain symbols that are unencodable by the +destination type. ```jldoctest julia> dna = dna"TTACGTAGACCG" 12nt DNA Sequence: @@ -201,7 +204,7 @@ CTTA ``` -However, this is not what happens, instead the following happens: +However, this is not what happens; instead, the following happens: ```@meta DocTestSetup = quote @@ -228,21 +231,21 @@ CTTAAA ``` -The reason for this is because the sequence literal is allocated only once -before the first time the function `foo` is called and run. Therefore, `s` in +The reason for this is that the sequence literal is allocated only once +before the first time the function `foo` is called. Therefore, `s` in `foo` is always a reference to that one sequence that was allocated. -So one sequence is created before `foo` is called, and then it is pushed to -every time `foo` is called. Thus, that one allocated sequence grows with every -call of `foo`. +So one sequence is created before `foo` is called, and then a nucleotide is appended +to it on every call to `foo`. Thus, that one allocated sequence grows with every +call to `foo`. -If you wanted `foo` to create a new sequence each time it is called, +If you want `foo` to create a new sequence each time it is called, then you can add a flag to the end of the sequence literal to dictate behaviour: -A flag of 's' means 'static': the sequence will be allocated before code is run, +A flag of `s` means 'static': the sequence will be allocated before the code is run, as is the default behaviour described above. -However providing 'd' flag changes the behaviour: 'd' means 'dynamic': +However, providing the `d` flag changes the behaviour: `d` means 'dynamic': the sequence will be allocated whilst the code is running, and not before. -So to change `foo` so as it creates a new sequence -each time it is called, simply add the 'd' flag to the sequence literal: +So to change `foo` so that it creates a new sequence +each time it is called, simply add the `d` flag to the sequence literal: ```@meta DocTestSetup = quote using BioSequences @@ -290,11 +293,11 @@ DocTestSetup = quote end ``` -So the take home message of sequence literals is this: +So the take-home message of sequence literals is this: Be careful when you are using sequence literals inside of functions, and inside the bodies of things like for loops. And if you use them and are unsure, use the - 's' and 'd' flags to ensure the behaviour you get is the behaviour you intend. +the `s` and `d` flags to ensure the behaviour you get is the behaviour you intend. ```@docs @dna_str @@ -303,7 +306,7 @@ the bodies of things like for loops. And if you use them and are unsure, use the ``` ## Loose parsing -As of version 3.2.0, BioSequences.jl provide the [`bioseq`](@ref) function, which can be used to build a `LongSequence` +As of version 3.2.0, BioSequences.jl provides the [`bioseq`](@ref) function, which can be used to build a `LongSequence` from a string (or an `AbstractVector{UInt8}`) without knowing the correct `Alphabet`. ```jldoctest @@ -312,7 +315,9 @@ julia> bioseq("ATGTGCTGA") ATGTGCTGA ``` -The function will prioritise 2-bit alphabets over 4-bit alphabets, and prefer smaller alphabets (like `DNAAlphabet{4}`) over larger (like `AminoAcidAlphabet`). +The function will prioritise 2-bit alphabets over 4-bit alphabets, and prefer +smaller alphabets (like `DNAAlphabet{4}`) over larger (like +`AminoAcidAlphabet`). If the input cannot be encoded by any of the built-in alphabets, an error is thrown: ```jldoctest @@ -321,8 +326,9 @@ ERROR: cannot encode 0x30 (Char '0') in AminoAcidAlphabet [...] ``` -Note that this function is only intended to be used for interactive, ephemeral work. -The function is necessarily type unstable, and the precise returned alphabet for a given input is a heuristic which is subject to change. +Note that this function is only intended to be used for interactive, ephemeral +work. The function is necessarily type unstable, and the precise returned +alphabet for a given input is a heuristic that is subject to change. ```@docs bioseq @@ -331,7 +337,7 @@ guess_alphabet ## Comparison to other sequence types Following Base standards, BioSequences do not compare equal to other containers even if they have the same elements. -To e.g. compare a BioSequence with a vector of DNA, compare the elements themselves: +To compare a BioSequence with a vector of DNA, for example, compare the elements themselves: ```jldoctest julia> seq = dna"GAGCTGA"; vec = collect(seq); @@ -339,5 +345,5 @@ julia> seq == vec, isequal(seq, vec) (false, false) julia> length(seq) == length(vec) && all(i == j for (i, j) in zip(seq, vec)) -true +true ``` diff --git a/docs/src/counting.md b/docs/src/counting.md index 58448530f..039047c8f 100644 --- a/docs/src/counting.md +++ b/docs/src/counting.md @@ -6,8 +6,8 @@ end ``` # Counting -`BioSequences` contains functionality to efficiently count biosymbols in a biosequence -that satisfies some predicate. +`BioSequences` contains functionality to efficiently count biosymbols in a +biosequence that satisfies some predicate. Consider a naive counting function like this: ```julia @@ -17,18 +17,18 @@ function count_Ns(seq::BioSequence{<:DNAAlphabet}) ns += (i == DNA_N)::Bool end ns -end +end ``` This function can be more efficiently implemented by exploiting the internal data layout of certain biosequences. Therefore, Julia provides optimised methods for `Base.count`, such that `count_Ns` -above can be more efficiently expressed `count(==(DNA_N), seq)`. +above can be more efficiently expressed as `count(==(DNA_N), seq)`. !!! note It is important to understand that this speed is achieved with custom methods of `Base.count`, and not by a generic mechanism that improves the speed of counting - symbols in `BioSequence `in general. + symbols in a `BioSequence` in general. Hence, while `count(==(DNA_N), seq)` may be optimised, `count(i -> i == DNA_N, seq)` is not, as this is a different method. @@ -44,7 +44,7 @@ have optimised methods. ## Matches and mismatches The methods `matches` and `mismatches` take two -sequences and count the number of positions where the sequences are unequal or equal, respectively. +sequences and count the number of positions where the sequences are equal or unequal, respectively. They are equivalent to `matches(a, b) = count(splat(==), zip(a, b))` (and with `!=`, respectively). @@ -62,7 +62,7 @@ gc_content ``` ## Deprecated aliases -Several of the optimised `count` methods have function names, which are deprecated: +Several of the optimised `count` methods have function names that are deprecated: | Deprecated function | Instead use | | :------------------- | :------------------------ | diff --git a/docs/src/index.md b/docs/src/index.md index 17237890f..92f567073 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -11,14 +11,14 @@ biological sequences, including DNA, RNA, and amino acid sequences. ## Installation -You can install BioSequences from the julia -REPL. Press `]` to enter pkg mode again, and enter the following: +You can install BioSequences from the Julia +REPL. Press `]` to enter pkg mode, and enter the following: ```julia add BioSequences ``` -If you are interested in the cutting edge of the development, please check out +If you are interested in the cutting edge of development, please check out the master branch to try new features before release. @@ -39,5 +39,5 @@ detailed contributor and maintainer guidelines, and code of conduct. ## Questions? If you have a question about contributing or using BioJulia software, come -on over and chat to us on [the #biology channel on the Julia SLack](https://julialang.org/slack/), or you can try the +on over and chat to us on [the #biology channel on the Julia Slack](https://julialang.org/slack/), or you can try the [Bio category of the Julia discourse site](https://discourse.julialang.org/c/domain/bio). diff --git a/docs/src/interfaces.md b/docs/src/interfaces.md index 6e69da626..e9f583c26 100644 --- a/docs/src/interfaces.md +++ b/docs/src/interfaces.md @@ -5,25 +5,25 @@ DocTestSetup = quote end ``` -# Custom BioSequences types +# Custom BioSequence types -If you're a developing your own Bioinformatics package or method, you may find -that the reference implementation of concrete `LongSequence` types provided in +If you're developing your own bioinformatics package or method, you may find +that the reference implementations of concrete `LongSequence` types provided in this package are not optimal for your purposes. -This page describes the interfaces for BioSequences' core types for -developers or other packages implementing their own sequence types or extending +This page describes the interfaces for BioSequences' core types, for developers +and other packages implementing their own sequence types or extending BioSequences functionality. ## Implementing custom Alphabets -Recall the required methods that define the [`Alphabet`](@ref) interface. +Recall the required methods that define the [`Alphabet`](@ref) interface. -To create an example custom alphabet, we need to create a singleton type, that +To create an example custom alphabet, we need to create a singleton type that implements a few methods in order to conform to the interface as described in the [`Alphabet`](@ref) documentation. -Let's do that for a restricted Amino Acid alphabet. We can test that it conforms +Let's do that for a restricted amino acid alphabet. We can test that it conforms to the interface with the [`BioSequences.has_interface`](@ref) function. ```jldoctest @@ -67,14 +67,14 @@ true ## Implementing custom BioSequences -Recall the required methods that define the [`BioSequence`](@ref) interface. +Recall the required methods that define the [`BioSequence`](@ref) interface. -To create an example custom alphabet, we need to create a singleton type, that +To create an example custom sequence type, we need to create a singleton type that implements a few methods in order to conform to the interface as described in the [`BioSequence`](@ref) documentation. Let's do that for a custom sequence type that is optimised to represent a small -sequence: A Codon. We can test that it conforms to the interface with the +sequence: a codon. We can test that it conforms to the interface with the [`BioSequences.has_interface`](@ref) function. ```jldoctest diff --git a/docs/src/io.md b/docs/src/io.md index 489b42e0b..e826ecaa4 100644 --- a/docs/src/io.md +++ b/docs/src/io.md @@ -3,10 +3,10 @@ Versions of BioSequences prior to v2.0 provided a FASTA, FASTQ, and 2Bit submodule for working with formatted sequence files. -After version v2.0, in order to neatly separate concerns, these submodules were +After v2.0, in order to neatly separate concerns, these submodules were removed. -Instead there will now be dedicated BioJulia packages for each format. Each +Instead, there are now dedicated BioJulia packages for each format. Each of these will be compatible with BioSequences. A list of all of the different formats and packages is provided below to help @@ -16,4 +16,4 @@ you find them quickly. |:------ |:------------------------------------------------ | | FASTA | [FASTX.jl](https://github.com/BioJulia/FASTX.jl) | | FASTQ | [FASTX.jl](https://github.com/BioJulia/FASTX.jl) | -| 2Bit | [TwoBit.jl](https://github.com/BioJulia/TwoBit.jl) | \ No newline at end of file +| 2Bit | [TwoBit.jl](https://github.com/BioJulia/TwoBit.jl) | diff --git a/docs/src/random.md b/docs/src/random.md index 7abc04c3e..382778263 100644 --- a/docs/src/random.md +++ b/docs/src/random.md @@ -8,8 +8,8 @@ end ## Long sequences -You can generate random long sequences using the `randdna` function and the -`Sampler`'s implemented in BioSequences: +You can generate random long sequences using the `randdnaseq` function and the +samplers implemented in BioSequences: ```@docs randseq @@ -18,4 +18,4 @@ randrnaseq randaaseq SamplerUniform SamplerWeighted -``` \ No newline at end of file +``` diff --git a/docs/src/recipes.md b/docs/src/recipes.md index 7a89dca79..00dc33a64 100644 --- a/docs/src/recipes.md +++ b/docs/src/recipes.md @@ -63,4 +63,3 @@ one_hot(dna"TGNTKCTW-T") 1 0 1 1 1 0 1 1 0 1 ``` - diff --git a/docs/src/sequence_search.md b/docs/src/sequence_search.md index d2210eef2..61e274af0 100644 --- a/docs/src/sequence_search.md +++ b/docs/src/sequence_search.md @@ -10,21 +10,22 @@ end There are many ways to search for particular motifs in biological sequences: 1. Exact searches, where you are looking for exact matches of a particular - character of substring. + character or substring. 2. Approximate searches, where you are looking for sequences that are sufficiently similar to a given sequence or family of sequences. 3. Searches where you are looking for sequences that conform to some sort of pattern. -Like other Julia sequences such as `Vector`, you can search a `BioSequence` with the `findfirst(predicate, collection)` method pattern. +Like other Julia sequences such as `Vector`, you can search a `BioSequence` with +the `findfirst(predicate, collection)` method pattern. All these kinds of searches are provided in BioSequences.jl, and they all -conform to the `findnext`, `findprev`, and `occursin` patterns established in `Base` for -`String` and collections like `Vector`. +conform to the `findnext`, `findprev`, and `occursin` patterns established in +`Base` for `String` and collections like `Vector`. The exception is searching using the specialised -regex provided in this package, which as you shall see, conforms to the `match` -pattern established in `Base` for pcre and `String`s. +regex provided in this package, which, as you shall see, conforms to the `match` +pattern established in `Base` for PCREs and `String`s. ## Symbol search @@ -70,22 +71,22 @@ Query patterns can be described in regular expressions. The syntax supports a subset of Perl and PROSITE's notation. Biological regexes can be constructed using the `BioRegex` constructor, for -example by doing `BioRegex{AminoAcid}("MV+")`. For bioregex literals, it is -instead recommended using the `@biore_str` macro: +example by doing `BioRegex{AminoAcid}("MV+")`. For BioRegex literals, it is +instead recommended to use the `@biore_str` macro: The Perl-like syntax starts with `biore` (BIOlogical REgular expression) -and ends with a symbol option: "dna", "rna" or "aa". For example, `biore"A+"dna` +and ends with a symbol option: `dna`, `rna`, or `aa`. For example, `biore"A+"dna` is a regular expression for DNA sequences and `biore"A+"aa` is for amino acid -sequences. The symbol options can be abbreviated to its first character: "d", -"r" or "a", respectively. +sequences. The symbol options can be abbreviated to their first character: `d`, +`r`, or `a`, respectively. Here are examples of using the regular expression for `BioSequence`s: ```jldoctest julia> match(biore"A+C*"dna, dna"AAAACC") -RegexMatch("AAAACC") +BioSequences.RE.RegexMatch{LongSequence{DNAAlphabet{4}}}(AAAACC, [1, 7]) julia> match(biore"A+C*"d, dna"AAAACC") -RegexMatch("AAAACC") +BioSequences.RE.RegexMatch{LongSequence{DNAAlphabet{4}}}(AAAACC, [1, 7]) julia> occursin(biore"A+C*"dna, dna"AAC") true @@ -95,25 +96,25 @@ false ``` -`match` will return a `RegexMatch` if a match is found, otherwise it will return `nothing` if no match is found. +`match` returns a `RegexMatch` if a match is found, and `nothing` otherwise. The table below summarizes available syntax elements. | Syntax | Description | Example | |:------:|:------------|:--------| | `\|` | alternation | `"A\|T"` matches `"A"` and `"T"` | -| `*` | zero or more times repeat | `"TA*"` matches `"T"`, `"TA"` and `"TAA"` | -| `+` | one or more times repeat | `"TA+"` matches `"TA"` and `"TAA"` | -| `?` | zero or one time | `"TA?"` matches `"T"` and `"TA"` | -| `{n,}` | `n` or more times repeat | `"A{3,}"` matches `"AAA"` and `"AAAA"` | -| `{n,m}` | `n`-`m` times repeat | `"A{3,5}"` matches `"AAA"`, `"AAAA"` and `"AAAAA"`| +| `*` | zero or more repetitions | `"TA*"` matches `"T"`, `"TA"`, and `"TAA"` | +| `+` | one or more repetitions | `"TA+"` matches `"TA"` and `"TAA"` | +| `?` | zero or one repetition | `"TA?"` matches `"T"` and `"TA"` | +| `{n,}` | `n` or more repetitions | `"A{3,}"` matches `"AAA"` and `"AAAA"` | +| `{n,m}` | `n`–`m` repetitions | `"A{3,5}"` matches `"AAA"`, `"AAAA"`, and `"AAAAA"` | | `^` | the start of the sequence | `"^TAN*"` matches `"TATGT"` | | `$` | the end of the sequence | `"N*TA$"` matches `"GCTA"` | | `(...)` | pattern grouping | `"(TA)+"` matches `"TA"` and `"TATA"` | -| `[...]` | one of symbols | `"[ACG]+"` matches `"AGGC"` | +| `[...]` | one of the symbols | `"[ACG]+"` matches `"AGGC"` | -`eachmatch` and `findfirst` are also defined, just like usual regex and strings -found in `Base`. +`eachmatch` and `findfirst` are also defined, just like the usual regex and string +methods found in `Base`. ```jldoctest julia> collect(matched(x) for x in eachmatch(biore"TATA*?"d, dna"TATTATAATTA")) # overlap @@ -139,27 +140,27 @@ julia> findfirst(biore"TATA*"d, dna"TATTATAATTA", 2) Noteworthy differences from strings are: * Ambiguous characters match any compatible characters (e.g. `biore"N"d` is equivalent to `biore"[ACGT]"d`). -* Whitespaces are ignored (e.g. `biore"A C G"d` is equivalent to `biore"ACG"d`). +* Whitespace is ignored (e.g. `biore"A C G"d` is equivalent to `biore"ACG"d`). The PROSITE notation is described in [ScanProsite - user manual](https://prosite.expasy.org/scanprosite/scanprosite_doc.html). The syntax supports almost all notations including the extended syntax. The PROSITE -notation starts with `prosite` prefix and no symbol option is needed because it +notation starts with the `prosite` prefix, and no symbol option is needed because it always describes patterns of amino acid sequences: ```jldoctest julia> match(prosite"[AC]-x-V-x(4)-{ED}", aa"CPVPQARG") -RegexMatch("CPVPQARG") +BioSequences.RE.RegexMatch{LongAA}(CPVPQARG, [1, 9]) julia> match(prosite"[AC]xVx(4){ED}", aa"CPVPQARG") -RegexMatch("CPVPQARG") +BioSequences.RE.RegexMatch{LongAA}(CPVPQARG, [1, 9]) ``` ### Position weight matrix search -A motif can be specified using [position weight +A motif can be specified using a [position weight matrix](https://en.wikipedia.org/wiki/Position_weight_matrix) (PWM) in a probabilistic way. This method searches for the first position in the sequence where a score @@ -231,7 +232,7 @@ probability ``p(s)`` as follows ([Wasserman2004]): ``` However, if you just want to quickly conduct a search, constructing the PFM and -PWM is done for you as a convenience if you build a `PWMSearchQuery`, using a +PWM is done for you as a convenience when you build a `PWMSearchQuery` using a collection of sequences: ```jldoctest diff --git a/docs/src/symbols.md b/docs/src/symbols.md index f63b2aa74..7a5d2b228 100644 --- a/docs/src/symbols.md +++ b/docs/src/symbols.md @@ -7,7 +7,7 @@ end # Biological symbols The `BioSequences` module reexports the biological symbol (character) types -that are provided by [BioSymbols.jl:](https://biojulia.dev/BioSymbols.jl/stable) +that are provided by [BioSymbols.jl](https://biojulia.dev/BioSymbols.jl/stable): | Type | Meaning | | :-------------- | :------------- | @@ -21,7 +21,7 @@ elements of strings. ## DNA and RNA nucleotides -Set of nucleotide symbols in BioSequences covers IUPAC nucleotide base plus +The set of nucleotide symbols in BioSequences covers the IUPAC nucleotide bases plus a gap symbol: | Symbol | Constant | Meaning | @@ -81,7 +81,7 @@ true ## Amino acids -Set of amino acid symbols also covers IUPAC amino acid symbols plus a gap symbol: +The set of amino acid symbols also covers the IUPAC amino acid symbols plus a gap symbol: | Symbol | Constant | Meaning | | :----------- | :-------------- | :-------------------------- | diff --git a/docs/src/transforms.md b/docs/src/transforms.md index 6d829f081..6fb5ca42c 100644 --- a/docs/src/transforms.md +++ b/docs/src/transforms.md @@ -9,7 +9,7 @@ end ## Indexing -Most `BioSequence` concrete subtypes for the most part behave like other vector +Most `BioSequence` concrete subtypes behave like other vector or string types. They can be indexed using integers or ranges: For example, with `LongSequence`s: @@ -29,7 +29,7 @@ TANAGTNNAGTACC ``` The biological symbol at a given locus in a biological sequence can be set using -setindex: +`setindex!`: ```jldoctest julia> seq = dna"ACGTTTANAGTNNAGTACC" @@ -42,9 +42,8 @@ DNA_A ``` !!! note - Some types such can be indexed using integers but not using ranges. + Some types can be indexed using integers but not using ranges. - For `LongSequence` types, indexing a sequence by range creates a copy of the original sequence, similar to `Array` in Julia's `Base` library. If you find yourself slowed down by the allocation of these subsequences, consider using a sequence view @@ -52,9 +51,9 @@ instead. ## Modifying sequences -In addition to `setindex`, many other modifying operations are possible for +In addition to `setindex!`, many other modifying operations are possible for biological sequences such as `push!`, `pop!`, and `insert!`, which should be -familiar to anyone used to editing arrays. +familiar to anyone who has used arrays. ```@docs push!(::BioSequences.BioSequence, ::Any) @@ -134,19 +133,19 @@ ACGTAT ``` -Many of these methods also have a version which makes a copy of the input +Many of these methods also have a version that makes a copy of the input sequence, so you get a modified copy, and don't alter the original sequence. Such methods are named the same, but without the exclamation mark. -E.g. `reverse` instead of `reverse!`, and `ungap` instead of `ungap!`. +E.g. `reverse` instead of `reverse!`, and `ungap` instead of `ungap!`. #### Translation -Translation is a slightly more complex transformation for RNA Sequences and so +Translation is a slightly more complex transformation for RNA sequences, and so we describe it here in more detail. -The [`translate`](@ref) function translates a sequence of codons in a RNA sequence -to a amino acid sequence based on a genetic code. The `BioSequences` package -provides all NCBI defined genetic codes and they are registered in +The [`translate`](@ref) function translates a sequence of codons in an RNA sequence +to an amino acid sequence based on a genetic code. The `BioSequences` package +provides all NCBI-defined genetic codes, and they are registered in [`ncbi_trans_table`](@ref). ```@docs diff --git a/docs/src/types.md b/docs/src/types.md index aaa599df2..5f5088bff 100644 --- a/docs/src/types.md +++ b/docs/src/types.md @@ -7,18 +7,18 @@ end # Abstract Types -BioSequences exports an abstract `BioSequence` type, and several concrete sequence -types which inherit from it. +BioSequences exports an abstract `BioSequence` type and several concrete sequence +types that inherit from it. ## The abstract BioSequence BioSequences provides an abstract type called a `BioSequence{A<:Alphabet}`. -This abstract type, and the methods and traits is supports, allows for +This abstract type, and the methods and traits it supports, allow for many algorithms in BioSequences to be written as generically as possible, thus reducing the amount of code to read and understand, whilst maintaining high performance when such code is compiled for a concrete BioSequence subtype. Additionally, it allows new types to be implemented that are fully compatible -with the rest of BioSequences, providing that key methods or traits are defined). +with the rest of BioSequences, provided that key methods or traits are defined. ```@docs BioSequence @@ -31,22 +31,21 @@ NucSeq AASeq ``` -Let's have a closer look at some of those methods that a subtype of `BioSequence` -must implement. Check out julia base library docs for `length`, `copy` and `resize!`. +Let's have a closer look at some of the methods that a subtype of `BioSequence` +must implement. Check out the Julia Base library docs for `length`, `copy`, and `resize!`. -```@docs +```@docs encoded_data_eltype extract_encoded_element encoded_setindex! ``` -A correctly defined subtype of `BioSequence` that satisfies the interface, will -find the vast majority of methods described in the rest of this manual -should work out of the box for that type. But they can always be overloaded if -needed. Indeed the `LongSequence` type overloads Indeed some of the generic -`BioSequence` methods, are overloaded for `LongSequence`, for example -for transformation and counting operations where efficiency gains can be made -due to the specific internal representation of a specific type. +A correctly defined subtype of `BioSequence` that satisfies the interface should +have the vast majority of methods described in the rest of this manual work out +of the box for that type. But they can always be overloaded if needed. Indeed, +the `LongSequence` type overloads some of the generic `BioSequence` methods, for +example transformation and counting operations where efficiency gains can be +made due to the specific internal representation of the type. ## The abstract Alphabet @@ -78,16 +77,16 @@ LongSequence ## Sequence views -Similar to how Base Julia offers views of array objects, BioSequences offers view of -`LongSequence`s - the `LongSubSeq{A<:Alphabet}`. +Similar to how Julia's `Base` offers views of array objects, BioSequences offers views of +`LongSequence`s—the `LongSubSeq{A<:Alphabet}`. Conceptually, a `LongSubSeq{A}` is similar to a `LongSequence{A}`, but instead of storing -their own data, they refer to the data of a `LongSequence`. Modiying the `LongSequence` +its own data, it refers to the data of a `LongSequence`. Modifying the `LongSequence` will be reflected in the view, and vice versa. If the underlying `LongSequence` is truncated, the behaviour of a view is undefined. For the same reason, some operations are not supported for views, such as resizing. -The purpose of `LongSubSeq` is that, since they only contain a pointer to the -underlying array, an offset and a length, they are much lighter than `LongSequences`, -and will be stack allocated on Julia 1.5 and newer. Thus, the user may construct +The purpose of `LongSubSeq` is that, since it only contains a pointer to the +underlying array, an offset, and a length, it is much lighter than `LongSequence`s, +and will be stack-allocated on Julia 1.5 and newer. Thus, the user may construct millions of views without major performance implications. diff --git a/src/BioSequences.jl b/src/BioSequences.jl index 9c52b20a2..e6441df8a 100644 --- a/src/BioSequences.jl +++ b/src/BioSequences.jl @@ -259,19 +259,19 @@ Base.eltype(::Type{<:Search}) = Int Base.IteratorSize(::Type{<:Search}) = Base.SizeUnknown() """ - findall(pattern, sequence::BioSequence[,rng::UnitRange{Int}]; overlap::Bool=true)::Vector + findall(pattern, sequence::BioSequence[, rng::UnitRange{Int}]; overlap::Bool=true)::Vector Find all occurrences of `pattern` in `sequence`. -The return value is a vector of ranges of indices where the matching sequences were found. +The return value is a vector of ranges of indices where matching sequences were found. If there are no matching sequences, the return value is an empty vector. The search is restricted to the specified range when `rng` is set. -With the keyword argument `overlap` set as `true`, the start index for the next search gets set to the start of the current match plus one; if set to `false`, the start index for the next search gets set to the end of the current match plus one. +With the keyword argument `overlap` set to `true`, the start index for the next search is set to the start of the current match plus one; if set to `false`, the start index for the next search is set to the end of the current match plus one. The default value for the keyword argument `overlap` is `true`. -The `pattern` can be a `Biosymbol` or a search query. +The `pattern` can be a `BioSymbol` or a search query. See also [`ExactSearchQuery`](@ref), [`ApproximateSearchQuery`](@ref), [`PWMSearchQuery`](@ref). diff --git a/src/alphabet.jl b/src/alphabet.jl index 3175ac72a..b7ca88140 100644 --- a/src/alphabet.jl +++ b/src/alphabet.jl @@ -11,15 +11,15 @@ Alphabet `Alphabet` is the most important type trait for `BioSequence`. An `Alphabet` -represents a set of biological symbols encoded by a sequence, e.g. A, C, G -and T for a DNA Alphabet that requires only 2 bits to represent each symbol. +represents a set of biological symbols encoded by a sequence, e.g. A, C, G, +and T for a DNA alphabet that requires only 2 bits to represent each symbol. # Extended help -* Subtypes of Alphabet are singleton structs that may or may not be parameterized. -* Alphabets span over a *finite* set of biological symbols. -* The alphabet controls the encoding from some internal "encoded data" to a BioSymbol +* Subtypes of `Alphabet` are singleton structs that may or may not be parameterized. +* Alphabets span a *finite* set of biological symbols. +* The alphabet controls the encoding from some internal "encoded data" to a `BioSymbol` of the alphabet's element type, as well as the decoding, the inverse process. -* An `Alphabet`'s `encode` method must not produce invalid data. +* An `Alphabet`'s `encode` method must not produce invalid data. ### Required methods Every subtype `A` of `Alphabet` must implement: @@ -27,13 +27,13 @@ Every subtype `A` of `Alphabet` must implement: * `symbols(::A)::Tuple{Vararg{S}}`. This gives tuples of all symbols in the set of `A`. * `encode(::A, ::S)::E` encodes a symbol to an internal data eltype `E`. * `decode(::A, ::E)::S` decodes an internal data eltype `E` to a symbol `S`. -* Except for `eltype` which must follow Base conventions, all functions operating +* Except for `eltype`, which must follow `Base` conventions, all functions operating on `Alphabet` should operate on instances of the alphabet, not the type. If you want interoperation with existing subtypes of `BioSequence`, the encoded representation `E` must be of type `UInt`, and you must also implement: -* `BitsPerSymbol(::A)::BitsPerSymbol{N}`, where the `N` must be zero - or a power of two in [1, 2, 4, 8, 16, 32, [64 for 64-bit systems]]. +* `BitsPerSymbol(::A)::BitsPerSymbol{N}`, where `N` must be zero + or a power of two in [1, 2, 4, 8, 16, 32] (or 64 on 64-bit systems). ### Optional methods * `BitsPerSymbol` for compatibility with existing `BioSequence`s @@ -74,10 +74,10 @@ Base.length(A::Alphabet) = length(symbols(A)) """ BitsPerSymbol{N} -A trait object specifying the number of bits it takes to encode a biosymbol in an `Alphabet` -Alphabets `A` should implement `BitsPerSymbol(::A)`. +A trait object specifying the number of bits it takes to encode a biosymbol in an `Alphabet`. +An alphabet `A` should implement `BitsPerSymbol(::A)`. For compatibility with existing BioSequences, the number of bits should be a power of two -between 1 and 32, both inclusive. +between 1 and 32, inclusive. See also: [`Alphabet`](@ref) """ struct BitsPerSymbol{N} end @@ -91,10 +91,10 @@ iscomplete(A::Alphabet) = Val(length(symbols(A)) === 1 << bits_per_symbol(A)) """ encode(::Alphabet, s::BioSymbol) -Internal function, do not use in user code. -Encode BioSymbol `s` to an internal representation using an [`Alphabet`](@ref). -This decoding is checked to enforce valid data element. -If `s` cannot be encoded to the given alphabet, throw an `EncodeError` +Internal function; do not use it in user code. +Encode the `BioSymbol` `s` to an internal representation using an [`Alphabet`](@ref). +This encoding is checked to enforce valid data elements. +If `s` cannot be encoded to the given alphabet, throw an `EncodeError`. """ @inline function encode(A::Alphabet, s::BioSymbol) @@ -106,10 +106,10 @@ tryencode(A::Alphabet, s::Any) = nothing """ tryencode(::Alphabet, x::S) -Try encoding BioSymbol `S` to the internal representation of [`Alphabet`](@ref), +Try to encode a `BioSymbol` `x` to the internal representation of [`Alphabet`](@ref), returning `nothing` if not successful. -See also: `encode`[@ref], `decode`[@ref] +See also: [`encode`](@ref), [`decode`](@ref) """ function tryencode end @@ -162,8 +162,8 @@ end ## Nucleic acid alphabets """ -Alphabet of nucleic acids. Parameterized by the number of bits per symbol, by -default only `2` or `4`-bit variants exists. +Alphabet of nucleic acids. Parameterized by the number of bits per symbol; by +default, only `2`- or `4`-bit variants exist. """ abstract type NucleicAcidAlphabet{N} <: Alphabet end @@ -278,8 +278,8 @@ abstract type AlphabetCode end """ AsciiAlphabet -Trait for alphabet using ASCII characters as String representation. -Define `codetype(A) = AsciiAlphabet()` for a user-defined `Alphabet` A to gain speed. +Trait for an alphabet that uses ASCII characters as its string representation. +Define `codetype(A) = AsciiAlphabet()` for a user-defined `Alphabet` `A` to gain speed. Methods needed: `BioSymbols.stringbyte(::eltype(A))` and `ascii_encode(A, ::UInt8)`. """ struct AsciiAlphabet <: AlphabetCode end @@ -299,7 +299,7 @@ codetype(::Alphabet) = UnicodeAlphabet() Encode the ASCII character represented by `b` to the internal alphabet encoding. For example, the input byte `UInt8('C')` is encoded to `0x01` and `0x02` for -2- and 4-bit DNA alphabets, reprectively. +2- and 4-bit DNA alphabets, respectively. This method is only needed if the `Alphabet` is an `AsciiAlphabet`. See also: [`BioSequences.AsciiAlphabet`](@ref) @@ -346,10 +346,10 @@ end """ guess_alphabet(s::Union{AbstractString, AbstractVector{UInt8}}) -> Union{Integer, Alphabet} -Pick an `Alphabet` that can encode input `s`. If no `Alphabet` can, return the index of the first -byte of the input which is not encodable in any alphabet. +Pick an `Alphabet` that can encode input `s`. If no `Alphabet` can encode it, return the index of the first +byte in the input that is not encodable in any alphabet. This function only knows about the alphabets listed below. If multiple alphabets are possible, -pick the first from the order below (i.e. `DNAAlphabet{2}()` if possible, otherwise `RNAAlphabet{2}()` etc). +pick the first from the order below (i.e. `DNAAlphabet{2}()` if possible, otherwise `RNAAlphabet{2}()`, etc.). 1. `DNAAlphabet{2}()` 2. `RNAAlphabet{2}()` 3. `DNAAlphabet{4}()` diff --git a/src/biosequence/biosequence.jl b/src/biosequence/biosequence.jl index bb75dc684..85c0abcdb 100644 --- a/src/biosequence/biosequence.jl +++ b/src/biosequence/biosequence.jl @@ -15,21 +15,21 @@ and is parameterized by an `Alphabet`, which controls the element type. Its subtypes are characterized by: * Being a linear container type with random access and indices `Base.OneTo(length(x))`. * Containing zero or more internal data elements of type `encoded_data_eltype(typeof(x))`. -* Being associated with an `Alphabet`, `A` by being a subtype of `BioSequence{A}`. +* Being associated with an `Alphabet`, `A`, by being a subtype of `BioSequence{A}`. -A `BioSequence{A}` is indexed by an integer. The biosequence subtype, the index +A `BioSequence{A}` is indexed by an integer. The `BioSequence` subtype, the index and the alphabet `A` determine how to extract the internal encoded data. The alphabet decides how to decode the data to the element type of the biosequence. Hence, the element type and container type of a `BioSequence` are separated. -Subtypes `T` of `BioSequence` must implement the following, with `E` begin an +Subtypes `T` of `BioSequence` must implement the following, with `E` being an encoded data type: * `Base.length(::T)::Int` * `encoded_data_eltype(::Type{T})::Type{E}` * `extract_encoded_element(::T, ::Integer)::E` * `copy(::T)` -* T must be able to be constructed from any iterable with `length` defined and +* `T` must be able to be constructed from any iterable with `length` defined and with a known, compatible element type. Furthermore, mutable sequences should implement @@ -58,8 +58,8 @@ end """ has_interface(::Type{BioSequence}, ::T, syms::Vector, mutable::Bool, compat::Bool=true) -Check if type `T` conforms to the `BioSequence` interface. A `T` is constructed from the vector -of element types `syms` which must not be empty. +Check if type `T` conforms to the `BioSequence` interface. A value of type `T` is +constructed from the vector of elements `syms`, which must not be empty. If the `mutable` flag is set, also check the mutable interface. If the `compat` flag is set, check for compatibility with existing alphabets. """ @@ -123,7 +123,7 @@ end """ join!(seq::BioSequence, iter) -Concatenate all biosequences/biosymbols in `iter` into `seq`, resizing it to fit. +Concatenate all biosequences or biosymbols in `iter` into `seq`, resizing it to fit. # Examples ``` @@ -132,7 +132,7 @@ julia> join(LongDNA(), [dna"TAG", dna"AAC"]) TAGAAC ``` -see also [`join`](@ref) +See also [`join`](@ref). """ join!(seq::BioSequence, it) = _join!(seq, it, Val(false)) @@ -159,7 +159,7 @@ end """ join(::Type{T <: BioSequence}, seqs) -Concatenate all the biosequences/biosymbols in `seqs` to a biosequence of type `T`. +Concatenate all the biosequences or biosymbols in `seqs` to a biosequence of type `T`. # Examples ``` @@ -168,7 +168,7 @@ julia> join(LongDNA, [dna"TAG", dna"AAC"]) TAGAAC ``` -see also [`join!`](@ref) +See also [`join!`](@ref). """ function Base.join(::Type{T}, it::Union{Vector, Tuple, Set}) where {T <: BioSequence} _join!(T(undef, reduce((a, b) -> a + joinlen(b), it, init=0)), it, Val(true)) @@ -197,7 +197,7 @@ Returns the element type of the encoded data of the `BioSequence`. This is the return type of `extract_encoded_element`, i.e. the data type that stores the biological symbols in the biosequence. -See also: [`BioSequence`](@ref) +See also: [`BioSequence`](@ref) """ function encoded_data_eltype end @@ -208,7 +208,7 @@ Returns the encoded element at position `i`. This data can be decoded using `decode(A(), data)` to yield the element type of the biosequence. -See also: [`BioSequence`](@ref) +See also: [`BioSequence`](@ref) """ function extract_encoded_element end @@ -217,9 +217,9 @@ function extract_encoded_element end encoded_setindex!(seq::BioSequence, x::E, i::Integer) Given encoded data `x` of type `encoded_data_eltype(typeof(seq))`, -sets the internal sequence data at the given index. +set the internal sequence data at the given index. -See also: [`BioSequence`](@ref) +See also: [`BioSequence`](@ref) """ function encoded_setindex! end diff --git a/src/biosequence/copying.jl b/src/biosequence/copying.jl index 2985cd15b..126cd4e7f 100644 --- a/src/biosequence/copying.jl +++ b/src/biosequence/copying.jl @@ -14,12 +14,12 @@ function Base.copyto!(dst::BioSequence, src::BioSequence) end """ - copyto!(dst::BioSequence, soff, src::BioSequence, doff, N) + copyto!(dst::BioSequence, doff, src::BioSequence, soff, N) -In-place copy `N` elements from `src` starting at `soff` to `dst`, starting at `doff`. +Copy `N` elements in place from `src`, starting at `soff`, to `dst`, starting at `doff`. The length of `dst` must be greater than or equal to `N + doff - 1`. -The first N elements of `dst` are overwritten, -the other elements are left untouched. The alphabets of `src` and `dst` must be compatible. +The `N` elements of `dst` starting at `doff` are overwritten; the other elements +are left untouched. The alphabets of `src` and `dst` must be compatible. # Examples ``` @@ -44,4 +44,4 @@ function Base.copyto!(dst::BioSequence{A}, doff::Integer, end return dst -end \ No newline at end of file +end diff --git a/src/biosequence/predicates.jl b/src/biosequence/predicates.jl index c8b3cfeaa..489725063 100644 --- a/src/biosequence/predicates.jl +++ b/src/biosequence/predicates.jl @@ -118,7 +118,7 @@ end Returns `true` if `seq` is canonical. For any sequence, there is a reverse complement, which is the same sequence, but -on the complimentary strand of DNA: +on the complementary strand of DNA: ``` -------> @@ -128,11 +128,11 @@ CGATCGAT ``` !!! note - Using the [`reverse_complement`](@ref) of a DNA sequence will give give this + Using the [`reverse_complement`](@ref) of a DNA sequence will give this reverse complement. -Of the two sequences, the *canonical* of the two sequences is the lesser of the -two i.e. `canonical_seq < other_seq`. +Of the two sequences, the *canonical* one is the lesser of the two, i.e. +`canonical_seq < other_seq`. """ function iscanonical(seq::NucleotideSeq) i = 1 diff --git a/src/biosequence/transformations.jl b/src/biosequence/transformations.jl index 35547cae6..96e8f826c 100644 --- a/src/biosequence/transformations.jl +++ b/src/biosequence/transformations.jl @@ -9,7 +9,7 @@ """ empty!(seq::BioSequence) -Completely empty a biological sequence `seq` of nucleotides. +Completely empty the biological sequence `seq`. """ Base.empty!(seq::BioSequence) = resize!(seq, 0) @@ -43,7 +43,7 @@ end """ insert!(seq::BioSequence, i, x) -Insert a biological symbol `x` into a biological sequence `seq`, at the given +Insert a biological symbol `x` into a biological sequence `seq` at the given index `i`. Returns the mutated `seq`. # Examples @@ -69,9 +69,9 @@ end """ spliceinto!(seq::BioSequence, i::Integer, x) -Insert the sequence `x` into a biological sequence `seq`, at the given index `i`. +Insert the sequence `x` into a biological sequence `seq` at the given index `i`. After splicing, the `seq`'s symbols at indices `i:i+length(x)-1` are equal to `x`, -and the the symbols that were previously there are moved to the right. +and the symbols that were previously there are moved to the right. # Examples ```jldoctest @@ -95,8 +95,8 @@ end """ spliceinto!(seq::BioSequence, span::UnitRange, x) -Delete the symbols at indices `span` in `seq`, and then copy `x` into the -first deleted position, then return `seq`. +Delete the symbols at indices `span` in `seq`, then copy `x` into the +first deleted position, and return `seq`. This is equivalent to `deleteat!(seq, span); spliceinto!(seq, first(span), x)`, but is more efficient. @@ -138,7 +138,7 @@ end """ deleteat!(seq::BioSequence, range::UnitRange{<:Integer}) -Deletes a defined `range` from a biological sequence `seq`. +Delete the specified `range` from a biological sequence `seq`. Modifies the input sequence. """ @@ -252,7 +252,7 @@ end """ complement(seq) -Make a complement sequence of `seq`. +Create a complement sequence of `seq`. """ function BioSymbols.complement(seq::NucleotideSeq) return complement!(copy(seq)) @@ -263,7 +263,7 @@ complement!(seq::NucleotideSeq) = map!(complement, seq) """ reverse_complement!(seq) -Make a reversed complement sequence of `seq` in place. +Create a reverse-complement sequence of `seq` in place. """ function reverse_complement!(seq::NucleotideSeq) return complement!(reverse!(seq)) @@ -272,7 +272,7 @@ end """ reverse_complement(seq) -Make a reversed complement sequence of `seq`. +Create a reverse-complement sequence of `seq`. """ function reverse_complement(seq::NucleotideSeq) return complement!(reverse(seq)) @@ -281,11 +281,11 @@ end """ canonical!(seq::NucleotideSeq) -Transforms the `seq` into its canonical form, if it is not already canonical. -Modifies the input sequence inplace. +Transform the `seq` into its canonical form if it is not already canonical. +Modify the input sequence in place. For any sequence, there is a reverse complement, which is the same sequence, but -on the complimentary strand of DNA: +on the complementary strand of DNA: ``` -------> @@ -295,13 +295,13 @@ CGATCGAT ``` !!! note - Using the [`reverse_complement`](@ref) of a DNA sequence will give give this + Using the [`reverse_complement`](@ref) of a DNA sequence will give this reverse complement. -Of the two sequences, the *canonical* of the two sequences is the lesser of the -two i.e. `canonical_seq < other_seq`. +Of the two sequences, the *canonical* one is the lesser of the two, i.e. +`canonical_seq < other_seq`. -Using this function on a `seq` will ensure it is the canonical version. +Using this function on `seq` will ensure it is the canonical version. """ function canonical!(seq::NucleotideSeq) if !iscanonical(seq) diff --git a/src/bit-manipulation/bitindex.jl b/src/bit-manipulation/bitindex.jl index 406de72bf..a332763af 100644 --- a/src/bit-manipulation/bitindex.jl +++ b/src/bit-manipulation/bitindex.jl @@ -14,7 +14,7 @@ """ BitIndex -`BitIndex` is an internal type used in BioSequences.It contains +`BitIndex` is an internal type used in BioSequences. It contains a bit offset. For biosequences with an internal array of coding units, it can be used to obtain the array index and element bit offset. diff --git a/src/counting.jl b/src/counting.jl index 5ccb11e9f..227dffbce 100644 --- a/src/counting.jl +++ b/src/counting.jl @@ -26,8 +26,8 @@ end """ gc_content(seq::BioSequence) -> Float64 -Calculate GC content of `seq`, i.e. the number of symbols that is `DNA_C`, `DNA_G`, -`DNA_C` or `DNA_G` divided by the length of the sequence. +Calculate GC content of `seq`, i.e. the number of symbols that are `DNA_C` or `DNA_G` +divided by the length of the sequence. # Examples ```jldoctest @@ -43,9 +43,9 @@ Base.count(::typeof(isGC), seq::NucleotideSeq) = _n_gc(seq) # Aliases """ - matches(a::BioSequence, b::BioSequences) -> Int + matches(a::BioSequence, b::BioSequence) -> Int -Count the number of positions in where `a` and `b` are equal. +Count the number of positions where `a` and `b` are equal. If `b` is given, and the length of `a` and `b` differ, look only at the indices of the shorter sequence. This function does not provide any special handling of ambiguous symbols, @@ -73,9 +73,9 @@ function matches end ) """ - mismatches(a::BioSequence, b::BioSequences) -> Int + mismatches(a::BioSequence, b::BioSequence) -> Int -Count the number of positions in where `a` and `b` differ. +Count the number of positions where `a` and `b` differ. If `b` is given, and the length of `a` and `b` differ, look only at the indices of the shorter sequence. This function does not provide any special handling of ambiguous symbols, @@ -111,7 +111,7 @@ of the shorter sequence. !!! warning Passing in two sequences is deprecated. In a future, breaking release of - BioSequences, this will throw a `MethodError` + BioSequences, this will throw a `MethodError`. # Examples ```jldoctest @@ -135,14 +135,14 @@ Base.count(::typeof(isgap), seq::BioSequence) = _n_gaps(seq) """ n_ambiguous(a::BioSequence, [b::BioSequence]) -> Int -Count the number of positions where `a` (or `b`, if present) have ambigious symbols. +Count the number of positions where `a` (or `b`, if present) have ambiguous symbols. If `b` is given, and the length of `a` and `b` differ, look only at the indices of the shorter sequence. -Gaps are not ambigous. +Gaps are not ambiguous. !!! warning Passing in two sequences is deprecated. In a future, breaking release of - BioSequences, this will throw a `MethodError` + BioSequences, this will throw a `MethodError`. # Examples @@ -167,7 +167,7 @@ Base.count(::typeof(isambiguous), seq::BioSequence) = _n_ambiguous(seq) """ n_certain(a::BioSequence, [b::BioSequence]) -> Int -Count the number of positions where `a` (and `b`, if present) have certain (i.e. non-ambigous +Count the number of positions where `a` (and `b`, if present) have certain (i.e. non-ambiguous and non-gap) symbols. If `b` is given, and the length of `a` and `b` differ, look only at the indices of the shorter sequence. @@ -175,7 +175,7 @@ Gaps are not certain. !!! warning Passing in two sequences is deprecated. In a future, breaking release of - BioSequences, this will throw a `MethodError` + BioSequences, this will throw a `MethodError`. # Examples @@ -395,4 +395,3 @@ for (sym, op, joiner) in [ @eval(@deprecate($(sym)(a::$(T)), count($(op), a))) end end - diff --git a/src/longsequences/copying.jl b/src/longsequences/copying.jl index b98e03606..3d1fcc339 100644 --- a/src/longsequences/copying.jl +++ b/src/longsequences/copying.jl @@ -10,7 +10,7 @@ """ copy!(dst::LongSequence, src::BioSequence) -In-place copy content of `src` to `dst`, resizing `dst` to fit. +Copy the content of `src` to `dst` in place, resizing `dst` to fit. The alphabets of `src` and `dst` must be compatible. # Examples @@ -220,12 +220,12 @@ function Base.copyto!(dst::SeqOrView{<:Alphabet}, src::SeqLike, C::AlphabetCode) end """ - copyto!(dst::LongSequence, soff, src, doff, N) + copyto!(dst::LongSequence, doff, src, soff, N) -In-place encode `N` elements from `src` starting at `soff` to `dst`, starting at `doff`. +Encode `N` elements in place from `src`, starting at `soff`, to `dst`, starting at `doff`. The length of `dst` must be greater than or equal to `N + doff - 1`. -The first N elements of `dst` are overwritten, -the other elements are left untouched. The content of `src` must be able to be encoded to +The `N` elements of `dst` starting at `doff` are overwritten; the other elements +are left untouched. The content of `src` must be able to be encoded to the alphabet of `dst`. # Examples diff --git a/src/longsequences/longsequence.jl b/src/longsequences/longsequence.jl index f46cd3bbe..cf3319d11 100644 --- a/src/longsequences/longsequence.jl +++ b/src/longsequences/longsequence.jl @@ -66,21 +66,21 @@ aliases for convenience. | `LongSequence{RNAAlphabet{N}}` | `RNA` | `LongRNA{N}` | | `LongSequence{AminoAcidAlphabet}` | `AminoAcid` | `LongAA` | -The `LongDNA` and `LongRNA` aliases use a DNAAlphabet{4}. +The `LongDNA` and `LongRNA` aliases use a `DNAAlphabet{4}`. `DNAAlphabet{4}` permits ambiguous nucleotides, and a sequence must use at least 4 bits to internally store each element (and indeed `LongSequence` does). If you are sure that you are working with sequences with no ambiguous nucleotides, -you can use `LongSequences` parameterised with `DNAAlphabet{2}` instead. +you can use `LongSequence` parameterised with `DNAAlphabet{2}` instead. `DNAAlphabet{2}` is an alphabet that uses two bits per base and limits to only unambiguous nucleotide symbols (A,C,G,T). -Changing this single parameter, is all you need to do in order to benefit from memory savings. +Changing this single parameter is all you need to do in order to benefit from memory savings. Some computations that use bitwise operations will also be dramatically faster. -The same applies with `LongSequence{RNAAlphabet{4}}`, simply replace the alphabet +The same applies to `LongSequence{RNAAlphabet{4}}`; simply replace the alphabet parameter with `RNAAlphabet{2}` in order to benefit. """ mutable struct LongSequence{A <: Alphabet} <: BioSequence{A} @@ -92,7 +92,7 @@ mutable struct LongSequence{A <: Alphabet} <: BioSequence{A} end end -"An alias for LongSequence{<:NucleicAcidAlpabet{N}}" +"An alias for LongSequence{<:NucleicAcidAlphabet{N}}" const LongNuc{N} = LongSequence{<:NucleicAcidAlphabet{N}} "An alias for LongSequence{DNAAlphabet{N}}" @@ -125,4 +125,4 @@ const FourBit = SeqOrView{<:NucleicAcidAlphabet{4}} const TwoBit = SeqOrView{<:NucleicAcidAlphabet{2}} const TwoBitSeq = BioSequence{<:NucleicAcidAlphabet{2}} -const FourBitSeq = BioSequence{<:NucleicAcidAlphabet{4}} \ No newline at end of file +const FourBitSeq = BioSequence{<:NucleicAcidAlphabet{4}} diff --git a/src/longsequences/operators.jl b/src/longsequences/operators.jl index d40c4c883..0d30713b3 100644 --- a/src/longsequences/operators.jl +++ b/src/longsequences/operators.jl @@ -7,17 +7,17 @@ Construct a matrix of nucleotides or amino acids from a vector of `BioSequence`s. -If parameter `major` is set to `:site`, the matrix is created such that one -nucleotide from each sequence is placed in each column i.e. the matrix is laid +If the parameter `major` is set to `:site`, the matrix is created such that one +nucleotide from each sequence is placed in each column, i.e. the matrix is laid out in site-major order. This means that iteration over one position of many sequences is efficient, -as julia arrays are laid out in column major order. +as Julia arrays are laid out in column-major order. If the parameter `major` is set to `:seq`, the matrix is created such that each -sequence is placed in one column i.e. the matrix is laid out in sequence-major +sequence is placed in one column, i.e. the matrix is laid out in sequence-major order. This means that iteration across each sequence in turn is efficient, -as julia arrays are laid out in column major order. +as Julia arrays are laid out in column-major order. # Examples ```julia @@ -75,17 +75,17 @@ end Construct a matrix of `T` from a vector of `BioSequence`s. -If parameter `major` is set to `:site`, the matrix is created such that one -nucleotide from each sequence is placed in each column i.e. the matrix is laid +If the parameter `major` is set to `:site`, the matrix is created such that one +nucleotide from each sequence is placed in each column, i.e. the matrix is laid out in site-major order. This means that iteration over one position of many sequences is efficient, -as julia arrays are laid out in column major order. +as Julia arrays are laid out in column-major order. If the parameter `major` is set to `:seq`, the matrix is created such that each -sequence is placed in one column i.e. the matrix is laid out in sequence-major +sequence is placed in one column, i.e. the matrix is laid out in sequence-major order. This means that iteration across each sequence in turn is efficient, -as julia arrays are laid out in column major order. +as Julia arrays are laid out in column-major order. # Examples ```julia @@ -146,13 +146,13 @@ end Construct a sequence that is a consensus of a vector of sequences. -The consensus is established by a simple majority vote rule, where ambiguous +The consensus is established by a simple majority-vote rule, where ambiguous nucleotides cast an equal vote for each of their possible states. -For each site a winner(s) out of A, T(U), C, or G is determined, in the cases -of ties the ambiguity symbol that unifies all the winners is returned. -E.g if A and T tie, then W is inserted in the consensus. If all A, T, C, and G -tie at a site, then N is inserted in the consensus. Note this means that if a -nucletide e.g. 'C' and a gap '-' draw, the nucleotide will always win over the +For each site, one or more winners from A, T(U), C, or G are determined; in cases +of ties, the ambiguity symbol that unifies all the winners is returned. +E.g., if A and T tie, then W is inserted in the consensus. If all A, T, C, and G +tie at a site, then N is inserted in the consensus. Note that this means that if a +nucleotide (e.g. 'C') and a gap ('-') draw, the nucleotide will always win over the gap, even though they tied. # Examples diff --git a/src/longsequences/randseq.jl b/src/longsequences/randseq.jl index 8938508bb..36515b183 100644 --- a/src/longsequences/randseq.jl +++ b/src/longsequences/randseq.jl @@ -12,7 +12,7 @@ import Random: Sampler, rand!, default_rng """ SamplerUniform{T} -Uniform sampler of type T. Instantiate with a collection of eltype T containing +Uniform sampler of type `T`. Instantiate with a collection with eltype `T` containing the elements to sample. # Examples @@ -41,8 +41,8 @@ const DefaultAASampler = SamplerUniform(aa"ACDEFGHIKLMNPQRSTVWY") """ SamplerWeighted{T} -Weighted sampler of type T. Instantiate with a collection of eltype T containing -the elements to sample, and an orderen collection of probabilities to sample +Weighted sampler of type `T`. Instantiate with a collection with eltype `T` containing +the elements to sample and an ordered collection of probabilities to sample each element except the last. The last probability is the remaining probability up to 1. @@ -106,17 +106,17 @@ randaaseq(len::Integer) = randaaseq(default_rng(), len) """ randseq([rng::AbstractRNG], A::Alphabet, len::Integer) -Generate a LongSequence{A} of length `len` from the specified alphabet, drawn +Generate a `LongSequence{A}` of length `len` from the specified alphabet, drawn from the default distribution. User-defined alphabets should implement this -method to implement random LongSequence generation. +method for random `LongSequence` generation. For RNA and DNA alphabets, the default distribution is uniform across A, C, G, and T/U. For AminoAcidAlphabet, it is uniform across the 20 standard amino acids. -For a user-defined alphabet A, default is uniform across all elements of +For a user-defined alphabet `A`, the default is uniform across all elements of `symbols(A)`. -# Example: +# Examples ``` julia> seq = randseq(AminoAcidAlphabet(), 50) 50aa Amino Acid Sequence: @@ -133,7 +133,7 @@ end Generate a LongSequence{A} of length `len` with elements drawn from the given sampler. -# Example: +# Examples ``` # Generate 1000-length RNA with 4% chance of N, 24% for A, C, G, or U julia> sp = SamplerWeighted(rna"ACGUN", fill(0.24, 4)) @@ -199,23 +199,23 @@ end """ randdnaseq([rng::AbstractRNG], len::Integer) -Generate a random LongSequence{DNAAlphabet{4}} sequence of length `len`, -with bases sampled uniformly from [A, C, G, T] +Generate a random `LongSequence{DNAAlphabet{4}}` sequence of length `len`, +with bases sampled uniformly from [`A`, `C`, `G`, `T`]. """ randdnaseq(rng::AbstractRNG, len::Integer) = randseq(rng, DNAAlphabet{4}(), len) """ randrnaseq([rng::AbstractRNG], len::Integer) -Generate a random LongSequence{RNAAlphabet{4}} sequence of length `len`, -with bases sampled uniformly from [A, C, G, U] +Generate a random `LongSequence{RNAAlphabet{4}}` sequence of length `len`, +with bases sampled uniformly from [`A`, `C`, `G`, `U`]. """ randrnaseq(rng::AbstractRNG, len::Integer) = randseq(rng, RNAAlphabet{4}(), len) """ randaaseq([rng::AbstractRNG], len::Integer) -Generate a random LongSequence{AminoAcidAlphabet} sequence of length `len`, +Generate a random `LongSequence{AminoAcidAlphabet}` sequence of length `len`, with amino acids sampled uniformly from the 20 standard amino acids. """ randaaseq(rng::AbstractRNG, len::Integer) = randseq(rng, AminoAcidAlphabet(), len) diff --git a/src/longsequences/stringliterals.jl b/src/longsequences/stringliterals.jl index 12624d822..03debb1b3 100644 --- a/src/longsequences/stringliterals.jl +++ b/src/longsequences/stringliterals.jl @@ -15,26 +15,26 @@ remove_newlines(s) = replace(s, r"\r|\n" => "") Create a `LongDNA{4}` sequence at parse time from string `seq`. If `flag` is `"s"` ('static', the default), the sequence is created at parse time, and inserted directly into the returned expression. -A static string ought not to be mutated +A static sequence ought not to be mutated. Alternatively, if `flag` is `"d"` (dynamic), a new sequence is parsed and created -whenever the code where is macro is placed is run. +whenever the code where the macro is placed is run. See also: [`@aa_str`](@ref), [`@rna_str`](@ref) # Examples In the example below, the static sequence is created once, at parse time, NOT -when the function `f` is run. This means it is the _same_ sequence that is -pushed to repeatedly. +when the function `f` is run. This means it is the _same_ sequence that is +repeatedly appended to. ```jldoctest julia> f() = dna"TAG"; -julia> string(push!(f(), DNA_A)) # NB: Mutates static string! +julia> string(push!(f(), DNA_A)) # NB: Mutates static sequence! "TAGA" julia> string(push!(f(), DNA_A)) "TAGAA" -julia> f() = dna"TAG"d; # dynamically make seq +julia> f() = dna"TAG"d; # dynamically make a sequence julia> string(push!(f(), DNA_A)) "TAGA" diff --git a/src/longsequences/transformations.jl b/src/longsequences/transformations.jl index e2ede26d2..3dcb891af 100644 --- a/src/longsequences/transformations.jl +++ b/src/longsequences/transformations.jl @@ -5,11 +5,11 @@ """ resize!(seq, size, [force::Bool=false]) -Resize a biological sequence `seq`, to a given `size`. Does not resize the underlying data -array unless the new size does not fit. If `force`, always resize underlying data array. +Resize a biological sequence `seq` to a given `size`. Does not resize the underlying data +array unless the new size does not fit. If `force`, always resize the underlying data array. -Note that resizing to a larger size, and then loading from uninitialized positions -is not allowed and may cause undefined behaviour. +Note that resizing to a larger size and then loading from uninitialized positions +is not allowed and may cause undefined behaviour. Make sure to always fill any uninitialized biosymbols after resizing. """ function Base.resize!(seq::LongSequence{A}, size::Integer, force::Bool=false) where {A} @@ -115,7 +115,7 @@ end """ complement!(seq) -Make a complement sequence of `seq` in place. +Complement the biological sequence `seq` in place. """ function complement!(seq::LongSequence{A}) where {A<:NucleicAcidAlphabet} seqdata = seq.data diff --git a/src/search/ApproxSearchQuery.jl b/src/search/ApproxSearchQuery.jl index 44d226b47..3e4e5b8dd 100644 --- a/src/search/ApproxSearchQuery.jl +++ b/src/search/ApproxSearchQuery.jl @@ -15,7 +15,7 @@ These queries are used as a predicate for the `Base.findnext`, `Base.findprev`, `Base.occursin`, `Base.findfirst`, and `Base.findlast` functions. Using these functions with these queries allows you to search a given sequence -for a sub-sequence, whilst allowing a specific number of errors. +for a subsequence, whilst allowing a specific number of errors. In other words they find a subsequence of the target sequence within a specific [Levenshtein distance](https://en.wikipedia.org/wiki/Levenshtein_distance) of the @@ -75,10 +75,10 @@ end """ ApproximateSearchQuery(pat::S, comparator::F = isequal) where {F<:Function,S<:BioSequence} -Construct an [`ApproximateSearchQuery`](@ref) predicate for use with Base find functions. +Construct an [`ApproximateSearchQuery`](@ref) predicate for use with `Base` find functions. # Arguments -- `pat`: A concrete BioSequence that is the sub-sequence you want to search for. +- `pat`: A concrete `BioSequence` that is the subsequence you want to search for. - `comparator`: A function used to compare the symbols between sequences. `isequal` by default. """ function ApproximateSearchQuery(pat::S, comparator::F = isequal) where {F<:Function,S<:BioSequence} @@ -107,7 +107,7 @@ end """ findnext(query, k, seq, start) -Return the range of the first occurrence of `pat` in `seq[start:stop]` allowing +Return the range of the first occurrence of `query` in `seq[start:stop]` allowing up to `k` errors. Symbol comparison is done using the predicate supplied to the query. diff --git a/src/search/ExactSearchQuery.jl b/src/search/ExactSearchQuery.jl index 2c1875d30..06a3834cf 100644 --- a/src/search/ExactSearchQuery.jl +++ b/src/search/ExactSearchQuery.jl @@ -3,8 +3,8 @@ Query type for exact sequence search. -An exact search, is one where are you are looking in some given sequence, for -exact instances of some given substring. +An exact search is one where you are looking in a given sequence for exact +instances of a given substring. These queries are used as a predicate for the `Base.findnext`, `Base.findprev`, `Base.occursin`, `Base.findfirst`, and `Base.findlast` functions. @@ -69,10 +69,10 @@ end """ ExactSearchQuery(pat::BioSequence, comparator::Function = isequal) -Construct an [`ExactSearchQuery`](@ref) predicate for use with Base find functions. +Construct an [`ExactSearchQuery`](@ref) predicate for use with `Base` find functions. # Arguments -- `pat`: A concrete BioSequence that is the sub-sequence you want to search for. +- `pat`: A concrete `BioSequence` that is the subsequence you want to search for. - `comparator`: A function used to compare the symbols between sequences. `isequal` by default. """ function ExactSearchQuery(pat::BioSequence, comparator::Function = isequal) @@ -221,7 +221,7 @@ end """ findnext(query::ExactSearchQuery, seq::BioSequence, start::Integer) -Return the index of the first occurrence of `query` in `seq`. +Return the range of the first occurrence of `query` in `seq`. Symbol comparison is done using the predicate supplied to the query. By default, `ExactSearchQuery`'s predicate is `isequal`. @@ -240,7 +240,7 @@ Base.findfirst(pat::ExactSearchQuery, seq::BioSequence) = findnext(pat, seq, fir """ findprev(query::ExactSearchQuery, seq::BioSequence, start::Integer) -Return the index of the last occurrence of `query` in `seq`. +Return the range of the last occurrence of `query` in `seq`. Symbol comparison is done using the predicate supplied to the query. By default, `ExactSearchQuery`'s predicate is `isequal`. @@ -259,6 +259,6 @@ Base.findlast(query::ExactSearchQuery, seq::BioSequence) = findprev(query, seq, """ occursin(x::ExactSearchQuery, y::BioSequence) -Return Bool indicating presence of exact match of x in y. +Return a `Bool` indicating the presence of an exact match of `x` in `y`. """ Base.occursin(x::ExactSearchQuery, y::BioSequence) = quicksearch(x, y, 1, lastindex(y)) != 0 diff --git a/src/search/pwm.jl b/src/search/pwm.jl index b15eefaa7..1b0c07334 100644 --- a/src/search/pwm.jl +++ b/src/search/pwm.jl @@ -140,7 +140,7 @@ end Create a position weight matrix from a position frequency matrix `pfm`. -The positive weight matrix will be `log2.((pfm ./ sum(pfm, 1)) ./ prior)`. +The position weight matrix will be `log2.((pfm ./ sum(pfm, 1)) ./ prior)`. """ function PWM(pfm::PFM{S}; prior = fill(1/4, 4)) where S <: Union{DNA,RNA} if !all(x -> x > 0, prior) @@ -388,4 +388,4 @@ end function Base.findlast(query::PWMSearchQuery, seq::BioSequence) return findprev(query, seq, lastindex(seq)) -end \ No newline at end of file +end diff --git a/src/search/re.jl b/src/search/re.jl index ef17e3e75..02f2980cb 100644 --- a/src/search/re.jl +++ b/src/search/re.jl @@ -581,9 +581,9 @@ end """ Regex{T}(pattern::AbstractString, syntax=:pcre) -Biological regular expression to seatch for `pattern` in sequences of type `T`, where -`T` can be `DNA`, `RNA`, and `AminoAcid`. `syntax` can be `:pcre` or `:prosite` for AminoAcid -acids. +Biological regular expression to search for `pattern` in sequences of type `T`, where +`T` can be `DNA`, `RNA`, or `AminoAcid`. `syntax` can be `:pcre` or `:prosite` for amino +acid sequences. """ struct Regex{T <: Union{BioSequences.DNA, BioSequences.RNA, BioSequences.AminoAcid}} pat::String # regular expression pattern (for printing) @@ -632,27 +632,19 @@ end Result of matching by `Regex`. +# Examples +```jldoctest julia> match(biore"A(C[TG])+N(CA)"d, dna"ACGACA") -RegexMatch("ACGACA", 1="CG", 2="", 3="CA") +BioSequences.RE.RegexMatch{LongSequence{DNAAlphabet{4}}}(ACGACA, [1, 7, 2, 4, 4, 0, 5, 7]) +``` """ struct RegexMatch{S} seq::S captured::Vector{Int} end -function Base.show(io::IO, m::RegexMatch) - print(io, "RegexMatch(") - for k in 1:div(length(m.captured), 2) - if k > 1 - print(io, ", ", k - 1, '=') - end - print(io, '"', m.seq[m.captured[2k-1]:m.captured[2k]-1], '"') - end - print(io, ')') -end - """ - matched(match::BioRegexMatch) + matched(match::RegexMatch) Return the matched pattern as a `BioSequence`. """ @@ -662,13 +654,13 @@ end """ - captured(match::BioRegexMatch) + captured(match::RegexMatch) Return a vector of the captured patterns, where a pattern not captured is `nothing`. -# Examples: +# Examples ``` -julia> captured(biore"A(C[TG])+N"d, dna"ACGAA"") +julia> captured(biore"A(C[TG])+N"d, dna"ACGAA") 2-element Vector{Union{Nothing, LongDNA{4}, LongNuc{4, DNAAlphabet{4}}}}: CG nothing diff --git a/test/search/regex.jl b/test/search/regex.jl index eec0c375e..5e65f333f 100644 --- a/test/search/regex.jl +++ b/test/search/regex.jl @@ -19,6 +19,10 @@ @test string(biore"A+"rna) == "biore\"A+\"rna" @test string(biore"A+"aa) == "biore\"A+\"aa" + m = match(biore"A+C*"dna, dna"AAAACC") + @test sprint(show, m) == sprint(Base.show_default, m) + @test startswith(sprint(show, m), "BioSequences.RE.RegexMatch") + @test occursin(biore"A"d, dna"A") @test !occursin(biore"A"d, dna"C") @test !occursin(biore"A"d, dna"G") From bf8b6ca1e29eadbdbf5273f40b9d96aa98eed38e Mon Sep 17 00:00:00 2001 From: Jakob Nybo Andersen Date: Sat, 29 Aug 2026 15:13:24 +0200 Subject: [PATCH 2/9] Fix: Boundscheck user-exposed getindex for GeneticCode GeneticCode is an AbstractDict{UInt64, AminoAcid}, and exposes a getindex function which does not boundscheck, thereby exposing UB. I think this issue is a smell of a deeper design issue with GeneticCode: If it is an AbstractDict, why is the key an undocumented and internal encoding? Why is GeneticDict private? How can users make their own? Here, do a bandaid solution. We fix it properly in a breaking release. --- src/geneticcode.jl | 38 ++++++++++++++++++++++---------------- test/translation.jl | 2 ++ 2 files changed, 24 insertions(+), 16 deletions(-) diff --git a/src/geneticcode.jl b/src/geneticcode.jl index 2695b8ca7..a6d3d6fd1 100644 --- a/src/geneticcode.jl +++ b/src/geneticcode.jl @@ -16,6 +16,8 @@ function unambiguous_codon(a::XNA, b::XNA, c::XNA) end # A genetic code is a table mapping RNA 3-mers (i.e. RNAKmer{3}) to AminoAcids. +# TODO: In the next breaking release, stop making GeneticCode an +# AbstractDict{UInt64, AminoAcid}; UInt64 codon encodings are an internal detail. "Type representing a Genetic Code" struct GeneticCode <: AbstractDict{UInt64, AminoAcid} name::String @@ -26,8 +28,13 @@ end ### Basic Functions ### +@inline function getindex_encoding(code::GeneticCode, encoding::UInt64) + return @inbounds code.tbl[encoding + one(UInt64)] +end + function Base.getindex(code::GeneticCode, codon::UInt64) - return @inbounds code.tbl[codon + one(UInt64)] + codon < 64 || throw(KeyError(codon)) + return getindex_encoding(code, codon) end Base.copy(code::GeneticCode) = code @@ -43,7 +50,7 @@ function Base.show(io::IO, ::MIME"text/plain", code::GeneticCode) print(io, " ") for z in rna codon = unambiguous_codon(x, y, z) - aa = code[codon] + aa = getindex_encoding(code, codon) print(io, x, y, z, ": ", aa) if z != RNA_U print(io, " ") @@ -61,7 +68,7 @@ function Base.iterate(code::GeneticCode, x = UInt64(0)) if x > UInt64(0b111111) return nothing else - return (x, @inbounds code[x]), x + 1 + return (x, getindex_encoding(code, x)), x + 1 end end @@ -93,11 +100,11 @@ function Base.show(io::IO, trans::TransTables) end """ -Genetic code list of NCBI. +List of genetic codes from NCBI. The standard genetic code is `ncbi_trans_table[1]` and others can be shown by `show(ncbi_trans_table)`. -For more details, consult the next link: +For more details, consult the following link: http://www.ncbi.nlm.nih.gov/Taxonomy/taxonomyhome.html/index.cgi?chapter=cgencodes. """ const ncbi_trans_table = TransTables() @@ -326,7 +333,7 @@ Base3 = TCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAG """ translate(seq; kwargs...)::LongAA -Same as [`translate!`](@ref), but allocate a new `LongAA` to hold the result. +Same as [`translate!`](@ref), but allocate a new `LongAA` to hold the result. """ function translate(ntseq::SeqOrView; code::GeneticCode = standard_genetic_code, @@ -372,7 +379,7 @@ MTAJBX # Extended help `seq` must be either `LongSequence` or `LongSubSeq` -of alphabets either `DNAAlphabet` or `RNAAlphabet`. +with an alphabet of either `DNAAlphabet` or `RNAAlphabet`. Translation is done using genetic code `code`, which defaults to the standard genetic code. @@ -381,14 +388,14 @@ If the length of `seq` is not divisible by 3, throw an error. If `seq` contains any gaps (i.e. `DNA_Gap` or `RNA_Gap`), throw an error. -If `allow_ambiguous_codons` is set (default), codons with ambiguous nucleotides -such as `DNA_W` will be translated to the most narrow amino acid which covers all +If `allow_ambiguous_codons` is set (the default), codons with ambiguous nucleotides +such as `DNA_W` will be translated to the narrowest amino acid that covers all non-ambiguous codons encompassed by the ambiguous codon. This process might slow translation down. -If not `allow_ambiguous_codons`, ambiguous codons will error. +If `allow_ambiguous_codons` is not set, ambiguous codons will cause an error. If `alternative_start` is set, the leading methionine will be set to `AA_M`, no -matter the starting codon, or the code. +matter the starting codon or the code. """ function translate!(aaseq::LongAA, ntseq::SeqOrView{<:NucleicAcidAlphabet{2}}; @@ -404,7 +411,7 @@ function translate!(aaseq::LongAA, b = ntseq[3i-1] c = ntseq[3i] codon = unambiguous_codon(a, b, c) - aaseq[i] = code[codon] + aaseq[i] = getindex_encoding(code, codon) end alternative_start && !isempty(aaseq) && (@inbounds aaseq[1] = AA_M) aaseq @@ -426,7 +433,7 @@ function translate!(aaseq::LongAA, if isgap(a) | isgap(b) | isgap(c) error("Cannot translate nucleotide sequences with gaps.") elseif iscertain(a) & iscertain(b) & iscertain(c) - aaseq[i] = code[unambiguous_codon(a, b, c)] + aaseq[i] = getindex_encoding(code, unambiguous_codon(a, b, c)) else aaseq[i] = try_translate_ambiguous_codon(code, a, b, c, allow_ambiguous_codons) end @@ -445,9 +452,9 @@ function try_translate_ambiguous_codon( ((a, b, c), unambigs) = Iterators.peel( Iterators.product(map(UnambiguousRNAs, (x, y, z))...) ) - aa = @inbounds code[unambiguous_codon(a, b, c)] + aa = getindex_encoding(code, unambiguous_codon(a, b, c)) @inbounds for (a, b, c) in unambigs - aa_new = code[unambiguous_codon(a, b, c)] + aa_new = getindex_encoding(code, unambiguous_codon(a, b, c)) aa_new == aa && continue allow_ambiguous || error("codon ", a, b, c, " cannot be unambiguously translated") aa = if aa_new in (AA_N, AA_D) && aa in (AA_N, AA_D, AA_B) @@ -475,4 +482,3 @@ function Base.iterate(x::UnambiguousRNAs, state=reinterpret(UInt8, x.x)) rna = reinterpret(RNA, 0x01 << (trailing_zeros(state) & 7)) (rna, state & (state - 0x01)) end - diff --git a/test/translation.jl b/test/translation.jl index ce4f19ba5..74cda4a39 100644 --- a/test/translation.jl +++ b/test/translation.jl @@ -64,6 +64,8 @@ # Basics @test length(BioSequences.standard_genetic_code) == 64 + @test_throws KeyError BioSequences.standard_genetic_code[UInt64(64)] + @test_throws KeyError BioSequences.standard_genetic_code[typemax(UInt64)] buf = IOBuffer() show(buf, BioSequences.standard_genetic_code) @test !iszero(length(take!(buf))) # just test it doesn't error From 309fc88967367dd4e3f6649cf1dc79f91e6fad98 Mon Sep 17 00:00:00 2001 From: Jakob Nybo Andersen Date: Sat, 29 Aug 2026 15:21:05 +0200 Subject: [PATCH 3/9] Fix: Avoid infinite loop when matching empty strings When matching empty strings in sequence search, an edge case meant that the search could not advance, resulting in an infinite loop. Always advance the search even for zero-length results. --- src/BioSequences.jl | 2 +- src/search/re.jl | 5 +++-- test/longsequences/find.jl | 1 + test/search/regex.jl | 1 + 4 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/BioSequences.jl b/src/BioSequences.jl index e6441df8a..e8a61f9b6 100644 --- a/src/BioSequences.jl +++ b/src/BioSequences.jl @@ -248,7 +248,7 @@ search(query, itr; overlap = DEFAULT_OVERLAP) = Search(query, itr, overlap) function Base.iterate(itr::Search, state=firstindex(itr.itr)) val = findnext(itr.query, itr.itr, state) val === nothing && return nothing - state = itr.overlap ? first(val) + 1 : last(val) + 1 + state = itr.overlap ? first(val) + 1 : max(first(val), last(val)) + 1 return val, state end diff --git a/src/search/re.jl b/src/search/re.jl index 02f2980cb..137a3bcfb 100644 --- a/src/search/re.jl +++ b/src/search/re.jl @@ -767,8 +767,9 @@ function advance!(threads, captured, s, re, seq, overlap) if runmatch!(threads, captured, re, seq) if !overlap empty!(threads) - s = captured[2] - if s <= lastindex(seq) + s = max(captured[1] + 1, captured[2]) + if s <= lastindex(seq) || + (captured[1] == captured[2] && s == lastindex(seq) + 1) push!(threads, (1, s)) end end diff --git a/test/longsequences/find.jl b/test/longsequences/find.jl index fe8ca96bb..2cdcdcd60 100644 --- a/test/longsequences/find.jl +++ b/test/longsequences/find.jl @@ -229,6 +229,7 @@ end # Check overlap key argument. @test findall(ExactSearchQuery(dna"GATC", iscompatible), seq; overlap = false) == [1:4, 6:9, 10:13] @test findall(ExactSearchQuery(dna"GATC", iscompatible), seq; overlap = true) == [1:4, 2:5, 6:9, 10:13] + @test findall(ExactSearchQuery(dna""), dna"AC"; overlap = false) == [1:0, 2:1, 3:2] # Check mapping of indices. @test findall(DNA_A, seq, 7:11) == [7, 11] diff --git a/test/search/regex.jl b/test/search/regex.jl index 5e65f333f..a03151965 100644 --- a/test/search/regex.jl +++ b/test/search/regex.jl @@ -88,6 +88,7 @@ for (i, m) in enumerate(collect(eachmatch(biore"GC|CG"d, dna"ACGTTATGCATGGCG", false))) @test matched(m) == matches[i] end + @test matched.(collect(eachmatch(BioRegex{DNA}(""), dna"AC", false))) == [dna"", dna"", dna""] # matchall matchall(pat, seq, overlap=true) = From 473831d8bf434038e51a5adf918113c8baa48eff Mon Sep 17 00:00:00 2001 From: Jakob Nybo Andersen Date: Sat, 29 Aug 2026 16:14:15 +0200 Subject: [PATCH 4/9] Fix: Refactor join of BioSequence Fix issue where join! on unknown-sized iterators did not shrink the target sequence after joining. Also redesign it: Assume AbstractArray, AbstractSet and Tuple are stateless, and for all other iterables, collect the iterable first. Then, precompute the full length, copy in the content and verify the size matches. --- src/biosequence/biosequence.jl | 34 ++++++++++++++++---------------- test/biosequences/biosequence.jl | 5 +++++ test/longsequences/basics.jl | 11 ++++++++++- 3 files changed, 32 insertions(+), 18 deletions(-) diff --git a/src/biosequence/biosequence.jl b/src/biosequence/biosequence.jl index 85c0abcdb..82fad5d42 100644 --- a/src/biosequence/biosequence.jl +++ b/src/biosequence/biosequence.jl @@ -115,9 +115,11 @@ function Base.similar(seq::BioSequence, len::Integer=length(seq)) return typeof(seq)(undef, len) end -# Fast path for iterables we know are stateless -function join!(seq::BioSequence, it::Union{Vector, Tuple, Set}) - _join!(resize!(seq, reduce((a, b) -> a + joinlen(b), it, init=0)), it, Val(true)) +# Fast path for collection types that can be iterated twice. +function join!(seq::BioSequence, it::Union{AbstractArray, AbstractSet, Tuple}) + len = sum(joinlen, it, init=0) + resize!(seq, len) + _join!(seq, it, len) end """ @@ -134,25 +136,23 @@ TAGAAC See also [`join`](@ref). """ -join!(seq::BioSequence, it) = _join!(seq, it, Val(false)) +join!(seq::BioSequence, it) = join!(seq, collect(it)::AbstractArray) -# B is whether the size of the destination seq is already -# known to be the final size -function _join!(seq::BioSequence, it, ::Val{B}) where B +function _join!(seq::BioSequence, it, target_len::Int) len = 0 - oldlen = length(seq) for i in it pluslen = joinlen(i) - if !B && oldlen < (len + pluslen) - resize!(seq, len + pluslen) + if len + pluslen > target_len + throw(ArgumentError("collection changed while joining")) end if i isa BioSymbol seq[len + 1] = i else - copyto!(seq, len + 1, i, 1, length(i)) + copyto!(seq, len + 1, i, 1, pluslen) end len += pluslen end + len == target_len || throw(ArgumentError("collection changed while joining")) seq end @@ -170,16 +170,16 @@ TAGAAC See also [`join!`](@ref). """ -function Base.join(::Type{T}, it::Union{Vector, Tuple, Set}) where {T <: BioSequence} - _join!(T(undef, reduce((a, b) -> a + joinlen(b), it, init=0)), it, Val(true)) +function Base.join(::Type{T}, it::Union{AbstractArray, AbstractSet, Tuple}) where {T <: BioSequence} + len = sum(joinlen, it, init=0) + _join!(T(undef, len), it, len) end # length is intentionally not implemented for BioSymbol -joinlen(x::Union{BioSequence, BioSymbol}) = x isa BioSymbol ? 1 : length(x) +joinlen(::BioSymbol) = 1 +joinlen(seq::BioSequence) = length(seq) -function Base.join(::Type{T}, it) where {T <: BioSequence} - _join!(empty(T), it, Val(false)) -end +Base.join(::Type{T}, it) where {T <: BioSequence} = join(T, collect(it)::AbstractArray) Base.repeat(chunk::BioSequence, n::Integer) = join(typeof(chunk), (chunk for i in 1:n)) Base.:^(x::BioSequence, n::Integer) = repeat(x, n) diff --git a/test/biosequences/biosequence.jl b/test/biosequences/biosequence.jl index 834c94d99..44a2ae677 100644 --- a/test/biosequences/biosequence.jl +++ b/test/biosequences/biosequence.jl @@ -19,6 +19,9 @@ BioSequences.encoded_setindex!(x::SimpleSeq, e::UInt, i::Integer) = x.x[i] = e SimpleSeq(::UndefInitializer, x::Integer) = SimpleSeq(Unsafe(), zeros(UInt, x)) Base.resize!(x::SimpleSeq, len::Int) = (resize!(x.x, len); x) +struct BadCollectIterator end +Base.collect(::BadCollectIterator) = BadCollectIterator() + # Not part of the API, just used for testing purposes random_simple(len::Integer) = SimpleSeq(rand([RNA_A, RNA_C, RNA_G, RNA_U], len)) @@ -63,6 +66,8 @@ random_simple(len::Integer) = SimpleSeq(rand([RNA_A, RNA_C, RNA_G, RNA_U], len)) @test join(SimpleSeq, [seq, seq2]) == join!(SimpleSeq([]), [seq, seq2]) @test join(SimpleSeq, gen) == join!(SimpleSeq([]), gen) @test join(SimpleSeq, [RNA_U, RNA_G, seq, RNA_U]) == SimpleSeq([RNA(i) for i in "UGCGUU"]) + @test_throws TypeError join!(SimpleSeq([]), BadCollectIterator()) + @test_throws TypeError join(SimpleSeq, BadCollectIterator()) @test copy!(SimpleSeq([]), seq) == seq seq3 = copy(seq2) diff --git a/test/longsequences/basics.jl b/test/longsequences/basics.jl index ca004f987..28677202f 100644 --- a/test/longsequences/basics.jl +++ b/test/longsequences/basics.jl @@ -277,6 +277,15 @@ end test_join(LongAA, [AA_G, AA_P, AA_L, aa"MNVWEED", AA_K], aa"GPLMNVWEEDK") test_join(LongRNA{4}, [RNA_M, RNA_U, RNA_S, rna"AGCGSK"], rna"MUSAGCGSK") test_join(LongDNA{2}, [dna"ATGCTTA", DNA_G, DNA_G, DNA_A, DNA_A, DNA_A], dna"ATGCTTAGGAAA") + + # Arbitrary iterators are collected before joining and resize an existing + # destination to their final length. + for (n, result) in [(0, dna""), (1, dna"C"), (4, dna"CCCC"), (6, dna"CCCCCC")] + @test join(LongDNA{4}, (dna"C" for _ in 1:n)) == result + seq = LongDNA{4}("AAAA") + @test join!(seq, (dna"C" for _ in 1:n)) == result + @test length(seq) == n + end end @testset "Length" begin @@ -407,4 +416,4 @@ end buf = IOBuffer() print(buf, LongSequence{ReducedAAAlphabet}("AGTDNWLE")) @test String(take!(buf)) == "AGTDNWLE" -end \ No newline at end of file +end From 1020db6cd25c54af3217067ca772038334dd2e71 Mon Sep 17 00:00:00 2001 From: Jakob Nybo Andersen Date: Sat, 29 Aug 2026 17:56:32 +0200 Subject: [PATCH 5/9] Fix: Fix findnext(!=, ...) for unencodable symbols E.g. `findnext(!=(DNA_A), rna"A")` should just return 1 without searching the RNA because DNA is unencodable in the RNA alphabet. --- src/longsequences/operators.jl | 10 ++++++---- test/longsequences/find.jl | 7 +++++++ 2 files changed, 13 insertions(+), 4 deletions(-) diff --git a/src/longsequences/operators.jl b/src/longsequences/operators.jl index 0d30713b3..47ad6c6b6 100644 --- a/src/longsequences/operators.jl +++ b/src/longsequences/operators.jl @@ -281,8 +281,9 @@ function _findnext( seq::SeqOrView{<:KNOWN_ALPHABETS}, i::Int, ) - enc = tryencode(Alphabet(seq), cmp.x) - enc === nothing && return nothing + x = cmp isa Base.Fix2 ? cmp.x : cmp.inner.x + enc = tryencode(Alphabet(seq), x) + enc === nothing && return i u_enc = enc * encoding_expansion(BitsPerSymbol(seq)) f = x -> x ⊻ u_enc _findnext_nonzero(f, seq, i) @@ -363,8 +364,9 @@ function _findprev( seq::SeqOrView{<:KNOWN_ALPHABETS}, i::Int, ) - enc = tryencode(Alphabet(seq), cmp.x) - enc === nothing && return nothing + x = cmp isa Base.Fix2 ? cmp.x : cmp.inner.x + enc = tryencode(Alphabet(seq), x) + enc === nothing && return i u_enc = enc * encoding_expansion(BitsPerSymbol(seq)) f = x -> x ⊻ u_enc _findprev_nonzero(f, seq, i) diff --git a/test/longsequences/find.jl b/test/longsequences/find.jl index 2cdcdcd60..c92e548a5 100644 --- a/test/longsequences/find.jl +++ b/test/longsequences/find.jl @@ -13,6 +13,13 @@ @test findfirst(==(AA_W), seq) === nothing @test findfirst(==('C'), seq) === nothing + # An unencodable symbol is unequal to every element of a two-bit sequence. + shortseq = LongDNA{2}("AC") + @test findnext(!=(DNA_N), shortseq, 1) === 1 + @test findprev(!=(DNA_N), shortseq, 2) === 2 + @test findnext(!isequal(DNA_N), shortseq, 1) === 1 + @test findprev(!isequal(DNA_N), shortseq, 2) === 2 + @test findlast(isequal(DNA_C), seq) == 14 @test findlast(isequal(DNA_G), seq) == 15 @test findlast(isequal(DNA_Y), seq) == 11 From 319d1690f30270a68766b9bc1ff6f77e3fbb42ba Mon Sep 17 00:00:00 2001 From: Jakob Nybo Andersen Date: Sat, 29 Aug 2026 18:05:20 +0200 Subject: [PATCH 6/9] Fix: Reject filter! for SubSeq Previously, this would mutate the sequence and then throw when attempting to resize!. Now, add a method that throws up front. This does not solve the problem generally for non-resizable BioSequences, but I don't think there is a good way to address that. At least we can add this safety check even while not solving it generically. --- src/longsequences/transformations.jl | 6 ++++++ test/longsequences/transformations.jl | 6 ++++++ 2 files changed, 12 insertions(+) diff --git a/src/longsequences/transformations.jl b/src/longsequences/transformations.jl index 3dcb891af..c09f0dfdd 100644 --- a/src/longsequences/transformations.jl +++ b/src/longsequences/transformations.jl @@ -24,6 +24,12 @@ function Base.resize!(seq::LongSequence{A}, size::Integer, force::Bool=false) wh end end +# A view may mutate its elements, but cannot change its length. Reject length- +# changing operations before the generic `filter!` implementation compacts data. +function Base.filter!(f, seq::LongSubSeq) + throw(ArgumentError("filter! cannot change the length of a LongSubSeq; use filter instead")) +end + """ reverse!(seq::LongSequence) diff --git a/test/longsequences/transformations.jl b/test/longsequences/transformations.jl index 33b8a6ee1..0b19bf31f 100644 --- a/test/longsequences/transformations.jl +++ b/test/longsequences/transformations.jl @@ -138,6 +138,12 @@ @test filter!(x -> x != DNA_N, seq) == seq @test seq == dna"ACGTACGT" + # A view cannot be resized, so filtering it must fail before compaction + # can overwrite its parent. + seq = dna"ACGT" + @test_throws ArgumentError filter!(x -> x == DNA_G, view(seq, 1:4)) + @test seq == dna"ACGT" + for len in [1, 2, 3, 5, 8, 9, 15, 19, 31, 32, 33, 50], _ in 1:5 str = random_dna(len) seq = LongDNA{4}(str) From 19fa4e1f54363f7103c1ad122531dd8fab7759bb Mon Sep 17 00:00:00 2001 From: Jakob Nybo Andersen Date: Sat, 29 Aug 2026 18:13:09 +0200 Subject: [PATCH 7/9] Fix: Yield Pair from GeneticCode, not Tuple GeneticCode subtypes AbstractDict, and so must yield Pairs. --- src/geneticcode.jl | 2 +- test/translation.jl | 6 ++++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/src/geneticcode.jl b/src/geneticcode.jl index a6d3d6fd1..5a6c60033 100644 --- a/src/geneticcode.jl +++ b/src/geneticcode.jl @@ -68,7 +68,7 @@ function Base.iterate(code::GeneticCode, x = UInt64(0)) if x > UInt64(0b111111) return nothing else - return (x, getindex_encoding(code, x)), x + 1 + return x => getindex_encoding(code, x), x + 1 end end diff --git a/test/translation.jl b/test/translation.jl index 74cda4a39..f0b5d63ec 100644 --- a/test/translation.jl +++ b/test/translation.jl @@ -66,6 +66,12 @@ @test length(BioSequences.standard_genetic_code) == 64 @test_throws KeyError BioSequences.standard_genetic_code[UInt64(64)] @test_throws KeyError BioSequences.standard_genetic_code[typemax(UInt64)] + code = BioSequences.standard_genetic_code + entries = collect(code) + @test entries == collect(pairs(code)) + @test eltype(entries) == Pair{UInt64,AminoAcid} + @test collect(keys(code)) == UInt64.(0:63) + @test collect(values(code)) == [code[key] for key in UInt64.(0:63)] buf = IOBuffer() show(buf, BioSequences.standard_genetic_code) @test !iszero(length(take!(buf))) # just test it doesn't error From 17c64612b575ac963cf476bb774b3d77783262ce Mon Sep 17 00:00:00 2001 From: Jakob Nybo Andersen Date: Sat, 29 Aug 2026 18:22:22 +0200 Subject: [PATCH 8/9] Fix: Actually implement findall for ApproxSearchQuery This was documented but not actually implemented. Implementation is trivial. --- src/search/ApproxSearchQuery.jl | 18 ++++++++++++++++++ test/search/ApproximateSearchQuery.jl | 9 +++++++++ 2 files changed, 27 insertions(+) diff --git a/src/search/ApproxSearchQuery.jl b/src/search/ApproxSearchQuery.jl index 3e4e5b8dd..fe0e74c47 100644 --- a/src/search/ApproxSearchQuery.jl +++ b/src/search/ApproxSearchQuery.jl @@ -134,6 +134,24 @@ Base.findfirst(query::ApproximateSearchQuery, k::Integer, seq::BioSequence) = fi Base.findlast(query::ApproximateSearchQuery, k::Integer, seq::BioSequence) = findprev(query, k, seq, lastindex(seq)) Base.occursin(query::ApproximateSearchQuery, k::Integer, seq::BioSequence) = !isnothing(_approxsearch(query, k, seq, 1, lastindex(seq), true)) +""" + findall(query, k, seq; overlap=true) + +Return all occurrences of `query` in `seq` allowing up to `k` errors. +""" +function Base.findall(query::ApproximateSearchQuery, k::Integer, seq::BioSequence; overlap::Bool = true) + matches = UnitRange{Int}[] + start = firstindex(seq) + stop = lastindex(seq) + 1 + while start <= stop + match = findnext(query, k, seq, start) + match === nothing && return matches + push!(matches, match) + start = overlap ? first(match) + 1 : max(first(match), last(match)) + 1 + end + return matches +end + function _approxsearch(query::ApproximateSearchQuery, k::Integer, seq::BioSequence, start::Integer, stop::Integer, forward::Bool) checkeltype(query.seq, seq) if k ≥ length(query.seq) diff --git a/test/search/ApproximateSearchQuery.jl b/test/search/ApproximateSearchQuery.jl index 7071832af..24089337c 100644 --- a/test/search/ApproximateSearchQuery.jl +++ b/test/search/ApproximateSearchQuery.jl @@ -1,6 +1,15 @@ @testset "Approximate" begin seq = dna"ACGTACG" + @testset "findall" begin + query = ApproximateSearchQuery(dna"ACA") + @test findall(query, 0, dna"ACACA") == [1:3, 3:5] + @test findall(query, 0, dna"ACACA"; overlap=false) == [1:3] + @test findall(ApproximateSearchQuery(dna"ACG"), 1, seq) == [1:2, 2:3, 5:6, 6:7] + @test findall(ApproximateSearchQuery(dna""), 0, dna"AC") == [1:0, 2:1, 3:2] + @test findall(query, 0, dna"GGGG") isa Vector{UnitRange{Int}} + end + @testset "iscompatible" begin @testset "forward" begin @test findnext(ApproximateSearchQuery(dna"", iscompatible), 0, seq, 1) === 1:0 From 540ce61649e4b705f12fc30f986e4122d24b2c38 Mon Sep 17 00:00:00 2001 From: Jakob Nybo Andersen Date: Sun, 30 Aug 2026 16:05:44 +0200 Subject: [PATCH 9/9] Bump version --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 11b929b57..5b3261a07 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "BioSequences" uuid = "7e6ae17a-c86d-528c-b3b9-7f778a29fe59" -authors = ["Sabrina Jaye Ward ", "Jakob Nissen "] -version = "3.5.1" +authors = ["Sabrina Jaye Ward ", "Jakob Nybo Andersen "] +version = "3.5.2" [deps] BioSymbols = "3c28c6f8-a34d-59c4-9654-267d177fcfa9"