## ----setup,include=FALSE------------------------------------------------------
library(BiocStyle)
library(Biostrings)

## -----------------------------------------------------------------------------
DNAString("ATGCATGC")

RNAString("AUGCAUGC")

AAString("ARNDC*")

BString("I am a BString")

## -----------------------------------------------------------------------------
DNAStringSet(c("ATGC", "CGTA", "AAAAA", "TT"))

## -----------------------------------------------------------------------------
filepath1 <- system.file("extdata", "someORF.fa", package="Biostrings")
readDNAStringSet(filepath1)

## -----------------------------------------------------------------------------
x <- DNAStringSet(c("ATGC", "CGTA", "AAAAA", "TT"))

# Subsetting to a DNAStringSet
x[1]

# Subsetting to a DNAString
x[[1]]

## -----------------------------------------------------------------------------
seq <- DNAString("AAATTTGGGCCC")

## Splitting into codons
Views(seq, 
      start = seq(1, length(seq), by = 3), 
      end = seq(3, length(seq), by = 3))

## All overlapping subsequences of length 3
Views(seq, 
      start = seq_len(length(seq) - 2), 
      width = 3)

## Views can also be "out of range":
Views(seq, 
      end = seq_along(seq), 
      width = 3)

## -----------------------------------------------------------------------------
seq <- DNAString("AAATTTGGGCCC")

## `codons` returns an XStringViews 
seq_codons <- codons(seq)

seq_codons

## We can also convert it to a DNAStringSet
DNAStringSet(seq_codons)

## -----------------------------------------------------------------------------
## the first string is the pattern,
## the second is the subject
neditAt(pattern = "AAAAA", subject = "AAAAB")

neditAt("ABCDEFG", "AAADDDG")

## Going from subject -> pattern
## AACGCA -> AACGTA -> AACCTA -> AAGCTA -> GAGCTA: 4 edits
neditAt("GAGCTA", "AACGCA")

## -----------------------------------------------------------------------------
## AACGCA -> GACGCA -> GAGCA -> GAGCTA: 3 edits
neditAt("GAGCTA", "AACGCA", with.indels = TRUE)

## -----------------------------------------------------------------------------
## 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")
neditAt("ABC", "ABCD", with.indels = TRUE)

## However, if we change the subject slightly
## (note just showing the first 3 characters of the subject)
neditAt("ABC", "DBCA") # DBCA -> ABCA: 1
neditAt("ABC", "DCBA") # DCBA -> DCCA -> DBCA -> ABCA : 3
neditAt("ABC", "DCBA", with.indels = TRUE) # DCBA -> ADCBA -> ABCBA: 2

## -----------------------------------------------------------------------------
neditAt("ABC", "ABCDE") == 0

isMatchingAt("ABC", "ABCDE")

## -----------------------------------------------------------------------------
neditAt("ABC", "BABCABC", at = 1)

neditAt("ABC", "BABCABC", at = 2)

neditAt("ABC", "BABCABC", at = 5)

## -----------------------------------------------------------------------------
subject <- "BABCABC"
edit_dists <- neditAt("ABC", subject, at = seq_len(nchar(subject)))
edit_dists

## These are the positions of the matches!
which(edit_dists == 0)

matches <- isMatchingAt("ABC", subject, at = seq_len(nchar(subject)))
which(matches)

## -----------------------------------------------------------------------------
pattern <- "AAB"
subject <- "AABBAAA"
s <- seq_len(nchar(subject))
neditAt(pattern, subject, at = s)

## Default is min.mismatch = max.mismatch = 0
isMatchingAt(pattern, subject, at = s)

## If max.mismatch = 1, we get more matches
isMatchingAt(pattern, subject, at = s, max.mismatch = 1)

## 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)

## -----------------------------------------------------------------------------
neditAt("GAGC", "AACG", with.indels = FALSE)
isMatchingAt("GAGCTA", "AACGCA", with.indels = FALSE, max.mismatch = 3)

neditAt("GAGC", "AACG", with.indels = TRUE)
isMatchingAt("GAGCTA", "AACGCA", with.indels = TRUE, max.mismatch = 3)

## -----------------------------------------------------------------------------
subject <- "BABCABC"
edit_dists <- neditAt("ABC", subject, at = seq_len(nchar(subject)))

