View source for Bioinformatics Introduction Sequence

<div id="BIO">
<div class="b1">
Sequence<br />
</div>

<table style="width:100%;"><tr>
<td style="height:30px; align:center; vertical-align:middle;">[[Bioinformatics_Introduction_Data|Data]]</td>
<td style="height:30px; align:center; vertical-align:middle;">''Sequence''</td>
<td style="height:30px; align:center; vertical-align:middle;">[[Bioinformatics_Introduction_Structure|Structure]]</td>
<td style="height:30px; align:center; vertical-align:middle;">[[Bioinformatics_Introduction_Phylogeny|Phylogeny]]</td>
<td style="height:30px; align:center; vertical-align:middle;">[[Bioinformatics_Introduction_Function|Function]]</td>
</tr></table>

{{Vspace}}

<!-- div class="alert">

Warning – this page is currently under construction (2016-12-26).

Do not use before this warning has been removed.

</div -->

{{Vspace}}

__TOC__


{{Vspace}}

==The Sequence Unit==

This Unit is part of a brief introduction to bioinformatics. The material is more or less interleaved with the <code>Sequence.R</code> Project File which is part of the RStudio project associated with this material. Refer to the course/workshop page for installation instructions.

This Unit introduces a target organism which has recently been genome-sequenced in which we will search for novel information about APSES domain transcription factors, it discusses pairwise- and multiple sequence alignment, and BLAST for fast database searches, it explores some methods for sequence analysis and introduces work with entire genomes.

{{Vspace}}

{{Vspace}}

==Introduction==

{{Vspace}}

<div class="colmask doublepage">
  <div class="colleft">
    <div class="col1">
      <!-- Column 1 start -->
Sequence alignment is a '''very''' large, and important topic.

