bch441-work-abc-units/scripts/ABC-makeMYSPElist.R

421 lines
15 KiB
R

# tocID <- "scripts/ABC-makeMYSPElist.R"
#
# Purpose: Create a list of genome sequenced fungi with protein annotations and
# Mbp1 homologues.
#
# Version: 1.3
#
# Date: 2016 09 - 2020 09
# Author: Boris Steipe (boris.steipe@utoronto.ca)
#
# Versions
# 1.3 Rewrite to change datasource. NCBI has not been updated
# since 2012. Use ensembl fungi as initial source.
# 1.2 Change from require() to requireNamespace()
# 1.1.2 Moved BLAST.R to ./scripts directory
# 1.1 Update 2017
# 1.0 First code 2016
#
# TODO:
#
# type out workflow
#
# ==============================================================================
#
# DO NOT source() THIS FILE!
#
# This file is code I provide for your deeper understanding of a process and
# to provide you with useful sample code. It is not actually necessary for
# you to run this code, but I encourage you to read it carefully and discuss
# if there are parts you don't understand.
#
# Run the commands that interact with the NCBI servers only if you want to
# experiment specifically with the code and/or parameters. I have commented out
# those parts. If you only want to study the general workflow, just load()
# the respective intermediate results.
#
#TOC> ==========================================================================
#TOC>
#TOC> Section Title Line
#TOC> --------------------------------------------------------
#TOC> 1 The strategy 56
#TOC> 2 PACKAGES AND INITIALIZATIONS 68
#TOC> 3 ENSEMBL FUNGI 76
#TOC> 3.1 Import 79
#TOC> 4 BLAST SEARCH 156
#TOC> 4.1 find homologous proteins 162
#TOC> 4.2 Identify species in "hits" 193
#TOC> 5 MERGE ENSEMBL AND BLAST RESULTS 283
#TOC> 6 STUDENT NUMBERS 366
#TOC>
#TOC> ==========================================================================
# = 1 The strategy ========================================================
# This script will create a list of "MYSPE" species and save it in an R object
# MYSPEspecies that is stored in the data subdirectory of this project from
# where it can be loaded. The strategy is as follows: we download a list of
# annotated fungal genomes from ensembl.fungi. All these are genome-sequenced
# species that have been annotated.
# Next we perform a BLAST search, to identify fungal species that have
# genes that are homologous to yeast MBP1.
#
# ...
# = 2 PACKAGES AND INITIALIZATIONS ========================================
# httr provides interfaces to Webservers on the Internet
if (! requireNamespace("httr", quietly = TRUE)) {
install.packages("httr")
}
# = 3 ENSEMBL FUNGI =======================================================
# == 3.1 Import ============================================================
# Navigate to https://fungi.ensembl.org and click on the link to the full
# list of all species: https://fungi.ensembl.org/species.html
# On the page, click on the spreadsheet symbol top right and choose
# "download whole table". The file will be named "Species.csv", in your
# usual downloads folder. Move it to the data folder, and read it.
sDat <- read.csv("./data/Species.csv")
str(sDat)
# The most obvious way to partition these is according to Classification ...
# (poking around a bit in the UniProt taxonomy database shows that the
# classification used here is the taxonomic rank of "order").
# how many classifications do we have?
length(unique(sDat$Classification)) # 66
# To have a good set for the class, we should have about 100.
# Let's see for which of these we can find Mbp1 homologues.
# First, we'll keep only the colums for name, classification, and taxID, and
# drop the rest ...
sDat <- sDat[ , c("Name", "Classification", "Taxon.ID")]
colnames(sDat) <- c("name", "order", "taxID")
# Next, we make an extra column: genus - the first part of the binomial name.
# We'll use the gsub() function, and for that we need a "regular expression"
# that matches to all characters from the first blank to the end of the string:
myPatt <- "\\s.*$" # one whitespace (\\s) ...
# followed by any character (.) 0..n times (*) ...
# until the end of the string
# using gsub() we substitue all matching characters with the empty string "" -
# this deletes the matching characters
# Test this:
gsub(myPatt, "", "Genus") # one word: unchanged
gsub(myPatt, "", "gEnus species") # two words: return only first
gsub(myPatt, "", "geNus species strain 123") # many words: return only first
# apply this to the "name" column and add the result as a separate column
# called "genus"
sDat$genus <- gsub(myPatt, "", sDat$name)
# what do we get?
c(head(unique(sDat$genus)),
tail(unique(sDat$genus))) # inspect the first and last few. Note that there
# is a problem that we have to keep in mind.
# (Always inspect your results!)
# Drop all rows for which the genus contains special chracters -
# like "[Candida]"
sDat <- sDat[ ! grepl("[^a-zA-Z]", sDat$genus) , ]
length(table(sDat$genus)) # how many genus?
hist(table(sDat$genus), col = "#E9F4FF") # Distribution ...
# most genus have very few, but
# some have very many species.
sort(table(sDat$genus), decreasing = TRUE)[1:10] # Top ten...
# We should have at least one species from each taxonomic order, but we can
# add a few genus until we have about 100 validated species.
# Let's add a column for species, by changing our regular expression a bit,
# using ^ (start of string), \\S (NOT a whitespace),
# and + (one or more matches), capturing the match (...), and returning
# it as the substitution (\\1) ...
myPatt <- "^(\\S+\\s\\S+)\\s.*$"
sDat$species <- gsub(myPatt, "\\1", sDat$name)
# And we reorder the columns, just for aesthetics:
sDat <- sDat[ , c("name", "species", "genus", "order", "taxID")]
# Final check:
any(grepl("[^a-zA-Z -]", sDat$species)) # FALSE means no special characters
#
# Now we check which of these have Mbp1 homologues ...
# = 4 BLAST SEARCH ========================================================
# We run a BLAST search to find all proteins related to yeast Mbp1 in any
# fungus. With the results, we'll annotate our sDat table.
# == 4.1 find homologous proteins ==========================================
#
# Use BLAST to fetch proteins related to Mbp1 and identify the species that
# contain them.
# Scripting against NCBI APIs is not exactly enjoyable - there is usually a fair
# amount of error handling involved that is not supported by the API in a
# principled way but requires rather ad hoc solutions. The code I threw together
# to make a BLAST interface (demo-quality, not research-quality) is in the file
# ./scripts/BLAST.R Feel encouraged to study how this works. It's a pretty
# standard task of communicating with servers and parsing responses - everyday
# fare in the bioinformatics lab. Surprisingly, there seems to be no good BLAST
# parser in currently available packages.
#
# DON'T use this for BLAST searches unless you have read the NCBI policy
# for automated tasks. If you indicriminately pound on the NCBI's BLAST
# server, they will blacklist your IP-address. See:
# https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=DeveloperInfo
#
# Use BLAST() to find yeast Mbp1 homologues in other fungi in refseq
# BLASThits <- BLAST("NP_010227", # Yeast Mbp1 RefSeq ID
# db = "refseq_protein", # database to search in
# nHits = 3000, # 945 hits in 2020
# E = 0.01, #
# limits = "txid4751[ORGN]") # = fungi
# saveRDS(BLASThits, file="data/BLASThits.rds")
#
# NO NEED TO ACTUALLY RUN THIS:you can load the results from the data directory
#
BLASThits <- readRDS(file = "data/BLASThits.rds")
# == 4.2 Identify species in "hits" ========================================
# This is a very big list that can't be usefully analyzed manually. Here
# we are only interested in the species names that it contains.
# How many hits in the list?
length(BLASThits$hits)
# Let's look at a hit somewhere down the list
str(BLASThits$hit[[277]])
# A fair amount of parsing has gone into the BLAST.R code to prepare the results
# in a useful way. The species information is in the $species element of every
# hit.
# Run a loop to extract all the species names into a vector. We subset ...
# Blasthits$hits ... the list of hits, from which we choose ...
# Blasthits$hits[[i]] ... the i-th hit, and get ...
# Blasthits$hits[[i]]$species ... the species element from that.
# Subsetting FTW.
BLASTspecies <- character()
for (i in seq_along(BLASThits$hits)) {
BLASTspecies[i] <- BLASThits$hits[[i]]$species
}
# You can confirm that BLASTspecies has the expected size.
length(BLASTspecies)
# if we delete some of these later on, we still want to remember which hit
# they came from. Thus we name() the elements with their index, which is the
# same as the index of the hit in BLASThits
names(BLASTspecies) <- 1:length(BLASTspecies)
# let's plot the distribution of E-values
eVals <- numeric()
for (i in seq_along(BLASThits$hits)) {
eVals[i] <- BLASThits$hits[[i]]$E
}
range(eVals)
sum(eVals == 0)
# let's plot the log of all values > 0 to see how they are distributed
# plotting only one vectyor of numbers plots their index as x, and
# their value as y ...
plot(log(eVals[eVals > 0]), col = "#CC0000")
# This is very informative: I would suspect that the first ten or so are
# virtually identical to the yeast protein, then we have about 700 hits with
# decreasing similarity, and then about 200 more that may actually be false
# positives. Also - we plotted them by index, that means the table is SORTED:
# Lower E-values strictly come before higher E-values.
# Again, some species appear more than once, e.g. ...
sum(BLASTspecies == "Saccharomyces cerevisiae")
# ... corresponding to the five homologous gene sequences (paralogues) of yeast.
# Therefore we remove duplicates. Removing duplicates will leave the FIRST
# in a list alone, and only remove the SUBSEQUENT ones. Which means, from each
# species, we will retain only the protein that has the highest similarity
# to yeast Mbp1, not any of its more distant paralogues.
sel <- ! duplicated(BLASTspecies)
BLASTspecies <- BLASTspecies[sel]
length(BLASTspecies)
# i.e. we got rid of about two thirds of the hits.
tail(BLASTspecies) # see how the names are useful!
# again - there are some special characters ...
# what are they?
BLASTspecies[grep("[^a-zA-Z ]", BLASTspecies)]
# remove the brackets ...
BLASTspecies <- gsub("\\[|\\]", "", BLASTspecies)
# drop any new duplicates ...
BLASTspecies <- BLASTspecies[ ! duplicated(BLASTspecies)]
# check the number again:
length(BLASTspecies)
# Think a bit about this: what may be the biological reason to find that
# on average, in 300 fungi across the entire phylogenetic tree, we have
# three sequences that are homologous to yeast Mbp1?
# Let's look at the distribution of E-values in this selection (Subsetting FTW):
# we plot all values that are TRUE in the vector "sel" that we created above,
# AND greater than 0
plot(log(eVals[sel & eVals > 0]), col = "#00CC00")
# = 5 MERGE ENSEMBL AND BLAST RESULTS =====================================
# Next we add the blast result to our sDat dataframe. We'll store the index,
# the E-value, and the Query-bounds from which we can estimate which domains
# of Mbp1 are actually covered by the hit. (True orthologues MUST align with
# Mbp1's N-terminal APSES domain.)
#
# First we pull the hits we wanted from the BLASTspecies:
iHits <- as.numeric(names(BLASTspecies))
length(iHits) # one index for each TRUE in sel
# add columns to sDat
l <- nrow(sDat)
sDat$iHit <- numeric(l) # index of the hit in the BLAST results
sDat$eVal <- numeric(l) # E-value of the hit
sDat$lAli <- numeric(l) # length of the aligned region
# extract and merge
for (iHit in iHits) {
thisSp <- BLASThits$hits[[iHit]]$species
sel <- sDat$species == thisSp
sDat$iHit[sel] <- iHit
sDat$eVal[sel] <- BLASThits$hits[[iHit]]$E
sDat$lAli[sel] <- BLASThits$hits[[iHit]]$lengthAli
}
# Are all reference species accounted for?
selA <- sDat$iHit != 0 # all rows which matched to a BLAST hit
REFspecies %in% sDat$species[selA] # yes, all there
selB <- sDat$species %in% REFspecies # all rows which have one of REF species
sum(selA & selB) # How many rows?
# sDat of course includes all duplicates. Some may be multiply sequenced, some
# may be different strains. We'll use the same strategy as before and keep
# only the best hit: order the rows by E-value, then drop all rows which
# are duplicated.
# drop all rows without BLAST hits ...
sDat <- sDat[ ! (sDat$iHit == 0) , ]
# order sDat by E-value ...
sDat <- sDat[order(sDat$eVal, decreasing = FALSE) , ]
# drop all rows with duplicated species ...
sDat <- sDat[ ! duplicated(sDat$species) , ]
# Lets look at the E-values ...
plot(log(sDat$eVal[sDat$eVal > 0]), col = "#00CC00")
# and alignment lengths ...
plot(sDat$lAli, col = "#00DDAA")
# How many ...
length(unique(sDat$name))
length(unique(sDat$species))
length(unique(sDat$genus))
length(unique(sDat$order))
# To get the final dataset, we remove the reference species with their
# entire orders ...
REForders <- unique(sDat$order[sDat$species %in% REFspecies])
sel <- sDat$order %in% REForders
REFdat <- sDat[sel , ]
sDat <- sDat[ ! sel , ]
# REFdat should now contain only the REFspecies ...
( REFdat <- REFdat[REFdat$species %in% REFspecies , ] )
# ... but all of them
sum(REFspecies %in% REFdat$species)
# ... and we have enough left in sDat to prune sDat to unique genus ...
sDat <- sDat[ ! duplicated(sDat$genus) , ]
# saveRDS(sDat, file = "data/sDat.rds")
# saveRDS(REFdat, file = "data/REFdat.rds")
# = 6 STUDENT NUMBERS =====================================================
#
# An asymmetric function to retrieve a MYSPE species
students <- read.csv("../BCH441-2020-students.csv")
sN <- students$Student.Number
range(sN)
any(duplicated(gsub(".+(.......)$", "\\1", sN)))
N <- 7
x <- numeric(N)
for (i in 1:N) {
x[i] <- H(substr(gsub(".+(.......)$", "\\1", sN), i, i))
}
plot(x, col = "#BB0000", type = "b")
keys <- as.numeric(gsub(".+(....).$", "\\1", sN))
any(duplicated(keys))
# =====
set.seed(112358)
names(sN) <- sample(1:nrow(sDat), length(sN))
MYSPEmap <- data.frame(keys = sprintf("%04d", 0:9999),
iMYSPE = sample(1:nrow(sDat), 10000, replace = TRUE))
rownames(MYSPEmap) <- MYSPEmap$keys
for (i in 1:length(sN)) {
rMap <- gsub(".+(....).$", "\\1", sN[i])
MYSPEmap[rMap, "iMYSPE"] <- as.integer(names(sN)[i])
}
# saveRDS(MYSPEmap, "./data/MYSPEmap.rds")
getMYSPE <- function(x) {
dat <- readRDS("./data/sDat.rds")
map <- readRDS("./data/MYSPEmap.rds")
key <- gsub(".+(....).$", "\\1", x)
return(dat$species[map[key, "iMYSPE"]])
}
# === validate
l <- length(sN)
sp <- character(l)
for(i in 1:l) {
sp[i] <- getMYSPE(sN[i])
}
any(duplicated(sp))
length(unique(sp))
which(! sDat$species %in% sp) # these can be assigned to late-comers
# Done.
# [END]