bch441-work-abc-units/RPR-Biostrings.R

# RPR-Biostrings.R
#
# Purpose:  A Bioinformatics Course:
#              R code accompanying the RPR-Biostrings unit.
#
# Version:  1.1
#
# Date:     2017  10  -  2019  01
# Author:   Boris Steipe (boris.steipe@utoronto.ca)
#
# Versions:
#           1.1    Change from require() to requireNamespace(),
#                      use <package>::<function>() idiom throughout,
#                      use Biocmanager:: not biocLite()
#           1.0    2017 Revisions
#           0.1    First code copied from 2016 material.
#
#
# TODO:
#
#
# == DO NOT SIMPLY  source()  THIS FILE! =======================================
#
# If there are portions you don't understand, use R's help system, Google for an
# answer, or ask your instructor. Don't continue if you don't understand what's
# going on. That's not how it works ...
#
# ==============================================================================


#TOC> ==========================================================================
#TOC>
#TOC>   Section  Title                                           Line
#TOC> ---------------------------------------------------------------
#TOC>   1        The Biostrings Package                            55
#TOC>   2        Getting Data into Biostrings Objects              86
#TOC>   3        Working with Biostrings Objects                  108
#TOC>   3.1        Properties                                     125
#TOC>   3.2        Subsetting                                     163
#TOC>   3.3        Operators                                      175
#TOC>   3.4        Transformations                                182
#TOC>   4        Getting Data out of Biostrings Objects           189
#TOC>   5        More                                             198
#TOC>   5.1        Views                                          200
#TOC>   5.2        Iranges                                        214
#TOC>   5.3        StringSets                                     220
#TOC>
#TOC> ==========================================================================


# This is a very brief introduction to the biostrings package, other units will
# be using more of the biostrings functions.


# =    1  The Biostrings Package  ==============================================


# First, we install and load the Biostrings package from bioconductor

if (! requireNamespace("BiocManager", quietly = TRUE)) {
  install.packages("BiocManager")
}
if (! requireNamespace("Biostrings", quietly = TRUE)) {
  BiocManager::install("Biostrings")
}
# Examine the package information:
library(help = Biostrings)       # basic information
browseVignettes("Biostrings")    # available vignettes
data(package = "Biostrings")     # available datasets


# At its core, Biostrings objects are "classes" of type XString (you can think
# of a "class" in R as a special kind of list), that can take on particular
# flavours for RNA, DNA or amino acid sequence information.

class(Biostrings::RNAString("AUG"))
class(Biostrings::DNAString("ATG"))
class(Biostrings::AAString("M"))

# An essential property of Biostrings objects is that they only allow letters
# from the applicable IUPAC alphabet:
Biostrings::RNAString("AUG")
Biostrings::DNAString("AUG")  # Error! No "U" in IUPAC DNA codes


# =    2  Getting Data into Biostrings Objects  ================================


# Example: read FASTA. Extract sequence. Convert to DNAString object.
rawSeq <- readLines("./data/S288C_YDL056W_MBP1_coding.fsa")
rawSeq <- dbSanitizeSequence(rawSeq)
biosDNAseq <- Biostrings::DNAString(rawSeq) # converts the nucleotide sequence
                                            # into an object of class DNAstring

# Multi FASTA files can be read directly as a "XStringSet) ...
rawMFAfile <- "./data/S288C_YDL056W_MBP1_coding.fsa"
(biosDNASet <- Biostrings::readDNAStringSet(rawMFAfile))

# ... and if you subset one sequence from the set, you get an XString object
# back again.
(Xseq <- biosDNASet[[1]])

biosDNAseq == Xseq           # the comparison evaluates to TRUE ...
identical(biosDNAseq, Xseq)  # ... and indeed the objects are deemed identical.



# =    3  Working with Biostrings Objects  =====================================

# Biostrings is a highly engineered package that is tightly integrated into
# the Bioconductor world - unfortunately that brings with it a somewhat
# undesirable level of computational overhead and dependencies. Using the
# package as we normally do - i.e. calling required functions with their
# explicit package prefix is therefore not advisable. There are generics
# that won't be propery dispatched. If you only need a small number of
# functions for a very specific context, you will probably get away with
# Biostrings::<function>() - but even in the demonstration code of this script
# not everything works out of the box. We'll therefore load the library,
# but we'll (redundantly) use the prefix anyway so as to emphasize where
# the functions come from.

