bch441-work-abc-units/BIN-ALI-Optimal_sequence_alignment.R

# BIN-ALI-Optimal_sequence_alignment.R
#
# Purpose:  A Bioinformatics Course:
#              R code accompanying the BIN-ALI-Optimal_sequence_alignment unit.
#
# Version:  1.3
#
# Date:     2017  09   -   2017  11
# Author:   Boris Steipe (boris.steipe@utoronto.ca)
#
# Versions:
#           1.3    Updated confirmation task with correct logic
#           1.2    Added missing load of seqinr package
#           1.1    Update annotation file logic - it could already have been
#                    prepared in the BIN-FUNC-Annotation unit.
#           1.0.1  bugfix
#           1.0    First 2017 live version.
#           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        Prepare                                                50
#TOC>   2        Biostrings Pairwise Alignment                          71
#TOC>   2.1      Optimal global alignment                               89
#TOC>   2.2      Optimal local alignment                               152
#TOC>   3        APSES Domain annotation by alignment                  176
#TOC>   4        Update your database script                           257
#TOC>   4.1      Preparing an annotation file ...                      263
#TOC>   4.1.1    If you HAVE NOT done the BIN-FUNC-Annotation unit     265
#TOC>   4.1.2    If you HAVE done the BIN-FUNC-Annotation unit         308
#TOC>   4.2      Execute and Validate                                  332
#TOC> 
#TOC> ==========================================================================


# =    1  Prepare  =============================================================

# To simplify code a bit, we will use seqinr's function s2c(x) to make
# character vectors from sequence strings below, rather than the lengthier
# base idiom unlist(strsplit(x, "").
if (!require(seqinr)) {
  install.packages("seqinr")
  library(seqinr)
}
# Package information:
#  library(help = seqinr)       # basic information
#  browseVignettes("seqinr")    # available vignettes
#  data(package = "seqinr")     # available datasets


# You need to recreate the protein database that you have constructed in the
# BIN-Storing_data unit.

source("makeProteinDB.R")


# =    2  Biostrings Pairwise Alignment  =======================================


if (!require(Biostrings, quietly=TRUE)) {
  if (! exists("biocLite")) {
    source("https://bioconductor.org/biocLite.R")
  }
  biocLite("Biostrings")
  library(Biostrings)
}
#  library(help = Biostrings)       # basic information
#  browseVignettes("Biostrings")    # available vignettes
#  data(package = "Biostrings")     # available datasets


# Biostrings stores sequences in "XString" objects. Once we have converted our
# target sequences to AAString objects, the alignment itself is straightforward.

# ==   2.1  Optimal global alignment  ==========================================

# The pairwiseAlignment() function was written to behave
# exactly like the functions you encountered on the EMBOSS server.

# First: make AAString objects ...
sel <- myDB$protein$name == "MBP1_SACCE"
aaMBP1_SACCE <- AAString(myDB$protein$sequence[sel])

sel <- myDB$protein$name == paste("MBP1_", biCode(MYSPE), sep = "")
aaMBP1_MYSPE <-   AAString(myDB$protein$sequence[sel])

?pairwiseAlignment
# ... and align.
# Global optimal alignment with end-gap penalties is default.
ali1 <-  pairwiseAlignment(
  aaMBP1_SACCE,
  aaMBP1_MYSPE,
  substitutionMatrix = "BLOSUM62",
  gapOpening = 10,
  gapExtension = 0.5)

str(ali1)  # ... it's complicated

# This is a Biostrings alignment object. But we can use Biostrings functions to
# tame it:
ali1
writePairwiseAlignments(ali1)   # That should look familiar

# And we can make the internal structure work for us  (@ is for classes as
# $ is for lists ...)
str(ali1@pattern)
ali1@pattern
ali1@pattern@range
ali1@pattern@indel
ali1@pattern@mismatch

# or work with "normal" R functions
# the alignment length
nchar(ali1@pattern)

# the number of identities
sum(s2c(as.character(ali1@pattern)) ==
    s2c(as.character(ali1@subject)))

