Ed's Big Plans

Computing for Science and Awesome

C# & Bioinformatics: Edit Strings and Local Alignment Subprofiles

with 2 comments

Alternatively: The awful finite string manipulation bug that took me forever to fix.

Update: Code last update 2010 Dec. 12 (branch logic fixed: gave default “else” case highest priority instead).

Every so often, a bug creeps in during the implementation of a utility function that halts an entire project until it can be fixed. I bumped into one that was so insidious — so retrospectively obvious — that I must share it with you. This particular bug has to do with the way Edit Strings behave when you use them to represent changes done during a Local Alignment rather than a Global Alignment.

Those of you who have looked into MUSCLE may have seen the 2004 BMC paper that goes into great detail about its implementation. In it, Dr. Edgar describes Edit Strings (e-strings) — these are integer strings which keep track of the edits performed on a sequence in order to insert gaps in the correct locations corresponding to an alignment.

I’ve previously discussed and explained e-strings and their use in sequence alignment here and also how to convert alignment strings (a recording of site-to-site insertions/deletions/substitutions) to edit strings here.

E-Strings in Local Alignments

E-strings weren’t originally designed with local alignments in mind — but that’s OK. We can make a few adjustments that will let us keep track of gaps the way we intend them. In a global alignment, all tokens (characters) in each of the aligned sequences are considered. In a local alignment, we need to be able to specify two aligned subsequences rather than the alignment of two full sequences. We can add a simple integer offset to indicate the number of positions since the beginning of aligned sequences to specify the start position of each subsequence.

These are two sequences s and t — (to keep things easy to visualize, I’ve just used characters from the English alphabet).

s = DEIJKLMNOPQRSTUVWXYZ
t = ABCDEFGHIJKLMNOUVWXY

Assume that the following is the optimal global alignment of s and t.

s(global) = ---DE---IJKLMNOPQRSTUVWXYZ
t(global) = ABCDEFGHIJKLMNO-----UVWXY-

The corresponding edit strings for s and t are es and et respectively.

e(s, global) = < -3, 2, -3, 18 >
e(t, global) = < 15, -5, 5, -1 >

Remember, we can take the sum of the absolute values in two edit strings corresponding to the same alignment to be sure that they contain the number of sites. This is true because a site in a given sequence or profile must be either an amino acid (in a sequence) or a distribution of amino acids (in a profile) or a gap (in either). In this case, the absolute values of both es and et correctly add up to 26.

Now let’s adapt the above to local alignments.

The corresponding optimal local alignment for s and t is as follows.

s(local) = DE---IJKLMNOPQRSTUVWXY
t(local) = DEFGHIJKLMNO-----UVWXY

The edit strings only need an offset to specify the correct start position in each sequence. Remember that our indexing is zero based, so the position for the character “A” was at position zero.

e(s, local) = 3 + < 2, -3, 17 >
e(t, local) = 3 + < 12, -5, 5 >

It’s a pretty straight forward fix. Notice that the absolute values of the integers inside the vectors now only sum up to 22 (corresponding to the sites stretching from “D” to “Y”).

Local Alignment E-String Subprofiles

The bug I faced wasn’t in the above however. It only occurred when I needed to specify subprofiles of local alignments. Imagine having a local alignment of two sequences in memory. Imagine now that there is a specific run of aligned sites therein that you are interested in and want to extract. That specific run should be specifiable by using the original sequences and a corresponding pair of e-strings and pair of offsets. Ideally, the new e-strings and offsets should be derived from the e-strings and offsets we had gathered from performing the original local alignment. We call this specific run of aligned sites a subprofile.

As you will see, this bug is actually a result of forgetting to count the number of gaps needed that occur after the original local alignment offset and before where the new subprofile should start. Don’t worry, this will become clearer with an example.

For the local alignment of s and t, I want to specify a subprofile stretching from the site at letter “I” to the site at letter “V”.

s(local) = DE---IJKLMNOPQRSTUVWXY → s(subprofile) = IJKLMNOPQRSTUV
t(local) = DEFGHIJKLMNO-----UVWXY → t(subprofile) = IJKLMNO-----UV

What do you suppose the edit string to derive would be?

Let’s puzzle this out — we want to take a substring from each s and t from position 5 to position 18 of the local alignment as shown below …

Position:       0    0            1  2
                0    5            8  1
     s(local) = DE---IJKLMNOPQRSTUVWXY →
s(subprofile) =      IJKLMNOPQRSTUV
                |    |            |  |
     t(local) = DEFGHIJKLMNO-----UVWXY →
t(subprofile) =      IJKLMNO-----UV
                     |            |
                     |<-  Want  ->|

Let’s break this down further and get first the vector part of the e-string and then the offset part of the e-string.

Getting the Vector part of the Edit String

Remember, we’re trying to get an edit string that will operate on the original sequences s, t and return ssubprofile, tsubprofile. Getting the vector part of the estring is straight forward, you literally find all of the sites specified in the subsequence and save them. Since we want positions [5, 18] in the local edit strings, we retrieve the parts of the edit strings which correspond to those sites (we’ll ignore offsets for now). To visualize this site-to-site mapping, let’s decompose our local estrings into something I call unit-estrings — an equivalent estring which only uses positive and negative ones (given as “+” and “-” below in the unit vectors us, ut for s, t respectively).

Position:                                0         0                         1     2
                                         0         5                         8     1
    s(local) =                           D E - - - I J K L M N O P Q R S T U V W X Y
