Ed's Big Plans

Computing for Science and Awesome

The Right Tool (Why I chose Java to do RSA)

without comments

Brief: I learned something valuable last week when working on this RSA encryption/decryption assignment for my Computer Security class. It’s important to be versatile when doing computer science — we must ensure we always use the most efficient tool. If we aren’t versatile, we risk taking tremendous amounts of time trying to reimplement something that already exists elsewhere.

So, what tool is the right tool to quickly throw together RSA encryption?

It turns out that Java does an excellent job. Its BigInteger class has all the ingredients you’ll ever need.

// This function generates a new probable prime ...
BigInteger p = BigInteger.probablePrime(bits, random);

// This function performs modulus inverse ...
BigInteger d = e.modInverse(phi_n);

// These functions can be used to check your work ...
BigInteger one = d.multiply(e).mod(phi_n);
BigInteger one = p.gcd(q);

// This function is at the heart of RSA ...
BigInteger C = M.modPow(d, n);

Before I looked at the Java documentation, I had plans to do this with Python and some of my classmates had plans with MATLAB. It’s not that these are inherently bad technologies — they’re just not the right tool.

As much as I have gripes with Java, I’ve got to say that this made me and a lot of my class very happy 😀

Eddie Ma

February 24th, 2011 at 6:26 pm

Posted in Pure Programming

Tagged with , ,

Move a Subdirectory to a Subdomain (Apache mod_alias)

without comments

In this post, I’m going to explain how to move resources located at a subdirectory of your domain to its own subdomain. The reason why I came across this was because I decided to move my MediaWiki installation to better separate the URLs corresponding to my blog and my notebook — this is paired with better separation of the actual files and directories corresponding to my blog and notebook as they exist on my server. I’ll assume that you’re running Apache; I’ll further assume that you’ve already gotten your A-RECORD set up with your domain registrar or name server provider (DNS) — if you don’t know what that means, ask me.

I’ll use my MediaWiki move as my running example for this post.

Here are some properties we want a redirect to satisfy — I’ll lead you through to making a redirect such that each of these conditions are met.

  1. Moving the subdirectory to a subdomain corresponds to moving files on the server filesystem.
    • (URL) http://eddiema.ca/wiki = ~/Sites/eddiema/notebook (server filesystem)
    • (URL) http://wiki.eddiema.ca = ~/Sites/notebook (server filesystem)
  2. Requests to the subdomain must work.
    • http://wiki.eddiema.ca
  3. Requests to the subdirectory must be forwarded to the subdomain.
    • http://eddiema.ca/wiki →
    • http://wiki.eddiema.ca
  4. Requests to subdirectories of the former subdirectory must be forwarded to subdirectories of the subdomain.
    • http://eddiema.ca/wiki/index.php/Special:AllPages
    • http://wiki.eddiema.ca/index.php/Special:AllPages
  5. URLs with GET requests to the subdirectory must be forwarded to the subdomain without losing the GET requests.
    • http://eddiema.ca/wiki/?search=Notes+CIS+6050
    • http://wiki.eddiema.ca/?search=Notes+CIS+6050

Let’s meet each of our above properties.

Changing the URL corresponds to a physical move (Property 1)

Move your subdirectory on the server to its target location on the server.

For my example, I moved my notebook out of my main site directory.

mv ~/Sites/eddiema/notebook ~/Sites/notebook

(If you’re moving a MediaWiki installation like me, there’s one more thing you have to do — this is explained at the end.)

Requests to the subdomain must work (Property 2)

Update your Apache httpd.conf configuration file with a new VirtualHost entry. This new entry stipulates a new document root for your target subdomain.

#<VirtualHost *:80>
#    ServerName {full URL to your new subdomain}
#    DocumentRoot {full path to your new document root}
#</VirtualHost>

# Example ...
<VirtualHost *:80>
    ServerName wiki.eddiema.ca
    DocumentRoot /Users/eddiema/Sites/notebook