# ... e.g. to calculate the percentage of identities
100 *
  sum(s2c(as.character(ali1@pattern)) ==
      s2c(as.character(ali1@subject))) /
  nchar(ali1@pattern)
# ... which should be the same as reported in the writePairwiseAlignments()
# output. Awkward to type? Then it calls for a function:
#
percentID <- function(al) {
  # returns the percent-identity of a Biostrings alignment object
  return(100 *
           sum(s2c(as.character(al@pattern)) ==
               s2c(as.character(al@subject))) /
           nchar(al@pattern))
}

percentID(ali1)

# ==   2.2  Optimal local alignment  ===========================================

# Compare with local optimal alignment (like EMBOSS Water)
ali2 <-  pairwiseAlignment(
  aaMBP1_SACCE,
  aaMBP1_MYSPE,
  type = "local",
  substitutionMatrix = "BLOSUM62",
  gapOpening = 50,
  gapExtension = 10)

writePairwiseAlignments(ali2)   # This has probably only aligned the N-terminal
                                # DNA binding domain - but that one has quite
                                # high sequence identity:
percentID(ali2)

# == TASK: ==

# Compare the two alignments. I have weighted the local alignment heavily
# towards an ungapped alignment by setting very high gap penalties. Try changing
# the gap penalties and see what happens: how does the number of indels change,
# how does the length of indels change...


# =    3  APSES Domain annotation by alignment  ================================

# In this section we define the MYSPE APSES sequence by performing a global,
# optimal sequence alignment of the yeast APSES domain with the full length
# protein sequence of the protein that was the most similar to the yeast APSES
# domain.
#

# I have annotated the yeast APSES domain as a feature in the
# database. To view the annotation, we can retrieve it via the proteinID and
# featureID. Here is the yeast protein ID:
(proID <- myDB$protein$ID[myDB$protein$name == "MBP1_SACCE"])


# ... and if you look at the feature table, you can identify the feature ID
(ftrID <- myDB$feature$ID[myDB$feature$name == "APSES fold"])

# ... and with the two annotations we can get the corresponding ID from the
# annotation table
(fanID <- myDB$annotation$ID[myDB$annotation$proteinID == proID &
                             myDB$annotation$featureID == ftrID])

myDB$annotation[myDB$annotation$ID == proID &
                myDB$annotation$ID == ftrID, ]

# The annotation record contains the start and end coordinates which we can use
# to define the APSES domain sequence with a substr() expression.

(start <- myDB$annotation$start[myDB$annotation$ID == fanID])
(end   <- myDB$annotation$end[myDB$annotation$ID == fanID])
(apses <- substr(myDB$protein$sequence[myDB$protein$ID == proID],
                 start,
                 end))

# Lots of code. But don't get lost. Let's recapitulate what we have done: we
# have selected from the sequence column of the protein table the sequence whose
# name is "MBP1_SACCE", and selected from the annotation table the start
# and end coordinates of the annotation that joins an "APSES fold" feature with
# the sequence, and used the start and end coordinates to extract a substring.

# Let's convert this to an AAstring and assign it:
aaMB1_SACCE_APSES <- AAString(apses)

# Now let's align these two sequences of very different length without end-gap
# penalties using the "overlap" type. "overlap" turns the
# end-gap penalties off and that is crucially important since
# the sequences have very different length.

aliApses <-  pairwiseAlignment(
  aaMB1_SACCE_APSES,
  aaMBP1_MYSPE,
  type = "overlap",
  substitutionMatrix = "BLOSUM62",
  gapOpening = 10,
  gapExtension = 0.5)

# Inspect the result. The aligned sequences should be clearly
# homologous, and have (almost) no indels. The entire "pattern"
# sequence from QIYSAR ... to ... KPLFDF  should be matched
# with the "query". Is this correct?
writePairwiseAlignments(aliApses)

# If this is correct, you can extract the matched sequence from
# the alignment object. The syntax is a bit different from what
# you have seen before: this is an "S4 object", not a list. No
# worries: as.character() returns a normal string.
as.character(aliApses@subject)