e(s, vector) = < 2, -3, 17 > -> u(s) = < + + - - - + + + + + + + + + + + + + + + + + >
e(s, newvec) = < 14 >        <-        <           + + + + + + + + + + + + + +       >
                                                   |<-        want!        ->|
e(t, newvec) = < 7, -5, 2 >  <-        <           + + + + + + + - - - - - + +       >
e(t, vector) = < 12, -5, 5 > -> u(t) = < + + + + + + + + + + + + - - - - - + + + + + >
    t(local) =                           D E F G H I J K L M N O - - - - - U V W X Y

We retrieve as our new vectors es, newvec = < 14 > and et, newvec = < 7, -5, 2 >. Notice that the integer vector elements are just a count of the runs of positive and negative ones in the subprofile region that we want.

Getting the Offset part of the Edit String (the Right way)

This part gets a bit tricky. The offset is a sum of three items.

  1. The offset from the original sequence to the local alignment.
  2. The offset from the local alignment to the subprofile.
  3. The negative of all of the gap characters that are inserted in the local alignment until the start of the subprofile.

The below diagram shows where each of these values come from.

            s =    DEIJKLMNOPQRSTUVWXYZ
     s(local) =    DE---IJKLMNOPQRSTUVWXY -> +0 -- offset from original to local
s(subprofile) =    .....IJKLMNOPQRSTUV    -> +5 -- offset from local to subprofile
                     ---                  -> -3 -- gaps: local start to subprofile start
                                          => +2 == TOTAL
                        |            |
            t = ABCDEFGHIJKLMNOUVWXY
     t(local) = ...DEFGHIJKLMNO-----UVWXY -> +3 -- offset from original to local
t(subprofile) =    .....IJKLMNO-----UV    -> +5 -- offset from local to subprofile
                                          -> -0 -- gaps: local start to subprofile start
                                          => +8 == TOTAL

You might be wondering about the weird gap subtraction for sequence s. This is the solution to my bug! We need this gap subtraction! When you forget to subtract these gaps, you end up with an alignment that is incorrect where the subprofile of s starts too late. But why is this so? When we transform a sequence from the original to the subprofile, the gaps that we’ve taken for granted in the local alignment don’t exist yet. This way, we would be mistakenly pushing our subprofile into a position where non-existent gapped sites are virtually aligned to existing amino acid sites in the alignment partner t. This was a very odd concept for me to grasp, so let yourself think on this for a bit if it still doesn’t make sense. … You could always take it for granted too once I show you below that this is indeed the correct solution 😛

The completed edit strings to specify the subprofile from “I” to “V” for the local alignment of s and t are thus:

e(s, subprofile) = 2 + < 14 >
e(t, subprofile) = 8 + < 7, -5, 2 >

So let’s make sure these work — let’s apply es, subprofile and et, subprofile onto s and t to see if we get back the correct subprofiles.

>>> apply offsets ...

s = DEIJKLMNOPQRSTUVWXYZ : +2 remaining ...
s = IJKLMNOPQRSTUVWXYZ

t = ABCDEFGHIJKLMNOUVWXY : +8 remaining ...
t = IJKLMNOUVWXY

>>> apply e-strings ...

s = IJKLMNOPQRSTUVWXYZ : < 14 > remaining ...
s = IJKLMNOPQRSTUV

t = IJKLMNOUVWXY : < 7, -5, 2 > remaining ...
t = IJKLMNO : UVWXY : < -5, 2 > remaining ...
t = IJKLMNO----- : UVWXY : < 2 > remaining ...
t = IJKLMNO-----UV

s = IJKLMNOPQRSTUV -- selected subprofile ok.
t = IJKLMNO-----UV -- selected subprofile ok.

The above application results in the desired subprofile.

The SubEString Function in C#

The below source code gives the correct implementation in order to find the e-string corresponding to a subprofile with the following arguments. It is an extension function for List<int> — a natural choice since I am treating lists of integer type as edit strings in my project.

  1. this List<int> u — the local alignment e-string from which to create the subprofile e-string.
  2. int start — the offset from the left specifying the location of the subprofile with respect to the local alignment.
  3. int end_before — the indexed site before which we stop; this is different from length since start and end_before are zero-based and are specified from the very start of the local alignment.
  4. out int lostGaps — returns a negative value to add to the offset — this is the number of gaps that occur between the start of the local alignment and the start of the subprofile.

This function returns a new List<int> e-string as you would expect.

