Search This Blog

Monday, July 26, 2010

Wrong about RALEE bug

After having email discussion with Sam Griffiths-Jones, the author of RALEE, I discovered that there was a difference interpretation, not a bug, and I decided that Sam's interpretation was better.

I had interpreted a contiguous run of base pairs to be a helix, but after looking at a few structures in VARNA, I realized that some of the "helices" could be considered one helix with a few bulges. Sam gave some good examples:

You have continuous on one side, but not the other:

.((((((...))).))).

Alignments where a helix is continuous in most sequences, but bulged occasionally:

AGCGCUUCG.CU
AGCGCUUCG.UU
AGCGCCGCG.CU
AGCGCUUCG.CU
AGGGCUUCC.CU
AGCGCUUCGGCU
((((...)).))

In the above case, what is a helix could be subject to arbitrary change, just by adding one sequence with a bulge.

Tuesday, July 13, 2010

Bug in RALEE code

I came across across a bug in the RALEE code that processes the secondary structure line in Stockholm files. The Stockholm file (of vault proteins) just happened to break the code that processes how many helices there are and associates a helix number with each base pair position. Here are the cases that I had to fix:

1. <<..<<..<<..<<<..>>>>>..>>>>

2. <<<<..<<<<<..>>>..>>..>>..>>

If there isn't another separate helix within the sequence, the helices that aren't separated by unpaired nucleotides can't be detected in the original code. I was able to fix the bug by adding more cases. Hopefully I haven't created another bug.

Monday, July 5, 2010

What happens when you open a Stockholm File?

Alot of things, and a journey through how Jalview works. I wrote out part of this before in my notebook, but I hadn't connected it to the Stockholm file parser before. I haven't described the different types of objects yet, but here is the data flow:


jalview.gui.Desktop calls file loading functions, Stockholm parse is called, and then a new window with the alignment is displayed.

In jalview.gui package
-Desktop.instance.inputLocalFileMenuItem_actionPerformed(viewport) is called
--->calls FileLoader.LoadFile() in jalview.io package
-------this starts a new thread, which calls FileLoader.run() (use start() method to do this)
-----------Checks file format
-----------Alignment object = FormatAdapter().readFromFile(source,format) in jalview.io package
-------------->calls AppletFormatAdapter.readFromFile
------------------>Calls StockholmFile()

In StockholmFile in jalview.io
-StockholmFile extends AlignFile which extends FileParse
--->Fileparse() resolves type, checks if valid file type, position in file
--->AlignFile() calls parse(), which is in StockholmFile()
------->StockholmFile.parse() called
-----------checks if file has Stockholm file features, e.g. starts with #STOCKHOLM \d.\d
-----------use regular expressions to figure out line type
--------------if not an annotation line, process sequence line and store in hashtable
--------------else, figure out annotation type
------------------>if GC annotation, call parseAnnotationRow() to parse secondary structure

StockholmFile.parseAnnotationRow(annotations member of AlignFile, ID, sequence)
-make Annotation array else
-for each position in sequence
---->make Annotation obj ann
-------->component secondary structure determined by getDssp3State (jalview.schemes.ResidueProperties)

(Basically this is creating an annotation for each position)

-create AlignmentAnnotation object annots
--->calls validateRangeAndDisplay()
------>calls labelsSecondaryStructure()

-returns annots


Finally, back in the Desktop instance, the viewport is changed and you see the alignment visualization.
--if viewport exists, update using firePropertyChange()
--else new AlignFrame object is created
-----new AlignViewport object is created and initialized

Regular expressions in Jalview

I was looking at this about a week ago, and I decided that it might be worthy of a blog post. Part of the parsing of Stockholm files in Jalview utilizes the power of regular expressions, and Jalview uses the com.stevesoft.pat package. It simplifies some of the syntax; a tutorial is online here.

With this package you can create Regex objects and compile and match regular expressions. The Stockholm file parser detects the end of the file, annotation and sequence lines, and parses lines using regular expressions.

Saturday, July 3, 2010

Bug tracker, svn

Jalview has a bug tracker which can be found here.

I found a bug with reading Stockholm files and I logged it there. There were a few Stockholm files from Rfam that I wasn't able to open. "\\" is used to denote the end of a Stockholm file, which coincidentally is found in urls! I fixed the bug by requiring that the line with "//" start with "//" or a space then "//".

After fixing this bug I decided to start a branch in my svn repository (just to try it out), but I found out that I did not do this correctly. The svn red book has a nice entry about branch philosophy. What you should do is copy the existing trunk into the new branch. This allows diffs to be created properly between different versions. For now I am not going to have branches since Jim is going to set up an official repository for Jalview later.

Thursday, June 3, 2010

How RALEE parses the secondary structure line

I've figured out how RALEE interprets the secondary structure line and figures out helices…I really don't like Lisp, haha. The GNU Emacs Lisp Reference Manual isn't too bad: http://www.gnu.org/software/emacs/manual/html_mono/elisp.html

Basically, this is what happens, in two separate functions.

Parse the secondary structure line from the Stockholm file (ralee-get-base-pairs)

Description: Keep a list of the positions of the open brackets and a list whose elements are paired positions. When you reach a close bracket, pair it with the last open bracket and store this pair in the "pairs" list.

  • Split the line into a "position" list (each character of the line is an element of the list)

  • Keep a "stack" list and a "pairs" list

  • Go through each element of the "position" list
    *if the base is an open bracket type, add the current position as the first element of the "stack" list

    *if the base is a close bracket type, create a 2 element list that consists of the first element of "stack" list and the current position. Add this list as the first element of the "pairs" list and remove the first element from the stack list.