library(Biostrings)


# ==   3.1  Properties  ========================================================
str(rawSeq)
str(biosDNAseq)

length(rawSeq)       # ... is 1: one string only. To get the number of
                     # characters in a string, you need nchar().
length(biosDNAseq)   # but the length of a "Bstring" is the number of elements
nchar(rawSeq)
nchar(biosDNAseq)    # ... but nchar() works too.

(uL <- Biostrings::uniqueLetters(biosDNAseq))

# Count frequencies - with strings, you would strsplit() into a character
# vector and then use table(). biost
Biostrings::alphabetFrequency(biosDNAseq)

# letterFrequency() works with a defined alphabet - such as what uniqueLetters()
# returns.
Biostrings::letterFrequency(biosDNAseq, uL)
sum(Biostrings::letterFrequency(biosDNAseq, c("G", "C"))) /
  length(biosDNAseq) # GC contents

Biostrings::dinucleotideFrequency(biosDNAseq)
barplot(sort(Biostrings::dinucleotideFrequency(biosDNAseq)), cex.names = 0.5)

(triNuc <- Biostrings::trinucleotideFrequency(biosDNAseq))
barplot(sort(triNuc), col="#4499EE33")
triNuc[triNuc == max(triNuc)]
triNuc[triNuc == min(triNuc)]
max(triNuc) / min(triNuc)  # AAA is more than 13 times as frequent as CGT

# compare to a shuffled sequence:
(triNuc <- Biostrings::trinucleotideFrequency(sample(biosDNAseq)))
barplot(sort(triNuc), col="#EEEE4433", add = TRUE)

# Interpret this plot.


# ==   3.2  Subsetting  ========================================================

# Subsetting any XString object works as expected:
biosDNAseq[4:15]

# ... well - maybe not expected, because rawSeq[4:15] would not work.

# Alternatively to the "[" operator, use the subseq() function - especially for
# long sequences. This is far more efficient.
Biostrings::subseq(biosDNAseq, start = 1, end = 30)


# ==   3.3  Operators  =========================================================

# RNAstring() and DNAstring() objects compare U and T as equals!
  Biostrings::RNAString("AUGUCUAACCAAAUAUACUCAGCGAGAUAU") ==
  Biostrings::DNAString("ATGTCTAACCAAATATACTCAGCGAGATAT")


# ==   3.4  Transformations  ===================================================

biosDNAseq[4:15]
Biostrings::reverseComplement(biosDNAseq[4:15])
Biostrings::translate(biosDNAseq[4:15])


# =    4  Getting Data out of Biostrings Objects  ==============================

# If you need a character object, use toString():

Biostrings::toString(biosDNAseq[4:15])

# save() and load() works like on all other R objects.


# =    5  More  ================================================================

# ==   5.1  Views  =============================================================

# Biostring "Views" are objects that store multiple substrings of one
# Biostring object.

(myView <- Biostrings::Views(biosDNAseq,
                             start = c(1, 19, 37),
                             end = c(15, 30, 45)))

# Views are convenient to store feature annotations
names(myView) <- c("Feature-A", "Feature-B", "Feature-C")
cat(sprintf("\n%s\t(%d)\t%s", names(myView), width(myView), myView ))


# ==   5.2  Iranges  ===========================================================

# Biostrings Iranges are like Views with a common start point. These can be
# useful for feature annotations. Instead of start/end you store start/width.


# ==   5.3  StringSets  ========================================================

# Biostring "StringSets" store multiple sequences.
#
ompA <- Biostrings::AAString("MKKTAIAIAVALAGFATVAQA")
sample(ompA) # sample can work directly on a Biostring object to shuffle it

x <- Biostrings::toString(ompA)
for (i in 2:10) {
  x[i] <- Biostrings::toString(sample(ompA))
}
shuffledPeptideSet <- Biostrings::AAStringSet(x)
names(shuffledPeptideSet) <- c("ompA", paste("shuffle.", 1:9, sep=""))
shuffledPeptideSet

length(shuffledPeptideSet)
Biostrings::width(shuffledPeptideSet)
Biostrings::alphabetFrequency(shuffledPeptideSet)


# [END]