bch441-work-abc-units/BIN-Storing_data.R

# BIN-Storing_data.R
#
# Purpose: A Bioinformatics Course:
#              R code accompanying the BIN-Storing_data unit
#
# Version: 1.1
#
# Date:    2017  10  08
# Author:  Boris Steipe (boris.steipe@utoronto.ca)
#
# V 1.1    Add instructions to retrieve UniProt ID from ID mapping service.
# V 1.0    First live version, complete rebuilt. Now using JSON data sources.
# V 0.1    First code copied from BCH441_A03_makeYFOlist.R
#
# TODO:
#
#
# == HOW TO WORK WITH LEARNING UNIT FILES ======================================
#
# DO NOT SIMPLY  source()  THESE FILES!
#
# 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        A Relational Datamodel in R: review                 62
#TOC>   1.1      Building a sample database structure               102
#TOC>   1.1.1    completing the database                            213
#TOC>   1.2      Querying the database                              248
#TOC>   1.3      Task: submit for credit (part 1/2)                 277
#TOC>   2        Implementing the protein datamodel                 289
#TOC>   2.1      JSON formatted source data                         315
#TOC>   2.2      "Sanitizing" sequence data                         355
#TOC>   2.3      Create a protein table for our data model          375
#TOC>   2.3.1    Initialize the database                            377
#TOC>   2.3.2    Add data                                           389
#TOC>   2.4      Complete the database                              409
#TOC>   2.4.1    Examples of navigating the database                436
#TOC>   2.5      Updating the database                              468
#TOC>   3        Add your own data                                  480
#TOC>   3.1      Find a protein                                     488
#TOC>   3.2      Put the information into JSON files                517
#TOC>   3.3      Create an R script to create your own database     540
#TOC>   3.3.1    Check and validate                                 560
#TOC>   3.4      Task: submit for credit (part 2/2)                 601
#TOC>
#TOC> ==========================================================================


# =    1  A Relational Datamodel in R: review  =================================

# A disclaimer at first: we are not building an industry-strength database at
# all here - but we are employing principles of such a database to keep common
# types of lab-data well organized. Don't think of this as emulating or even
# replacing a "real" database, but think of this as improving the ad hoc
# approaches we normally employ to store data in the lab. That does not mean
# such ad hoc approaches are necessarily bad - the best solution always depends
# on your objectives, the details of your tasks, and the context in which you
# are working.

# The principle we follow in implementing a relational data model is to build a
# list of dataframes . This list is our "database":
#  - Each _entity_ of the datamodel is a dataframe. In an SQL database, these
#    would also be called "tables". In a spreadsheet this would be a "sheet".
#  - Each instance of an entity, i.e. one stored _item_, is a row of the data
#    frame. In an SQL database this would be a record. In a spreadsheet this is
#    a row.
#  - Each _attribute_ of an entity is is a column of the dataframe. In an SQL
#    database this is a column, in a spreadsheet too.
#  - This doesn't necessarily solve the question of how we will store and curate
#    our source data - we will defer that to later. At first we talk only about
#    data representation internal to our R session, where we need it for
#    processing and analysis.

# Lets review syntax for creating and accessing such a structure, a list of data
# frames. You'll have to be absolutely confident with this, or you'll get lost
# in all the later learning units. We'll start from a compact example, a tiny
# database of philosophers to keep things brief. That database will have three
# tables: person, works and book. Person stores biographical data, book stores
# books, and works is a join table associating persons with their work. You
# should already be familiar with "join tables" and why we need them. This is
# the structure:
#
# person: id, name, born, died, school
# book:   id, title, published
# works:  id, person$id, book$id

# Perhaps draw out this schema to make things more clear.

# ==   1.1  Building a sample database structure  ==============================

# Let's build this structure.

philDB <- list()  # This is an empty list

# This is a data frame that we initialize with two philosophers
x <- data.frame(id = c(1,2),
                name = c("Laozi", "Martin Heidegger"),
                born = c(NA, "1889"),
                died = c("531 BCE", "1976"),
                school = c("Daoism", "Phenomenology"),
                stringsAsFactors = FALSE)
str(x)

# Lets add the dataframe to the philDB list and call it "person" there.
philDB[["person"]] <- x
str(philDB)