One of the foundations of bioinformatics is the empirical observation that related sequences conserve structure, and often function. Much of what we know about a protein's physiological function is based on the '''conservation''' of that function as the species evolves. Indeed, conservation is a defining aspect of what can rightly be said to be a protein's "function" in the first place. 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 multiple effects on a species' reproductive fitness. 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.

We assess conservation by comparing sequences between related proteins. This is the basis on which we can make inferences from well-studied model organisms for species that have not been studied as deeply. The foundation is to measure protein sequence similarity. If two sequences are much more similar than we could expect from chance, we hypothesize that their similarity comes from shared ancestry plus conservation. The measurement of sequence similarity however requires sequence alignment<ref>This is not strictly true in all cases: some algorithms measure similarity through an alignment-free approach, for example by comparing structural features, or domain annotations. These methods are less sensitive, but important when sequences are so highly diverged that no meaningful sequence alignment can be produced.</ref>.

     <!-- Column 1 end -->
    </div>
    <div class="col2">
      <!-- Column 2 start -->

A carefully done sequence alignment is a cornerstone for the annotation of the essential properties a gene or protein. It can already tell us a lot about which proteins we expect to have similar functions in different species. 

Multiple sequence alignments ('''MSAs''') are further useful to resolve ambiguities in the precise placement of "indels"<ref>"indel": '''in'''sertion / '''del'''etion – a difference in sequence length between two aligned sequences that is accommodated by gaps in the alignment. Since we can't tell from the comparison of two sequences whether such a change was introduced by ''insertion into'' or ''deletion from'' the ancestral sequence, we join both into a {{WP|Portmanteau|''portmanteau''}}.</ref> and to ensure that columns in alignments actually contain amino acids that evolve in a similar context. MSAs serve as input for 
* functional annotation;
* protein homology modelling; 
* phylogenetic analyses, and
* sensitive homology searches in databases.

Below we will explore the essentials of 

<div class="emphasis-box">
* optimal global and local pairwise alignment;
* Fast BLAST searches to determine best matches in large databases, and reciprocal best matches;
* PSI BLAST searches for exhaustive matches; 
* Domain annotation by sequence alignment to statistical models;
* Multiple sequence alignments; and
* Genomes.
</div>


This is the scenario: the genome sequence of early-diverging fungus ''Spizellomyces punctatus'' has been recently published<ref>{{#pmid: 27540072}}</ref>. Does it contain proteins that are related to yeast Mbp1? If there are more, which one is the most closely related protein? Is its DNA binding domain conserved? How can we identify '''all''' related genes in ''S. punctatus''? And, what can we learn from that collection of sequences?
      <!-- Column 2 end -->
    </div>
  </div>
</div>

{{Vspace}}



==An ''S. punctatus'' homologue to Mbp1==

{{Vspace}}

BLAST is a fast algorithm for searching related sequences in large databases. We will explore it in more detail below, for now, we just briefly use it to find the sequence in ''S. punctatus'' that is most closely related to <code>MBP1_SACCE</code>.

{{Vspace}}

{{task|1=

# Return to the [http://www.ncbi.nlm.nih.gov/protein/NP_010227 Mbp1 protein page at the NCBI] and follow the link to [https://blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE=Proteins&PROGRAM=blastp&BLAST_PROGRAMS=blastp&QUERY=NP_010227.1&LINK_LOC=protein&PAGE_TYPE=BlastSearch Run BLAST] under "Analyze this sequence" in the right-hand column.
# This allows you to perform a sequence similarity search. You need to set two parameters:
## As '''Database''', select '''Reference proteins (refseq_protein)''' from the drop down menu;
## In the '''Organism''' field, type <code>spizellomyces punctatus</code>. The field will autocomplete to two choices, choose the DAOM BR117 strain (taxid:645134), as that is the strain which has been sequenced.
# Click on '''BLAST''' to start the search. This should find a handful of genes, all of them in ''S. punctatus''. If you find none, or hundreds, or they are not all in the same species, you did something wrong. Ask on the mailing list and make sure to fix the problem.
# Look at the top "hit" in the '''Alignments section'''. It should be the ''hypothetical protein SPPG_00022'', and it should have two matching '''Ranges''' reported. 
# In the header information for each hit is a link to its database entry, right next to '''Sequence ID'''.  It says <code>[https://www.ncbi.nlm.nih.gov/protein/1026947531 XP_016612325.1]</code> ... follow that link.
# Open the [https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=645134 Taxonmoy Browser] page in a new tab (it's linked from the <code>ORGANISM</code> line of the sequernce record), explore it, and note the taxonomy ID.
# Note the RefSeq ID, and the search results %ID, E-value, whether one or more similar regions were found etc. in your Journal. We will refer to this sequence as "''S. punctatus'' Mbp1" or <code>MBP1_SPIPU</code> or similar in the future.
# Finally access the [http://www.uniprot.org/uploadlists/ UniProt ID mapping service] to retrieve the UniProt ID for the protein. Paste the RefSeq ID and choose '''RefSeq Protein''' as the '''From:''' option and '''UniProtKB''' as the '''To:''' option.

:<small>If the mapping works, the UniProt ID will be in the '''Entry:''' column of the table that is being returned. Click the link and have a look at the UniProt entry page while you're there.</small>

<!-- What could go wrong? Sometimes the mapping does not work: I.e. UniProt contains the sequence, but the mapping service does not know.
I don't know why the mapping for some sequences is not available.
If this happens, you can work around the problem as follows.

1. Load the RefSeq protein page
2. View the protein as FASTA and copy the sequence.
3. Open the UniProt BLAST page http://www.uniprot.org/blast/
  (Yes, UniProt runs its own BLAST version, and that searches UniProt databases, not Genbank)
4. Paste the sequence into the search form and run BLAST.

... if the sequence is in UniProt, you will get the top hit with 100% sequence identity.

-->

}}

{{Vspace}}



===Add MBP1_SPIPU to the database===

{{Vspace}}

{{task|1=

* (Re)-open the RStudio project and open the <code>Sequence.R</code> file.
* Work through <code>PART ONE: ADDING DATA</code>.

}}

{{Vspace}}

{{Vspace}}

==Analyze==

Let us perform a few simple sequence analyses using the online EMBOSS tools. EMBOSS (the European Molecular Biology laboratory Open Software Suite) combines a large number of simple but fundamental sequence analysis tools. The tools can be [http://emboss.sourceforge.net/download/ installed locally on your own machine], or run via a public Web interface. Google for [http://www.google.ca/search?q=EMBOSS+explorer EMBOSS explorer], public access points include http://emboss.bioinformatics.nl/ .

Access an EMBOSS Explorer service and familiarize yourself with the offerings in the EMBOSS package. I find some of the nucleic acid tools indispensable in the lab, such as restriction-site mapping tools, and I frequently use the alignment tools <code>Needle</code> and <code>Water</code>, but by and large the utility of many of the components&ndash;while fast, efficient and straightforward to use&ndash; suffers from lack of reference and comparison and from terse output. The routines show their conceptual origin in the 1970s and 1980s. We will encounter alternatives in later assignments.

{{Vspace}}

{{task|1=
;Local composition
# Find <code>pepinfo</code> under the '''PROTEIN COMPOSITION''' heading.
# Retrieve the <code>MBP1_SACCE</code> and <code>MBP1_SPIPU</code> sequences from your '''R''' database, e.g. with something like <br /><code>   cat(db$protein[db$protein$name == "MBP1_SACCE"), "sequence"]</code>
# Do the following for both sequences in separate windows:
## Copy and paste the sequence into the input field.
## Run with default parameters.
## Scroll to the figures all the way at the bottom.
# Try to compare the result ... <small>(kind of hard without reference, overlay and expectation, isn't it?)</small>
}}


{{task|1=
;Global composition
# Find <code>pepstats</code> under the '''PROTEIN COMPOSITION''' heading.
# Do the following for both sequences in separate windows:
## Paste the <code>MBP1_SACCE</code> resp. the <code>MBP1_SPIPU</code> sequence into the input field.
## Run with default parameters.
# Try to compare ... are there significant and unexpected differences?
}}


{{task|1=
;Transmembrane sequences
# Find <code>tmap</code>. Also find <code>shuffleseq</code>.
# Use the <code>MBP1_SPIPU</code> sequence to annotate transmembrane helices in the protein and at least five shuffled sequences. <code>MBP1_SPIPU</code> is not expected to have TM helices, nor are the shuffled sequences expected to have any. If you '''do''' find some, these are most likely "''false positives''".

# Also compare the following positive control: Gef1 - a yeast chloride channel with 10 trans-membrane helices and outward localized N-terminus:
<source lang="text">
>gi|6322500|ref|NP_012574.1| Gef1p [Saccharomyces cerevisiae S288c]
MPTTYVPINQPIGDGEDVIDTNRFTNIPETQNFDQFVTIDKIAEENRPLSVDSDREFLNSKYRHYREVIW
DRAKTFITLSSTAIVIGCIAGFLQVFTETLVNWKTGHCQRNWLLNKSFCCNGVVNEVTSTSNLLLKRQEF
ECEAQGLWIAWKGHVSPFIIFMLLSVLFALISTLLVKYVAPMATGSGISEIKVWVSGFEYNKEFLGFLTL
VIKSVALPLAISSGLSVGKEGPSVHYATCCGYLLTKWLLRDTLTYSSQYEYITAASGAGVAVAFGAPIGG
VLFGLEEIASANRFNSSTLWKSYYVALVAITTLKYIDPFRNGRVILFNVTYDRDWKVQEIPIFIALGIFG
GLYGKYISKWNINFIHFRKMYLSSWPVQEVLFLATLTALISYFNEFLKLDMTESMGILFHECVKNDNTST
FSHRLCQLDENTHAFEFLKIFTSLCFATVIRALLVVVSYGARVPAGIFVPSMAVGATFGRAVSLLVERFI
SGPSVITPGAYAFLGAAATLSGITNLTLTVVVIMFELTGAFMYIIPLMIVVAITRIILSTSGISGGIADQ
MIMVNGFPYLEDEQDEEEEETLEKYTAEQLMSSKLITINETIYLSELESLLYDSASEYSVHGFPITKDED
KFEKEKRCIGYVLKRHLASKIMMQSVNSTKAQTTLVYFNKSNEELGHRENCIGFKDIMNESPISVKKAVP
VTLLFRMFKELGCKTIIVEESGILKGLVTAKDILRFKRIKYREVHGAKFTYNEALDRRCWSVIHFIIKRF
TTNRNGNVI
</source>

# Evaluate the output: does the algorithm (wrongly) predict TM-helices in <code>MBP1_SPIPU</code>? In the shuffled sequences? Does it find all ten TM-helices in Gef1?
}}


{{task|1=
;Motifs
# Find <code>pepcoil</code>, an algorithm to detect {{WP|coiled coil}} motifs.
# Run this with the YFO Mbp1 sequence and yeast Mbp1.
# Try to compare ... do both sequences have coiled-coil motif predictions? Are they annotated in approximately comparable regions of the respective sequence?
}}


A textbook case of a {{WP|Coiled_coil|coiled coil}} is the so-called {{WP|Leucine_zipper|leucine zipper}} domain. It is a protein-protein interaction module that has a number of leucine residues in i, i+7 position. We don't know what the intervening residues may be, thus we need to define the sequence we are looking for as a generic pattern. This is efficiently done using {{WP|Regular_expression|'''regular expressions'''}}. 

{{Vspace}} 

====A Brief Digression to Regular Expressions====

{{Vspace}}

;Regular expressions are a concise description language to define patterns for pattern-matching in strings.

Truth be told, many programmers have a love-hate relationship with regular expressions. The syntax of regular expressions is very powerful and expressive, but also terse, not always intuitive, and sometimes hard to understand. I'll introduce you to a few principles here that are quite straightforward and they will probably cover 99% of the cases you will encounter.

Here is our test-case: the sequence of Mbp1, copied from the [https://www.ncbi.nlm.nih.gov/protein/NP_010227 NCBI Protein database page for yeast Mbp1].

        1 msnqiysary sgvdvyefih stgsimkrkk ddwvnathil kaanfakakr trilekevlk
       61 ethekvqggf gkyqgtwvpl niakqlaekf svydqlkplf dftqtdgsas pppapkhhha
      121 skvdrkkair sastsaimet krnnkkaeen qfqsskilgn ptaaprkrgr pvgstrgsrr
      181 klgvnlqrsq sdmgfprpai pnssisttql psirstmgpq sptlgileee rhdsrqqqpq
      241 qnnsaqfkei dledglssdv epsqqlqqvf nqntgfvpqq qssliqtqqt esmatsvsss
      301 pslptspgdf adsnpfeerf pgggtspiis miprypvtsr pqtsdindkv nkylsklvdy
      361 fisnemksnk slpqvllhpp phsapyidap idpelhtafh wacsmgnlpi aealyeagts
      421 irstnsqgqt plmrsslfhn sytrrtfpri fqllhetvfd idsqsqtvih hivkrksttp
      481 savyyldvvl skikdfspqy rielllntqd kngdtalhia skngdvvffn tlvkmgaltt
      541 isnkegltan eimnqqyeqm miqngtnqhv nssntdlnih vntnnietkn dvnsmvimsp
      601 vspsdyityp sqiatnisrn ipnvvnsmkq masiyndlhe qhdneikslq ktlksisktk
      661 iqvslktlev lkesskdeng eaqtnddfei lsrlqeqntk klrkrliryk rlikqkleyr
      721 qtvllnklie detqattnnt vekdnntler lelaqeltml qlqrknklss lvkkfednak
      781 ihkyrriire gtemnieevd ssldvilqtl iannnknkga eqiitisnan sha
 //


{{task|1=

Navigate to http://regexpal.com and paste the sequence into the '''lower''' box. This site is one of a number of online regular expression testers; their immediate, visual feedback is invaluable when you are developing regular expression patterns. Set the dialect to '''PCRE''' in the drop-down menu, it is '''Javascript''' by default.

Lets try some expressions:

;Most characters are matched literally.
:Type "<code>a</code>" in to the '''upper''' box and you will see all "<code>a</code>" characters matched. Then replace <code>a</code> with <code>q</code>. 
: Now type "<code>aa</code>" instead. Then <code>krnnkk</code>. ''Sequences'' of characters are also matched literally.

;The pipe character {{pipe}} that symbolizes logical OR can be used to define that more than one character should match:
:<code>i(s{{pipe}}m{{pipe}}q)n</code> matches <code>isn</code> OR <code>imn</code> OR <code>iqn</code>. Note how we can group with parentheses, and try what would happen without them.

;We can more conveniently specify more than one character to match if we place it in square brackets.
:<code>[lq]</code> matches <code>l</code> OR <code>q</code>. <code>[familyvw]</code> matches hydrophobic amino acids. 

;Within square brackets, we can specify "ranges".
:<code>[1-5]</code> matches digits from 1 to 5.

;Within square brackets, we can specify characters that should NOT be matched, with the caret, <code>^</code>.
:<code>[^0-9]</code> matches everything EXCEPT digits. <code>[^a-z]</code> matches everything that is not a lower-case letter. That's what we need (try it).

}}


{{Vspace}}

{{task|1 =
* Work through <code>PART TWO: REGULAR EXPRESSIONS</code> section of the <code>Sequence.R</code> '''R''' script for a deeper exposure to regular expressions.
}}

{{Vspace}}

{{Vspace}}

=='''R''' Sequence Analysis Tools==

It's interesting to see this collection of tools that were carefully designed some twenty years ago, as an open source replacement for a set of software tools - the GCG package - that was indispensable for molecular biology labs in the 80s and 90s, but whose cost had become prohibitive. Fundamentally this is a building block approach, and the field has turned to programming solutions instead.

As for functionality, much more sophisticated functions are available on the Web: do take a few minutes and browse the [http://bioinformatics.ca/links_directory/tag/protein-sequence-analysis curated Web services directory of bioinformatics.ca].

As for versatility, '''R''' certainly has the edge. Let's explore some of the functions available in the <code>seqinr</code> package that you already encountered in the introductory [[R tutorial]]. They are comparatively basic - but it is easy to prepare our own analysis.


{{Vspace}}

{{task|1 =
* Study the code in the <code>PART THREE: SEQUENCE ANALYSIS</code> section of the <code>Sequence.R</code> script
}}

{{Vspace}}

{{Vspace}}

==Alignment==

{{Vspace}}

;Sequence alignments are the single most important task of bioinformatics algorithms.

{{Vspace}}

=== DotPlots and the Mutation Data Matrix ===

Before we start calculating alignments, we should get a better sense of the underlying sequence similarity. A Dotplot is a perfect tool for that, because it displays alignment-free similarity information. Let's make a dotplot that uses the BLOSUM62 Mutation Data Matrix to measure pairwise amino acid similarity. The NCBI makes its alignment matrices available by ftp. They are located at  ftp://ftp.ncbi.nih.gov/blast/matrices - for example here is a link to the [ftp://ftp.ncbi.nih.gov/blast/matrices/BLOSUM62 '''BLOSUM62 matrix''']<ref>That directory also contains sourcecode to generate the PAM matrices. This may be of interest if you ever want to produce scoring matrices from your own datasets.</ref>. 

{{Vspace}}

Scoring matrices are also available in the Bioconductor Biostrings package.

<source lang="text">

BLOSUM62

   A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V  B  J  Z  X  *
A  4 -1 -2 -2  0 -1 -1  0 -2 -1 -1 -1 -1 -2 -1  1  0 -3 -2  0 -2 -1 -1 -1 -4
R -1  5  0 -2 -3  1  0 -2  0 -3 -2  2 -1 -3 -2 -1 -1 -3 -2 -3 -1 -2  0 -1 -4
N -2  0  6  1 -3  0  0  0  1 -3 -3  0 -2 -3 -2  1  0 -4 -2 -3  4 -3  0 -1 -4
D -2 -2  1  6 -3  0  2 -1 -1 -3 -4 -1 -3 -3 -1  0 -1 -4 -3 -3  4 -3  1 -1 -4
C  0 -3 -3 -3  9 -3 -4 -3 -3 -1 -1 -3 -1 -2 -3 -1 -1 -2 -2 -1 -3 -1 -3 -1 -4
Q -1  1  0  0 -3  5  2 -2  0 -3 -2  1  0 -3 -1  0 -1 -2 -1 -2  0 -2  4 -1 -4
E -1  0  0  2 -4  2  5 -2  0 -3 -3  1 -2 -3 -1  0 -1 -3 -2 -2  1 -3  4 -1 -4
G  0 -2  0 -1 -3 -2 -2  6 -2 -4 -4 -2 -3 -3 -2  0 -2 -2 -3 -3 -1 -4 -2 -1 -4
H -2  0  1 -1 -3  0  0 -2  8 -3 -3 -1 -2 -1 -2 -1 -2 -2  2 -3  0 -3  0 -1 -4
I -1 -3 -3 -3 -1 -3 -3 -4 -3  4  2 -3  1  0 -3 -2 -1 -3 -1  3 -3  3 -3 -1 -4
L -1 -2 -3 -4 -1 -2 -3 -4 -3  2  4 -2  2  0 -3 -2 -1 -2 -1  1 -4  3 -3 -1 -4
K -1  2  0 -1 -3  1  1 -2 -1 -3 -2  5 -1 -3 -1  0 -1 -3 -2 -2  0 -3  1 -1 -4
M -1 -1 -2 -3 -1  0 -2 -3 -2  1  2 -1  5  0 -2 -1 -1 -1 -1  1 -3  2 -1 -1 -4
F -2 -3 -3 -3 -2 -3 -3 -3 -1  0  0 -3  0  6 -4 -2 -2  1  3 -1 -3  0 -3 -1 -4
P -1 -2 -2 -1 -3 -1 -1 -2 -2 -3 -3 -1 -2 -4  7 -1 -1 -4 -3 -2 -2 -3 -1 -1 -4
S  1 -1  1  0 -1  0  0  0 -1 -2 -2  0 -1 -2 -1  4  1 -3 -2 -2  0 -2  0 -1 -4
T  0 -1  0 -1 -1 -1 -1 -2 -2 -1 -1 -1 -1 -2 -1  1  5 -2 -2  0 -1 -1 -1 -1 -4
W -3 -3 -4 -4 -2 -2 -3 -2 -2 -3 -2 -3 -1  1 -4 -3 -2 11  2 -3 -4 -2 -2 -1 -4
Y -2 -2 -2 -3 -2 -1 -2 -3  2 -1 -1 -2 -1  3 -3 -2 -2  2  7 -1 -3 -1 -2 -1 -4
V  0 -3 -3 -3 -1 -2 -2 -3 -3  3  1 -2  1 -1 -2 -2  0 -3 -1  4 -3  2 -2 -1 -4
B -2 -1  4  4 -3  0  1 -1  0 -3 -4  0 -3 -3 -2  0 -1 -4 -3 -3  4 -3  0 -1 -4
J -1 -2 -3 -3 -1 -2 -3 -4 -3  3  3 -3  2  0 -3 -2 -1 -2 -1  2 -3  3 -3 -1 -4
Z -1  0  0  1 -3  4  4 -2  0 -3 -3  1 -1 -3 -1  0 -1 -2 -2 -2  0 -3  4 -1 -4
X -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -4
* -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4  1

</source>


{{task|
* Study this table and make sure you understand what this information represents, how it can be used, and what a reasonable range of values for identities and pairscores for non-identical, similar and dissimilar residues is. Ask on the mailing list in case you have questions. '''This piece of data is the foundation of any sequence alignment. without it, no sensible alignment could be produced!'''
* Figure out the following values:
** Compare an identical match of histidine with an identical match of serine. What does this mean?
** How similar are lysine and leucine, as compared to leucine and isoleucine? Is this what you expect?
** PAM matrices are sensitive to an interesting artefact. Since W and R can be interchanged with a single point mutation, the probability of observing W&rarr;R and R&rarr;W exchanges in '''closely related sequences''' is much higher than one would expect from the two amino acid's biophysical properties. (Why?) PAM matrices were compiled from hypothetical point exchanges and then extrapolated. Therefore these matrices assign a relatively high degree of similarity to (W, R), that is not warranted considering what actually happens in nature. Do you see this problem in the BLOSUM matrix? If BLOSUM does not have this issue, why not? 
}}

{{Vspace}}

Next, let's apply the scoring matrix for actual comparison:

{{Vspace}}

{{task|1 =

* Return to your RStudio session.
* Study and work through the code in the <code>PART FOUR: DOTPLOT AND MDM</code> section of the <code>Sequence.R</code> script

}}

{{Vspace}}

{{Vspace}}

== Pairwise Alignments: Optimal ==

{{Vspace}}

Optimal pairwise sequence alignment is the mainstay of sequence comparison. To consider such alignments in practice, we'll align the same sequences that we have just mapped in the dotplot exercise: <code>MBP1_SACCE</code> and <code>MBP1_SPIPU</code>. Your dotplots should have shown you two regions of similarity: a highly similar region focussed somewhere around the N-terminal 100 amino acids, and a more extended, but somewhat less similar region in the middle of the sequences. You can think of the sequence alignment algorithm as building the similarity matrix, and then discovering the best path along high-scoring diagonals.

{{Vspace}}

=== Optimal Sequence Alignment: EMBOSS online tools===

{{Vspace}}

Online programs for optimal sequence alignment are part of the EMBOSS tools. The programs take FASTA files or raw text files as input. 

'''Local''' optimal sequence alignment using "water" 

{{task|1=
# Fetch the sequences for <code>MBP1_SACCE</code> and <code>MBP1_SPIPU</code> from your database.
# Access the [http://emboss.bioinformatics.nl/ EMBOSS Explorer site] (if you haven't done so yet, you might want to bookmark it.)
# Look for '''ALIGNMENT LOCAL''', click on '''water''', paste your sequences and run the program with default parameters.
# Study the results. You will probably find that the alignment extends over most of the protein, but does not include the termini.
# Change the '''Gap opening''' and '''Gap extension''' parameters to high values (e.g. 30 and 5). Then run the alignment again.
# Note what is different.
}}


'''Global''' optimal sequence alignment using "needle" 
{{task|1=
# Look for '''ALIGNMENT GLOBAL''', click on '''needle''', paste the <code>MBP1_SACCE</code> and <code>MBP1_SPIPU</code> sequences again and run the program with default parameters.
# Study the results. You will find that the alignment extends over the entire protein, likely with long ''indels'' at the termini.
}}

{{Vspace}}

=== Optimal Sequence Alignment with '''R''': Biostrings ===

{{Vspace}}

Biostrings has extensive functions for sequence alignments. They are generally well written and tightly integrated with the rest of Bioconductor's functions. There are a few quirks however: for example alignments won't work with lower-case sequences<ref>While this seems like an unnecessary limitation, given that we could easily write such code to transform to-upper when looking up values in the MDM, perhaps it is meant as an additional sanity check that we haven't inadvertently included text in the sequence that does not belong there, such as the FASTA header line perhaps.</ref>.


{{Vspace}}

{{task|1 =

* Return to your RStudio session.
* Study and work through the code in the <code>PART FIVE: BIOSTRINGS PAIRWISE ALIGNMENT</code> section of the <code>Sequence.R</code> script

}}

{{Vspace}}

{{Vspace}}

==Heuristic pairwise alignments: BLAST==

{{Vspace}}


<div class="colmask doublepage">
  <div class="colleft">
    <div class="col1">
      <!-- Column 1 start -->
[http://www.ncbi.nlm.nih.gov/blast '''BLAST'''] is by a margin the most important computational tool of molecular biology. It is so important, that we have already used BLAST above, even before properly introducing the algorithm and the principles, to find the most similar sequence to <code>MBP1_SACCE</code> in SPIPU. 

In this part of the unit we will use BLAST to perform '''Reciprocal Best Matches'''. 

One of the important questions of ''model-organism based inference'' is: which genes perform the same function in two different organisms. In the absence of other information, our best guess is that these are the two genes that are '''mutually''' most similar. The keyword here is '''mutually'''. If <code>MBP1_SACCE</code> from ''S. cerevisiae'' is the best match to <code>RES2_SCHPO</code> in ''S. pombe'', the two proteins are only mutually most similar if  <code>RES2_SCHPO</code> is more similar to <code>MBP1_SACCE</code> than to any other ''S. cerevisiae'' protein. We call this a '''Reciprocal Best Match''', or "RBM"<ref>Note that RBMs are usually orthologues, but the definition of orthologue and RBM is not the same. Most importantly, many orthologues are not RBMs. We will explore this more when we discuss phylogenetic inference.</ref>.

The argument is summarized in the figure on the right: genes that evolve under continuos selective pressure on their function have relatively lower mutation rates and are thus more similar to each other, than genes that undergo neo- or sub-functionalization after duplication.

However, there is a catch: proteins are often composed of multiple domains that implement distinct roles of their function. Under the assumptions above we could hypothesize:
* a gene in SPIPU that has the "same" function as the Mbp1 cell-cycle checkpoint switch in yeast should be an RBM to Mbp1;
* a gene that binds to the same DNA sites as Mbp1 should have a DNA-binding domain that is an RBM to the DNA binding domain of Mbp1.

Thus we'll compare RBMs in SPIPU for full-length <code>Mbp1_SACCE</code> and its DNA-binding domain, and see if the results are the same.


      <!-- Column 1 end -->
    </div>
    <div class="col2">
      <!-- Column 2 start -->
[[Image:RBM.jpg|frame|A hypothetical phylogenetic gene tree. "S" is a speciation in the tree, "D" is a duplication within a species. The duplicated gene (teal triangle) evolves towards a different function and thus acquires more mutations than its paralogue (teal circle). If an RBM search start from the blue triangle, it finds the red circle. However the reciprocal match finds the teal circle. The red and teal circles fulfill the '''RBM''' criterion.]]

      <!-- Column 2 end -->
    </div>
  </div>
</div>


{{Vspace}}

===Full-length RBM===

{{Vspace}}

You have already performed the first half of the experiment: matching from ''S. cerevisiae'' to ''S. punctatus'' . The backward match is simple.

{{task|1=
# Access [http://www.ncbi.nlm.nih.gov/blast '''BLAST'''] and follow the link to the '''protein blast''' program.
# Enter the RefSeq ID for <code>MBP1_SPIPU</code> in the '''Query sequence''' field.
# Select <code>refseq_protein</code> as the '''database''' to search in, and enter <code>Saccharomyces cerevisiae (taxid:4932)</code> to restrict the '''organism''' for which hits are reported.
# Run BLAST. Examine the results.

If your top-hit is <code>NP_010227</code>, you have confirmed the RBM between <code>Mbp1_SACCE</code> and <code>Mbp1_SPIPU</code>. If it is not, let me know. I expect this to be the same and would like to verify your results if it is not.

}}

{{Vspace}}

===RBM for the DNA binding domain===

{{Vspace}}

The DNA-binding domain of  <code>Mbp1_SACCE</code> is called an '''APSES''' domain. If the RBM between ''Saccharomyces cerevisiae'' Mbp1 and ''Spizellomyces punctatus'' is truly an orthologue, we expect all of the protein's respective domains to have the RBM property as well. But let's not simply assume what we can easily test. We'll define the sequence of the APSES domain in MBP1_SACCE and YFO and see how these definitions reflect in a BLAST search.

{{Vspace}}

====Defining the range of the APSES domain annotation====

{{#lst:Reference annotation yeast Mbp1|CDD_APSES}}


{{Vspace}}

====Executing the forward search====

{{Vspace}}

{{task|1=
# Access [http://www.ncbi.nlm.nih.gov/blast '''BLAST'''] and follow the link to the '''protein blast''' program.
# '''Forward search:'''
## Paste only the APSES domain sequence for <code>MBP1_SACCE</code> in the '''Query sequence''' field (copy the sequence from above).
## Select <code>refseq_protein</code> as the '''database''' to search in, and enter the correct taxonomy ID for ''Spizellomyces punctatus'' in the '''Organism''' field.
## Run BLAST. Examine the results.
## If the top hit is the same protein you have already seen, oK. If it's not, something went wrong.

}}

With this we have confirmed the sequence with the most highly conserved APSES domain in ''Spizellomyces punctatus''. Can we take the sequence for the reverse search from the alignment that BLAST returns? Actually, that is not a good idea. The BLAST alignment is not guaranteed to be optimal. We should do an optimal sequence alignment instead. That is: we use two different tools here for two different purposes: we use BLAST to identify one protein as the most similar among many alternatives and we use optimal sequence alignment to determine the best alignment between two sequences. That best alignment is what we will annotate as the ''Spizellomyces punctatus'' APSES domain.

{{Vspace}}

====Alignment to define the YFO APSES domain for the reverse search====

{{Vspace}}


{{task|1 =

* Return to your RStudio session.
* Study and work through the code in the <code>PART SIX: APSES DOMAIN ANNOTATION BY ALIGNMENT</code> section of the <code>Sequence.R</code> script.

}}

{{Vspace}}

====Executing the reverse search====

{{Vspace}}

{{task|1=
# Paste the the APSES domain sequence for the ''Spizellomyces punctatus'' protein and enter it into '''Query sequence''' field of the BLAST form.
## Select <code>refseq_protein</code> as the '''database''' to search in, and enter <code>Saccharomyces cerevisiae (taxid:4932)</code> to restrict the '''organism''' for which hits are reported.
## Run BLAST. Examine the results.

If your top-hit is again <code>NP_010227</code>, you have confirmed the RBM between the APSES domain of <code>Mbp1_SACCE</code> and <code>Mbp1_SPIPU</code>. If it is not, something went wrong.
}}


{{Vspace}}

{{Vspace}}

==Heuristic profile-based alignment: PSI BLAST==

{{Vspace}}

<div class="colmask doublepage">
  <div class="colleft">
    <div class="col1">
      <!-- Column 1 start -->
It is (deceptively) easy to perform BLAST searches via the Web interface, but to use such powerful computational tools to their greatest advantage takes a considerable amount of care, caution and consideration.

PSI-BLAST allows to perform very sensitive searches for homologues that have diverged so far that their pairwise sequence similarity has become insignificant. It achieves this by establishing a profile of sequences to align with the database, rather than searching with individual sequences. This deemphasizes parts of the sequence that are variable and inconsequential, and focusses on the parts of greater structural and functional importance. As a consequence, the signal to noise ratio is greatly enhanced. 

In this part of the 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 for multiple alignments, calculation of conservation ''etc''. 

The first methodical problem we have to address is what sequence to search with. The full-length Mbp1 sequence from ''Saccharomyces cerevisiae'' or its RBM from ''Spizellomyces punctatus'' are not suitable: They contain multiple domains (in particular the ubiquitous Ankyrin domains) and would create broad, non-'''specific''' profiles. The APSES domain sequence by contrast is structurally well defined. The KilA-N domain, being shorter, is less likely to make a '''sensitive''' profile. Indeed one of the results of our analysis will be  to find whether APSES domains in fungi all have the same length as the Mbp1 domain, or whether some are indeed much shorter, like the KILA-N domain, as suggested by the Pfam alignment.

The second methodical problem we must address is how to perform a sensitive PSI-BLAST search '''in one organism'''. We need to balance two conflicting objectives:

* If we restrict the PSI-BLAST search to ''Spizellomyces punctatus'', PSI-BLAST has little chance of building a meaningful profile - the number of homologues that actually are '''in''' ''S. punctatus'' is too small. Thus the search will not become very sensitive.

      <!-- Column 1 end -->
    </div>
    <div class="col2">
      <!-- Column 2 start -->

* If we don't restrict our search, but search in all species, the number of hits may become unwieldily large. It becomes increasingly difficult to closely check all hits as to whether they have good coverage. Also we need to evaluate the fringe cases of marginal E-value: should a new sequence be added to the profile, or should we hold off on it for one or two iterations, to see whether its E-value drops significantly. By all means, we need to avoid profile corruption. 

Perhaps this is still be manageable when we are searching in fungi, but imagine you are working with a bacterial protein, or a protein that is conserved across the entire tree of life: your search may find tens of thousands of sequences. And by next year, thousands more will have been added. 

Therefore we have to find a middle ground: add enough organisms (sequences) to compile a sensitive profile, but not so many that we can no longer individually assess the sequences that contribute to the profile. We need to define a broadly representative but manageable set of species - to exploit the transitivity of homology - even if we are interested only in matches in one species: YFO. Please reflect on this and make sure you understand why we include sequences in a PSI-BLAST search that we are not actually interested in.

We need a subset of species
# that represent as large a '''range''' as possible on the evolutionary tree; 
# that are as well '''distributed''' as possible on the tree; and 
# whose '''genomes''' are fully sequenced.



      <!-- Column 2 end -->
    </div>
  </div>
</div>


&nbsp;




{{Vspace}}

===Selecting species for a PSI-BLAST search===

{{Vspace}}


<div class="colmask doublepage">
  <div class="colleft">
    <div class="col1">
      <!-- Column 1 start -->

To select species, we will use an approach that is conceptually simple: select 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. 

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 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, 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.)

<table>

<tr class="sh">
<td>Rank</td>
<td>Suffix</td>
<td>Example</td>
</tr>

<tr class="s1">
<td>Domain</td>
<td>&nbsp;</td>
<td>Eukaryota (Eukarya)</td>
</tr>

<tr class="s2">
<td>&nbsp;&nbsp;Subdomain</td>
<td>&nbsp;</td>
<td>Opisthokonta</td>
</tr>

<tr class="s1">
<td>'''Kingdom'''</td>
<td>&nbsp;</td>
<td>Fungi</td>
</tr>

<tr class="s2">
<td>&nbsp;&nbsp;Subkingdom</td>
<td>&nbsp;</td>
<td>Dikarya</td>
</tr>

<tr class="s1">
<td>'''Phylum'''</td>
<td>&nbsp;</td>
<td>Ascomycota</td>
</tr>

<tr class="s2">
<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">
<td>&nbsp;&nbsp;Subphylum</td>
<td>-mycotina</td>
<td>Saccharomycotina</td>
</tr>

<tr class="s2">
<td>'''Class'''</td>
<td>-mycetes</td>
<td>Saccharomycetes</td>
</tr>

<tr class="s1">
<td>&nbsp;&nbsp;Subclass</td>
<td>-mycetidae</td>
<td>&nbsp;</td>
</tr>

<tr class="s2">
<td>'''Order'''</td>
<td>-ales</td>
<td>Saccharomycetales</td>
</tr>

<tr class="s1">
<td>'''Family'''</td>
<td>-aceae</td>
<td>Saccharomycetaceae</td>
</tr>

<tr class="s2">
<td>&nbsp;&nbsp;Subfamily</td>
<td>-oideae</td>
<td>&nbsp;</td>
</tr>

<tr class="s1">
<td>&nbsp;&nbsp;Tribe</td>
<td>-eae</td>
<td>&nbsp;</td>
</tr>

<tr class="s2">
<td>&nbsp;&nbsp;Subtribe</td>
<td>-ineae</td>
<td>&nbsp;</td>
</tr>

<tr class="s1">
<td>'''Genus'''</td>
<td>&nbsp;</td>
<td>Saccharomyces</td>
</tr>

<tr class="s2">
<td>'''Species'''</td>
<td>&nbsp;</td>
<td>''Saccharomyces cerevisiae''</td>
</tr>

</table>


      <!-- Column 1 end -->
    </div>
    <div class="col2">
      <!-- Column 2 start -->

You can see that there is no common mapping between the yeast lineage listed at the NCBI 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 
subclass-, order-, or family rank (hover over the names to see their taxonomic rank.)

I have chosen the 10 species below to define a well-distributed search-space for PSI-BLAST. Of course '''we must also include ''Spizellomyces punctatus'' in the selection'''.

To enter these 10 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>)

<source lang="text">

"Wallemia mellicola"[organism] OR
"Puccinia Graminis"[organism] OR
"Ustilago maydis"[organism] OR
"Cryptococcus neoformans"[organism] OR
"Coprinopsis cinerea"[organism] OR
"Schizosaccharomyces pombe"[organism] OR
"Aspergillus nidulans"[organism] OR
"Neurospora crassa"[organism] OR
"Bipolaris oryzae"[organism] OR
"Saccharomyces cerevisiae"[organism] OR
"Spizellomyces punctatus"[organism]
</source>


      <!-- Column 2 end -->
    </div>
  </div>
</div>


{{Vspace}}

===Executing the PSI-BLAST search===

{{Vspace}}


We have a list of species. Good. Next up: how do we '''use''' it.

{{task|1=


# Navigate to the [http://www.ncbi.nlm/nih.gov/blast/ BLAST homepage].
# Select '''protein BLAST'''.
# Paste the MBP1_SACCE APSES domain sequence into the search field:
 >APSES_MBP1 Residues 4-102 of S. cerevisiae Mbp1
 QIYSARYSGVDVYEFIHSTGSIMKRKKDDWVNATHILKAANFAKAKRTRI
 LEKEVLKETHEKVQGGFGKYQGTWVPLNIAKQLAEKFSVYDQLKPLFDF
# Select '''refseq''' as the database.
# Copy the Entrez restrictions from above. Paste the list into the '''Entrez Query''' field.
# In the '''Algorithm''' section, select PSI-BLAST.
#Click on '''BLAST'''.
}}


Evaluate the results carefully. Since we did not change the algorithm parameters, the threshold for inclusion was set at an '''E-value''' of 0.005 by default, and that may be a bit too lenient, i.e. it might include sequences that are not homologous. 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, they 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 will drop out of the profile and not corrupt it.


{{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 '''Query cover''' to sort the table by coverage, not by score. 
# '''Deselect''' the check mark next to these sequences in the second-to-rightmost column '''Select for PSI blast'''. 
# Then scroll to '''Run PSI-BLAST iteration 2 ...''' and click on '''<code>Go</code>'''.
}}


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:
# we are explicitly removing sequences with poor coverage; and
# we are requiring a more stringent minimum E-value for each sequence.


{{task|1=
#Again, study the table of hits. Sequences highlighted in yellow have met the search criteria in the second iteration and are proposed for inclusion in the next iteration. Note that the coverage of (some) of the previously excluded sequences is now above 80%. These are the ones you need to check carefully: do you agree that they should be included? If there is any doubt, perhaps because of a really marginal E-value, poor coverage or a function annotation that is not compatible with your query, it is safer to exclude a sequence than to risk profile corruption. If the sequence is a true positive, it will return to the list in later iterations, usually with a better E-value as the profile improves. It's a good idea to note such sequences in your journal so you can keep track of how their E-values change.
# Let's exclude partial matches one more time. Again, deselect all sequences with less than 80% coverage. Then run the third iteration.
# Iterate the search in this way, successively relaxing the coverage threshold, until no more "New" sequences are added to the profile. The search has converged. Obviously the result depends on your data, but it would be unusual if the search had not converged after 6 iterations or so, and there is probably a mistake if there are more than 70 hits or so. 
# Now look at the list of excluded hits (if any), the hits that are reasonable but didn't quite make the cut. Are there any from ''Spizellomyces punctatus'' that seem like they should actually be included? Perhaps their E-value is only marginally above the threshold? If that's the case, try returning the E-value threshold to the default 0.005 and see what happens...
}}


Once no "new" sequences have been added, we would always get the same result on additional iterations because there are no more changes to the profile. We say that the search has '''converged'''. Time to harvest.


{{task|1=
# In the header section of the BLAST report, click on '''Taxonomy reports''' and find ''Spizellomyces punctatus'' in the '''Organism Report''' section. These are '''your APSES domain homologs'''. All of them. There is a link to the alignment, the BLAST score, the E-value, and a link to the entry in RefSeq.
# From the report copy the sequence identifiers from ''Spizellomyces punctatus'', with E-values above your defined threshold to your notebook.
}}

For example, the list of ''Saccharomyces'' genes contains the following information:

<code>
[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 \>
4e-37   [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]<br \>
2e-30   [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]<br \>
4e-27   [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]<br \>
4e-27   [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]<br \>
7e-06   [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]<br \>
</code>


<small>[[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.</small>


{{task|1=

# To add the sequences to your database, open each of the links to the RefSeq record for ''Spizellomyces punctatus'' organism into a separate tab.
# Find the UniProt IDs
# Go through the (short) section <code>add PSI BLAST results</code> in the Assignment 04 R-script.

}}

<!--

 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

-->


So much for using PSI-BLAST. The last step seems a bit tedious, adding all this information by hand. There's got to be a better way, right?

But for now, we'll have a look at what the sequences tell us.

{{Vspace}}

{{Vspace}}

==Model Based Alignments: PSSMs and HMMs==

{{Vspace}}

<div class="colmask doublepage">
  <div class="colleft">
    <div class="col1">
      <!-- Column 1 start -->

;Position Specific Scoring Matrices (PSSMs)

The sensitivity of PSI-BLAST is based on the alignment of profiles of related sequences. The profiles are represented as position specific scoring matrices compiled from the alignment of hits, first to the original sequence and then to the profile. Incidentally, this process can also be turned around, and a collection of pre-compiled PSSMs can be used to annotate protein sequence: this is the principle employed by RPS-BLAST, the tool that identifies conserved domains at the beginning of every BLAST search, and has been used to build the CDD database of conserved domains (for a very informative help-page on CDD [https://www.ncbi.nlm.nih.gov/Structure/cdd/cdd_help.shtml '''see here'''].

<!--
=== CDD domain annotation ===

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.


{{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.
}}
-->



      <!-- Column 1 end -->
    </div>
    <div class="col2">
      <!-- Column 2 start -->

; Hidden Markov Models (HMMs)

An approach to represent such profile information that is more general than PSSMs is a {{WP|Hidden Markov model|'''Hidden Markov model (HMM)'''}} and the standard tool to use HMMs in Bioinformatics is [http://hmmer.org/ '''HMMER'''], written by Sean Eddy. HMMER has allowed to represent the entirety of protein sequences as a collection of profiles, stored in databases such as [http://pfam.xfam.org/ '''Pfam'''], [https://www.ebi.ac.uk/interpro/ '''Interpro'''], and [http://smart.embl-heidelberg.de/ '''SMART'''].  While the details are slightly different, all of these services allow to scan sequences for the presence of domains. Importantly thus, the alignment results are not collections of full-length protein families, but annotate to domain families, ''i.e.'' full-length proteins are decomposed into their homologous domains. This is a very powerful approach towards the functional annotation of unknown sequences.

In this section, we will annotate the ''Spizellomyces punctatus'' sequence with the domains it contains, using the database of domain HMMs curated by SMART in Heidelberg and Pfam at the EMBL. We will then compare these annotations with those determined for the orthologues in the reference species. In this way we can enhance the information about one protein by determining how its features are conserved.

      <!-- Column 2 end -->
    </div>
  </div>
</div>


{{Vspace}}

{{Vspace}}

==Suboptimal alignments==

Many sequences contain internal duplications - our's do too, in the ankyrin domains. To detect these, programs exist that will detect suboptimal alignments. 

{{task|1=

* Retrieve the sequence of <code>MBP1_SACCE</code> from <code>myDB</code>.
* Access the [http://www.ebi.ac.uk/Tools/pfa/radar/ '''RADAR''' server at the EBI] and paste the sequence into the form.
* Review the results. Keep the window open, to compare the regions flagged as self-similar by RADAR with the SMART annotations you preduce below.

}}


{{Vspace}}

{{Vspace}}

== SMART domain annotation ==


The [http://smart.embl-heidelberg.de/ SMART database] at the EMBL in Heidelberg integrates a number of feature detection tools including Pfam domain annotation and its own, HMM based SMART domain database. You can search by sequence, or by accession number and retrieve domain annotations and more.


===SMART search===

{{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 <code>MBP1_SPIPU</code> UniProtID into the '''Sequence ID or ACC''' field. Note: you will need to use the <code>chunkSeq()</code> utility function to get its entire length printed to the console so you can copy/paste it - RStudio will truncate long sequences.
# Check all the boxes for:
## '''outlier homologues''' (also including homologues in the PDB structure database)
## '''PFAM domains''' (domains defined by sequence similarity in the PFAM database)
## '''signal peptides''' (using 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)
# 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 ...", try again with the actual sequence instead of the accession number.)</small>

Study the results. 

# Note down the following information so you can enter the annotation in the protein database for ''Spizellomyces punctatus'':
## From the section on "Confidently predicted domains ..."
### The start and end coordinates of the '''KilA-N''' domain <small>(...according to SMART, not Pfam, in case the two differ)</small>.
### All start and end coordinates of '''low complexity segments'''
### All start and end coordinates of '''ANK''' (Ankyrin) domains. 
### Start and end coordinates of '''coiled coil''' domain(s).
### Start and end coordinates of '''AT hook''' domain(s) <small>I expect at most one - not all Mbp1 orthologues have one.</small>
## From the section on "Features NOT shown ..."
### Start and end coordinates of '''low complexity segments''' for which the ''Reason'' is "overlap".
### Start and end coordinates of overlapping '''Ankyrin''' (ANK_...) domains. Don't enter them all individually, but consolidate th contiguous/overlapping domains into a single annotation.
### Start and end coordinates of overlapping '''coiled coil''' segments.
### <small>If there are other features I have not mentioned here, feel encouraged to let me know.</small>
## From the section on "Outlier homologues ..."
### Start and end coordinates of a '''PDB:1SW6{{!}}B''' annotation (if you have one): this is a region of sequence similarity to a protein for which the 3D structural coordinate are known.
### <small>Of course there should also be annotations to the structure of 1BM8 / 1MB1 and/or 1L3G - all of which are structures of the Mbp1 APSES domain that we have already annotated as an "APSES fold" feature previously. And there may be BLAST annotations to Ankyrin domains. We will not annotate these separately either.</small>
# Follow the links to the database entries for the annotated domains so you know what these domains and features are and understand their biological implications.

Finally,  we'll enter the features into our database, so we can compare them with the annotations that I have prepared from SMART annotations of Mbp1 orthologues from the ten reference fungi.

* Complete <code>PART SIX: SMART DOMAIN ANNOTATIONS</code> of the <code>Sequence.R</code> script.

}}

{{Vspace}}

After you execute the code, your plot should look similar to this one but include the <code>MBP1_SPIPU</code> annotations:

[[Image:DomainAnnotations.jpg|frame|none|SMART domain annotations for Mbp1 proteins for the ten reference fungi.
]]

A note on the '''R''' code up to this point: You will find that we have been writing a lot of nested expressions for selections that join data from multiple tables of our data model. When I teach '''R''' workshops for graduate students, postdocs and research fellows, I find that the single greatest barrier in their actual research work is the preparation of data for analysis: filtering, selecting, cross-referencing, and integrating data from different sources. By now, I hope you will have acquired a somewhat robust sense for achieving this. You can imagine that there are ways to simplify those tasks with functions you write, or special resources from a variety of different packages you cab install. But the "pedestrian" approach we have been taking in our scripts has the advantage of working from a very small number of principles, with very few syntactic elements.

{{Vspace}}

{{Vspace}}

==Multiple Sequence Alignment==

{{Vspace}}

In order to perform a multiple sequence alignment, we obviously need a set of homologous sequences. This is not trivial. All interpretation of MSA results depends '''absolutely''' on how the input sequences were chosen. Should we include only orthologues, or paralogues 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: 
*orthologues are expected to be functionally and structurally conserved; 
*paralogues 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.

&nbsp;<br>


{{Vspace}}

===Computing an MSA in R===

{{Vspace}}


We shall use the Bioconductor '''msa''' package to align the sequences we have. Study and run the following code


{{task|1 =

* Return to your RStudio session.
* Make sure you have saved <code>myDB</code> as instructed previously. 
* Bring code and data resources up to date by pulling the most recent version from GitHub.
* Type <code>init()</code> to re-load any updated, supporting scripts.
* Study and work through the code in the <code>PART SEVEN: MULTIPLE SEQUENCE ALIGNMENT</code> section of the <code>Sequence.R</code> script.
* Note that the final task asks you to print out some results, I may ask you to hand these in for credit at a later point.

}}

{{Vspace}}

===Sequence alignment editors===

{{Vspace}}

Really excellent software tools have been written that help you visualize and manually curate multiple sequence alignments. If anything, I think they tend to do too much. Past versions of the course have used Jalview, but I have heard good things of AliView <small>(and if you are on a Mac [https://github.com/4ment/seqotron seqotron] might interest you, but I only cover software that is free and runs on all three major platforms)</small>.

Right now, I am just mentioning the two alignment editors. If you have experience with comparing them, let us know.

* [[http://www.jalview.org/ '''Jalview''']] an integrated MSA editor and sequence annotation workbench from the Barton lab in Dundee. Lots of functions.
* [[http://www.ormbunkar.se/aliview/ '''AliView''']] from Uppsala: fast, lean, looks to be very practical.


{{Vspace}}

<!--

====Jalview: alignment editor====


Geoff Barton's lab in Dundee has developed an integrated MSA editor and sequence annotation workbench with a number of very useful functions. It is written in Java and should run on Mac, Linux and Windows platforms without modifications. 


{{#pmid: 19151095}}


<!--
We will quickly install Jalview and explore its features in other assignments. 

{{task|1=
#Navigate to the [http://www.jalview.org/ Jalview homepage] click on the '''Download''' link, and install Jalview on your computer. For Mac OS X, use the '''Install Jalview Only''' link.
# Start Jalview. A number of windows that showcase the program's abilities will load, you can close these.
#Select File &rarr; Input Alignment &rarr; from File and open the <code>APSES_proteins.mfa</code> file you have prepared above. An alignment window with sequences should appear.
# Choose '''Web Service''' &rarr; '''Alignment''' &rarr; '''Tcoffee with Defaults''' to run a Tcoffee MSA remotely at the Barton lab. The program should execute remotely and download the aligned results into a new window. Scroll along the window to get a sense of what has and hasn't been aligned.
#Select File &rarr; Input Alignment &rarr; from File and open the <code>APSES_proteins_muscle.mfa</code> file you have prepared above. An alignment window with your Muscle alignment should appear.
#Compare the two alignments and get a sense for how similar or different they are.

}}


&nbsp;
-->

===Computing alignments online===


The EBI has a very convenient [http://www.ebi.ac.uk/Tools/msa/ page to access a number of MSA algorithms]. This is especially convenient when you want to compare, e.g. T-Coffee and Muscle and MAFFT results to see which regions of your alignment are robust. 

{{task|1=
# Open the <code>.mfa</code> file <code>APSES_proteins.mfa</code> that you hve created in the previous RStudio session. 
# Copy the sequences.
# Navigate to the [http://www.ebi.ac.uk/Tools/msa '''EBI MSA tools page'''], continue to the MAFFT page.
# Paste your sequences into the form.
# Click on '''Submit'''.
# In separate windows, do the same for '''Clustal Omega''', '''Kalign''', and '''TCoffee'''.
# Compare and interpret the results. Did any of the algorithms seem to get a difficult region more correct than the others? Which of the annotated regions does this correspond to?
# Post a summary of your comparison / interpretation on the mailing list.
}}

{{Vspace}}

{{Vspace}}

==Sequence annotation==

{{Vspace}}

In this assignment we will first download a number of APSES domain containing sequences into our database - and we will automate the process. Then we will annotate them with domain data. First manually, and then again, we will automate this. Next we will extract the APSES domains from our database according to the annotations. And finally we will align them, and visualize domain conservation in the 3D model to study parts of the protein that are conserved.

{{Vspace}}

===Downloading Protein Data From the Web===

{{task|1=

* Work through <code>PART EIGHT: DOWNLOADING DATA FROM THE WEB</code> of the <code>Sequence.R</code> script.
* Part of this code makes use of the <code>reutils</code> package: open the [https://cran.r-project.org/web/packages/reutils/reutils.pdf package vignette on CRAN] so you can read in more detail about the functions we cover.
}}

{{Vspace}}

===Computing annotations over MSAs===

{{task|1=

* Work through <code>PART NINE: COMPUTING OVER MSAs</code> of the <code>Sequence.R</code> script.

}}


[[Image:InformationPlot.jpg|frame|none|Plot of information vs. sequence position produced by the '''R''' script above, for an alignment of Mbp1 ortholog APSES domains.]]

{{Vspace}}

{{Vspace}}

==Genomics==

{{Vspace}}


Large scale genome sequencing and annotation has made a wealth of information available that is all related to the same biological objects: the DNA. The information however can be of very different types, it includes:
* the actual sequence
* sequence variants (SNPs and CNVs)
* conservation between related species
* genes (with introns and exons)
* mRNAs
* expression levels
* regulatory features such as transcription factor bindings sites
and much more.

Since all of this information relates to specific positions or ranges on the chromosome, displaying it alongside the chromosomal coordinates is a useful way to integrate and visualize it. We call such strips of annotation ''tracts'' and display them in ''genome browsers''. Quite a number of such browsers exist and most work on the same principle: server hosted databases are queried through a Web interface; the resulting data is displayed graphically in a Web browser window. The large data centres each have their own browsers, but arguably the best engineered, most informative and mostly widely used one is provided by the University of California Santa Cruz (UCSC) Genome Browser Project. 

Compiling the data requires a massive annotation effort, which has not been completed for all genome-sequenced species. In particular, ''Spizellomyces punctatus'' has not been included in the major model-organism annotation efforts. The general strategy for analysis of a gene in ''Spizellomyces punctatus'' is thus to map it to homologous genes in {{WP|Model organism|model organisms}}. In this assignment you will explore the UCSC genome browser and we will go through an exercise that relates fungal replication genes to human genes. We have previously focused a lot on Mbp1 homologs, but these have no clear equivalences in "higher" eukaryotes. However one of the key target genes of Mbp1 is the cell cycle protein {{WP|Cdc6}}, which is well conserved in fungi and other eukaryotes eukaryotes and has a {{WP|CDC6|human homolog}}. Since generally speaking the annotation level for human genes is the highest, we will have a closer look at that gene.

{{Vspace}}

===The UCSC genome browser===
{{smallvspace}}

The University of California Santa Cruz (UCSC) Genome Browser Project has the largest offering of annotation information. However it is strictly model-organism oriented and you will probably not find YFO among its curated genomes. Nevertheless, if you are studying eg. human genes, or yeast, the UCSC browser will probably be your first choice.

{{task|1=

In this task you will access the UCSC genome browser view of the <!-- yeast Cdc6 gene and its human orthologue --> human Cdc6 gene. You will explore some of the very large number of tracks that are available and study the transcription factor binding region.

# Navigate to the [http://genome.ucsc.edu/ '''UCSC''' Genome Bioinformatics entry page] and follow the link to the '''Genome Browser''' in the "Our tools" section.
<!--
# From the available menus, access the ''S. cerevisiae'' information ('''group &rarr; other''') and enter Cdc6 as the '''search term'''.
# Click on the link to the [http://genome.ucsc.edu/cgi-bin/hgTracks?position=chrX:69338-70879&hgsid=311433759&sgdGene=pack&hgFind.matches=YJL194W, Cdc6 gene] on chromosome X.
# Click on the button to zoom out '''3x''' - we want to see the upstream regulatory region.
# In the subsection for '''Expression and Regulation''', find the menu for '''Regulatory Code''' and select '''full'''; select '''hide''' for all other expression tracks. Click '''refresh'''.

Up to now, this looks very similar to the SGD genome browser.

# Open a second window, and access the UCSC Genome browser for the '''human genome'''. Search for CDC6 and click the link to [http://genome.ucsc.edu/cgi-bin/hgTracks?position=chr17:38444146-38459413&hgsid=394751891_WVshDjZBOw5nRfbXOotacA9pGJn5&knownGene=pack&hgFind.matches=uc002huj.1, <code>Homo sapiens cell division cycle 6 (CDC6), mRNA</code>] on chromosome 17.
-->
# Click on the link to humans. Note that this is the hg38 assembly.
# Enter CDC6 into the "Position/Search Term" field and click "Go". You should get a list of entries, click on the top link, the gene on chromosome 17: <tt>[http://genome.ucsc.edu/cgi-bin/hgTracks?position=chr17:40287633-40304657&hgsid=570479629_xD9YY3QMJ4u2xrTagkgV7xJMqEen&knownGene=pack&hgFind.matches=uc002huj.2, CDC6 (uc002huj.2) at chr17:40287633-40304657]</tt>

# Zoom out '''1.5x''' to view the upstream regulatory region: the end of the adjacent WIPF2 gene should have just come into view on the left.
# Study the Genome Browser view of the human CDC6 homolog.
## In particular, note the extensive functional annotations of DNA and the alignments of vertebrate syntenic regions that allow detailed genomic comparisons.
## Distinguish between exon and intron sequence.
## Note that the mammal Conservation track has high values for all of the exons, but not only for exons.
## Find more information on the "Layered H3K27Ac" tract.

# Note the '''large''' number of available tracks that have been integrated into this view. Most of them are switched off. Find the '''Regulation''' section, and follow the link to the "ORegAnno" information to see what that is about. Note that you can switch individual annotations on or off on this page, as well as set the display format for all of the results. Select the check-box '''only''' for "transcription factor binding site" to be on, select the "Display mode" to '''full''' and click '''submit'''.
# Study this information and note:
## There is a cluster of TFBS just upstream of the transcription initiation site.
## This cluster coincides with the highest H3K27Ac density.
## If you &lt;control&gt;-click (right-click?) on the top orange bar of this cluster, a contextual menu opens from which you can access the details page for OREG1791811 in a new window. Follow the link to the RBL2 transcription factor via [http://useast.ensembl.org/Homo_sapiens/Transcript/Summary?g=ENSG00000103479;r=16:53445781-53491648;t=ENST00000379935 ENST00000379935] ... from where you can access transcript and gene and expression and protein family and GO and all other information.
# Go back to the Genome Browser and set the ORegAnno tract to "pack" and click "refresh".
# Slide the SNP track to just beneath the RefSeq genes track that contains the introns and exons. You will notice that one of the SNPs is green, and two are red. Why? Set the "Common SNPs" track display mode to "pack" and click "refresh".
}}


Based on this kind of information, it should be straightforward to identify human transcription factors that potentially regulate human Cdc6 and determine - via sequence comparisons - whether any of them are homologous to any of the yeast transcription factors or factors in YFO. Through a detailed analysis of existing systems, their regulatory components and the conservation of regulation, one can in principle establish functional equivalences across large evolutionary distances.

<!--
The UCSC browser has a sometimes bewildering amount of information available. But its curators are aware of the need for educating users regarding the utility of their tools.

{{task|1=

In this task you will access some of the tutorial information that UCSC provides.
# Return to the [http://genome.ucsc.edu/ '''UCSC''' Genome Bioinformatics entry page] and follow the link to '''Training''' in the left-hand menu.
# Follow the link to the [http://www.openhelix.com/ucsc '''OpenHelix UCSC tutorials'''].
# Download the Hands-on exercise PDF file and work through '''Exercise 2''' (the rat leptin exercise).
}}

This exercise includes a number of interesting options to work with the UCSC data - the BLAT tool for genomic region alignment and the selective display of SNP annotations.

; Optional
* Work through exercise one and three of the OpenHelix UCSC introduction.
* Access the [http://www.openhelix.com/ENCODE2 OpenHelix '''ENCODE''' tutorial], download the '''Hands-on Exercises''' pdf and work through the exercises. Exercise 3 is particularly valuable, as it teaches you how to create results from complex intersections of queries.
* You can also work through the [http://www.nature.com/scitable/ebooks/guide-to-the-ucsc-genome-browser-16569863 Guide to the UCSC Genome Browser at "nature"] which gives an excellent, in-depth overview.
* Study the ''User's guide to ENCODE'' paper linked below.

{{task|1=

Finally: 
# Print this page, but print the first page only.
# With a red pen, mark and label the following four items on your print-out:
## The first exon of CDC6.
## The chromosomal coordinates of the current view.
## The binding sites for the transcription factors that bind to the CDC6 promoter.
## The locations of the missense-variant SNPs. 
# Write your name and Student number on this page and bring it to class to hand it in on Tuesday.
}}

-->


&nbsp;

== Links and resources ==
* [[Regular Expressions|More on Regular Expressions]]

{{#pmid: 22764121}}
{{smallvspace}}
{{#pmid: 26527727}}
{{#pmid: 25762420}}
{{#pmid: 21526222}}
{{smallvspace}}
{{#pmid: 25645873}}

*{{#pmid: 15286655}}

*{{PDFlink|[https://www.bioconductor.org/packages/release/bioc/vignettes/Biostrings/inst/doc/BiostringsQuickOverview.pdf  Biostrings Quick Overview]| summary of Biostrings functions (PDF)}}

*[[Reference_species_for_fungi|Our 10 "reference" fungi]]



{{#lst:Bioinformatics_Introduction_Data|intro_footer}}

{{Vspace}}


<table style="width:100%;"><tr>
<td style="height:30px; align:center; vertical-align:middle;">[[Bioinformatics_Introduction_Data|Data]]</td>
<td style="height:30px; align:center; vertical-align:middle;">''Sequence''</td>
<td style="height:30px; align:center; vertical-align:middle;">[[Bioinformatics_Introduction_Structure|Structure]]</td>
<td style="height:30px; align:center; vertical-align:middle;">[[Bioinformatics_Introduction_Phylogeny|Phylogeny]]</td>
<td style="height:30px; align:center; vertical-align:middle;">[[Bioinformatics_Introduction_Function|Function]]</td>
</tr></table>


[[Category:Bioinformatics]]
</div>