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

321 lines
12 KiB
R
Raw Normal View History

# tocID <- "RPR-FASTA.R"
#
2017-10-14 22:25:46 +00:00
#
# Purpose: A Bioinformatics Course:
# R code accompanying the RPR-FASTA unit.
#
2020-11-13 05:38:08 +00:00
# Version: 1.1.1
2017-10-14 22:25:46 +00:00
#
# Date: 2017-10 - 2020-09
2017-10-14 22:25:46 +00:00
# Author: Boris Steipe (boris.steipe@utoronto.ca)
#
# Versions:
2020-11-13 05:38:08 +00:00
# 1.1.1 bugfix - wrong function name
# 1.1 2020 Maintenance. Rewrite validation logic. Add data
# to utilities. Define AACOLS
2017-10-14 22:25:46 +00:00
# 1.0 New unit.
#
#
# TODO: Make a simple solution first, then extend it to error checking, and
# to handle .mfa files.
2017-10-14 22:25:46 +00:00
#
#
# == 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 ...
#
# ==============================================================================
2017-10-14 22:25:46 +00:00
#TOC> ==========================================================================
#TOC>
#TOC> Section Title Line
#TOC> -----------------------------------------------------
2020-11-13 05:38:08 +00:00
#TOC> 1 Reading and validating FASTA 45
#TOC> 1.1 Validating FASTA 81
#TOC> 2 Parsing FASTA 227
#TOC> 3 Interpreting FASTA 247
#TOC> 4 Writing FASTA 274
#TOC>
2017-10-14 22:25:46 +00:00
#TOC> ==========================================================================
2017-10-14 22:25:46 +00:00
# = 1 Reading and validating FASTA ========================================
2017-10-14 22:25:46 +00:00
# FASTA is a text based format, structured in lines that are separated by
# line-feed or paragraph-break characters. Which one of these is used, depends
# on your operating system. But R's readLines() function knows how to handle
2017-10-14 22:25:46 +00:00
# these correctly, accross platforms. Don't try to read such files "by hand".
# Here is the yeast Mbp1 gene, via SGD.
file.show("./data/S288C_YDL056W_MBP1_coding.fsa")
faMBP1 <- readLines("./data/S288C_YDL056W_MBP1_coding.fsa")
2017-10-14 22:25:46 +00:00
# The warning is generated because the programmer at the NCBI who implemented
# the code to write this FASTA file neglected to place a line-break character
# after the last sequence character. While this is not technically incorrect,
# it is poor practice.
2017-10-14 22:25:46 +00:00
head(faMBP1)
2017-10-14 22:25:46 +00:00
# Note that there are NO line-break characters ("\n") at the end of these
# strings, even though they were present in the original file. readLines()
# has "consumed" these characters while reading - but every single line is in
# a vector of its own.
2017-10-14 22:25:46 +00:00
tail(faMBP1)
2017-10-14 22:25:46 +00:00
# Also note that the last line has fewer characters - this means readLines()
# imported the whole line, despite it not being terminated by "\n".
2017-10-14 22:25:46 +00:00
# It's very straightforward to work with such data, for example by collapsing
# everything except the first line into a single string ...
2017-10-14 22:25:46 +00:00
f <- c(faMBP1[1], paste(faMBP1[-1], sep = "", collapse = ""))
2017-10-14 22:25:46 +00:00
f[1]
nchar(f[2])
# == 1.1 Validating FASTA ==================================================
# The code above is making the assumption that everything from line 2 until
# the end IS sequence, the whole sequence and nothing but sequence.
# That assumption can break down in many ways:
2017-10-14 22:25:46 +00:00
#
# - there could be more than one header line. The specification says otherwise,
# but some older files use multiple, consecutive header lines. You don't
# want that to end up in your sequence.
# - this could be not a FASTA file at all. It could be raw sequence, a
# different sequence file format, or a wholly different file altogether.
# If you look at the file, you can immediately tell, but if you are
# reading the file in a complex workflow, your could easily import wrong
# data into your analysis.
# - there could be more than one sequence in the file. Such Multi-FASTA files
# occur commonly, as downloads of ORFs from genome regions or other
# sets of genes or proteins, or as the input / output for multiple
# sequence alignment programs.
#
# Data "from the wild" can (and usually does) have the most unexpected
# variations and it is really, really important to be clear about the
# assumptions that you are making. It is possible to "fix" things, according
# to the "Robustness Principle" :
# "Be conservative in what you send,
# be liberal in what you accept".
# (cf. https://en.wikipedia.org/wiki/Robustness_principle )
# ... but if you think about this, that's actually a really poor idea,
# which is much more likely to dilute standards, make unwarranted
# assumptions, and allow errors to pass silently and corrupt data.
2017-10-14 22:25:46 +00:00
#
# Let's discard this principle on the trash-heap of
# things-that-sound-like-a-good-idea-but-aren't. What we do instead is test,
# identify problems, and follow the principle: "crash early, crash often". Of
# course I can write code that would reformat any possible input as a FASTA
# file - but what good will it do me if it parses the file I receive
# from a server into FASTA format like:
2017-10-14 22:25:46 +00:00
#
# >404- Page Not Found</title</head>
# dyh-PagentfndhpThepageyreqesteddesnteistnthisserverCheckthe
# spellingrcntacttheadministratrsdyhtml
#
# Therefore, we write ourselves a FASTA checker that will enforce the following:
# (1) a FASTA file contains one or more sequences separated by zero or
# more empty lines
# (2) a sequence contains one header line followed by
# one or more sequence lines
# (3) a sequence line contains one or more uppercase or lowercase single
# letter amino acid codes, hyphens (gap character), or * (stop).
#
# Anything else should generate an error.
# (Case 1): Header(s) exist
fX <- c("ABC",
"defghi",
"klmnpq")
sel <- grepl("^>", fX) # "^>" is a regular expression that
# means: the exact character ">" at the
# beginning ("^") of the line.
if ( ! any(sel) ) { stop("no header lines in input.") }
# (Case 2) No adjacent header lines
fX <- c(">ABC",
">123",
"defghi",
"klmnpq")
sel <- grepl("^>", fX)
sel <- sel[- length(sel)] & sel[-1] # comparing shifted vectors
if ( any(sel)) { stop("adjacent header lines in input.") }
# (Case 3.1) all sequence lines contain only valid characters
# (constants for valid characters AAVALID, NUCVALID, and NUCAMBIG
# are defined with the .utilities.R script)
AAVALID
fX <- c(">ABC",
"def ;-) ghi",
"klmnpq")
myRegex <- sprintf("[^%s]", AAVALID) # NOT a valid character
sel <- ! grepl("^>", fX) # NOT headers
if (any(grepl(myRegex, fX[sel]))) {
stop("invalid chracter(s) outside of header lines.")
}
2017-10-14 22:25:46 +00:00
# (Case 3.2) all headers are followed directly by
# at least one letter of sequence
fX <- c(">ABC",
"",
">123",
"defghi",
"klmnpq")
sel <- grep("^>", fX) + 1 # indexes of headers + 1
myRegex <- sprintf("[%s]+", AAVALID) # at least one valid character
if (! all(grepl(myRegex, fX[sel]))) {
stop("a header has no adjacent sequence.")
2017-10-14 22:25:46 +00:00
}
# Ah, you might ask - couldn't we just have dropped all empty lines, and
# then caught this in Case 2? No - for two reasons: we would still miss headers
# at the end of file, and, we would have changed the line numbering - and
# ideally our "production" function will create information about where the
# error is to be found.
2017-10-14 22:25:46 +00:00
# Now combine this into a function ...
2017-10-14 22:25:46 +00:00
val <- function(fa) {
2017-10-14 22:25:46 +00:00
if ( ! any(grepl("^>", fa)) ) {
stop("no header lines in input.")
}
2017-10-14 22:25:46 +00:00
sel <- grepl("^>", fa)
if ( any(sel[- length(sel)] & sel[-1])) {
stop("adjacent header lines in input.")
}
2017-10-14 22:25:46 +00:00
sel <- ! grepl("^>", fa)
if ( any(grepl(sprintf("[^%s]", AAVALID), fa[sel]))) {
stop("invalid chracter(s) outside of header lines.")
}
2017-10-14 22:25:46 +00:00
sel <- grep("^>", fa) + 1
if (! all(grepl(sprintf("[%s]+", AAVALID), fa[sel]))) {
stop("a header has no adjacent sequence.")
2017-10-14 22:25:46 +00:00
}
return(invisible(NULL))
2017-10-14 22:25:46 +00:00
}
# Here is an example
FA <- c(">head1",
"acdef",
"ghi",
"",
">head2",
"kl",
">head3",
"mn",
"pqrs")
2020-11-13 05:38:08 +00:00
val(FA) # ... should not create an error
# A somewhat more elaborate validateFA() function was loaded with the
# ./utilities.R script. It needs a bit more bookkeeping, since NCBI multi-
2020-11-13 05:38:08 +00:00
# fasta files have space-characters in their spacer lines. Try it ...
validateFA(FA)
# = 2 Parsing FASTA =======================================================
# Once we have validated our assumptions about our input, it's quite
# painless to parse it. I have put this together as a function and the function
# gets loaded from ./.utilities.R
#
2017-10-14 22:25:46 +00:00
# Lets try this:
# - the first 3 elements of faMBP1:
readFASTA(faMBP1[1:3])
2017-10-14 22:25:46 +00:00
# - a multi FASTA file of aligned APSES domain sequences:
2017-10-14 22:25:46 +00:00
refAPSES <- readFASTA("./data/refAPSES.mfa")
2017-10-14 22:25:46 +00:00
# Subset the sequence with "P39678" in the header
refAPSES[grep("P39678", refAPSES$head) ,]
2017-10-14 22:25:46 +00:00
# = 3 Interpreting FASTA ==================================================
2017-10-14 22:25:46 +00:00
# FASTA files are straightforward to interpret - just one thing may be of note:
# when working with strings, we can use substr(<string>, <start>, <stop>) to
# extract substrings, but more often we expand the string into a vector of
# single characters with strsplit(<string>, ""). strsplit() returns a list,
# to accommodate that <string> could be a vector of many elements, therefore
2017-10-14 22:25:46 +00:00
# we usually unlist() the result if we use it only on a single string.
# Example: How many positive charged residues in "MBP1_SACCE"?
s <- unlist(strsplit(refAPSES$seq[grep("MBP1_SACCE", refAPSES$head)], ""))
s
2017-10-14 22:25:46 +00:00
sum(grepl("[HKR]", s)) # 20 (+) charged residues. grepl() returns TRUE and FALSE
# for the characters, sum() coerces to 1 and 0
# respectively, and that gives us the result.
100 * sum(grepl("[HKR]", s)) / length(s) # in percent: 20.2 %
# residue distribution
x <- factor(s, levels = names(AACOLS))
pie(table(x)[names(AACOLS)], col = AACOLS)
2017-10-14 22:25:46 +00:00
# = 4 Writing FASTA =======================================================
2017-10-14 22:25:46 +00:00
# Writing FASTA files is mostly just the reverse of reading, with one
2017-10-14 22:25:46 +00:00
# twist: we need to break the long sequence string into chunks of the desired
# width. The FASTA specification calls for a maximum of 120 characters per line,
# but writing out much less than that is common, since it allows to comfortably
2017-10-14 22:25:46 +00:00
# view lines on the console, or printing them on a sheet of paper (do we still
# do that actually?). How do we break a string into chunks? A combination of
# seq(<from>, <to>, <by>) with substring(<string>, <start>, <stop>) will work
# nicely. (Note that substring() is vectorized, whereas substr() is not!) As we
# loop through our FASTA object in memory, we can build the output by c()'ing
# blocks of header + sequence to each other. For VERY large objects this might
# be slow - in that case, we might want to precalculate the size of the output
# object. But that's more of a hypothetical consideration.
( s <- refAPSES$seq[2] )
2017-10-14 22:25:46 +00:00
nchar(s)
w <- 30 # width of chunk
(starts <- seq(1, nchar(s), by = w)) # starting index of chunk
(ends <- c((starts - 1)[-1], nchar(s))) # ending index of chunk
# Task: Is this safe? What happens if nchar(s) is shorter than w?
# What happens if nchar(s) is an exact multiple of w?
substring(s, starts, ends)
# confirm that the output contains the first and last residue, and both
# residues adjacent to the breaks
2017-10-14 22:25:46 +00:00
# As always, the function has been defined in ".utilities.R" for to use
# any time... type writeFASTA to examine it.
2017-10-14 22:25:46 +00:00
# Let's try this...
2017-10-14 22:25:46 +00:00
writeFASTA(refAPSES, width = 40)
2017-10-14 22:25:46 +00:00
# roundtrip for validation: write refAPSES with a different format,
# read it back in - the new dataframe must be identical
# to the original dataframe.
fname <- tempfile()
writeFASTA(refAPSES, fn = fname, width = 30)
identical(refAPSES, readFASTA(fname))
2017-10-14 22:25:46 +00:00
# ...works for me :-)
2017-10-14 22:25:46 +00:00
# [END]