Difference between revisions of "BIO Assignment Week 6"

From "A B C"
Jump to navigation Jump to search
m
m
Line 2: Line 2:
 
<div class="b1">
 
<div class="b1">
 
Assignment for Week 6<br />
 
Assignment for Week 6<br />
<span style="font-size: 70%">Sensitive database searches with PSI-BLAST</span>
+
<span style="font-size: 70%">Function</span>
 
</div>
 
</div>
 
<table style="width:100%;"><tr>
 
<table style="width:100%;"><tr>
Line 15: Line 15:
  
 
__TOC__
 
__TOC__
 +
 +
  
  
Line 20: Line 22:
 
==Introduction==
 
==Introduction==
  
<div style="padding: 2px; background: #F0F1F7border:solid 1px #AAAAAA; font-size:125%;color:#444444">
+
&nbsp;
 +
 
 +
In this assignment we will perform a few computations with coordinate data in PDB files, look more in depth at domain annotations, and compare them across different proteins related to yeast Mbp1. We will write a function in '''R''' to help us plot a graphic of the comparison and collaborate to share data for our plots. But first we need to go through a few more '''R''' concepts.
  
  
&nbsp;<br>
+
&nbsp;
  
;Take care of things, and they will take care of you.
+
== Programming '''R''' code  ==
:''Shunryu Suzuki''
 
</div>
 
  
 +
