Yet another bioinformatics blog

It has been a while since posting anything.  A long time back, I posted about Neural Network magic.  I finally got around to posting code implementing these tips and tricks for classification and ranking on github.  Check it out:  https://github.com/mkayala/nnutils

If you are pushing/pulling through ssh with git, git really does not like your shell to write anything to stdout.

Today, I was debugging a problem with the proper $PATH not getting set when submitting jobs to the Sun Grid Engine.  After debugging, I left an errant 'echo' sitting in my '~/.bashrc' file on the remote server.

20 minutes later, I try to push a recent local git commit to the same remote server that had the messed up '.bashrc', and I get this error:

$ git push fatal: protocol error: bad line length character: 006

A similar error occurs with a pull.  A quick google perusal shows that this is something that comes up often, but there aren't any well-explained solutions.  There are a number of bug-reports on the git mailing list.  

However, its simply a matter of making sure that your remote host doesn't write anything to stdout.  The simplest way to do this (assuming a bash login shell of course) is to just add this to the beginning of your ~/.bashrc:

if test "${PS1:-blah}" == "blah" ; then exit fi

Basically, don't do anything if this isn't an interactive shell (i.e, the $PS1 variable is not set).  I don't have this setup because I use my ~/.bashrc to dynamically load up environments depending on the machine I login to (where all machines I could log into have my home directory mounted on the same NFS server), and things like the Sun Grid Engine use non-interactive shells. So, in this case, I just make sure there are no echo debug statements hidden in there.

One project I work on involves using Django connected to a backend PostgreSQL database.  The glue between these two (specifically Python and PostgreSQL) is made up of the python Psycopg2 package.  

It turns out that there is a clash between the latest versions of Django (1.3) and Psycopg2 (2.4.2).  Here is the bug report with the Django devs.  The solution is "Don't use psycopg2 2.4.2, just use 2.4.1".  Simple enough.

It would be nice to work with everything in the standard Ubuntu repositories, but especially for dev packages, these often are either extremely out of date, or extremely up to date, but with no option of installing older versions.  However this was not that big of a deal to fix.  Simply uninstall Ubuntu's psycopg2 package, then install the PyPI installer: 

sudo apt-get install python-pip

Then install the specific pyscopg2 version with PyPI:

sudo pip install psycopg2==2.4.1

I broke down and decided to ditch my MacBook Pro recently.  I ended up buying a new Lenovo ThinkPad T420i. I bought the base model (with some minor exceptions).

Of course, the first thing I did was ditch Windows and install a flavor of linux.  Before getting the laptop, I tried all three main flavors of Ubuntu (straight Ubuntu, Kubuntu, and Xubuntu) in VirtualBox on the MacBook.  As recent updates to Ubuntu have a new strange desktop, I ended up picking Xubuntu.

For the most part, everything "Just Worked".  None of these are problems per se, but here are some notes to self (or others) about what I did to set up things to my liking. 

Installing Sun Java

For some reason, Chemaxon's java applets don't support OpenJDK Java, only Sun Java. Also, it seems like Ubuntu no longer supports Sun Java.  So following the instructions there, I had to:

sudo add-apt-repository ppa:ferramroberto/java sudo apt-get update sudo apt-get install sun-java6-jre sun-java6-plugin

And because I actually have OpenJDK java installed, I had to run this:

sudo update-alternatives --config java

Setting up dual monitors

While this didn't "just work", this is pretty easy.  I just run the following in a shell script whenever I plug/unplug my monitor.

xrandr --output VGA1 --auto --right-of LVDS1

To make this even easier, I set the key (Fn + F7) that does normal monitor switching to run this script.  This is done by opening Settings -> Settings Editor -> xfce4-keyboad-shortcuts.  Then I change the entry XF86Display to run this shell script.

Develop environment setup

This is a big week for me!  Not only did my first first-author journal paper just get accepted, but I'm also getting married on Friday!

We've put up a webserver to allow playing around with our machine learning based reaction predictor.  There is a pre-print of the journal paper on there as well:  http://reactions.ics.uci.edu/rxnpred/

I'm putting together a paper for NIPS 2011 using LaTeX.  Style guidelines for NIPS tell you to make the "References" header to be \subsubsection font size.  However, if you use Bibtex, this defaults to \section font size.  

After much googling I found this gem of a link.  It is a bit of a hack, but it worked like a charm!

About a year ago, I was exploring how to run quantum mechanics calculations to build descriptors for my reaction prediction work.  I ended up not using any of the calculations in a high-throughput manner, but it was an interesting exercise.  