public static List<int> SubEString(
     this List<int> u, int start, int end_before, out int lostGaps) {

    if(end_before < start) {
        Console.Error.WriteLine("Error: end_before < start?! WHY?");
        throw new ValueOutOfRange();
    }

    var retObj = new List<int>();
    var consumed = 0; // the amount "consumed" so far.
    var remaining_start = start; // subtract "remaining_start" to zero.
    var ui = 0; // the current index in the e-string
    var uu = 0; // the current element of the e-string
    var uu_abs = 0; // the absolute value of the current element
    lostGaps = 0; // the number of gaps that is subtracted from the offset.

    //START -- dealing with the "start" argument ...
    for(; ui < u.Count; ui ++) {
        uu = u[ui];
        uu_abs = Math.Abs(uu);
        if(uu_abs < remaining_start) {
            if(uu < 0) lostGaps += uu; //LOSTGAPS
            // add uu this time
            // it's the fraction of remaining_start we care about.
            consumed += uu_abs; // uu smaller than start amount, eat uu.
            remaining_start -= uu_abs; // start shrinks by uu.
            uu_abs = 0; // no uu left.
            uu = 0; // no uu left.
            // not done -- go back and get more uu for me to eat.
        }
        else if(consumed + uu_abs == remaining_start) {
            if(uu < 0) lostGaps += -remaining_start; //LOSTGAPS
            // -- if thus uu is equal to the start amount,
            // count them into LOSTGAPS
            // we want to get everything occurring before the subestring start
            consumed += remaining_start; // same size -- consume start amount.
            uu_abs = 0; // MAD (mutually assured destruction 0 -> 0)
            remaining_start = 0; // MAD
            uu = 0; // MAD
            ui ++; break; // done.
        }
        else if(consumed + uu_abs > remaining_start) {
            if(consumed + uu_abs > end_before && consumed < remaining_start) {
                // easier to say ui == u.Count-1?
                if(uu < 0) lostGaps += -remaining_start; //LOSTGAPS
                uu_abs = (end_before - (consumed + remaining_start));
                uu = uu > 0 ? uu_abs : -uu_abs;
                retObj.Add(uu);
                return retObj;
            } else if(consumed + uu_abs < end_before) {
                consumed += uu_abs;
                // Seemingly wrong -- however:
                // push all the way to the end; subtract the beginning of uu_abs
                // ensures correct calculation of remaining string
                uu_abs -= remaining_start; // remove only start amount.
                // we're getting the end half of uu_abs!
                if(uu < 0) lostGaps += -remaining_start; //LOSTGAPS
                remaining_start = 0; // no start left.
                uu = uu > 0 ? uu_abs : -uu_abs; // remaining uu at this ui.
                if(uu != 0) retObj.Add(uu); // fixes a zero-entry bug
                ui ++; break; // done.
            } else {
                if(uu < 0) lostGaps += -remaining_start; //LOSTGAPS
                uu_abs = (end_before - remaining_start);
                uu = uu > 0 ? uu_abs : -uu_abs;
                retObj.Add(uu);
                return retObj;
            }
        }
    }

    //END_BEFORE -- dealing with the "end_before" argument ...
    for(; ui < u.Count; ui ++) {
        uu = u[ui];
        uu_abs = Math.Abs(uu);
        if((consumed + uu_abs) == end_before) {
            retObj.Add(uu);
            break;
        } else if((consumed + uu_abs) < end_before) {
            retObj.Add(uu);
            consumed += uu_abs;
        } else /*if((consumed + uu_abs) > end_before)*/ {
            uu_abs = end_before - consumed;
            uu = uu > 0 ? uu_abs : -uu_abs;
            retObj.Add(uu);
            break;
        }
    }
    return retObj;
}

I’ve really sent the above code through unit tests and also in practical deployment so it should be relatively bug free. I hope that by posting this here, you can tease out the logic from both the code and the comments and port it to your preferred language. In retrospect, maybe this bug wasn’t so obvious after all. When you think about it, each of the pieces of code and logic are doing exactly what we asked them to. It was only in chaining together an e-string subprofile function along with a local alignment that this thing ever even popped out. Either way, the fix is here. Hopefully it’ll be helpful.

Two Important Notes

  • Since global alignments represent the aligned sites of all sequences involved, no offset is used and hence this bug would never occur when specifying a subprofile of a global alignment.
  • Gaps that occur in the original sequence don’t affect retrieving the subprofile e-string calculation at all because those gaps DO exist and can be authentically counted into an offset.

Next

I haven’t had a need yet for porting the multiply function (earlier post) from global alignment logic to local alignment logic. If there are any insidious logical flaws like the ones I encountered this time, I’ll be sure to make a post.

Interactive Diagrams Using Inline JS+SVG+XHTML

featured post

without comments

>>> Attached: ( inline_svg_demo.xhtml — demo using the below technique for a silly drawing ) <<<

Compatibility: The following method doesn’t work for Internet Explorer 8, so if interoperability is really important to you, skip this series! Don’t worry though — Firefox 4.0, Internet Explorer 9 and WebKit (Safari / Chrome) will all natively support HTML5. HTML5 includes inline SVG for HTML so this ad hoc version can be retired when these three products take hold. The below method is compatible with Firefox 3.x, Opera 9.x and the current version of WebKit (Safari 5.x, Chrome 7.x). For display reasons, I’ll show off bitmap files inline while displaying inline SVG in linked pages.

Part One of Two

Part one is for everyone! I describe how to display an SVG inline on a webpage (XHTML) along with some JavaScript driven behaviour (clicking changes a displayed message and CSS defined colours).

Part two is for the bioinformaticians and phylogenists that survive and stick around — I utilize a parsed Newick format traversal to build up the drawing of a tree recursively — then I finish off by reintroducing JavaScript driven interactivity — tree nodes are highlighted via CSS when clicked, and information is displayed for the selected node (perhaps a multiple sequence alignment, a global alignment score, branch lengths — etc.).

( See part two here: Simple Interactive Phylogenetic Tree Sketches JS+SVG+XHTML )

Wait a sec — what do you mean by inline SVG?

Scalable Vector Graphics (SVG) is an XML vector image format. All browsers allow you to specify SVG files to be included into a webpage as an image. What we’ll be discussing here is how to include an SVG image as an XML element directly in a XHTML file — that is, the single webpage file will contain all of the markup for itself and its included SVG images without linking to any external files.

