<![CDATA[Benchling Engineering]]>Benchling.com]]>http://benchling.engineering/Ghost 0.11Mon, 27 Feb 2017 05:37:22 GMT60<![CDATA[Building a Regex Search Engine for DNA]]>This winter, Tony interned at Benchling and built the foundation of our DNA search feature. Somak and I brought it into production shortly after his internship ended, and we're finally writing it up!

We recently launched DNA search on Benchling - you can have a library of thousands of plasmids

http://benchling.engineering/dna-regex-search/74d5fc59-861d-44f3-aa39-b52bfddf4da1Wed, 20 Apr 2016 05:59:00 GMTThis winter, Tony interned at Benchling and built the foundation of our DNA search feature. Somak and I brought it into production shortly after his internship ended, and we're finally writing it up!

We recently launched DNA search on Benchling - you can have a library of thousands of plasmids and primers, and we'll search through them in less than 100ms end to end. It was a fun project, and through it we ended up building a regex search engine (on top of Elasticsearch) to support the various search needs of scientists. Here's how we did it.


Some rough biology: each DNA strand is modeled as a string of A, T, C, and G characters, and each character is called a base. For a string like ATG AAA TTT the DNA gets read 3 bases at a time. 3 bases together are called a codon, and each codon eventually gets translated into a series of amino acids: ATG becomes methionine, AAA becomes lysine, and TTT becomes phenylalanine. Amino acids then are assembled into proteins, which then perform an astounding number of functions inside an organism. In this way, DNA stores the "code" of life.

DNA search is essentially substring search across a database of strings (albeit with a smaller alphabet). Some common use cases: one, scientists will search for certain genes that they've used in engineered plasmids. Two, since multiple codons can translate to the same amino acid, a process called codon optimization might replace codons with equivalent ones that work better given a certain kind of organism - your CAA (glutamine) might become CAG (still glutamine) instead. You might want to search by amino acid ("glutamine") so that you find it regardless of whether CAA or CAG was used.

As scientists design and document their research on Benchling, the constructs that they create build up over time. DNA search should help scientists cut through the noise and find the relevant sequences they're looking for.

General Approach

Scientists have a few ways of describing sequences, all of which map well to regular expressions:

  • IUPAC notation describes DNA bases. On top of A, C, T, and G, you have letters like N (any base: A, C, G, or T) and W (A or T). You can search for TAGNGG, which becomes TAG(A|C|T|G)GG. Here's a chart from Wikipedia:
  • Amino acids have their own set of letters: glutamine is Q, for example. Scientists can search by amino acids, and multiple codons will translate to the same amino acid - Q maps to the regex (CAA|CAG). Searching for QWD becomes (CAA|CAG)(TGG)(GAT|GAC)
  • Finally, amino acids have a regex-like syntax introduced by the protein database Prosite - Q(3)-x-x denotes three glutamines Q followed by any two amino acids x, or (CAA|CAG)(CAA|CAG)(CAA|CAG)(NNN)(NNN).

To address these three cases, we built a "regex search engine" - users can use IUPAC notation, amino acid notation, or prosite notation and we'll translate it into a regex to feed into the engine.

Our approach uses a combination of regex parsing and Elasticsearch indexing:

  • We build an inverted trigram index of DNA bases
  • For each regex, we convert it into a set of constraints that we can use as filter. The filter can allow for false positives (but not false negatives) - it helps narrow the search space: (A|C|G|T)GG becomes AGG OR CGG OR GGG OR TGG
  • That filter is sent to Elasticsearch, which uses the trigram index to quickly find potentially matching bases. We use a post-search regex filter to eliminate false positives
  • That's it! The result is sub 100ms searches across a user's entire collection of sequences.

We used Elasticsearch because we already used it for indexing and searching the various file types we have on Benchling. Searching for bases is just another way to match items and combines well with our existing filters and scoring functions. (We considered building it using pg_trgm / Postgres instead, but Elasticsearch has Postgres beat in terms of built-in analyzers and query features.)

We're definitely not the first to build regex search, and we won't be the last - Russ Cox has a great writeup of how it was done for Google code search. What differentiates our approach is that we built it for DNA search using relatively standard software bits.

Finally, a quick note on BLAST (https://en.wikipedia.org/wiki/BLAST): it and several other similar algorithms have been developed for quickly finding similar sequences given an input sequence. The exact use case here is different, though: BLAST is based on a heuristic algorithm to approximately match against a large set of sequences, whereas here we're focused on precise matching of a scientist's personal library. The narrower scope allows us to be more precise and support more complex regexes, but as you'll see it still requires some work to be performant.

Parsing Regex

It turned out finding a Python library to parse regexes (and not parse using regexes) was surprisingly difficult. You could use the debug flag (re.compile(regex, re.DEBUG)) to see the AST it parses, but 1) it prints to stdout and 2) we'd need to reparse it again to get our own AST, so it seemed too hacky of an approach:

>>> re.compile('a(b.)*', re.DEBUG)
literal 97  
max_repeat 0 4294967295  
  subpattern 1
    literal 98
    any None