which(edit_dists == 0)

matchPattern("ABC", subject)

## -----------------------------------------------------------------------------
pattern <- DNAString("AG")
subject <- DNAString("AGTAGCAGA")

matchPattern(pattern, subject)

## -----------------------------------------------------------------------------
pattern <- DNAString("AGY") # Y = T or G
subject <- DNAString("AGTAGCAGAAGY")

## By default, ambiguities only match ambiguities
## (i.e., Y = Y)
matchPattern(pattern, subject)

## But with fixed = FALSE, it can match other letters
## (i.e., Y = Y, T, or G)
matchPattern(pattern, subject, fixed = FALSE)

## -----------------------------------------------------------------------------
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)

## -----------------------------------------------------------------------------
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

second_matches <- matchPattern(pattern2, first_matches, fixed = FALSE)
second_matches

## -----------------------------------------------------------------------------
pattern <- DNAString("ANN")
subject <- DNAString("ATGTTTAATTTATAATATTTATAAAT")

length(matchPattern(pattern, subject, fixed = FALSE))

countPattern(pattern, subject, fixed = FALSE)

## -----------------------------------------------------------------------------
pattern <- DNAString("ATNG")
subjects <- DNAStringSet(c("ATTGATTAATGGATCG", "AGTGAGACAGATTG"))

vcountPattern(pattern, subjects, fixed = FALSE)
vcountPattern(pattern, subjects, fixed = FALSE, max.mismatch = 2)

vmatchPattern(pattern, subjects, fixed = FALSE, max.mismatch = 1)

## -----------------------------------------------------------------------------
pattern <- DNAString("ATNG")
subjects <- DNAStringSet(c("ATTGATTAATGGATCG", "AGTGAGACAGATTG"))

matches <- vmatchPattern(pattern, subjects, fixed = FALSE)

## 3 matches at indices 1, 9, and 13
matches[[1]]

## Note this is the same as in matchPattern
matchPattern(pattern, subjects[[1]], fixed = FALSE)

## ...and same for the second match
matches[[2]]
matchPattern(pattern, subjects[[2]], fixed = FALSE)

## -----------------------------------------------------------------------------
patterns <- DNAStringSet(c("ATGC", "GGGG", "CCCC", "ATAT"))
pdict <- PDict(patterns)
pdict

## -----------------------------------------------------------------------------
pdict[[1]]

## -----------------------------------------------------------------------------
patterns <- c(seq1 = "ATGC", seq2 = "CGTA")
pdict <- PDict(patterns)
names(pdict)

## -----------------------------------------------------------------------------
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)

## Where are the matches?
matchPDict(pdict, subject)

## -----------------------------------------------------------------------------
## only the first and second patterns have matches
whichPDict(pdict, subject)

## -----------------------------------------------------------------------------
subjects <- DNAStringSet(c(subj1 = "ATGCAAAATTTAATTGC", subj2 = "CCAATCCTACTATGC"))

vwhichPDict(pdict, subjects)

## -----------------------------------------------------------------------------
vcountPDict(pdict, subjects)

## -----------------------------------------------------------------------------
vmatchPDict(pdict, subjects)

## -----------------------------------------------------------------------------
all_matches <- vmatchPDict(pdict, subjects)

## first entry is the matches for all patterns against the first subject:
all_matches[[1]]

## -----------------------------------------------------------------------------
patterns <- DNAStringSet(c("TTGCTAMKT", "AGCTAGTTAAG"))
pdict1 <- PDict(patterns, tb.start = 2, tb.width = 4)
pdict1

pdict2 <- PDict(patterns, tb.end = -4, tb.width = 5)
pdict2

## -----------------------------------------------------------------------------
head(pdict2)
tb(pdict2)
tail(pdict2)

## -----------------------------------------------------------------------------
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)

## -----------------------------------------------------------------------------
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]]

## Two matches; N in tail can match any nucleotide
matchPDict(pdict, subject, fixed=FALSE)[[1]]

## Since N matches anything, this matches any AT** pattern
matchPDict(pdict, subject, fixed=FALSE, max.mismatch = 1)[[1]]