I won’t go into detail about how to make SVG images. You may create them with software like InkScape or by manual manipulation of the raw XML using a text or code editor (which I’ve done in this post). If you do end up using drawing software, simply copy out all of the XML elements with a text editor afterward and inject them into the XHTML page you create.

XHTML File Nuances

Include the following as a attribute of your html tag — this ensures that the browser reading your page knows which document type definition (DTD) to use.

<html xmlns="http://www.w3.org/1999/xhtml">
...
</html>

You must include your SVG with the following tags.

<svg xmlns="http://www.w3.org/2000/svg" width="300" height="300"
    xmlns:xlink="http://www.w3.org/1999/xlink">
...
</svg>

In the above, we parameterize the “svg” tag with the correct xml namespace (xmlns) as an attribute. We specify the width and height of the image we want to inline and also specify that we want to be able to use hyperlinks within the SVG image with “xlink”. Notice that both namespace attributes are meant to enable a client browser to correctly render your page.

SmileAgain Thumbnail

Smile Again!

A Silly Drawing and Target Behaviour

In this post, we’ll be manipulating a picture of a silly face I’ve named Smile Again — if you can render the attached file, the face looks like the one to the right.

The eventual target behaviour we want is as follows. Clicking on different parts of Smile Again’s face causes (1) the displayed message to change and (2) the clicked part to change colours. Notice that I’ll be using JavaScript (JS) to perform both of these functions. We’ll get into the specifics soon — if you’re feeling up to it, you can always go to the attached demo xhtml file now and look at the source to reverse engineer the thing.

We can create this face manually with the following inline SVG code — we’ll start with pure SVG and worry about the JavaScript behaviour afterward.

<svg xmlns="http://www.w3.org/2000/svg" width="300" height="300"
 xmlns:xlink="http://www.w3.org/1999/xlink">

<circle cx="100" cy="100" r="40"
 style="fill:white;stroke:black;stroke-width:2" />

<circle cx="100" cy="100" r="20"
 style="fill:black;stroke:black;stroke-width:2" />

<circle cx="200" cy="100" r="60"
 style="fill:white;stroke:black;stroke-width:2" />

<circle cx="200" cy="100" r="10"
 style="fill:black;stroke:black;stroke-width:2" />

<path d="M60 200 C80 300 220 300 240 200 L60 200"
 style="fill:white;stroke:black;stroke-width:2" />

<path d="M100 235 L200 235"
 style="fill:white;stroke:black;stroke-width:2" />

<path d="M130 224 L130 246"
 style="fill:white;stroke:black;stroke-width:2" />

<path d="M150 220 L150 250"
 style="fill:white;stroke:black;stroke-width:2" />

<path d="M170 224 L170 246"
 style="fill:white;stroke:black;stroke-width:2" />

</svg>

Each of the XML elements above correspond to a different piece of the SVG image. The first four circles draw the eyes, next the complicated path draws the mouth and the four simpler paths draw the stitches of the teeth.

Changing the Displayed Message with JavaScript

Now that we’ve drawn Smile Again, we want to have it react to mouse clicks. I want to have both the clicked features highlight as well as the displayed message to change. This would be an important feature of actual practical diagrams. To make this task a bit easier to explain, I will break down the behaviour into smaller components and functions and build them back up as I go along. This will also make it a lot easier to read.

Let us focus on just the two circles that makes up Smile Again’s right eye (on the left side of the diagram).

We want our clicks to call a function — in order for that to happen, we need to use anchor tags (hyperlinks) within the SVG. This is made possible with our careful xmlns (namespace) attributes we declared earlier!

<svg xmlns="http://www.w3.org/2000/svg" width="300" height="300"
    xmlns:xlink="http://www.w3.org/1999/xlink">

<a xlink:href="#" onclick="javascript:change_thoughts('Oww! My right cornea!')">
    <circle cx="100" cy="100" r="40"
        style="fill:white;stroke:black;stroke-width:2" />
</a>

<a xlink:href="#" onclick="javascript:change_thoughts('Ouch! My right iris!')">
    <circle cx="100" cy="100" r="20"
        style="fill:black;stroke:black;stroke-width:2" />
</a>
...
</svg>

Notice that we’ve now enclosed entire clickable elements in anchor tags and have used the “xlink:href” hyperlink attribute. In reality, we can have these links point to actual other webpages. Instead, these point to the same page and call the JavaScript function to change the text in the description called “change_thoughts”. In this case, the output is an expression of pain after being poked in the eye (“Oww! My right cornea!”). We can now define the “change_thoughts” function and also the place where we want our messages to appear.

<head>
...
function change_thoughts(say_this) {
    document.getElementById("thoughts").innerHTML = "\"" + say_this + "\""
}
...
</head>
<body>
...
<h2 id = "thoughts" style="font-family:monospace;">
</h2>
...
</body>

The “change_thoughts” function can go in the head or somewhere early in the body of the XHTML. It doesn’t really matter as long as it occurs before the SVG. This function gets an element with the id “thoughts” and changes the HTML inside of it. In the above, I add a pair of double quotes to the string to show. The element with the id “thoughts” is a h2 element — this really could been any valid HTML element. You do however have to place the “thoughts” element somewhere after  you declare the “change_thoughts” function.

Changing the Highlighting with JavaScript and CSS

The last thing we want to do is to change the highlighting by making our JavaScript function touch the style (CSS) of the inline SVG element that was clicked. We could use a separate CSS file, but in the spirit of making everything one inlined file, the CSS is specified with the “style” html attribute.

