Jermdemo Raised to the Law

AGBT: digesting diposable MinIONs in diaspora

2012-02-20T21:44:00.000-05:00

Despite my current ranking of 15th in Biostar, myriad page views of my BAS™ post (albeit mostly misdirected perverts), and positive response for my celebrated campaign against more microarray papers, for some reason I was not "comped" an all-expenses paid trip as honorary blog journalist to this year's Advances in Biology and Genome Technology, which is kind of like CES for sequencing people, except AGBT is still worth attending. Normally the oversight would not bother me, as bioinformatics itself is not the focus of this meeting, but the flood of #AGBT tweets would not let me forget this fact and I was forced to stew and blog in envy.

The first game changing disruptive revolutionary thing from England since 1964

Even from my distant perch it was obvious all the scientific presentations at AGBT were overshadowed by a 17-minute showstopping demo from Clive Brown of Oxford Nanopore, a company that by all appearances would either die, focus on some minor stuff, or bring it. They chose the third option, and in so doing boosted the "Clive index" to unprecedented levels. OxN's recent decision to enlist famed geneticist and serial startup advisor George Church struck me as a huge gamble, as the string of Route 128 flameouts touting his name lead me to assume long ago that Church had stowed away some cursed Tiki idol in his luggage like Bobby in that episode of the Brady Bunch. However, after reading up on OxN, I had to admit I was just bitter about Dr. Church's refusal to invest in my chain of Polonator-based paternity testing clinics, Yo Po'lonatizz!™

Two new sequencer platforms were announced:

A MinION. Forget to hit eject before removing this
and you will instantly lose $900.

The MinION, a $900 "disposable" USB drive which detects minute changes in voltage incurred by the passage of DNA through a robust and delicious lipid bilayer. Finally a device capable of sequencing filthy rabbit blood right on the spot!
The GridION system, a scalable rackmounted sequencer, which despite some lack of pricing clarity, should produce an actual $1000 15-minute human genome by 2013.

These exotic machines must be truly game-changing because they made properly expanding Albert Vilella's NGS sequencer spreadsheet quite difficult. The MinION, in particular, could be viewed as a free device with $900 of consumables. This effectively lowers the bar to getting high-throughput sequence in the doctor's office to a 100% unamortized billable transaction. These things also claim fucking unlimited read lengths.

Expression microarrays, SAGE, 454, ABI SOLiD, and now Pacific Biosciences have all left bad tastes of uncertainty and dissatisfaction in the mouths of scientists. It is easy to disappoint people on a grand scale with a $700,000 machine, but $900 worth of chemicals in a USB drive is a different animal, and it seems likely this invention will find a following if it even delivers on a fraction of what it promises.

The GridION - put it in a rack or right on the floor.

Good information on this sequencer-on-a-stick is to be found at Nick Loman's blog, Genomes Unzipped, and official press releases. An excellent discussion of the nanopores themselves can be found at Omically Speaking.

More cringeworthy marketing from the West coast

The Oxford Nanopore machines are so jaw-dropping, in fact, that Jonathan Rothberg is already crying vaporware. His complaints do seem warranted, given disappointments from past year's announcements and the lack of publicly available sequence from these devices.

Unfortunately Ion Torrent has spent all of its goodwill on an inane and hamfisted advertising war against Illumina's MiSeq, an intentionally crippled opponent. Seemingly orchestrated by castoffs from the Celebrity Apprentice, this assault began with cringe-inducing derivations of Apple commercials, and has expanded to include a sort of "feature combover." Through some convoluted logic involving consensus, a professional whiteboard artist attempts to convince the public how the homopolymer error rate is actually lower using Ion Torrent PGM than MiSeq. This is the sequencing equivalent of having your mom try to convince you two apples is better than one devil dog, or some such utter nonsense.

My response was predictably measured and cerebral:

This is not the first time I have tweet-confronted Ion Torrent over its odious approach. All this is rather unnecessary because overall, and despite the homopolymer issues, the utility of the PGM has been more or less within expectations. The MiSeq is also exactly within expectations, since it is basically a transparent, measly 1/50th slice of a HiSeq. The same cannot really be said for the RS, whose error rate is clearly far above what was expected at the outset. So if anyone requires an aggressive smokescreen-type marketing campaign (or a new machine) it is Pacific Biosciences.

When can we expect the last damn microarray paper?

2012-01-18T21:15:00.000-05:00

With bonus R code

It came as a shock to learn from PubMed that almost 900 papers were published with the word "microarray" in their titles last year alone, just 12 shy of the 2010 count. More alarming, many of these papers were not of the innocuous "Microarray study of gene expression in dog scrotal tissue" variety, but dry rehashings along the lines of "Statistical approaches to normalizing microarrays to the reference brightness of Ursa Minor".

It's an ugly truth we must face: people aren't just using microarrays, they're still writing about them.

See for yourself:

getCount<-function(term){function(year){
  nihUrl<-concat("http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=pubmed&term=",term,"+",year,"[pdat]")
  #cleanurl<-gsub('\\]','%5D',gsub('\\[','%5B',x=url))
  #http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=pubmed&term=microarray%5btitle%5d+2003%5bpdat%5d
  xml<-xmlTreeParse(URLencode(nihUrl),isURL=TRUE)
  #Data Mashups in R, pg17
  as.numeric(xmlValue(xml$doc$children$eSearchResult$children$Count$children$text))
}}

years<-1995:2011
df<-data.frame(type="obs",year=years,
    mic=sapply(years,function(x){do.call(getCount('microarray[title]'),list(x))}),
    ngs=sapply(years,function(x){do.call(getCount('"next generation sequencing"[title] OR "high-throughput sequencing"[title]'),list(x))})
)
#papers with "microarray" in title
> df[,c("year","mic")]
   year  mic
1  1995    2
2  1996    4
3  1997    0
4  1998    7
5  1999   28
6  2000  108
7  2001  273
8  2002  553
9  2003  770
10 2004 1032
11 2005 1135
12 2006 1216
13 2007 1107
14 2008 1055
15 2009  981
16 2010  909
17 2011  897

Reading another treatise on microarray normalization in 2012 would be just tragic. Who still reads these? Who still writes these papers? Can we stop them? If not, when can we expect NGS to wipe them off the map?

#97 is a fair start
df<-subset(df,year>=1997)
mdf<-melt(df,id.vars=c("type","year"),variable_name="citation")

c<-ggplot(mdf,aes(x=year))
p<-c+geom_point(aes(y=value,color=citation)) +
  ylab("papers") +
  stat_smooth(aes(y=value,color=citation),data=subset(mdf,citation=="mic"),method="loess") +
  scale_x_continuous(breaks=seq(from=1997,to=2011,by=2))
print(p)

Here I plot both microarray and next-generation sequencing papers (in title). We see kurtosis is working in our favor, and LOESS seems to agree!

But when will the pain end? Let us extrapolate, wildly.

#Return 0 for negative elements
# noNeg(c(3,2,1,0,-1,-2,2))
# [1] 3 2 1 0 0 0 2
noNeg<-function(v){sapply(v,function(x){max(x,0)})}

#Return up to the first negative/zero element inclusive
# toZeroNoNeg(c(3,2,1,0,-1,-2,2))
# [1] 3 2 1 0
toZeroNoNeg<-function(v){noNeg(v)[1:firstZero(noNeg(v))]}

#return index of first zero
firstZero<-function(v){which(noNeg(v)==0)[1]}

#let's peer into the future
df.lo.mic<-loess(mic ~ year,df,control=loess.control(surface="direct"))

#when will it stop?
mic_predict<-as.integer(predict(df.lo.mic,data.frame(year=2012:2020),se=FALSE))
zero_year<-2011+firstZero(mic_predict)
cat(concat("LOESS projects ",sum(toZeroNoNeg(mic_predict))," more microarray papers."))
cat(concat("The last damn microarray paper is projected to be in ",(zero_year-1),"."))

#predict ngs growth
df.lo.ngs<-loess(ngs ~ year,df,control=loess.control(surface="direct"))
ngs_predict<-as.integer(predict(df.lo.ngs,data.frame(year=2012:zero_year),se=FALSE))

pred_df<-data.frame(type="pred",year=c(2012:zero_year),mic=toZeroNoNeg(mic_predict),ngs=ngs_predict)
df2<-rbind(df,pred_df)

mdf2<-melt(df2,id.vars=c("type","year"),variable_name="citation")

c2<-ggplot(mdf2,aes(x=year))
p2<-c2+geom_point(aes(y=value,color=citation,shape=type),size=3) +
    ylab("papers") +
    scale_y_continuous(breaks=seq(from=0,to=1600,by=200))+
    scale_x_continuous(breaks=seq(from=1997,to=zero_year,by=2))
print(p2)

LOESS projects 2038 more microarray papers.
The last damn microarray paper is projected to be published in 2016.

Yeah, right...

Full R code here: https://gist.github.com/1637248

Making R's paste act more like CONCAT

2011-10-31T13:17:00.003-04:00

While vector-friendly, R's paste function has a few behaviors I don't particularly like.

One is using a space as the default separator:

> adjectives<-c("lean","fast","strong")
> paste(adjectives,"er")
> paste(adjectives,"er")
[1] "lean er"   "fast er"   "strong er"  #d'oh
> paste(adjectives,"er",sep="")
[1] "leaner"   "faster"   "stronger"

Empty vectors get an undeserved first class treatment:

> paste(indelPositions,"i",sep="")
[1] "i"
> indelPositions<-c(5,6,7)
> paste(indelPositions,"i",sep="")
[1] "5i" "6i" "7i" #good

> indelPositions<-c()
> paste(indelPositions,"i",sep="")
[1] "i"  #not so good

And perhaps worst of all, NA values get replaced with a string called "NA":

> placing<-"1"
> paste(placing,"st",sep="")
[1] "1st" #awesome

> placing<-NA_integer_
> paste(placing,"st",sep="")
[1] "NAst" #ugh

This is inconvenient in situations where I don't know a priori if I will get a value, a vector of length 0, or an NA.

Working from Hadley Wickham's str_c function in the stringr package, I decided to write a paste function that behaves more like CONCAT in SQL:

library(stringr)
concat<-CONCAT<-function(...,sep="",collapse=NULL){
  strings<-list(...)
  #catch NULLs, NAs
  if(
    all(unlist(llply(strings,length))>0)
    &&
    all(!is.na(unlist(strings)))
    ){
    do.call("paste", c(strings, list(sep = sep, collapse = collapse)))
  }else{
    NULL
  }
}

This function has the behaviors I expect:

> concat(adjectives,"er")
[1] "leaner"   "faster"   "stronger"

> concat(indelPositions,"i")
NULL

> concat(placing,"st")
NULL

That's more like it!

SELinux for enhanced headaches

2011-10-06T13:10:00.000-04:00

Security Enhanced Linux (SELinux) is a new extra hidden layer of permissions that makes configuring things more difficult, without ever identifying itself as the culprit - kind of like ACLs but more cryptic. Though it may be more secure, it is not an enhancing experience to deal with, and probably not worth it for the average user.