(Update: Andrew Dalke had a great tip - there's actually an undocumented sre_parse module that returns a parse tree - definitely worth looking into.)

We ended up writing a simple recursive descent parser to parse our subset of supported regex. It sounds fancier than it is: Matt Might has a great article on writing one that came in handy. The nice thing about recursive descent parsers is that they're pretty simple to write for a language like regex. I wrote the core parser on a plane ride - it weights in at about 200 lines of code and easy to test. So we started off with that as a proof of concept.

parse(...) takes in the regex string to parse and outputs a data structure commonly referred to as an abstract syntax tree (AST):

>>> parse('a(b.)*')
ArrayNode([CharNode('a'), StarNode(ArrayNode([CharNode('b'), PeriodNode()]), False)])  

Our AST generally consists of ArrayNode, which corresponds to a sequence of terms, and AlternateNode, which corresponds to | in regex. (There's a bunch of other nodes for specific regex features like PlusNode and StarNode.) So A(T|C)G becomes ArrayNode('A', AlternateNode('T', 'C'), 'G'). This AST provides an easier-to-use representation than a regex string, and it'll be very useful for further analysis (see Constraint Generation below).

It wasn't completely without hitches. Python doesn't do tail recursion elimination (Guido doesn't like it), and we quickly hit recursion limits. Our recursive code looked something like parse(ATG) => Sequence(A, parse(TG)) => Sequence(A, Sequence(T, parse(G))). The code is dead simple this way, but you quickly run out of stack space for a longer regex. We tokenized consecutive characters, so ATG is treated as a single token instead of three tokens, which addresses the most common use case of consecutive, normal bases. For genuinely long regexes (e.g. a long amino acid chain: (AAA|TTT){1000 more}), we simply bumped the recursion limit. We have some ideas for improving the parser, but more on that later.

Constraint Generation

Now that we have a parsed regex, we need to extract as much information as we can for the search. The idea here is to map the regex to a filter consisting of boolean expression - for a regex, what are filters we can generate to help narrow down potential matches? We have a function generate_constraints(...) that takes in a regex and returns a tree of AND and OR nodes (with strings at the leaves). We'll use this as an example:


After parsing, we have

=> ArrayNode(ATTAAA, StarNode(DotNode), AlternateNode(TTT, GGG))

Our generate_constraints function takes in regex in the form of a parsed AST and recursively calls itself: generate(ArrayNode(X, Y)) becomes AND(generate(X), generate(Y)) and generate(AlternateNode(X, Y)) becomes OR(generate(X), generate(Y)). Any nodes that are optional - StarNode, for example - are replaced with an ANY() node. Any string literals are left alone, as strings.
Our example becomes

=> ArrayNode(ATTAAA, ANY(), OR(TTT, GGG))

We then simplify the tree: ANY() inside of AND() is removed, and OR() that contains ANY() is replaced with ANY(). Any single-term AND or OR clause is unwrapped - AND(X) becomes X.


Finally, we flatten constraints a bit - you get a lot of AND(X, AND(Y, Z)) due to how we parsed things, which we replace with AND(X, Y, Z) since AND is associative. The same goes for OR. This is actually really important - we found that Elasticsearch doesn't do much algebraic simplification, so if you search for AND(AND(X, Y)) a query could take 4 seconds while AND(X, Y) would take 10 milliseconds.

Now that we have a set of boolean constraints, we can pass it into Elasticsearch.

Elasticsearch Indexing

If you're not familiar, Elasticsearch is a distributed document store that is really good at creating inverted indices. Given an input source, you configure how it extracts "terms" out of the source, and Elasticsearch then provides a rich API to query for which documents contain your terms.

In our case, we used Elasticsearch to index DNA bases. Elasticsearch's out-of-the-box text indexing won't work here: Elasticsearch splits by whitespace and normalizes to lowercase, which works great for e.g. English documents. However, DNA is one long string - a common approach, especially for substring search, is to index by trigrams instead. (An n-gram is just a fancy word for a substring of length n. Given abcd, the 3-grams, or trigrams, are abc and bcd.) We pass it a DNA string and have it extract all trigrams for indexing.

We then search using trigrams - if our query is ATTG, we'll search for any DNA that has ATT and TTG. This can have false positive matches (such as a sequence TTGGATT), but the goal here is to quickly eliminate sequences that definitely don't match. For the remaining matches we can run the actual query to verify - as long as our false positive rate isn't too high, this works well.

It gets one level deeper since we're using boolean filters from above. Elasticsearch uses a JSON-formatted DSL, so our constraint from above looks something like:

  "and": [
    {"query": {"match": "ATTAAA"}},
      "or": [
        {"query": {"match": "TTT"}},
        {"query": {"match": "GGG"}}

Elasticsearch will (basically) perform each query independently - each query gets the same trigram treatment, so ATTAAA actually becomes ATT, TTA, TAA, AAA and a sequence that has all the trigrams will match that clause. Those results are intersected (due to the and) with the sequences that match the or clause. After running it through the set of boolean clauses, the sequences that match have been narrowed down and are good candidates for the original regex. In our experience, this is super fast (a few milliseconds), even with millions of sequences.

Inspired by a thread on the Elasticsearch mailing list, we use a groovy script to run through each of the final matching sequences with the original regex. Originally we tried using the built-in Regexp Filter, but Elasticsearch would attempt to run it first, before any filtering from the boolean query, resulting in a slow scan over the entire database. Elasticsearch will run custom scripts last, which alleviates this problem.

The groovy script is run against each matching document and filters by the regex. This is where false positives are removed, leaving us with the DNA we were looking for!

(Side note: configuring a groovy script was surprisingly annoying: Elasticsearch allows both dynamic scripts (where the query itself has the script to run as a string) and stored scripts (where the script is stored on the node). Running dynamic scripts on Elasticsearch looks like a security vulnerability (they're not sandboxed and are discouraged), so we had to use Chef to upload a .groovy file to each node in the cluster.)

Gotcha: Indexing

We found indexing to be the most problematic step. We already had an indexing infrastructure in place to re-index sequences on change, so that wasn't problematic. However, out of the box Elasticsearch choked on longer sequences, and we haven't had too much time to tune the cluster.

We originally indexed not only trigrams, but 4-grams and 5-grams as well - the thinking there being that the larger n-grams would have less false positives. We measured that the indexing time increased linearly with the number of n-gram types, so indexing 3-grams took half as long as indexing both 3-grams and 4-grams. (This is not so surprising if you assume it's just reading it over twice, but ideally it could do both in one pass.) Luckily, only 3-grams still provides us with fast enough searching (not too many false positives), so we stuck with that.

It still takes 30 seconds for a 200,000 basepair sequence. This is on the higher end of what our indexing processes can bear, so we currently only index sequences shorter than 200kb (200 kilobases). A ton of use cases are addressed with indexing sequences under 200kb, but we'd obviously like to support larger sequences as well. More on that later.

Gotcha: Circular Sequences

In bacteria, DNA can come in a circular form (called a plasmid) - this means that for circular DNA, our regexes can wrap around the end of a string: T.*A matches AAATTT if AAATTT is circular! This is a twist on the age-old interview question of checking if one string is a rotation of the other - we store the DNA, if circular, as a doubled string. So AAATTT is stored as AAATTTAAATTT - we run the same search as before against it, but in our groovy script we also check that the length of the match is less than the length of the original DNA string. (Otherwise, something like AAA.*AAA would incorrectly match AAATTTAAATTT.)

We need one last change - a regex like AAA.*T will match greedily and now breaks: it matches AAATTTAAATTT and is filtered out by the length filter we just added. Instead, during our parsing step we transform the AST to a non-greedy version - AAA.*T becomes AAA.*?T. Now, the match is just AAAT, and we're all set.

Next Steps

As mentioned earlier, our recursive descent parser works but isn't a great fit for Python due to lack of tail recursion optimization. Moving forward, we'd probably use a parser generator like ANTLR to generate an LL parser, which won't have the same recursion issues. ANTLR is heavier and requires some more tooling to work, but it's a lot more efficient, pretty standard in the parsing community, and less likely to be buggy than writing one ourselves. Our recursive descent parser works in the meanwhile, though!

Matching longer sequences is going to take more infrastructure work. Our current thinking is that we can cover a large percentage of remaining use cases if we can support indexing long sequences but still cap the length of the allowed match (to about 100kb). We could index sequences in 200kb chunks, each one overlapping by 100kb, and parallelize across multiple Elasticsearch shards to speed up indexing as well as search. A search would look across all chunks (any matches than 100kb would guarantee to be found) and then de-duplicate the resulting sequences.

Comments or ideas? I'd love to hear them: josh at benchling.com.

Finally, as always, we're hiring engineers who are excited to make research faster by solving these sorts of problems.

Discuss on Hacker News

<![CDATA[A Deep Dive into React Perf Debugging]]>This is the second half of a 2-part series on performance engineering in React. In Part 1, we took a look at how to use React's Perf tools, common React rendering bottlenecks, and a few tips for debugging. Take a look if you haven't already!

In Part 2, we'll now

http://benchling.engineering/deep-dive-react-perf-debugging/40c61515-bf20-4881-b744-264ec055e5c5Wed, 24 Feb 2016 21:23:08 GMTThis is the second half of a 2-part series on performance engineering in React. In Part 1, we took a look at how to use React's Perf tools, common React rendering bottlenecks, and a few tips for debugging. Take a look if you haven't already!

In Part 2, we'll now dive deeper into debugging workflows - given all these ideas, how do they manifest in practice? We'll walk through a few real-life-inspired examples and use Chrome devtools to diagnose and fix performance. (And if you end up with any suggestions or additions, let me know!)

We'll refer to the sample code below - you can see that it renders a simple todo list in React. In the JS fiddle snippets that follow, you can click "Result" for an interactive version, complete with performance repros. We'll post updated JS fiddles as we go.

Case Study #1: TodoList

Let's start from the TodoList above. Try typing fairly quickly into the non-optimized code, and you'll notice how slow it is.

Let's start the Timeline profiler in Chrome's dev tools, which gives a detailed profile of what the browser is doing: handling user events, running JS, rendering and painting. Type in a single character in the input, and then stop the timeline profiler. The slowness won’t be as apparent, since we typed just a single character, but it is the fastest way to produce the minimum amount of information needed for profiling.

Note the long bar for Event (textInput) and that Aggregated Time shows 121.10ms in Scripting (Children). The timeline profiler reveals that the slowness is a scripting issue and not a style / recalculate issue.

So we dive into scripting. Switch to the Profiles tab - while Timeline gives you an overview across the browser (and supports JS Profile), Profiles lets us do a deeper dive into JS-land and has some different visualization tools. Recording another profile indicates that the slowness isn't from our application code:

Looking at the profile by Heavy (Bottom Up), descending sort by Total indicates that most of the time is spent in React's batchedUpdates call, which is a pretty clear hint that it's in React-land. In contrast, Self measures time spent in the function, excluding time in child functions - sort by Self to see if there are any specific expensive functions. It looks like there's no obvious bottleneck in a user-land function, so let's try React's Perf tool instead.

To generate a measurement profile for our slow action, in the console, we call React.addons.Perf.start(), perform the slow action by typing a character, and then finish by calling React.addons.Perf.stop(). We can then see time spent on unnecessary renders with React.addons.Perf.printWasted():

The first item shows TodoItem was rendered by Todos; however Perf.printWasted() found that there we could have saved 100ms if we avoided re-building the render tree. This looks like a prime candidate for optimization.

To diagnose why TodoItem wastes so much time, we've created a custom mixin, cleverly named WhyDidYouUpdateMixin. It hooks into a component and logs information about which updates happened and why. Here's the code; feel free to adapt it to your own needs

Once we add this mixin to TodoItem, we can see what's going on:

Aha! We see that tags is the same before and after - the mixin says it's avoidable if two objects are deeply equal but not strictly equal. On the other hand, it has a harder time figuring out if two functions are equal, since Function.bind produces a new function even with the same bound args. These are useful clues, though - we look back at how we're passing in tags and deleteItem, and it looks like we're passing in new values each time we construct a TodoItem.

If we instead pass in the unbound function to the TodoItem, and we store the tags as a constant, we avoid the issue:

WhyDidYouUpdateMixin now shows that the prevProps is shallow equal to the new props. We can use PureRenderMixin, which will skip updates if props (and state) are shallow equal.

When we run the profiler again, we see that it only takes about 35ms (4x faster than before):

This is better, but still not ideal. Typing in an input should not result in this much overhead. It turns out we're still doing O(number of items in list) work. We've merely reduced the constant, but we still need to do a shallow compare for each item.

At this point, you might decide that 1000 items in a todo list is already an extreme case, and 30ms is reasonable for your application. If you're expecting to support a few thousand children, though, this isn't quite yet the ideal 60fps (16ms per frame - any slower and you get noticeable lag).

Breaking the component into multiple components is a reasonable next step (arguably a valid first step as well). We observe that the Todos component really consists of two disjoint subcomponents, an AddTaskForm component containing the input and a button, and a TodoItems component containing the list of items.

Each of these refactors provides substantial performance gains:

  • If we create a TodoItems component that uses PureRenderMixin, then we avoid doing the O(n) work since it would skip re-rendering each item, as prevProps.items === this.props.items.
  • If we create an AddTaskForm component, the state of what text has been entered can now live there. When that text changes, the Todos component no longer re-renders (avoiding the O(n) render).

Combined, we end up at 10ms per keypress!

Case Study #2:

Scenario: We want to conditionally render a warning if the user has too many tasks (> 3000), and we also want to style the todo items so that every other item has a background color.


  • We have a similar todo list example, with TodoItems implemented - for this example we'll store the text of the input on the top-level component state.
  • We create a TaskWarning component that chooses to render the message based on the number of tasks. To encapsulate the logic within the component, we have it return null if it shouldn't render.
  • We add some CSS to style div:nth-child(even) with a gray background.

Observation: Type quickly in the input, the page lags quite a bit (with no more than 3000 tasks). If we first add 1 more item to our todo list (> 3000 tasks), the lag goes away when button mashing. Surprisingly, adding more tasks fixes the problem!

Debugging: The timeline profile shows something very interesting:

For some reason, typing a single character causes a large Recalculate style that takes 30ms (which is why if we type faster than 30ms/character, we observe jank).

Take a look at the First invalidated section towards the bottom of the above image. This is indicating that Danger.dangerouslyReplaceNodeWithMarkup is causing layout invalidation, which then leads to style recalculation. Here's react-with-addons.js:2301:

oldChild.parentNode.replaceChild(newChild, oldChild);

For some reason, React is replacing a DOM node with an entirely new DOM node! Recall that DOM manipulations can be expensive. Using Perf.printDOM(), we can see what DOM manipulations React is performing:

The update attributes reflect typing abc into the input while TaskWarning is not visible. However, the replace indicates that React is deciding to touch the DOM for the TaskWarning component, even though it seems like it should have an identical virtual DOM.

As it turns out, React (<= v0.13) uses a noscript tag to render "no component" but incorrectly treats two noscript tags as not equal: the noscript ends up needlessly replaced by another noscript. Furthermore, recall we styled every other div with a gray background. Because of the CSS, an individual render of any of the 3000 item nodes is dependent on the sibling nodes before it. Each time the noscript tag is replaced, the subsequent DOM nodes have their styles recalculated.

To fix this issue, we can either:

  • Have TaskWarning return an empty div
  • Move the TaskWarning component inside a div as it won't be able to affect CSS selectors of nodes that follow it.
  • Upgrade React :-)

But that's beside the point. The key takeaway here is that we were able to diagnose this ourselves, just from the timeline profiler!


I hope this was useful in showing how React performance issues can appear in dev tools - a combination of Timeline, Profiles, and React's Perf tools goes a long way.

Todo lists with thousands of items and arbitrary coloring may seem contrived, but there were very similar problems we encountered when rendering large documents and spreadsheets, when building an electronic lab notebook. And yes, we're still growing our team - if building complex React apps excites you, get in touch.

Any suggestions or comments? Let me know! You can reach me via saif at benchling.com.

Discuss on Hacker News

<![CDATA[Performance Engineering with React]]>This is the first half of a 2-part series on performance engineering in React. Part 2 - A Deep Dive into React Perf Debugging is now up!

This post is for those of you with a complex React application. If you're building something smaller, you might not need to focus

http://benchling.engineering/performance-engineering-with-react/886c29d9-b501-42b9-9c1c-1593b2e2896bThu, 04 Feb 2016 18:18:37 GMTThis is the first half of a 2-part series on performance engineering in React. Part 2 - A Deep Dive into React Perf Debugging is now up!

This post is for those of you with a complex React application. If you're building something smaller, you might not need to focus on performance yet. Don't prematurely optimize! Go build things!

However, if you're building DNA design tools, gel image analysis software, a rich-text editor, or full-feature spreadsheets, you're going to hit performance bottlenecks, and you're going to need to solve them. We hit our fair share of this at Benchling, so this post attempts to share some of what we've learned - it's targeted towards folks out there on the Internet and fellow Benchlings alike. (And yes, if you like these kinds of problems, we're hiring!)

In this post, I'll cover the basics of using React's Perf tools, some common issues that can lead to React rendering bottlenecks, and tips to keep in mind while debugging.

Baseline React

To recap browser performance in 3 sentences: you ideally want to render 60 frames per second in the browser, leaving you 16.7ms per frame. When your app is slow, it's often taking too long to respond to user events, taking too long to process the data or to re-render the new data. In a majority of cases, you're not doing complex data processing onClock - you're wasting time re-rendering.

By using React, you instantly get some performance gains without any extra work:

Because React handles all DOM manipulations, you largely avoid issues regarding DOM parsing and layout. Behind the scenes, React maintains a virtual DOM in JavaScript, which it can use to quickly determine the minimal changes needed to bring the document to the desired state.

Because a React component's state is stored in JavaScript, we avoid accessing the DOM. A classic performance issue is accessing the DOM at inopportune moments, which can result in issues like forced synchronous layouts (in order to read e.g. someNode.style.left, the browser is forced to render a frame). Instead of doing:

someNode.style.left = parseInt(someNode.style.left) + 10 + "px";  

we would declaratively call <SomeComponent style={{left: this.state.left}} />, and to animate the component we would simply update the state without reading from the DOM:

this.setState({left: this.state.left + 10}).  

To be clear, these optimizations are possible without React - I'm simply pointing out React tends to solve these problems ahead of time.

For simple applications, these performance optimizations that come with React are sufficient - I think of them as the minimum work needed for its declarative framework to be feasible. However, as you develop more complicated views, maintaining and comparing virtual DOMs can become an expensive operation. Fortunately, React provides some tools to detect where performance issues exist and means for avoiding them.

Performance issues caused by debugging

Watch out - there's some overhead that debugging itself can cause, leading to confusing debugging sessions that seem to go away in production.

Elements pane

The elements pane is a nice, simple way to see what is getting re-rendered - it flashes a color when an attribute changes or a DOM node is updated/inserted/replaced. However, the flashing/re-rendering of the element panel will affect performance! I often switch away from the Elements pane to the Console to get a more accurate sense of the FPS.


In development builds of React, PropType validation occurs when a component is rendered - the props that components receive are checked to help with debugging and development. Using the Chrome's JS Profiler, you may observe that your React component spends a majority of the time in the validate method:

While the development build raises warnings useful for debugging, they incur costs that don't reflect production. I sometimes switch to using React's production build to ignore this false sense of sluggishness. (To enable the production build, set NODE_ENV to production: https://facebook.github.io/react/downloads.html#npm.)

Identifying Perf issues with React.addons.Perf

Before we dive into common fixes, it's important to emphasize that you should only spend time fixing issues that you were able to measure. It's very easy to end up shooting in the dark if you're not disciplined - again, focus on building things and only invest time on fixing the key performance bottlenecks.

Identifying bottlenecks using standard debugging tools still works, but it is often difficult to interpret the data because application code can result in more time spent in React-land code (e.g. a complicated render method you wrote runs quickly, but the resulting virtual DOM computations are much more expensive). It becomes difficult to identify what application code triggered the visible bottlenecks in React-land.

Fortunately, React is bundled with some perf tools that can be used in the non-production build of React (docs). You'll find it as React.addons.Perf via react/addons in <= 0.13, and in its own react-addons-perf package in 0.14 onwards.


To use Perf, simply call Perf.start() from the console, perform the actions to record, and then call Perf.stop(). You can then call one of the methods below to view sets of useful measurements.

When I'm in Perf-Debugging Mode, I like to use a simple record button to start/stop recording easily. (The code is super simple - it's just a React component in the corner of the screen that calls through to React.addons.Perf.) It looks like the React DevTools are also looking into adding one officially: facebook/react-devtools#71. Jeff had a great tip where he binds keyboard shortcuts to start/stop, useful for when you need to profile with the mouse.


Perf.printWasted() is easily where most of the usefulness of React.addons.Perf comes from. It tells you how much time was wasted doing render tree construction and virtual DOM comparisons that did not result in a DOM modification. The components that are surfaced here are prime candidates to be improved via PureRenderMixin / other techniques that I listed above.

Perf.printInclusive() / Perf.printExclusive()

Perf.printInclusive() and Perf.printExclusive() print the times it took to render your components. I haven't found this to be too useful, as the bottleneck in rendering is often solved by not rendering (via Perf.printWasted() analysis and PureRenderMixin, more on that later) vs rendering faster. However, it can help highlight which components that perform expensive computations in lifecycle methods. I have typically found that after resolving printWasted issues it is my application code that is expensive. At that point it's easier to use the standard Chrome DevTool's JS Profiler and look directly at the most expensive function calls.


Perf.printDOM() returns all DOM operations that took place when rendering the React trees. In my experience, this is often difficult to interpret / visualize as it is a long list of entries describing exactly what happened, e.g. each attribute change and each DOM insertion, and if your app is sufficiently complicated this can represent a fairly large changeset.

After your initial component is rendered, it is expected that future re-renders should be re-using / updating existing DOM nodes and not creating new ones - after all, this is the optimization that React's virtual DOM affords us.

I've used this occasionally to discover weird, expensive browser quirks or notice large, unexpected amounts of DOM modifications.

Avoiding renders with shouldComponentUpdate

While React does wonders by maintaining a virtual DOM representation to avoid expensive DOM operations, the maintenance of that virtual DOM can also be expensive. Imagine a very large, complicated render tree. If you update the props of any node, React needs to recompute the render tree from that node all the way down to the leaves, where it finally does its virtual DOM comparisons. Fortunately, React provides a mechanism for avoiding this, namely shouldComponentUpdate - return false from this method, and the component's entire subtree does not bother re-rendering. We just need to figure out how/when to return false.

The simplest way to take advantage of this is to keep your render methods pure - your component should render only based on state and props (as opposed to reading from the DOM, cookies, or anything else.) This "pure rendering" technique is mentioned quite often, but it's to stress as it's a good habit that leads to components that are easier to reason about. You'll have external state once in a while - try to isolate that to a few components and keep the rest pure.

Once you do this, you can use PureRenderMixin in your components. From the source, the mixin just calls through to shallowCompare. (If you're using ES6 classes, you'll want to use shallowCompare directly as well.)

var ReactComponentWithPureRenderMixin = {  
  shouldComponentUpdate: function(nextProps, nextState) {
    return shallowCompare(this, nextProps, nextState);

If no changes in props/state are detected, it doesn't re-render. For a component to still behave correctly, a component should ensure that:

  • render() is strictly dependent on props and state, i.e. it shouldn't be reading values from some global state.
  • props and state should never be mutated - any changes should result in a new variable, since shallowCompare only checks for strict equality between top-level props. react-addons-update will help you perform immutable updates, and we use Object.assign / _.extend for simple cases as well. ImmutableJS is more significant switch, but once you use immutable data structures you can use PureRenderMixin easily. Beware of the urge to do things like this.state.myItem.stars++ - it's easy to forget that you're mutating the state directly, especially if things happen to work due to some other state changing as well.

If you adhere to pure components while developing, when you notice bottlenecks it becomes much easier to toss in a PureRenderMixin.

A small gotcha

You may get a false sense of performance increase if you use PureRenderMixin - it will avoid propType validation in child components, which seems faster in development but will be skipped anyways in production builds.

A much larger gotcha!

Even if you adhere to this stricter policy, you may not immediately reap the benefits of PureRenderMixin. As described earlier, React does a shallow-equal comparison, not a deep comparison, to decide if it needs to re-render. There are surprisingly many easy ways to accidentally create deep-equal props that are not shallow-equal (more on that later).

One quick way to address this is to use deep comparison, e.g. _.isEqual:

shouldComponentUpdate(nextProps, nextState) {
  return !_.isEqual(this.props, nextProps) ||
    !_.isEqual(this.state, nextState);

If you mostly re-use props, _.isEqual first does a shallow comparison, so comparing the props that haven't changed is cheap and performance is fine. In practice, where _.isEqual suffices, we haven't found the deep compare to be a performance issue.

You could also write a custom shouldComponentUpdate method that is specifically tailored to your component, but I'd only do this on simpler components. If the custom method isn't properly maintained, you will run into issues with your component not updating when it really should.

Optimizing for shallow-equal props

Often, it turns out using best practices to avoid creating new objects when possible helps naturally with rendering optimizations.

Function.bind() / inline (anonymous) functions

Function.bind is a convenient way to expose a contextual call to a component. Unfortunately, each call to Function.bind produces a new function:

> console.log.bind(null, 'hi') === console.log.bind(null, 'hi')
> function(){console.log(‘hi');} === function(){console.log(‘hi');}

// New function each time
render() {
  return <MyComponent onClick={() => this.setState(...)} />

No amount of prop checking will help, and your component will always re-render. (You can turn on the react/jsx-no-bind eslint rule to disallow using bind and arrow functions in jsx props.)

The simplest solution we've found is to pass the unbound function and the desired args to the child component, and have the child use an instance method, e.g:

const TodoItem = React.createClass({
  deleteItem() {

It feels odd to expose a more general method to a subcomponent, with the contract that it will pass back in the index. To get better encapsulation for more complex scenarios, we sometimes use an IntermediateBinder whose only purpose is to act as a binding context for e.g. the id argument. It takes the id as a prop, defines its own methods bound to itself, and passes those methods to the child component.

This allows us to write:

  {(boundProps) => <TodoItem deleteItem={boundProps.deleteItem} />}

(Another possibility we've explored is using a custom bind function that stores metadata on the function itself, which in combination with a more advanced check function, could detect bound functions that haven't actually changed. This didn't seem explicit enough for our tastes.)

Literal array/object construction

It's simple but often overlooked. Array literals will break PureRenderMixin:

> ['important', 'starred'] === ['important', 'starred']

If this object is never expected to change, you can move it into a module constant / component static variable:

const TAGS = ['important', 'starred'];


Defining content boundaries between a component and its subcomponent often lends itself to easy performance optimizations - well-encapsulated component interfaces lend naturally to performant updates. Refactoring out intermediate components can help improve where you can use PureRenderMixin and save updates:

  <ComplexForm props={this.props.complexFormProps} />
    <li prop={this.props.items[0]}>item A</li>
    ...1000 items...

In this case, if complexFormProps and items come from the same store, typing in the ComplexForm might lead to store updates, and each store update leads to re-rendering the entire <ul>. Virtual DOM diffing is great, but it still has to check every <li>. Instead, refactor out <ul> into its own subcomponent that takes in this.props.items, and only update if this.props.items changes:

  <CustomList items={this.props.items} />
  <ComplexForm props={this.props.complexFormProps} />

Cache expensive computations

This goes against the "single source of state" principle, but if computations on a prop are expensive you can cache them on the component. Instead of directly using doExpensiveComputation(this.prop.someProp) in the render method, we can wrap the call that caches the value if the prop is unchanged:

getCachedExpensiveComputation() {
  if (this._cachedSomeProp !== this.prop.someProp) {
    this._cachedSomeProp = this.prop.someProp;
    this._cachedComputation = doExpensiveComputation(this.prop.someProp);
  return this._cachedComputation;

Candidates for this optimization would be best discovered using the JS Profiler.

Link State

React's Two Way Binding Helpers can be very useful for simple inversion of control, allowing a child component to communicate new state to the parent. If only used with valueLink for a React form component, it isn't so bad as the React form inputs are very simple. But if you start threading it through more components like we were doing, you may run into issues. linkState is implemented as follows:

linkState(key) {
  return new ReactLink(
    ReactStateSetters.createStateKeySetter(this, key)

Every call to linkState returns a new object, even if the the state hasn't changed! This means shallowCompare will never work. Our workaround is unfortunately simply not to use a linkState. If you instead flatten the linkState into a getter prop and a setter prop, we avoid creating a new object, e.g. nameLink={this.linkState(‘name')} could be replaced with name={this.state.name} setName={this.setName}. (We've considered writing a linkState that caches itself...)

Compiler Optimizations

Newer versions of Babel and React support inlining React elements and automatically hoisting constant React elements. We haven't played too much with this yet, unfortunately, but they will help with reducing calls to React.createElement and in speeding up DOM reconciliation, respectively.

Wrapping Up

We went through a lot just now (you should've seen the original list!), but the key point to take away is that you should 1) get comfortable with profiling and 2) shouldComponentUpdate will get you a long way. We hope this has been useful!

Any suggestions, comments, or things we've missed? Let us know - saif at benchling.com.

Stay tuned for part 2, where we'll discuss our React debugging workflows, dive into real examples of non-performant code, and subsequently fix them.

Discuss on Hacker News

Update: Part 2 is out! Check it out - A Deep Dive into React Perf Debugging.

As always, we're looking for product-loving folks to join the team. :)

Thanks to Jeff Chan, Victoria Sun, Harry Yu, and Mark Zhang for reading drafts of this.

<![CDATA[Powering CRISPR with AWS Lambda]]>In our last post, we talked about the algorithm we use to do sub-second CRISPR searches on billions of DNA bases. Our initial search infrastructure was too expensive and unscalable, so one of our summer interns, Harini, focused on rebuilding this infrastructure using AWS Lambda. In this post, we'll focus

http://benchling.engineering/crispr-aws-lambda/34027238-5049-49d0-9d3f-8e5766922a0cFri, 25 Sep 2015 09:44:12 GMTIn our last post, we talked about the algorithm we use to do sub-second CRISPR searches on billions of DNA bases. Our initial search infrastructure was too expensive and unscalable, so one of our summer interns, Harini, focused on rebuilding this infrastructure using AWS Lambda. In this post, we'll focus on how we use AWS Lambda to parallelize our CRISPR searches.

The Problem

CRISPR Cas9 is a revolutionary genome editing technique that allows scientists to modify parts of a genome with extreme precision. It's already being used by thousands of scientists worldwide to build disease models, cure genetic disorders, and more1.

A CRISPR search involves finding a 20 base string of DNA (a "guide") that occurs in the target region of the genome and nowhere else. Given a particular guide, we need to find all locations in the genome that are "close enough" to the guide.

We want to be able to do this search across dozens of genomes in a way that's fast, cheap, and dynamic. In our case, dynamic means it should be easy to add new genomes. Since our algorithm doesn't require any indexing or preprocessing of the genome, we should also be able to support genomes that users upload.

When analyzing our infrastructure, we'll keep the following in mind:

  • We currently run a few hundred thousand CRISPR searches a month. These searches happen in bursts, and are not uniformly distributed across the day.
  • The speed should not depend on how many concurrent searches are running.
  • We currently have 30 genomes, each 3GB in size; we should be able to support an arbitrary number of genomes.
  • We have an algorithm (implemented in C++) that takes in a string of DNA and list of queries, which runs at about 1 second per 100MB.
  • To reduce complexity, we'd like this to run easily in our current AWS architecture.

From these requirements, we have an idea of what our infrastructure needs:

  • Dynamic server allocation - we want to serve requests quickly even at peak times, and it would cost too much to allocate for the maximum load.
  • Shared file system or separate servers per genome - we have 90GB (and growing) of genome data, so storing this on the disk of each (dynamically allocated) server is not feasible.

Original Infrastructure

Our old CRISPR infrastructure used several servers to process CRISPR search tasks. Each server stored all of the genomes on disk. Using our existing Python Celery infrastructure, a single CRISPR search was split into several subtasks and parallelized across the servers. Each subtask read the corresponding genome region from disk, performed the search, and returned the results. The results were combined by a master task and returned to the user.

There were several issues with this approach:

  • All genomes were stored on each of the servers, making maintenance difficult.
  • We couldn't spin up new servers fast enough if demand was high.
  • We were paying for the servers even when they were not being used.

These are all common problems with scaling, and Amazon provides several possible solutions. Auto Scaling is one such service, but is too slow for our spiky requests. Each CRISPR search takes around 1 second, so waiting 1 minute for an EC2 instance to start is too much. Instead, we went with Amazon's new serverless platform, AWS Lambda. In the following sections, we'll discuss the infrastructure we built on AWS Lambda, and the various complications we ran into.

AWS Lambda


AWS Lambda provides all of the infrastructure around dynamic server allocation. To get started, you upload a ZIP file with the code you need to run, and specify the entry point for the request. For any request, Lambda will run the code, creating new servers if necessary. You pay for the amount of time used and memory allocated for each request. We use Lambda to offload long-running and parallelizable work from our web servers.

New Infrastructure

Our goal is to split up a CRISPR search across several Lambda "tasks" to reduce costs and increase scalability. The CRISPR search problem is easily parallelized by splitting the genome into smaller regions that can be processed separately. The general approach is fairly straightforward:

  1. A web server receives a request to do a CRISPR search on a specific genome.
  2. The web server splits up the genome into smaller regions.
  3. The web server invokes the Lambda function for each region.
  4. Once all of the Lambda tasks complete, the results are combined and returned to the user.

This approach completely solves our dynamic scaling issue. We don't need to maintain several servers to perform searches, as AWS will handle all of the allocation for us. That was easy!

But, as always, there were...


AWS Lambda imposes several limitations on the size and types of requests it can process. Here's a quick list of the main issues we ran into, and how we solved them.

1. Running C++ code with AWS Lambda

Lambda only supports Javascript (Node.js) and Java, so we needed a way to run our C++ code. Porting the code over to Javascript or Java would have resulted in a significant performance hit.

Node.js provides a way to call into C++ libraries directly with addons. We used node-gyp to compile our C++ code into a binary that Node.js can call directly.

Lambda also only runs on Amazon Linux, so we needed to compile our C++ code on an Amazon Linux machine. To deploy a new AWS Lambda function, our process is:

  1. Spin up an Amazon Linux machine
  2. Compile the search code using node-gyp
  3. Upload the compiled code as a new AWS Lambda function.

Although this process is fairly cumbersome, it hasn't been a problem, as the core search functionality has remained relatively unchanged. We can automate this process if necessary.

2. Passing genome data to the Lambda function

The Lambda function needed to easily access regions of the 3GB genome. Amazon has several size limitations that made this hard:

  • Maximum size of the uploaded zipped code: 50MB
  • Maximum size of the unzipped code: 250MB
  • Maximum request size: 6MB
  • Maximum response size: 6MB

We considered several options:

  1. Pass the genome bases as an argument to the Lambda function
  2. Include the genome region in the zipped code
  3. Download the genome region from S3 on each Lambda invocation

Option 1 was too slow, and option 2 would have been hard to maintain with dozens of genomes. We went with option 3.

We stored each of the genomes as its own file in S3. We then specified the genome region to the Lambda function by passing in the genome name with start and end indices.

function CRISPRSearch(genome, startIndex, endIndex) {  
    # Download genome region from startIndex to endIndex
    # Run the C++ search algorithm
    # Return the results

The S3 API supports downloading specific regions of a file, so we only download the bases we need. The bandwidth between S3 buckets and EC2 instances in the same region is around 50MB/s2, so this was fast with small enough genome regions.

A large issue we ran into was the Lambda Resource Model.


Each Lambda function has an allocated memory size. The only way to increase CPU speed is to pay for more memory. Download speeds with EC2 are directly proportional to the CPU power, and thus the memory. This meant that even though our algorithm required very little memory, we had to pay for 1GB to get 50MB/s download speeds.


We wanted our infrastructure to be fast, cheap, and dynamic. Let's look at how our new Lambda infrastructure compares with the old server infrastructure.

1. Speed and Scalability

Our old server infrastructure was built to run CRISPR searches in around 1 second. It unfortunately did not scale well as demand fluctuated, and would have costed a lot as overall demand increased.

Our new Lambda infrastructure can theoretically run at arbitrary speeds, since we can choose the amount of parallelization. In reality, there are other problems as the number of Lambda tasks gets large, so we currently run it at around 1 second per CRISPR search. It also scales perfectly with demand, meaning that even at peak times we can perform CRISPR searches quickly.

2. Cost

Our old server infrastructure cost thousands of dollars each month just for server costs.

Using the new Lambda infrastructure, we pay for the number of Lambda invocations, the total duration of the requests, and the number of S3 requests3. This comes out to $60/month for hundreds of thousands of CRISPR searches!

3. Dynamic User Uploads

Our old server infrastructure required us to add new genomes to the disk of each server. It also made analyzing user uploaded genomes difficult, as we would have needed to copy the uploaded genome to several different servers.

The Lambda infrastructure gets around this problem easily, since the genome is only stored in one place. We can now support user-specified genomes by uploading them directly to S3. As CRISPR becomes more widespread, this will become more useful for companies and users that have custom genomes to analyze.


This was our first foray into using AWS Lambda. We were able to use the serverless architecture to parallelize long running CRISPR searches. It was significantly cheaper, and will scale well as our CRISPR usage grows. We're super excited to be building tools that scientists use every day to do their research.

We're Hiring!

And yes, Benchling is hiring engineers to join our small-but-growing team - if you're interested in building high-quality tools to power a new generation of research, see our jobs page or contact us.

  1. Hsu, Lander, Zhang 2014: http://www.cell.com/cell/pdf/S0092-8674(14)00604-7.pdf

  2. This was measured in the us-west-2 region. Bandwidth will vary across regions.

  3. Data transfer between S3 buckets and EC2 instances in the same region is free, so we only need to pay for the number of GET requests.

<![CDATA[Optimizing CRISPR: Sub-second searches on a 3 billion base genome]]>http://benchling.engineering/optimizing-crispr/5e990d4c-c275-4246-bbde-1ee395efd748Wed, 22 Jul 2015 17:23:00 GMTCRISPR is a genome-editing technique that has rapidly gained popularity in the world of synthetic biology over the last two years.

CRISPR allows scientists to edit genomes with unprecedented precision, efficiency, and flexibility. Gizmodo

It is completely transforming how medical research is done - scientists are already applying CRISPR towards curing genetic disorders, creating disease models, accelerating drug development, and more1. We believe software will play a crucial role in these advancements, and we've released the Benchling CRISPR tool to make working with CRISPR faster and more accurate.

This post focuses on how we've improved 15-minute batch jobs to sub-second searches on a 3-billion base genome. Faster analysis means faster feedback, and faster feedback means faster iteration. Without going into too much detail about the biology of CRISPR (see Andrew Gibiansky's post for a great primer), let's go over a few basics to first understand the motivation. Then we'll dive into the code.


We model DNA as a string of base molecules, where each base is either adenine (A), thymine (T), cytosine (C), or guanine (G). So the human genome is a long string of ACGT characters - that's it2.

At its core, to work with CRISPR scientists design a 20 base string of DNA called a "guide." This guide is integrated into a special protein called Cas9, which then cuts the genome at all occurrences of the guide. Cutting the genome at a position means actually breaking the DNA strands apart. This is where CRISPR gets its powers - by designing the right 20 bases, you can break apart specific parts of the genome.

For example, to break a gene:

  1. Find a 20 base guide that occurs near the beginning of the gene
  2. Use the CRISPR/Cas9 system to cut at that location
  3. The cell's natural repair mechanisms will attempt to repair the cut. Since it's an error-prone process, continual cutting by Cas9 eventually results in an incorrect repair. The gene is essentially "corrupted" around that region, and the cell is no longer able to express that gene as intended.

However, this is also the hard part of CRISPR - given a 20 base guide, it's possible for it to occur elsewhere in the genome. You'd like to pick your guides to target a certain gene but avoid breaking other things - while CRISPR makes this much easier, it still has to be carefully done. Additionally, a perfect match isn't necessary for CRISPR to cut at a location - empirically, up to 4 mismatches can occur and still result in cutting.

So we want to find the optimal guide that matches near a gene of interest, without having too many false positives. This sounds like something a computer would be great at doing, and it's one of our core CRISPR features.

We help scientists find a set of guides and, for each guide, determine where else it binds in a 3 billion base genome.

Our Problem

Let's define a function match(guide, genome_region) which returns true if guide and genome_region are close enough that Cas9 would successfully cut at the genome_region. In other words, it will return true if they are 1) the same length, and 2) differ in at most 4 characters.


Given a genome G and a set of guides S (all with the same length), we'd like to find all locations in G that match at least one guide in S3.

Let's optimize our algorithm for the following constraints:

  • A genome length of 3 billion bases (the length of the human genome)
  • A guide length of 20
  • 200 guides in S

Our algorithm should be both parallelizable and cost efficient. We'd like it to run, parallelized, in under 1 second.

(For readability, I've used Python "pseudocode" for any code snippets - all CPU times shown here are based on C++ implementations.)

Naive Solution

Let's start with a simple solution as a baseline. We can iterate over the entire genome, and check each index to see if any guide is a match.

CPU Time: 3.25 hours

Ouch. Now, what are the constraints of our problem, and how can we take advantage of them?

  • We have the genome in advance. Each species has only one or two reference genomes that scientists use.
  • DNA strings are made up of only 4 different characters (A, C, G, and T).
  • The length of a guide is 20 bases.
  • The maximum edit distance (replaces only) is 4.

Let's start with the first constraint.

Solution #1: Index the Genome

Since we have the genome beforehand, why don't we just index the entire genome? We can create a hash_table that maps every 20 character string in the genome to the list of locations that it occurs.

For a 23 base genome AAAAAAAAAAAAAAAAAAAAATG, hash_table would be:


Since there are only 4 possible characters, we can represent each one with just two bits, using a simple encoding: A => 0b00, C => 0b01, G => 0b10, and T => 0b11. For example, TGGA would be represented as 0b11101000. This means that a 20 base string can fit into one 64-bit int. That changes hash_table to:

0b0000000000000000000000000000000000000000: [0, 1]  
0b0000000000000000000000000000000000000011: [2]  
0b0000000000000000000000000000000000001110: [3]  

Given this map, how do we find all of the matches for a given guide? We have to allow an edit distance of up to 4, so we can't just do a single lookup for the guide. For a given guide, there are 424,996 strings4 within edit distance 4, so we could just search for each one...

CPU Time: 400 seconds

Unfortunately, the pseudocode hides a very important detail - each lookup in hash_table is actually a disk seek. The table is too large to store in memory (around 36GB). Doing 200 * 424996 disk reads takes quite a long time, even on an SSD.

Solution #2: Preprocess everything

If the problem is the number of disk reads, maybe we should try preprocessing some more information from the genome.

One idea would be to store the matches for every possible guide - essentially, a full map from all possible inputs to their corresponding answer. That way, given a guide, we could do a single query into the hash table to find all of its matches.

CPU Time: Unknown

Assuming the genome is random, there should be about 1000 (32-bit int) matches for each guide. This would result in a storage requirement of at least 420 * 1000 * 4 bytes/match, or more than 4 petabytes!

Even after playing around with this for a while, it was hard to come up with a solution that was both fast and used a reasonable amount of storage.

Let's take a step back, and look at the problem as a whole. We are trying to search for 200 guides over a 3 billion base genome. We can roughly divide the solution space into two groups:

  • Process the genome, and do a separate search for each guide
  • Process the guides, and do a separate search for each position in the genome

My intuition tells me that the second scheme should be much worse. We can't process the guides in advance (whereas we can process the genome). We will have to do at least 3 billion "operations," 1 per position in the genome (versus ~200 "operations," 1 per guide, by processing the genome). But let's give it a shot anyway.

Suppose we constrained the problem - let's search for an exact match for a single guide. How long does it actually take to naively iterate over the entire genome?

CPU Time (single guide): 15 seconds

Now, let's take advantage of the limits of DNA. Since we can store an entire 20 base string in a single 64-bit integer, we should be able to check equality in a single CPU instruction. As we go through the genome, we can keep a rolling hash of the last 20 bases we've seen. To add the next character, we can use simple bitwise operations to remove the highest 2 bits, and append on 2 bits to the end. Our "hash" is a bit of a misnomer - it's just the encoded bases.

Now, to check for an exact match, we can just compare our two 64-bit ints, the guide and the rolling hash.

CPU Time (single guide): 2.61 seconds

That's fast! It only takes 3 seconds to scan 3 billion bases and check equality at every index. Of course, we're not checking for exact matches. We're looking for the 424,996 guides that are within an edit distance of 4 from each of our 200 guides. Let's alter the loop to look for all of the guides simultaneously.

Solution #3: Hash set for guides

Let's create a hash set possible_matches with an entry for every possible match, for each guide. Remembering that we're using 200 guides, there will be 424,996 * 200 = 84,999,2005 such matches in possible_matches. Then, we can simply iterate over the genome, and check at every index if a corresponding hash set entry exists.

CPU time: 773.8 seconds
Creating hash set: 41.8 seconds
Searching over genome: 732 seconds

That's weird. We just found that iterating over the entire genome to check equality only took 2.61 seconds. So why does this method take so much longer?

If we examine our loop, we are checking rolling_hash in possible_matches in each iteration. Checking existence in a hash set requires computing the hash of rolling_hash (pardon the naming). This means we need to do 3 billion hash operations as we iterate over the genome. We need to decrease the time taken at each position of the genome.

Solution #4: Two boolean arrays

Since evaluating the hash function was the bottleneck, let's replace our hash set with a boolean array. array[key] = True if and only if key in hash_set. This means a set of {1, 2, 5} would be an array of [False, True, True, False, False, True]. We can avoid hashing, at the expense of additional memory usage.

Remember that we represent each guide as a number, using the two-bit encoding. In this case, we need a slot in the array for every possible guide. For length 20 guides, we need about 1 terabyte of memory (420 1-byte booleans)6.

Instead, let's try using two boolean arrays (this is actually the key insight to solving our problem!)7. The first array first_half_matches will correspond to matches for the first 10 bases of any guide, and second_half_matches will correspond to matches for the last 10 bases. Our storage needs drop from 1TB (420 1-byte booleans) to 2MB (2 * 410 1-byte booleans).

Remember we'd like to decide if rolling_hash is within edit distance 4 of any guide. If rolling_hash's first half is not in first_half_matches, then its first half is already more than edit distance 4 from any guide and we can abort. The second half follows similarly.

If both halves of rolling_hash are present in their respective arrays, we still need to verify that rolling_hash matches overall. For example, AAAAAAAAAAAAAAAAAAAA matches AAAAAAACCC in the first half and AAAAAAAGGG in the second half but does not match AAAAAAACCCAAAAAAAGGG (edit distance of 6). For now, we'll simply check the rolling_hash against all guides to verify if it's a match.

CPU time: 25.43 minutes
Creating arrays: 0.149 seconds
Searching over genome: 25.43 minutes

We're almost there! Preprocessing the queries is now extremely fast, and uses only 2MB of memory. But searching over the genome takes more than 25 minutes. Iterating over the genome to check equality only took 2.6 seconds, so what's taking so long now?

It turns out that our first if statement (checking first_half_matches and second_half_matches) has too many false positives. On some test data, 2,898,421,200 indices pass the first if check, but only 222,000 are real matches. That means that 96.6% of the 3 billion locations pass our first filter - not much of a filter, is it? The time taken to perform the expensive check takes up more than 99.999% of our CPU time.

Let's alter our approach slightly to fix this.

Final Solution: Two Integer Arrays

Instead of using a boolean array to denote whether a 10 character string matches the first half or second half of the guide, let's also keep track of the minimum edit distance between the string and any guide. Now, we can change our filter to be more aggressive - if the summed min edit distances of the two halves exceeds 4, we can safely rule out a match.

Our false positive count drops from 2,898,421,200 to 34,706,100 (-98.8%), and our times reflect that:

CPU time: 28.5 seconds
Creating arrays: 0.16 seconds
Searching genome: 28.4 seconds

We now have a serial solution that takes only 29 seconds to complete, down from 3.25 hours!


After even more optimizations, we were able to bring the serial version down to 8 seconds. This algorithm is also "embarrassingly parallelizable" - we chop the genome into 10 or so chunks, run the search over each chunk, and merge the results together. Empirically, we were able to easily achieve that final sub-second runtime.

Here's a quick overview of the key optimizations we made:

  • Splitting an exponential space (all possible guides) into two halves
  • Replacing hash tables with boolean arrays for faster lookup
  • Operating on 64-bit ints instead of strings
  • Using a rolling hash to remember the last 20 bases of the genome
  • Parallelizing across machines

Our solution has one unexpected benefit as well - we don't need to process the genome beforehand. While not a huge priority, this gives us some flexibility in how we build our infrastructure and makes it easy to support more genomes.


So there you have it - an analysis that used to take 10-15 minutes now takes a second on Benchling. We're pretty excited to be able to contribute these kinds of improvements to the scientific community, for free.

As far as infrastructure goes, we currently use a bunch of cores on AWS c4.2xlarge boxes for the computations, but one of our interns is in the process of moving onto AWS Lambda for a super scalable (and cheap) setup. More on that in a future post!

We're Hiring!



  2. There are technically still some regions of the human genome which we are unsure about, and these regions are represented with a different set of letters. We ignore these regions when searching for matches.

  3. Technically, we also need to deal with PAM sites. This is just a minor addition to the final algorithm, though, so I've left it out of this blog post.

  4. The number of possible guides is computed as: (20 choose 4) * 34 + (20 choose 3) * 33 + (20 choose 2) * 32 + (20 choose 1) * 31 + (20 choose 0) * 30 = 424996

  5. We could represent this as a bitmap instead of a boolean array. This would cut our memory needs in 8, but we would still need 128GB. In addition, we would need at least 4 operations to check if a particular value exists - one to get the value in the array, two to extract the correct bit, and one to check if it was 0 or 1.

  6. This number actually depends on the guides themselves. Some values may overlap, so 84,999,200 is actually an upper bound on the size of possible_matches. In practice, the actual size will be very close to 84,999,200.

  7. I played around with the number and size of the arrays, and two arrays of 10 seems to be optimal for a 20 base pair guide. The tradeoffs are 1) amount of time to preprocess, 2) memory requirements, 3) number of CPU instructions in the filter, and 4) number of false positives. One interesting note is: even for a size 30 guide, we can ignore the last ten bases and use two arrays of size 10, since we have to check for false positives regardless.