# and let's remove x so we don't mix up things later.
rm(x)

# We can address elements with the usual subsetting operators. I will use
# the $ operator for tables and columns, the [] operator for elements in
# columns. For example ...

philDB$person$name[1]   # Laozi

# task: Write an expression that returns all "school" entries from the
#       person table.


# Let's now add another person. There are several ways to do this, the
# conceptually cleanest is to create a one-row dataframe with the data, and
# rbind() it to the existing dataframe. Doing this, we must take care that
# the data frame column names are identical. What happens if they are not?
# Let's find out:
(x <- data.frame(a=1:4, b=11:14))
(y <- data.frame(a=6, c=17))
rbind(x, y)

(y <- data.frame(a=6, b=17))
rbind(x, y)

# All clear? That's good - this behaviour provides us with a sanity check on the
# operation.

(x <- data.frame(id = 2,
                 name = "Zhuangzi",
                 born = "369 BCE",
                 died = "286 BCE",
                 school = "Daoism"))

# Add this to the "person" table in our database with rbind() ...
philDB$person <- rbind(philDB$person, x)

# ... and examine the result:
str(philDB)

# Now one thing you should note is that we had forgotten to declare
# stringsAsFactors = FALSE when we created x - but this did not damage
# the database. This is because the existing columns had type chr and the
# implicit coercion, e.g. ...
as.character(x$name)
# happened to do the right thing. Don't rely on that. The Right Way is to
# turn factors off, even when you are making just a single row.

# But we made a serious error in our data! Did you spot it?
#
# If not, look at ...
philDB$person$id
# ... does that look oK?
#
# Absolutely not! id is the Primary Key in the table, and it has to be
# unique. How can we guarantee it to be unique? Certainly not when we
# enter it by hand. We need a function that generates a unique key. Here's
# a simple version, without any error-checking. It assumes that a column
# named "id" exists in the table, and that it holds the Primary Keys:
autoincrement <- function(table) {
  return(max(table$id) + 1)
}

#Try it:
autoincrement(philDB$person)

# Once that is clear, let's remove the Zhuangzi entry and recreate it correctly.
# Many ways to remove, here we use a logical expression to select matching record(s),
# apply the results to subset the data frame, and overwrite the existing table
# with the new one.

sel <- !(philDB$person$name == "Zhuangzi")   # select ...
philDB$person <- philDB$person[sel, ]        # ... and replace

str(philDB)

# Now let's add Zhuangzi with correct data. Note how we use the autoincrement
# function for the id

x <- data.frame(id = autoincrement(philDB$person),
                name = "Zhuangzi",
                born = "369 BCE",
                died = "286 BCE",
                school = "Daoism",
                stringsAsFactors = FALSE)
philDB$person <- rbind(philDB$person, x)
str(philDB)

# So far so good. Be honest with yourself. If you didn't follow any of this,
# go back, re-read, play with it, and ask for help. This is essential.


# ===  1.1.1  completing the database


# Next I'll add one more person, and create the other two tables:

x <- data.frame(id = autoincrement(philDB$person),
                name = "Kongzi",
                born = "551 BCE",
                died = "479 BCE",
                school = "Confucianism",
                stringsAsFactors = FALSE)
philDB$person <- rbind(philDB$person, x)


philDB[["books"]] <- data.frame(id = 1:5,
                                title = c("Zhuangzi",
                                          "Analects",
                                          "Being and Time",
                                          "Daodejing",
                                          "On the Way to Language"),
                                published = c("300 BCE",
                                              "220 BCE",
                                              "1927",
                                              "530 BCE",
                                              "1959"),
                                stringsAsFactors = FALSE)

philDB[["works"]] <- data.frame(id = 1:5,
                                personID = c(3, 4, 2, 1, 2),
                                bookID = c(1, 2, 3, 4, 5),
                                stringsAsFactors = FALSE)

str(philDB)


# ==   1.2  Querying the database  =============================================


# To retrieve data, we need to subset tables, possibly based on conditions we
# find in other tables. Sometimes we can simply get the information, e.g.
# all names ...
philDB$person$name
# ... or all book titles ...
philDB$books$title