# Now, what are the aligned start and end coordinates? You can read them from
# the output of writePairwiseAlignments(), or you can get them from the range of
# the match.

str(aliApses@subject@range)

# start is:
aliApses@subject@range@start

# ... and end is:
aliApses@subject@range@start + aliApses@subject@range@width - 1


# =    4  Update your database script  =========================================


# Since we have this feature defined now, we can create a feature annotation
# right away and store it in myDB.

# ==   4.1  Preparing an annotation file ...  ==================================
#
# ===  4.1.1  If you HAVE NOT done the BIN-FUNC-Annotation unit
#
#
#   You DON'T already have a file called "<MYSPE>-Annotations.json" in the
#   ./data/ directory:
#
#   - Make a copy of the file "./data/refAnnotations.json" and put it in your
#     project directory.
#
#   - Give it a name that is structured like "<MYSPE>-Annotations.json" - e.g.
#     if MYSPE is called "Crptycoccus neoformans", your file should be called
#     "CRYNE-Annotations.json" (and the "name" of your Mbp1 orthologue is
#     "MBP1_CRYNE").
#
#   - Open the file in the RStudio editor and delete all blocks for
#     the Mbp1 protein annotations except the first one.
#
#   - From that block, delete all lines except for the line that says:
#
# {"pName" : "MBP1_SACCE", "fName" : "APSES fold", "start" : "4", "end" : "102"},
#
#   - Then delete the comma at the end of the line (your file will just have
#     this one annotation).
#
#   - Edit that annotation: change MBP1_SACCE  to MBP1_<MYSPE> and change the
#     "start" and "end" features to the coordinates you just discovered for the
#     APSES domain in your sequence.
#
#   - Save your file.
#
##   - Validate your file online at https://jsonlint.com/
#
#   - Update your "makeProteinDB.R" script to load your new
#     annotation when you recreate the database. Open the script in the
#     RStudio editor, and add the following command at the end:
#
#     myDB <- dbAddAnnotation(myDB, fromJSON("<MYSPE>-Annotations.json"))
#
#   - save and close the file.
#
# Then SKIP the next section.
#
#
# ===  4.1.2  If you HAVE done the BIN-FUNC-Annotation unit    
#
#
#   You DO already have a file called "<MYSPE>-Annotations.json" in the
#   ./data/ directory:
#
#   - Open the file in the RStudio editor.
#
#   - Below the last feature lines (but before the closing "]") add the
#     following feature line (without the "#")
#
# {"pName" : "MBP1_SACCE", "fName" : "APSES fold", "start" : "4", "end" : "102"}
#
#   - Edit that annotation: change MBP1_SACCE  to MBP1_<MYSPE> and change the
#     "start" and "end" features to the coordinates you just discovered for the
#     APSES domain in your sequence.
#
#   - Add a comma after the preceding feature line.
#
#   - Save your file.
#
#   - Validate your file online at https://jsonlint.com/
#
#
# ==   4.2  Execute and Validate  ==============================================
#
#   - source() your database creation script:
#
#     source("makeProteinDB.R")
#
#     This should run without errors or warnings. If it doesn't work and you
#     can't figure out quickly what's happening, ask on the mailing list for
#     help.
#
#   - Confirm
#     The following commands should retrieve the correct start and end
#     coordinates and sequence of the MBP1_MYSPE APSES domain:

sel <- which(myDB$protein$name == paste("MBP1_", biCode(MYSPE), sep = ""))

(proID <- myDB$protein$ID[sel])
(ftrID <- myDB$feature$ID[myDB$feature$name == "APSES fold"])
(fanID <- myDB$annotation$ID[myDB$annotation$proteinID == proID &
                             myDB$annotation$featureID == ftrID])
(start <- myDB$annotation$start[myDB$annotation$ID == fanID])
(end   <- myDB$annotation$end[myDB$annotation$ID == fanID])
(apses <- substr(myDB$protein$sequence[myDB$protein$ID == proID],
                 start,
                 end))


# [END]