BIN-ALI-Optimal sequence alignment

From "A B C"
Jump to navigation Jump to search

Optimal global and local sequence alignment

(NWS (optimal global) and SW (optimal local) algorithms, alignment via EMBOSS tools in practice, interpretation of alignments)


 


Abstract:

This unit covers the concepts and algorithms for optimal pairwise sequence alignments.


Objectives:
This unit will ...

  • ... discuss how homology is inferred from optimal sequence alignments, by using scoring matrices that represent an evolutionary relationship;
  • ... introduce the principle of dynamic programming alignment works by optimizing the sum of (context independent) pairscores, using an affine gap model for indels, and backtracking to reconstruct an alignment from contributing cells in the path-matrix;
  • ... point out problems associated with affine gap functions and how parameter choice influences size and distribution of indels;
  • ... teach the difference between global and local optimal alignment and in which situation these algorithms are appropriately used;
  • ... demonstrate how to calculate optimal sequence alignments with online EMBOSS tools, and in R code with the Biostrings package.;

Outcomes:
After working through this unit you ...

  • ... can produce and interpret optimal sequence alignments, online, and in R code.

Deliverables:

  • Time management: Before you begin, estimate how long it will take you to complete this unit. Then, record in your course journal: the number of hours you estimated, the number of hours you worked on the unit, and the amount of time that passed between start and completion of this unit.
  • Journal: Document your progress in your Course Journal. Some tasks may ask you to include specific items in your journal. Don't overlook these.
  • Insights: If you find something particularly noteworthy about this unit, make a note in your insights! page.
  • Your protein database
    Add APSES domain annotations for MBP1_MYSPE proteins to your database.

    Prerequisites:
    This unit builds on material covered in the following prerequisite units:


     



     



     


    Evaluation

    Evaluation: NA

    This unit is not evaluated for course marks.

    Contents

    Pairwise Alignments: Optimal

     

    Task:


     

    Optimal pairwise sequence alignment is the mainstay of sequence comparison. To try our first alignments in practice, we will start with aligning Mbp1 and its MYSPE relative. For simplicity, I will call the two proteins MBP1_SACCE and MBP1_MYSPE through the remainder of the unit.


     

    Optimal Sequence Alignment: EMBOSS online tools

     

    EMBOSS tools are a collection of standard sequence analysis programs. The most important ones are hosted at the EBI, but the EMBOSS explorer site hosts many more. They offer Needlman-Wunsch and Smith-Waterman alignments.


    Task:

    • Fetch the sequences for MBP1_SACCE and MBP1_MYSPE from your database that you have prepared in the BIN-Storing_data unit. Open the RStudio project and enter the code below - substituting the proper name for MYSPE where appropriate.
    source("makeProteinDB.R")
    
    # Print the MBP1_SACCE sequence
    sel <- myDB$protein$name == "MBP1_SACCE"
    myDB$protein$sequence[sel]
    
    # Print the MBP1_MYSPE sequence
    sel <- myDB$protein$name == paste0("MBP1_", biCode(MYSPE))
    myDB$protein$RefSeqID[sel]
    
    

    (If this didn't work, fix the problem. Did you give your sequence the right name in your database?)

    1. Access the EMBOSS tools page at the EBI.
    2. Look for Water, click on protein, paste your sequences and run the program with default parameters.
    3. Study the results. You will probably find that the alignment extends over most of the protein, but does not include the termini.
    4. Considering the sequence identity cutoff we discussed in class (25% over the length of a domain), do you believe that the N-terminal domains (the APSES domains) are homologous?
    5. Change the Gap opening and Gap extension parameters to high values (e.g. 25 and 5). Then run the alignment again.
    6. Note what is different.


    Global optimal sequence alignment using "needle"

    Task:

    1. Look for Needle, click on protein, paste the MBP1_SACCE and MBP1_MYSPE sequences again and run the program with default parameters.
    2. Study the results. You will find that the alignment extends over the entire protein, likely with significant indels at the termini.


     


    Optimal Sequence Alignment with R: Biostrings

     

    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[1].


     

    Task:

     
    • Open RStudio and load the ABC-units R project. If you have loaded it before, choose FileRecent projectsABC-Units. If you have not loaded it before, follow the instructions in the RPR-Introduction unit.
    • Choose ToolsVersion ControlPull Branches to fetch the most recent version of the project from its GitHub repository with all changes and bug fixes included.
    • Type init() if requested.
    • Open the file BIN-ALI-Optimal_sequence_alignment.R and follow the instructions.


     

    Note: take care that you understand all of the code in the script. Evaluation in this course is cumulative and you may be asked to explain any part of code.


     


     

    Further reading, links and resources

    Fitch (2000) Homology a personal view on some of the problems. Trends Genet 16:227-31. (pmid: 10782117)

    PubMed ] [ DOI ] There are many problems relating to defining the terminology used to describe various biological relationships and getting agreement on which definitions are best. Here, I examine 15 terminological problems, all of which are current, and all of which relate to the usage of homology and its associated terms. I suggest a set of definitions that are intended to be totally consistent among themselves and also as consistent as possible with most current usage.

    Notes

    1. 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.


     


    About ...
     
    Author:

    Boris Steipe <boris.steipe@utoronto.ca>

    Created:

    2017-08-05

    Modified:

    2020-09-24

    Version:

    1.1

    Version history:

    • 1.1 2020 Updates
    • 1.0 First live
    • 0.1 First stub

    CreativeCommonsBy.png This copyrighted material is licensed under a Creative Commons Attribution 4.0 International License. Follow the link to learn more.