First, we will cover essentials of '''R''' programming: the fundamental statements that are needed to write programs&ndash;conditional expressions and loops, and how to define ''functions'' that allow us to use programming code. But let's start with a few more data types of '''R''' so we can use the concepts later on: matrices, lists and data frames.
  
Anyone can click buttons on a Web page, but to use the powerful sequence database search tools ''right'' often takes considerable more care, caution and consideration.
 
  
Much of what we know about a protein's physiological function is based on the '''conservation''' of that function as the species evolves. We assess conservation by comparing sequences between related proteins. Conservation - or its opposite: ''variation'' - is a consequence of '''selection under constraints''': protein sequences change as a consequence of DNA mutations, this changes the protein's structure, this in turn changes functions and that has the multiple effects on a species' fitness function. Detrimental variants may be removed. Variation that is tolerated is largely neutral and therefore found only in positions that are neither structurally nor functionally critical. Conservation patterns can thus provide evidence for many different questions: structural conservation among proteins with similar 3D-structures, functional conservation among homologues with comparable roles, or amino acid propensities as predictors for protein engineering and design tasks.
+
{{task|
 +
Please begin by working through the short [http://biochemistry.utoronto.ca/steipe/abc/index.php/R_tutorial#Matrices '''R''' - tutorial: matrices] section and the following sections on "Lists" and "Data frames".
  
Measuring conservation requires alignment. Therefore a carefully done multiple sequence alignment ('''MSA''') is a cornerstone for the annotation of the essential properties a gene or protein. MSAs are also useful to resolve ambiguities in the precise placement of indels and to ensure that columns in alignments actually contain amino acids that evolve in a similar context. MSAs serve as input for
+
Note that what we have done here is just the bare minimum on vectors, matrices and lists. The concepts are very generally useful, and there are many useful ways to extract subsets of values. We'll come across these in various parts of '''R''' sample code. But unless you read the provided code examples line by line, make sure you understand every single statement and '''ask''' if you are not clear about the syntax, this will all be quickly forgotten. Use it or lose it!
* functional annotation;
+
}}
* protein homology modeling;
 
* phylogenetic analyses, and
 
* sensitive homology searches in databases.
 
  
In order to perform a multiple sequence alignment, we obviously need a set of homologous sequences. This is where the trouble begins. All interpretation of MSA results depends '''absolutely''' on how the input sequences were chosen. Should we include only orthologs, or paralogs as well? Should we include only species with fully sequenced genomes, or can we tolerate that some orthologous genes are possibly missing for a species? Should we include all sequences we can lay our hands on, or should we restrict the selection to a manageable number of ''representative'' sequences? All of these choices influence our interpretation:
 
*orthologs are expected to be functionally and structurally conserved;
 
*paralogs may have divergent function but have similar structure;
 
*missing genes may make paralogs look like orthologs; and
 
*selection bias may weight our results toward sequences that are over-represented and do not provide a fair representation of evolutionary divergence.
 
  
 +
'''R''' is a full featured programming language and in order to be that, it needs the ability to manipulate '''Control Flow''', i.e. the order in which instructions are executed. By far the most important control flow statement is a '''conditional branch''' i.e. the instruction to execute code at alternative points of the program, depending on the presence or absence of a condition. Using such conditional branch instructions, two main types of programming constructs can be realized: ''conditional expressions'' and ''loops''.
  
In this assignment, we will set ourselves the task to use PSI-BLAST and '''find all orthologs and paralogs of the APSES domain containing transcription factors in YFO'''. We will use these sequences later for multiple alignments, calculation of conservation ''etc''. The methodical problem we will address is: how do we perform a sensitive PSI-BLAST search '''in one organism'''. There is an issue to consider:
+
=== Conditional expressions ===
* If we restrict the PSI-BLAST search to YFO, PSI-BLAST has little chance of building a meaningful profile - the number of homologues that actually are '''in''' YFO is too small. Thus the search will not become very sensitive.
 
* If we don't restrict our search, but search in all species, the number of hits may become too large. It becomes increasingly difficult to closely check all hits as to whether they have good coverage, and how will we evaluate the fringe cases of marginal E-value, where we need to decide whether to include a new sequence in the profile, or whether to hold off on it for one or two iterations, to see whether the E-value drops significantly. Profile corruption would make the search useless. This is maybe still manageable if we restrict our search to fungi, but imagine you are working with a bacterial protein, or a protein that is conserved across the entire tree of life: your search will find thousands of sequences. And by next year, thousands more will have been added.
 
  
Therefore we have to find a middle ground: add enough species (sequences) to compile a sensitive profile, but not so many that we can no longer individually assess the sequences that contribute to the profile.
+
The template for conditional expression in R is:
  
 +
<source lang="rsplus">
 +
if( <expression 1> ) {
 +
  <statement 1>
 +
}
 +
else if ( <expresssion 2> ) {
 +
  <statement 2>
 +
}
 +
else {
 +
  <statement 3>
 +
}
  
Thus in practice, a sensitive PSI-BLAST search needs to address two issues before we begin:
+
</source>
# We need to define the sequence we are searching with; and
 
# We need to define the dataset we are searching in.
 
  
 +
...where both the <code>else if (...) { ... }</code> and the <code>else (...)</code> block are optional. We have encountered this construct previously, when we assigned the appropriate colors for amino acids in the frequency plot:
  
 +
<source lang="rsplus">
 +
if      (names(logRatio[i]) == "F") { barColors[i] <- hydrophobic }
 +
else if (names(logRatio[i]) == "G") { barColors[i] <- plain }
 +
  [... etc ...]
 +
else                                { barColors[i] <- plain }
  
 +
</source>
 +
 +
==== Logical expressions ====
  
 +
We have to consider the <code>&lt;expression&gt;</code> in a bit more detail: anything that is, or produces, or can be interpreted as, a Boolean TRUE or FALSE value can serve as an expression in a conditional statement.
  
==Defining the sequence to search with==
+
{{task|1=
 +
Here are some examples. Copy the code to an '''R''' script, predict what will happen on in each line and try it out:
  
 +
<source lang="rsplus">
 +
# A boolean constant is interpreted as is:
 +
if (TRUE)    {print("true")} else {print("false")}
 +
if (FALSE)  {print("true")} else {print("false")}
  
Consider again the task we set out from: '''find all orthologs and paralogs of the APSES domain containing transcription factors in YFO'''.  
+
# Strings of "true" and "false" are coerced to their
 +
# Boolean equivalent, but contrary to other programming languages
 +
# arbitrary, non-empty or empty strings are not interpreted.
 +
if ("true") {print("true")} else {print("false")}
 +
if ("false") {print("true")} else {print("false")}
 +
if ("widdershins")    {print("true")} else {print("false")}
 +
if ("") {print("true")} else {print("false")}
  
 +
# All non-zero, defined numbers are TRUE
 +
if (1)      {print("true")} else {print("false")}
 +
if (0)      {print("true")} else {print("false")}
 +
if (-1)    {print("true")} else {print("false")}
 +
if (pi)    {print("true")} else {print("false")}
 +
if (NULL)  {print("true")} else {print("false")}
 +
if (NA)    {print("true")} else {print("false")}
 +
if (NaN)    {print("true")} else {print("false")}
 +
if (Inf)    {print("true")} else {print("false")}
  
{{task|1=
+
# functions can return Boolean values
What query sequence should you use? Should you ...
+
affirm <- function() { return(TRUE) }
 +
deny <- function() { return(FALSE) }
 +
if (affirm())    {print("true")} else {print("false")}
 +
if (deny())    {print("true")} else {print("false")}
  
 +
# N.B. coercion of Booleans into numbers can be done as well
 +
and is sometimes useful: consider ...
 +
a <- c(TRUE, TRUE, FALSE, TRUE, FALSE)
 +
a
 +
as.numeric(a)
 +
sum(a)
  
# Search with the full-length Mbp1 sequence from ''Saccharomyces cerevisiae''?
+
# ... or coercing the other way ...
# Search with the full-length Mbp1 homolog that you found in YFO?
+
as.logical(-1:1)
# Search with the structurally defined ''S. cerevisiae'' APSES domain sequence?
+
</source>
# Search with the APSES domain sequence from the YFO homolog, that you have defined by sequence alignment with the yeast protein?
+
}}
# Search with the KilA-N domain sequence?
 
  
 +
==== Logical operators ====
  
<div class="mw-collapsible mw-collapsed" data-expandtext="Expand" data-collapsetext="Collapse" style="border:#000000 solid 1px; padding: 10px; margin-left:25px; margin-right:25px;">Reflect on this (pretend this is a quiz question) and come up with a reasoned answer. Then click on "Expand" to read my opinion on this question.
+
To actually write a conditional statement, we have to be able to '''test a condition''' and this is what logical operators do. Is something '''equal''' to something else? Is it less? Does something exist? Is it a number?
<div  class="mw-collapsible-content">
 
;The full-length Mbp1 sequence from ''Saccharomyces cerevisiae''
 
:Since this sequence contains multiple domains (in particular the ubiquitous Ankyrin domains) it is not suitable for BLAST database searches. You must restrict your search to the domain of greatest interest for your question. That would be the APSES domain.
 
  
;The full-length Mbp1 homolog that you found in YFO
+
{{task|1=
:What organism the search sequence comes from does not make a difference. Since you aim to find '''all''' homologs in YFO, it is not necessary to have your search sequence more or less similar to '''any particular''' homologs. In fact '''any''' APSES sequence should give you the same result, since they are '''all''' homologous. But the full-length sequence in YFO has the same problem as the ''Saccharomyces'' sequence.
 
  
;The structurally defined ''S. cerevisiae'' APSES domain sequence?
+
Here are some examples. Again, predict what will happen ...
:That would be my first choice, just because it is structurally well defined as a complete domain, and the sequence is easy to obtain from the <code>1BM8</code> PDB entry. (<code>1MB1</code> would also work, but you would need to edit out the penta-Histidine tag at the C-terminus that was engineered into the sequence to help purify the recombinantly expressed protein.)
 
  
;The APSES domain sequence from the YFO homolog, that you have defined by sequence alignment with the yeast protein?
+
<source lang="rsplus">
:As argued above: since they are all homologs, any of them should lead to the same set of results.
+
TRUE            # Just a statement.
  
;The KilA-N domain sequence?
+
#  unary operator
:This is a shorter sequence and a more distant homolog to the domain we are interested in. It would not be my first choice: the fact that it is more distantly related might make the search '''more sensitive'''. The fact that it is shorter might make the search '''less specific'''. The effect of this tradeoff would need to be compared and considered. By the way: the same holds for the even shorter subdomain 50-74 we discussed in the last assignment. However: one of the results of our analysis will be '''whether APSES domains in fungi all have the same length as the Mbp1 domain, or whether some are indeed much shorter, as suggested by the Pfam alignment.'''
+
! TRUE          # NOT ...
  
 +
# binary operators
 +
FALSE > TRUE    # GREATER THAN ...
 +
FALSE < TRUE    # ... these are coerced to numbers
 +
FALSE < -1       
 +
0 == FALSE      # Careful! == compares, = assigns!!!
  
So in my opinion, you should search with the yeast Mbp1 APSES domain, i.e. the sequence which you have previously studied in the crystal structure. Where is that? Well, you might have saved it in your journal, or you can get it again from the [http://www.pdb.org/pdb/explore/explore.do?structureId=1BM8 '''PDB'''] (i.e. [http://www.pdb.org/pdb/files/fasta.txt?structureIdList=1BM8 here], or from [[BIO_Assignment_Week_3#Search input|Assignment 3]].
+
"x" == "u"      # using lexical sort order ...
 +
"x" >= "u"
 +
"x" <= "u"
 +
"x" != "u"
 +
"aa" > "u"      # ... not just length, if different.
 +
"abc" < "u" 
  
</div>
+
TRUE | FALSE    # OR: TRUE if either is true
</div>
+
TRUE & FALSE    # AND: TRUE if both are TRUE
}}
 
  
&nbsp;
+
# equality and identity
 +
?identical
 +
a <- c(TRUE)
 +
b <- c(TRUE)
 +
a; b
 +
a == b
 +
identical(a, b)
  
==Selecting species for a PSI-BLAST search==
+
b <- 1
 +
a; b
 +
a == b
 +
identical(a, b)  # Aha: equal, but not identical
  
  
As discussed in the introduction, in order to use our sequence set for studying structural and functional features and conservation patterns of our APSES domain proteins, we should start with a well selected dataset of APSES domain containing homologs in YFO. Since these may be quite divergent, we can't rely on '''BLAST''' to find all of them, we need to use the much more sensitive search of '''PSI-BLAST''' instead. But even though you are interested only in YFO's genes, it would be a mistake to restrict the PSI-BLAST search to YFO. PSI-BLAST becomes more sensitive if the profile represents more diverged homologs. Therefore we should always search with a broadly representative set of species, even if we are interested only in the results for one of the species. This is important. Please reflect on this for a bit and make sure you understand the rationale why we include sequences in the search that we are not actually interested in.
+
# some other useful tests for conditional expressions
 +
?all
 +
?any
 +
?duplicated
 +
?exists
 +
?is.character
 +
?is.factor
 +
?is.integer
 +
?is.null
 +
?is.numeric
 +
?is.unsorted
 +
?is.vector
 +
</source>
 +
}}
  
  
But you can also search with '''too many''' species: if the number of species is large and PSI-BLAST finds a large number of results:
+
=== Loops ===
# it becomes unwieldy to check the newly included sequences at each iteration, inclusion of false-positive hits may result, profile corruption and loss of specificity. The search will fail.
 
# since genomes from some parts of the Tree Of Life are over represented, the inclusion of all sequences leads to selection bias and loss of sensitivity.
 
  
 +
Loops allow you to repeat tasks many times over. The template is:
  
We should therefore try to find a subset of species
+
<source lang="rsplus">
# that represent as large a '''range''' as possible on the evolutionary tree;
+
for (<name> in <vector>) {
# that are as well '''distributed''' as possible on the tree; and
+
  <statement>
# whose '''genomes''' are fully sequenced.
+
}
 +
</source>
  
These criteria are important. Again, reflect on them and understand their justification. Choosing your species well for a PSI-BLAST search can be crucial to obtain results that are robust and meaningful.
+
{{task|1=
 +
Consider the following: Again, copy the code to a script, study it, predict what will happen and then run it.
  
How can we '''define''' a list of such species, and how can we '''use''' the list?
+
<source lang="rsplus">
 +
# simple for loop
 +
for (i in 1:10) {
 +
print(c(i, i^2, i^3))
 +
}
  
The definition is a rather typical bioinformatics task for integrating datasources: "retrieve a list of representative fungi with fully sequenced genomes". Unfortunately, to do this in a principled way requires tools that you can't (yet) program: we would need to use a list of genome sequenced fungi, estimate their evolutionary distance and select a well-distributed sample. Regrettably you can't combine such information easily with the resources available from the NCBI.
+
# Compare excution times: one million square roots from a vector ...
 +
n <- 1000000
 +
x <- 1:n
 +
y <- sqrt(x)
  
We will use an approach that is conceptually similar: selecting a set of species according to their shared taxonomic rank in the tree of life. {{WP|Biological classification|'''Biological classification'''}} provides a hierarchical system that describes evolutionary relatedness for all living entities. The levels of this hierarchy are so called {{WP|Taxonomic rank|'''taxonomic ranks'''}}. These ranks are defined in ''Codes of Nomenclature'' that are curated by the self-governed international associations of scientists working in the field. The number of ranks is not specified: there is a general consensus on seven principal ranks (see below, in bold) but many subcategories exist and may be newly introduced. It is desired&ndash;but not mandated&ndash;that ranks represent ''clades'' (a group of related species, or a "branch" of a phylogeny), and it is desired&ndash;but not madated&ndash;that the rank is sharply defined. The system is based on subjective dissimilarity. Needless to say that it is in flux.
+
# ... or done explicitly in a for-loop
 +
for (i in 1:n) {
 +
  y[i] <- sqrt (x[i])
 +
}
  
If we follow a link to an entry in the NCBI's Taxonomy database, eg. [http://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=559292 ''Saccharomyces cerevisiae S228c''], the strain from which the original "yeast genome" was sequenced in the late 1990s, we see the following specification of its taxonomic lineage:
+
</source>
  
 +
''If'' you can achieve your result with an '''R''' vector expression, it will be faster than using a loop. But sometimes you need to do things explicitly, for example if you need to access intermediate results.
  
<source lang="text">
+
}}
cellular organisms; Eukaryota; Opisthokonta; Fungi; Dikarya;
 
Ascomycota; Saccharomyceta; Saccharomycotina; Saccharomycetes;
 
Saccharomycetales; Saccharomycetaceae; Saccharomyces; Saccharomyces cerevisiae
 
</source>
 
  
  
These names can be mapped into taxonomic ranks ranks, since the suffixes of these names e.g. ''-mycotina'', ''-mycetaceae'' are specific to defined ranks. (NCBI does not provide this mapping, but {{WP|Taxonomic rank|Wikipedia}} is helpful here.)
+
Here is an example to play with loops: a password generator. Passwords are a '''pain'''. We need them everywhere, they are frustrating to type, to remember and since the cracking programs are getting smarter they become more and more likely to be broken. Here is a simple password generator that creates random strings with consonant/vowel alterations. These are melodic and easy to memorize, but actually as '''strong''' as an 8-character, fully random password that uses all characters of the keyboard such as <code>)He.{2jJ</code> or <code>#h$bB2X^</code> (which is pretty much unmemorizable). The former is taken from 20<sup>7</sup> * 7<sup>7</sup> 10<sup>15</sup> possibilities, the latter is from 94<sup>8</sup> ~ 6*10<sup>15</sup> possibilities. HIgh-end GPU supported {{WP|Password cracking|password crackers}} can test about 10<sup>9</sup> passwords a second, the passwords generated by this little algorithm would thus take on the order of 10<sup>6</sup> seconds or eleven days to crack<ref>That's assuming the worst case in that the attacker needs to know the pattern with which the password is formed, i.e. the number of characters and the alphabet that we chose from. But note that there is an even worse case: if the attacker had access to our code and the seed to our random number generator. When the random number generator starts up, a new seed is generated from system time, thus the possible space of seeds can be devastatingly small. But even if a seed is set explicitly with the <code>set.seed()</code> function, that seed is a 32-bit integer and thus can take only a bit more than 4*10<sup>9</sup> values, six orders of magnitude less than the 10<sup>15</sup> password complexity we thought we had. It turns out that the code may be a much greater vulnerability than the password itself. Keep that in mind. <small>Keep it secret. <small>Keep it safe.</small></small></ref>. This is probably good enough to deter a casual attack.
  
<table>
+
{{task|1=
 +
Copy, study and run ...
 +
<source lang="rsplus">
 +
# Suggest memorizable passwords
 +
# Below we use the functions:
 +
?nchar
 +
?sample
 +
?substr
 +
?paste
 +
?print
  
<tr class="sh">
+
#define a string of  consonants ...
<td>Rank</td>
+
con <- "bcdfghjklmnpqrstvwxz"
<td>Suffix</td>
+
# ... and a string of of vowels
<td>Example</td>
+
vow <- "aeiouy"
</tr>
 
  
<tr class="s1">
+
for (i in 1:10) {  # ten sample passwords to choose from ...
<td>Domain</td>
+
    pass = rep("", 14)  # make an empty character vector
<td></td>
+
    for (j in 1:7) {    # seven consonant/vowel pairs to be created ...
<td>Eukaryota (Eukarya)</td>
+
        k  <- sample(1:nchar(con), 1)  # pick a random index for consonants ...
</tr>
+
        ch  <- substr(con,k,k)          #  ... get the corresponding character ...
 +
        idx <- (2*j)-1                  # ... compute the position (index) of where to put the consonant ...
 +
        pass[idx] <- ch                # ...  and put it in the right spot
  
<tr class="s2">
+
        # same thing for the vowel, but coded with fewer intermediate assignments
<td>&nbsp;&nbsp;Subdomain</td>
+
        # of results to variables
<td>&nbsp;</td>
+
        k <- sample(1:nchar(vow), 1)
<td>Opisthokonta</td>
+
        pass[(2*j)] <- substr(vow,k,k)
</tr>
+
    }
 +
    print( paste(pass, collapse="") )  # collapse the vector in to a string and print
 +
}
 +
</source>
 +
}}
  
<tr class="s1">
+
=== Functions ===
<td>'''Kingdom'''</td>
 
<td>&nbsp;</td>
 
<td>Fungi</td>
 
</tr>
 
  
<tr class="s2">
+
Finally: functions. Functions look very much like the statements we have seen above. the template looks like:
<td>&nbsp;&nbsp;Subkingdom</td>
 
<td>&nbsp;</td>
 
<td>Dikarya</td>
 
</tr>
 
  
<tr class="s1">
+
<source lang="rsplus">
<td>'''Phylum'''</td>
+
<name> <- function (<parameters>) {
<td>&nbsp;</td>
+
  <statements>
<td>Ascomycota</td>
+
}
</tr>
+
</source>
  
<tr class="s2">
+
In this statement, the function is assigned to the ''name'' - any valid name in '''R'''. Once it is assigned, it the function can be invoked with <code>name()</code>. The parameter list (the values we write into the parentheses followin the function name) can be empty, or hold a list of variable names. If variable names are present, you need to enter the corresponding parameters when you execute the function. These assigned variables are available inside the function, and can be used for computations. This is called "passing the variable into the function".
<td>&nbsp;&nbsp;''rankless taxon''<ref>The -myceta are well supported groups above the Class rank. See {{WP|Leotiomyceta|''Leotiomyceta''}} for details and references.</ref></td>
 
<td>-myceta</td>
 
<td>Saccharomyceta</td>
 
</tr>
 
  
<tr class="s1">
+
You have encountered a function to choose YFO names. In this function, your Student ID was the parameter. Here is another example to play with: a function that calculates how old you are. In days. This is neat - you can celebrate your 10,000 birth'''day''' - or so.
<td>&nbsp;&nbsp;Subphylum</td>
 
<td>-mycotina</td>
 
<td>Saccharomycotina</td>
 
</tr>
 
  
<tr class="s2">
+
{{task|1=
<td>'''Class'''</td>
 
<td>-mycetes</td>
 
<td>Saccharomycetes</td>
 
</tr>
 
  
<tr class="s1">
+
Copy, explore and run ...
<td>&nbsp;&nbsp;Subclass</td>
 
<td>-mycetidae</td>
 
<td>&nbsp;</td>
 
</tr>
 
  
<tr class="s2">
+
;Define the function ...
<td>'''Order'''</td>
+
<source lang = "rsplus">
<td>-ales</td>
+
# A lifedays calculator function
<td>Saccharomycetales</td>
 
</tr>
 
  
<tr class="s1">
+
myLifeDays <- function(date = NULL) { # give "date" a default value so we can test whether it has been set
<td>'''Family'''</td>
+
    if (is.null(date)) {
<td>-aceae</td>
+
        print ("Enter your birthday as a string in \"YYYY-MM-DD\" format.")
<td>Saccharomycetaceae</td>
+
        return()
</tr>
+
    }
 +
    x <- strptime(date, "%Y-%m-%d") # convert string to time
 +
    y <- format(Sys.time(), "%Y-%m-%d") # convert "now" to time
 +
    diff <- round(as.numeric(difftime(y, x, unit="days")))
 +
    print(paste("This date was ", diff, " days ago."))
 +
}
 +
</source>
  
<tr class="s2">
+
;Use the function (example):
<td>&nbsp;&nbsp;Subfamily</td>
+
<source lang = "rsplus">
<td>-oideae</td>
+
  myLifeDays("1932-09-25")  # Glenn Gould's birthday
<td>&nbsp;</td>
+
</source>
</tr>
+
}}
  
<tr class="s1">
+
Here is a good opportunity to play and practice programming: modify this function to accept a second argument. When a second argument is present (e.g. 10000) the function should print the calendar date on which the input date will be that number of days ago. Then you could use it to know when to celebrate your 10,000<sup>th</sup> lifeDay, or your 777<sup>th</sup> anniversary day or whatever.
<td>&nbsp;&nbsp;Tribe</td>
 
<td>-eae</td>
 
<td>&nbsp;</td>
 
</tr>
 
  
<tr class="s2">
+
Enjoy.
<td>&nbsp;&nbsp;Subtribe</td>
 
<td>-ineae</td>
 
<td>&nbsp;</td>
 
</tr>
 
  
<tr class="s1">
+
== Protein structure annotation ==
<td>'''Genus'''</td>
 
<td>&nbsp;</td>
 
<td>Saccharomyces</td>
 
</tr>
 
  
<tr class="s2">
+
Let's do something very quick in Chimera: calculate a Ramachandran plot.
<td>'''Species'''</td>
 
<td>&nbsp;</td>
 
<td>''Saccharomyces cerevisiae''</td>
 
</tr>
 
  
<table>
+
{{task|1=
 +
# Open Chimera and load <code>1BM8</code>.
 +
# Choose '''Favorites''' &rarr; '''Model Panel'''
 +
# Look for the Option '''Ramachandran plot...''' in the choices on the right.
 +
# Click the button and study the result. What do you see? What does this mean? What options are there? What do they do? What significance does this have?
 +
}}
  
  
You can see that there is not a common mapping between the yeast lineage and the commonly recognized categories - not all ranks are represented. Nor is this consistent across species in the taxonomic database: some have subfamily ranks and some don't. And the tree is in no way normalized - some of the ranks have thousands of members, and for some, only a single extant member may be known, or it may be a rank that only relates to the fossil record. But the ranks do provide some guidance to evolutionary divergence. Say you want to choose four species across the tree of life for a study, you should choose one from each of the major '''domains''' of life: Eubacteria, Euryarchaeota, Crenarchaeota-Eocytes, and Eukaryotes. Or you want to study a gene that is specific to mammals. Then you could choose from the clades listed in the NCBI taxonomy database under [http://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=40674&lvl=4 '''Mammalia'''] (a {{WP|Mammal_classification|'''class rank'''}}, and depending how many species you would want to include, use the
+
&nbsp;
subclass-, order-, or family rank (hover over the names to see their taxonomic rank.)
 
  
There will still be quite a bit of manual work involved and an exploration of different options on the Web may be useful. For our purposes here we can retrieve a good set of organisms from the [http://fungi.ensembl.org/info/website/species.html '''ensembl fungal genomes page'''] - maintained by the EBI's genome annotation group - that lists species grouped by taxonomic ''order''. All of these organisms are genome-sequenced, we can pick a set of representatives:
+
== CDD domain annotation ==
  
# Capnodiales&nbsp;&nbsp;&nbsp;''Zymoseptoria tritici''
+
In the last assignment, you followed a link to '''CDD Search Results''' from the [http://www.ncbi.nlm.nih.gov/protein/NP_010227 RefSeq record for yeast Mbp1] and briefly looked at the information offered by the NCBI's Conserved Domain Database, a database of ''Position Specific Scoring Matrices'' that embody domain definitions. Rather than access precomputed results, you can also search CDD with sequences: assuming you have saved the YFO Mbp1 sequence in FASTA format, this is straightforward. If you did not save this sequence, return to [[BIO_Assignment_Week_3|Assignment 3]] and retrieve it again.
# Erysiphales&nbsp;&nbsp;&nbsp;''Blumeria graminis''
 
# Eurotiales&nbsp;&nbsp;&nbsp;''Aspergillus nidulans''
 
# Glomerellales&nbsp;&nbsp;&nbsp;''Glomerella graminicola''
 
# Hypocreales&nbsp;&nbsp;&nbsp;''Trichoderma reesei''
 
# Magnaporthales&nbsp;&nbsp;&nbsp;''Magnaporthe oryzae''
 
# Microbotryales&nbsp;&nbsp;&nbsp;''Microbotryum violaceum''
 
# Pezizales&nbsp;&nbsp;&nbsp;''Tuber melanosporum''
 
# Pleosporales&nbsp;&nbsp;&nbsp;''Phaeosphaeria nodorum''
 
# Pucciniales&nbsp;&nbsp;&nbsp;''Puccinia graminis''
 
# Saccharomycetales&nbsp;&nbsp;&nbsp;''Saccharomyces cerevisiae''
 
# Schizosaccharomycetales&nbsp;&nbsp;&nbsp;''Schizosaccharomyces pombe''
 
# Sclerotiniaceae&nbsp;&nbsp;&nbsp;''Sclerotinia sclerotiorum''
 
# Sordariales&nbsp;&nbsp;&nbsp;''Neurospora crassa''
 
# Tremellales&nbsp;&nbsp;&nbsp;''Cryptococcus neoformans''
 
# Ustilaginales&nbsp;&nbsp;&nbsp;''Ustilago maydis''
 
  
This set of organisms thus can be used to generate a PSI-BLAST search in a well-distributed set of species. Of course '''you must also include YFO''' (<small>if YFO is not in this list already</small>).
 
  
To enter these 16 species as an Entrez restriction, they need to be formatted as below. (<small>One could also enter species one by one, by pressing the '''(+)''' button after the organism list</small>)
+
{{task|1=
 +
# Access the [http://www.ncbi.nlm.nih.gov/Structure/cdd/cdd.shtml '''CDD database'''] at http://www.ncbi.nlm.nih.gov/Structure/cdd/cdd.shtml
 +
# Read the information. CDD is a superset of various other database domain annotations as well as NCBI-curated domain definitions.
 +
# Copy the YFO Mbp1 FASTA sequence, paste it into the search form and click '''Submit'''.
 +
## On the result page, clik on '''View full result'''
 +
## Note that there are a number of partially overlapping ankyrin domain modules. We will study ankyrin domains in a later assignment.
 +
## Also note that there may be blocks of sequence colored cyan in the sequence bar. Hover your mouse over the blocks to see what these blocks signify.
 +
## Open the link to '''Search for similar domain architecture''' in a separate window and study it. This is the '''CDART''' database. Think about what these results may be useful for.
 +
## Click on one of the ANK superfamily graphics and see what the associated information looks like: there is a summary of structure and function, links to specific literature and a tree of the relationship of related sequences.
 +
}}
  
 +
== SMART domain annotation ==
  
<source lang="text">
 
Aspergillus nidulans[orgn]
 
OR Blumeria graminis[orgn]
 
OR Cryptococcus neoformans[orgn]
 
OR Glomerella graminicola[orgn]
 
OR Magnaporthe oryzae[orgn]
 
OR Microbotryum violaceum[orgn]
 
OR Neurospora crassa[orgn]
 
OR Phaeosphaeria nodorum[orgn]
 
OR Puccinia graminis[orgn]
 
OR Sclerotinia sclerotiorum[orgn]
 
OR Trichoderma reesei[orgn]
 
OR Tuber melanosporum[orgn]
 
OR Saccharomyces cerevisiae[orgn]
 
OR Schizosaccharomyces pombe[orgn]
 
OR Ustilago maydis[orgn]
 
OR Zymoseptoria tritici[orgn]
 
  
</source>
+
The [http://smart.embl-heidelberg.de/ SMART database] at the EMBL in Heidelberg offers an alternative view on domain architectures. I personally find it more useful for annotations because it integrates a number of additional features. You can search by sequence, or by accession number and that raises the question of how to retrieve a database cross-reference from an NCBI sequence identifier to a UniProt sequence ID:
  
  
 +
===ID mapping===
  
&nbsp;
 
  
==Executing the PSI-BLAST search==
+
{{task|
 +
<!-- Yeast:  NP_010227 ... P39678 -->
 +
# Access the [http://www.uniprot.org/mapping/ UniProt ID mapping service] at http://www.uniprot.org/mapping/
 +
# Paste the RefSeq identifier for YFO Mbp1 into the search field.
 +
# Use the menu to choose '''From''' ''RefSeq Protein'' and '''To''' ''UniprotKB AC''&ndash;the UniProt Knowledge Base Accession number.
 +
# Click on '''Map''' to execute the search.
 +
# Note the ID - it probably starts with a letter, followed by numbers and letters. Here are some examples for fungal Mbp1-like proteins: <code>P39678 Q5B8H6 Q5ANP5 P41412</code> ''etc.''
 +
# Click on the link, and explore how the UniProt sequence page is similar or different from the RefSeq page.
 +
}}
  
We have a list of species. Good. Next up: how do we '''use''' it.
+
===SMART search===
  
 
{{task|1=
 
{{task|1=
 +
# Access the [http://smart.embl-heidelberg.de/ '''SMART database'''] at http://smart.embl-heidelberg.de/
 +
# Click the lick to access SMART in the '''normal''' mode.
 +
# Paste the YFO Mbp1 UniProtKB Accession number into the '''Sequence ID or ACC''' field.
 +
# Check the boxes for:
 +
## '''PFAM domains''' (domains defined by sequence similarity in the PFAM database)
 +
## '''signal peptides''' (using the Gunnar von Heijne's SignalP 4.0 server at the Technical University in Lyngby, Denmark)
 +
## '''internal repeats''' (using the programs ''ariadne'' and ''prospero'' at the Wellcome Trust Centre for Human Genetics at Oxford University, England)
 +
## '''intrinsic protein disorder''' (using Rune Linding's DisEMBL program at the EMBL in Heidelberg, Germany)
 +
# Click on '''Sequence SMART''' to run the search and annotation. <small>(In case you get an error like: "Sorry, your entry seems to have no SMART domain ...", let me know and repeat the search with the actual FASTA sequence instead of the accession number.)</small>
  
 +
Study the results. Specifically, have a look at the proteins with similar domain '''ORGANISATION''' and '''COMPOSITION'''. This is similar to the NCBI's CDART.
  
# Navigate to the BLAST homepage.
 
# Select '''protein BLAST'''.
 
# Paste the APSES domain sequence into the search field.
 
# Select '''refseq''' as the database.
 
# Copy the organism restriction list from above '''and enter the correct name for YFO''' into the list if it is not there already. Obviously, you can't find sequences in YFO if YFO is not included in your search space. Paste the list into the '''Entrez Query''' field.
 
# In the '''Algorithm''' section, select PSI-BLAST.
 
#Click on '''BLAST'''.
 
 
}}
 
}}
 +
 +
=== Visual comparison of domain annotations in '''R''' ===
  
  
Evaluate the results carefully. Since we used default parameters, the threshold for inclusion was set at an '''E-value''' of 0.005 by default, and that may be a bit too lenient. If you look at the table of your hits&ndash; in the '''Sequences producing significant alignments...''' section&ndash; there may also be a few sequences that have a low query coverage of less than 80%. Let's exclude these from the profile initially: not to worry, if they are true positives, the will come back with improved E-values and greater coverage in subsequent iterations. But if they were false positives, their E-values will rise and they should drop out of the profile and not contaminate it.
+
The versatile plotting functions of '''R''' allow us to compare domain annotations. The distribution of segments that are annotated as being "low-complexity" or "disordered is particulalry interesting: these are functional features of the amino acid sequence that are often not associated with sequence similarity.
  
 +
In the following code tutorial, we create a plot similar to the CDD and SMART displays. It is based on the SMART domain annotations of the six-fungal reference species for the course.
  
 
{{task|1=
 
{{task|1=
#In the header section, click on '''Formatting options''' and in the line "Format for..." set the '''with inclusion threshold''' to <code>0.001</code> (This means E-values can't be above 10<sup>-03</sup> for the sequence to be included.)
 
# Click on the '''Reformat''' button (top right).
 
# In the table of sequence descriptions (not alignments!), click on the '''Query coverage''' to sort the table by coverage, not by score.
 
# Copy the rows with a coverage of less than 80% and paste them into some text editor so you can compare what happens with these sequences in the next iteration.
 
# '''Deselect''' the check mark next to these sequences in the right-hand column '''Select for PSI blast'''. (For me these are six sequences, but with YFO included that may be a bit different.)
 
# Then scroll to '''Run PSI-BLAST iteration 2 ...''' and click on '''<code>Go</code>'''.
 
}}
 
  
 +
Copy the code to an '''R''' script, study and execute it.
 +
<source lang="R">
 +
 +
# plotDomains
 +
# tutorial and functions to plot a colored rectangle from a list of domain annotations
 +
 +
 +
# First task: create a list structure for the annotations: this is a list of lists
 +
# As you see below, we need to mix strings, numbers and vectors of numbers. In R
 +
# such mixed data types must go into a list.
 +
 +
Mbp1Domains <- list()  # start with an empty list
  
This is now the "real" PSI-BLAST at work: it constructs a profile from all the full-length sequences and searches with the '''profile''', not with any individual sequence. Note that we are controlling what goes into the profile in two ways:
+
# For each species annotation, compile the SMART domain annotations in a list.
# we are explicitly removing sequences with poor coverage; and
+
Mbp1Domains <- rbind(Mbp1Domains, list(  # rbind() appends the list to the existing
# we are requiring a more stringent minimum E-value for each sequence.
+
    species = "Saccharomyces cerevisiae",
 +
    code    = "SACCE",
 +
    ACC    = "P39678",
 +
    length  = 833,
 +
    KilAN  = c(18,102),  # Note: Vector of (start, end) pairs
 +
    AThook  = NULL,      # Note: NULL, because this annotation was not observed in this sequence
 +
    Seg    = c(108,122,236,241,279,307,700,717),
 +
    DisEMBL = NULL,
 +
    Ankyrin = c(394,423,427,463,512,541),  # Note: Merge overlapping domains, if present
 +
    Coils  = c(633, 655)
 +
    ))
  
 +
Mbp1Domains <- rbind(Mbp1Domains, list(
 +
    species = "Emericella nidulans",
 +
    code    = "ASPNI",
 +
    ACC    = "Q5B8H6",
 +
    length  = 695,
 +
    KilAN  = c(23,94),
 +
    AThook  = NULL,
 +
    Seg    = c(529,543),
 +
    DisEMBL = NULL,
 +
    Ankyrin = c(260,289,381,413),
 +
    Coils  = c(509,572)
 +
    ))
  
{{task|1=
+
Mbp1Domains <- rbind(Mbp1Domains, list(
#Again, study the table of hits. Sequences highlighted in yellow have met the search criteria in the second iteration. Note that the coverage of (some) of the previously excluded sequences is now above 80%.
+
    species = "Candida albicans",
# Let's exclude partial matches one more time. Again, deselect all sequences with less than 80% coverage. Then run the third iteration.
+
    code    = "CANAL",
# Iterate the search in this way until no more "New" sequences are added to the profile. Then scan the list of excluded hits ... are there any from YFO that seem like they could potentially make the list? Marginal E-value perhaps, or reasonable E-value but less coverage? If that's the case, try returning the E-value threshold to the default 0.005 and see what happens...
+
    ACC    = "Q5ANP5",
}}
+
    length  = 852,
 +
    KilAN  = c(19,102),
 +
    AThook  = NULL,
 +
    Seg    = c(351,365,678,692),
 +
    DisEMBL = NULL,
 +
    Ankyrin = c(376,408,412,448,516,545),
 +
    Coils  = c(665,692)
 +
    ))
  
 +
Mbp1Domains <- rbind(Mbp1Domains, list(
 +
    species = "Neurospora crassa",
 +
    code    = "NEUCR",
 +
    ACC    = "Q7RW59",
 +
    length  = 833,
 +
    KilAN  = c(31,110),
 +
    AThook  = NULL,
 +
    Seg    = c(130,141,253,266,514,525,554,564,601,618,620,629,636,652,658,672,725,735,752,771),
 +
    DisEMBL = NULL,
 +
    Ankyrin = c(268,297,390,419),
 +
    Coils  = c(500,550)
 +
    ))
  
Once no "new" sequences have been added, if we were to repeat the process again and again, we would always get the same result because the profile stays the same. We say that the search has '''converged'''. Good. Time to harvest.
+
Mbp1Domains <- rbind(Mbp1Domains, list(
 +
    species = "Schizosaccharomyces pombe",
 +
    code    = "SCHPO",
 +
    ACC    = "P41412",
 +
    length  = 657,
 +
    KilAN  = c(21,97),
 +
    AThook  = NULL,
 +
    Seg    = c(111,125,136,145,176,191,422,447),
 +
    DisEMBL = NULL,
 +
    Ankyrin = c(247,276,368,397),
 +
    Coils  = c(457,538)
 +
    ))
  
 +
Mbp1Domains <- rbind(Mbp1Domains, list(
 +
    species = "Ustilago maydis",
 +
    code    = "USTMA",
 +
    ACC    = "Q4P117",
 +
    length  = 956,
 +
    KilAN  = c(21,98),
 +
    AThook  = NULL,
 +
    Seg    = c(106,116,161,183,657,672,776,796),
 +
    DisEMBL = NULL,
 +
    Ankyrin = c(245,274,355,384),
 +
    Coils  = c(581,609)
 +
    ))
  
{{task|1=
 
# At the header, click on '''Taxonomy reports''' and find YFO in the '''Organism Report''' section. These are your APSES domain homologs. All of them. Actually, perhaps more than all: the report may also include sequences with E-values above the inclusion threshold.
 
# From the report copy the sequence identifiers
 
## from YFO,
 
## with E-values above your defined threshold.
 
}}
 
  
For example, the list of ''Saccharomyces'' genes is the following:
+
# Working with data in lists and dataframes can be awkward, since the result
 +
# of filters and slices are themselves lists, not vectors.
 +
# Therefore we need to use the unlist() function a lot. When in doubt: unlist()
  
<code>
+
#### Boxes #####
<b>[http://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=559292 Saccharomyces cerevisiae S288c]</b> [http://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=4890 [ascomycetes]] taxid 559292<br \>
+
# Define a function to draw colored boxes, given input of a vector with
[http://www.ncbi.nlm.nih.gov:80/entrez/query.fcgi?cmd=Retrieve&db=Protein&list_uids=6320147&dopt=GenPept ref|NP_010227.1|] Mbp1p [Saccharomyces cerevisiae S288c]          [ 131] 1e-38<br \>
+
# (start,end) pairs, a color, and the height where the box should go.
[http://www.ncbi.nlm.nih.gov:80/entrez/query.fcgi?cmd=Retrieve&db=Protein&list_uids=6320957&dopt=GenPept ref|NP_011036.1|] Swi4p [Saccharomyces cerevisiae S288c]         [ 123]  1e-35<br \>
+
drawBoxes <- function(v, c, h) {  # vector of xleft, xright pairs; color; height
[http://www.ncbi.nlm.nih.gov:80/entrez/query.fcgi?cmd=Retrieve&db=Protein&list_uids=6322808&dopt=GenPept ref|NP_012881.1|] Phd1p [Saccharomyces cerevisiae S288c]          [  91]  1e-25<br \>
+
    if (is.null(v)) { return() }
[http://www.ncbi.nlm.nih.gov:80/entrez/query.fcgi?cmd=Retrieve&db=Protein&list_uids=6323658&dopt=GenPept ref|NP_013729.1|] Sok2p [Saccharomyces cerevisiae S288c]          [  93]  3e-25<br \>
+
    for (i in seq(1,length(v),by=2)) {
[http://www.ncbi.nlm.nih.gov:80/entrez/query.fcgi?cmd=Retrieve&db=Protein&list_uids=6322090&dopt=GenPept ref|NP_012165.1|] Xbp1p [Saccharomyces cerevisiae S288c]          [  40]  5e-07<br \>
+
        rect(v[i], h-0.1, v[i+1], h+0.1, border="black", col=c)
</code>
+
    }
 +
}
  
[[Saccharomyces cerevisiae Xbp1|Xbp1]] is a special case. It has only very low coverage, but that is because it has a long domain insertion and the N-terminal match often is not recognized by alignment because the gap scores for long indels are unrealistically large. For now, I keep that sequence with the others.
+
#### Annotations ####
 +
# Define a function to write the species code, draw a grey
 +
# horizontal line and call drawBoxes() for every annotation type
 +
# in the list
 +
drawGene <- function(rIndex) {
 +
    # define colors:
 +
    kil <- "#32344F"
 +
    ank <- "#691A2C"
 +
    seg <- "#598C9E"
 +
    coi <- "#8B9998"
 +
    xxx <- "#EDF7F7"
 +
   
 +
    text (-30, rIndex, adj=1, labels=unlist(Mbp1Domains[rIndex,"code"]), cex=0.75 )
 +
    lines(c(0, unlist(Mbp1Domains[rIndex,"length"])), c(rIndex, rIndex), lwd=3, col="#999999")
  
 +
    drawBoxes(unlist(Mbp1Domains[rIndex,"KilAN"]),  kil, rIndex)
 +
    drawBoxes(unlist(Mbp1Domains[rIndex,"AThook"]),  xxx, rIndex)
 +
    drawBoxes(unlist(Mbp1Domains[rIndex,"Seg"]),    seg, rIndex)
 +
    drawBoxes(unlist(Mbp1Domains[rIndex,"DisEMBL"]), xxx, rIndex)
 +
    drawBoxes(unlist(Mbp1Domains[rIndex,"Ankyrin"]), ank, rIndex)
 +
    drawBoxes(unlist(Mbp1Domains[rIndex,"Coils"]),  coi, rIndex)
 +
}
  
Next we need to retrieve the sequences. Tedious to retrieve them one by one, but we can get them all at the same time:
+
#### Plot ####
 +
# define the size of the plot-frame
 +
yMax <- length(Mbp1Domains[,1])  # number of domains in list
 +
xMax <- max(unlist(Mbp1Domains[,"length"]))  # largest sequence length
  
 +
# plot an empty frame
 +
plot(1,1, xlim=c(-100,xMax), ylim=c(0, yMax) , type="n", yaxt="n", bty="n", xlab="sequence number", ylab="")
  
{{task|1=
+
# Finally, iterate over all species and call drawGene()
 +
for (i in 1:length(Mbp1Domains[,1])) {
 +
    drawGene(i)
 +
}
  
# Return to the BLAST results page and again open the '''Formatting options'''.
+
# end
# Find the '''Limit results''' section and enter YFO's name into the field. For example <code>Saccharomyces cerevisiae [ORGN]</code>
 
# Click on '''Reformat'''
 
# Scroll to the '''Descriptions''' section, check the box at the left-hand margin, next to each sequence you want to keep. Then click on '''Download &rarr; FASTA complete sequence &rarr; Continue'''.
 
  
 +
</source>
 +
}}
  
  
<div class="mw-collapsible mw-collapsed" data-expandtext="Expand" data-collapsetext="Collapse" style="border:#000000 solid 1px; padding: 10px; margin-left:25px; margin-right:25px;">There are actually several ways to download lists of sequences. Using the results page utility is only one. But if you know the GIs of the sequences you need, you can get them more directly by putting them into the URL...
+
When you execute the code, your plot should look similar to this one:
<div  class="mw-collapsible-content">
 
  
* http://www.ncbi.nlm.nih.gov/protein/6320147,6320957,6322808,6323658,6322090?report=docsum  - The default report
+
[[Image:DomainAnnotations.jpg|frame|none|SMART domain annotations for Mbp1 proteins from six fungal species.
* http://www.ncbi.nlm.nih.gov/protein/6320147,6320957,6322808,6323658,6322090?report=fasta - FASTA sequences with NCBI HTML markup
+
]]
  
Even more flexible is the [http://www.ncbi.nlm.nih.gov/books/NBK25500/#chapter1.Downloading_Full_Records '''eUtils'''] interface to the NCBI databases. For example you can download the dataset in text format by clicking below.
+
{{task|1=
  
* http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=protein&id=6320147,6320957,6322808,6323658,6322090&rettype=fasta&retmode=text
+
On the Student Wiki, add the annotations for YFO to the plot:
  
Note that this utility does not '''show''' anything, but downloads the (multi) fasta file to your default download directory.
+
# Copy one of the list definitions for Mbp1 domains and edit it with the appropriate values for your own annotations.
   
+
# Test that you can add the YFO annotation to the plot.
</div>
+
# Submit your validated code block to the [http://biochemistry.utoronto.ca/steipe/abc/students/index.php/BCH441_2014_Assignment_4_domain_annotations '''Student Wiki here''']. The goal is to compile an overview of all species we are studying in class.   
</div>
+
# If your working annotation block is in the Wiki before noontime on Wednesday, you will be awarded a 10% bonus on the quiz.
 
}}
 
}}
 +
  
  
Line 405: Line 546:
 
<!-- <div class="reference-box">[http://www.ncbi.nlm.nih.gov]</div> -->
 
<!-- <div class="reference-box">[http://www.ncbi.nlm.nih.gov]</div> -->
  
<!--
 
 
Add to this assignment:
 
- study the BLAST output format, links, tools, scores ...
 
- compare the improvement in E-values to standard BLAST
 
- examine this in terms of sensitivity and specificity
 
 
-->
 
  
 
&nbsp;
 
&nbsp;

Revision as of 14:32, 2 October 2015

Assignment for Week 6
Function

< Assignment 5 Assignment 7 >

Note! This assignment is currently inactive. Major and minor unannounced changes may be made at any time.

 
 

Concepts and activities (and reading, if applicable) for this assignment will be topics on next week's quiz.




 

Introduction

 

In this assignment we will perform a few computations with coordinate data in PDB files, look more in depth at domain annotations, and compare them across different proteins related to yeast Mbp1. We will write a function in R to help us plot a graphic of the comparison and collaborate to share data for our plots. But first we need to go through a few more R concepts.


 

Programming R code

First, we will cover essentials of R programming: the fundamental statements that are needed to write programs–conditional expressions and loops, and how to define functions that allow us to use programming code. But let's start with a few more data types of R so we can use the concepts later on: matrices, lists and data frames.


Task:

Please begin by working through the short R - tutorial: matrices section and the following sections on "Lists" and "Data frames".

Note that what we have done here is just the bare minimum on vectors, matrices and lists. The concepts are very generally useful, and there are many useful ways to extract subsets of values. We'll come across these in various parts of R sample code. But unless you read the provided code examples line by line, make sure you understand every single statement and ask if you are not clear about the syntax, this will all be quickly forgotten. Use it or lose it!


R is a full featured programming language and in order to be that, it needs the ability to manipulate Control Flow, i.e. the order in which instructions are executed. By far the most important control flow statement is a conditional branch i.e. the instruction to execute code at alternative points of the program, depending on the presence or absence of a condition. Using such conditional branch instructions, two main types of programming constructs can be realized: conditional expressions and loops.

Conditional expressions

The template for conditional expression in R is:

if( <expression 1> ) {
  <statement 1>
} 
else if ( <expresssion 2> ) {
  <statement 2>
} 
else {
  <statement 3>
}

...where both the else if (...) { ... } and the else (...) block are optional. We have encountered this construct previously, when we assigned the appropriate colors for amino acids in the frequency plot:

if      (names(logRatio[i]) == "F") { barColors[i] <- hydrophobic }
else if (names(logRatio[i]) == "G") { barColors[i] <- plain }
   [... etc ...]
else                                { barColors[i] <- plain }

Logical expressions

We have to consider the <expression> in a bit more detail: anything that is, or produces, or can be interpreted as, a Boolean TRUE or FALSE value can serve as an expression in a conditional statement.

Task:
Here are some examples. Copy the code to an R script, predict what will happen on in each line and try it out:

# A boolean constant is interpreted as is:
if (TRUE)    {print("true")} else {print("false")}
if (FALSE)   {print("true")} else {print("false")}

# Strings of "true" and "false" are coerced to their 
# Boolean equivalent, but contrary to other programming languages
# arbitrary, non-empty or empty strings are not interpreted.
if ("true") {print("true")} else {print("false")}
if ("false") {print("true")} else {print("false")}
if ("widdershins")    {print("true")} else {print("false")}
if ("") {print("true")} else {print("false")}

# All non-zero, defined numbers are TRUE
if (1)      {print("true")} else {print("false")}
if (0)      {print("true")} else {print("false")}
if (-1)     {print("true")} else {print("false")}
if (pi)     {print("true")} else {print("false")}
if (NULL)   {print("true")} else {print("false")}
if (NA)     {print("true")} else {print("false")}
if (NaN)    {print("true")} else {print("false")}
if (Inf)    {print("true")} else {print("false")}

# functions can return Boolean values
affirm <- function() { return(TRUE) }
deny <- function() { return(FALSE) }
if (affirm())    {print("true")} else {print("false")}
if (deny())    {print("true")} else {print("false")}

# N.B. coercion of Booleans into numbers can be done as well
and is sometimes useful: consider ...
a <- c(TRUE, TRUE, FALSE, TRUE, FALSE)
a
as.numeric(a)
sum(a)

# ... or coercing the other way ...
as.logical(-1:1)

Logical operators

To actually write a conditional statement, we have to be able to test a condition and this is what logical operators do. Is something equal to something else? Is it less? Does something exist? Is it a number?

Task:
Here are some examples. Again, predict what will happen ...

TRUE             # Just a statement.

#   unary operator
! TRUE           # NOT ...

# binary operators
FALSE > TRUE     # GREATER THAN ...
FALSE < TRUE     # ... these are coerced to numbers
FALSE < -1        
0 == FALSE       # Careful! == compares, = assigns!!!

"x" == "u"       # using lexical sort order ...
"x" >= "u"
"x" <= "u"
"x" != "u"
"aa" > "u"       # ... not just length, if different.
"abc" < "u"  

TRUE | FALSE    # OR: TRUE if either is true
TRUE & FALSE    # AND: TRUE if both are TRUE

# equality and identity
?identical
a <- c(TRUE)
b <- c(TRUE)
a; b
a == b
identical(a, b)

b <- 1
a; b
a == b
identical(a, b)   # Aha: equal, but not identical


# some other useful tests for conditional expressions
?all
?any
?duplicated
?exists
?is.character
?is.factor
?is.integer
?is.null
?is.numeric
?is.unsorted
?is.vector


Loops

Loops allow you to repeat tasks many times over. The template is:

for (<name> in <vector>) {
   <statement>
}

Task:
Consider the following: Again, copy the code to a script, study it, predict what will happen and then run it.

# simple for loop
for (i in 1:10) {
	print(c(i, i^2, i^3))
}

# Compare excution times: one million square roots from a vector ...
n <- 1000000
x <- 1:n
y <- sqrt(x)

# ... or done explicitly in a for-loop
for (i in 1:n) {
  y[i] <- sqrt (x[i])
}

If you can achieve your result with an R vector expression, it will be faster than using a loop. But sometimes you need to do things explicitly, for example if you need to access intermediate results.


Here is an example to play with loops: a password generator. Passwords are a pain. We need them everywhere, they are frustrating to type, to remember and since the cracking programs are getting smarter they become more and more likely to be broken. Here is a simple password generator that creates random strings with consonant/vowel alterations. These are melodic and easy to memorize, but actually as strong as an 8-character, fully random password that uses all characters of the keyboard such as )He.{2jJ or #h$bB2X^ (which is pretty much unmemorizable). The former is taken from 207 * 77 1015 possibilities, the latter is from 948 ~ 6*1015 possibilities. HIgh-end GPU supported password crackers can test about 109 passwords a second, the passwords generated by this little algorithm would thus take on the order of 106 seconds or eleven days to crack[1]. This is probably good enough to deter a casual attack.

Task:
Copy, study and run ...

# Suggest memorizable passwords
# Below we use the functions: 
?nchar
?sample
?substr
?paste
?print

#define a string of  consonants ...
con <- "bcdfghjklmnpqrstvwxz"
# ... and a string of of vowels
vow <- "aeiouy"

for (i in 1:10) {  # ten sample passwords to choose from ...
    pass = rep("", 14)  # make an empty character vector
    for (j in 1:7) {    # seven consonant/vowel pairs to be created ...
        k   <- sample(1:nchar(con), 1)  # pick a random index for consonants ...
        ch  <- substr(con,k,k)          #  ... get the corresponding character ...
        idx <- (2*j)-1                  # ... compute the position (index) of where to put the consonant ...
        pass[idx] <- ch                 # ...  and put it in the right spot

        # same thing for the vowel, but coded with fewer intermediate assignments
        # of results to variables
        k <- sample(1:nchar(vow), 1)
        pass[(2*j)] <- substr(vow,k,k) 
    }
    print( paste(pass, collapse="") )  # collapse the vector in to a string and print
}

Functions

Finally: functions. Functions look very much like the statements we have seen above. the template looks like:

<name> <- function (<parameters>) {
   <statements>
}

In this statement, the function is assigned to the name - any valid name in R. Once it is assigned, it the function can be invoked with name(). The parameter list (the values we write into the parentheses followin the function name) can be empty, or hold a list of variable names. If variable names are present, you need to enter the corresponding parameters when you execute the function. These assigned variables are available inside the function, and can be used for computations. This is called "passing the variable into the function".

You have encountered a function to choose YFO names. In this function, your Student ID was the parameter. Here is another example to play with: a function that calculates how old you are. In days. This is neat - you can celebrate your 10,000 birthday - or so.

Task:
Copy, explore and run ...

Define the function ...
# A lifedays calculator function

myLifeDays <- function(date = NULL) { # give "date" a default value so we can test whether it has been set
    if (is.null(date)) {
        print ("Enter your birthday as a string in \"YYYY-MM-DD\" format.")
        return()
    }
    x <- strptime(date, "%Y-%m-%d") # convert string to time
    y <- format(Sys.time(), "%Y-%m-%d") # convert "now" to time
    diff <- round(as.numeric(difftime(y, x, unit="days")))
    print(paste("This date was ", diff, " days ago."))
}
Use the function (example)
   myLifeDays("1932-09-25")  # Glenn Gould's birthday

Here is a good opportunity to play and practice programming: modify this function to accept a second argument. When a second argument is present (e.g. 10000) the function should print the calendar date on which the input date will be that number of days ago. Then you could use it to know when to celebrate your 10,000th lifeDay, or your 777th anniversary day or whatever.

Enjoy.

Protein structure annotation

Let's do something very quick in Chimera: calculate a Ramachandran plot.

Task:

  1. Open Chimera and load 1BM8.
  2. Choose FavoritesModel Panel
  3. Look for the Option Ramachandran plot... in the choices on the right.
  4. Click the button and study the result. What do you see? What does this mean? What options are there? What do they do? What significance does this have?


 

CDD domain annotation

In the last assignment, you followed a link to CDD Search Results from the RefSeq record for yeast Mbp1 and briefly looked at the information offered by the NCBI's Conserved Domain Database, a database of Position Specific Scoring Matrices that embody domain definitions. Rather than access precomputed results, you can also search CDD with sequences: assuming you have saved the YFO Mbp1 sequence in FASTA format, this is straightforward. If you did not save this sequence, return to Assignment 3 and retrieve it again.


Task:

  1. Access the CDD database at http://www.ncbi.nlm.nih.gov/Structure/cdd/cdd.shtml
  2. Read the information. CDD is a superset of various other database domain annotations as well as NCBI-curated domain definitions.
  3. Copy the YFO Mbp1 FASTA sequence, paste it into the search form and click Submit.
    1. On the result page, clik on View full result
    2. Note that there are a number of partially overlapping ankyrin domain modules. We will study ankyrin domains in a later assignment.
    3. Also note that there may be blocks of sequence colored cyan in the sequence bar. Hover your mouse over the blocks to see what these blocks signify.
    4. Open the link to Search for similar domain architecture in a separate window and study it. This is the CDART database. Think about what these results may be useful for.
    5. Click on one of the ANK superfamily graphics and see what the associated information looks like: there is a summary of structure and function, links to specific literature and a tree of the relationship of related sequences.

SMART domain annotation

The SMART database at the EMBL in Heidelberg offers an alternative view on domain architectures. I personally find it more useful for annotations because it integrates a number of additional features. You can search by sequence, or by accession number and that raises the question of how to retrieve a database cross-reference from an NCBI sequence identifier to a UniProt sequence ID:


ID mapping

Task:

  1. Access the UniProt ID mapping service at http://www.uniprot.org/mapping/
  2. Paste the RefSeq identifier for YFO Mbp1 into the search field.
  3. Use the menu to choose From RefSeq Protein and To UniprotKB AC–the UniProt Knowledge Base Accession number.
  4. Click on Map to execute the search.
  5. Note the ID - it probably starts with a letter, followed by numbers and letters. Here are some examples for fungal Mbp1-like proteins: P39678 Q5B8H6 Q5ANP5 P41412 etc.
  6. Click on the link, and explore how the UniProt sequence page is similar or different from the RefSeq page.

SMART search

Task:

  1. Access the SMART database at http://smart.embl-heidelberg.de/
  2. Click the lick to access SMART in the normal mode.
  3. Paste the YFO Mbp1 UniProtKB Accession number into the Sequence ID or ACC field.
  4. Check the boxes for:
    1. PFAM domains (domains defined by sequence similarity in the PFAM database)
    2. signal peptides (using the Gunnar von Heijne's SignalP 4.0 server at the Technical University in Lyngby, Denmark)
    3. internal repeats (using the programs ariadne and prospero at the Wellcome Trust Centre for Human Genetics at Oxford University, England)
    4. intrinsic protein disorder (using Rune Linding's DisEMBL program at the EMBL in Heidelberg, Germany)
  5. Click on Sequence SMART to run the search and annotation. (In case you get an error like: "Sorry, your entry seems to have no SMART domain ...", let me know and repeat the search with the actual FASTA sequence instead of the accession number.)

Study the results. Specifically, have a look at the proteins with similar domain ORGANISATION and COMPOSITION. This is similar to the NCBI's CDART.

Visual comparison of domain annotations in R

The versatile plotting functions of R allow us to compare domain annotations. The distribution of segments that are annotated as being "low-complexity" or "disordered is particulalry interesting: these are functional features of the amino acid sequence that are often not associated with sequence similarity.

In the following code tutorial, we create a plot similar to the CDD and SMART displays. It is based on the SMART domain annotations of the six-fungal reference species for the course.

Task:
Copy the code to an R script, study and execute it.

# plotDomains
# tutorial and functions to plot a colored rectangle from a list of domain annotations


# First task: create a list structure for the annotations: this is a list of lists
# As you see below, we need to mix strings, numbers and vectors of numbers. In R
# such mixed data types must go into a list.

Mbp1Domains <- list()   # start with an empty list

# For each species annotation, compile the SMART domain annotations in a list.
Mbp1Domains <- rbind(Mbp1Domains, list(  # rbind() appends the list to the existing
    species = "Saccharomyces cerevisiae",
    code    = "SACCE",
    ACC     = "P39678",
    length  = 833,
    KilAN   = c(18,102),  # Note: Vector of (start, end) pairs
    AThook  = NULL,       # Note: NULL, because this annotation was not observed in this sequence
    Seg     = c(108,122,236,241,279,307,700,717),
    DisEMBL = NULL,
    Ankyrin = c(394,423,427,463,512,541),   # Note: Merge overlapping domains, if present
    Coils   = c(633, 655)
    ))

Mbp1Domains <- rbind(Mbp1Domains, list(
    species = "Emericella nidulans",
    code    = "ASPNI",
    ACC     = "Q5B8H6",
    length  = 695,
    KilAN   = c(23,94),
    AThook  = NULL,
    Seg     = c(529,543),
    DisEMBL = NULL,
    Ankyrin = c(260,289,381,413),
    Coils   = c(509,572)
    ))

Mbp1Domains <- rbind(Mbp1Domains, list(
    species = "Candida albicans",
    code    = "CANAL",
    ACC     = "Q5ANP5",
    length  = 852,
    KilAN   = c(19,102),
    AThook  = NULL,
    Seg     = c(351,365,678,692),
    DisEMBL = NULL,
    Ankyrin = c(376,408,412,448,516,545),
    Coils   = c(665,692)
    ))

Mbp1Domains <- rbind(Mbp1Domains, list(
    species = "Neurospora crassa",
    code    = "NEUCR",
    ACC     = "Q7RW59",
    length  = 833,
    KilAN   = c(31,110),
    AThook  = NULL,
    Seg     = c(130,141,253,266,514,525,554,564,601,618,620,629,636,652,658,672,725,735,752,771),
    DisEMBL = NULL,
    Ankyrin = c(268,297,390,419),
    Coils   = c(500,550)
    ))

Mbp1Domains <- rbind(Mbp1Domains, list(
    species = "Schizosaccharomyces pombe",
    code    = "SCHPO",
    ACC     = "P41412",
    length  = 657,
    KilAN   = c(21,97),
    AThook  = NULL,
    Seg     = c(111,125,136,145,176,191,422,447),
    DisEMBL = NULL,
    Ankyrin = c(247,276,368,397),
    Coils   = c(457,538)
    ))

Mbp1Domains <- rbind(Mbp1Domains, list(
    species = "Ustilago maydis",
    code    = "USTMA",
    ACC     = "Q4P117",
    length  = 956,
    KilAN   = c(21,98),
    AThook  = NULL,
    Seg     = c(106,116,161,183,657,672,776,796),
    DisEMBL = NULL,
    Ankyrin = c(245,274,355,384),
    Coils   = c(581,609)
    ))


# Working with data in lists and dataframes can be awkward, since the result
# of filters and slices are themselves lists, not vectors.
# Therefore we need to use the unlist() function a lot. When in doubt: unlist()

#### Boxes #####
# Define a function to draw colored boxes, given input of a vector with
# (start,end) pairs, a color, and the height where the box should go.
drawBoxes <- function(v, c, h) {   # vector of xleft, xright pairs; color; height
    if (is.null(v)) { return() }
    for (i in seq(1,length(v),by=2)) {
        rect(v[i], h-0.1, v[i+1], h+0.1, border="black", col=c)
    }
}

#### Annotations ####
# Define a function to write the species code, draw a grey
# horizontal line and call drawBoxes() for every annotation type
# in the list
drawGene <- function(rIndex) {
    # define colors:
    kil <- "#32344F"
    ank <- "#691A2C"
    seg <- "#598C9E"
    coi <- "#8B9998"
    xxx <- "#EDF7F7"
    
    text (-30, rIndex, adj=1, labels=unlist(Mbp1Domains[rIndex,"code"]), cex=0.75 )
    lines(c(0, unlist(Mbp1Domains[rIndex,"length"])), c(rIndex, rIndex), lwd=3, col="#999999")

    drawBoxes(unlist(Mbp1Domains[rIndex,"KilAN"]),   kil, rIndex)
    drawBoxes(unlist(Mbp1Domains[rIndex,"AThook"]),  xxx, rIndex)
    drawBoxes(unlist(Mbp1Domains[rIndex,"Seg"]),     seg, rIndex)
    drawBoxes(unlist(Mbp1Domains[rIndex,"DisEMBL"]), xxx, rIndex)
    drawBoxes(unlist(Mbp1Domains[rIndex,"Ankyrin"]), ank, rIndex)
    drawBoxes(unlist(Mbp1Domains[rIndex,"Coils"]),   coi, rIndex)
}

#### Plot ####
# define the size of the plot-frame
yMax <- length(Mbp1Domains[,1])   # number of domains in list
xMax <- max(unlist(Mbp1Domains[,"length"]))  # largest sequence length

# plot an empty frame
plot(1,1, xlim=c(-100,xMax), ylim=c(0, yMax) , type="n", yaxt="n", bty="n", xlab="sequence number", ylab="")

# Finally, iterate over all species and call drawGene()
for (i in 1:length(Mbp1Domains[,1])) {
    drawGene(i)
}

# end


When you execute the code, your plot should look similar to this one:

SMART domain annotations for Mbp1 proteins from six fungal species.

Task:
On the Student Wiki, add the annotations for YFO to the plot:

  1. Copy one of the list definitions for Mbp1 domains and edit it with the appropriate values for your own annotations.
  2. Test that you can add the YFO annotation to the plot.
  3. Submit your validated code block to the Student Wiki here. The goal is to compile an overview of all species we are studying in class.
  4. If your working annotation block is in the Wiki before noontime on Wednesday, you will be awarded a 10% bonus on the quiz.


That is all.


 

Links and resources

 


Footnotes and references

  1. That's assuming the worst case in that the attacker needs to know the pattern with which the password is formed, i.e. the number of characters and the alphabet that we chose from. But note that there is an even worse case: if the attacker had access to our code and the seed to our random number generator. When the random number generator starts up, a new seed is generated from system time, thus the possible space of seeds can be devastatingly small. But even if a seed is set explicitly with the set.seed() function, that seed is a 32-bit integer and thus can take only a bit more than 4*109 values, six orders of magnitude less than the 1015 password complexity we thought we had. It turns out that the code may be a much greater vulnerability than the password itself. Keep that in mind. Keep it secret. Keep it safe.


 

Ask, if things don't work for you!

If anything about the assignment is not clear to you, please ask on the mailing list. You can be certain that others will have had similar problems. Success comes from joining the conversation.



< Assignment 5 Assignment 7 >