<?xml version="1.0" encoding="iso-8859-1"?>

<feed xmlns="http://www.w3.org/2005/Atom" xml:lang="en">
<title type="text">Daily Life in an Ivory Basement</title>
<subtitle type="html"><![CDATA[
This Space Intentionally Left Empty
]]></subtitle>
<id>http://ivory.idyll.org/blog/tags/python</id>
<link rel="alternate" type="text/html" href="http://ivory.idyll.org/blog" />
<link rel="self" type="text/xml" href="http://ivory.idyll.org/blog/tags/python" />

<author>
<name>Titus Brown</name>
<uri>http://ivory.idyll.org/blog/tags/python</uri>
<email>titus+blog1@idyll.org</email>
</author>
<rights>Copyright 2004-2009, C. Titus Brown</rights>
<generator uri="http://pyblosxom.sourceforge.net/" version="1.3.2 2/13/2006">
PyBlosxom http://pyblosxom.sourceforge.net/ 1.3.2 2/13/2006
</generator>

<updated>2010-07-07T15:09:20Z</updated>
<!-- icon?  logo?  -->

<entry>
<title type="html">A memory efficient way to remove low-abundance k-mers from large (metagenomic?) DNA data sets</title>
<category term="/u/t/blog/entries/jul-10" />
<id>http://ivory.idyll.org/blog/2010/07/07/kmer-filtering</id>
<updated>2010-07-07T15:09:20Z</updated>
<published>2010-07-07T15:09:20Z</published>
<link rel="alternate" type="text/html" href="http://ivory.idyll.org/blog/jul-10/kmer-filtering" />
<content type="html">&lt;div class=&quot;document&quot;&gt;
&lt;p&gt;I&apos;ve spent the last few weeks working on a simple solution to a
challenging problem in DNA sequence assembly, and I think we&apos;ve got a
nice simple theoretical solution with an actual implementation.  I&apos;d
be interested in comments!&lt;/p&gt;
&lt;div class=&quot;section&quot;&gt;
&lt;h1&gt;&lt;a id=&quot;introduction&quot; name=&quot;introduction&quot;&gt;Introduction&lt;/a&gt;&lt;/h1&gt;
&lt;p&gt;Briefly, the algorithmic challenge is this:&lt;/p&gt;
&lt;p&gt;We have a bunch of DNA sequences in (let&apos;s say) FASTA format,&lt;/p&gt;
&lt;pre class=&quot;literal-block&quot;&gt;
&amp;gt;850:2:1:1145:4509/1
CCGAGTCGTTTCGGAGGGACCCCGCCATCATACTCGGGGAATTCATCTGAAAGCATGATCATAGAGTCACCGAGCA
&amp;gt;850:2:1:1145:4509/2
AGCCAAGAGCACCCCAGCTATTCCTCCCGGACTTCATAACGTAACGGCCTACCTCGCCATTAAGACTGCGGCGGAG
&amp;gt;850:2:1:1145:14575/1
GACGCAAAAGTAATCGTTTTTTGGGAACATGTTTTCATCCTGATCATGTTCCTGCCGATTCTGATCTCGCGACTGG
&amp;gt;850:2:1:1145:14575/2
TAACGGGCTGAATGTTCAGGACCACGTTTACTACCGTCGGGTTGCCATACTTCAACGCCAGCGTGAAAAAGAACGT
&amp;gt;850:2:1:1145:2219/1
GAAGACAGAGTGGTCGAAAGCTATCAGACGCGATGCCTAACGGCATTTTGTAGGGCGTCTGCGTCAGACGCCAACC
&amp;gt;850:2:1:1145:2219/2
GAAGGCGGGTAATGTCCGACAAACATTTGACGTCAAAGCCGGCTTTTTTGTAGTGGGTTTGACTCTTTCGCTTCCG
&amp;gt;850:2:1:1145:5660/1
GATGGCGTCGTCCGGGTGCCCTGGTGGAAGTTGCGGGGATGCGGATTCATCCGGGACGCGCAGACGCAGGCGGTGG
&lt;/pre&gt;
&lt;p&gt;and we want to pre-filter these sequences so that only sequences
containing high-abundance DNA words of length k (&amp;quot;k-mers&amp;quot;),
remain. For example, given a set of DNA sequences, we might want to
remove any sequence that does not contain a k-mer present at least 5
times in the entire data set.&lt;/p&gt;
&lt;p&gt;This is very straightforward to do for small numbers of sequences, or
for small k.  Unfortunately, we are increasingly confronted by data sets
containing 250 million sequences (or more), and we would like to be
able to do this for large k (k &amp;gt; 20, at least).&lt;/p&gt;
&lt;p&gt;You can break the problem down into two basic steps: first,
counting k-mers; and second, filtering sequences based on those k-mer
counts.  It&apos;s not immediately obvious how you would parallelize this
task: the counting should be very quick (basically, it&apos;s I/O
bound) while the filtering is going to rely on wide-reaching lookups
that aren&apos;t localized to any subset of k-mer space.&lt;/p&gt;
&lt;p&gt;tl; dr? we&apos;ve developed a way to do this for arbitrary k, in linear
time and constant memory, efficiently utilizing as many computers as
you have available.  It&apos;s open source and works today, but, uhh, could
use some documentation...&lt;/p&gt;
&lt;/div&gt;
&lt;div class=&quot;section&quot;&gt;
&lt;h1&gt;&lt;a id=&quot;digression-the-bioinformatics-motivation&quot; name=&quot;digression-the-bioinformatics-motivation&quot;&gt;Digression: the bioinformatics motivation&lt;/a&gt;&lt;/h1&gt;
&lt;p&gt;(You can skip this if you&apos;re not interested in the biological motivation ;)&lt;/p&gt;
&lt;p&gt;What we &lt;em&gt;really&lt;/em&gt; want to do with this kind of data is assemble it,
using a &lt;a class=&quot;reference&quot; href=&quot;http://en.wikipedia.org/wiki/De_Bruijn_graph&quot;&gt;De Bruijn graph approach&lt;/a&gt; a la &lt;a class=&quot;reference&quot; href=&quot;http://genome.cshlp.org/content/18/5/821.full&quot;&gt;Velvet&lt;/a&gt;, &lt;a class=&quot;reference&quot; href=&quot;http://www.bcgsc.ca/platform/bioinfo/software/abyss&quot;&gt;ABySS&lt;/a&gt;, or
&lt;a class=&quot;reference&quot; href=&quot;http://soap.genomics.org.cn/soapdenovo.html&quot;&gt;SOAPdenovo&lt;/a&gt;.
However, De Bruijn graphs all rely on building a graph of overlapping
k-mers in memory, which means that their memory usage scales as a
function of the number of unique k-mers.  This is good in general, but
Bad in certain circumstances -- in particular, whenever the data set
you&apos;re trying to assemble has a lot of genomic novelty.  (See &lt;a class=&quot;reference&quot; href=&quot;http://www.ncbi.nlm.nih.gov/pubmed/20211242&quot;&gt;this
fantastic review&lt;/a&gt; and
my &lt;a class=&quot;reference&quot; href=&quot;http://ged.msu.edu/angus/files/lecture5-assembly.pdf&quot;&gt;assembly lecture&lt;/a&gt; for some
background here.)&lt;/p&gt;
&lt;p&gt;One particular circumstance in which De Bruijn graph-based assemblers
flail is in &lt;a class=&quot;reference&quot; href=&quot;http://en.wikipedia.org/wiki/Metagenomics&quot;&gt;metagenomics&lt;/a&gt;, the isolation and
sequencing of genetic material from &amp;quot;the wild&amp;quot;, e.g.  soil or the
human gut.  This is largely because the diversity of bacteria present
in soil is so huge (although the relatively high error rate of
next-gen platforms doesn&apos;t help).&lt;/p&gt;
&lt;p&gt;Prefiltering can help reduce this memory usage by removing erroneous
sequences as well as not-so-useful sequences.  First, any sequence
arising as the result of a sequencing error is going to be extremely
rare, and abundance filtering will remove those.  Second, genuinely
&amp;quot;rare&amp;quot; (shallowly-sequenced) sequences will generally not contribute
much to the assembly, and so removing them is a good heuristic for
reducing memory usage.&lt;/p&gt;
&lt;p&gt;We are currently playing with dozens (probably hundreds, soon) of gigabytes
of metagenomic data, and we really need to do this prefiltering in order
to have a chance at getting a useful assembly out of it.&lt;/p&gt;
&lt;p&gt;It&apos;s worth noting that this is in no way an original thought: in
particular, the Beijing Genome Institute (BGI) did this kind of
prefiltering in their landmark Human Microbiome paper (&lt;a class=&quot;reference&quot; href=&quot;http://www.nature.com/nature/journal/v464/n7285/full/nature08821.html&quot;&gt;ref&lt;/a&gt;).
We want to do it, too, and the BGI wasn&apos;t obliging enough to give
us source code (booooooo hisssssssssssssss).&lt;/p&gt;
&lt;/div&gt;
&lt;div class=&quot;section&quot;&gt;
&lt;h1&gt;&lt;a id=&quot;existing-approaches&quot; name=&quot;existing-approaches&quot;&gt;Existing approaches&lt;/a&gt;&lt;/h1&gt;
&lt;p&gt;Existing approaches are inadequate to our needs, as far as we can tell
from a literature survey and private conversations.  Everyone seems
to rely on big-RAM machines, which are nice if you have them, but shouldn&apos;t
be necessary.&lt;/p&gt;
&lt;p&gt;There are two basic approaches.&lt;/p&gt;
&lt;p&gt;First, you can build a large table in memory, and then map k-mers into
it.  This involves writing a simple hash function that converts DNA
words into numbers.  However, this approach scales poorly with k: for
example, there are 4**k unique DNA sequences of length k (or roughly
(4**k) / 2 + (4**(k/2))/2, considering reverse complements).  So the table
for k=17 needs 4**17 entries -- that&apos;s 17 gb at 1 byte per counter,
which stretches the memory of cheaply available commodity hardware.
And we want to use bigger k than 17 -- most assemblers will be more
effective for longer k, because it&apos;s more specific.  (We&apos;ve been using
k values between 30 and 70 for our assemblies, and we&apos;d like to filter
on the same k.)&lt;/p&gt;
&lt;p&gt;Second, you can observe that k-mer space (for sufficiently large k) is
likely to be quite sparsely occupied -- after all, there&apos;s only so
many actual 30-mers present in a 100gb data set! So, you can do
something clever like use a tree representation of k-mers and then
attach counters to the end nodes of the tree (see, for example,
&lt;a class=&quot;reference&quot; href=&quot;http://www.ncbi.nlm.nih.gov/pubmed/18976482&quot;&gt;tallymer&lt;/a&gt;.  The
problem here that you need to use pointers to connect nodes in the
tree, which means that while the tree size is going to scale with the
amount of novel k-mers -- ok! -- it&apos;s going to have a big constant in
front of it -- bad!.  In our experience this has been prohibitive for
real metagenomic data filtering.&lt;/p&gt;
&lt;p&gt;These seem to be the two dominant approaches, although if you don&apos;t
need to actually &lt;em&gt;count&lt;/em&gt; the k-mers but only assess presence or
absence, you can use something like a &lt;a class=&quot;reference&quot; href=&quot;http://en.wikipedia.org/wiki/Bloom_filter&quot;&gt;Bloom filter&lt;/a&gt; -- for example, see
&lt;a class=&quot;reference&quot; href=&quot;http://bioinformatics.oxfordjournals.org/cgi/content/full/26/13/1595&quot;&gt;Classification of DNA sequences using a Bloom filter&lt;/a&gt;,
which uses Bloom filters to look for novel sequences (the exact
opposite of what we want to do here!).  References to other approaches
welcome...&lt;/p&gt;
&lt;p&gt;Note that you really, really, really want to avoid disk access by
keeping everything in memory.  These are ginormous data sets and we
would like to be able to quickly explore different parameters and
methods of sequence selection.  So we want to come up with a good counting
scheme for k-mers that involves relatively small amounts of memory and
as little disk access as possible.&lt;/p&gt;
&lt;p&gt;I think this is a really fun kind of problem, actually.  The more
clever you are, the more likely you are to come up with a completely
inappropriate data structure, given the amount of data and the basic
algorithmic requirements.  It demands KISS!  Read on...&lt;/p&gt;
&lt;/div&gt;
&lt;div class=&quot;section&quot;&gt;
&lt;h1&gt;&lt;a id=&quot;an-approximate-approach-to-counting&quot; name=&quot;an-approximate-approach-to-counting&quot;&gt;An approximate approach to counting&lt;/a&gt;&lt;/h1&gt;
&lt;p&gt;So, we&apos;ve come up with a solution that scales with the amount of
genomic novelty, and efficiently uses your available memory.  It can
also make effective use of multiple computers (although not multiple
processors).  What is this magic approach?&lt;/p&gt;
&lt;p&gt;&lt;a class=&quot;reference&quot; href=&quot;http://en.wikipedia.org/wiki/Hash_table&quot;&gt;Hash tables&lt;/a&gt;.  Yep.  Map
k-mers into a fixed-size table (presumably one about as big as your
available memory), and have the table entries be counters for k-mer
abundance.&lt;/p&gt;
&lt;p&gt;This is an obvious solution, except for one problem: collisions.  The
big problem with hash tables is that you&apos;re going to have collisions,
wherein multiple k-mers are mapped into a single counting bin.  Oh
noes!  The traditional way to deal with this is to keep track of each
k-mer individually -- e.g. when there&apos;s a collision, use some sort of
data structure (like a linked list) to track the actual k-mers that
mapped to a particular bin.  But now you&apos;re stuck with using gobs of
memory to keep these structures around, because each collision
requires a new pointer of some sort.  It may be possible to get around
this efficiently, but I&apos;m not smart enough to figure out how.&lt;/p&gt;
&lt;p&gt;Instead of becoming smarter, I reconfigured my brain to look at the problem
differently.  (Think Different?)&lt;/p&gt;
&lt;p&gt;The big realization here is that collisions &lt;strong&gt;may not matter&lt;/strong&gt; all
that much.  Consider the situation where you&apos;re filtering on a maximum
abundance of 5 -- that is, you want at least one k-mer in each
sequence to be present five or more times across the data set.  Well,
if you look at the hash bin for a specific k-mer and see the number
&lt;strong&gt;4&lt;/strong&gt;, you immediately know that whether or not there are any
collisions, that particular k-mer isn&apos;t present five or more times,
and can be discarded!  That is, the count for a particular hash bin is
the sum of the (non-negative) individual counts for k-mers mapping to
that bin, and hence that sum is an upper bound on any individual
k-mer&apos;s abundance.&lt;/p&gt;
&lt;img alt=&quot;http://ivory.idyll.org/permanent/kmer-hashtable.png&quot; src=&quot;http://ivory.idyll.org/permanent/kmer-hashtable.png&quot; style=&quot;width: 20%;&quot; /&gt;
&lt;p&gt;The tradeoff is false positives: as k-mer space fills up, the hash
table is going to be hit by more and more collisions.  In turn, more
of the k-mers are going to be falsely reported as being over the
minimum abundance, and more of the sequences will be kept.  You can
deal with this partly by doing iterative filtering with different
prime hash table sizes, but that will be of limited utility if you
have a really saturated hash table.&lt;/p&gt;
&lt;p&gt;Note that the counting and filtering is still O(N), while the false
positives grow as a function of k-mer space occupancy -- which is to
say that the more diversity you have in your data, the more trouble
you&apos;re in.  That&apos;s going to be a problem no matter the approach, however.&lt;/p&gt;
&lt;p&gt;You can see a simple example of approximate and iterative filtering
here, for k=15 (a k-mer space of approximately 1 billion) and hash
table sizes ranging from 50m to 100m.  The curves all approach the
correct final number of reads (which we can calculate exactly, for
this data set) but some take longer than others.  (The code for this
is &lt;a class=&quot;reference&quot; href=&quot;http://github.com/ctb/khmer/blob/master/scripts/ctb-iterative-bench-2.py&quot;&gt;here&lt;/a&gt;.)&lt;/p&gt;
&lt;img alt=&quot;http://ivory.idyll.org/permanent/kmer-filtering-iterative.png&quot; src=&quot;http://ivory.idyll.org/permanent/kmer-filtering-iterative.png&quot; style=&quot;width: 50%;&quot; /&gt;
&lt;/div&gt;
&lt;div class=&quot;section&quot;&gt;
&lt;h1&gt;&lt;a id=&quot;making-use-of-multiple-computers&quot; name=&quot;making-use-of-multiple-computers&quot;&gt;Making use of multiple computers&lt;/a&gt;&lt;/h1&gt;
&lt;p&gt;While I was working out the details of the (really simple) approximate
counting approach, I was bugged by my inability to make effective use
of all the juicy computers to which I have access.  I don&apos;t have many
&lt;em&gt;big&lt;/em&gt; computers, but I do have lots of medium sized computers with
memory in the ~10-20 gb range.  How can I use them?&lt;/p&gt;
&lt;p&gt;This is not a pleasantly parallel problem, for at least two reasons.
First, it&apos;s I/O bound -- reading the DNA sequences in takes more time
than breaking them down into k-mers and counting them.  And since it&apos;s
really memory bound -- you want to use the biggest hash table possible
to minimize collisions
-- it doesn&apos;t seem like using multiple processors on a single machine
will help.  Second, the hash table is going to be too big to
effectively share between computers: 10-20 gb of data per computer is
too much to send over the network.  So what do we do?&lt;/p&gt;
&lt;p&gt;I was busy explaining to &lt;a class=&quot;reference&quot; href=&quot;http://en.wikipedia.org/wiki/Charles_Ofria&quot;&gt;a colleague&lt;/a&gt; that this was
impossible -- always a useful way for me to solve problems ;) -- when
it hit me that you could use &lt;em&gt;multiple chassis&lt;/em&gt; (RAM + CPU + disk) to
decrease the false positive rate with only a small amount of
communication overhead.  Basically, my approach is to partition k-mer
space into Z subsets (one for each chassis) and have each computer count
only the k-mers that fall into their partition.  Then, after the
counting stage, each chassis records a max k-mer abundance per
partition per sequence, and communicates &lt;em&gt;that&lt;/em&gt; to a central
node.  This central node in turn calculates the max k-mer abundance
across all partitions.&lt;/p&gt;
&lt;p&gt;The partitioning trick is a more general form of the specific &apos;prefix&apos;
approach -- that is, separately count and get max abundances/sequence
for all k-mers starting with &apos;A&apos;, then &apos;C&apos;, then &apos;G&apos;, and then &apos;T&apos;.
For each sequence you will then have four values (the max
abundance/sequence for k-mers start with &apos;A&apos;, &apos;C&apos;, &apos;G&apos;, and &apos;T&apos;),
which can be cheaply stored on disk or in memory.  Now you can do a
single-pass integration and figure out what sequences to keep.&lt;/p&gt;
&lt;p&gt;This approach effectively multiplies your working
memory by a factor of Z, decreasing your false positive rate
equivalently - no mean feat!&lt;/p&gt;
&lt;img alt=&quot;http://ivory.idyll.org/permanent/kmer-hashtable-par.png&quot; src=&quot;http://ivory.idyll.org/permanent/kmer-hashtable-par.png&quot; style=&quot;width: 20%;&quot; /&gt;
&lt;img alt=&quot;http://ivory.idyll.org/permanent/kmer-filter-process-2.png&quot; src=&quot;http://ivory.idyll.org/permanent/kmer-filter-process-2.png&quot; style=&quot;width: 40%;&quot; /&gt;
&lt;p&gt;The communication load is significant but not prohibitive: replicate a
read-only sequence data set across Z computers, and then communicate
max values (1 byte for each sequence) back -- 50-500 mb of data per
filtering round.  The result of each filtering round can be returned
to each node as a readmask against the already-replicated initial
sequence set, with one bit per sequence for ignore or process.  You can
even do it on a single computer, with a local hard drive, in multiple
iterations.&lt;/p&gt;
&lt;p&gt;You can see a simple in-memory implementation of this approach &lt;a class=&quot;reference&quot; href=&quot;http://github.com/ctb/khmer/blob/master/python/khmer/split.py&quot;&gt;here&lt;/a&gt;,
and some tests &lt;a class=&quot;reference&quot; href=&quot;http://github.com/ctb/khmer/blob/master/python/test_split.py&quot;&gt;here&lt;/a&gt;.
I&apos;ve implemented this using readmask tables and min/max tables (serializable
data structures) more generally, too; see &amp;quot;the actual code&amp;quot;, below.&lt;/p&gt;
&lt;/div&gt;
&lt;div class=&quot;section&quot;&gt;
&lt;h1&gt;&lt;a id=&quot;similar-approaches&quot; name=&quot;similar-approaches&quot;&gt;Similar approaches&lt;/a&gt;&lt;/h1&gt;
&lt;p&gt;By allowing for false positives, I&apos;ve effectively turned the hash
table into a probabilistic data structure.  The closest analog I&apos;ve
seen is the &lt;a class=&quot;reference&quot; href=&quot;http://en.wikipedia.org/wiki/Bloom_filter&quot;&gt;Bloom filter&lt;/a&gt; which records
presence/absence information using multiple hash functions for
arbitrary k.  The hash approach outlined above devolves into a
maximally efficient Bloom filter for fixed k when only
presence/absence information is recorded.&lt;/p&gt;
&lt;/div&gt;
&lt;div class=&quot;section&quot;&gt;
&lt;h1&gt;&lt;a id=&quot;the-actual-code&quot; name=&quot;the-actual-code&quot;&gt;The actual code&lt;/a&gt;&lt;/h1&gt;
&lt;p&gt;Theory and practice are the same, in theory.  In practice, of course,
they differ.  A whole host of minor interface and implementation
design decisions needed to be made.  The result can be seen at the
&apos;khmer&apos; project, here: &lt;a class=&quot;reference&quot; href=&quot;http://github.com/ctb/khmer&quot;&gt;http://github.com/ctb/khmer&lt;/a&gt;.  It&apos;s slim on
documentation, but there are some example scripts and a reasonable
amount of tests.  It requires nothing but Python 2.6 and a compiler;
nose if you want to run the tests.&lt;/p&gt;
&lt;p&gt;The implementation is in C++ with a Python wrapper, much like the
paircomp and motility software packages.&lt;/p&gt;
&lt;p&gt;It will filter 1m 70-mer sequences in about 45 seconds, and 80 million sequences
in less than an hour, on a 3 GHz Xeon with 16 gbs of RAM.&lt;/p&gt;
&lt;p&gt;Right now it&apos;s limited to k &amp;lt;= 32, because we encode each DNA k-mer in
a 64-bit &apos;long long&apos;.&lt;/p&gt;
&lt;p&gt;You can see an exact filtering script here:
&lt;a class=&quot;reference&quot; href=&quot;http://github.com/ctb/khmer/blob/master/scripts/filter-exact.py&quot;&gt;http://github.com/ctb/khmer/blob/master/scripts/filter-exact.py&lt;/a&gt; .  By
varying the hash table size (second arg to new_hashtable) you can turn
it into an &lt;em&gt;inexact&lt;/em&gt; filtering script quite easily.&lt;/p&gt;
&lt;p&gt;Drop me a note if you want help using the code, or a specific example.
We&apos;re planning to write documentation, doctests, examples, robust
command line scripts, etc. prior to publication, but for now we&apos;re
primarily trying to use it.&lt;/p&gt;
&lt;/div&gt;
&lt;div class=&quot;section&quot;&gt;
&lt;h1&gt;&lt;a id=&quot;other-uses&quot; name=&quot;other-uses&quot;&gt;Other uses&lt;/a&gt;&lt;/h1&gt;
&lt;p&gt;It has not escaped our notice that you can use this approach for a bunch of
other k-mer based problems, such as repeat discovery and calculating abundance
distributions... but right now we&apos;re focused on actually using it for
filtering metagenomic data sets prior to assembly.&lt;/p&gt;
&lt;/div&gt;
&lt;div class=&quot;section&quot;&gt;
&lt;h1&gt;&lt;a id=&quot;acks&quot; name=&quot;acks&quot;&gt;Acks&lt;/a&gt;&lt;/h1&gt;
&lt;p&gt;I talked a fair bit with Prof. Rich Enbody about this approach, and he
did a wonderful job of double-checking my intuition.  Jason Pell and
Qingpeng Zhang are co-authors on this project; in particular, Jason
helped write the C++ code, and Qingpeng has been working with k-mers
in general and tallymer in specific on some other projects, and
provided a lot of background help.&lt;/p&gt;
&lt;/div&gt;
&lt;div class=&quot;section&quot;&gt;
&lt;h1&gt;&lt;a id=&quot;conclusions&quot; name=&quot;conclusions&quot;&gt;Conclusions&lt;/a&gt;&lt;/h1&gt;
&lt;p&gt;We&apos;ve taken a previously hard problem and made it practically
solvable: we can filter ~88m sequences in a few hours on a
single-processor computer with only 16gb of RAM.  This seems useful.&lt;/p&gt;
&lt;p&gt;Our main challenge now is to come up with a better hashing function.
Our current hash function is not uniform, even for a
uniform distribution of k-mers, because of the way it handles reverse
complements.&lt;/p&gt;
&lt;p&gt;The approach scales reasonably nicely.  Doubling the amount of data
doubles the compute time.  However, if you have double the novelty,
you&apos;ll need to do double the partitions to keep the same false
positive rate, in which case you quadruple the compute time.  So it&apos;s
O(N^2) for the worst case (unending novelty) but should be something
better for real-world cases.  That&apos;s what we&apos;ll be looking at over
the next few months.&lt;/p&gt;
&lt;p&gt;I haven&apos;t done enough background reading to figure out if our approach
is particularly novel, although in the space of bioinformatics it seems
to be reasonably so.  That&apos;s less important than actually solving our
problem, but it would be nice to punch the &amp;quot;publication&amp;quot; ticket if possible.
We&apos;re thinking of writing it up and sending it to BMC Bioinformatics,
although suggestions are welcome.&lt;/p&gt;
&lt;p&gt;It would be particularly ironic if the first publication from my lab
was this computer science-y, given that I have no degrees in CS and
am in the CS department by kind of a fluke of the hiring process ;).&lt;/p&gt;
&lt;/div&gt;
&lt;/div&gt;
</content>
</entry>

<entry>
<title type="html">Teaching scientists how to use computers - hub &amp; spokes</title>
<category term="/u/t/blog/entries/jul-10" />
<id>http://ivory.idyll.org/blog/2010/07/05/swc-hub-spokes</id>
<updated>2010-07-06T03:12:24Z</updated>
<published>2010-07-06T03:12:24Z</published>
<link rel="alternate" type="text/html" href="http://ivory.idyll.org/blog/jul-10/swc-hub-spokes" />
<content type="html">&lt;div class=&quot;document&quot;&gt;
&lt;p&gt;After my recent &lt;a class=&quot;reference&quot; href=&quot;http://ivory.idyll.org/blog/jun-10/ngs-course-postmortem&quot;&gt;next-gen sequencing course&lt;/a&gt;, which
was supposed to tie into the whole &lt;a class=&quot;reference&quot; href=&quot;http://software-carpentry.org/&quot;&gt;software carpentry (SWC) effort&lt;/a&gt; but didn&apos;t really succeed in doing
so the first time through, I started thinking about the Right Way to
tie in the SWC material.  In particular, how do you both motivate
scientists to look at the SWC material, and (re)direct people to the
appropriate places?&lt;/p&gt;
&lt;p&gt;It&apos;s not clear that a Plan is in place.  Greg Wilson seems to assume
that scientists will find at least some of the material immediately
obviously usable, but I think he&apos;s targetted at a more sophisticated
population of users -- physicists and the like.  My experience with
bioinformaticians, however, is that they either come from straight
biology backgrounds (with little or no computational background and
rather limited on-the-job training), straight computation backgrounds
(with very little biology), or physics (gonzo programming skills, but
no biology).  The latter fit neatly into the SWC fold, but they (we ;)
are rare in biology.  I think computer scientists and biologists are
going to need guidance to dive into SWC at an early enough time for it
to be the most rewarding.&lt;/p&gt;
&lt;p&gt;So, what&apos;s a good model for SWC to guide scientists from multiple
disciplines into the appropriate material?  It&apos;s obviously not going
to be possible to have Greg et al. tailor the SWC material to individual
subgroups -- he doesn&apos;t know much (any ;) biology, for example.  I don&apos;t
have the time, patience, or skillset to integrate my next-gen notes
into his SWC material, either.  So, instead, I propose the hub &amp;amp; spokes
model!&lt;/p&gt;
&lt;img alt=&quot;http://ivory.idyll.org/permanent/hub-spokes.png&quot; src=&quot;http://ivory.idyll.org/permanent/hub-spokes.png&quot; /&gt;
&lt;p&gt;Here, the &amp;quot;hub&amp;quot; is the SWC material, and the spokes are all of the
individual disciplines.&lt;/p&gt;
&lt;p&gt;Basically, the idea is that individual sites (like my own ANGUS site
on next-gen sequencing, &lt;a class=&quot;reference&quot; href=&quot;http://ged.msu.edu/angus/&quot;&gt;http://ged.msu.edu/angus/&lt;/a&gt;) will develop their
own field-specific content, and then link from that content into the
SWC notes.  This way the experts with feet in both fields can link
appropriately, and Greg only has to worry about making the central
content general -- which he&apos;s already doing quite well, I think.  Yes,
It&apos;s more work than asking Greg to do it, but frankly I&apos;m going to be
happy with a kick-ass central SWC site to which I can link -- right
now it&apos;s dismayingly challenging to teach students why this stuff
matters and how to learn it.&lt;/p&gt;
&lt;p&gt;From the psychosocial perspective, it&apos;s a great fit.  Students can get
hands on tutorials on how to do X, Y, and Z in their own field -- and
then connect into the SWC material to learn the background, or
additional computational techniques in support of it.  Motivation first!&lt;/p&gt;
&lt;p&gt;What do we need SWC to do to support this?  Not much -- basically, the
central SWC notes need to be stable enough (with permalinks) that I
can link into them from my own site(s) and not have to worry about the
links becoming broken or (worse) silently migrating in topic.  There
are other solutions (wholesale incorporation of SWC into my own notes,
for example) but I think the permalink idea is the most
straightforward.  Oh, and we should have a Greg-gets-hit-by-a-bus plan,
too; at some point he&apos;s going to move on from SWC (perhaps when his
lovely wife decide she&apos;s had enough and he needs to stop obsessing over
it, or perhaps under more dire circumstances ;( and it would be good to
know who holds the domain and site keys.&lt;/p&gt;
&lt;p&gt;Thoughts?  Comments?&lt;/p&gt;
&lt;p&gt;--titus&lt;/p&gt;
&lt;/div&gt;
</content>
</entry>

<entry>
<title type="html">Which functional programming language(s) should we teach?</title>
<category term="/u/t/blog/entries/jun-10" />
<id>http://ivory.idyll.org/blog/2010/06/24/functional-programming-languages</id>
<updated>2010-06-24T18:31:48Z</updated>
<published>2010-06-24T18:31:48Z</published>
<link rel="alternate" type="text/html" href="http://ivory.idyll.org/blog/jun-10/functional-programming-languages" />
<content type="html">&lt;div class=&quot;document&quot;&gt;
&lt;p&gt;Laurie Dillon just posted the SIGPLAN eduction board article on &lt;a class=&quot;reference&quot; href=&quot;http://mt4.acm.org/educationboard/2010/06/why-undergraduates-should-learn-the-principles-of-programming-languages.html&quot;&gt;Why
Undergraduates Should Learn the Principles of Programming Languages&lt;/a&gt;
to our faculty mailing list at the &lt;a class=&quot;reference&quot; href=&quot;http://www.cse.msu.edu&quot;&gt;MSU Computer Science department&lt;/a&gt;.  One question that came up in the ensuing
conversation was: what functional programming language(s) would/should
we teach?&lt;/p&gt;
&lt;p&gt;I mentioned OCaml, Haskell, and Erlang as reasonably pure but still
pragmatic FP languages.  Anything else that&apos;s both &amp;quot;truly&amp;quot; functional
and used somewhat broadly in the real world?&lt;/p&gt;
&lt;p&gt;thanks!&lt;/p&gt;
&lt;p&gt;--titus&lt;/p&gt;
&lt;/div&gt;
</content>
</entry>

<entry>
<title type="html">Teaching next-gen sequencing data analysis to biologists</title>
<category term="/u/t/blog/entries/jun-10" />
<id>http://ivory.idyll.org/blog/2010/06/14/ngs-course-postmortem</id>
<updated>2010-06-14T15:38:31Z</updated>
<published>2010-06-14T15:38:31Z</published>
<link rel="alternate" type="text/html" href="http://ivory.idyll.org/blog/jun-10/ngs-course-postmortem" />
<content type="html">&lt;div class=&quot;document&quot;&gt;
&lt;p&gt;Our &lt;a class=&quot;reference&quot; href=&quot;http://bioinformatics.msu.edu/ngs-summer-course-2010&quot;&gt;sequencing analysis course&lt;/a&gt; ended last
Friday, with an overwhelmingly positive response from the students.
The few negative comments that I got were largely about organizational
issues, and could be reshaped as suggestions for next time rather than
as condemnations of this year&apos;s course.&lt;/p&gt;
&lt;img alt=&quot;http://ivory.idyll.org/permanent/ngs-2010-group.png&quot; src=&quot;http://ivory.idyll.org/permanent/ngs-2010-group.png&quot; style=&quot;width: 65%;&quot; /&gt;
&lt;p&gt;The 23 students -- most with no prior command-line experience -- spent
two weeks experiencing at first hand the challenges of dealing with
dozens of gigabytes of sequencing data.  Each of the students went
through genome-scale mapping, genome assembly, mRNAseq analysis on an
&amp;quot;emerging model organism&amp;quot; (a.k.a &amp;quot;one with a crappy genome&amp;quot;, lamprey),
resequencing analysis on E. coli, and ChIP-seq analysis on Myxococcus
xanthus.  By the beginning of the second week, many students were working
with their own data -- a real victory.  Python programming competency
may take a bit longer, but many of them seem motivated.&lt;/p&gt;
&lt;p&gt;If you had told me three weeks ago that we could pull this off, I
would have told you that you were crazy.  This does beg the question
of what I was thinking when I proposed the course -- but don&apos;t dwell
on that, please...&lt;/p&gt;
&lt;p&gt;The locale was great, as you can see:&lt;/p&gt;
&lt;img alt=&quot;http://ivory.idyll.org/permanent/ngs-2010-beach.png&quot; src=&quot;http://ivory.idyll.org/permanent/ngs-2010-beach.png&quot; style=&quot;width: 65%;&quot; /&gt;
&lt;p&gt;One of the most important lessons of the course for me is that &lt;a class=&quot;reference&quot; href=&quot;http://ivory.idyll.org/blog/jun-10/ngs-course-with-aws.html&quot;&gt;cloud
computing works well to backstop this kind of course&lt;/a&gt;.  I
was very worried about the suitabiliy and reliability and ease of use,
but AWS did a great job, providing an easy-to-use Web interface and a
good range of machine images.  I have little doubt that this course
would have been nearly impossible (and either completely ineffective
or much more expensive) without it.&lt;/p&gt;
&lt;p&gt;In the end, we spent more on beer than on computational power.  That says
something important to me :)&lt;/p&gt;
&lt;p&gt;The &lt;a class=&quot;reference&quot; href=&quot;http://ged.msu.edu/angus/&quot;&gt;course notes&lt;/a&gt; are available under a
CC license although they need to be reworked to use publicly available
data sets before they become truly useful.  At that point I expect them
to become awesomely useful, though.&lt;/p&gt;
&lt;p&gt;From the scientific perspective, the students derived a number of
significant benefits from the course.  One that I had not really
expected was that some students had no idea what went in to
computational &amp;quot;sausage&amp;quot;, and were kind of shocked to see what kinds of
assumptions us comp bio people made on their behalf.  This was
especially true in the case of students from companies, who have
pipelines that are run on their data.  One student lamented that &amp;quot;we
used to look at the raw traces... now all we get are spreadsheet
summaries!&amp;quot;  Another student came to me in a panic because they didn&apos;t
realize that there &lt;em&gt;was&lt;/em&gt; no one true answer -- that that was in fact
part of the &amp;quot;fun&amp;quot; of &lt;em&gt;all&lt;/em&gt; biology, not just experimental biology.
These reactions alone made teaching the course worthwhile.&lt;/p&gt;
&lt;p&gt;Of course, the main point is that many of the students seem to be
capable of at least starting their own analyses now.  I was surprised
at the practical power of our cut-and-paste approach -- for example,
if you look at the &lt;a class=&quot;reference&quot; href=&quot;http://ged.msu.edu/angus/tutorials/short-read-assembly.html&quot;&gt;Short-read assembly with ABySS tutorial&lt;/a&gt;, it
turns out to be relatively straightforward to adapt this to doing
assemblies of your own genomic or transcriptomic data.  I based our
approach on Greg Wilson&apos;s post on &lt;a class=&quot;reference&quot; href=&quot;http://pyre.third-bit.com/blog/archives/3761.html&quot;&gt;the failure of inquiry-based
teaching&lt;/a&gt; and so
far I like it.&lt;/p&gt;
&lt;p&gt;I am particularly amused that we have now documented, in replicable
detail, the Kroos Lab MrpC ChIP analysis.  We also have the best
documentation for Jeff Barrick&apos;s breseq software, I think; this is
what is used to analyze the &lt;a class=&quot;reference&quot; href=&quot;http://en.wikipedia.org/wiki/E._coli_long-term_evolution_experiment&quot;&gt;Long Term Evolution Experiment&lt;/a&gt;
lines -- and I can&apos;t wait for the anti-evolutionists to pounce on
that...  &amp;quot;Titus Brown -- making evolution experiments accessible to
creationists.&amp;quot;  Yay?&lt;/p&gt;
&lt;p&gt;There were a number of problems and mistakes that we had to
steamroller through.  In particular, more background and more advanced
tutorials would have be great, but we just didn&apos;t have time to write
them.  Some 454, Helicos, and SOLiD data sets (and next year, PacBio?)
would be a good addition.  We had a general lack of multiplexing data,
which is becoming a Big Thing now that sequencing is so ridiculously
deep. I would also like to introduce additional real data analyses
next year, reprising things like the &lt;a class=&quot;reference&quot; href=&quot;http://www.nature.com/nbt/journal/v28/n5/abs/nbt.1621.html&quot;&gt;Cufflinks analysis&lt;/a&gt; and
whole-vertebrate-genome ChIP-seq/mRNAseq &lt;a class=&quot;reference&quot; href=&quot;http://www.nature.com/nmeth/journal/v6/n11s/abs/nmeth.1371.html&quot;&gt;a la the Wold Lab&lt;/a&gt;.
I&apos;m weighing adding metagenomics data analysis in for a day, although
it&apos;s a pretty separate field of inquiry (and frankly much harder in
terms of &amp;quot;unknown unknowns&amp;quot;).  We also desperately need some plant
genomics expertise, because frankly I know nothing about plant
genomes; my last-minute plant genomics TA fell through due to lack of
planning on my part.  (Conveniently, plant genomics is something MSU
is particularly good at, so I&apos;m sure I can find someone next year.)&lt;/p&gt;
&lt;p&gt;Oops, did I say next year?  Well, yes.  &lt;em&gt;If&lt;/em&gt; I can find funding for my
princely salary, &lt;em&gt;then&lt;/em&gt; I will almost certainly run the course again
next year.  I can cover TAs and my own room/board and speakers with
workshop fees, but if I&apos;m going to keep room+board+fees under
$1000/student -- a practical necessity for most -- there&apos;s no way I
can pay myself, too.  And while this year I relied on my lovely,
patient, and frankly long-suffering wife to hold down the home fort
while I was away for two weeks, I simply can&apos;t put her through that
again, so I will need to pay for a nanny next year.  So doing it for
free is not an option.&lt;/p&gt;
&lt;p&gt;In other words, &lt;strong&gt;if you are a sequencing company, or an NIH/NSF/USDA
program director, interested in keeping this going, please get in
touch&lt;/strong&gt;.  I plan to apply for this &lt;a class=&quot;reference&quot; href=&quot;http://grants.nih.gov/grants/guide/pa-files/PAR-09-245.html&quot;&gt;Initiative to Maximize Research
Education in Genomics&lt;/a&gt; in
September, but I am not confident of getting that on the first try,
and in any case I will need letters of support from interested folks.
So &lt;a class=&quot;reference&quot; href=&quot;mailto:ctb&amp;#64;msu.edu&quot;&gt;drop me a note at ctb&amp;#64;msu.edu&lt;/a&gt;.&lt;/p&gt;
&lt;p&gt;Course development this year was sponsored by the &lt;a class=&quot;reference&quot; href=&quot;http://www.bch.msu.edu/GEDD/&quot;&gt;MSU Gene Expression
in Disease and Development&lt;/a&gt;, to whom
I am truly grateful.  The course would simply not have been possible
without their support.&lt;/p&gt;
&lt;p&gt;My overall conclusion is that it is possible to teach bench biologists
with no prior computational experience to achieve at least minimal
competency in real-world data analysis of next-generation sequencing
data.  I can&apos;t conclusively &lt;em&gt;demonstrate&lt;/em&gt; this without doing a better
job of course evaluation, and of course only time will tell if it
sticks for any of the students, but right now I&apos;m feeling pretty good
about the course overall.  Not to mention massively relieved.&lt;/p&gt;
&lt;p&gt;--titus&lt;/p&gt;
&lt;p&gt;p.s. Update from one student -- &amp;quot;It&apos;s not even 12 o&apos;clock Monday
morning and I&apos;m already getting people asking me how to run assemblies
and analyze data.&amp;quot;  Heh.&lt;/p&gt;
&lt;/div&gt;
</content>
</entry>

<entry>
<title type="html">Running a next-gen sequence analysis course using Amazon Web Services</title>
<category term="/u/t/blog/entries/jun-10" />
<id>http://ivory.idyll.org/blog/2010/06/08/ngs-course-with-aws</id>
<updated>2010-06-08T14:52:29Z</updated>
<published>2010-06-08T14:52:29Z</published>
<link rel="alternate" type="text/html" href="http://ivory.idyll.org/blog/jun-10/ngs-course-with-aws" />
<content type="html">&lt;div class=&quot;document&quot;&gt;
&lt;p&gt;So, I&apos;ve been teaching a &lt;a class=&quot;reference&quot; href=&quot;http://bioinformatics.msu.edu/ngs-summer-course-2010&quot;&gt;course&lt;/a&gt; on
next-generation sequence analysis for the last week, and one of the
issues I had to deal with before I proposed the course was how to deal
with the &lt;strong&gt;volume of data&lt;/strong&gt; and the required computation.&lt;/p&gt;
&lt;p&gt;You see, next-generation sequence analysis involves analyzing not just
entire genomes (which are, after all, only 3gb or so in size) but data
sets that are 100x or 1000x as big!  We want to not just map these
data sets (which is CPU-intensive), but also perform memory-intensive
steps like assembly.  If you have a class with 20+ students in it, you
need to worry about a lot of things:&lt;/p&gt;
&lt;blockquote&gt;
&lt;ul class=&quot;simple&quot;&gt;
&lt;li&gt;computational power: how do you provide 24 &amp;quot;good&amp;quot; workstations&lt;/li&gt;
&lt;li&gt;memory&lt;/li&gt;
&lt;li&gt;disk space&lt;/li&gt;
&lt;li&gt;bandwidth&lt;/li&gt;
&lt;li&gt;&amp;quot;take home&amp;quot; ability&lt;/li&gt;
&lt;/ul&gt;
&lt;/blockquote&gt;
&lt;p&gt;One strategy would be to simply provide some Linux or Mac
workstations, with cut-down data sets.  But then you wouldn&apos;t be
teaching reality -- you&apos;d be teaching a cut-down version of reality.
This would make the course particularly irrelevant given that one of
the extra-fun things about next-gen sequence analysis is how hard it
is to deal with the volume of data.  You also have to worry that the
course would be made even &lt;em&gt;more&lt;/em&gt; irrelevant because the students would
leave the course and be unable to use the information without finding
infrastructure and installing a bunch of software and then administering
the machine.&lt;/p&gt;
&lt;p&gt;While I enjoy setting up computers and installing software and
managing users, I&apos;m clearly masochistic.  It&apos;s also &lt;em&gt;entirely besides
the point&lt;/em&gt; for bioinformaticians and biologists - they just want to
analyze data!&lt;/p&gt;
&lt;p&gt;The solution I came up with was to use Amazon Web Services and rent
some EC2 machines.  There&apos;s a large variety of hardware configurations
available (see &lt;a class=&quot;reference&quot; href=&quot;http://aws.amazon.com/ec2/#instance&quot;&gt;instance types&lt;/a&gt;) and they&apos;re not that
expensive per hour (see &lt;a class=&quot;reference&quot; href=&quot;http://aws.amazon.com/ec2/#pricing&quot;&gt;pricing&lt;/a&gt;).&lt;/p&gt;
&lt;p&gt;&lt;strong&gt;This has worked out really, really well.&lt;/strong&gt;&lt;/p&gt;
&lt;p&gt;It&apos;s hard to enumerate the benefits, because there have been so many
of them ;).  A few of the obvious ones --&lt;/p&gt;
&lt;p&gt;We&apos;ve been able to write tutorials (temporary home here:
&lt;a class=&quot;reference&quot; href=&quot;http://ged.msu.edu/angus/&quot;&gt;http://ged.msu.edu/angus/&lt;/a&gt;) that make use of specific images and should
be as future-proof as they can be.  We&apos;ve given students cut and paste
command lines that Just Work, and that they can tweak and modify as
they want.  If it borks, they always just throw it away and start from
a clean install.&lt;/p&gt;
&lt;p&gt;It&apos;s dirt cheap.  We spent less than $50 the first week, for ~30
people using an average of 8 hours of CPU time.  The second week will
increase to an average of 8 hours of CPU time a day, and for larger
instances -- so probably about $300 total, or maybe even $500 -- but
that&apos;s ridiculously cheap, frankly, when you consider that there are
no hardware issues or OS re-install problems to deal with!&lt;/p&gt;
&lt;p&gt;Students can choose whatever machine specs they need in order to do
their analysis.  More memory?  Easy.  Faster CPU needed?  No problem.&lt;/p&gt;
&lt;p&gt;All of the data analysis takes place off-site.  As long as we can provide
the data sets somewhere else (I&apos;ve been using S3, of course) the students
don&apos;t need to transfer multi-gigabyte files around.&lt;/p&gt;
&lt;p&gt;The students can go home, rent EC2 machines, and do their own analyses
-- without their labs buying any required infrastructure.&lt;/p&gt;
&lt;p&gt;Home institution computer admins can use the EC2 tutorials as
documentation to figure out what needs to be installed (and
potentially, maintained) in order for their researchers to do next-gen
sequence analysis.&lt;/p&gt;
&lt;p&gt;The documentation should even serve as a general set of tutorials,
once I go through and remove the dependence on private data sets!
There won&apos;t be any need for students to do difficult or tricky configurations
on their home machines in order to make use of the tutorial info.&lt;/p&gt;
&lt;p&gt;So, truly awesome.  I&apos;m going to be using it for all my courses from now
on, I think.&lt;/p&gt;
&lt;p&gt;There have been only two minor hitches.&lt;/p&gt;
&lt;p&gt;First, I&apos;m using Consolidated Billing to pay for all of the students&apos;
computer use during the class, and Amazon has some rules in place to
prevent abuse of this.  They&apos;re limiting me to 20 consolidated billing
accounts per AWS account, which means that I&apos;ve needed to get a second
AWS account in order to add all 30 students, TAs, and visiting
instructors.  I wouldn&apos;t even mention it as a serious issue but for
the fact that &lt;em&gt;they don&apos;t document it anywhere&lt;/em&gt;, so I ran into this on
the first day of class and then had to wait for them to get back to me
to explain what was going on and how to work around it.  Grr.&lt;/p&gt;
&lt;p&gt;Second, we had some trouble starting up enough Large instances
simultaneously on the day we were doing assembly.  Not sure what that
was about.&lt;/p&gt;
&lt;p&gt;Anyway, so I give a strong +1 on Amazon EC2 for large-ish style data
analysis.  Good stuff.&lt;/p&gt;
&lt;p&gt;cheers,
--titus&lt;/p&gt;
&lt;/div&gt;
</content>
</entry>
</feed>