I just stumbled across a nice little tutorial on using some open-source software to run QM calculations.  This uses Avogadro, GAMESS, and MacMolPlt on Windows.  When I was hacking around, I tried a number of different of configurations, but ultimately settled on GAMESS, pyMOL, and MacMolPlt.  I tried to get Avogadro to work as a pyMOL replacement, but there were (at least at the time) a number of bugs in the Mac/Linux versions.  It is totally possible that these have been fixed (I was on a pre-1.0 release, current = 1.01).  I would be interested in finding out because pyMOL was recently scooped up by Schrodinger, so although a version of it is still open-source, the future of the project is murky.  

1  Introduction

CyberT works by calculating a bayesian regularized t-statistic for each gene. This is a lot of language to basically say that we calculate a t-statistic for each gene, but we tweak the variance of each group by using genes with similar expression values.

The Mann-Whitney (or Wilcoxon Rank-Sum, or Mann-Whitney-Wilcoxon (MWW)) test is the non-parametric version of the t-test. Non-parametric in this case means that the statistic works on the ranks of the data, rather than assuming some parameterized underlying distribution. In contrast, proper analysis of the t-statistic is dependent on assuming that each group is from a normal distribution, i.e., is from a distribution parameterized by a mean and variance. The MWW statistic is very closely related to the area under the receiver operating characteristic curve, i.e., the area under the ROC curve, or AUC. In fact, if x is a vector of predictions on a set of positive examples, and y is a vector of predictions on a set of negative examples, then:

AUC(x,y) = MWW(x, y)/(|x| * |y|)

In this post, I'm going to demonstrate how to calculate the AUC for each gene in a microarray dataset. This basically says, if we used a single gene's expression value as a predictor for two different phenotypes, what would the AUC of this predictor be? Or another way to think of it is as a MWW statistic that is normalized to be between 0 and 1, or really between 0.5 and 1. One can also calculate p-values for these AUC's based on the p-values of the MWW stats.

2  The Function

Here is a function to do this calculation (in R of course). This actually calculates a one-vs-rest for all groups.

> ###################################################################### > #   calcAUCArray > # > #   Function to calculate the AUC by row for a matrix of predictions > #    > #   Input: > #       x           -   the data.frame or array  > #       repPerCond  -   a vector of the number per condition (assumes  > #                       the data.frame is in order). > ###################################################################### > calcAUCArray <- function(x, repPerCond, calcPVals = T){ +     classFactor <- factor(rep(1:length(repPerCond), repPerCond)) +      +     #The Auc is equiv to the Mann-Whitney or Wilcox test.   +     #For generalizations beyond 2 class, do a one-vs-rest scheme. +     classIndices <- tapply(1:ncol(x), classFactor, list) +      +     theRes <- sapply(classIndices, function(inds){ +         sapply(as.data.frame(t(x)), function(y){ +             maxW <- length(inds)*length(y[-inds]); +             theTest <-  wilcox.test(y[inds], y[-inds]); +             res <- theTest$statistic/maxW; +              +         }) +     }) +      +     ## This line flips the prediction (AUC < 0.5 is worse than random) +     theRes <- apply(theRes, 2, function(x){ifelse(x > 0.5, x, 1-x)}) +     colnames(theRes) <- paste("AUC_Cond_", 1:length(repPerCond), sep=''); +     if (calcPVals){ +         thePVals <- sapply(classIndices, function(inds){ +             sapply(as.data.frame(t(x)), function(y){ +                 theTest <-  wilcox.test(y[inds], y[-inds]); +                 theTest$p.value; +             }) +         }) +         colnames(thePVals) <- paste("AUC_PVal_Cond_", 1:length(repPerCond), sep=''); +         theRes <- cbind(theRes, thePVals) +     }   +     theRes; + }   

3  An Example

Here, let's fake up some data and show how the MWW could be useful in finding some genes of interest. First, make some fake completely normal data:

> data <- data.frame(matrix(rnorm(100), nrow=10)) > colnames(data) <- c(sprintf('Ctrl%d', 1:5),  +                     sprintf('Expt%d', 1:5)) > head(data)  

        Ctrl1       Ctrl2       Ctrl3      Ctrl4      Ctrl5      Expt1 1  1.05075025 -1.30088148  2.04442575  0.2396374 -0.1640639 0.46255951 2 -0.04274524  0.14674350  0.19292496 -3.0934266  1.3382232 1.44391446 3  0.32498054 -0.38788724 -1.18156093  0.4588947  0.4427220 1.74483569 4  1.61588825 -0.78319953 -0.07895495 -0.3512329  1.7633533 0.08020991 5 -0.15815873 -0.71218677 -0.54358175  1.1996652 -0.9436978 0.64061365 6 -1.68811312 -0.01021523  0.81717540  1.0148944 -1.6237105 1.39887543        Expt2      Expt3       Expt4        Expt5 1  0.1977868  1.2951139 -0.94278682 -0.458105214 2  1.1318220 -0.3654534  0.45353342 -0.514351083 3 -0.9689614 -0.7129198  0.07177276 -1.175522844 4 -1.7997450 -0.4597961  0.56408353  1.044657360 5 -0.4004646  0.5248907  0.49722154  0.009011202 6  0.5671171  0.9390662 -0.12485235  0.362181267  

