SequenceExtraction

In BRGenomics: Tools for the Efficient Analysis of High-Resolution Genomics Data

In this section, we'll give one more example showing the benefit of Bioconductor integration by using the Biostrings package to extract sequences given by GRanges objects.

library(BRGenomics)
library(Biostrings)

We've included a twobit file of sequences, although users can use fasta files, as well.

# get path to included 2bit file
sfile <- system.file("extdata", "dm6_chr4chrM.2bit",
                     package = "BRGenomics")

We could import the entire sequence using the rtracklayer import() function, which will figure out the file format and import a suitable object. In this case, a DNAStringSet:

seq_chr4 <- import(sfile)
seq_chr4

We included mitochondrial DNA to demonstrate how the DNAStringSet treats multiple chromosomes.

However, we don't need to import all the sequences. Instead, we can make a TwoBitFile object that points to the file, and extract desired sequences from it directly using the getSeq() function:

data("txs_dm6_chr4")
txs_pr <- promoters(txs_dm6_chr4, 0, 100)

seq_txs_pr <- getSeq(TwoBitFile(sfile), txs_pr)
seq_txs_pr

The sequences are stranded as well, such that if a plus and minus strand gene overlapped perfectly, the minus strand sequence would be the reverse complement of the plus strand sequence.

The Biostrings package itself is richly featured, and we'll demonstrate only a couple functions below. This functionality is extended by packages like ggseqlogo, for example, which plots sequence logos directly from DNAStringSets.

RNAStringSet(seq_txs_pr)
suppressWarnings(translate(seq_txs_pr))
oligonucleotideFrequency(seq_txs_pr[1:5], width = 1)
oligonucleotideFrequency(seq_txs_pr[1:5], width = 2)

tss_seq <- getSeq(TwoBitFile(sfile), promoters(txs_dm6_chr4, 4, 4))
tsspwm <- PWM(tss_seq)
tsspwm

BRGenomics documentation built on Nov. 8, 2020, 8:03 p.m.