The package Biostrings
is hosted on Bioconductor:
# source("https://bioconductor.org/biocLite.R")
# biocLite("Biostrings")
suppressMessages(library(Biostrings)) # load this library silently
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
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