For example to have Apache serve personal websites (i.e. http://server/~leipzig) it is no longer enough to alter httpd.conf, because you will get mysterious 403 errors until you do this (as others have experienced):

chcon -R -t httpd_sys_content_t /home/leipzig

You forget about this change until xauth starts complaining about stuff for no apparent reason:

/usr/bin/xauth:  timeout in locking authority file /home/leipzig/.Xauthority

so of course you need to do this (thanks Madhav Diwan for this post):

chcon unconfined_u:object_r:user_home_dir_t:s0 /home/leipzig

I have no idea what these things actually mean, nor any real interest in learning. I'm sure this stuff is great for sysadmin cocktail chat but at least for private servers it is just another the brake on the wheel of getting things done. For the time being I have set the level to "permissive", which means it displays warnings but does not interfere, but am leaning toward "disabled" or maybe something else:

# This file controls the state of SELinux on the system.
# SELINUX= can take one of these three values:
#     enforcing - SELinux security policy is enforced.
#     permissive - SELinux prints warnings instead of enforcing.
#     disabled - No SELinux policy is loaded.
SELINUX=excoriated
# SELINUXTYPE= can take one of these two values:
#     targeted - Targeted processes are protected,
#     mls - Multi Level Security protection.
SELINUXTYPE=targeted

More on the pros and cons:
http://unix.stackexchange.com/questions/9163/does-selinux-provide-enough-extra-security-to-be-worth-the-hassle-of-learning-set

Installing RStudio Server on Scientific Linux 6: My bash notebook

2011-08-18T10:03:00.008-04:00

Granted, not a brilliant sysadmin mind at work here, but this might help someone someday.
Scientific Linux (SL) is built from Red Hat Enterprise Linux

See installation instructions here:
http://rstudio.org/download/server

[leipzig@localhost ~]$ sudo rpm -Uvh
http://download.fedoraproject.org/pub/epel/6/x86_64/epel-release-6-5.noarch.rpm
[sudo] password for leipzig: 
Retrieving
http://download.fedoraproject.org/pub/epel/6/x86_64/epel-release-6-5.noarch.rpm
warning: /var/tmp/rpm-tmp.S2RQAH: Header V3 RSA/SHA256 Signature, key ID
0608b895: NOKEY
Preparing...                ########################################### [100%]
   1:epel-release           ########################################### [100%]
[leipzig@localhost ~]$ rpm -qa | grep epel
epel-release-6-5.noarch
[leipzig@localhost ~]$ which R
/usr/local/bin/R
[leipzig@localhost ~]$ R

R version 2.13.0 (2011-04-13)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-unknown-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

&gt; q()
Save workspace image? [y/n/c]: n
[leipzig@localhost ~]$ wget
https://s3.amazonaws.com/rstudio-server/rstudio-server-0.94.92-x86_64.rpm
--2011-08-17 13:06:36-- 
https://s3.amazonaws.com/rstudio-server/rstudio-server-0.94.92-x86_64.rpm
Resolving s3.amazonaws.com... 72.21.211.170
Connecting to s3.amazonaws.com|72.21.211.170|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 11373769 (11M) [application/x-redhat-package-manager]
Saving to: “rstudio-server-0.94.92-x86_64.rpm”

100%[===========================================================================
=============================================&gt;] 11,373,769  7.89M/s   in 1.4s   


2011-08-17 13:06:37 (7.89 MB/s) - “rstudio-server-0.94.92-x86_64.rpm” saved
[11373769/11373769]

[leipzig@localhost ~]$ sudo rpm -Uvh rstudio-server-0.94.92-x86_64.rpm
error: Failed dependencies:
	libR.so()(64bit) is needed by rstudio-server-0.94.92-1.x86_64
	libRblas.so()(64bit) is needed by rstudio-server-0.94.92-1.x86_64
	libRlapack.so()(64bit) is needed by rstudio-server-0.94.92-1.x86_64
	libcrypto.so.6()(64bit) is needed by rstudio-server-0.94.92-1.x86_64
	libgfortran.so.1()(64bit) is needed by rstudio-server-0.94.92-1.x86_64
	libssl.so.6()(64bit) is needed by rstudio-server-0.94.92-1.x86_64
[leipzig@localhost ~]$ sudo yum install R  
epel/metalink                                                                   
                                                          |  14 kB     00:00    

epel                                                                            
                                                          | 4.3 kB     00:00    

epel/primary_db                                                                 
                                                          | 4.0 MB     00:00    

sl                                                                              
                                                          | 3.2 kB     00:00    

sl-security                                                                     
                                                          | 1.9 kB     00:00    

Setting up Install Process
Resolving Dependencies
--&gt; Running transaction check
---&gt; Package R.x86_64 0:2.13.1-1.el6 set to be updated
--&gt; Processing Dependency: libRmath-devel = 2.13.1-1.el6 for package:
R-2.13.1-1.el6.x86_64
--&gt; Processing Dependency: R-devel = 2.13.1-1.el6 for package:
R-2.13.1-1.el6.x86_64
--&gt; Running transaction check
---&gt; Package R-devel.x86_64 0:2.13.1-1.el6 set to be updated
--&gt; Processing Dependency: R-core = 2.13.1-1.el6 for package:
R-devel-2.13.1-1.el6.x86_64
--&gt; Processing Dependency: bzip2-devel for package: R-devel-2.13.1-1.el6.x86_64
--&gt; Processing Dependency: gcc-gfortran for package: R-devel-2.13.1-1.el6.x86_64
--&gt; Processing Dependency: tk-devel for package: R-devel-2.13.1-1.el6.x86_64
--&gt; Processing Dependency: pcre-devel for package: R-devel-2.13.1-1.el6.x86_64
--&gt; Processing Dependency: tcl-devel for package: R-devel-2.13.1-1.el6.x86_64
---&gt; Package libRmath-devel.x86_64 0:2.13.1-1.el6 set to be updated
--&gt; Processing Dependency: libRmath = 2.13.1-1.el6 for package:
libRmath-devel-2.13.1-1.el6.x86_64
--&gt; Running transaction check
---&gt; Package R-core.x86_64 0:2.13.1-1.el6 set to be updated
--&gt; Processing Dependency: cups for package: R-core-2.13.1-1.el6.x86_64
--&gt; Processing Dependency: libtk8.5.so()(64bit) for package:
R-core-2.13.1-1.el6.x86_64
---&gt; Package bzip2-devel.x86_64 0:1.0.5-7.el6_0 set to be updated
---&gt; Package gcc-gfortran.x86_64 0:4.4.4-13.el6 set to be updated
---&gt; Package libRmath.x86_64 0:2.13.1-1.el6 set to be updated
---&gt; Package pcre-devel.x86_64 0:7.8-3.1.el6 set to be updated
---&gt; Package tcl-devel.x86_64 1:8.5.7-6.el6 set to be updated
---&gt; Package tk-devel.x86_64 1:8.5.7-5.el6 set to be updated
--&gt; Running transaction check
---&gt; Package cups.x86_64 1:1.4.2-35.el6_0.1 set to be updated
--&gt; Processing Dependency: portreserve for package:
1:cups-1.4.2-35.el6_0.1.x86_64
--&gt; Processing Dependency: poppler-utils for package:
1:cups-1.4.2-35.el6_0.1.x86_64
---&gt; Package tk.x86_64 1:8.5.7-5.el6 set to be updated
--&gt; Running transaction check
---&gt; Package poppler-utils.x86_64 0:0.12.4-3.el6_0.1 set to be updated
---&gt; Package portreserve.x86_64 0:0.0.4-4.el6 set to be updated
--&gt; Finished Dependency Resolution

Dependencies Resolved

================================================================================
================================================================================
==
 Package                                 Arch                            Version
                                      Repository                            Size
================================================================================
================================================================================
==
Installing:
 R                                       x86_64                         
2.13.1-1.el6                                  epel                              
   17 k
Installing for dependencies:
 R-core                                  x86_64                         
2.13.1-1.el6                                  epel                              
   33 M
 R-devel                                 x86_64                         
2.13.1-1.el6                                  epel                              
   88 k
 bzip2-devel                             x86_64                         
1.0.5-7.el6_0                                 sl-security                       
  249 k
 cups                                    x86_64                         
1:1.4.2-35.el6_0.1                            sl-security                       
  2.3 M
 gcc-gfortran                            x86_64                         
4.4.4-13.el6                                  sl                                
  4.7 M
 libRmath                                x86_64                         
2.13.1-1.el6                                  epel                              
  111 k
 libRmath-devel                          x86_64                         
2.13.1-1.el6                                  epel                              
   21 k
 pcre-devel                              x86_64                         
7.8-3.1.el6                                   sl                                
  317 k
 poppler-utils                           x86_64                         
0.12.4-3.el6_0.1                              sl-security                       
   72 k
 portreserve                             x86_64                         
0.0.4-4.el6                                   sl                                
   21 k
 tcl-devel                               x86_64                         
1:8.5.7-6.el6                                 sl                                
  161 k
 tk                                      x86_64                         
1:8.5.7-5.el6                                 sl                                
  1.4 M
 tk-devel                                x86_64                         
1:8.5.7-5.el6                                 sl                                
  495 k

Transaction Summary
================================================================================
================================================================================
==
Install      14 Package(s)
Upgrade       0 Package(s)

Total download size: 43 M
Installed size: 89 M
Is this ok [y/N]: y
Downloading Packages:
(1/14): R-2.13.1-1.el6.x86_64.rpm                                               
                                                          |  17 kB     00:00    

(2/14): R-core-2.13.1-1.el6.x86_64.rpm                                          
                                                          |  33 MB     00:05    

(3/14): R-devel-2.13.1-1.el6.x86_64.rpm                                         
                                                          |  88 kB     00:00    

(4/14): bzip2-devel-1.0.5-7.el6_0.x86_64.rpm                                    
                                                          | 249 kB     00:00    

(5/14): cups-1.4.2-35.el6_0.1.x86_64.rpm                                        
                                                          | 2.3 MB     00:01    

(6/14): gcc-gfortran-4.4.4-13.el6.x86_64.rpm                                    
                                                          | 4.7 MB     00:02    

(7/14): libRmath-2.13.1-1.el6.x86_64.rpm                                        
                                                          | 111 kB     00:00    

(8/14): libRmath-devel-2.13.1-1.el6.x86_64.rpm                                  
                                                          |  21 kB     00:00    

(9/14): pcre-devel-7.8-3.1.el6.x86_64.rpm                                       
                                                          | 317 kB     00:00    

(10/14): poppler-utils-0.12.4-3.el6_0.1.x86_64.rpm                              
                                                          |  72 kB     00:00    

(11/14): portreserve-0.0.4-4.el6.x86_64.rpm                                     
                                                          |  21 kB     00:00    

(12/14): tcl-devel-8.5.7-6.el6.x86_64.rpm                                       
                                                          | 161 kB     00:00    

(13/14): tk-8.5.7-5.el6.x86_64.rpm                                              
                                                          | 1.4 MB     00:00    

(14/14): tk-devel-8.5.7-5.el6.x86_64.rpm                                        
                                                          | 495 kB     00:00    

--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
--
Total                                                                           
                                                 3.1 MB/s |  43 MB     00:13    

warning: rpmts_HdrFromFdno: Header V3 RSA/SHA256 Signature, key ID 0608b895:
NOKEY
epel/gpgkey                                                                     
                                                          | 3.2 kB     00:00 ...

Importing GPG key 0x0608B895 "EPEL (6) epel@fedoraproject.org" from
/etc/pki/rpm-gpg/RPM-GPG-KEY-EPEL-6
Is this ok [y/N]: y
Running rpm_check_debug
Running Transaction Test
Transaction Test Succeeded
Running Transaction
Warning: RPMDB altered outside of yum.
  Installing     : 1:tk-8.5.7-5.el6.x86_64                                      
                                                                            1/14

  Installing     : portreserve-0.0.4-4.el6.x86_64                               
                                                                            2/14

  Installing     : poppler-utils-0.12.4-3.el6_0.1.x86_64                        
                                                                            3/14

  Installing     : 1:cups-1.4.2-35.el6_0.1.x86_64                               
                                                                            4/14

  Installing     : R-core-2.13.1-1.el6.x86_64                                   
                                                                            5/14

  Installing     : gcc-gfortran-4.4.4-13.el6.x86_64                             
                                                                            6/14

  Installing     : libRmath-2.13.1-1.el6.x86_64                                 
                                                                            7/14

  Installing     : 1:tcl-devel-8.5.7-6.el6.x86_64                               
                                                                            8/14

  Installing     : 1:tk-devel-8.5.7-5.el6.x86_64                                
                                                                            9/14

  Installing     : libRmath-devel-2.13.1-1.el6.x86_64                           
                                                                           10/14

  Installing     : bzip2-devel-1.0.5-7.el6_0.x86_64                             
                                                                           11/14

  Installing     : pcre-devel-7.8-3.1.el6.x86_64                                
                                                                           12/14

  Installing     : R-devel-2.13.1-1.el6.x86_64                                  
                                                                           13/14

  Installing     : R-2.13.1-1.el6.x86_64                                        
                                                                           14/14


Installed:
  R.x86_64 0:2.13.1-1.el6                                                       
                                                                                


Dependency Installed:
  R-core.x86_64 0:2.13.1-1.el6                R-devel.x86_64 0:2.13.1-1.el6     
  bzip2-devel.x86_64 0:1.0.5-7.el6_0       cups.x86_64 1:1.4.2-35.el6_0.1     
  gcc-gfortran.x86_64 0:4.4.4-13.el6          libRmath.x86_64 0:2.13.1-1.el6    
  libRmath-devel.x86_64 0:2.13.1-1.el6     pcre-devel.x86_64 0:7.8-3.1.el6    
  poppler-utils.x86_64 0:0.12.4-3.el6_0.1     portreserve.x86_64 0:0.0.4-4.el6  
  tcl-devel.x86_64 1:8.5.7-6.el6           tk.x86_64 1:8.5.7-5.el6            
  tk-devel.x86_64 1:8.5.7-5.el6              

Complete!
[leipzig@localhost ~]$ sudo rpm -Uvh rstudio-server-0.94.92-x86_64.rpm
error: Failed dependencies:
	libcrypto.so.6()(64bit) is needed by rstudio-server-0.94.92-1.x86_64
	libgfortran.so.1()(64bit) is needed by rstudio-server-0.94.92-1.x86_64
	libssl.so.6()(64bit) is needed by rstudio-server-0.94.92-1.x86_64
[leipzig@localhost ~]$ sudo yum install libcrypto.so.6
Setting up Install Process
Resolving Dependencies
--&gt; Running transaction check
---&gt; Package openssl098e.i686 0:0.9.8e-17.el6 set to be updated
--&gt; Processing Dependency: libc.so.6(GLIBC_2.3.4) for package:
openssl098e-0.9.8e-17.el6.i686
--&gt; Processing Dependency: libkrb5.so.3(krb5_3_MIT) for package:
openssl098e-0.9.8e-17.el6.i686
--&gt; Processing Dependency: libc.so.6(GLIBC_2.1) for package:
openssl098e-0.9.8e-17.el6.i686
--&gt; Processing Dependency: libcom_err.so.2 for package:
openssl098e-0.9.8e-17.el6.i686
--&gt; Processing Dependency: libc.so.6(GLIBC_2.0) for package:
openssl098e-0.9.8e-17.el6.i686
--&gt; Processing Dependency: libk5crypto.so.3 for package:
openssl098e-0.9.8e-17.el6.i686
--&gt; Processing Dependency: libk5crypto.so.3(k5crypto_3_MIT) for package:
openssl098e-0.9.8e-17.el6.i686
--&gt; Processing Dependency: libdl.so.2(GLIBC_2.0) for package:
openssl098e-0.9.8e-17.el6.i686
--&gt; Processing Dependency: libc.so.6(GLIBC_2.7) for package:
openssl098e-0.9.8e-17.el6.i686
--&gt; Processing Dependency: libkrb5.so.3 for package:
openssl098e-0.9.8e-17.el6.i686
--&gt; Processing Dependency: libc.so.6(GLIBC_2.4) for package:
openssl098e-0.9.8e-17.el6.i686
--&gt; Processing Dependency: libgssapi_krb5.so.2 for package:
openssl098e-0.9.8e-17.el6.i686
--&gt; Processing Dependency: libdl.so.2(GLIBC_2.1) for package:
openssl098e-0.9.8e-17.el6.i686
--&gt; Processing Dependency: libc.so.6(GLIBC_2.1.3) for package:
openssl098e-0.9.8e-17.el6.i686
--&gt; Processing Dependency: libresolv.so.2 for package:
openssl098e-0.9.8e-17.el6.i686
--&gt; Processing Dependency: libz.so.1 for package: openssl098e-0.9.8e-17.el6.i686
--&gt; Processing Dependency: libc.so.6 for package: openssl098e-0.9.8e-17.el6.i686
--&gt; Processing Dependency: libdl.so.2 for package:
openssl098e-0.9.8e-17.el6.i686
--&gt; Processing Dependency: libc.so.6(GLIBC_2.3) for package:
openssl098e-0.9.8e-17.el6.i686
--&gt; Running transaction check
---&gt; Package glibc.i686 0:2.12-1.7.el6_0.5 set to be updated
--&gt; Processing Dependency: libfreebl3.so for package:
glibc-2.12-1.7.el6_0.5.i686
--&gt; Processing Dependency: libfreebl3.so(NSSRAWHASH_3.12.3) for package:
glibc-2.12-1.7.el6_0.5.i686
---&gt; Package krb5-libs.i686 0:1.9-9.el6_1.1 set to be updated
--&gt; Processing Dependency: libkeyutils.so.1(KEYUTILS_0.3) for package:
krb5-libs-1.9-9.el6_1.1.i686
--&gt; Processing Dependency: libkeyutils.so.1 for package:
krb5-libs-1.9-9.el6_1.1.i686
--&gt; Processing Dependency: libselinux.so.1 for package:
krb5-libs-1.9-9.el6_1.1.i686
---&gt; Package libcom_err.i686 0:1.41.12-3.el6 set to be updated
---&gt; Package zlib.i686 0:1.2.3-25.el6 set to be updated
--&gt; Running transaction check
---&gt; Package keyutils-libs.i686 0:1.4-1.el6 set to be updated
---&gt; Package libselinux.i686 0:2.0.94-2.el6 set to be updated
---&gt; Package nss-softokn-freebl.i686 0:3.12.8-1.el6_0 set to be updated
--&gt; Finished Dependency Resolution

Dependencies Resolved

================================================================================
================================================================================
==
 Package                                     Arch                         
Version                                     Repository                          
 Size
================================================================================
================================================================================
==
Installing:
 openssl098e                                 i686                         
0.9.8e-17.el6                               sl                                  
772 k
Installing for dependencies:
 glibc                                       i686                         
2.12-1.7.el6_0.5                            sl-security                         
4.3 M
 keyutils-libs                               i686                         
1.4-1.el6                                   sl                                  
 19 k
 krb5-libs                                   i686                         
1.9-9.el6_1.1                               sl-security                         
711 k
 libcom_err                                  i686                         
1.41.12-3.el6                               sl                                  
 33 k
 libselinux                                  i686                         
2.0.94-2.el6                                sl                                  
106 k
 nss-softokn-freebl                          i686                         
3.12.8-1.el6_0                              sl-security                         
108 k
 zlib                                        i686                         
1.2.3-25.el6                                sl                                  
 71 k

Transaction Summary
================================================================================
================================================================================
==
Install       8 Package(s)
Upgrade       0 Package(s)

Total download size: 6.0 M
Installed size: 18 M
Is this ok [y/N]: y
Downloading Packages:
(1/8): glibc-2.12-1.7.el6_0.5.i686.rpm                                          
                                                          | 4.3 MB     00:02    

(2/8): keyutils-libs-1.4-1.el6.i686.rpm                                         
                                                          |  19 kB     00:00    

(3/8): krb5-libs-1.9-9.el6_1.1.i686.rpm                                         
                                                          | 711 kB     00:00    

(4/8): libcom_err-1.41.12-3.el6.i686.rpm                                        
                                                          |  33 kB     00:00    

(5/8): libselinux-2.0.94-2.el6.i686.rpm                                         
                                                          | 106 kB     00:00    

(6/8): nss-softokn-freebl-3.12.8-1.el6_0.i686.rpm                               
                                                          | 108 kB     00:00    

(7/8): openssl098e-0.9.8e-17.el6.i686.rpm                                       
                                                          | 772 kB     00:00    

(8/8): zlib-1.2.3-25.el6.i686.rpm                                               
                                                          |  71 kB     00:00    

--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
--
Total                                                                           
                                                 1.2 MB/s | 6.0 MB     00:04    

Running rpm_check_debug
Running Transaction Test
Transaction Test Succeeded
Running Transaction
  Installing     : nss-softokn-freebl-3.12.8-1.el6_0.i686                       
                                                                             1/8

  Installing     : glibc-2.12-1.7.el6_0.5.i686                                  
                                                                             2/8

  Installing     : libcom_err-1.41.12-3.el6.i686                                
                                                                             3/8

  Installing     : zlib-1.2.3-25.el6.i686                                       
                                                                             4/8

  Installing     : libselinux-2.0.94-2.el6.i686                                 
                                                                             5/8

  Installing     : keyutils-libs-1.4-1.el6.i686                                 
                                                                             6/8

  Installing     : krb5-libs-1.9-9.el6_1.1.i686                                 
                                                                             7/8

  Installing     : openssl098e-0.9.8e-17.el6.i686                               
                                                                             8/8


Installed:
  openssl098e.i686 0:0.9.8e-17.el6                                              
                                                                                


Dependency Installed:
  glibc.i686 0:2.12-1.7.el6_0.5        keyutils-libs.i686 0:1.4-1.el6           
     krb5-libs.i686 0:1.9-9.el6_1.1       libcom_err.i686 0:1.41.12-3.el6      
  libselinux.i686 0:2.0.94-2.el6       nss-softokn-freebl.i686 0:3.12.8-1.el6_0 
     zlib.i686 0:1.2.3-25.el6            

Complete!
[leipzig@localhost ~]$ sudo rpm -Uvh rstudio-server-0.94.92-x86_64.rpm
error: Failed dependencies:
	libcrypto.so.6()(64bit) is needed by rstudio-server-0.94.92-1.x86_64
	libgfortran.so.1()(64bit) is needed by rstudio-server-0.94.92-1.x86_64
	libssl.so.6()(64bit) is needed by rstudio-server-0.94.92-1.x86_64
[leipzig@localhost ~]$ sudo yum install libcrypto.so.6
Setting up Install Process
Package openssl098e-0.9.8e-17.el6.i686 already installed and latest version
Nothing to do
[leipzig@localhost ~]$ sudo yum install libgfortran.so.1
Setting up Install Process
Resolving Dependencies
--&gt; Running transaction check
---&gt; Package compat-libgfortran-41.i686 0:4.1.2-39.el6 set to be updated
--&gt; Finished Dependency Resolution

Dependencies Resolved

================================================================================
================================================================================
==
 Package                                           Arch                         
   Version                                  Repository                      Size
================================================================================
================================================================================
==
Installing:
 compat-libgfortran-41                             i686                         
   4.1.2-39.el6                             sl                              99 k

Transaction Summary
================================================================================
================================================================================
==
Install       1 Package(s)
Upgrade       0 Package(s)

Total download size: 99 k
Installed size: 488 k
Is this ok [y/N]: y
Downloading Packages:
compat-libgfortran-41-4.1.2-39.el6.i686.rpm                                     
                                                          |  99 kB     00:00    

Running rpm_check_debug
Running Transaction Test
Transaction Test Succeeded
Running Transaction
  Installing     : compat-libgfortran-41-4.1.2-39.el6.i686                      
                                                                             1/1


Installed:
  compat-libgfortran-41.i686 0:4.1.2-39.el6                                     
                                                                                


Complete!
[leipzig@localhost ~]$ sudo rpm -Uvh rstudio-server-0.94.92-x86_64.rpm
error: Failed dependencies:
	libcrypto.so.6()(64bit) is needed by rstudio-server-0.94.92-1.x86_64
	libgfortran.so.1()(64bit) is needed by rstudio-server-0.94.92-1.x86_64
	libssl.so.6()(64bit) is needed by rstudio-server-0.94.92-1.x86_64
[leipzig@localhost ~]$ sudo yum install libssl.so.6
Setting up Install Process
Package openssl098e-0.9.8e-17.el6.i686 already installed and latest version
Nothing to do

[leipzig@localhost ~]$ sudo rpm -Uvh --nodeps rstudio-server-0.94.92-x86_64.rpm
Preparing...                ########################################### [100%]
   1:rstudio-server         ########################################### [100%]
rsession: no process killed
Starting rstudio-server: /usr/lib/rstudio-server/bin/rserver: error while
loading shared libraries: libssl.so.6: cannot open shared object file: No such
file or directory
[FAILED]

#trying some stuff recommended here
#http://support.rstudio.org/help/discussions/problems/839-installing-rstudio-
from-source-after-installing-r-from-source

[leipzig@localhost ~]$ sudo yum install openssl098e-0.9.8e
Setting up Install Process
Resolving Dependencies
--&gt; Running transaction check
---&gt; Package openssl098e.x86_64 0:0.9.8e-17.el6 set to be updated
--&gt; Finished Dependency Resolution

Dependencies Resolved

================================================================================
================================================================================
==
 Package                                  Arch                               
Version                                       Repository                      
Size
================================================================================
================================================================================
==
Installing:
 openssl098e                              x86_64                             
0.9.8e-17.el6                                 sl                             
762 k

Transaction Summary
================================================================================
================================================================================
==
Install       1 Package(s)
Upgrade       0 Package(s)

Total download size: 762 k
Installed size: 2.2 M
Is this ok [y/N]: y
Downloading Packages:
openssl098e-0.9.8e-17.el6.x86_64.rpm                                            
                                                          | 762 kB     00:00    

Running rpm_check_debug
Running Transaction Test
Transaction Test Succeeded
Running Transaction
Warning: RPMDB altered outside of yum.
rstudio-server-0.94.92-1.x86_64 has missing requires of libcrypto.so.6()(64bit)
rstudio-server-0.94.92-1.x86_64 has missing requires of
libgfortran.so.1()(64bit)
rstudio-server-0.94.92-1.x86_64 has missing requires of libssl.so.6()(64bit)
  Installing     : openssl098e-0.9.8e-17.el6.x86_64                             
                                                                             1/1


Installed:
  openssl098e.x86_64 0:0.9.8e-17.el6                                            
                                                                                


Complete!
[leipzig@localhost ~]$ sudo yum install gcc41-libgfortran-4.1.2
Setting up Install Process
No package gcc41-libgfortran-4.1.2 available.
Error: Nothing to do
[leipzig@localhost ~]$ sudo yum install pango-1.28.1
Setting up Install Process
Package pango-1.28.1-3.el6_0.5.x86_64 already installed and latest version
Nothing to do
[leipzig@localhost ~]$ sudo rpm -Uvh --nodeps rstudio-server-0.94.92-x86_64.rpm
Preparing...                ########################################### [100%]
	package rstudio-server-0.94.92-1.x86_64 is already installed

[leipzig@localhost ~]$ sudo rstudio-server start
[leipzig@localhost ~]$ sudo rstudio-server verify-installation
Stopping rstudio-server:                                   [  OK  ]
/usr/lib/rstudio-server/bin/rsession: error while loading shared libraries:
libgfortran.so.1: wrong ELF class: ELFCLASS32
Starting rstudio-server:                                   [  OK  ]
[leipzig@localhost ~]$ sudo yum install libgfortran.so.1
Setting up Install Process
Package compat-libgfortran-41-4.1.2-39.el6.i686 already installed and latest
version
Nothing to do
[leipzig@localhost ~]$ sudo rpm -Uvh
ftp.scientificlinux.org/linux/scientific/6.0/x86_64/os/Packages/compat-
libgfortran-41-4.1.2-39.el6.x86_64.rpm
error: open of
ftp.scientificlinux.org/linux/scientific/6.0/x86_64/os/Packages/compat-
libgfortran-41-4.1.2-39.el6.x86_64.rpm failed: No such file or directory
[leipzig@localhost ~]$ wget
ftp.scientificlinux.org/linux/scientific/6.0/x86_64/os/Packages/compat-
libgfortran-41-4.1.2-39.el6.x86_64.rpm
--2011-08-18 04:39:39-- 
http://ftp.scientificlinux.org/linux/scientific/6.0/x86_64/os/Packages/compat-
libgfortran-41-4.1.2-39.el6.x86_64.rpm
Resolving ftp.scientificlinux.org... 131.225.110.147
Connecting to ftp.scientificlinux.org|131.225.110.147|:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: 128080 (125K) [application/x-rpm]
Saving to: “compat-libgfortran-41-4.1.2-39.el6.x86_64.rpm”

100%[===========================================================================
=============================================&gt;] 128,080      488K/s   in 0.3s   


2011-08-18 04:39:39 (488 KB/s) - “compat-libgfortran-41-4.1.2-39.el6.x86_64.rpm”
saved [128080/128080]

[leipzig@localhost ~]$ sudo rpm -Uvh
compat-libgfortran-41-4.1.2-39.el6.x86_64.rpm 
Preparing...                ########################################### [100%]
   1:compat-libgfortran-41  ########################################### [100%]
[leipzig@localhost ~]$ sudo rstudio-server verify-installation
Stopping rstudio-server:                                   [  OK  ]
Starting rstudio-server:                                   [  OK  ]

Big-Ass Servers™ and the myths of clusters in bioinformatics

2011-06-23T17:14:00.000-04:00

Spending $55k for a 512GB machine (Big-Ass Server™ or BAS™) can be a tough sell for a bioinformatics researcher to pitch to a department head.

Dell PowerEdge r900, available in orange and lemon-lime

Speaking as someone who keeps his copy of CLR safely stored in the basement, ready to help rebuild society after a nuclear holocaust, I am painfully aware of the importance of algorithm development in the history of computing, and the possibilities for parallel computing to make problems tractable.

Having recently spent 3 years in industry, however, I am now more inclined to just throw money at problems. In the case of hardware, I think this approach is more effective than clever programming for many of the current problems posed by NGS.

From an economic and productivity perspective, I believe most bioinformatics shops doing basic research would benefit more from having access to a BAS™ than a cluster. Here's why:

The development of multicore/multiprocessor machines and memory capacity has outpaced the speed of networks. NGS analyses tends to be more memory-bound and IO-bound rather than CPU-bound, so relying on a cluster of smaller machines can quickly overwhelm a network.
NGS has forced the number of high-performance applications from BLAST and protein structure prediction to doing dozens of different little analyses, with tools that change on a monthly basis, or are homegrown to deal with special circumstances. There isn't time or ability to write each of these for parallel architectures.

If those don't sound very convincing, here is my layman's guide to dealing with the myths you might encounter concerning NGS and clusters:

Myth: Google uses server farms. We should too.

Google has to focus on doing one thing very well: search.

Bioinformatics programmers have to explore a number of different questions for any given experiment. There is not time to develop a parallel solution to many of these questions as they will lead to dead ends.

Many bioinformatic problems, de-novo assembly being a prime example, are notoriously difficult to divide among several machines without being overwhelmed with messaging. You can imagine trying to divide a jigsaw puzzle among friends sitting several tables, you would spend more time talking about the pieces than fitting them together.

Myth: Our development setup should mimic our production setup

An experimental computing structure with a BAS™ allows for researchers to freely explore big data without having to think about how to divide it efficiently. If an experiment is successful and there is the need to scale-up to a clinical or industrial platform, that can happen later.

Myth: Clusters have been around a long time so there is a lot of shell-based infrastructure to distribute workflows

There are tools for queueing jobs, but those are often quite helpless to assist in managing workflows that are written as parallel and serial steps - for example, waiting for steps to finish before merging results.

Various programming languages have features to take advantage of clusters. For example, R has SNOW. But Rsamtools requires you to load BAM files into memory, so a BAS™ is not just preferable for NGS analysis with R, it's required.

Myth: The rise of cloud computing and Hadoop means that homegrown clusters are irrelevant that but also means we don't need a BAS™

The popularity of cloud computing in bioinformatics is also driven by the newfound ability to rent time on a BAS™. The main problem with cloud computing is the bottleneck posed by transferring GBs data to the cloud.

Myth: Crossbow and Myrna are based on Hadoop, we can develop similar tools

Ben Langmead, Cole Trapnell, and Michael Schatz, alums of Steven Salzberg's group at UMD, have developed NGS solutions using the Hadoop MapReduce framework.

Crossbow is a Hadoop-based implementation of Bowtie.
Myrna is an RNA-Seq pipeline.
Contrail is a de novo short read assembler.

These are difficult programs to develop, and these examples are also somewhat limited experimental proofs of concept or are married to components that may be undesirable for certain analyses. The Bowtie stack (Bowtie, Tophat, Cufflinks), while revolutionary in its implementation of Burroughs-Wheeler algorithm, is itself is built around the limitations of computers in the year 2008. For many it lacks the sensitivity to deal with, for example, 1000 Genomes data.

The dynamic scripting languages used most bioinformatics programmers are not as well suited to Hadoop as Java. To imply we can all develop similar tools of this sophistication is unrealistic. Many bioinformatics programs are not even threaded, much less designed to work amongst several machines.

Myth: embarrassingly parallel problems imply a cluster is needed

A server with 4 quad-core processors is often adequate for handling these embarrassing problems. Dividing the work just tends to lead to further embarrassments.

Here is a particularly telling quote from Biohaskell developer Ketil Malde on Biostar:

In general, I think HPC are doing the wrong thing for bioinformatics. It's okay to spend six weeks to rewrite your meteorology program to take advantage of the latest supercomputer (all of which tend to be just a huge stack of small PCs these days) if the program is going to run continously for the next three years. It is not okay to spend six weeks on a script that's going to run for a couple of days.

In short, I keep asking for a big PC with a bunch of the latest Intel or AMD core, and as much RAM as we can afford.

Myth: We don't have money for a BAS™ because we need a new cluster to handle things like BLAST

IBM System x3850 X5 expandable to 1536GB, mouse not included

Even the BLAST setup we think of as being the essence of parallel (a segmented genome index - every node gets a part of the genome) is often not the one that many institutions have settled on. Many rely on farming out queries to a cluster in which every node has the full genome index in memory.

Secondly, the mpiBLAST appears to be more suited to dividing an index among older machines than today's, which typically have >32GB RAM. Here is a telling FAQ entry:

I benchmarked mpiBLAST but I don't see super-linear speedup! Why?!

mpiBLAST only yields super-linear speedup when the database being searched is significantly larger than the core memory on an individual node. The super-linear speedup results published in the ClusterWorld 2003 paper describing mpiBLAST are measurements of mpiBLAST v0.9 searching a 1.2GB (compressed) database on a cluster where each node has 640MB of RAM. A single node search results in heavy disk I/O and a long search time.
http://www.mpiblast.org/Docs/FAQ#super-linear

Your comments on this topic are welcome!

RStudio: My thoughts

2011-03-15T13:17:00.001-04:00

Let me get this out of the way: I just love RStudio.

Created by a team lead by JJ Allaire, a name that should ring a bell if you were involved in web development during the Clinton administration, RStudio is an R IDE that is actually designed for R from the ground up. RStudio works on Linux, Mac, and Windows platforms, and can even run over the web.

While borrowing many of the best features from ESS, the Mac R-GUI, and maybe Anup Parikh's Red-R, RStudio provides solutions to several long-standing barriers that have hampered R code development. For instance, to do Sweave->tex->pdf (then view the pdf) in ESS was a frustrating, arthritic (M-n s M-n P)process that flummoxed even the greatest minds of our generation. RStudio has a handy button (Compile PDF) that brings you all the way from .Rnw to Acrobat. Although this command appears to run in its own session, leading to some unexpected behavior compared to running Sweave from the command line, the fact that this IDE is already geared for Sweave bodes well for future development.

This is fucking genius

The movement of commands back and forth from console to editor is another task that other editors made unnecessarily difficult - the old Mac R GUI console would not let you copy-and-paste a subset of the history, ESS was always geared to having users write code in the editor then executing lines but never writing code in the console then committing to the script. RStudio provides means of easily going in either direction. Control over multiple plots (solving both the overwritten X-Window and the annoying type=Cairo PNG problem on OS X) is a welcome relief.

RStudio offers very good autocompletion for such a relatively weird language - in addition to package methods it is aware of data frame columns and user-defined functions, for instance.

RStudio has already garnered a good number of suggestions. Here's personally what I would like to see:

More support for LaTeX markup, including menu driven formatting options so users don't have to memorize stuff like \textbf{}
More built-in aesthetic support for ggplot2, something where users are given a WYSIWYG manipulating an existing plot similar to Jeroen Ooms' ggplot2 web application
A non-sudo Linux binary and a method of specifying different R and TeX installations kicking around a server without re-installing from source.
Better control over the working directory (already reported and a likely future feature)
A means of quickly seeing where source files are actually located without mouseover
Integration with version control.
Code cleanup and indenting

Chromosome bias in R, my notebook

2010-12-23T15:42:00.000-05:00

My goal is to develop a means of detecting chromosome bias from a human BAM file.

Because I've been working with proprietary and novel plant genomes for the last three years, I haven't had the chance to use any of the awesome UCSC-based annotational features that have been introduced and refined in Bioconductor until now. I've returned to biomedical research and I have some catching up to do.

BSgenome might sound like horsecrap, but each Biostrings-based genome data package is actually a huge digested version of a UCSC/NCBI genome freeze and basic sequence annotation compiled into R objects.

BSgenome at Bioconductor

Be careful with googling bioconductor help - often the results point to older versions. Make sure your link has "release" in the url.

Here are the BSgenomes available today:

> available.genomes(type=getOption("pkgType"))
BioC_mirror = http://www.bioconductor.org
Change using chooseBioCmirror().
 [1] "BSgenome.Amellifera.BeeBase.assembly4" "BSgenome.Amellifera.UCSC.apiMel2"     
 [3] "BSgenome.Athaliana.TAIR.01222004"      "BSgenome.Athaliana.TAIR.04232008"     
 [5] "BSgenome.Btaurus.UCSC.bosTau3"         "BSgenome.Btaurus.UCSC.bosTau4"        
 [7] "BSgenome.Celegans.UCSC.ce2"            "BSgenome.Celegans.UCSC.ce6"           
 [9] "BSgenome.Cfamiliaris.UCSC.canFam2"     "BSgenome.Dmelanogaster.UCSC.dm2"      
[11] "BSgenome.Dmelanogaster.UCSC.dm3"       "BSgenome.Drerio.UCSC.danRer5"         
[13] "BSgenome.Drerio.UCSC.danRer6"          "BSgenome.Ecoli.NCBI.20080805"         
[15] "BSgenome.Ggallus.UCSC.galGal3"         "BSgenome.Hsapiens.UCSC.hg17"          
[17] "BSgenome.Hsapiens.UCSC.hg18"           "BSgenome.Hsapiens.UCSC.hg19"          
[19] "BSgenome.Mmusculus.UCSC.mm8"           "BSgenome.Mmusculus.UCSC.mm9"          
[21] "BSgenome.Ptroglodytes.UCSC.panTro2"    "BSgenome.Rnorvegicus.UCSC.rn4"        
[23] "BSgenome.Scerevisiae.UCSC.sacCer1"     "BSgenome.Scerevisiae.UCSC.sacCer2"

Select and load hg19

biocLite("BSgenome.Hsapiens.UCSC.hg19")
library("BSgenome.Hsapiens.UCSC.hg19")

When we get an alignment file one of the first things we want to do is look for red flags that might indicate something went awry in the lab or downstream. An example is chromosome bias - are we seeing more reads aligned to certain chromosomes than would be expected on size alone? A sticky question, since any experiment will introduce confounds based on the inherent uneven distribution of interesting genomic features, not to mention mapability. And yet I think this is still a worthwhile exercise and should be part of any ngs sequencing pipeline.

What we don't want to do is ignore that 7.6% of the GRCh37 freeze is sequence that looks like "NNNNNNN" - gaps representing unsequencable regions such as centromeres, scaffold gap delinations, and the like. We especially don't want to ignore gaps because they are not evenly distributed across the chromosomes (chrY is 56% gaps).

Raw chromosome length can be obtained from the BAM file header, but for this chromosome bias analysis I need the "non-gappy" length, the portion eligible for alignment. This is one of the "masks" turned on by default for BSgenomes in order to allow various functions to work properly (see MaskCollection in the IRanges package for more information).

> masks(Hsapiens)
Error in function (classes, fdef, mtable)  : 
  unable to find an inherited method for function masks, for signature "BSgenome"
#oops I see masks are a member of MaskedDNAString objects (i.e. chromosomes) not BSgenome objects
> masks(Hsapiens$chrY)
MaskCollection of length 4 and width 59373566
masks:
  maskedwidth maskedratio active names                               desc
1    33720000 0.567929506   TRUE AGAPS                      assembly gaps
2           0 0.000000000   TRUE   AMB   intra-contig ambiguities (empty)
3    16024357 0.269890426  FALSE    RM                       RepeatMasker
4      587815 0.009900281  FALSE   TRF Tandem Repeats Finder [period<=12]
all masks together:
  maskedwidth maskedratio
     49783032   0.8384713
all active masks together:
  maskedwidth maskedratio
     33720000   0.5679295
#I think the maskedwidth should reveal sum of actively masked nucleotides
> maskedwidth(Hsapiens$chrY)
[1] 33720000
#can we mess with the masks?
> active(masks(Hsapiens$chrY))["RM"]<-TRUE
Error in `$<-`(`*tmp*`, "chrY", value = < S4 object of class "MaskedDNAString">) : 
  no method for assigning subsets of this S4 class
#oops I can't manipulate a BSgenome this way - it is behaving like a class instead of an instance of a class
> chrY<-Hsapiens$chrY
> active(masks(chrY))["RM"]<-TRUE
> maskedwidth(chrY)
[1] 49744357
# ok maskedwidth is working as I figured, but i need unmasked width
> unmaskedWidth<-function(chr){length(chr)-maskedwidth(chr)}
> unmaskedWidth(Hsapiens$chrY)
[1] 25653566
#how can I iterate over something with a $ operator? let's try [[]]
> unmaskedWidth(Hsapiens[["chrY"]])
[1] 25653566

Now I want to create a data frame of with sequence names and unmaskedWidths to go with some read counts from a BAM file. Whenever I want to go from a list, through a function, to a data frame I think plyr, specifically ldply (list to data frame).

# let's take chr 1-22,X,Y, skipping the unscaffolded sequences and mitochondrial chr
> maskedSizes<-ldply(.data=seqnames(Hsapiens)[1:24],
  .fun=function(x){
    data.frame(chr=x,seqlength=length(Hsapiens[[x]]),
    unmaskedWidth=unmaskedWidth(Hsapiens[[x]]))},
  .progress="text",
  .parallel=TRUE)
> maskedSizes
                     chr seqlength unmaskedWidth
1                   chr1 249250621     225280621
2                   chr2 243199373     238204518
3                   chr3 198022430     194797135
4                   chr4 191154276     187661676
5                   chr5 180915260     177695260
6                   chr6 171115067     167395066
7                   chr7 159138663     155353663
8                   chr8 146364022     142888922
9                   chr9 141213431     120143431
10                 chr10 135534747     131314738
11                 chr11 135006516     131129516
12                 chr12 133851895     130481393
13                 chr13 115169878      95589878
14                 chr14 107349540      88289540
15                 chr15 102531392      81694766
16                 chr16  90354753      78884753
17                 chr17  81195210      77795210
18                 chr18  78077248      74657229
19                 chr19  59128983      55808983
20                 chr20  63025520      59505520
21                 chr21  48129895      35106642
22                 chr22  51304566      34894545
23                  chrX 155270560     151100560
24                  chrY  59373566      25653566

Load the BAM file and get read counts in a data frame.

#other methods include scanBam and readAligned
bamFile<-readBamGappedAlignments("myIndexedSortedBamFile.bam")
> levels(rname(bamFile))
 [1] "1"          "2"          "3"          "4"          "5"         
 [6] "6"          "7"          "8"          "9"          "10"        
[11] "11"         "12"         "13"         "14"         "15"        
[16] "16"         "17"         "18"         "19"         "20"        
[21] "21"         "22"         "X"          "Y"          "MT"        
[26] "GL000207.1" "GL000226.1" "GL000229.1" "GL000231.1" "GL000210.1"
[31] "GL000239.1" "GL000235.1" "GL000201.1" "GL000247.1" "GL000245.1"
[36] "GL000197.1" "GL000203.1" "GL000246.1" "GL000249.1" "GL000196.1"
[41] "GL000248.1" "GL000244.1" "GL000238.1" "GL000202.1" "GL000234.1"
[46] "GL000232.1" "GL000206.1" "GL000240.1" "GL000236.1" "GL000241.1"
[51] "GL000243.1" "GL000242.1" "GL000230.1" "GL000237.1" "GL000233.1"
[56] "GL000204.1" "GL000198.1" "GL000208.1" "GL000191.1" "GL000227.1"
[61] "GL000228.1" "GL000214.1" "GL000221.1" "GL000209.1" "GL000218.1"
[66] "GL000220.1" "GL000213.1" "GL000211.1" "GL000199.1" "GL000217.1"
[71] "GL000216.1" "GL000215.1" "GL000205.1" "GL000219.1" "GL000224.1"
[76] "GL000223.1" "GL000195.1" "GL000212.1" "GL000222.1" "GL000200.1"
[81] "GL000193.1" "GL000194.1" "GL000225.1" "GL000192.1"
#the deflines in my reference do not match the BSgenome names, must fix at least the chromosomes of interest
levels(rname(bamFile))[1:25]<-paste('chr',c(1:22,'X','Y','M'),sep='')

#run length encoded read counts per chromosome
readRle<-rname(bamFile)

#get a data frame with chromosome and read counts
allReadsDf<-ldply(runValue(readRle),function(x){data.frame(chr=levels(runValue(readRle))[x],reads=runLength(readRle)[x])})
> head(allReadsDf)
   chr   reads
1 chr1 3616909
2 chr2 3642052
3 chr3 2843019
4 chr4 2636141
5 chr5 2590352
6 chr6 2497123

Merge the read counts with unmasked chromosome lengths and plot their relationship.

chrSizesReads<-merge(maskedSizes,readCounts,sort=FALSE)
library(ggplot2)
p<-ggplot(data=chrSizesReads, aes(x=unmaskedWidth, y=reads, label=chr)) + 
  geom_point() +
  geom_text(vjust=2,size=3) +
  stat_smooth(method="lm", se=TRUE,level=0.95) +
  ylab("Reads aligned") +
  xlab("Unmasked chromosome size") +
  opts(title = "Reads vs Chromosome Size")
print(p)

There should be a strong linear relationship between read count and chromosome size. We can test this using a linear regression model, the null hypothesis being the number of reads aligned to a chromosome is independent of its size.

> mylm<-lm(reads~unmaskedWidth,data=chrSizesReads)
> mysummary<-summary(mylm)
> mysummary

Call:
lm(formula = reads ~ unmaskedWidth, data = chrSizesReads)

Residuals:
    Min      1Q  Median      3Q     Max 
-271816 -108122  -43984   42826  676284 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.774e+05  9.505e+04   1.866   0.0754 .  
unmaskedWidth 1.455e-02  7.145e-04  20.365 9.12e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 206600 on 22 degrees of freedom
Multiple R-squared: 0.9496, Adjusted R-squared: 0.9473 
F-statistic: 414.8 on 1 and 22 DF,  p-value: 9.123e-16

The low p-value (that chr size has no influence) and R-squared (predictive value of the linear model) suggest this model is sound.

The following plot is obtained from the standardized residuals (the standardized difference between data observed and values expected) of the linear model described earlier.

Chromosome bias refers to uneven read alignment distribution across various chromosomes. We can expect some chromosome bias in treatment sets because of the inherient nature any experimental conditions - recovered fragments will not be evenly distributed among chromosomes because regions of affect are not evenly distributed. Other possible factors of chromosome bias include heterochromatin, uneven repeat content, and the potential for aligning the against an incorrect set of sex chromsomes. Aligners will typically randomly, evenly, assign discrete positions to reads which map ambiguously to multiple locations.

> p<-qplot(chrSizesReads$chr,rstandard(mylm))+
   aes(label=chrSizesReads$chr)+
   geom_text(vjust=2,size=3)+
   xlab("Chromosome")+
   ylab("Std Residual from lm (reads)")+
   geom_abline(slope=0,intercept=0)+
   opts(axis.text.x = theme_text(angle=45,hjust=1))+
   opts(title = "Linear Regression Residuals")
> print(p)

Fortunately, there is no clear pattern to these residual values, which would indicate some model problems, but with a Z-score of 3.36, chrX appears to be an outlier. With 46M total alignments this is certainly not due to sampling error, but we can still test our observation with a Lund statistic.

#http://stackoverflow.com/questions/1444306/how-to-use-outlier-tests-in-r-code
>lundcrit<-function(a, n, q) {
 F<-qf(c(1-(a/n)),df1=1,df2=n-q-1,lower.tail=TRUE)
 crit<-((n-q)*F/(n-q-1+F))^0.5
 crit
}
> n<-nrow(chrSizesReads)
> q<-length(mylm$coefficients)
> crit<-lundcrit(0.05,n,q)
> chrSizesReads[which(rstandard(mylm)>crit),"chr"]
[1] chrX

Happy holidays!

Directory-based bash histories

2010-12-09T17:41:00.000-05:00

Using a directory-based bash history allows for a record of shell actions on a directory basis, so a group of developers have some record of what was done while in a directory, when, and by whom. This can be helpful when trying to reconstruct history with limited documentation.

I know this setup will be of some benefit to my successor at my previous job because he has access to everything I ever did in any project directory.

Place this code in your ~/.bash_profile or ~/.bashrc

(type source ~/.bash_profile (or .bashrc) to load this for your current session)

NGS viewers reviewed

2010-08-30T14:12:00.000-04:00

I gathered up some of the recent free next generation sequence viewers that were capabale of viewing BAM files - and put each through the motions with a few BAM files and reference sequences of various sizes. While there are some great ideas and several choices to be found along the feature spectrum, I think we are still in the dark ages with this stuff. No viewer has really been able to entirely combine usability with performance and analysis capabilities, let alone extensibility and web connectivity.

tview
My take: tview is the barebones, text-rendered viewer that is included with Samtools. People who favor this as their BAM viewer probably think Vim is too polished. Even the very limited navigation is remarkably unituitive (goto a coordinate requires chr:position even if you have just one chromosome, no errors are displayed if you forget chr).
Most resembles which video game: Oregon Trail
Good standout feature: command line access
Bad feature: text display

BamView
My take: Bamview is a wicked fast simple BAM file viewer. It doesn't have much in the way of features, but for cursory examination of BAM files it is more palatable than tview.
Most resembles which video game: Burgertime
Good standout feature: strand split screen
Bad feature: drag selecting a region turns it red, and umm... that's all it does

GenoViewer
My take: The only NGS viewer endorsed by Speak the Hungarian rapper, unfortunately this recent entrant leaves a lot to be desired in terms of performance with large files. GenoViewer is very hard on the eyes - the indiscriminate use of primary colors looks like a kid somehow vomited up the ball pit at Chuck E. Cheese.
Most resembles which video game: Centipede
Good standout feature: promises that "You will not get lost in the details, and can easily figure out the true meaning behind the data. Guaranteed."
Bad feature: graphics

MagicViewer
My take: MagicViewer has come a long way since its initial release last December. With its interesting pie-chart icon renderings of SNP purity, decent treatment of annotation tracks, and improved performance, MagicViewer might soon be a contender to Tablet in the midweight category. The navigation is workable but takes some getting used to - it's unclear when to use scroll bars vs. arrows
Most resembles which video game: Moon Patrol
Good standout feature: primer design tool
Bad feature: some regions are simply not visible - no reference, no reads

Tablet
My take: Hands down the most attractive of the viewers, Tablet is aimed at fostering a delicate balance between performance, features, and aesthetics. Tablet comes with a suite of read views - Packed, Stacked, and Classic - to suit both young children and elderly scientists alike. GFF feature files can be loaded but they appear to merely serve as position indices.
Most resembles which video game: SimCity
Good standout feature: interface
Bad feature: read insertions not displayed correctly

IGV
My take: The Integrative Genomics Viewer is a serious tool for exploring and analyzing large datasets. In addition to viewing, IGV is designed to allow users to extract the kind of hard publishable data that has typically been the domain of Bio* scripts. Like the true product of an Ivy League education, IGV can appear aloof and arrogant to newcomers. The viewer will let you load BAM files and other annotation tracks that have nothing to do with the reference without comment or guidance, then require an extra unitutive click to actually generate a view. While not the easiest or even the best in performance, in terms of generating real queries on your data from the viewer there is nothing comparable.
Most resembles which video game: Dig Dug
Good standout feature: analysis tools
Bad feature: throws a hissy fit if it cannot connect to home server

gbrowse2
My take: gbrowse2 is the AJAXified protege to the venerable generic genome browser. Samtools integration is a recent addition to this highly extensible platform, which has been used for years to display everything from large genomes to small sequencing projects. Of all the viewers here, GB2 provides the best set of visualization tools, such that virtually any biological information that can be rendered linearly has been done so as gbrowse tracks. The lone true web application, gbrowse genomic positions can be hyperlinked or even snapshot-embedded on the web. The web provides the best platform to share visual genomic data among several users. However, a gbrowse2 installation with BAM tracks can be a massive pain to install, configure, and debug ("landmark chrI not found" is the most popular google search in all of bioinformatics). Novices can expect a minefield of historical gotchas and arcane conventions ("Name=" not "name=" field in gff3, bp_seqfeature_load instead of bp_load_gff for gff3, no validator for conf files), and even experienced users are often baffled by cryptic errors that pop up in server logs.
Most resembles which video game: Sissyfight 2000
Good standout feature: hyperlinks
Bad feature: setup

My thoughts on the acquisition of Ion Torrent by Life Technologies

2010-08-18T12:42:00.000-04:00

Yesterday it was announced that Ion Torrent, makers of the Ion Personal Genome Machine (PGM) sequencer, would be acquired by ~~Invitrogen~~ ~~ABI~~ Life Technologies for $375M in cash and stock, with the possibility of another $350M if various milestones are met.

Gregory Lucier, chairman and CEO of Life Technologies, said, "We believe Ion Torrent's technology will represent a profound change for the life sciences industry, as fundamental as the one we saw with the introduction of qPCR." That analogy might not sound earth-shattering, but I suspect he is using qPCR as an example because it is one of the few molecular biology (as opposed to biochemical) techniques regularly used in clinical settings.

The PGM is a unique machine because it is the first second generation sequencer (tentative specs: ~3 million 200bp reads at $500 1hr run, $50k machine) that could be conceivably leased by a small clinical testing facility, like the ones fed by those LabCorp boxes you see scattered all over strip malls. Bear in mind even a capillary-based sanger sequencer like the ABI3730xl costs a whopping $375,000, which comes as a shock to those of us who mostly work with next-generation sequencing data. You can see not only why the PGM is really a game changer, but how it might fit into Life Technologies offerings.

Exactly which clinical diagnostics are, or will be, suited for this machine are unclear. With the right foolproof software this could replace a lot of PCR and microarray-based tests. Read lengths and throughput are bound to increase just as they did with Solexa. Future applications like tumor sequencing and 16S rRNA microbiome sequencing have not even entered into medical practice yet.

The key to Life Technology succeeding in these areas will likely come down to non-technical challenges:

Getting FDA approval for the PGM as a medical device and for various protocols based on the PGM as approved clinical tests. Life Technology has experience with this process. Ion Torrent clearly does not.
Convincing insurance companies and HMOs that these tests are cost-effective diagnostics. When you consider the exorbitant cost of other tests - e.g. several MRIs over a course of chemotherapy - this might not be such hard sell.
Developing a business model that will allow small clinics to lease the machine for little or even no cost provided they agree to purchase a minimum amount of consumables. Life Technologies will likely market the PGM to smaller labs who themselves will inevitably face competition from larger dedicated sequencing centers trying to centralize this type of work using bigger machines from Illumina or PacBio (or even LT's SOLiD).

Although these obstacles seem daunting, I have sat through academic departmental meetings where very knowledgeable sequencer salespeople were invited back multiple times only to be jerked around as the faculty hemmed and hawed over whether this was the right $400k machine to buy at the time. That was undoubtedly an expensive and frustrating sales process for Solexa and 454. With the PGM, Life Technologies can pitch this much cheaper machine to an albeit different group of skeptics, clinical lab managers, but one with a clear bottom line in mind.

Use local instead of my in Perl when using $$ (double dollar sign) inside an enclosed block

2010-03-10T13:21:00.000-05:00

I thought I would be all clever and initialize several hash indices at once using the $$ notation in Perl - evaluating strings as variables. Unfortunately using "my" screws this up bigtime.

This must be another circumstance spelled out in the incomprehensible "my" treatise:
http://perldoc.perl.org/perlsub.html#Private-Variables-via-my()


my %foo;
$foo{'bar'}=1;
foreach $fooDex(qw(foo)){
    $$fooDex{'bar'}=2;
}
print $foo{'bar'}."\n"; #prints 1


local %foo;
$foo{'bar'}=1;
foreach $fooDex(qw(foo)){
    $$fooDex{'bar'}=2;
}
print $foo{'bar'}."\n"; #prints 2

This is one of those solutions that gives you a feeling of dread instead of accomplishment.

Getting the basics from readAligned

2010-03-09T16:57:00.000-05:00

The UCR guide is a little sparse with regard to getting basic information from readAligned.

I'd like to add to the general cookbook. If some bioc people out there can contribute some alignment recipes can fill me in on some more basics please comment:


alignedReads <- readAligned("./", pattern="output.bowtie", type="Bowtie")

#how many reads did I attempt to align
#i don't think you can't get this from alignedReads

#how many reads aligned (one or more times)
length(unique(id(alignedReads)))

#how many hits were there?
length(alignedReads)

#how many reads produced multiple hits
length(unique(id(alignedReads[srduplicated(id(alignedReads))])))

#how many reads produced multiple hits at the best strata?
#please fill me in on this one

#how many reads aligned uniquely (with exactly one hit)
length(unique(id(alignedReads)))-length(unique(id(alignedReads[srduplicated(id(alignedReads))])))

#how many reads aligned uniquely at the best strata (the other hits were not as good)
#please fill me in on this one

#how many unique positions were hit? what if I ignore strand?
#please fill me in on this one

#how many converging hits were there (two query sequences aligned to the same genomic position)
#please fill me in on this one

Quality trimming in R using ShortRead and Biostrings

2010-03-03T16:48:00.000-05:00

I wrote an R function to do soft-trimming, right clipping FastQ reads based on quality.

This function has the option of leaving out sequences trimmed to extinction and will do left-side fixed trimming as well.

#softTrim
#trim first position lower than minQuality and all subsequent positions
#omit sequences that after trimming are shorter than minLength
#left trim to firstBase, (1 implies no left trim)
#input: ShortReadQ reads
#       integer minQuality
#       integer firstBase
#       integer minLength
#output: ShortReadQ trimmed reads
library("ShortRead")
softTrim<-function(reads,minQuality,firstBase=1,minLength=5){
qualMat<-as(FastqQuality(quality(quality(reads))),'matrix')
qualList<-split(qualMat,row(qualMat))
ends<-as.integer(lapply(qualList,
function(x){which(x < minQuality)[1]-1}))
#length=end-start+1, so set start to no more than length+1 to avoid negative-length
starts<-as.integer(lapply(ends,function(x){min(x+1,firstBase)}))
#use whatever QualityScore subclass is sent
newQ<-ShortReadQ(sread=subseq(sread(reads),start=starts,end=ends),
quality=new(Class=class(quality(reads)),
quality=subseq(quality(quality(reads)),
start=starts,end=ends)),
id=id(reads))

#apply minLength using srFilter
lengthCutoff <- srFilter(function(x) {
width(x)>=minLength
},name="length cutoff")
newQ[lengthCutoff(newQ)]
}

To use:

library("ShortRead")
source("softTrimFunction.R") #or whatever you want to name this
reads<-readFastq("myreads.fq") trimmedReads<-softTrim(reads=reads,minQuality=5,firstBase=4,minLength=3) writeFastq(trimmedReads,file="trimmed.fq")

I strongly recommend reading the excellent UC Riverside HT-Sequencing Wiki cookbook and tutorial if you wish to venture into using R for NGS handling. Among other things, it will explain how to perform casting if you have Solexa scaled (base 64) fastq files. The function should respect that. http://manuals.bioinformatics.ucr.edu/home/ht-seq

Using Vmatch to combine assemblies

2009-11-19T13:56:00.000-05:00

In praise of Vmatch

If I could take only one bioinformatics tool with me to a desert island it would be Vmatch.

In addition to being a versatile alignment tool, Vmatch has many lesser-known features which also leverage suffix trees. The dbcluster option allows Vmatch to cluster similar sequences using Vmatch index files. This method is about 1000x faster than attempting to align all sequences against each other, no matter how clever the algorithm.

Originally presented as means to cluster EST sequences, this option is useful for processing the output from de novo short read assemblers like Velvet and ABYSS. Vmatch -dbcluster allows us to easily create a non-redundant set of contig sequences originating from several assemblies.

Why would we want to combine assemblies?

Every kmer is sacred

Why not just take the one with the highest N50? At the recent Genome Informatics meeting at CSHL, Inanc Birol (ABYSS) said "every kmer is sacred". Our experience with the de novo transcriptome assembly of plants has been that the best set of contigs is spread out all over the parameter space. Longer kmer and cvCut settings can produce a longer set of elite contigs at the expense of omitting lowly expressed (or spurious?) contigs that appear at less stringent settings.

Reads held hostage

To the point above, I would add "every cvCut is also sacred". That doesn't exactly roll off the tongue, but we have seen some instances in which an assembly will have a higher read usage at a cvCut of 10 than at 5. This paradox suggests there is a "contig read hostage" situation, in which reads critical to the formation of longer contigs can be held captive in low coverage short contigs. Raising the cvCut threshold frees up these critical reads to extend more plausible contigs and allowing additional reads from the pool to be recruited into the assembly.

What is a non-redundant set?

The following sequences have some redundancies. NODE2 and NODE3 do not add any information.
The non-redundant set will have only NODE1 and NODE4.

>NODE1
AATTCAGTTGAAGTAATAGAGGCAGCTGCTGTTAGAACTTCGCTACGACTAGCAACGCTATGTCGAGTTGTACCTTCCCACCTCGTATACAAGGAGCATGAAGTCATCAGCC
CTTTCTACAGAATCTAGGTCCGAAATGATGAAATTAAGAGAAAGACAACTGAAGGTTTTAGGAGGCAAGGTTTACAGGGTAATGAACACGGAAAAACCCACAAGCTAGGAAC
AGTGTGTCTTGGAGTTTAAACTGATTTGGTAGTAGTTCGAAAACAAATTGGAAGGGACATTTAAAGTCCGAGTTGACGTTATCTGAGACAACTTTGTCTTTAACCGACAGGG
AGTTGAGGTAGGAGAGAGTGTCCACATATTTAACGTTGTTCAGATATGGGATGTAGCAGTTGTAACCGAAGCATGTAGGAAGGTTAAAGGGGTCCATACCTCTATTCTAGTC
CCGAAGGTTGGTAGCTAGACTCGGTGACCTAAAATGAGAACGGAAGAAACGGAGGTGACATCCATGGGGCTTGTCGTATATCCAATCATACCTTTGGAGAAGGAAGTTAAAG
GTCAAAACTTTAAAAACCATGAGGACCATTTTATCCTCGTCACTGTCGACTATGGTGAAGTGACCTAGCAGGTTTGTGTAGTTACTGTTTTGTAAGTGTAAAGTGCGTTGCT
GGCTATAGGAGTTTCGCATGAAACATGCTCGCTCTTGTGACCCATCGGTTACCAAGTTTACAGACTAGGTGAGGACTAGGTCACCTCTGTTTGGTAACCAAAAAGTAGAGAA
AGAATATCAACAAGGTATATACAACAACTTAGAGTAGCAAACCTTAAAAGAGTTGTCGTCAGTCACAAAGGGAGAGTTGTACTAAAGGGAACTTTTGTTCACAGAGCTAAAA
CTCATTAAGACCATTTTAATGTGGTCTACCAAGTCTGTCGAAAAGAAGTGGATGCTACGGCAGA
>NODE2 a substring of NODE1
AATTCAGTTGAAGTAATAGAGGCAGCTGCTGTTAGAACTTCGCTACGACTAGCAACGCTATGTCGAGTTGTACCTTCCCACCTCGTATACAAGGAGCATGAAGTCATCAGCC
CTTTCTACAGAATCTAGGTCCGAAATGATGAAATTAAGAGAAAGACAACTGAAGGTTTTAGGAGGCAAGGTTTACAGGGTAATGAACACGGAAAAACCCACAAGCTAGGAAC
AGTGTGTCTTGGAGTTTAAACTGATTTGGTAGTAGTTCGAAAACAAATTGGAAGGGACATTTAAAGTCCGAGTTGACGTTATCTGAGACAACTTTGTCTTTAACCGACAGGG
AGTTGAGGTAGGAGAGAGTGTCCACATATTTAACGTTGTTCAGATATGGGATGTAGCAGTTGTAACCGAAGCATGTAGGAAGGT
>NODE3 a reverse complement substring of NODE1
CGACAGACTTGGTAGACCACATTAAAATGGTCTTAATGAGTTTTAGCTCTGTGAACAAAAGTTCCCTTTAGTACAACTCTCCCTTTGTGACTGACGACAACTCTTTTAAGGT
TTGCTACTCTAAGTTGTTGTATATACCTTGTTGATATTCTTTCTCTACTTTTTGGTTACCAAACAGAGGTGACCTAGTCCTCACCTAGTCTGTAAACTTGGTAACCGATGGG
TCACAAGAGCGAGCATGTTTCATGCGAAACTCCTATAGCCAGCAACGCACTT
>NODE4 a new sequence
TTACGAACGATAGCATCGATCGAAAACGCTACGCGCATCCGCTAAGCACTAGCATAATGCATCGATCGATCGACTACGCCTACGATCGACTAGCTAGCATCGAGCATCGATC
AGCATGCATCGATCGATCGAT

How do we create a non-redundant set?

Do not use -nonredundant

No, really. There is an option called -nonredundant which should presumably do what we want, but unfortunately that writes a "representative" member of each cluster to a file, which may or may not be the longest contig. I'm not sure what makes a sequence representative of a cluster, but for this application we want the longest member of each cluster.

On April 29th 2010, Vmatch 2.1.3 was released. The most important change is that the option -nonredundant now delivers the longest sequence from the corresponding cluster (instead
of an unspecified representative). This should make the longestSeq.pl approach unnecessary.

To create a non-redundant set we will produce cluster files and then extract the longest sequence from each file. Use the following commands to produce your cluster files from an index:

mkvtree -allout -pl -db contigs1.fa contigs2.fa -dna -indexname myIndex
#mkvtree will accept multiple fasta or gzipped fasta files

#if using Vmatch 2.1.3 or later:
vmatch -d -p -l 25 -dbcluster 100 0 -v -nonredundant nonredundantset.fa myIndex > mySeqs.rpt
# thx to Hamid Ashrafi for debugging this syntax


#if using older Vmatch
mkdir sequenceMatches 
#this is where vmatch will put each cluster

vmatch -d -p -l 25 -dbcluster 100 0 mySeqs "(1,0)" -v -s myIndex > mySeqs.rpt
mv mySeqs*.match mySeqs*.fna sequenceMatches

What do these options do?

-d direct matches (forward strand)
-p palindromic matches (reverse strand)
-l search length (set this below your shortest sequence)
-dbcluster queryPerc targetPerc - runs an internal alignment of the index created by mkvtree.
- The two numeric arguments specify what percentage of your query sequence (the smaller) is involved in an alignment to the cluster sentinel/target superstring. For our purposes we require 100% of our query sequence substring to match the target. We don't care what percentage of the target is aligned, so set the second parameter to 0.
- The third argument to dbcluster is the index prefix name that vmatch will give to the THOUSANDS of cluster .fna and .match files it will create.
- The fourth argument "(1,0)" specifies that we want to keep singletons (in a file called mySeqs.single.fna) and that there is no limit to the number of sequences in an acceptable cluster.
-v verbose report (redirected to a the .rpt)
-s create the individual cluster fasta files and match reports

Consult the vmatch manual for fuzzy matches and more examples:
http://www.zbh.uni-hamburg.de/vmatch/virtman.pdf

The standard output details the clusters and singlets.

0:761437: NODE_83_length_31_cov_22.451612857642: 
NODE_106_length_31_cov_22.451612621409: NODE_152_length_29_cov_27.758621749981: 
NODE_185_length_29_cov_27.7586211:761861: NODE_531_length_424_cov_26.851416805590: 
NODE_1413_length_320_cov_28.187500837480: NODE_1236_length_320_cov_28.187500858407: 
NODE_937_length_320_cov_28.187500765510: NODE_1542_length_108_cov_34.870369621979: 
NODE_786_length_425_cov_32.602352750915: NODE_1290_length_321_cov_34.392525

Extract the longest sequence from your cluster files

If using Vmatch 2.1.3 or later this is unnecessary
Here is a perl script, which we will call longestSeq.pl, to do that

#!/usr/bin/perl
#print the longest sequence in a fasta file
use strict;my %seqs;my $defline;
while (<>) {
  chomp;
  if ( /^#/ || /^\n/ ) {
    #comment line or empty line do nothing
  }elsif (/>/) {
    s/>//g;$defline = $_;
    $seqs{$defline} = "";
  }elsif ($defline) {
    $seqs{$defline} .= $_;
  }else{
    die('invalid FASTA file');
  }
}
my $max = $defline;
foreach my $def ( keys %seqs ) {
  $max = ( length( $seqs{$def} ) > length( $seqs{$max} ) ) ? $def : $max;
}
print ">" . $max . "\n" . $seqs{$max} . "\n";

We want the longest member of each cluster and all sequences in the singletons file:

for f in sequenceMatches/*fna;
do
if [ "$f" = "sequenceMatches/mySeqs.single.fna" ];
then 
cat $f >> mySeqs_longest_seqs.fa;
else perl longestSeq.pl $f >> mySeqs_longest_seqs.fa;
fi;
done

That's it - now you have a comprehensive and non-redundant set of the longest contigs from a number of assemblies.

R's xtabs for total weighted read coverage

2009-11-04T12:54:00.000-05:00

Samtools and its BioPerl wrapper Bio::DB:Sam prefer to give read coverage on a depth per base pair basis. This is typically an array of depths, one for every position that has at least one read aligned. OK, works for me. But how can we quickly see which targets (in my case transcripts) have the greatest total weighted read coverage (i.e. sum every base pair of every read that aligned)?

My solution is to take that target-pos-depth information and import a table into R with at least the following columns:
targetName
depth

I added the pos column here to emphasize the base-pair granularity

         tx pos depth
1   tx500090 227     1
2   tx500090 228     1
3   tx500090 229     1
4   tx500090 230     1
5   tx500090 231     1
...
66  tx500123 184     1
67  tx500123 185     1
68  tx500123 186     1
69  tx500123 187     2
70  tx500123 188     2
71  tx500123 189     2

In R

myCoverage<-read.table("myCoverage.txt",header=TRUE)
myxTab<-xtabs(depth ~ tx,data=myCoverage)

xtabs will sum up depth-weighted positions by default (i suppose this is what tabulated contigency really means) and return an unsorted list of transcripts and their weighted coverage (total base pair read coverage)

> myxTab[1:100]
tx
tx500090 tx500123 tx500134 tx500155 tx500170 tx500178 tx500203 tx500207
      38       92      610       46      176       46       92      130
tx500273 tx500441 tx500481 tx500482 tx500501 tx500507 tx500667 tx500684
      76     2390      114    71228      762      222      542      442
tx500945 tx500955 tx501016 tx501120 tx501127 tx501169 tx501190 tx501192
    1378     3604       46       46      420      854      130      352
tx501206 tx501226 tx501229 tx501245 tx501270 tx501297 tx501390 tx501405
     244     1204      206    15926      214       46      168       46
tx501406 tx501438 tx501504 tx501694 tx501702 tx501877 tx501902 tx502238
      38     2572     7768     3274      314      298       84      198
tx502320 tx502364 tx502403 tx502414 tx502462 tx502515 tx502517 tx502519
     122       38      588       46       46       38       38      466
tx502610 tx502624 tx502680 tx502841 tx502882 tx503090 tx503192 tx503204
     206       38      168     3750       38      122       76       92
tx503416 tx503468 tx503523 tx503536 tx503571 tx503578 tx503623 tx503700
     260       38      168       38       46       46       84       38
tx503720 tx503721 tx503722 tx503788 tx503872 tx503892 tx503930 tx503970
   97112       38       38     4708       38       38     1290       84
tx503995 tx504107 tx504115 tx504346 tx504353 tx504355 tx504357 tx504398
      46      152      206       46     3416     1402      122      290
tx504434 tx504483 tx504523 tx504589 tx504612 tx504711 tx504751 tx504827
     290     8728      176       46       46       76     5644     1308
tx504828 tx504834 tx504882 tx504931 tx504952 tx505017 tx505029 tx505078
    2336      328       46    34138     1000     1838       46      474
tx505123 tx505146 tx505159 tx505184
      38   123344      160      588

this is approximately 10000x faster than using a formula like:

by(myCoverage,myCoverage$tx,function(x){sum(x$depth)}

Installing Bio::DB::Sam from CPAN

2009-10-23T15:22:00.000-04:00

Bio::DB::Sam is Lincoln Stein's BioPerl API to the SamTools package.

Installing via CPAN might skip a necessary question that will cause it to fail.


cpan> install Bio::DB::Sam
...
...
DIED. FAILED tests 1-93
  Failed 93/93 tests, 0.00% okay
Failed Test Stat Wstat Total Fail  Failed  List of Failed
-------------------------------------------------------------------------------
t/01sam.t      2   512    93  186 200.00%  1-93
Failed 1/1 test scripts, 0.00% okay. 93/93 subtests failed, 0.00% okay.
make: *** [test_dynamic] Error 2
/usr/bin/make test -- NOT OK
Running make install
make test had returned bad status, won't install without force

Navigate to where CPAN has downloaded the gz ffile
A closer examination reveals that Build.PL file wants to know where the SamTools header files are located.


Please enter the location of the bam.h and compiled libbam.a files:

I have no idea how to pass these arguments using CPAN. I would just avoid this method of installation. Do the local build instead.

Disable file locking in Eclipse for OS X

2009-10-12T16:46:00.000-04:00

Eclipse will refuse to use a workspace on an automounted OS X Server home directory.

Workspace in use or cannot be created

To remedy this problem do the following:

Right click the Eclipse application and select "Show Package Contents"
Contents->MacOS
Edit the eclipse.ini file in a text editor
Add -Dosgi.locking=none to the line below -vmargs

My 2009 Bridge-to-Bridge Experience

2009-09-23T20:18:00.000-04:00

Bridge-to-Bridge is a 105-mile ride up to the top of Grandfather Mountain, one of the highest peaks in the Blue Ridge mountains.

Although B2B is not an officially sanctioned race, the organizers conduct it just as professionally (with the exception of neutral support). Cops manage the major crossings, volunteers provide hand-offs at the dozen feed stations, and the event is officially timed with the aid of some magnetic shoe things.

I last did this ride in 1999. Now 10 years older but about 10 pounds lighter I had somehow forgotten how much suffering was involved and figured this was a good time to tackle the challenge, despite falling ill a couple of weeks beforehand.

Due to constant rain and very heavy fog, this year was utter torture for the 299 finishers and 371 non-finishers who braved the elements. I believe there may have been another 130 non-starters who stayed in their hotel rooms enjoying the Golden Girls marathon on tv. While I'm sure I would have felt some sense of accomplishment doing that, I was obligated to finish this ride as we had already driven down from Philadelphia (en-route to a wedding in Nashville the following weekend).

The Grandfather Mountain staff said riders could not enter the park before 3pm. To accommodate both the riders and this odd rule two start times were offered - 10 and 11 am, with slower riders encouraged to start first. To me this was a welcome change from the ungodly pre-dawn start times of most big rides. Riders were advised to start at the later time if they estimated they would be pushing the 3pm threshold.

I was on the fence about which group to join. I saw several very fit looking riders and expensive bikes in the 10am pack. I was still riding the same Litespeed I used in 1999. The word going around was that more rain was on the way (this turned out to be true for everyone). In the end I felt the risk of having an inexperienced rider fall in front of me to be the deciding factor to go with the second group. I knew I would have to pass several of that first group anyway but it would be on the later climbs instead of the early rollers. Ironically I almost got clipped by some idiot drifting carelessly in our pack. At the end there was considerable overlap in finishing times between the two groups.

The two times I did this ride ('98 and '99) I got dropped by the leaders on the 13 mile climb up NC181, then spent the rest of the ride either riding alone or with a couple other guys. Despite my best efforts this year turned out roughly the same except I stuck with that front group for about half the climb instead of just the first couple miles. Remind me to buy a compact crank.

This climb is very difficult psychologically - a relentless slog up a roughly-paved 4 lane highway. I was not very familiar with the profile and prematurely thought I had crested three times - each time putting in a kick over the "top". The fog would clear and I would see yet another rise.

Feeling dispirited and exhausted, I nearly froze to death on the descents of 181 and the Blue Ridge Parkway and was eventually caught by a small group that stuck together until the ascent of Grandfather. I was amazed by how few words were exchanged in that group during the hour or so we traded pulls through the fog and rain, which only got worse as we neared the finish. One guy did say something that stuck with me, "It's like we're riding through a horror movie."

After crawling to the finish in a 39x27, I was very fortunate that Mary Ellen had the foresight to drive up to the summit well in advance to meet me with a warm car. I thanked her by singing the Golden Girls theme song all the way to the hotel.

B2B '09 Results:
http://www.caldwellcochamber.org/support/pagepics/09Bridge.txt
An account by Bruce Humphries (1st place)
http://dieseldiaries.com/hom/?p=232
Some more blog posts on the '09 ride:

Printing a specific line from bash_history

2009-09-16T12:20:00.000-04:00

Often I want to archive specific commands out of my immediate bash history to a batch script that I can run later.
Unfortunately I could find no way of redirecting !# (where # is the bash history line I wish to save, e.g. !58 executes line 58) to a file. There is the "colon p" option - where !#:p will print the command instead executing it, but I could not redirect or pipe that output either.

So I added this one-liner to a bin directory in my path:


#!/bin/bash
cat $HISTFILE | sed -n "$1p"

I call it getHistLine. So now to save line 58 to a batch script I can just type:
getHistLine 58 >> myBatchScript.sh

To save a range of lines I can type
getHistLine 50,58 >> myBatchScript.sh

Charles Dickens de Bruijn Graph

2009-08-26T23:27:00.000-04:00

Standardized Velvet Assembly Report

2009-08-25T10:17:00.001-04:00

http://code.google.com/p/standardized-velvet-assembly-report/

I finally got my Velvet Assembler report script up on google code. This "program" consists of some short scripts and a Sweave report designed to help Velvet users identify the optimal kmer and cvCut parameters.

GregorianCalendar foolishness

2009-08-13T15:33:00.000-04:00

I wanted to add queries to my Grails application to find tasks that needed to be completed today, or were delinquent (due date < last midnight).
My application kept thinking things were delinquent the afternoon of the due date.

The problem was that I neglected to read how HOUR_OF_DAY differed from HOUR and absentmindedly mixed the two.
http://java.sun.com/j2se/1.4.2/docs/api/java/util/Calendar.html

You can try to set HOUR to 0 but it will default to 12. If it is the afternoon when you request the Calendar object then it will assume you mean 12 noon.



Calendar lastMidnight = Calendar.getInstance();
//DO NOT DO THIS!!
  lastMidnight.set( Calendar.HOUR, lastMidnight.getMinimum(Calendar.HOUR_OF_DAY ));
  ...snip...


  //this is ok
  lastMidnight.set( Calendar.HOUR_OF_DAY, lastMidnight.getMinimum(Calendar.HOUR_OF_DAY ));
  ...snip...


  //this is also ok
  lastMidnight.set( Calendar.AM_PM, Calendar.AM );
  lastMidnight.set( Calendar.HOUR, lastMidnight.getMinimum(Calendar.HOUR ));
  ...snip...

The other settings are pretty self-explanatory


  lastMidnight.set( Calendar.MINUTE, lastMidnight.getMinimum(Calendar.MINUTE));
  lastMidnight.set( Calendar.SECOND, lastMidnight.getMinimum(Calendar.SECOND));
  lastMidnight.set( Calendar.MILLISECOND,lastMidnight.getMinimum(Calendar.MILLISECOND));

Beware! Groovy split and tokenize don't treat empty elements the same

2009-07-28T14:47:00.000-04:00

Groovy's tokenize, which returns a List, will ignore empty elements (when a delimiter appears twice in succession). Split keeps such elements and returns an Array. If you want to use List functions but you don't want to lose your empty elements, then just use split and convert your Array into a List in a separate step.

This might be important if you are parsing CSV files with empty cells.


import groovy.util.GroovyTestCase


class StringTests extends GroovyTestCase {

  protected void setUp() {
    super.setUp()
  }

  protected void tearDown() {
    super.tearDown()
  }

  void testSplitAndTokenize() {
    assertEquals("This,,should,have,five,items".tokenize(',').size(),5)
    assertEquals("This, ,should,have,six,items".tokenize(',').size(),6)

    assertEquals("This, ,should,have,six,items".split(',').size(),6)
    assertEquals("This,,should,have,six,items".split(',').size(),6)

    //convert array to List and re-evaluate
     def fieldArray = "This,,should,have,six,items".split(',')
     def fields=fieldArray.collect{it}
     assert fields instanceof java.util.List
     assertEquals(fields.size(),6)
  }
}

Loading files in grails bootstrap

2009-05-05T12:45:00.000-04:00

I've had to do this for two separate projects.

A web application always has weird ideas of where it is in the path. I tend to look for examples that work without assumptions. This works for the bootstrap scripts I use to start up grails apps. Locating stuff in /web-app in the deployment context is a different story.

I think the best way is as follows:
create a directory: grails-app/conf/resources

For a tab delimited-file I used the following routine to load and create domain objects


class BootStrap {

  def init = {servletContext ->
    //AC204211.3    c0189C22    ACTIVEFIN    UNKNOWN    153071100    153252400    ctg708    9800    191100

    def filePath = "resources/fpc_report.txt"

    def text = ApplicationHolder.application.parentContext.getResource("classpath:$filePath").inputStream.text
    text.eachLine {
      println it
      def BacFields = it.split("\t")
      new Bacs(accession: BacFields[0],
              cloneName: BacFields[1],
              status: BacFields[2],
              chr: BacFields[3],
              chrStart: BacFields[4],
              chrEnd: BacFields[5],
              contig: BacFields[6],
              contigStart: BacFields[7],
              contigEnd: BacFields[8]).save()
    }
     }
}

For a CSV file I used the opencsv library
http://opencsv.sourceforge.net/
I added to my maven pom the following dependency


 
    net.sf.opencsv
    opencsv
    1.8

and used something similar to get the grails appHolder object to give up a Java File


import org.codehaus.groovy.grails.commons.ApplicationHolder

class BootStrap {

  def init = { servletContext ->

    def filePath = "resources/E357Lims.CSV"
    def appHolder=ApplicationHolder.application.parentContext.getResource("classpath:$filePath")

    def reader = new CSVReader(new FileReader(appHolder.getFile()))
    def header = true
    reader.readAll().each{ csvrow ->
      if(!header){
        new FlatReport(name:csvrow[0].trim()).save()
      }
      header = false
    }
}