Before diving into searching and matching in Biostrings, let’s first briefly review the central data structures and concepts in the package. If you’re already familiar with XString objects and the various containers to hold them, feel free to skip forward to the next section.
The focus of Biostrings is providing containers to efficiently work
with biological data. Operations in Biostrings revolve around XString1 objects.
XString isn’t actually a type; it’s the collective name for a set of
objects of the form *String:
Aside from BString, all the XString types also support standard IUPAC ambiguity codes.
The simplest way to create XString objects is directly from a character object:
## 8-letter DNAString object
## seq: ATGCATGC## 8-letter RNAString object
## seq: AUGCAUGC## 6-letter AAString object
## seq: ARNDC*## 14-letter BString object
## seq: I am a BStringSingle XString objects are nice, but many biological sequence analyses need to work with multiple sequences. Our primary container to achieve this functionality is the XStringSet2 object. This container holds a set of XString objects that are all the same type (all DNA, all RNA, etc.). There is one XStringSet type for each XString type (i.e., DNA, RNA, AA, B).
It’s also possible to create XStringSet objects from character vectors:
## DNAStringSet object of length 4:
## width seq
## [1] 4 ATGC
## [2] 4 CGTA
## [3] 5 AAAAA
## [4] 2 TTHowever, they are most often used to read in sequence sets from files, such as FASTA files:
## DNAStringSet object of length 7:
## width seq names
## [1] 5573 ACTTGTAAATATATCTTTTATTT...CTTATCGACCTTATTGTTGATAT YAL001C TFC3 SGDI...
## [2] 5825 TTCCAAGGCCGATGAATTCGACT...AGTAAATTTTTTTCTATTCTCTT YAL002W VPS8 SGDI...
## [3] 2987 CTTCATGTCAGCCTGCACTTCTG...TGGTACTCATGTAGCTGCCTCAT YAL003W EFB1 SGDI...
## [4] 3929 CACTCATATCGGGGGTCTTACTT...TGTCCCGAAACACGAAAAAGTAC YAL005C SSA1 SGDI...
## [5] 2648 AGAGAAAGAGTTTCACTTCTTGA...ATATAATTTATGTGTGAACATAG YAL007C ERP2 SGDI...
## [6] 2597 GTGTCCGGGCCTCGCAGGCGTTC...AAGTTTTGGCAGAATGTACTTTT YAL008W FUN14 SGD...
## [7] 2780 CAAGATAATGTCAAAGTTAGTGG...GCTAAGGAAGAAAAAAAAATCAC YAL009W SPO7 SGDI...Once we have an XStringSet, we can use single brackets
([]) to get a subset of the XStringSet, and double brackets
([[]]) to get the component XStrings:
## DNAStringSet object of length 1:
## width seq
## [1] 4 ATGC## 4-letter DNAString object
## seq: ATGCWhen we search for patterns in a sequence, we’re looking for specific subsequences within a single sequence that match whatever we’re looking for. To implement that functionality, we have the XStringViews3 classes, which store “views” (subsequences) of a single XString object. These will be very important for searching and matching patterns later on.
We can directly construct an XStringViews from an XString object with
the Views method by by either specifying start/end indices,
or using the width of the subsequences:
seq <- DNAString("AAATTTGGGCCC")
## Splitting into codons
Views(seq,
start = seq(1, length(seq), by = 3),
end = seq(3, length(seq), by = 3))## Views on a 12-letter DNAString subject
## subject: AAATTTGGGCCC
## views:
## start end width
## [1] 1 3 3 [AAA]
## [2] 4 6 3 [TTT]
## [3] 7 9 3 [GGG]
## [4] 10 12 3 [CCC]## Views on a 12-letter DNAString subject
## subject: AAATTTGGGCCC
## views:
## start end width
## [1] 1 3 3 [AAA]
## [2] 2 4 3 [AAT]
## [3] 3 5 3 [ATT]
## [4] 4 6 3 [TTT]
## [5] 5 7 3 [TTG]
## [6] 6 8 3 [TGG]
## [7] 7 9 3 [GGG]
## [8] 8 10 3 [GGC]
## [9] 9 11 3 [GCC]
## [10] 10 12 3 [CCC]## Views on a 12-letter DNAString subject
## subject: AAATTTGGGCCC
## views:
## start end width
## [1] -1 1 3 [ A]
## [2] 0 2 3 [ AA]
## [3] 1 3 3 [AAA]
## [4] 2 4 3 [AAT]
## [5] 3 5 3 [ATT]
## ... ... ... ... ...
## [8] 6 8 3 [TGG]
## [9] 7 9 3 [GGG]
## [10] 8 10 3 [GGC]
## [11] 9 11 3 [GCC]
## [12] 10 12 3 [CCC]This class will be commonly returned from the pattern matching
functions we’ll discuss later, as well as some other internal functions
(e.g., codons, which returns all the codons in a
DNAString)
It’s worth noting that we can initialize an XStringSet from an XStringViews, in case we want to save our subsequences as a new XStringSet object.
seq <- DNAString("AAATTTGGGCCC")
## `codons` returns an XStringViews
seq_codons <- codons(seq)
seq_codons## Views on a 12-letter DNAString subject
## subject: AAATTTGGGCCC
## views:
## start end width
## [1] 1 3 3 [AAA]
## [2] 4 6 3 [TTT]
## [3] 7 9 3 [GGG]
## [4] 10 12 3 [CCC]## DNAStringSet object of length 4:
## width seq
## [1] 3 AAA
## [2] 3 TTT
## [3] 3 GGG
## [4] 3 CCCKnowledge of XStrings, XStringSets, and XStringViews are all that
you’ll need to be able to understand the rest of this vignette. However,
there exist a number of additional containers and utility functions to
work with these objects that I haven’t mentioned here. For more
information on these utility functions, look at ?XString.
Other containers that may be of interest include XStringSetList,
MaskedXString, and QualityScaledXStringSet.
Now we can start to look at how to perform some basic matching
operations in Biostrings. Matching in Biostrings is done in terms of
“patterns” and “subjects”. A pattern is what we’re searching
for, and a subject is where we’re searching. For example, if we
wanted to find all ATG codons in an E. coli
genome, ATG would be the pattern, and the E. coli
genome is the subject. A match is a subsequence of the pattern
that is equal to the pattern (or similar enough, if using things like
ambiguity codes).
Before we can get to matching, we first need a way to quantify the
edit distance between two strings. The edit distance is the
minimum number of operations it takes to transform one string into
another. For instance, the edit distance between AAA and
AAB is 1, since the only letter that differs between the
two strings is the final character.
## the first string is the pattern,
## the second is the subject
neditAt(pattern = "AAAAA", subject = "AAAAB")## [1] 1## [1] 4## Going from subject -> pattern
## AACGCA -> AACGTA -> AACCTA -> AAGCTA -> GAGCTA: 4 edits
neditAt("GAGCTA", "AACGCA")## [1] 4However, this only works for strings of the same length. If we wanted
to get the edit distance between two strings of unequal length, we’d
also have to consider insertions and deletions (e.g., AAA
and AC would have a distance of two via AAA
-> AA -> AC). Considering insertions,
deletions, and substitutions produces the Levenstein distance, which is
also supported by neditAt if we set
with.indels=TRUE:
## [1] 3Note that this tries to match the pattern against the subject, so it
only looks at the first n characters of the subject, where
n is the length of the pattern. This can lead to some
unintuitive results:
## Both of these will return 0, because nchar("ABC") = 3,
## so we only look at the first 3 characters of the pattern ("ABC")
neditAt("ABC", "ABCD")## [1] 0## [1] 0## However, if we change the subject slightly
## (note just showing the first 3 characters of the subject)
neditAt("ABC", "DBCA") # DBCA -> ABCA: 1## [1] 1## [1] 3## [1] 2The previous result for neditAt probably seems
unintuitive. Why is this function set up the way it is?
It turns out that we can use this neditAt function to
find matches. We can define a “match” as a pattern that has an edit
distance of 0 with the subject.
## [1] TRUE## [1] TRUEThese functions have an “at” suffix because we can specify where we
should start searching in the subject. By default, we start searching at
the first character (1), but we can start anywhere:
## [1] 3## [1] 0## [1] 0Both functions are vectorized to calculate for multiple starting positions at the same time:
## [1] 3 0 3 3 0 3 3## [1] 2 5## [1] 2 5What if we want to be a little more lenient than exact matching?
Previously, we defined a match as any location where
neditAt returns 0. We can relax these conditions with the
min.mismatch and max.mismatch arguments. With
these, the definition of a “match” becomes anywhere where
min.mismatch <= neditAt <= max.mismatch.
pattern <- "AAB"
subject <- "AABBAAA"
s <- seq_len(nchar(subject))
neditAt(pattern, subject, at = s)## [1] 0 1 3 2 1 1 2## [1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE## If max.mismatch = 1, we get more matches
isMatchingAt(pattern, subject, at = s, max.mismatch = 1)## [1] TRUE TRUE FALSE FALSE TRUE TRUE FALSE## We can also restrict the lower bound
## This returns all patterns with neditAt == 2
isMatchingAt(pattern, subject, at = s, min.mismatch = 2, max.mismatch = 2)## [1] FALSE FALSE FALSE TRUE FALSE FALSE TRUENote that since we’re using neditAt to define a match,
the we can also use the Levenshtein distance as our matcher:
## [1] 3## [1] FALSE## [1] 2## [1] TRUEWith isMatchingAt, we have a basic way to get matches in
pattern/subject syntax. All future sections will build off this
foundational function.
Note that, while we used character vectors for the examples in this section, all the functions shown also work on XString and XStringSet objects.
In the previous section, we talked about how to check if a particular pattern has a match at a specific location in a sequence. While this is theoretically sufficient to do any arbitrary pattern matching, the API is pretty clunky. Ideally, we could just call a single function to automatically find a pattern throughout the subject, without having to worry about where to look.
matchPatternThe matchPattern function does exactly what we’re
looking for. This function has the same input arguments as
isMatchingAt, but searches in the entire subject for the
pattern rather than at specific points.
subject <- "BABCABC"
edit_dists <- neditAt("ABC", subject, at = seq_len(nchar(subject)))
which(edit_dists == 0)## [1] 2 5## Views on a 7-letter BString subject
## subject: BABCABC
## views:
## start end width
## [1] 2 4 3 [ABC]
## [2] 5 7 3 [ABC]The output is almost the same as running isMatchingAt
over the entire subject, except that the return type is an XStringViews
object.
XStringViews objects are perfect for pattern matching, because they can hold all the subsequences that match in the subject along with their start and end positions. Let’s see what that looks like using DNAStrings:
## Views on a 9-letter DNAString subject
## subject: AGTAGCAGA
## views:
## start end width
## [1] 1 2 2 [AG]
## [2] 4 5 2 [AG]
## [3] 7 8 2 [AG]Here we can see that the pattern AG is found starting at
indices 1, 4, and 7. In this case, all of our matches are identical to
the pattern, but this may not always be the case.
For instance, if we set fixed = FALSE, then we can use
IUPAC ambiguity codes:
pattern <- DNAString("AGY") # Y = T or G
subject <- DNAString("AGTAGCAGAAGY")
## By default, ambiguities only match ambiguities
## (i.e., Y = Y)
matchPattern(pattern, subject)## Views on a 12-letter DNAString subject
## subject: AGTAGCAGAAGY
## views:
## start end width
## [1] 10 12 3 [AGY]## But with fixed = FALSE, it can match other letters
## (i.e., Y = Y, T, or G)
matchPattern(pattern, subject, fixed = FALSE)## Views on a 12-letter DNAString subject
## subject: AGTAGCAGAAGY
## views:
## start end width
## [1] 1 3 3 [AGT]
## [2] 4 6 3 [AGC]
## [3] 10 12 3 [AGY]Similarly, if we use the max.mismatch arguments, we can
get variable length matches:
pattern <- DNAString("AGY") # Y = T or G
subject <- DNAString("AGTAGCAGAAGY")
## But with fixed = FALSE, it can match other letters
## (i.e., Y = Y, T, or G)
matchPattern(pattern, subject, fixed = FALSE, max.mismatch = 2)## Views on a 12-letter DNAString subject
## subject: AGTAGCAGAAGY
## views:
## start end width
## [1] 1 3 3 [AGT]
## [2] 4 6 3 [AGC]
## [3] 7 9 3 [AGA]
## [4] 9 11 3 [AAG]
## [5] 10 12 3 [AGY]Sometimes you may get too many matches to deal with. However, we can also match a pattern against an XStringViews subject. In other words, we can find matches in our matches to filter them out:
pattern1 <- DNAString("ANN") # Any 3-mer beginning with A
pattern2 <- DNAString("NNT") # Any 3-mer ending with T
subject <- DNAString("ATGTTTAATTTATAATATTTATAAAT")
first_matches <- matchPattern(pattern1, subject, fixed = FALSE)
first_matches## Views on a 26-letter DNAString subject
## subject: ATGTTTAATTTATAATATTTATAAAT
## views:
## start end width
## [1] 1 3 3 [ATG]
## [2] 7 9 3 [AAT]
## [3] 8 10 3 [ATT]
## [4] 12 14 3 [ATA]
## [5] 14 16 3 [AAT]
## [6] 15 17 3 [ATA]
## [7] 17 19 3 [ATT]
## [8] 21 23 3 [ATA]
## [9] 23 25 3 [AAA]
## [10] 24 26 3 [AAT]## Views on a 26-letter DNAString subject
## subject: ATGTTTAATTTATAATATTTATAAAT
## views:
## start end width
## [1] 7 9 3 [AAT]
## [2] 8 10 3 [ATT]
## [3] 14 16 3 [AAT]
## [4] 17 19 3 [ATT]
## [5] 24 26 3 [AAT]countPatternSometimes we don’t care specifically about what the matches are, we
just want to know how many matches exist. For this, we have
countPattern. countPattern returns the number
of matches that would be found by matchPattern.
pattern <- DNAString("ANN")
subject <- DNAString("ATGTTTAATTTATAATATTTATAAAT")
length(matchPattern(pattern, subject, fixed = FALSE))## [1] 10## [1] 10What if we want to search for a single pattern in multiple subjects? A common example of this is searching for a motif across multiple chromosomes.
We can accomplish this workflow using the vcountPattern
and vmatchPattern. The “v” stands for “vectorized”, since
it operates across all the strings in the subject like a vectorized
function.
While we can only use a single pattern (for now), the subject of these functions can have multiple sequences, like a character vector or XStringSet:
pattern <- DNAString("ATNG")
subjects <- DNAStringSet(c("ATTGATTAATGGATCG", "AGTGAGACAGATTG"))
vcountPattern(pattern, subjects, fixed = FALSE)## [1] 3 1## [1] 7 9## MIndex object of length 2
## [[1]]
## IRanges object with 5 ranges and 0 metadata columns:
## start end width
## <integer> <integer> <integer>
## [1] 1 4 4
## [2] 5 8 4
## [3] 8 11 4
## [4] 9 12 4
## [5] 13 16 4
##
## [[2]]
## IRanges object with 3 ranges and 0 metadata columns:
## start end width
## <integer> <integer> <integer>
## [1] 1 4 4
## [2] 7 10 4
## [3] 11 14 4vmatchPattern returns an object we haven’t seen before:
an MIndex. This container stores matches for a set of a
patterns against one or more subject sequences. Each entry in the
MIndex is an IRanges object, which is
basically an XStringViews without the substring of the
match attached (i.e., just the start/end position of each match and its
width).
pattern <- DNAString("ATNG")
subjects <- DNAStringSet(c("ATTGATTAATGGATCG", "AGTGAGACAGATTG"))
matches <- vmatchPattern(pattern, subjects, fixed = FALSE)
## 3 matches at indices 1, 9, and 13
matches[[1]]## IRanges object with 3 ranges and 0 metadata columns:
## start end width
## <integer> <integer> <integer>
## [1] 1 4 4
## [2] 9 12 4
## [3] 13 16 4## Views on a 16-letter DNAString subject
## subject: ATTGATTAATGGATCG
## views:
## start end width
## [1] 1 4 4 [ATTG]
## [2] 9 12 4 [ATGG]
## [3] 13 16 4 [ATCG]## IRanges object with 1 range and 0 metadata columns:
## start end width
## <integer> <integer> <integer>
## [1] 11 14 4## Views on a 14-letter DNAString subject
## subject: AGTGAGACAGATTG
## views:
## start end width
## [1] 11 14 4 [ATTG]matchPattern and countPattern are the first
“real” matching functions we’ve seen. While it’s possible to do complex
matching using isMatchingAt and neditAt, the
*Pattern functions provide true search functionality that
is useful in real scenarios.
One point of discussion we’ll skip in this vignette is the
algorithm parameter of matchPattern. The best
practice is to just leave it set to its default value of
"auto", which will automatically pick the best pattern for
your search. If a particular algorithm can be used for your search, it
will perform identically to all other usable algorithms, so you don’t
need to worry too much about this parameter.
We unfortunately won’t be able to do “multiple pattern, multiple
sequence” searching quite yet. For that, we’ll need PDict
objects, which we’ll talk about later.
PDict is short for “Preprocessed Dictionary”. A PDict object is essentially a container that holds a collection of DNA patterns to search with. The container preprocesses these patterns to allow fast matching against one or more subject sequences.
In the standard configuration, PDicts have the following limitations:
ATGC, no
ambiguity codes allowed)We’ll relax these limitations later.
## TB_PDict object of length 4 and width 4 (preprocessing algo="ACtree2")The PDict object shows how many patterns it holds (“length 4”) and how long the sequences are (“width 4”). We can also recover the sequences inside of it:
## 4-letter DNAString object
## seq: ATGCWe can also provide named sequences, and the PDict will keep track of them:
## [1] "seq1" "seq2"It also mentions the preprocessing algorithm, which must be one of the following:
ACtree2: The default processing algorithm. Supports
matching against subjects that have ambiguity codes (i.e.,
A in the PDict could match N in a
sequence).Twobit: Much faster than ACtree2, but
can’t match against ambiguity codes in subjects.Regardless of algorithm, we can’t include ambiguity codes in the patterns stored in the PDict. However, we can get around this and unlock the full power of PDicts by using a “Trusted Band”.
Once we have a PDict, we can start matching. These functions are very
similar to the *Pattern functions:
matchPDict is equivalent to
matchPatterncountPDict is equivalent to
countPatternvmatchPDict is similar to
vmatchPatternvcountPDict is equivalent to
vcountPatternThese enable us to do “many pattern, single subject” matching:
patterns <- DNAStringSet(c(seq1 = "ATGC", seq2 = "AATT", seq3 = "CCTA"))
pdict <- PDict(patterns)
subject <- DNAString("ATGCAAAATTTAATTGC")
## How many matches does each pattern have?
countPDict(pdict, subject)## [1] 1 2 0## MIndex object of length 3
## $seq1
## IRanges object with 1 range and 0 metadata columns:
## start end width
## <integer> <integer> <integer>
## [1] 1 4 4
##
## $seq2
## IRanges object with 2 ranges and 0 metadata columns:
## start end width
## <integer> <integer> <integer>
## [1] 7 10 4
## [2] 12 15 4
##
## $seq3
## IRanges object with 0 ranges and 0 metadata columns:
## start end width
## <integer> <integer> <integer>We also have whichPDict, which returns a vector
indicating which patterns have at least one match:
## [1] 1 2We can also do “many pattern, many subject” matching with the
v*PDict functions. However, the output types of these are
different than previously. Let’s break them down one by one.
First, vwhichPDict:
subjects <- DNAStringSet(c(subj1 = "ATGCAAAATTTAATTGC", subj2 = "CCAATCCTACTATGC"))
vwhichPDict(pdict, subjects)## [[1]]
## [1] 1 2
##
## [[2]]
## [1] 1 3This returns a list with the output of calling
whichPDict on each sequence in the subject. For this input,
we find that patterns 1 and 2 have matches in the first sequence, and
patterns 1 and 3 have matches in the second sequence.
If we wanted to find out how many matches each pattern has in each
sequence, we’d use vcountPDict:
## [,1] [,2]
## [1,] 1 1
## [2,] 2 0
## [3,] 0 1This function returns a matrix of counts. The rows correspond to
patterns, and the columns correspond to subjects. Thus, since
vcountPDict(pdict, subjects)[2,1] = 2, we know that pattern
2 has two matches in subject 1.
If we wanted to find out where those matches were, we’d use
vmatchPDict:
## MIndexList of length 2
## [["subj1"]] 3 patterns, 3 matches
## [["subj2"]] 3 patterns, 2 matchesThe result of this is an MIndexList class. This is a similar concept
to what we saw with the MIndex class when we discussed
matchPattern. The MIndex class is esssentially a list
storing offsets of matches for a set of patterns against a subject, and
the MIndexList class is a list storing an MIndex object for each
subject.
all_matches <- vmatchPDict(pdict, subjects)
## first entry is the matches for all patterns against the first subject:
all_matches[[1]]## MIndex object of length 3
## $seq1
## IRanges object with 1 range and 0 metadata columns:
## start end width
## <integer> <integer> <integer>
## [1] 1 4 4
##
## $seq2
## IRanges object with 2 ranges and 0 metadata columns:
## start end width
## <integer> <integer> <integer>
## [1] 7 10 4
## [2] 12 15 4
##
## $seq3
## IRanges object with 0 ranges and 0 metadata columns:
## start end width
## <integer> <integer> <integer>Note that the last entry of all_matches[[1]] is empty,
because the third pattern had no matches against the first subject.
Note: vmatchPDict is the only
function that supports “many pattern, many subject” matching.
vmatchPattern can only match a single pattern against one
or more sequences.
PDicts in their default configuration are very limited. However, this changes when we define a “Trusted Band”. PDicts with a trusted band support ambiguity codes, variable length sequences, and fuzzy/inexact matching.
The trusted band is a fixed size subsequence of the patterns that obeys the prior limitations (no ambiguities, fixed size). When we do matching, the PDict can first search for matches in the trusted band to limit the search space, and then expand to the regions of each pattern outside the trusted band.
The trusted bands do not need to all be in the same position, but they all need to be the same length. For example, we could imagine the following sequences with trusted bands:
Sequences:
----------
TTGCTAMKT
AGCTAGTTAAG
Some Possible Trusted Bands:
----------------------------
T | TGCT | AMKT
A | GCTA | GTTAAG
T | TGCTA | MKT
AGC | TAGTT | AAGNotice how the bands are all the same width and contain no ambiguity codes, but they don’t have to all be in the same position in each sequence. In R, we could do this with the following:
patterns <- DNAStringSet(c("TTGCTAMKT", "AGCTAGTTAAG"))
pdict1 <- PDict(patterns, tb.start = 2, tb.width = 4)
pdict1## TB_PDict object of length 2 and variable width (preprocessing algo="ACtree2"):
## - with a head of width 1
## - with a Trusted Band of width 4
## - with a tail of variable width (min=4 / max=6)## TB_PDict object of length 2 and variable width (preprocessing algo="ACtree2"):
## - with a head of variable width (min=1 / max=3)
## - with a Trusted Band of width 5
## - with a tail of width 3We can either define the trusted band relative to the start
(tb.start), or relative to the end (tb.end).
Negative numbers are indices relative to the end of the sequence, so
-4 means “the fourth base counting from the end of the
sequence”.
The trusted band splits the sequences into three parts: the head (before the band), the band, and the tail (part after the band). We can extract these to look at them:
## DNAStringSet object of length 2:
## width seq
## [1] 1 T
## [2] 3 AGC## DNAStringSet object of length 2:
## width seq
## [1] 5 TGCTA
## [2] 5 TAGTT## DNAStringSet object of length 2:
## width seq
## [1] 3 MKT
## [2] 3 AAGGiven a PDict with a trusted band, we can start to do inexact matching. Inexact matching with a PDict that has a trusted band works in the following steps:
This means that an exact match has to be found to the trusted band before any inexact matching takes place. This is a key distinction, and can lead to some unintuitive behavior. For example, suppose we make two trusted bands from a given sequence:
Seq: AATTGGC
Trusted Band 1: AA | TTGG | C
Trusted Band 2: A | ATTG | GCThe trusted band is the center section, and the left and right regions are the head / tail. Now suppose we want to match against a target sequence with at most one mismatch:
Subject: ACTTGGC
Goal: Match with at most 1 mismatch
Trusted Band 1 finds a match:
AA | TTGG | C
AC TTGG C
Trusted Band 2 finds NO match.
A | ATTG | GC
A CTTG GC
(One mismatch, but no exact match in trusted band)Because of this, the trusted band should be defined as the minimal match that must exist for a match to be correct.
To do this in R:
seq <- DNAStringSet(c("AATTGGC"))
# A | ATTG | GC
pdict <- PDict(seq, tb.start=2, tb.end=5)
subjects <- DNAStringSet(c("AATTGGCC", "CATTGGCT", "ACTTGGCG"))
# Pattern: A | ATTG | GC
# Seq 1: A ATTG GCC (matches)
# Seq 2: C ATTG GCT (matches, one mismatch)
# Seq 3: A CTTG GCG (no exact match in TB)
vcountPDict(pdict, subjects, max.mismatch = 1L)## [,1] [,2] [,3]
## [1,] 1 1 0The use of the trusted band allows us to use the other fuzzy matching
techniques we used with the *Pattern functions.
We can match ambiguity codes using fixed:
pattern <- DNAStringSet("ATNC")
# | AT | NC
pdict <- PDict(pattern, tb.start = 1, tb.end = 2)
subject <- DNAString("CCAATATACATNC")
## Only one match, since ambiguities are treated literally
matchPDict(pdict, subject)[[1]]## IRanges object with 1 range and 0 metadata columns:
## start end width
## <integer> <integer> <integer>
## [1] 10 13 4## IRanges object with 2 ranges and 0 metadata columns:
## start end width
## <integer> <integer> <integer>
## [1] 6 9 4
## [2] 10 13 4## Since N matches anything, this matches any AT** pattern
matchPDict(pdict, subject, fixed=FALSE, max.mismatch = 1)[[1]]## IRanges object with 3 ranges and 0 metadata columns:
## start end width
## <integer> <integer> <integer>
## [1] 4 7 4
## [2] 6 9 4
## [3] 10 13 4It’s possible to use the *PDict functions without a
PDict, by just passing in a DNAStringSet. This is slightly slower than
using a preprocessed dictionary, but will be slightly faster than using
matchPattern in a loop, since matchPDict loops
at the C level. Additionally, the preprocessing restrictions for PDicts
only apply to initializing the object, so it’s possible to call
matchPDict for a variable width sequence set if the
sequences aren’t preprocessed:
## Patterns are all different lengths!
patterns <- DNAStringSet(c("AA", "ATGC", "TTGGCC", "T"))
subject <- DNAString("AATTGGCCAATGC")
countPDict(patterns, subject)## [1] 2 1 1 3This is also the only way to use the with.indels = TRUE
argument.
Inexact matching can become extremely slow for large numbers of subjects and patterns, especially if the trusted band is small. The most time consuming part is the inexact matching of the head/tail with the flanking regions of the subject’s matches. Larger trusted bands mean fewer exact matches, which means better performance. The Biostrings documentation recommends aiming for \(\frac{L}{4^w} < 10\), where \(L\) is the number of letters in the subject, and \(w\) is the width of the trusted band. Larger trusted bands will always result in fewer, faster matches, so prioritize making it as large as you confidently can to get the matches you’re looking for.
Biostrings also includes a couple general-purpose helper functions
that operate on matches returned by matching functions (e.g.,
matchPattern, matchPDict).
Fuzzy / inexact matching can match with some mismatches. It’s often
useful to figure out which letters were mismatching in an inexact match
call. For this, we have the mismatch function, which
provides the position of mismatching letters of a given pattern relative
to its matches in the subject:
subject <- DNAString("ATGCCGTA")
pattern <- DNAString("GCGGAA")
matches <- matchPattern(pattern, subject, max.mismatch = 2)
matches## Views on a 8-letter DNAString subject
## subject: ATGCCGTA
## views:
## start end width
## [1] 3 8 6 [GCCGTA]## the pattern matched positions 3-8 of the subject
## pattern: GCGGAA
## subject: ATGCCGTA
## mismatches: | | (positions 3 and 5 of the match)
mismatch(pattern = pattern, x = matches)## [[1]]
## [1] 3 5## [1] 2If we’re matching using ambiguity codes, we can also use the
fixed argument to treat ambiguities in the subject as
ambiguities (and not literal matches)
subject <- DNAString("ATNGTNCTN")
pattern <- DNAString("ATC")
matches <- matchPattern(pattern, subject, fixed = FALSE, max.mismatch = 1L)
matches## Views on a 9-letter DNAString subject
## subject: ATNGTNCTN
## views:
## start end width
## [1] 1 3 3 [ATN]
## [2] 4 6 3 [GTN]
## [3] 5 7 3 [TNC]
## [4] 7 9 3 [CTN]## [[1]]
## [1] 3
##
## [[2]]
## [1] 1 3
##
## [[3]]
## [1] 1 2
##
## [[4]]
## [1] 1 3## [[1]]
## integer(0)
##
## [[2]]
## [1] 1
##
## [[3]]
## [1] 1
##
## [[4]]
## [1] 1Given a set of matches on a subject, we may also be interested in determining how many matches we got at each position. In other words, how well does the pattern “cover” the sequence?
Let’s use the same example from the previous section, including ambiguities:
subject <- DNAString("ATNGTNCTN")
pattern <- DNAString("ATC")
matches <- matchPattern(pattern, subject, fixed = FALSE, max.mismatch = 1L)
matches## Views on a 9-letter DNAString subject
## subject: ATNGTNCTN
## views:
## start end width
## [1] 1 3 3 [ATN]
## [2] 4 6 3 [GTN]
## [3] 5 7 3 [TNC]
## [4] 7 9 3 [CTN]We can imagine overlaying the regions that match the pattern:
Subject: A T N G T N C T N
Match 1: - - -
Match 2: - - -
Match 3: - - -
Match 4: - - -If we then count how many matches we have for each position, we get the coverage:
Subject: A T N G T N C T N
Match 1: - - -
Match 2: - - -
Match 3: - - -
Match 4: - - -
Coverage: 1 1 1 1 2 2 2 1 1In code, this looks like:
subject <- DNAString("ATNGTNCTN")
pattern <- DNAString("ATC")
matches <- matchPattern(pattern, subject, fixed = FALSE, max.mismatch = 1L)
## This is a run-length encoding
coverage(matches)## integer-Rle of length 9 with 3 runs
## Lengths: 4 3 2
## Values : 1 2 1## [1] 1 1 1 1 2 2 2 1 1If we’re matching multiple patterns against a subject,
coverage will give us the number of matches from
any pattern for each position:
patterns <- DNAStringSet(c("ATC", "ATG"))
pd <- PDict(patterns, tb.start = 1, tb.end = 2)
subject <- DNAString("ATATATATATGAC")
countPDict(pd, subject, max.mismatch = 1L)## [1] 5 5## [1] 2 2 4 2 4 2 4 2 4 2 2The previous sections have covered everything you need to know for matching sequences with Biostrings. In this last section, we’ll cover some of the more application-specific functions you can use for biological analyses.
Finding palindromic sequences is a critical analysis for many
analyses. Biostrings provides the findPalindromes function
to identify palindromic regions in a sequence.
By palindromic sequences, we mean two regions of a sequence that are reverse complements of each other, potentially separated by a non-matching region. These palindromic sequences can form hairpin loops, like so:
Palindromic Sequence: GAATTCAATCTAAGAATTC
This forms a hairpin:
A A
G A A T T C T
| | | | | | C
C T T A A G T
A AThis palindrome has three regions: a left arm, a right arm, and a loop. The arms are the two regions that match, with left / right corresponding to which comes first / last in the sequence (respectively). The loop is the middle region that doesn’t bind to itself.
In the previous example, we’d see the following regions:
Palindromic Sequence: GAATTCAATCTAAGAATTC
Separated into groups: GAATTC AATCTAA GAATTC
/ | \
left arm loop right armWe can find this with the findPalindrome function in
Biostrings:
## Views on a 19-letter DNAString subject
## subject: GAATTCAATCTAAGAATTC
## views:
## start end width
## [1] 1 19 19 [GAATTCAATCTAAGAATTC]In this case, we get a single view on the whole string, since the
entire sequence is a palindrome. Note that we set
max.looplength to control the maximum size of loop we look
for. We can also set min.armlength and
max.armlength to control the size of the arms, and
min.looplength to control the size of the loop.
If we have a sequence with multiple palindromes, we can find them all at once:
subject <- DNAString("GAAGTTCTTCATGCATGTAACATG")
# Default parameters:
# - min.armlength = 4
# - max.looplength = 1 (strict palindromes only)
findPalindromes(subject, max.looplength = 3L)## Views on a 24-letter DNAString subject
## subject: GAAGTTCTTCATGCATGTAACATG
## views:
## start end width
## [1] 1 10 10 [GAAGTTCTTC]
## [2] 10 17 8 [CATGCATG]
## [3] 14 24 11 [CATGTAACATG]Now we have three palindromes! However, it’s a little difficult to tell what the palindromic regions are. We can identify the left/right arms with some helper functions:
subject <- DNAString("GAAGTTCTTCATGCATGTAACATG")
palindromes <- findPalindromes(subject, max.looplength = 3L)
## lengths of all arms
palindromeArmLength(palindromes)## [1] 4 4 5## Views on a 24-letter DNAString subject
## subject: GAAGTTCTTCATGCATGTAACATG
## views:
## start end width
## [1] 1 4 4 [GAAG]
## [2] 10 13 4 [CATG]
## [3] 14 18 5 [CATGT]## Views on a 24-letter DNAString subject
## subject: GAAGTTCTTCATGCATGTAACATG
## views:
## start end width
## [1] 7 10 4 [CTTC]
## [2] 14 17 4 [CATG]
## [3] 20 24 5 [ACATG]We can also set allow.wobble to TRUE to
treat G-T or G-U pairings as matches:
subject <- DNAString("ATGCGTAT")
## No matches, because G-T pairings are ignored
findPalindromes(subject)## Views on a 8-letter DNAString subject
## subject: ATGCGTAT
## views: NONE## Views on a 8-letter DNAString subject
## subject: ATGCGTAT
## views:
## start end width
## [1] 1 8 8 [ATGCGTAT]## We can also just allow mismatches directly, which aren't limited to G-T pairs
findPalindromes(subject, max.mismatch = 1L)## Views on a 8-letter DNAString subject
## subject: ATGCGTAT
## views:
## start end width
## [1] 1 8 8 [ATGCGTAT]findPalindromes will search for reverse complement
sequences if the input is a DNAString or RNAString. If the input is a
BString or AAString, it’ll just search for regular palindromes:
## Views on a 39-letter BString subject
## subject: amanaplanacanalpanama racecar momomomom
## views:
## start end width
## [1] 1 21 21 [amanaplanacanalpanama]
## [2] 22 30 9 [ racecar ]
## [3] 31 39 9 [momomomom]Another common workflow is identifying amplicons that match to a
given pair of primers. For this, we have the matchProbePair
function, which takes as input a pair (forward and reverse) of primers
and a target sequence, and returns all the amplicons mapped to it.
Let’s look at an example, using a primer pair derived directly from a yeast chromosome.
options(width = 70) # limiting the width of XString output
data("yeastSEQCHR1") # bundled with the Biostrings package
yeast_chr1 <- DNAString(yeastSEQCHR1)
## We'll target 500 bases starting at position 50,000
region_start <- 50000L
region_width <- 500L
probe_size <- 20L
## Getting some subsequences to act as primers
forward_probe <- subseq(yeast_chr1, start = region_start, width = probe_size)
reverse_probe <- subseq(
yeast_chr1,
end = region_start + region_width - 1,
width = probe_size
)
## The reverse probe should be the reverse complement
reverse_probe <- reverseComplement(reverse_probe)
## Now we can identify the amplicons
matchProbePair(forward_probe, reverse_probe, subject = yeast_chr1)## Views on a 230208-letter DNAString subject
## subject: CCACACCACACCCACACACCCACACACCA...GGGTGTGGTGTGGGTGTGGTGTGTGTGGG
## views:
## start end width
## [1] 50000 50499 500 [AGAGGGTGATTTCTTTCTCGA...CTTTCCCATAGTACCATTAA]Under the hood, matchProbePair uses
matchPattern, which is why it returns an XStringViews
object. That means it supports all the features of
matchPattern:
options(width = 70) # limiting the width of XString output
## Including mismatches
matchProbePair(
forward_probe,
reverse_probe,
yeast_chr1,
max.mismatch = 5
)## Views on a 230208-letter DNAString subject
## subject: CCACACCACACCCACACACCCACACACCA...GGGTGTGGTGTGGGTGTGGTGTGTGTGGG
## views:
## start end width
## [1] 50000 50499 500 [AGAGGGTGATTTCTTTCTCG...CTTTCCCATAGTACCATTAA]
## [2] 85724 87210 1487 [AGAGCCTTATCTCTTTCTAG...CTTTCCCACAGCATCTATAA]
## [3] 91467 134389 42923 [ATATGGTAATTTCTTTCTGT...CCTTCCCATAATAGTTTTAA]## Not sure why you'd want to do this, but you can!
matchProbePair(
forward_probe,
reverse_probe,
yeast_chr1,
min.mismatch = 6,
max.mismatch = 8
)## Views on a 230208-letter DNAString subject
## subject: CCACACCACACCCACACACCCACACACCA...GGGTGTGGTGTGGGTGTGGTGTGTGTGGG
## views:
## start end width
## [1] 603 834 232 [ACATTTTGATATCTATATC...CTATCATGGTATCACTAA]
## [2] 3085 3419 335 [AGACTGCCAAATTTTTCTT...TTTTCAAAGAAACTTCAA]
## [3] 4140 4257 118 [AGAGAGTACTTTTTTTGGA...CTAAAAAATTCACTGTAT]
## [4] 5089 5337 249 [TGAAGATTTGTATATGAAA...TCTAATCAGTTCCATTAA]
## [5] 6956 7097 142 [TTTATAGTAGTAGTGTGAA...AGAAAGATATAAAAATTA]
## ... ... ... ... ...
## [310] 227170 227460 291 [GTAAAAGTATGAAGTGAAA...TTTCCCATTTAACATTTA]
## [311] 227687 227742 56 [TCAAATATAATATGTGCAA...CGAAAGCAATGGCCGCCA]
## [312] 228251 228324 74 [TGAGGACAACTCTTTACTC...CACAAGGTATCACATTAT]
## [313] 229121 229121 1 [A]
## [314] 229186 229328 143 [TTAACTTTACTAAAGGAAA...GTTAAATGGGATCATTAA]We can also find matches in sequences using position weight matrices
(PWMs). This functions exactly like all the other matching functions
we’ve covered, except we’ll use a PWM as the pattern rather than an
XString directly. We can create a PWM with the PWM
function:
## Some possible TATA box binding motifs
motifs <- DNAStringSet(c("TATAAA", "TATATA", "TATAAT", "TATATT"))
## Converting into a PWM
pwm <- PWM(motifs)
pwm## [,1] [,2] [,3] [,4] [,5] [,6]
## A 0.0000000 0.1801462 0.0000000 0.1801462 0.1397077 0.1397077
## C 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## G 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## T 0.1801462 0.0000000 0.1801462 0.0000000 0.1397077 0.1397077Now we can use our countPWM and matchPWM
just like we would with a pattern, we’re just going to use the PWM as
our pattern:
options(width = 70) # limiting the width of XString output
data("yeastSEQCHR1") # bundled with the Biostrings package
yeast_chr1 <- DNAString(yeastSEQCHR1)
## Getting the number of matching motifs
countPWM(pwm, yeast_chr1)## [1] 7085## Views on a 230208-letter DNAString subject
## subject: CCACACCACACCCACACACCCACACACCA...GGGTGTGGTGTGGGTGTGGTGTGTGTGGG
## views:
## start end width
## [1] 316 321 6 [CATATT]
## [2] 342 347 6 [TAAATA]
## [3] 372 377 6 [TTTATA]
## [4] 408 413 6 [TGTATA]
## [5] 410 415 6 [TATACT]
## ... ... ... ... ...
## [7081] 229813 229818 6 [TATCAA]
## [7082] 229826 229831 6 [TATGTA]
## [7083] 229908 229913 6 [GATATA]
## [7084] 229910 229915 6 [TATATT]
## [7085] 230071 230076 6 [TAGAAT]As you can see, there are tons of matches. First, let’s get an idea of how well each of these matches scored:
## finding matches with score included
matches <- matchPWM(pwm, yeast_chr1, with.score = TRUE)
## `mcols` references the metadata columns
all_scores <- mcols(matches)$score
table(round(all_scores, 2))##
## 0.82 0.86 1
## 5867 637 581The high scores are because we have such a small motif, but we do have some subsequences that scored higher than others. Let’s trim out the lower scoring sequences:
## [1] 7085## [1] 1218## [1] 581There are a lot of features with the *PWM functions. You
can use logged probability ratios (the default) or raw probabilities for
scores, set priors for the estimated probability of bases at each
position, and renormalize your PWMs to a unit scale. Check out
?PWM for more details on these features.
If you’re interested in diving more into the details of matching in Biostrings, check out these man pages:
lowlevel-matching: low-level matchingmatchPDict-inexactmatch-utilsMatching in Biostrings is still a work in progress, and there are many features we’re still working on. Some of them are even noted in the man pages. If you’re interested in contributing to Biostrings, feel free to file a pull request!
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets
## [7] methods base
##
## other attached packages:
## [1] affydata_1.61.0 affy_1.91.0 hgu95av2cdf_2.18.0
## [4] hgu95av2probe_2.18.0 AnnotationDbi_1.75.0 Biobase_2.73.1
## [7] Biostrings_2.81.5 Seqinfo_1.3.0 XVector_0.53.0
## [10] IRanges_2.47.2 S4Vectors_0.51.5 BiocGenerics_0.59.10
## [13] generics_0.1.4 BiocStyle_2.41.0
##
## loaded via a namespace (and not attached):
## [1] bit_4.6.0 preprocessCore_1.75.0 jsonlite_2.0.0
## [4] compiler_4.6.1 BiocManager_1.30.27 crayon_1.5.3
## [7] blob_1.3.0 jquerylib_0.1.4 png_0.1-9
## [10] yaml_2.3.12 fastmap_1.2.0 R6_2.6.1
## [13] knitr_1.51 maketools_1.3.2 DBI_1.3.0
## [16] bslib_0.11.0 affyio_1.83.0 rlang_1.3.0
## [19] KEGGREST_1.53.5 cachem_1.1.0 xfun_0.60
## [22] sass_0.4.10 sys_3.4.3 bit64_4.8.2
## [25] otel_0.2.0 RSQLite_3.53.3 memoise_2.0.1
## [28] cli_3.6.6 digest_0.6.39 lifecycle_1.0.5
## [31] vctrs_0.7.3 evaluate_1.0.5 buildtools_1.0.0
## [34] httr_1.4.8 rmarkdown_2.31 tools_4.6.1
## [37] htmltools_0.5.9