## -----------------------------------------------------------------------------
## Patterns are all different lengths!
patterns <- DNAStringSet(c("AA", "ATGC", "TTGGCC", "T"))
subject <- DNAString("AATTGGCCAATGC")
countPDict(patterns, subject)

## -----------------------------------------------------------------------------
subject <- DNAString("ATGCCGTA")
pattern <- DNAString("GCGGAA")

matches <- matchPattern(pattern, subject, max.mismatch = 2)
matches

## 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)

## How many mismatches were there overall?
nmismatch(pattern = pattern, x = matches)

## -----------------------------------------------------------------------------
subject <- DNAString("ATNGTNCTN")
pattern <- DNAString("ATC")

matches <- matchPattern(pattern, subject, fixed = FALSE, max.mismatch = 1L)
matches

## mismatches if we ignored ambiguities in subject
mismatch(pattern, matches)

## mismatches if we use ambiguities in subject
mismatch(pattern, matches, fixed = FALSE)

## -----------------------------------------------------------------------------
subject <- DNAString("ATNGTNCTN")
pattern <- DNAString("ATC")

matches <- matchPattern(pattern, subject, fixed = FALSE, max.mismatch = 1L)
matches

## -----------------------------------------------------------------------------
subject <- DNAString("ATNGTNCTN")
pattern <- DNAString("ATC")

matches <- matchPattern(pattern, subject, fixed = FALSE, max.mismatch = 1L)

## This is a run-length encoding
coverage(matches)

## converting to a vector, we get:
as.vector(coverage(matches))

## -----------------------------------------------------------------------------
patterns <- DNAStringSet(c("ATC", "ATG"))
pd <- PDict(patterns, tb.start = 1, tb.end = 2)
subject <- DNAString("ATATATATATGAC")

countPDict(pd, subject, max.mismatch = 1L)

matches <- matchPDict(pd, subject, max.mismatch = 1L)
as.vector(coverage(matches))

## -----------------------------------------------------------------------------
subject <- DNAString("GAATTCAATCTAAGAATTC")
findPalindromes(subject, max.looplength = 7)

## -----------------------------------------------------------------------------
subject <- DNAString("GAAGTTCTTCATGCATGTAACATG")

# Default parameters:
# - min.armlength = 4
# - max.looplength = 1 (strict palindromes only)
findPalindromes(subject, max.looplength = 3L)

## -----------------------------------------------------------------------------
subject <- DNAString("GAAGTTCTTCATGCATGTAACATG")
palindromes <- findPalindromes(subject, max.looplength = 3L)

## lengths of all arms
palindromeArmLength(palindromes)

## Get only the left arm
palindromeLeftArm(palindromes)

## Get only the right arm
palindromeRightArm(palindromes)

## -----------------------------------------------------------------------------
subject <- DNAString("ATGCGTAT")

## No matches, because G-T pairings are ignored
findPalindromes(subject)

## Allowing wobble base pairs, we find a match
findPalindromes(subject, allow.wobble = TRUE)

## We can also just allow mismatches directly, which aren't limited to G-T pairs
findPalindromes(subject, max.mismatch = 1L)

## -----------------------------------------------------------------------------
subject <- BString("amanaplanacanalpanama racecar momomomom")
findPalindromes(subject)

## -------------------------------------------------------------------
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)

## -------------------------------------------------------------------
options(width = 70) # limiting the width of XString output

## Including mismatches
matchProbePair(
  forward_probe,
  reverse_probe,
  yeast_chr1,
  max.mismatch = 5
)

## 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
)

## -----------------------------------------------------------------------------
## Some possible TATA box binding motifs
motifs <- DNAStringSet(c("TATAAA", "TATATA", "TATAAT", "TATATT"))

## Converting into a PWM
pwm <- PWM(motifs)
pwm

## -------------------------------------------------------------------
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)

## Finding the matches
matchPWM(pwm, yeast_chr1)

## -----------------------------------------------------------------------------
## 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))

## -----------------------------------------------------------------------------
## Way too many matches
countPWM(pwm, yeast_chr1)

## ...but we can limit the matches with min.score
countPWM(pwm, yeast_chr1, min.score = "86%")

## Score can also be expressed as a decimal value
countPWM(pwm, yeast_chr1, min.score = 0.95)

## ----sessionInfo, echo=FALSE--------------------------------------------------
sessionInfo()