R Lesson #14 - Introduction to Biostrings
For this lesson we are going to be introduced to the Biostrings package, which is a popular way to work with biological sequences in R. First, we must ensure Biostrings is installed before we can load the library.Show output
1library(Biostrings) # load the Biostrings package
The Biostrings package introduces XStringSets to store biological sequences. Amino acid sequences are stored in the class AAStringSet, nucleotide sequences are stored in the class DNAStringSet or RNAStringSet. XStringSets can be imported from FASTA/FASTQ files or created directly with the appropriate constructor.
2dna <- DNAStringSet(c("ACTG", "NSYR", "CG"))
An XStringSet is simply a collection of XStrings. Individual XStrings (i.e., single sequences) can be accessed with `[[`, whereas `[` returns an XStringSet.
345dna[[1]]dna[1]dna[1:2]
Biostrings uses a byte encoding to store nucleotide sequences. The letter A is stored as 0000 0001, C as 0000 0010, G as 0000 0100, and T/U as 0000 1000. Since RNA and DNA sequences are stored in the same manner, they can be equal even though they represent different letters.
678rna <- RNAStringSet(dna)rnadna==rna
Byte encoding allows ambiguity codes to be supported as combinations of letters. For example, the letter N, meaning any nucleotide, would be stored as 0000 1111. The letter S, meaning C or G, is stored as 0000 0110. Biostrings uses the IUPAC standard to represent ambiguity codes.
9IUPAC_CODE_MAP
The advantage of storing nucleotides in this manner is that they can more quickly be compared with binary operations. For example, the AND operator could be used to test whether N matches C (i.e., when fixed = FALSE), whereas the equality operator could be use to test whether N is C (i.e., when fixed = TRUE).
1011neditAt(dna[[1]], dna[[2]], fixed=FALSE) # 0 mismatchesneditAt(dna[[1]], dna[[2]], fixed=TRUE) # 4 mismatches
Other positions of the byte are used to store gap ("-": 0001 0000), mask ("+": 0010 0000), or unknown characters (".": 0100 0000). This allows XStringSets to contain the full set of symbols commonly used in sequences and alignments.
1213align <- DNAStringSet(c(".+AC-TG.", ".+ACCTG."))align
Many standard functions work with XStringSet inputs.
14151617length(align)unique(align)names(align) <- c("one", "two")align
There are also many functions specific to XStringSets.
181920alphabet(align)width(align) # analogous to nchar()alphabetFrequency(align)
Biostrings has the ability to read sequences from popular file formats such as FASTA and FASTQ. It will also read (gzip) compressed files. For example, we can import the sequences from an example set of fly DNA sequences that are included in the Biostrings package.
21222324-252627-28fas1 <- system.file("extdata", "dm3_upstream2000.fa.gz", package="Biostrings")flyDNA <- readDNAStringSet(fas1) # This set contains 26,454 sequences that originated fromthe 2,000 nucleotides upstream of annotated transcriptionstart sites in the reference D. melanogaster genome. flyDNA
One of the advantages of working with XStringSets is that they are very efficient to manipulate relative to character vectors. For example, we could count all of the 3-mer subsequences in the XStringSet nearly instantaneously.
293031oligonucleotideFrequency(flyDNA, width=3, simplify.as="collapse")
The way that XStringSets are internally stored allows them to be subset with almost no memory footprint. For example, we could select the first 100 nucleotides from a random subset of the sequences, and none of the nucleotide information would actually be replicated in memory! This operation is analogous to using substring() but far more memory efficient.
3233s <- sample(length(flyDNA), 1000)subseq(flyDNA[s], 1, 100)
Biostrings is built on top of the IRanges package, which is used to define ranges in a sequence. To effectively use Biostrings, it is necessary to understand how to create and access IRanges objects. Here we will extract the first 100 nucleotides using IRanges.
3435at1 <- IRanges(start=1, end=100)extractAt(flyDNA, at1) # DNAStringSetList
This command returned a DNAStringSetList, which is a list of XStringSets. Like other lists, it can be unlisted to merge its contents into a single XStringSet. Note that XStringSets can be unlisted to become a single XString.
3637unlist(extractAt(flyDNA, at1)) # DNAStringSetunlist(unlist(extractAt(flyDNA, at1))) # DNAString
If we wished to extract different positions from every sequence, we would need to provide an IRangesList rather than a single set of IRanges. For example, we could extract a random nucleotide from every sequence.
3839404142at2 <- sample(100, length(flyDNA), replace=TRUE)at2 <- lapply(at2, IRanges, width=1)at2 <- as(at2, "IRangesList")at2unlist(extractAt(flyDNA, at2))
The Biostrings package is particularly useful for rapidly searching sequences for a specific pattern. The following command finds the pattern "ACTG" in every sequence and returns a list of IRanges.
434445v1 <- vmatchPattern("ACTG", flyDNA)v1unlist(extractAt(flyDNA, v1))
Patterns can match ambiguities (when fixed = FALSE), and allow for mismatches (when max.mismatches > 0). However, this also matches ambiguities in the subject (flyDNA) such that runs of Ns appear as matches. We can change this behavior by only allow ambiguities in the pattern sequence (fixed = "subject").
4647484950v2 <- vmatchPattern(pattern="AACTGNCATCG", subject=flyDNA, fixed="subject", max.mismatch=1)unlist(extractAt(flyDNA, v2))
An analogous function, replaceAt, exist to replace ranges in an XStringSet. For example, we could delete the matches by replacing them with nothing (the default), or we could replace them with another sequence.
5152replaceAt(flyDNA, v2) # deletes the matchesreplaceAt(flyDNA, v2, "+") # replaces the matches
Several Biostrings functions make use of its excellent searching abilities. The findPalindromes function will search for palindromic sequences in an XString. Unfortunately, it does not work with XStringSets, so we must first unlist our XStringSet to convert it into an XString. Then we can try to find potential stem-loop structures in the regions upstream of transcribed regions in the fly genome.
5354555657unlist(flyDNA)pal <- findPalindromes(unlist(flyDNA), min.armlength=8, max.looplength=10, max.mismatch=1)
Now we must map each palindrome back to its respective sequence positions in the XStringSet. Here we can tally the number of times that we saw an upstream palindrome to see if there are any palindrome hotspots.
58596061626364656667686970717273bounds <- c(1, cumsum(width(flyDNA)))head(bounds) # boundaries of each sequencestarts <- start(pal)ends <- end(pal)cbind(head(starts), head(ends)) # pairs of positionsseq <- 1L # current sequencecounts <- integer(2000) # counts for upstream positionsfor (i in seq_along(starts)) { while (ends[i] >= bounds[seq + 1L]) { seq <- seq + 1L } if (starts[i] >= bounds[seq]) { range <- starts[i]:ends[i] - bounds[seq] + 1L counts[range] <- counts[range] + 1L }}
Now we can plot the palindrome counts to see that there are about twice as many palindromes as usual about 200 nucleotides upstream of transcript start sites in the fly genome.
7475767778plot(-rev(seq_along(counts)), counts, type="l", xlab="Positions from transcription start site", ylab="Number of palindromes covering position")
Biostrings also includes a number of functions for working with amino acid sequences, i.e. AAStringSets. As an example, we can import some Yeast open reading frames (ORFs) from a FASTA file included with Biostrings.
7980818283fas2 <- system.file("extdata", "someorf.fa.gz", package="Biostrings")ORFs <- readDNAStringSet(fas2)ORFs
These seven sequences contain the open reading frame plus 1,000 nucleotides to either side. We can remove these flanking nucleotides with subseq.
8485ORFs <- subseq(ORFs, 1001, -1001)ORFs
We can now translate the open reading frames to convert the nucleotides into amino acid sequences.
86aa <- translate(ORFs)
Open reading frames number 1 and 3 contain introns, so we will find more than the expected stop codon at the very end of the sequence.
87vmatchPattern("*", aa)
< Previous Lesson