Again, we’ll be looking only at the two circles corresponding to Smile Again’s right eye (I’ve changed the way the lines are broken in the code listing below so that it will actually fit on this page).

<a xlink:href="#"
    onclick=
"javascript:change_thoughts('Oww! My right cornea!','right_cornea', '#F00')">
    <circle cx="100" cy="100" r="40"
        style="fill:white;stroke:black;stroke-width:2"
        id="right_cornea" />
</a>

<a xlink:href="#"
    onclick=
"javascript:change_thoughts('Ouch! My right iris!', 'right_iris', '#00F')">
    <circle cx="100" cy="100" r="20"
        style="fill:black;stroke:black;stroke-width:2"
        id="right_iris" />
</a>

Here are our two changes. First — each of the SVG elements now have a unique id (“right_cornea”, “right_iris”) so that we can use JS to manipulate each shape independently. Second — we’re now calling a new version of “change_thoughts” — this time, we’re giving the function the id of the shape to highlight and the colour (in hexadecimal) to highlight it with.

Finally, I can show you the new version of “change_thoughts” which will accept the new arguments and change the colours of the shapes.

function change_thoughts(say_this, highlight_this, colour) {

    // Say this ...
    document.getElementById("thoughts").innerHTML = "\"" + say_this + "\""

    // Highlight this ...
    var change = document.getElementById(highlight_this)
    change.setAttribute("style",
        "fill:" + colour + ";stroke:black;stroke-width:5")
}

In the above, we’ve added the code after “// Highlight this …”. We select the element “highlight_this” — remember, the id and colour are given by the function called with mouse clicks in the SVG. We then change the style by overwriting the inline CSS with the new specified fill colour (and a fat black stroke).

We now have a problem — each new mouseclick on a different shape will cause that new shape to highlight. The previous highlighted shape doesn’t automatically revert, and eventually all of the shapes are highlighted. What we need to do now is to improve the “change_thoughts” function one last time with some way to “undo” the previous change, so that new calls to this function will not only highlight the new shape, but will also remove the highlight from the previous shape.

Finishing with an Undo Object to remove previous changes

In the following version, we add an object to behave like a hash in the JavaScript. This hash will take html id attributes as keys, and the previous style as values. This behaves like an “undo” stack since we apply the styles to the elements with the corresponding id before we commit the changes to the new shape.

var _undo_render_hash = new Object

function change_thoughts(say_this, highlight_this, colour) {

    // Say this ...
    document.getElementById("thoughts").innerHTML = "\"" + say_this + "\""

    // Undo previous highlight ...
    for(var i in _undo_render_hash) {
        document.getElementById(i).setAttribute("style", _undo_render_hash[i])
    }

    // Reset undo stack ...
    _undo_render_hash = new Object

    // Highlight this ...
    var change = document.getElementById(highlight_this)
    _undo_render_hash[highlight_this] = change.getAttribute("style")
    change.setAttribute("style",
        "fill:" + colour + ";stroke:black;stroke-width:5")
}

Walking through the code, we start with a blank object as the undo stack — this is OK. On the first call of this function, since the undo object is empty, simply nothing is done before the new shape is highlighted. After the first call, all of the changes that were committed in the last call of the function are undone as their original styles are reapplied. We just add the next object to change to the undo object before applying changes to it.

The undo object is treated like a hash here, and in fact, we iterate all of its elements even though we only expect to have captured one id in the last call of the function. This means we can use this undo technique for functions that will modify many shapes with many changes since each of those modifications is iteratively undone in the next call. We must of course ensure that we save the original style of each of these shapes.

Next

At some point, I’ll want to revisit this and combine some of the other JavaScript stuff I’ve shown in this blog before. Displaying a rendered clickable phylogenetic tree is a task I figured out since I needed to quickly visualize some data for my team in my thesis work. Since I’ve already covered parsing a newick tree in this blog, it makes sense to complete the discussion with how I ended up with my phylogeny visualizer.

Eddie Ma

October 31st, 2010 at 9:23 pm

BioCompiler might start life as BioCOBOL

with 3 comments

Update: Matthew has found an even more thorough review paper that discusses computer assisted synthetic biology approaches — it can be found at doi:10.1016/j.copbio.2009.08.007.

The iGEM competition year is running to a close. The teams are headed into November 2010 and have roughly one month left before attending the conference. I’m personally not attending the conference this year — I think the undergraduates will get more out of the experience.

The current year sees our continued efforts to precipitate a dedicated software team — a functioning autonomous unit that will serve to supplement and enhance Waterloo’s impact in the iGEM competition. More importantly, we’re going to do some very interesting science. We’ve had some success talking with other student groups across campus — notably, we should probably talk with the student IEEE/CUBE chapter when we have more work completed. We had involved BIC as well, however, it was early on and we had even less ground to stand upon ( — the primordial software team was mostly interested in BioMortar and BrickLayer at the time — the later project having been taken up by the iGEM Coop students as in Python this year).

This whole BioCompiler business started when Matthew uncovered a nice candidate problem: the compilation of a schematic for behaviour to a fully functioning synthetic biological circuit. Let’s be as precise as we can be here. I mean to say, we will take a description (which could resemble a piece of formal language source code) — and have it compiled into the sequence of BioBricks that will produce the desired behaviour.