# ... but sometimes we need to cross-reference information via join tables. Here
# is an example where we list authors and their works, sorted alphabetically by
# author:
(sel <- order(philDB$person$name))   # check out ?order and describe to
                                     # someone you know what it does, so that
                                     # you are sure you understand it.
(pID <- philDB$person$id[sel])
sel <- numeric()   # initialize the vector
for (ID in pID) {
  sel <- which(philDB$works$personID == ID)          # get all rows for which
                                                     #   the condition is TRUE
  cat(sprintf("%s: ", philDB$person$name[ID]))       # output the person
  cat(sprintf("\"%s\"  ", philDB$books$title[sel]))  # output the book
  cat("\n")
}

# Examine the intermediate results and trace the logic until this is clear.


# ==   1.3  Task: submit for credit (part 1/2)  ================================


#    Write and submit code that adds another philosopher to the datamodel:
#       Immanuel Kant, (1724 - 1804), Enlightenment Philosophy.
#       Works: Critique of Pure Reason (1781), Critique of Judgement (1790)
#    Write and submit code that lists the books in alphabetical order,
#    followed by the author and the year of publishing. Format your output like:
#    "Analects" - Kongzi (220 BCE)
#    Show the result.


# =    2  Implementing the protein datamodel  ==================================


# Working with the code above has probably illustrated a few concerns about
# curating data and storing it for analysis. In particular the join tables
# seem problematic - figuring out the correct IDs, it's easy to make
# mistakes.
#  - Data needs to be captured in a human-readable form so it can be verified
#      and validated;
#  - Some aspects of the database should _never_ be done by hand because they
#      errors are easy to make and hard to see. That essentially includes
#      every operation that has to do with abstract, primary keys;
#  - Elementary operations we need to support are: adding data, selecting
#      data, modifying data and deleting data.

# We will therefore construct our database in the following way:
#  - For each table, we will keep the primary information in JSON files. There
#      it is easy to read, edit if needed, and modify it.
#  - We will use simple scripts to read the JSON data and assemble it in
#      our database for further analysis.
#  - I have constructed initial files for yeast Mbp1 and nine other reference
#      species.
#  - I have written a small number of utility functions to read those files
#      and assemble them into a database.


# ==   2.1  JSON formatted source data  ========================================


# Have a look at the structure of the yeast Mbp1 protein data:
file.show("./data/MBP1_SACCE.json")

# - The whole thing is an array: [ ... ]. This is not necessary for a single
#     object, but we will have more objects in other files. And it's perfectly
#     legal to have an array with a single element.
# - The data is formatted as "key" : "value" pairs inside an object { ... }.
#     This keeps the association between data items and their semantics
#     explicit.
# - All keys are strings and they are unique in the object.
# - Values are mostly single strings and integers ...
# - ... except for "sequence". That one is an array of strings. Why? This is to
#     make it easier to format and maintain the data. JSON does not allow line
#     breaks within strings, but the strings we copy/paste from Genbank or other
#     sources might have line breaks, sequence numbers etc. So we need to
#     sanitize the sequence at some point. But since we need to do that
#     anyway, it is easier to see the whole sequence if we store it in chunks.

# Let's load the "jsonlite" package and have a look at how it reads this data.

if (! require(jsonlite, quietly=TRUE)) {
  install.packages("jsonlite")
  library(jsonlite)
}
# Package information:
#  library(help = jsonlite)       # basic information
#  browseVignettes("jsonlite")    # available vignettes
#  data(package = "jsonlite")     # available datasets


x <- fromJSON("./data/MBP1_SACCE.json")
str(x)

x$name
unlist(x$sequence)


# ==   2.2  "Sanitizing" sequence data  ========================================


# Examine the dbSanitizeSequence() function:

dbSanitizeSequence

# Try:

dbSanitizeSequence(c("GAA", "ttc"))
dbSanitizeSequence("MsnQ00%0 I@#>YSary    S
                     G1 V2DV3Y>")
x <- "        1 msnqiysary sgvdvyefih stgsimkrkk ddwvnathil kaanfakakr trilekevlk
       61 ethekvqggf gkyqgtwvpl niakqlaekf svydqlkplf dftqtdgsas pppapkhhha
      121 skvdrkkair sastsaimet krnnkkaeen qfqsskilgn ptaaprkrgr pvgstrgsrr
..." # copy/paste from Genbank

dbSanitizeSequence(x)


# ==   2.3  Create a protein table for our data model  =========================

# ===  2.3.1  Initialize the database


# The function dbInit contains all the code to return a list of empty
# data frames for our data model.

dbInit

myDB <- dbInit()
str(myDB)


# ===  2.3.2  Add data


# fromJSON() returns a dataframe that we can readily process to add data
# to our table. Have a look at the function to add protein entries:

dbAddProtein

myDB <- dbAddProtein(myDB, fromJSON("./data/MBP1_SACCE.json"))
str(myDB)

# Lets check that the 833 amino acids of the yeast MBP1 sequence have
# safely arrived. Note the genral idiom we use here to retrieve the data:
# we define a boolean vector that satisfies a condition, then we subset
# a column with that vector.

sel <- myDB$protein$name == "MBP1_SACCE"
nchar(myDB$protein$sequence[sel])


# ==   2.4  Complete the database  =============================================


# Completing the database with Mbp1 data and data for 9 other "reference"
# species is more of the same. I have assembled the code in a script
# "./scripts/ABC-createRefDB.R" - open it, check it out, and then source it.
# It's really very simple, just reading some prepared files of data I have
# formatted with JSON, and assembling the data in our data model.
#
# The code is also very simple and in particular there is no checking for errors
# or inconsistencies. Have a look:

# Totally straightforward ...
dbAddTaxonomy
dbAddFeature

# Just slightly more complex, since we need to match the protein or feature
# name in the JSON file with its internal ID, and, when doing that confirm
# that it CAN be matched and that the match is UNIQUE
dbAddAnnotation

# Now: create the database
source("./scripts/ABC-createRefDB.R")

str(myDB)


# ===  2.4.1  Examples of navigating the database


# You can look at the contents of the tables in the usual way we access
# elements from lists and dataframes. Here are some examples:
myDB$protein
myDB$protein$RefSeqID
myDB$protein[,"name"]
myDB$taxonomy
myDB$taxonomy$species
biCode(myDB$taxonomy$species)

# Comparing two tables:
# Are all of the taxonomyIDs in the protein table present in the
# taxonomy table? We ought to check, because the way we imported the
# data from JSON objects, we could have omitted or forgotten some. But we can
# check this with one simple expression. Unravel it and study its components.

all(myDB$protein$taxonomyID %in% myDB$taxonomy$ID)

# If this is not TRUE, you MUST fix the problem before continuing.

# Cross-referencing information:
# What is the species name of the protein whose name is "MBP1_COPCI"?

sel <- myDB$protein$name == "MBP1_COPCI"
x <- myDB$protein$taxonomyID[sel]

sel <- myDB$taxonomy$ID == x
myDB$taxonomy$species[sel]


# ==   2.5  Updating the database  =============================================


# Basic tasks for databases include retrieving data, selecting data, updating
# and deleting data. Here we will take a simple, pedestrian approach:
#
#     In case we need to modify any of the data, we modify it in the JSON file
#     save that, and recreate the database. The myDB database will only be
#     used for analysis.
#


# =    3  Add your own data  ===================================================


# You have chosen an organism as "MYSPE", and you final task will be to find the
# protein in MYSPE that is most similar to yeast Mbp1 and enter its information
# into the database.


# ==   3.1  Find a protein  ====================================================


# The BLAST algorithm will be properly introduced in a later learning unit -
# for now just use it in the following way:
#
# - Navigate to https://blast.ncbi.nlm.nih.gov/Blast.cgi and click on
#   Protein BLAST.
# - Enter NP_010227 into the "Query Sequence" field.
# - Choose "Reference proteins (refseq_protein)" as the "Database".
# - Paste the MYSPE species name into the "Organism" field.
#
# - Click "BLAST".

# You will probably get more than one result. If you get dozens of results or
# more, or if you get no results, something went wrong. Reconsider whether the
# problem was with your input, try something different, or ask for help.

# Otherwise, look for the top-hit in the "Alignments" section. In some cases
# there will be more than one hit with nearly similar E-values. If this is the
# case for MYSPE, choose the one with the higher degree of similarity (more
# identities) with the N-terminus of the query - i.e. the Query sequence of
# the first ~ 100 amino acids.

# -  Follow the link to the protein data page, linked from "Sequence ID".
# -  From there, in a separate tab, open the link to the taxonomy database page
#      for MYSPE which is linked from the "ORGANISM" record.


# ==   3.2  Put the information into JSON files  ===============================


# - Next make a copy of the file "./data/MBP1_SACCE.json" in your project
#     directory and give it a new name that corresponds to MYSPE - e.g. if
#     MYSPE is called "Crptycoccus neoformans", your file should be called
#     "MBP1_CRYNE.json"; in that case "MBP1_CRYNE" would also be the
#     "name" of your protein. Open the file in the RStudio editor and replace
#     all of the MBP1_SACCE data with the corresponding data of your protein.
#
#     The UniProt ID may not be discoverable from the NCBI page. To retrieve
#     it, navigate to http://www.uniprot.org/mapping/ , paste your RefSeq ID
#     into the query field, make sure "RefSeqProtein" is selected for "From"
#     and "UniProtKB" is selected for "To", and click "Go". In case this does
#     not retrieve a single UniProt ID, contact me.
#
# - Do a similar thing for the MYSPE taxonomy entry. Copy
#     "./data/refTaxonomy.json" and make a new file named "MYSPEtaxonomy.json".
#     Create a valid JSON file with only one single entry - that of MYSPE.
#
# - Validate your two files online at https://jsonlint.com/


# ==   3.3  Create an R script to create your own database  ====================


# Next: to create your own database.
# - Make a new R script, call it "makeProteinDB.R"
# - enter the following expression as the first command:
#     source("./scripts/ABC-createRefDB.R")
# - than add the two commands that add your protein and taxonomy data,
#     they should look like:
#     myDB <- dbAddProtein(    myDB, fromJSON("MBP1_<code>.json"))
#     myDB <- dbAddTaxonomy(   myDB, fromJSON("MYSPEtaxonomy.json"))
#
# - save the file and source() it:
#     source("makeProteinDB.R")

# This command needs to be executed whenever you recreate
# the database. In particular, whenver you have added or modified data
# in any of the JSON files. Later you will add more information ...


# ===  3.3.1  Check and validate


# Is your protein named according to the pattern "MBP1_MYSPE"? It should be.
# And does the taxonomy table contain the systematic name? It should be the same
# that you get when you type MYSPE into the console.

# Let's compute sequence lengths on the fly (with the function nchar() ), and
# open this with the table viewer function View()

View(cbind(myDB$protein[ , c("ID", "name", "RefSeqID")],
           length = nchar(myDB$protein$sequence)))

# Does your protein appear in the last row of this table? Where does your
# protein's length fall relative to the reference proteins? About the same? Much
# shorter? Much longer? If it is less then 500 amino acids long, I would suspect
# an error. Contact me for advice.

# Is that the right sequence? Is it the same as the one on the NCBI protein
# database page?
myDB$protein$sequence[nrow(myDB$protein)]

# If not, don't continue! Fix the problem first.
# Let me repeat: If this does not give you the right sequence of the MYSPE
#                Mbp1 homologue, DO NOT CONTINUE. Fix the problem.

# Is that the right taxonomy ID and binomial name for MYSPE?
sel <- myDB$taxonomy$species == MYSPE
myDB$taxonomy[sel, ]

# If not, or if the result was "<0 rows> ... " then DO NOT CONTINUE.
# Fix the problem first.

# Does this give you the right refseq ID for MBP1_MYSPE?
sel <- myDB$protein$name == paste0("MBP1_", biCode(MYSPE))
myDB$protein$RefSeqID[sel]

# If not, or if the result was "<0 rows> ... " then DO NOT CONTINUE.
# Fix the problem first.


# ==   3.4  Task: submit for credit (part 2/2)  ================================


# - On your submission page, note the E-value of your protein and link
#     to its NCBI protein database page.
# - Copy and paste the contents of your two JSON files on your submission
#     page on the Student Wiki
# - Execute the two commands below and show the result on your submission page

biCode(myDB$taxonomy$species) %in% biCode(MYSPE)
myDB$protein$taxonomyID %in% myDB$taxonomy$ID[(myDB$taxonomy$species == MYSPE)]

# That is all.


# [END]