</VirtualHost>

Requests to the subdirectory must be forwarded to the subdomain (Properties 3~5)

There are two ways to accomplish this task — you can either modify Apache’s httpd.conf file, or you can create a new .htaccess file. I chose the latter because I wanted to keep all of the changes related to this move localized.

To accomplish this task, you must perform two steps.

First, you must recreate the subdirectory on your server filesystem that you have just moved elsewhere. We do this so that a request for that directory causes Apache to look there for a .htaccess file in the correct place. This .htaccess file will tell Apache to redirect to your new subdomain instead.

For my example, I recreated the subdirectory with…

mkdir ~/Sites/eddiema/notebook

Second, Create a new plaintext .htaccess file in the recreated subdirectory with a RedirectMatch rule.

#RedirectMatch 301 ^/{subdirectory}/(.*)$ http://{full subdomain URL}/$1

# Example ...
RedirectMatch 301 ^/wiki/(.*)$ http://wiki.eddiema.ca/$1

RedirectMatch is part of mod_alias. If you’re doing something more sophisticated than just forwarding URLs, then you will need to use mod_redirect which has the ability to manipulate query strings.

The number 301 tells browsers and search engines that the URL pointing at your subdirectory is no longer valid and has permanently moved to the one pointing at your new subdomain. The next part is a PERL style regular expression. Let’s break it apart.

  • ^ and $ match the start of, and end of the string respectively.
  • .* means we want to match zero or more characters in a string.
  • (.*) means we want to save the matched string as an autovariable $1.

The last string is appended with $1 to push the tail of any URL of your subdirectory onto the tail of your subdomain.

A final note for moving a MediaWiki installation