This idea has been approached by several groups before — but each time, a different subproblem was considered (this is a reorganization of Matthew’s very nice list here).

  • Synthetic biology programming language: Genetic Engineering of Living Cells (GEC) (Microsoft Research) is a project that hones in on a formal language specification.
  • Synthetic biology computer assisted design (CAD): Berkeley Software (iGEM 2009) created a suite of items — Eugene (formal language for synthetic biology), Spectacles (visualizer integrating parts with their behaviours), Kepler (a dataflow broker). As well, Berkeley is responsible for the award winning Clotho (iGEM 2008 – Best Software Tool) which is a workbench that connects with the parts registry database (amongst other possible resources).
  • Systems biology pathway reaction simulation: Systems Biology Markup Language (SBML), Systems Biology Workbench (SBW) and Jarnac are a set of tools that perform systems biology analysis (which we consider to be an output of synthetic biology). SBML is the formal language, Jarnac is a reaction network simulator (which utilizes JDesigner as a front end) and SBML is a dataflow broker between SBML and Jarnac.
  • BioBrick specific pathway simulation: Minnesota’s Team Comparator (iGEM 2008) created SynBioSS — a tool which estimates the concentration of reactants and products given the appropriate BioBricks on a simulated circuit.

Update — here are a few more items thanks to Matthew Gingerich, George Zarubin and Andre Masella.

  • Molecular biology and bioinformatics analysis: The European Molecular Biology Open Software Suite (EMBOSS) is a toolkit developed by the European Molecular Biology Network (EMBnet) for bioinformatics. This might not be immediately relevant, but it is interesting. EMBOSS is actually relatively complete. More distantly along this vein, there’s also Bioconductor which focuses mostly on microarray analysis and is implemented in R.
  • More synthetic biology computer aided design (CAD): TinkerCell is a GUI-driven piece of software that supports extension with C++ and Python. TinkerCell along with the suite created by Berkeley Software are the two most promising target systems in which to integrate BioCompiler. Finally, there’s GenoCAD which appears to be in an early phase of development — this software looks to emphasize construction with correct syntax and attribute grammars — I’ll have to read more about this.
  • Formal laboratory protocol description: BioCoder is another Microsoft research project — the designed language aims to be both human readable and complete for automation. It reminds me of standard operating procedures (SOP) with greater precision. While BioCoder compiles from protocol to automation, we’ll be compiling to circuitry and protocol. BioCoder will give us some insights about the kinds of protocols that others are thinking of.

There is of course more software, but these are the items that we have become most familiar with — that we like — and that we consider to be standards toward which our own work should strive.

Two interesting problems arise when we think about these subproblems. First, the programming languages specified (including Eugene from Berkeley’s CAD suite) are exactly what they claim to be. Formal specifications. This is possible because of how concrete they are. They literally document what a synthetic biology circuit is. But this isn’t too different from what humans have been doing in iGEM all along. Second, the synthetic biology items don’t really seem to talk to the systems biology items — whereas we expect that the two — being input and output — to be inextricably linked. I explain my thoughts on the two below.

Why a programming language?

Andre enlightened me to this the other day. Humans invented programming languages to do two things. On an arrow running from the concrete toward the esoteric, we have the practical concern of compressing the amount of code that we want to write while retaining our programs’ expressiveness. This is the origin of macro systems such as COBOL. On an arrow pointing in the reverse — from the cerebral abstract down to the literal — we have the theoretical concern of mathematical beauty, of completeness. This is the origin of such functional languages like Lisp. For humans to have any hope of creating such a Lisp-like language for synthetic biology — we would need to understand all of the reactionary nuances about it, the system with which we tamper — at least inasmuch as a painful heuristic approximation. This is a feat we are no where near completing though footholds are managed with systems biology.

So here we are, BioCOBOL is the first step. A developing simple — though complete macro-like system that is in league in terms of abstraction with the programming language / CAD -like projects we’ve seen thus far. Only, we aim to increase the efficiency of circuit design; so that the programming language is not a literal mapping of the human document — rather, it is an explanation of behaviour. We will abstract it ever so slightly with each iteration of development — departing further and further away from CAD. A subteam headed by Brandon is currently developing the syntax and search algorithms required for the job — Brandon suggests that a weighted traversal of valid circuits should form our algorithmic primitive. My subteam is attempting to characterize as many known circuits in iGEM as possible — analyzing what pieces of input (stimuli: chemical concentrations, gradients, quora, oscillators) may be compressed for their common usage — what function prototypes already exist ( — reacts to an input: promoters) for compression — what the standard outputs are (analogy to printing error messages: GFP-family) etc. — again in the effort to realize what is losslessly compressible. Our software should eventually provide the correct laboratory protocol as well.

In this sense, we are respecting what a compiler is before we even approach more sophisticated compilation: it transforms a document in one language to another.

Incidentally, I should mention that Jordan and George are working on a modelling problem with the lab team — I’m not clear on the specifics, but I take it they want to pull out some differential equations on a set of promoters.

Why do we care that Synthetic Biology logic should talk to Systems Biology logic?

Finally, it is clear that we will eventually want to become even more abstract — even more mathematically complete, even more expressive. While we may never know enough about systems biology to create BioLISP (in our lifetime), we expect there to be sufficient research for us to discover — and perhaps research we can conduct ourselves to come ever closer. Systems biology allows us to think about synthetic biology in terms of reaction concentrations; free energy etc.. It gives the notion of compilation its own ground; the ground we want to cover. Imagine the perfect BioCompiler — stating the a problem to be compiled in terms of the input and output of the system. Let’s be precise here: I mean to say, the products and reactants or behaviour of our circuit. Let us describe what our circuit will do instead of what our circuit is made of. This — the missing link, this compiler — is the logical final step of BioCompiler.

