What can `Biostrings’ do?


Loading the library

The package Biostrings is hosted on Bioconductor:

# source("https://bioconductor.org/biocLite.R")
# biocLite("Biostrings")
suppressMessages(library(Biostrings))   # load this library silently  


Working with DNA sequences

fastafile = "seqDNA.fasta"
dna_seqs <- readDNAStringSet(fastafile)
class(dna_seqs)
## [1] "DNAStringSet"
## attr(,"package")
## [1] "Biostrings"
# methods(class = "DNAStringSet")
names(dna_seqs)
## [1] "AK002358.PE1            501 residues"
## [2] "HSU78678.PE1            501 residues"
## [3] "RNU73525.PE1            501 residues"
nchar(dna_seqs)
## [1] 501 501 501
sum(nchar(dna_seqs))   # total number of bases 
## [1] 1503
dna_seqs[1] == dna_seqs[2] 
## [1] FALSE
# consensusMatrix(dna_seqs) 
toString(dna_seqs[1])
## [1] "ATGGCTCAGCGGCTCCTCCTGGGGAGGTTCCTGACCTCAGTCATCTCCAGGAAGCCTCCTCAGGGTGTGTGGGCTTCCCTCACCTCTAAGACCCTGCAGACCCCTCAGTACAATGCTGGTGGTCTAACAGTAATGCCCAGCCCAGCCCGGACAGTACACACCACCAGAGTCTGTTTGACGACCTTTAACGTCCAGGATGGACCTGACTTTCAAGACAGAGTTGTCAACAGTGAGACACCAGTTGTTGTGGACTTTCATGCACAGTGGTGTGGCCCCTGCAAGATCCTAGGACCGCGGCTAGAGAAGATGGTCGCCAAGCAGCACGGGAAGGTGGTCATGGCCAAAGTGGACATTGACGATCACACAGACCTTGCCATTGAATATGAGGTGTCAGCTGTGCCTACCGTGCTAGCCATCAAGAACGGGGACGTGGTGGACAAGTTTGTGGGGATCAAGGACGAGGACCAGCTAGAAGCCTTCCTGAAGAAGCTGATTGGCTGA"
consensusString(dna_seqs)
## [1] "ATGGCTCAGCGRCTYCTYCTGRGGAGGTTCCTGRCCTCWGTCATCTCCAGGAAGCCYYCTCAGGGTSWGTGGSCWYCCCTCACYTCYAMRRSCCTGCAGACCCCWCMRTRCARTSCTGGTGGYCTRACWGKAAYRCCCARCCCWGCCCGGACADTWYACACCACSAGRRTCTSYTYRACRACCTTTAAYRTCCAGGATGGACCTGACTTTCAAGACMGAGTKGTCAACAGTGAGACACCAGTKGTYGTGGAYTTYCAYGCACAGTGGTGTGGMCCCTGCAAGATCCTRGGRCCKMGGYTAGAGAAGATGGTVGCCAARCAGCACGGGAAGGTGGTSATGGCCAARGTGGAYATTGAYGAYCACACAGACCTYGCCATTGARTAYGAGGTGTCWGCKGTGCCYACYGTGCTRGCCATSAAGAAYGGGGACGTGGTGGACAAGTTTGTGGGSATCAAGGAYGARGAYCAGYTRGARGCCTTCCTGAAGAAGCTRATTGGCTGA"
trinucleotideFrequency(dna_seqs)  
##      AAA AAC AAG AAT ACA ACC ACG ACT AGA AGC AGG AGT ATA ATC ATG ATT CAA
## [1,]   1   4  14   3  13  13   7   2  14  11   9  11   1   5   8   3   9
## [2,]   1   4  13   4  11  12   4   5  10   8  13   9   4   6  10   4  12
## [3,]   2   8  13   2  14  13   8   3  13   9   9   9   0   5   6   4  11
##      CAC CAG CAT CCA CCC CCG CCT CGA CGC CGG CGT CTA CTC CTG CTT GAA GAC
## [1,]   9  20   6  12   9   3  18   3   2   5   3   7  10  11   6   8  19
## [2,]  13  15   4  16  11   3  15   4   2   3   1   0   8  15   7   7  14
## [3,]  10  15   6  11   8   3  19   4   0   5   5   5  10  15   7   9  16
##      GAG GAT GCA GCC GCG GCT GGA GGC GGG GGT GTA GTC GTG GTT TAA TAC TAG
## [1,]   7   6   5  11   2  10  14   7   9  10   3   8  17   6   4   3   4
## [2,]  11  11   7  13   2   4  15  11   7  14   2   5  19   5   2   1   1
## [3,]   9   6   4  11   1   9  15   7   7  11   3   7  18   6   3   4   3
##      TAT TCA TCC TCG TCT TGA TGC TGG TGT TTA TTC TTG TTT
## [1,]   1  14   9   1   4  10   8  17  10   1   5   9   5
## [2,]   4  10   9   1   6  15   5  24   7   2   7   7   4
## [3,]   0  13   9   2   6   9   9  19   9   2   8   7   5
reverseComplement(dna_seqs)
##   A DNAStringSet instance of length 3
##     width seq                                          names               
## [1]   501 TCAGCCAATCAGCTTCTTCAG...AGGAGGAGCCGCTGAGCCAT AK002358.PE1     ...
## [2]   501 TCAGCCAATCAGCTTCTTCAG...AGAAGAAGTCGCTGAGCCAT HSU78678.PE1     ...
## [3]   501 TCAGCCAATTAGCTTCTTCAG...AGGAGAAGCCGCTGAGCCAT RNU73525.PE1     ...
translate(dna_seqs)
##   A AAStringSet instance of length 3
##     width seq                                          names               
## [1]   167 MAQRLLLGRFLTSVISRKPPQ...FVGIKDEDQLEAFLKKLIG* AK002358.PE1     ...
## [2]   167 MAQRLLLRRFLASVISRKPSQ...FVGIKDEDQLEAFLKKLIG* HSU78678.PE1     ...
## [3]   167 MAQRLLLRRFLTSVISRKPPQ...FVGIKDEDQLEAFLKKLIG* RNU73525.PE1     ...
pairwiseAlignment(dna_seqs[1], dna_seqs[2])
## Global PairwiseAlignmentsSingleSubject (1 of 1)
## pattern: [1] ATGGCTCAGCGGCTCCTCCTGGGGAGGTTCCT...TAGAAGCCTTCCTGAAGAAGCTGATTGGCTGA 
## subject: [1] ATGGCTCAGCGACTTCTTCTGAGGAGGTTCCT...TGGAGGCCTTCCTGAAGAAGCTGATTGGCTGA 
## score: 449.0658


Working with amino acid sequences

aa_seqs <- readAAStringSet("seqAA.fasta")
class(aa_seqs)
## [1] "AAStringSet"
## attr(,"package")
## [1] "Biostrings"
# methods(class = "AAStringSet")   # see what you can do with objects of this class 
alphabetFrequency(aa_seqs)
##      A  R N D C Q  E G H I  L  K M F  P  S T W Y V U O B J Z X * - + .
## [1,] 8 13 7 6 6 9 18 8 1 9 29 14 5 6 10 16 7 3 1 6 0 0 0 0 0 0 1 0 0 0
##      other
## [1,]     0
names(aa_seqs)
## [1] "A06852                  183 residues"
nchar(aa_seqs)
## [1] 183
toString(aa_seqs)
## [1] "MPRLFSYLLGVWLLLSQLPREIPGQSTNDFIKACGRELVRLWVEICGSVSWGRTALSLEEPQLETGPPAETMPSSITKDAEILKMMLEFVPNLPQELKATLSERQPSLRELQQSASKDSNLNFEEFKKIILNRQNEAEDKSLLELKNLGLDKHSRKKRLFRMTLSEKCCQVGCIRKDIARLC*"
vmatchPattern("SYLL", aa_seqs)
## MIndex object of length 1
## $`A06852                  183 residues`
## IRanges object with 1 range and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         6         9         4