In order to stop MediaWiki from appending the subdirectory “wiki” to the tail of your new subdomain (“http://wiki.{your domain}.*/wiki“), you can change the variable $wgScriptPath to the empty string in LocalSettings.php found in the root directory of your MediaWiki installation.

#$wgScriptPath       = "/wiki";
$wgScriptPath       = "";

Note that MediaWiki doesn’t seem to like this idea — but they don’t really say explain why.

After this step, you can use Short URLs (same as above link) to get rid of index.php appearing in the URL, but I won’t go into that today.

That’s everything 😀

Eddie Ma

February 23rd, 2011 at 12:36 pm

C & Math: Sieve of Eratosthenes with Wheel Factorization

featured post

without comments

In the first assignment of Computer Security, we were to implement The Sieve of Eratosthenes. The instructor gives a student the failing grade of 6/13 for a naive implementation, and as we increase the efficiency of the sieve, we get more marks. There are the three standard optimizations: (1) for the current prime being considered, start the indexer at the square of the current prime; (2) consider only even numbers; (3) end crossing out numbers at the square root of the last value of the sieve.

Since the assignment has been handed in, I’ve decided to post my solution here as I haven’t seen C code on the web which implements wheel factorization.

We can think of wheel factorization as an extension to skipping all even numbers. Since we know that all even numbers are multiples of two, we can just skip them all and save half the work. By the same token, if we know a pattern of repeating multiples corresponding to the first few primes, then we can skip all of those guaranteed multiples and save some work.

The wheel I implement skips all multiples of 2, 3 and 5. In Melissa O’Neill’s The Genuine Sieve of Erastothenes, an implementation of the sieve with a priority queue optimization is shown in Haskell while wheel factorization with the primes 2, 3, 5 and 7 is discussed. The implementation of that wheel (and other wheels) is left as an exercise for her readers 😛

But first, let’s take a look at the savings of implementing this wheel. Consider the block of numbers in modulo thirty below corresponding to the wheel for primes 2, 3 and 5 …

0 1 2 3 4 5 6 7 8 9
10 11 12 13 14 15 16 17 18 19
20 21 22 23 24 25 26 27 28 29

Only the highlighted numbers need to be checked to be crossed out during sieving since the remaining values are guaranteed to be multiples of 2, 3 or 5. This pattern repeats every thirty numbers which is why I say that it is in modulo thirty. We hence skip 22/30 of all cells by using the wheel of thirty — a savings of 73%. If we implemented the wheel O’Neill mentioned, we would skip 77% of cells using a wheel of 210 (for primes 2, 3, 5 and 7).

(Note that the highlighted numbers in the above block also correspond to the multiplicative identity one and numbers which are coprime to 30.)

Below is the final code that I used.

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

const unsigned int SIEVE = 15319000;
const unsigned int PRIME = 990000;

int main(void) {
    unsigned char* sieve = calloc(SIEVE + 30, 1); // +30 gives us incr padding
    unsigned int thisprime = 7;
    unsigned int iprime = 4;

    unsigned int sieveroot = (int)sqrt(SIEVE) +1;

    // Update: don't need to zero the sieve - using calloc() not malloc()

    sieve[7] = 1;

    for(; iprime < PRIME; iprime ++) {
        // ENHANCEMENT 3: only cross off until square root of |seive|.
        if(thisprime < sieveroot) {
            // ENHANCEMENT 1: Increment by 30 -- 4/15 the work.
            // ENHANCEMENT 2: start crossing off at prime * prime.
            int i = (thisprime * thisprime);
            switch (i % 30) { // new squared prime -- get equivalence class.
                case 1:
                    if(!sieve[i] && !(i % thisprime)) {sieve[i] = 1;}
                    i += 6;
                case 7:
                    if(!sieve[i] && !(i % thisprime)) {sieve[i] = 1;}
                    i += 4;
                case 11:
                    if(!sieve[i] && !(i % thisprime)) {sieve[i] = 1;}
                    i += 2;
                case 13:
                    if(!sieve[i] && !(i % thisprime)) {sieve[i] = 1;}
                    i += 4;
                case 17:
                    if(!sieve[i] && !(i % thisprime)) {sieve[i] = 1;}
                    i += 2;
                case 19:
                    if(!sieve[i] && !(i % thisprime)) {sieve[i] = 1;}
                    i += 4;
                case 23:
                    if(!sieve[i] && !(i % thisprime)) {sieve[i] = 1;}
                    i += 6;
                case 29:
                    if(!sieve[i] && !(i % thisprime)) {sieve[i] = 1;}
                    i += 1; // 29 + 1 (mod 30) = 0 -- just in step
            }
            for(; i < SIEVE; i += 30) {
                if(!sieve[i+1] && !((i+1) % thisprime)) sieve[i+1] = 1;
                if(!sieve[i+7] && !((i+7) % thisprime)) sieve[i+7] = 1;
                if(!sieve[i+11] && !((i+11) % thisprime)) sieve[i+11] = 1;
                if(!sieve[i+13] && !((i+13) % thisprime)) sieve[i+13] = 1;
                if(!sieve[i+17] && !((i+17) % thisprime)) sieve[i+17] = 1;
                if(!sieve[i+19] && !((i+19) % thisprime)) sieve[i+19] = 1;
                if(!sieve[i+23] && !((i+23) % thisprime)) sieve[i+23] = 1;
                if(!sieve[i+29] && !((i+29) % thisprime)) sieve[i+29] = 1;
            }
        }

        {
            int i = thisprime;
            switch (i % 30) { // write down the next prime in 'thisprime'.
                case 1:
                    if(!sieve[i]) {thisprime = i; sieve[i] = 1; goto done;}
                    i += 6;
                case 7:
                    if(!sieve[i]) {thisprime = i; sieve[i] = 1; goto done;}
                    i += 4;
                case 11:
                    if(!sieve[i]) {thisprime = i; sieve[i] = 1; goto done;}
                    i += 2;
                case 13:
                    if(!sieve[i]) {thisprime = i; sieve[i] = 1; goto done;}
                    i += 4;
                case 17:
                    if(!sieve[i]) {thisprime = i; sieve[i] = 1; goto done;}
                    i += 2;
                case 19:
                    if(!sieve[i]) {thisprime = i; sieve[i] = 1; goto done;}
                    i += 4;
                case 23:
                    if(!sieve[i]) {thisprime = i; sieve[i] = 1; goto done;}
                    i += 6;
                case 29:
                    if(!sieve[i]) {thisprime = i; sieve[i] = 1; goto done;}
                    i += 1;
            }
            for(; i < SIEVE; i += 30) {
                if(!sieve[i+1]) {thisprime = i+1; sieve[i+1] = 1; goto done;}
                if(!sieve[i+7]) {thisprime = i+7; sieve[i+7] = 1; goto done;}
                if(!sieve[i+11]) {thisprime = i+11; sieve[i+11] = 1; goto done;}
                if(!sieve[i+13]) {thisprime = i+13; sieve[i+13] = 1; goto done;}
                if(!sieve[i+17]) {thisprime = i+17; sieve[i+17] = 1; goto done;}
                if(!sieve[i+19]) {thisprime = i+19; sieve[i+19] = 1; goto done;}
                if(!sieve[i+23]) {thisprime = i+23; sieve[i+23] = 1; goto done;}
                if(!sieve[i+29]) {thisprime = i+29; sieve[i+29] = 1; goto done;}
            }
            done:;
        }
    }
    printf("%d\n", thisprime);
    free(sieve);
    return 0;
}

Notice that there is a switch construct — this is necessary because we aren’t guaranteed that the first value to sieve for a new prime (or squared prime) is going to be an even multiple of thirty. Consider sieving seven — the very first prime to consider. We start by considering 72 = 49. Notice 49 (mod 30) is congruent to 19. The switch statement incrementally moves the cursor from the 19th equivalence class to the 23rd, to the 29th before pushing it one integer more to 30 — 30 (mod 30) is zero — and so we are able to continue incrementing by thirty from then on in the loop.

The code listed is rigged to find the 990 000th prime as per the assignment and uses a sieve of predetermined size. Note that if you want to use my sieve code above to find whichever prime you like, you must also change the size of the sieve. If you take a look at How many primes are there? written by Chris K. Caldwell, you’ll notice a few equations that allow you to highball the nth prime of your choosing, thereby letting you calculate that prime with the overshot sieve size.

Note also that this sieve is not the most efficient. A classmate of mine implemented The Sieve of Atkin which is magnitudes faster than this implementation.

Eddie Ma

February 3rd, 2011 at 11:08 pm

Add ‘Day/Week/Month’ to WordPress.com Stats dashboard widget

without comments

Brief: The WordPress.com Stats plug-in does exactly what it’s supposed to. In the Dashboard, it shows you your top posts, top search queries and your most active posts. The problem here is that it won’t say for what period of time it’s showing you these lists. In this post, I’ll show you a quick hack to do just that so that each of those list titles will now be accompanied by the text Day, Week or Month.

In case you didn’t know, you can change over which period you want each of those lists to be shown by mousing over the top-right corner of the WordPress.com Stats Dashboard widget and clicking the revealed text, configure. Below are two images that illustrate our objective today.

The screenshot on the left (before) has the default list headers and the one on the right (after) has period labels.

The PHP file you have to edit is wp-content/plugins/stats/stats.php — you can open this file on your filesystem — or if your WordPress installation is hosted elsewhere, you can use the built-in editor. To get to this file using the built-in editor, go to your dashboard (wp-admin), then click on PluginsEditor — then in the drop-down labelled Select plugin to edit, choose WordPress.com Stats.

Once you have access to the file, look for a block of HTML that looks like the following. There’s only one div with the id stats-info — so you can use that to string search over the file. The ellipses (“…”) means that I’m skipping several irrelevant lines.

<div id="stats-info">
	<div id="top-posts" class='stats-section'>
		<div class="stats-section-inner">
		<h4 class="heading"><?php _e( 'Top Posts' , 'stats'); ?></h4>
		...
		<h4 class="heading"><?php _e( 'Top Searches' , 'stats'); ?></h4>
		...
		<h4 class="heading"><?php _e( 'Most Active' , 'stats'); ?></h4>
		...
		</div>
	</div>
</div>

You’ll next need to add the following line directly above <div id=”stats-info”> . What it does is rebuild a dictionary (associative array) so that a period of one means “Day”, seven means “Week”, and thirty-one means “Month”. If you’ve discovered a different dictionary that has these values keyed in reverse, let me know since I couldn’t find it.

<?php $re_periods = array('1' => 'Day', '7' => 'Week', '31' => 'Month'); ?>

Finally, replace the h4 elements in the original code with the h4 elements in the below code. This completes the customization.

<div id="stats-info">
	<div id="top-posts" class='stats-section'>
		<div class="stats-section-inner">
		<h4 class="heading">
	<?php _e( 'Top Posts &middot; '.$re_periods[$options[top]] , 'stats'); ?>
		</h4>
		...
		<h4 class="heading">
	<?php _e( 'Top Searches &middot; '.$re_periods[$options[search]] , 'stats'); ?>
		</h4>
		...
		<h4 class="heading">
	<?php _e( 'Most Active &middot; '.$re_periods[$options[active]] , 'stats'); ?>
		</h4>
		...
		</div>
	</div>
</div>

The above code uses the $options dictionary to look for the correct period, then uses that period in our new $re_periods dictionary to look for the correct name. That name is then concatenated to the string inside the h4 element and serves as our new title or list header inside the dashboard.

Enjoy!

Eddie Ma

January 28th, 2011 at 12:52 pm

C# & Science: A 2D Heatmap Visualization Function (PNG)

featured post

without comments

Today, let’s take advantage of the bitmap writing ability of C# and output some heatmaps. Heatmaps are a nice visualization tool as they allow you to summarize numeric continuous values as colour intensities or hue spectra. It’s far easier for the human to grasp a general trend in data via a 2D image than it is to interpret a giant rectangular matrix of numbers. I’ll use two examples to demonstrate my heatmap code. The function provided is generalizable to all 2D arrays of doubles so I welcome you to download it and try it yourself.

>>> Attached: ( tgz | zip — C# source: demo main, example data used below, makefile ) <<<

Let’s start by specifying the two target colour schemes.

  • RGB-faux-colour {blue = cold, green = medium, red = hot} for web display
  • greyscale {white = cold, black = hot} suitable for print

Next, let’s take a look at the two examples.

Example Heat Maps 1: Bioinformatics Sequence Alignments — Backtrace Matrix

Here’s an application of heat maps to sequence alignments. We’ve visualized the alignment traces in the dynamic programming score matrices of local alignments. Here, a pair of tryptophan-aspartate-40 protein sequences are aligned so that you can quickly pick out two prominent traces with your eyes. The highest scoring spot is in red — the highest scoring trace carries backward and upward on a diagonal from that spot.

The left heatmap is shown in RGB-faux-colour decorated with a white diagonal. The center heatmap is shown in greyscale with no decorations. The right heatmap is decorated with a major diagonal, borders, lines every 100 values and ticks every ten values. Values correspond to amino acid positions where the top left is (0, 0) of the alignment.

Example Heat Maps 2: Feed Forward Back Propagation Network — Weight Training

Here’s another application of heat maps. This time we’ve visualized the training of neural network weight values. The weights are trained over 200 epochs to their final values. This visualization allows us to see the movement of weights from pseudorandom noise (top) to their final values (bottom).

Shown above are the weight and bias values for a 2-input, 4-hidden-node, 4-hidden-node, 2-output back-propagation feed-forward ANN trained on the toy XOR-EQ problem. The maps are rendered as RGB-faux-colour (left); greyscale (center); and RGB-faux-colour decorated with horizontal lines every 25 values (=50px) (right). Without going into detail, the leftmost 12 columns bridge inputs to the first hidden layer, the next 20 columns belong to the next weight layer, and the last 10 columns bridge to the output layer.

Update: I changed the height of the images of this example to half their original height — the horizontal lines now occur every 25 values (=25px).

C# Functions

This code makes use of two C# features — (1) function delegates and (2) optional arguments. If you also crack open the main file attached, you’ll notice I also make use of a third feature — (3) named arguments. A function delegate is used so that we can define and select which heat function to use to transform a three-tuple of double values (real, min, max) into a three-tuple of byte values (R, G, B). Optional arguments are used because the function signature has a total of twelve arguments. Leaving some out with sensible default values makes a lot of sense. Finally, I use named arguments in the main function because they allow me to (1) specify my optional arguments in a logical order and (2) read which arguments have been assigned without looking at the function signature.

Heat Functions

For this code to work, heat functions must match this delegate signature. Arguments: val is the current value, min is the lowest value in the matrix and max is the highest value in the matrix; we use min and max for normalization of val.

public delegate byte[] ValueToPixel(double val, double min, double max);

This is the RGB-faux-colour function that breaks apart the domain of heats and assigns it some amount of blue, green and red.

public static byte[] FauxColourRGB(double val, double min, double max) {
    byte r = 0;
    byte g = 0;
    byte b = 0;
    val = (val - min) / (max - min);
           if(               val <= 0.2) {
        b = (byte)((val / 0.2) * 255);
    } else if(val >  0.2 &&  val <= 0.7) {
        b = (byte)((1.0 - ((val - 0.2) / 0.5)) * 255);
    }
           if(val >= 0.2 &&  val <= 0.6) {
        g = (byte)(((val - 0.2) / 0.4) * 255);
    } else if(val >  0.6 &&  val <= 0.9) {
        g = (byte)((1.0 - ((val - 0.6) / 0.3)) * 255);
    }
           if(val >= 0.5               ) {
        r = (byte)(((val - 0.5) / 0.5) * 255);
    }
    return new byte[]{r, g, b};
}

This is a far simpler greyscale function.

public static byte[] Greyscale(double val, double min, double max) {
    byte y = 0;
    val = (val - min) / (max - min);
    y = (byte)((1.0 - val) * 255);
    return new byte[]{y, y, y};
}

Heatmap Writer

The function below is split into a few logical parts: (1) we get the minimum and maximum heat values to normalize intensities against; (2) we set the pixels to the colours we want; (3) we add the decorations (borders, ticks etc.); (4) we save the file.

public static void SaveHeatmap(
    string fileName, double[,] matrix, ValueToPixel vtp,
    int pixw = 1, int pixh = 1,
    Color? decorationColour = null,
    bool drawBorder = false, bool drawDiag = false,
    int hLines = 0, int vLines = 0, int hTicks = 0, int vTicks = 0)
{
    var rows = matrix.GetLength(0);
    var cols = matrix.GetLength(1);
    var bitmap = new Bitmap(rows * pixw, cols * pixh);

    //Get min and max values ...
    var min = Double.PositiveInfinity;
    var max = Double.NegativeInfinity;
    for(int i = 0; i < matrix.GetLength(0); i ++) {
        for(int j = 0; j < matrix.GetLength(1); j ++) {
            max = matrix[i,j] > max? matrix[i,j]: max;
            min = matrix[i,j] < min? matrix[i,j]: min;
        }
    }
    //Set pixels ...
    for(int j = 0; j < bitmap.Height; j ++) {
        for(int i = 0; i < bitmap.Width; i ++) {
            var triplet = vtp(matrix[i / pixw, j / pixh], min, max);
            var color = Color.FromArgb(triplet[0], triplet[1], triplet[2]);
            bitmap.SetPixel(i, j, color);
        }
    }
    //Decorations ...
    var dc = decorationColour ?? Color.Black;
    for(int i = 0; i < bitmap.Height; i ++) {
        for(int j = 0; j < bitmap.Width; j ++) {
            if(drawBorder) {
                if(i == 0 || i == bitmap.Height -1) {
                    // Top and Bottom Borders ...
                    bitmap.SetPixel(j, i, dc);
                } else if(j == 0 || j == bitmap.Width -1) {
                    // Left and Right Borders ...
                    bitmap.SetPixel(j, i, dc);
                }
            }
            if(bitmap.Width == bitmap.Height && drawDiag && (i % 2 == 0)) {
                // Major Diagonal ... (only draw if square image +explicit)
                bitmap.SetPixel(i, i, dc);
            }
            //Zeroed lines and zeroed ticks are turned off.
            if(hLines != 0 && i % (hLines*pixh) == 0) {
                if(j % (2*pixw) == 0) {
                    //Horizontal Bars ...
                    bitmap.SetPixel(j, i, dc);
                }
            } else if(hTicks != 0 && i % (hTicks*pixh) == 0) {
                // Dots: H-Spacing
                if(vTicks != 0 && j % (vTicks*pixw) == 0) {
                    // Dots: V-Spacing
                    bitmap.SetPixel(j, i, dc);
                }
            } else if(i % (2*pixh) == 0) {
                if(vLines != 0 && j % (vLines*pixw) == 0) {
                    //Vertical Bars
                    bitmap.SetPixel(j, i, dc);
                }
            }
        }
    }

    //Save file...
    bitmap.Save(fileName, ImageFormat.Png);
    bitmap.Dispose();
}

Happy mapping 😀

Compatibility Notes: The C# code discussed was developed with Mono using Monodevelop. All code is compatible with the C#.NET 4.0 specification and is forward compatible. I’ve linked in two assemblies — (1) -r:System.Core for general C#.Net 3.5~4.0 features, and (2) r:System.Drawing to write out PNG files. If you need an overview of how to set up Monodevelop to use the System.Drawing assembly, see C Sharp PNG Bitmap Writing with System.Drawing in my notebook.

Eddie Ma

January 21st, 2011 at 1:20 pm

Simple Interactive Phylogenetic Tree Sketches JS+SVG+XHTML

featured post

without comments

>>> Attached: ( parse_newick_xhtml.js | draw_phylogeny.js — implementation in JS | index.html — demo with IL-1 ) <<<
>>> Attached: ( auto_phylo_diagram.tgz — includes all above, a python script, a demo makefile and IL-1 data ) <<<

In this post, I discuss a Python script I wrote that will convert a Newick Traversal and a FASTA file into a browser-viewable interactive diagram using SVG+JS+XHTML. This method is compatible with WebKit (Safari, Chrome) and Mozilla (Firefox) renderers but not Opera. XHTML doesn’t work with IE, so we’ll have to wait for massive adoption of HTML5 before I can make this work for everyone — don’t worry, it’s coming.

To the left is what the rendered tree looks like (shown as a screen captured PNG).

I recommend looking at the demo index.xhtml above if your browser supports it now. The demo will do a far better job than text in explaining exactly what features are possible using this method.

Try clicking on the different nodes and hyperlinks in the demo — notice that the displayed alignment strips away sites (columns) that are 100% gaps so that each node displays only a profile relevant for itself.

I originally intended to explain all the details of the JavaScripts that render and drive the diagram, but I think it’d be more useful to first focus on the Python script and the data input and output from the script.

The Attached Python Script to_xhtml.py found in auto_phylo_diagram.tgz

I’ll explain how to invoke to_xhtml.py with an example below.

python to_xhtml.py IL1fasta.txt IL1tree2.txt "IL-1 group from NCBI CDD (15)" > index.xhtml

Arguments …

  1. Plain Text FASTA file — IL1fasta.txt
  2. Plain Text Newick Traversal file — IL1tree2.txt
  3. Title of the generated XHTML — “IL-1 group from NCBI CDD (15)”

The data I used was generated with MUSCLE with default parameters. I specified the output FASTA with -out and the output second iteration tree with -tree2.

The Attached makefile found in auto_phylo_diagram.tgz

The makefile has two actual rules. If index.html does not exist, it will be regenerated with the above Python script and IL1fasta.txt and IL1tree2.txt. If either IL1fasta.txt or IL1tree2.txt or both do not exist, these files are regenerated using MUSCLE with default parameters on the included example data file ncbi.IL1.cl00094.15.fasta.txt. Running make on the included package shouldn’t do anything since all of the files are included.

The Example Data

The example data comes from the Interleukin-1 group (cd00100) from the NCBI Conserved Domains Database (CDD).

Finally, I’ll discuss briefly the nitty gritty of the JavaScripts found as stand-alones above that are also included in the package.

Part Two of Two

As is the nature of programatically drawing things, this particular task comes with its share of prerequisites.

In a previous — Part 1, I covered Interactive Diagrams Using Inline JS+SVG+XHTML — see that post for compatibility notes and the general method (remember, this technique is a holdover that will eventually be replaced with HTML5 and that the scripts in this post are not compatible with Opera anymore).

In this current post — Part 2, I’ll also assume you’ve seen Parsing a Newick Tree with recursive descent (you don’t need to know it inside and out, since we’ll only use these functions) — if you look at the example index.html attached, you’ll see that we actually parse the Newick Traversal using two functions in parse_newick_xhtml.jstree_builder() then traverse_tree(). The former converts the traversal to an in-memory binary tree while the latter accepts the root node of said tree and produces an array. This array is read in draw_phylogeny.js by fdraw_tree_svg() which does the actual rendering along with attaching all of the dynamic behaviour.

Parameters ( What are we drawing? )

To keep things simple, let’s focus on three requirements: (1) the tree will me sketched top-down with the root on top and its leaves below; (2) the diagram must make good use of the window width and avoid making horizontal scroll bars (unless it’s absolutely necessary); (3) the drawing function must make the diagram compatible with the overall document object model (DOM) so that we can highlight nodes and know what nodes were clicked to show the appropriate data.

Because of the second condition, all nodes in a given depth of the tree will be drawn at the same time. A horizontal repulsion will be applied so that each level will appear to be centre-justified.

General Strategy

Because we are centre-justifying the tree, we will be drawing each level of the tree simultaneously to calculate the offsets we want for each node. We perform this iteratively for each level of depth of the tree in the function fdraw_tree_svg(). Everything that is written dynamically by JavaScript goes to a div with a unique id. All of the SVG goes to id=”tree_display” and all of the textual output including raster images, text and alignments goes to id=”thoughts”. To attach behaviour to hyperlinks, whether it’s within the SVG or using standard anchor tags, the “onclick” property is defined. In this case, we always call the same function “onclick=javascript:select_node()“. This function accepts a single parameter: the integer serial number that has been clicked.

Additional Behaviour — Images and Links to Databases

After I wrote the first drafts of the JavaScripts, I decided to add two more behaviours in the current version. First, if a particular node name has the letters “GI” in it, the script will attempt to hyperlink it with the proceeding serial number against NCBI. Second, if a particular node name has the letters “PDB” in it, the script will attempt to hyperlink the following identifier against RSCB Protein Database and also pull the first thumbnail with an img tag for the 3D protein structure.

Enjoy!

This was done because I wanted to show a custom alignment algorithm to my advisors. In this demo, I’ve stripped away the custom alignment leaving behind a more general yet pleasing visualizer. I hope that the Python script and JavaScripts are useful to you. Enjoy!

Eddie Ma

January 1st, 2011 at 11:55 pm