Eddie Ma

October 26th, 2010 at 10:57 am

Virus capsids are pretty

without comments

Brief: The majority of virus protein coats (capsids) are in the shape of an icosahedron — a figure with twenty equilateral triangles. The first time I saw this rendered was in a paper by David S. Goodsell. In it, Goodsell describes proteins with structural symmetries. Four viruses are used as examples — they are tobacco necrosis virus (2BUK), tomato bushy stunt virus (2TBV), bluetongue virus (3IYK) and simian virus 40 (1SVA) linked here to their RSCB PDB entries.

Pretty, aren’t they? Very pretty.

If you follow the PDB links, you can take a look at how a single tessellation unit appears, how long a chain is and how massive the capsid is.

Notice: The icosahedron (twenty equilateral triangles) must not be confused with the dodecahedron (twenty points).

Eddie Ma

September 27th, 2010 at 1:00 pm

C# & Bioinformatics: Indexers & Substitution Matrices

without comments

I’ve recently come to appreciate the convenience of C# indexers. What indexers allow you to do is to subscript an object using the familiar bracket notation. I’ve used them for substitution matrices as part of my phylogeny project. Indexers are essentially syntactical sugar that obey the rules of method overloading. I first describe what I think are useful substitution matrix indexers and then a bare bones substitution matrix class (you could use your own). The indexer notation implementation is discussed last, so feel free to skip the preamble if you’re able to deduce what I’m doing.

Note: I’ve only discussed accessors (getters) and not mutators (setters) today.

Some Reasonable Substitution Matrix Indexers

This is the notation you might expect from the indexer notation in C#.

// Let there be a class called SubMatrix which contains data from BLOSUM62.

var  sm = new SubMatrix( ... );
     // I'll assume you already have some constructors.

int  index_of_proline = sm['P'];
     // Returns the row or column that corresponds to proline, 14.

char token_at_three = sm[3];
     // Returns the amino acid at position three, aspartate.

int  score_proline_to_aspartate = sm['P', 'D'];
     // Returns the score for a mutation from proline to aspartate, -1.

int  score_aspartate_to_proline = sm[3, 14];
     // Returns the score for a mutation from aspartate to proline, -1.

An Example Bare Bones Substitution Matrix Class

Let’s say you’ve loaded up the BLOSUM62 and are representing it internally in some 2D array…

// We've keying the rows and columns in the order given by BLOSUM62:
// ARNDCQEGHILKMFPSTWYVBZX* (24 rows, columns)

int[,] imatrix;

For convenience, let’s say you’ll also keep a dictionary to map at which array position one finds each amino acid…

// Keys = amino acid letters, Values = row or column index

Dictionary<char, int> indexOfAA;

Finally, we’ll put these two elements into a class and assume that you’ve already written your own constructors that will take care of the above two items — either from accepting arrays and strings as arguments or by reading from raw files. If this isn’t true and you need more help, feel free to leave a comment and I’ll extend this bare bones example.

// Bare bones class ...
public partial class SubMatrix {

    // Properties ...
    private int[,] imatrix;
    private Dictionary<char, int> indexOfAA;

    // Automatic Properties ...
    public int Width { // Returns number of rows of the matrix.
        get {
            return imatrix.GetLength(0);
        }
    }

    // Constructors ...
    ...
}

I’ve added the automatic property “Width” above — automatic properties are C# members that provide encapsulation: a public face to some arbitrary backend data — you’ve been using these all along when you’ve called “List.Count” or “Array.Length“.

Substitution Matrix Indexer Implementation

You can implement the example substitution matrix indexers as follows. Notice the use of the “this” keyword and square[] brackets[] to specify[] arguments[].

//And finally ...
public partial class SubMatrix {

    // Indexers ...
    // Give me the row or column index where I can find the token "aa".
    public int this[char aa] {
        get {
            return this.indexOfAA[aa];
        }
    }
    // Give me the amino acid at the row or column index "index".
    public char this[int index] {
        get {
            return this.key[index];
        }
    }
    // Give me the score for mutating token "row" to token "column".
    public double this[char row, char column] {
        get {
            return this.imatrix[this.indexOfAA[row], this.indexOfAA[column]];
        }
    }
    // Give me the score for mutating token at index "row" to index "column".
    public double this[int row, int column] {
        get {
            return this.imatrix[row, column];
        }
    }
}

Similar constructs are available in Python and Ruby but not Java. I’ll likely cover those later as well as how to set values too.

Eddie Ma

September 1st, 2010 at 11:15 am

C# & Bioinformatics: Align Strings to Edit Strings

without comments

This post follows roughly from the e-strings (R. Edgar, 2004) topic that I posted about here. The previous source code was listed in Ruby and JS, but this time I’ve used C#.

In this post, I discuss alignment strings (a-strings), then why and how to convert them into edit strings (e-strings). Incidentally, I can’t seem to recover where I first saw alignment strings so tell me if you know.

Alignment strings

When you perform a pairwise sequence alignment, the transformation that you perform on the two sequences is finite. You can record precisely where insertions, deletions and substitutions (or matches) are made. This is useful if you want to retain the original sequences, or later on build a multiple sequence alignment while keeping the history of modifications. There’s a datastructure I’ve seen described called alignment strings, and in it you basically list out the characters ‘I’, ‘D’ and ‘M’ to describe the pairwise alignment.