Numbers indicate the order in which the pairs are added to the "pairs" list

<<..<<..>>..<<..>>>>
65..21..12..43..3456

Create a Hash of positions and which helix they belong to (ralee-helix-map)
Description: This function takes the pairs list generated by the above function as its input. It goes through the pairs and figures out how many helices there are and creates a hash where the keys are positions and the values are which helix they belong to.

Variables:
helix = the current helix, an int
helices = hash of the positions and which helix they belong to
lastclose = position of the last close bracket reviewed
lastopen = position of the last open bracket reviewed
open = position of current open bracket
close = position of current close bracket
pair = current pair (consists of open and close)
pairs = list of base pairs (input, elements are 2 element lists)
i = current pair (in pairs list)

  • Go through each item in pairs list

    *for the current pair, open is the first element and close is the second element
    *Check for an inner helix:
     catch things like
    ; <<..>>..<<..>>
    ; *
    if the lastclose comes before open, increment the current helix
    *Check for bulges:
     catch things like
    ; <<..<<..>>..<<..>>>>
    ; *

    <<..<<..>>..<<..>>>>
    cc..aa..aa..bb..bbcc

    compare current pair to all other pairs (this works because of the way that the pairs are stored)
    **if open of a pair comes before lastopen and also after current open
    **then find which helix the open belongs to. If it belongs to the current helix, do nothing. Otherwise increment the number of helices.

  • add the current open and close to the helices hash with the helix as the value

  • set lastopen and lastclose to open and close

Wednesday, June 2, 2010

Parsing WUSS notation of RNA secondary structure annotation

A key part of this project is to parse the secondary structure line of Stockholm files so that it can be interpreted for coloring schemes. I have been adding mini-goals as appropriate. I will probably also need to add code to check that the sequence length and secondary structure length are the same, as well as the same number of open and closed parentheses.

WUSS notation is used in RNA stockholm files to indicate secondary structure. WUSS notation can support more characters than I thought, but Rfam uses the simplified version that the covariance modeling program Infernal uses. The description of Rfam on the Janelia Farm page is


Rfam is a collection of multiple sequence alignments and covariance models covering many common non-coding RNA families. The main use of Rfam is as a source of RNA multiple alignments with consensus secondary structure annotation in a consistent format. In conjunction with the Infernal software package, Rfam covariance models (CMs) can be used to search genomes or other DNA sequence databases for homologs to known structural RNA families.


WUSS notation uses <>, (), [], and {} to indicate base pairs and ':', ',', '_', '.', and '~' as single stranded columns. Each type of symbols has subtle meaning, but for Infernal the structure annotation line only needs to indicate which columns are base paired to each other. Thus, full WUSS notation is not necessary and a simple minimal annotation uses <> to indicate base pairs and '.' for single stranded positions of the alignment.

In more detail taken from the Infernal user guide:

Base pairs: the different symbols indicate different depth
*<> for simple terminal stems
*() for "internal" helices enclosing a multifunction of all terminal stems
*[] for internal helices enclosing a multifunction that includes at least one annotated () stem already
*{} for all internal helices enclosing deeper multifurcations

Hairpin loops
*indicated by underscores '_'
*Simple stem loops example: <<<____>>>

Bulge, interior loops
*indicated by dashes '-'

Multifurcation loops
*indicated by commas ','
*example: <<<___>>>,,<<<__>>>

External residues completely outside structure
*indicated by colons ':'

Insertions
* . to a known structure
* ~ used to indicate that a local structural alignment left regions of target and query unaligned.

Pseudoknots
* pairs of upper case/lower case letters
* example: <<<<_AAAA____>>>>aaaa


Things that I am thinking about:

-I need to interpret WUSS notation in a general way. It shouldn't be too difficult, but it is necessary since the same structure can be written in multiple ways. An example from the Infernal user guide is : <<<<....>>>> and ((((____)))) and <(<(._._)>)> all indicate a four base stem with a four base loop

-How should I store the secondary structure line so that it will be easily interpreted to implement coloring schemes?

Potentially I can store pairs of positions like how the disulfide bond positions are stored as annotations (Jim pointed this one out). I also need to keep in mind that bulges might exist, so I can't just interpret a run of the same type of bracket as part of the same stem. VARNA interprets bulges just fine, so I don't have to worry about that. An example of a complicated structure with a bulge:

<<<<……<<<< <<<<…..>>>>..>>>>……<<<<…>>>>….>>>>

-How can I make sure that there are 4 stems instead of 3? I can't simply scan through from left to right or eat away at both ends at the same time. It looks like the RALEE mode for Emacs handles bulges just fine based on this example in the readme




0123456789012345678901234
.<<<<<...>>.<<...>>..>>>.


Column 1 pairs with 23
2 with 22
3 with 21
4 with 10
5 with 9
12 with 18
13 with 17



The image is from VARNA. Note the numbering in the image starts at 1 instead of 0.

RALEE is written in Emacs Lisp, so I need to look up some basics in Lisp before I can feel confident that I'm interpreting the code correctly! I think that this code will cut down on my thinking time, however.

To do
-check that secondary structure line and sequence are the same length. Does Jalview already do this?
-change all bracket types to () for VARNA (I just noticed that VARNA only likes (), not <> for base pairing! )
-convert all WUSS symbols to something simple, like how Jalview already does for protein secondary structure (simple helices and sheets)
-Need to figure out how to detect pseudoknots
-Add support for error checking when a user adds a base pair annotation. Make sure same number of column groups are selected
-How will colors cycle for different numbers of stems?