Then, let's randomly pick a gene and just right shift the expression of just the Expt cols here by 2.

> fakeGeneId <- sample(1:10, 1) > fakeGeneId  

[1] 10  

> data[fakeGeneId, 6:10] <- data[fakeGeneId, 6:10] + 2;  

Now, let's calculate the AUCs and p-values here.

> repPerCond <- c(5,5) # 5 in each group > aucs <- calcAUCArray(data, repPerCond) > aucs  

      AUC_Cond_1 AUC_Cond_2 AUC_PVal_Cond_1 AUC_PVal_Cond_2 V1.W        0.56       0.56     0.841269841     0.841269841 V2.W        0.60       0.60     0.690476190     0.690476190 V3.W        0.60       0.60     0.690476190     0.690476190 V4.W        0.60       0.60     0.690476190     0.690476190 V5.W        0.76       0.76     0.222222222     0.222222222 V6.W        0.68       0.68     0.420634921     0.420634921 V7.W        0.52       0.52     1.000000000     1.000000000 V8.W        0.80       0.80     0.150793651     0.150793651 V9.W        0.56       0.56     0.841269841     0.841269841 V10.W       1.00       1.00     0.007936508     0.007936508  

> aucs[fakeGeneId,]  

     AUC_Cond_1      AUC_Cond_2 AUC_PVal_Cond_1 AUC_PVal_Cond_2      1.000000000     1.000000000     0.007936508     0.007936508   

There you go! AUCs and you should see a much lower p-value for the fakeGeneId!

File translated from TEX by TTH, version 3.85. On 7 Apr 2011, 17:54.

In the ML world, neural networks can be considered old hat, not sexy, etc, etc.  It's also common knowledge that it takes a magician to get a neural network to converge quickly, not overfit, etc.

Two things however:  

1.  My advisor is a neural networks guru, he made his name proving mathematical properties about them.  So, he likes them.   Again, he's my advisor, ergo I like them.

2.  With a few tricks, they actually work pretty well.  Not magic per se, just a couple of tricks.  

Here is a great tutorial about training neural networks.  Particularly useful are the sections about preconditioning the network and learning rate adaptation.  The one thing the tutorial lacks is a good reference to cite.  It took me a while to find it, but all of the ideas above seem to have been originally written up (or at least reviewed) here.

Yesterday was the first day of the ACS Meeting in Anaheim.  I was very interested in attending a symposium in honor of James Hendrickson.  Hendrickson was one of the pioneers of using computers in chemistry.  In particular, he was one of the first people to work on retro-synthesis, creating a program called SynGen.  The symposium was filled with talks about the history of this research in addition to current retro-synthesis work by other academics and software companies.  While, I don't work on retro-synthesis per se, this is very similar to what I work on.  

This day long talk of research close to what I do was absolutely great, especially because I presented a poster at the CINF Scholarship for Scientific Excellence poster session yesterday evening.  It was a great experience to have long, meaningful discussions with some of the leading researchers in the field.  I have business cards galore, several offers of collaboration, and many new ideas.  Here is a copy of my poster.  

Two types of files mess up my subversion working directories:  

*.pyc files - compiled python byte code files

*graffle files

OmniGraffle is my favorite diagramming tool for the Mac (and is tied with sheer laziness for why I don't just wipe my machine and use a real Linux distro).  It is great diagramming software, in that it's easy to use, the final images look great and are exportable as vector graphic, etc.  However, their design decisions about file storage pretty much suck.  *.graffle files are 'files' to Mac OSX, but they're really directories with a ton of other files hidden inside with data in non-standard formats that cmd line tools (like svn) choke on.  I am obsessive about version control.  So I commit everything (including papers I'm writing) to a subversion repo.  Unfortunately, without first archiving the graffle files, svn can't handle them.  So, I want svn to ignore them (and I version control an exported eps file).

*.pyc files are just created by the python interpreter when you run python scripts.  They are completely unnecessary to version control.

As these are global tweaks to subversion, the easiest way for me to 'ignore' these is to edit my subversion global-ignores in the config file. To do this:

emacs ~/.subversion/config

and edit the line starting with global-ignores to have something like:

global-ignores = *.o *.lo *.la #*# .*.rej *.rej .*~ *~ .#* .DS_Store *graffle* *.pyc

Machine learning and LaTeX, what's not to like? Check out the hand drawn LaTeX symbol recognizer.