Consider the example two protein subsequences below.

gi|115372949: GQAGDIRRLQSFNFQTYFVRHYN
gi|29827656:  SLSTGVSRSFQSVNYPTRYWQ

A global alignment of the two using the BLOSUM62 substitution matrix with the gap penalties -12 to open, -1 to extend yields the following.

gi|115372949: GQAGDIR-----RLQSFNFQTYFVRHYN
gi|29827656:  SL----STGVSRSFQSVNYPT---RYWQ

The corresponding alignment string looks like this…

gi|115372949: GQAGDIR-----RLQSFNFQTYFVRHYN
gi|29827656:  SL----STGVSRSFQSVNYPT---RYWQ
Align String: MMDDDDMIIIIIMMMMMMMMMDDDMMMM

Remember, the alignment is described such that the top string is modelled as occurring earlier than the bottom string which occurs later — this is why a gap in the top string is an insertion (that new material is inserted in the later, bottom string) while a gap in the bottom string is a deletion (that material is deleted to make the later string). Notice in reality, it doesn’t really matter what’s on top and what’s on bottom– the important thing is that the alignment now contains gaps.

Why we should probably use e-strings instead

An alignment string describes the relationship between two strings and their gaps — we are actually recording some redundant information if we only want to take one string at a time into consideration and the path needed to construct profiles it participates in. The top and bottom sequences are also treated differently, where both ‘M’ and ‘D’ indicates retention of the original characters for the top string and both ‘M’ and ‘I’ indicates retention for the bottom string; the remaining character of course implies the inclusion of a gap character. A pair of e-strings would give each of these sequences their own data structure and allow us to cleanly render the sequences as they appear in the deeper nodes of a phylogenetic tree using the multiply operation described last time.

Here are the corresponding e-strings for the above examples.

gi|115372949: GQAGDIR-----RLQSFNFQTYFVRHYN
e-string:     < 7, -5, 16 >
gi|29827656:  SL----STGVSRSFQSVNYPT---RYWQ
e-string:     < 2, -4, 16, -3, 4 >

Recall that an e-string is a list of alternating positive and negative integers; positive integers mean to retain a substring of the given length from the originating sequence, and negative integers mean to place in a gap of the given length.

Converting a-strings to e-strings

Below is a C# code listing for an implementation I used in my project to convert from a-strings to the more versatile e-strings. The thing is — I really don’t use a-strings to begin with anymore. In earlier versions of my project, I used to keep track of my movement across the score matrix using an a-string by dumping down a ‘D’ for a vertical hop, a ‘I’ for a horizontal hop and a ‘M’ for a diagonal hop. I now just count the number of relevant steps to generate the matching pair of e-strings.

The function below, astring_to_estring takes an a-string (string u) as an argument and returns two e-strings in a list (List<List<int>>) — don’t let that type confuse you, it simply breaks down to a list with two elements in it: the e-string for the top sequence (List<int> v), and the e-string for the bottom sequence (List<int> w).

public static List<List<int>> astring_to_estring(string u) {
    /* Defined elsewhere are the constant characters ...
       EDITSUB = 'M';
       EDITINS = 'I';
       EDITDEL = 'D';
    */
    var v = new List<int>(); // Top e-string
    var w = new List<int>(); // Bottom e-string
    foreach(var uu in u) {
        if(uu == EDITSUB) { // If we receive a 'M' ...
            if(v.Count == 0) { // Working with e-string v (top, keep)
                v.Add(1);
            } else if(v[v.Count -1] <= 0) {
                v.Add(1);
            } else {
                v[v.Count -1] += 1;
            }
            if(w.Count == 0) { // Working with e-string w (bottom, gap)
                w.Add(1);
            } else if(w[w.Count -1] <= 0) {
                w.Add(1);
            } else {
                w[w.Count -1] += 1;
            }
        } else if(uu == EDITINS) { // If we receive a 'I' ...
            if(v.Count == 0) { // Working with e-string v (top, gap)
                v.Add(-1);
            } else if(v[v.Count -1] >= 0) {
                v.Add(-1);
            } else {
                v[v.Count -1] -= 1;
            }
            if(w.Count == 0) { // Working with e-string w (bottom, keep)
                w.Add(1);
            } else if(w[w.Count -1] <= 0) {
                w.Add(1);
            } else {
                w[w.Count -1] += 1;
            }
        } else if(uu == EDITDEL) { // If we receive a 'D' ...
            if(v.Count == 0) { // Working with e-string v (top, keep)
                v.Add(1);
            } else if(v[v.Count -1] >= 0) {
                v[v.Count -1] += 1;
            } else {
                v.Add(1);
            }
            if(w.Count == 0) { // Working with e-string w (bottom, keep)
                w.Add(-1);
            } else if(w[w.Count -1] >= 0) {
                w.Add(-1);
            } else {
                w[w.Count -1] -= 1;
            }
        }
    }
    var vw = new List<List<int>>(); // Set up return list ...
    vw.Add(v); // Top e-string v added ...
    vw.Add(w); // Bottom e-string w added ...
    return vw;
}

The conversion back from e-strings to a-strings is also easy, but I don’t cover that today. Enjoy and happy coding 😀

Eddie Ma

August 21st, 2010 at 9:03 pm