Exploring random forest hyper-parameters using a League of Legends dataset

Over the last few blog posts, I have used random forests to investigate data from the game League of Legends. In this last post, I will explore model optimization. Specifically I will look at how hyper-parameters like forest size, and node size can influence classification accuracy, show that dimensionality reduction doesn't help random forests, and compare random forest performance to Naive Bayes. For details of this post, please look at this Jupyter notebook.

Background

As a refresher, League of Legends is an online multiplayer game where teams try to destroy each others' base. The aim of this project is to predict the eventual winner of a game at early timepoints, using game conditions; for example, I would like to be able to say the Blue Team with a 10,000 gold lead has an X% chance of winning, and determine X by machine learning. I gathered data from Riot Games's API, and created pandas dataframes containing information from 30,000 games at 5 minute intervals. For examples, you can see previous posts. In this post, as mentioned above, I will explore how I optimized my random forest model.

Forest size

The basic task my random forest is trying to do is predict the winner of a game of League of Legends, a classification task. Random forests work by generating a large number of decision trees, and averaging the results across trees. Thus one of the most obvious random forest parameters to optimize is the number of trees in the forest. To look at this, I ran the prediction algorithm many times, with different number of trees for each run.

As you increase the number of trees in the forest, the accuracy improves, until it plateaus around 25 trees. I am using around a dozen features in my dataset, which makes me hypothesize that your forest should have a number of trees equal to twice the number of features. It would be interesting to look at datasets with more features to see how they perform, and see how general this is.

Node size

The software package I am using for building random forests, scikit-learn, by default creates decision trees where each leaf is pure (that is, each group contains only wins or losses, and no mixture of wins and losses). This can lead to overfitting of individual trees, and by extension the random forest. However, you can set a parameter in the random forest to increase the minimum leaf size, or the minimum size of nodes for splitting. To see how this influenced prediction accuracy, I ran the prediction algorithm over a range of minimum node sizes.

The model accuracy goes up as the minimum sample size increases, plateauing around a 200 sample minimum. The model then maintains its accuracy for larger minimums, before decreasing once the minimum size is 5,000, or approximately 15% of the data. Once the minimum is that large, the trees are restricted in their depth, and cannot make enough splits to separate the data.

Rather than making graphs like this for every parameter, I ran a grid search over tree depth, leaf size, and node size. The best parameters for my dataset with 30,000 games, and a dozen features was a maximum depth of 10 levels, minimum leaf size of 10, and minimum split size of 1000. Optimizing these parameters yielded a 2-5% increase in accuracy over the default parameters (over a range of timepoints).

Dimensionality reduction

There is a lot of collinearity in my dataset, which can be a problem for certain machine learning algorithms; theoretically, random forests are resilient to collinearity. To see this for myself, I performed dimensionality reduction, specifically PCA, on my features, and then ran the random forest algorithm again. Here I am presenting the data slightly differently, in terms of prediction accuracy at specific times within the game.

The PCA model is actually a little worse! In the past I have used PCA on data before performing regression, and which improved regression performance by 5%. It's good to know that random forests work as advertised, and don't require dimentionality reduction beforehand.

Naïve Bayes

Part of the reason I did this project was to gain experience with random forests, and compare their performance with other algorithms. Naïve Bayes is a decent, fast default algorithm which is often used for benchmarking. I hadn't used it previously, as there is a little bit of manipulation necessary to combine categorical and continuous features in Naïve Bayes (like 5 lines of code!). As a final experiment, I wanted to see how Naïve Bayes compared to Random Forests.

Naïve Bayes performs almost as well! And runs a whole lot faster! What is probably happening here is that the dataset is just not deep enough for the strength of random forests to shine.

Influence of region, skill and patch on predictability in League of Legends

Last month, I used random forests to predict the winner of League of Legends games. Since then I have downloaded datasets from different regions, different ELOs, and on different patches, and checked how these factors influence predictability. This post will summarize the results, but for details of the analyses, please see this notebook.

Differences between regions

Korea is famed for its great players, and an "open mid" philosophy (in a normal game of League of Legends, you cannot surrender until 20 minutes into the game; in Korea, players accelerate this process by letting the other team win the game quickly through an "open mid"). Are we able to see those attributes in our model? As a first way to look at this, we can plot the length of games in different regions.

The EUW and NA have similar distributions of game lengths, with around 8% of games being surrendered at 20 minutes. The Korean games, on the other hand, have a few percent more games ending before 20 minutes, due to open mids, and a 14% of games ending right after 20 minutes. Since so many games are ending earlier, can we see this in the model's predictions?

Indeed we can, with games being a few percentage points more predictable for games that last 20 minutes or less. However, once the games last 25 minutes, all the regions are similarly predictable. This means that the Koreans are surrendering winnable games.

Skill and predictability

We can run the same type of analysis for skill. In the notebook I compare Bronze V games to challenger solo queue, and show that low skill games are slower, and less predictable than high skill games. Here, I will compare high skill solo queue games to master tier team ranked. The team ranked dataset is quite a bit smaller than the solo queue data (1,000 games vs 30,000), but in my tuning, I've found there is not a big difference in accuracy when the dataset is over 1,000 games.

We can start again by plotting game length.

Team ranked games are about two minutes faster than solo queue games. We can then run the prediction algorithm again on the team ranked games.

Team ranked games are much more predictable than solo queue games! (This is in fact the largest difference in predictability I have seen between different datasets). What is likely happening here is that coordinated teams are more able to capitalize on their advantages, and press them around the map. With fewer mistakes to allow comebacks, the games are faster, and more predictable. It would be interesting to further look at professional games to see how they compare.

Patches and predictability

In recent patches, Riot has made significant changes to the game, which makes it feel like the games are faster, and comebacks more difficult. I scraped challenger and master ranked games from the most recent patch, and compared them to ranked games from Season 5. We can again start with game length.

If you thought games were faster, your intuition was right, as games are around three minutes faster in preseason 2016. What about predictability?

As you'd expect given the previous relationships between speed and predictability, preseason 2016 is a few percentage easier to predict. In a previous blog post, I investigated the importance of different features on predictability, and found that dragons were not important in helping teams win. Putting on my amateur designer hat, this means that dragon is even less important, since games are less likely to last long enough for the accumulated advantages to matter.

In my next post, I'm going to go more in depth on the modeling, investigating random forest hyper-parameters, and maybe trying a few different models.

Visualizing champion difficulty in Jupyter

I've been playing around some more with League of Legends analysis and Jupyter. This time I scraped champion information from League of Graphs, and calculated which champions were the easiest / hardest to play, and which champions are best in the early / late game. Unfortunately, Blogger does not interact with Jupyter, so I have posted the notebook on nbviewer, if you are interested in taking a look!

Playing in random forests in League of Legends

I've wanted to learn more about machine learning, specifically python's scikit-learn module. I'm an avid League of Legends player (summoner names lemmingo and Umiy), and Riot Games provides a thorough API for querying game data, so I decided to explore machine learning using League of Legends (LoL). Specifically, I wanted to see if I could predict the eventual winner of a game long before the game finishes. In this post, I'm going to explore some features of the LoL dataset, describe what features are important in predicting the winner, show how the predictions change over time, and investigate how winnable surrendered games are. If you are interested in the details of how I did my analysis, go to the github page for this project, or this Jupyter notebook.

For those who don't know what League of Legends is, it's an online competitive game where teams of five players try to destroy each other's base. Some key things players do during the game are: 1. get gold; 2. destroy towers; 3. kill other players; and 4. kill boss monsters called dragons. Games usually last 25-40 minutes, and each game has a clear winner (no ties). Teams can surrender after twenty minutes if the the game is a rout.

Loading the data

The details of how I loaded the data are in the Jupyter notebook, but I will give a summary here. To get a list of players, I first queried the API for the list of featured games; from those feature games I extracted 50 player names.  Then I created a list of every ranked game those players were in during the year 2015, approximately 30,000 games in total. These should all be high skill-level games. I then queried the API for the first 6,000 of these games. This took a few hours because the API is limited to ~0.8 requests per second. A JSON with all of the data for 6,000 games takes up ~1GB (this is due to the JSON being text, and containing a lot of redundant information). Finally, I extracted features from the JSON. Now that I have my scripts more established, I can query games, extract their features, and save them directly to a (smaller) csv.

Exploring the features

As with any good data science project, the first thing I did was explore the features. I started by plotting a histogram of the game lengths:

As Riot says, the average game lasts between 25-35 minutes. Few games last less than 20 minutes, with a large spike of games ending right after 20 minutes, due to surrender. There are others spikes just after 30 and 40 minutes, which I don't have a clear explanation for. Perhaps large team fights break out around then.

The second thing I plotted was the average gold difference at 20 minutes.

This is about as normal a distribution as one could hope for, with a standard deviation of 6,300. Rather than plot each feature individually, we can plot how pairs of features are related to each other:

The features shown here are the gold difference between the teams, the kill difference, difference in number of towers destroyed, and the difference in dragons killed. Apologies for the ugly axes.

As you might expect, there is a large degree of collinearity between gold difference, kill difference, tower difference, and dragon difference. This would be a problem for linear regression, but I am going to use random forests, which are less sensitive to collinearity. If I were interested purely in prediction accuracy and generality, I would perform PCA on these components, but for now I'll just use them raw.

Feature selection

I was curious about which features are most predictive of the eventual game winner. One of the cool things about random forests is that they can identify which features are important based on how many trees they are found in. I created a random forest classifier and used it to predict the winner of the game at five minute intervals. Then for each timepoint, I extracted the feature importance, and plotted a heatmap:

Some features: Gold difference, kill difference, tower difference, dragon difference, and firsts of bunch of featues. Carry share is the max kills on a team divided by total team kills. 

As one might expect, gold difference is the most predictive feature, at every time point. Next most important are the kill and tower difference. Surprisingly, dragon difference is relatively unimportant. This is probably because killing dragons yields only a small advantage. If I were to put on my amateur designer hat, I'd recommend Riot decrease the number of dragons required to get the double dragon buff. Another note is that first blood, first dragon, and first tower are mostly uninformative.

To get a better look at how feature importances changes with time, let's take slices of feature importance at 20 and 35 minutes.

There are some significant changes in importance between these time points. Notably, barons and inhibitors gain importance at 35 minutes. To compensate for this, there is a loss of importance in the rest of the features. Once again, dragons are relatively unimportant, even at the late timepoints. Inhibitors may be ever so slightly more predictive than barons.

Accuracy over time

Given the above feature importance data, I reduced the number of features in the classifier to avoid overfitting the data, and only used the features for the differences, and the baron and inhibitor features. With these features, I then created a small random forest at each time point, and performed cross-validated predictions for the eventual winner of the game.

At the beginning of the game, the classifier isn't successful because there isn't enough information to work with. As the game goes on, the classifier improves, peaking at ~80% accuracy between 20-35 minutes. However, the accuracy drops for longest games. One possibility for this is that there are fewer games that last 40 minutes (less than 1000), which means there is less data to train the classifier. Another possibility is that games that last a long time are fairly close, and close games are harder to predict.

How winnable are surrendered games?

As one last fun thing to do, I wanted to investigate the games that were surrendered early. I separated the games into "early surrenders," and "close" games. To explore the surrendered games, I plotted the gold difference at twenty minutes. The gold difference is bimodal, with one team normally having a large gold lead:

I then created used a random forest to predict the winner of these games, and it did so with 99% accuracy. To see whether these surrendered games were really unwinnable, I trained a random forest on the close games, and then used that forest to predict the surrendered games. Here, it was only 94% accurate. My interpretation of this is that around 6% of surrendered games are winnable!

I have lots of ideas for how to expand on these initial results. I would like to train the model on low ELO games, and see how predictable those are in comparison to high ELO games. Riot has an API for professional games as well, and it would be fun to see if pro games are more (or less) predictable than amateur games. I could scrape games from different regions and compare them. One long term goal is to create a live predictor for LCS games as they happen. If you have any ideas, let me know in the comments!

Walk Along the Paper Trail: Garfield Gap

It's been three years since I did a Walk Along paper summary! Wow!

Recently in our journal club, we discussed a paper by Garfield et al from the Lowell lab. In discussion, some unusually interesting points were raised, and I'd like to think about them here.

Background

I've written about this before, but here is a quick refresher on hunger in the brain. Many types of neurons control metabolism, but the most famous are AgRP neurons in the arcuate nucleus of the hypothalamus. If you stimulate AgRP neurons with channelrhodopsin, or chemogenetically, you can make mice "voraciously feed," aka stuff food in their face. Notably, AgRP neurons are GABAergic, which means they are shutting down neurons in other areas.

While I've written about AgRP neurons before, I've generally ignored what AgRP itself is: a neuropeptide. AgRP does not bind to AgRP receptors, but is actually an endogenous antagonist for the melanocortin-4 receptor (MC4R), whose principle ligand is α-MSH. AgRP neurons cohabitate in the arcuate nucleus with POMC neurons that produce α-MSH; stimulation of POMC neurons can reduce feeding (albeit on a longer timescale than AgRP neurons).

AgRP neurons project to a lot of different brain areas. To see whether all of these projections can induce feeding, the Sternson lab stimulated axon terminals in a bunch of different brain regions, and found that stimulation of only some of these terminals (namely the PVH, aBNST, and LH) could induce feeding. One question Garfield et al wanted to answer was, "what is the molecular identity of these downstream targets?"

Is occlusion anything?

Since AgRP does not have any agonistic receptor, Garfield and colleagues investigated neurons expressing MC4R in various brain regions. They started by performing channelrhodopsin assisted circuit mapping (CRACM!) to see if AgRP neurons connect to MC4R neuons. To do this, they infected AgRP^Cre :: MC4R^Cre mice with channelrhodopsin in the AgRP neurons, and GFP in the MC4R neurons in the PVH. They sliced brains, patched onto PVH^MC4R neurons, and photostimulated the AgRP axon terminals (see diagram below). They found that 25/30 MC4R neurons received IPSCs following photostimulation, showing they were connected to AgRP neurons (panel A). They also patched onto non-MC4R neurons in the PVH, and found that only 2/10 neurons received AgRP input, showing that the AgRP-> PVH connection was fairly specific for MC4R neurons (panel B).

Patch recordings of MC4R neurons downstream from AgRP. AgRP neurons express ChR2-mCherry (red), and MC4R neurons express GFP (green). A. When you photostimulate AgRP neuron terminals, most MC4R neurons receive GABAergic inputs. B. When you photostimulate AgRP input to Non-MC4R neurons, most neurons do not receive input.

A previous paper by Atasoy and Sternson had claimed that AgRP neurons project to oxytocin or SIM1 neurons in the PVH, so Garfield investigated a few other neuron types in the PVH as well, but found no connections to any of them.

After showing the connection in-vitro, they wanted to show it functionally in-vivo using behaviour. They performed an occlusion study where they infected both AgRP and PVH^MC4R neurons with ChR2, then put an optic fibre over the PVH to stimulate the AgRP fibre terminals simultaneously alongside PVH^MC4R cell bodies (panel g, below). When they did this, they found the food intake was lower than AgRP neuron stimulation alone (panel h).

g. Diagram of experiment. AgRP neurons express ChR2. In different experiments, MC4R or OXT neurons also express ChR2. The optic fibre is placed over the PVH to stimulate cell bodies and AgRP terminals. h. Stimulation of AgRP terminals increases feeding. This is reduced ("occluded") by simultaneous stimulation of MC4R neurons. It is NOT reduced by simultaneous stimulation of OXT neurons.

I originally liked the idea of behavioural occlusion, perhaps because it reminded me of LTP occlusion. However, I'm not sure that it's informative in circuit mapping. First, my intuition is that direct excitation beats indirect inhibition. So if you stimulate AgRP terminals and MC4R neurons at same time, and direct excitation wins, it doesn't really tell you anything. Second, if you know two brain regions oppositely modulate behaviour, stimulating both of them does not tell you whether they directly interact. For example, stimulation of PKC-delta neurons in the CeA reduces feeding, and PKC-delta neurons are not connected to AgRP neurons. If you simultaneously stimulate AgRP and PKC-delta neurons, and one "occludes" the other, it won't mean they are directly connected. It only means one is stronger than the other!

In fact, I think the term "occlussion" is misleading, and not used in the same way as it was in LTP. In LTP, two protocols "occlude" each other if they both induce LTP alone, but stimulating both of them does not yield additional LTP. They are said to occlude each other because they use the same signaling pathway. In behaviour, "occlusion" has referred to stimulation of opposing pathways where one wins out. This experimental paradigm seems to be catching on, but I'm not sure it actually means anything.

Do all neuropeptides synapse on receptors?

Garfield next started looking for other AgRP-MC4R connections in other brain regions. As before, they infected AgRP neurons with ChR2-mCherry, and patched onto GFP-infected MC4R neurons in the anterior BNST, and the lateral hypothalamus (LH). Unlike before, there were no connections between AgRP neurons and MC4R neurons in these other brain regions.

Whole-cell patch onto neurons MC4R downstream from AgRP neurons expressing ChR2. No neurons, either in the aBNST or LH, were connected to AgRP neurons. Note the n's are not connected.

This raises question to me, how often do peptide neurons synapse onto peptide receptor neurons? (AgRP is certainly a strange case insofar as it doesn't have a natural receptor). There are hundreds of peptide-receptor pairs, and from the brief literature I've looked at, people don't always verify these cells are connected. For example, in the recent Dietrich paper about NPY5R (NPY is another neurotransmitter for AgRP neurons), they stimulated AgRP neurons and applied NPY5R antagonists, but never actually patched onto neurons to see if they were connected.

All neuropeptide neurons express multiple neurotransmitters, including classic ones and neuropeptides, and it's possible each transmitter work at different temporal and spatial scales. Fast neurotransmitters like glutamate or GABA are reuptaken quickly, and so cannot diffuse; in contrast, neuropeptides can last for minutes or hours. This would allow, for example, AgRP neurons to form GABAergic synapses on specific targets, and let their neuropeptides diffuse and act as paracrine [?] signaling.

In any case, I hope that as people explore these systems, they verify that the neurons we assume are connected, in fact are.

Collateral

Once they identified PVH^MC4R neurons as important for feeding, they wanted to know where they projected, and specifically if PVH^MC4R projected to multiple targets or single ones. They identified PBN as a major target of PVH^MC4R neurons, and used rabies virus to retrogradely label PVH^MC4R that project to the PBN. When they investigated other brain regions that PVH^MC4R projects to, they did not see axons, showing these neurons do not send collaterals.

(Top) Diagram of experiment. PVH^MC4R neurons are labeled in red. PVH^MC4R that project to the PBN are labeled in green. (Images) There are red and green cells in the PVH, and red and green terminals in the PBN. However, there are only red terminals in the vlPAG and NTS/DMV.

Previous research has shown that AgRP neurons do not have collaterals, and made me wonder whether this is a common feature of mid- or hind-brain neurons. However, a recent paper from Luqin Luo's lab mapped axons of locus coeruleus norepinephrine neurons, and found that those neurons do send collaterals widely, albeit biased towards specific areas. As we get more information about cell types, and perform more extensive tracing studies, we will get a better sense of what parts of the brain have discrete pathways versus overlapping projections.

The catch-22 of recording from identified cell populations

Recording from individual cells in genetically identified populations is the hottest technique in systems neuroscience right now (I am, of course, totally biased since that's what I'm trying to do). To record from identified populations, you first choose a mouse line that expresses Cre driven by a cell-specific marker like D1R. Then you transduce those cells with floxed ChR2 or GCaMP so you can phototag or image them. Finally, you run the mice through behaviours while recording to see how the identified neurons respond. If everything works out, you can directly link identified neurons' activity with specific behaviours.

There is, though, a catch-22 in interpreting these results. If all the cells respond the same way, then you have confirmed that your population is unitary, but the single cell nature of your recording yields no new information. On the other hand, if the cells you record from respond to a bunch of different things, they're not really a unified population, but rather a conglomeration of different neuron types. To illustrate these issue, I'd like to briefly show figures from two papers.

All together now

If you stimulate AgRP neurons in the hypothalamus, you can get a mouse to shove chow in its face. What remained unclear until this year, however, was how these neurons fire naturally. Some evidence gave hints; cFos staining had shown that these neurons are more active in hungry mice; in vitro recordings have shown that in hungry mice, these neurons receive more excitatory input, a cool form of short-term plasticity (Yang, ..., Sternson, 2011); and imaging has shown that these neurons undergo rapid spinogenesis and pruning when mice are hungry and fed (Liu, ...,  Lowell, 2012). In general, the working hypothesis was that AgRP neurons fire at a high rate when a mouse is hungry, which causes a mouse to seek food, or eat; when a mouse is sated, AgRP neurons turn off.

Given that basic model, there were many unanswered questions. How fast do AgRP neurons turn on and off? Do they turn off when you start eating, or do they take time to integrate enteric (gut) signals? What rate do they fire at? To answer these types of questions, we needed the development of easier in-vivo recording techniques for deep brain areas.

Earlier this year, the Knight lab at UCSF answered many of these questions by doing fibre photometry of AgRP and POMC neurons expressing GCaMP6s. (Chen, ..., Knight, 2015). They found that the activity of AgRP neurons in hungry mice actually decreases before mice start eating, when the mice first sense food (see below). In addition to receiving gut information, AgRP neurons receive fast input from brain areas that can identify food, which was unexpected. These results also question whether AgRP neurons are "hunger" neurons, or something slightly different like food seeking neurons.

AgRP neurons decrease their activity when mice see food. B. When a hungry mouse sees food (red trace), AgRP neuron fluorescence decreases. When it see an object (black trace), the neuron does not change activity. C. Summary of AgRP neuron fluorescence changes in fasted and fed mice responding to objects and chow. Only fasted mice that see chow decrease fluorescence.
From Chen, et. al., 2015.

Even knowing how AgRP neurons respond on average, it is possible that individual AgRP neurons respond differentially due to differences in protein expression or projections. AgRP neurons express variable levels of metabolic receptors like insulin or ghrelin, which could influence firing. Subsets of AgRP neurons project to different brain regions without collaterals; stimulation of some of these projections induces feeding while others do not, implying functional differences. Given these differences, it's possible some AgRP neurons respond to food cues, while other respond to enteric signals. To answer these questions, you would need to record from single cells.

Fortuitously, the Sternson lab did just that, using an in-vivo endomicroscope to image individual AgRP neurons expressing GCaMP6 (Betley, ..., Sternson, 2015). They found that all AgRP neurons act pretty much the same. They quantified fluorescence changes from AgRP neurons when a mouse was well fed, or food deprived; 54/61 neurons had brighter fluorescence when the mouse was food deprived (panel e, below). Like Chen et. al., they found that AgRP neurons decreased their activity before the mice started to eat (panel f); 96% of them to be precise (panel i).

AgRP neurons have higher activity when a mouse is hungry. e. Magnitude of fluorescence events for ad-libitum (AL) and food restricted (FR) mice. f. Single trial fluorescence traces from neurons during consumption. Most neurons decrease their activity before food is consumed. i. Magnitude of fluorescence before eating, after eating for one trial, and during satiety. 96% of traces decrease between the first trial base and first trial food.
From Betley, et. al., 2015.

I can tell you from personal experience these experiments are not trivial. The pessimistic interpretation of these results is that the technical prowess did not yield any new information. In talking to colleagues, however, that is probably unfair. Since we know that stimulation of different subsets of AgRP neurons can elicit different behaviour, the uniformity of their responses is surprising. This creates the question of how subsets of neurons which express the same protein and respond to the environment in the same way can elicit different behaviour.

How am I different?

In contrast to all the neurons acting the same, there is the possibility that all the neurons act differently. To illustrate this, I've selected a recent paper from the Stuber lab at UNC which investigated GABA neurons in the lateral hypothalamus (LH; Jennings, ..., Stuber, 2015). They used an in-vivo endoscope to image Vgat-Cre neurons expressing GCaMP6m. They had previously shown that these neurons are involved in consummatory behaviour.

During imaging, they ran the mice through two sets of behaviours. First, they let the mice eat in a cage with food located in two corners. If neurons had increased activity in the food zone, they were categorized as food zone excited (FZe); if they decreased activity, they were food zone inhibited (FZi). Around 10-15% of neurons were FZi or FZe (panel F, below). In a second set of behaviours, the mice were taught a progressive ratio task (PR3) where they could nose poke for food. Here they found some neurons responded during the nose poke (23%), while others responded during the consumption of food (10%; panel H, below). Finally, one might imagine that FZe or FZi neurons were correlated with nose poke or consumption activity, so they investigated the overlap between these populations. 28/40 FZe neurons responded during the PR3 task, split between consumption and nose poke; 12/40 FZi neurons responded during the PR3 task, again split between consumption and nose poke (panel J).

LH GABA neurons do all sorts of stuff. F. Histogram of change in neuronal activity when a mouse enters a food zone. H. Venn diagram of cell responses during PR3 task. J. Overlap of FZe and FZi responses with PR3 responses. There does not appear to be a consistent pattern of activity either in the tasks individually, or across tasks.

The authors highlight the diversity of responses in their discussion:

Our data suggest that separate subsets of appetitive-coding and consumption-coding ensembles exist within the LH GABAergic network. Thus, the LH GABAergic network can be viewed as a mosaic of functionally and computationally distinct cell types, requiring further definition. Nevertheless, these important computational differences among individual LH GABAergic neurons would have gone unnoticed if only bulk neuromodulatory approaches were employed, further underscoring the necessity of identifying the natural activity dynamics within a network to better understand the precise neural underpinnings of complex behavioral states.

The pessimist in me is frustrated by these results. Now we have to track down the neuronal subtypes, and repeat the experiments for each subtype (and lord knows I hate repeating experiments)! On the positive side, this could be a building block for future experiments. There is now a lower bound for the number of cell types to look for (at least five). I think single cell sequencing is the only way to identify these cell types reliably, and without recording.

Anyway, that is a long way of getting at what I see as the catch-22 of single cell identified neuronal recording: if they're all the same, you didn't need single cell resolution; and if they're all different, you don't have an strongly identified cell type.

References

Betley, J., Xu, S., Cao, Z., Gong, R., Magnus, C., Yu, Y., & Sternson, S. (2015). Neurons for hunger and thirst transmit a negative-valence teaching signal Nature, 521 (7551), 180-185 DOI: 10.1038/nature14416

Chen Y, Lin YC, Kuo TW, & Knight ZA (2015). Sensory detection of food rapidly modulates arcuate feeding circuits. Cell, 160 (5), 829-41 PMID: 25703096

Jennings, J., Ung, R., Resendez, S., Stamatakis, A., Taylor, J., Huang, J., Veleta, K., Kantak, P., Aita, M., Shilling-Scrivo, K., Ramakrishnan, C., Deisseroth, K., Otte, S., & Stuber, G. (2015). Visualizing Hypothalamic Network Dynamics for Appetitive and Consummatory Behaviors Cell, 160 (3), 516-527 DOI: 10.1016/j.cell.2014.12.026

Liu, T., Kong, D., Shah, B., Ye, C., Koda, S., Saunders, A., Ding, J., Yang, Z., Sabatini, B., & Lowell, B. (2012). Fasting Activation of AgRP Neurons Requires NMDA Receptors and Involves Spinogenesis and Increased Excitatory Tone Neuron, 73 (3), 511-522 DOI: 10.1016/j.neuron.2011.11.027

Yang, Y., Atasoy, D., Su, H., & Sternson, S. (2011). Hunger States Switch a Flip-Flop Memory Circuit via a Synaptic AMPK-Dependent Positive Feedback Loop Cell, 146 (6), 992-1003 DOI: 10.1016/j.cell.2011.07.039

Playing with deconvolution and GCaMP6 imaging data

The Palmiter lab recently got an Inscopix microscope. We are still troubleshooting our surgeries and recordings right now, so we don't have any imaging data yet. Given that, I wanted to set up our analysis pipeline ahead of time. Specifically, I wanted to see how we can identify calcium events. In this post I will explore how well deconvolution works on calcium imaging data from the Svoboda lab.

Why deconvolution

Inscopix provides image analysis software, Mosaic, which is pretty good at motion correction, and identifying ROIs. For Mosaic's signal detection, it identifies events using the simple criterion:

ΔF / F > F₀ + 3SD

This works adequately for low frequency events on a low background, but completely ignores the magnitude of events, or what happens when events occur in quick succession. While I've done imaging before, I've never done deconvolution, so I thought it would be fun to try it out.

To get sample GCaMP6 imaging data, I turned to the CRCNS, which has a bunch of neuroscience datasets. The most relevant dataset was from the original GCaMP6 paper wherein they recorded calcium fluorescence in parallel with loose seal electrophysiology, to allow comparisons between calcium and spikes. For this example, I loaded the processed data from the GCaMP6s imaging dataset. If you don't care about loading data in python, you can skip the next section and head straight to Deconvolution time!

Loading Svoboda lab data into python

The data is stored in .mat files using the Svoboda lab data format. Since understanding someone else's data structure is as difficult as understanding their musical taste, I'll walk through how to load the data in python. To load the data, you can simply use the scipy.io.loadmat function, with the following arguments:

import scipy.io as scio

svoboda_data = scio.loadmat( 'data_20120521_cell5_003', squeeze_me = True, struct_as_record = False)

It is critical to use the squeeze_me and struct_as_record flags! If you don't, python will load a weird object that will crash your python instance. This took me an embarrassingly long time to figure out.

svoboda_data is a dictionary with four fields; the only one we care about is the 'obj' field, so we can extract it:

svoboda_obj = svoboda_data['obj']

svoboda_obj has a bunch of fields as well; the most important one is timeSeriesArrayHash.value, which contains the data. The value field is an array of five structs, number 0-4, which have:

0: average fluorescence of an ROI
1: average fluorescence of the surrounding neuropil
2: raw e-phys trace from the cell
3: high-pass filtered e-phys
4: detected spikes

The value struct has two fields we care about, .time (which has the... time) and .valueMatrix (which has the values of the fluorescence or voltage). Knowing this, we can load the fluorescence and spike data:

ca_data = svoboda_obj.timeSeriesArrayHash.value[0].valueMatrix

ca_time = svoboda_obj.timeSeriesArrayHash.value[0].time

spike_data = svoboda_obj.timeSeriesArrayHash.value[4].valueMatrix

spike_time = svoboda_obj.timeSeriesArrayHash.value[4].time

I wrote a small wrapper function in python to load full Svoboda lab experiments. You can then plot the calcium fluorescence and spike times on each other:

Fluorescence (blue trace, DF / F) and spike data (green) from the GCaMP6s data set, recorded May 15, 2012, cell 1.

For this cell, most spikes generated calcium events, and multiple spikes in a short time period generated larger calcium events. However, some spikes did not yield much (if any) calcium response. Calcium events, in general, had mostly exponential decays.

Deconvolution time!

Now that the data is loaded, we can deconvolve it. Thankfully, Alistair Muldal has implemented the fast non-negative deconvolution described in Vogelstein, 2010 (PyFNND) . To be honest, I only partially understand the deconvolution algorithm. It starts by assuming that calcium can be modeled by spikes that return to baseline with exponential decay, but how they turn that generative model into a fit went beyond me.

In any case, we can run the PyFNND deconvolution on our calcium trace (I've calculated the ΔF / F separately), and get both a fit of the calcium trace, and an imputed spike train:

from pyfnnd import deconvolve, plotting 

n_best, c_best, LL, theta_best = deconvolve(
df_f.reshape, dt=0.0166, verbosity=1, learn_theta=(0, 1, 1, 1,0) )

plotting.plot_fit(df_f.reshape(-1, 1).T, n_best, c_best, theta_best, 0.0166)

Top: Original fluoresence (blue) and fit fluorescence (red), given the spike heuristics shown below.
Bottom: Imputed spike heuristics for the calcium trace shown above. Note that n_hat does not give a spike train, but rather a spike heuristic.

I played around with the theta parameters to see if I could get a better fit, but the best results I got were with the default ones above.

To get a spike train, you need to threshold the spike heuristic n_hat. Here I used a threshold of 0.1 (to avoid getting false positive spikes), and plotted the imputed spike train, real spike train, and calcium signal.

The imputed spike train (red) for a threshold of 0.1, and the real spikes (green).

Given how often the real spikes did not generate clean calcium events, I feel like the deconvolution did a pretty good job identifying spikes. The deconvolved spikes were often slightly delayed compared to the actual ones (although this could probably be fixed).

To get a sense for how accurately the deconvoluted spike train matched the real spike train, I calculated the Victor-Purpura distance (VP distance) between the two. The VP distance calculates how many times you would have to insert, delete, or move a spike to turn one spike train into another. VP was handily implemented by Nicolas Jimenez in the python module fit_neuron. (Note, if you want to use this module you should download the .tar.gz from Github, as the pip install has a bug at the moment.) I also wanted to use this metric to get a better sense of the optimal threshold for n_hat, so I ran the VP distance for thresholds of n_hat ranging from 0 to 1. I used a cost of q=5, so that spikes would be inserted or deleted only if there was not a nearby spike within 200 ms.

from fit_neuron.evaluate import spkd_lib as spkd

vp_cost = 5

spk_times = ca_data.spike_time[ca_data.spike_train>0.5]

vp_dist = spkd.victor_purpura_dist(ca_data.fluor_time[n_best >thresh], spk_times, vp_cost)

Similarity between the imputed spike train and real spike train. The normalized VP distance is the VP distance divided by the total number of real spikes. The threshold is the threshold for n_hat. 

A normalized VP distance of 0 means all the spikes were the same for each spike trains; and a distance of 2 means all the spikes from one train had to be deleted and re-inserted to recreate the other spike train. I think a distance of 0.8 means that just over half the real spikes were shared with the imputed spikes. I ran the deconvolution for all 9 neurons in the GCaMP6s dataset, and the minimum normalized distance ranged from 0.73-1. Some cells were undermined by bursts of firing, which caused a lot of false positives in the deconvolution. Others were difficult due to random non-spike related calcium fluctuations. I tried some filtering or baseline subtraction to make the results better but could not do significantly better than the default settings

A few concluding thoughts. Python is sweet, and people have implemented a lot of useful algorithms in it. This allowed me to try things out that would have been otherwise impossible. In the future I will try to clean up python posts by using iPython notebooks.

Deconvolution works pretty well at turning large calcium events into number of spikes. The timing won't be exact, but if you're doing calcium imaging, precise timing (<50ms) doesn't matter anyway.

[The following comments reflect my understanding after a few days of coding / playing. If I make any mistakes here, hopefully my commenters can correct me!] Finally, GCaMP6s was sold as having single spike resolution, but in my hands it does not, at least not for single trials. In Fig. 3F of the GCaMP6 paper they claim to detect 99% of single spikes in visual cortex with a 1% false positive rate. When I read that, I thought they had done something similar to what I just did: given a calcium fluorescence trace, and no prior knowledge, try to identify as many action potentials as you can, and compare to ground truth. What they actually did was answer the question: if you know what a single spike's fluorescent trace is, can you tell the difference between known spike events and random noise? In real imaging, of course, you don't know what spike calcium events look like (it's certainly not a simple exponential!). Also, that false positive rate, while sounding stringent, is quite permissive: if you image at 60Hz, with a 1% false positive rate, you will detect false spikes at a rate of 0.6Hz; if your cell fires at 0.5Hz, as their example cells do, that's a problem.

Having said that, using deconvolution with GCaMP6 works well depending on the type of data you're interested in. If you are recording from low frequency neurons where each action potential can cause spiking in downstream neurons, missing half the spikes is important. But if your neuron fires has a high baseline firing rate, and individual spikes aren't too important, the combination works well.

A cheap source for 230 μm ferrules

UPDATE: Somewhat embarrassingly, and fortuitously, I hadn't researched enough ferrule suppliers before making this post. After contacting a couple more, I found a supplier that sells 230 μm ID ferrules for $1.5 / pc. It is the Shenzhen Han Xin Hardware Mold Co. They do not have the 230 μm ID ferrules listed on their Alibaba page, but you can contact them for a custom order. Our lab ordered a set of 10 from them, and they worked well. If you use a lot of ferrules, I suggest ordering a test set yourself (shipping is $50).

If you want an American supplier, you can use Kientec. They sell 225 uμm ferrules (FZI-LC-225) for $3.25 / pc with an order of 100. These ferrules fit the fibre a little more snugly.

(The below is now outdated.)

Tl;dr: I am trying to get a bunch of labs to pool money to get create a new source for 230 μm ferrules.

The Palmiter lab makes its own fibre optic cannulae for implantation, using the protocol described in Sparta and Stuber. This has saved us a lot of money, making each cannula cost less than $10. As we have ramped up our optogentics, however, the monthly cost has risen enough that I wanted to find a way to save more money. The biggest cost in cannula construction is the zirconia ferrule (1.25 mm OD, 230 μm ID) we use, which cost $5 a piece. This is in contrast to ferrules with smaller 100 μm bores, which cost $1 a piece.

After hearing about Alibaba in the news, I looked there, and found lots of manufacturers for 100 μm ID ferrules, but none for 230 μm ID ferrules. Curious, I contacted Huangshi Sunshine Photoelectric, Co. to see if they could make 230 μm ID ferrules. They responded that they could, but that they would have to make a new "model." Making the new model would entail a one-time fixed cost of $7500, but once the model is made, the ferrules will cost $0.95 / piece for orders of 1000. (I sent them the specs of the Precision Fiber ferrules. They recommended a slight reduction of the diameter of the epoxy concavity, from 0.9 to 0.8mm; otherwise they should be the same. Once we have the money, they will make a confirmation diagram.)

I would like to get a group of labs together to cover the cost of the $7500 model (labs at UW should be able to cover $1500). Once the model is made, if your lab is big enough, you will be able to order directly for Sunshine Photoelectric in batches of 1000. If you can't use that many, the Palmiter lab should be able to pool together a bunch of smaller orders, and then distribute them.

If your lab uses these ferrules, and would be interested in contributing to getting the model made, please e-mail me at map222 at uw.edu.

UPDATE: The Zweifel lab here is on board, so we only need $6000 more!

Where's Bregma?

OR: A fun game for every surgery!

If you're interested in any of the important parts of the brain (i.e. not cortex), you're going to need do to stereotaxic surgery to target the area you're interested in. And to do stereotaxic surgery, you need fiducial coordinates on the skull, which canonically are bregma and lambda.

You would think that the definitions of bregma and lambda are well known, but I found out last year that I had a completely wrong definition of lambda. So in the interest of clarity, here's a diagram of a mouse skull with four possible locations for bregma and lambda (nose to the bottom). Without looking at the caption, can you, dear reader, identify the correct locations?

Diagram of mouse skull (not to scale). The bregma, for mice, is where the coronal suture intersects the midline, point B. Lambda is defined as the intersection of the midpoint and a curve fitting the "lamboid structure" (whatever that is). Basically, it's where the rear sinus would intersect the midline if the rear sinus didn't curve, point D. It is NOT point C, as any reasonable person would guess.

Now that you know the actual definitions of bregma and lambda, you can impress your professors at happy hour, and do accurate surgeries! Except bregma, is almost NEVER that simple in practice. The coronal suture curves this way and that; the left and right sides almost never meet at the same place on the midline. Sometimes the midline curves too!

Just for fun, I've diagrammed some common bregmas you might encounter during surgery. Which of the points listed below do you think are bregma? My opinion is in the caption.

Realistic bregmas (not to scale).
Top left: bregma where the two sides do not meet at the midline. The bregma here is at point B.
Top right: bregma where the mouse's right coronal suture veers. Here I'd say Bregmas is at C.
Bottom left: bregma where the midline curves. This is a bit trickier, but I'd say point B.
Bottom right: bregma where the midline shifts. Here there is actually not enough information! I would follow the midline all the way to lambda to see if it shifts back.

Disagree? Tell me in the comments!

The (Near) Future of Cell Type Specificity

I have been growing more interested in genomics, so this quarter I took a class on new techniques in genomics.* What I learned is that the most important aspect of modern genomics is cost. Sequencing gets exponentially cheaper every year, which makes answering old questions less expensive, and opens up new experimental possibilities. For example for cancer, in the past doctors might have genotyped a tumor once, while doctors now can sequence tumors in real time to identify drug resistance.

* If you miss the days of scientific tables, I suggest you check out the supplement of any genomics paper.

So how can genomics be useful for neuroscience? Sequencing is most famously used on DNA, but can be used on RNA too (RNA-seq). This makes it possible to characterize transcriptomes using sequencing rather than use microarrays. Recently, with cost reductions and new techniques, this has enabled single neuron transcriptomics. If you believe than a neuron’s mRNA and proteins define it, this means we can characterize single neurons with RNA-seq, then look at populations of neurons to figure out how neurons can be grouped into cell types.

So today I'd like to write about what cell specificity is, review two papers on single neuron transcriptomics, and the implications that this has on how systems neuroscience looks at cell-specific circuits.

What is a cell type?

The brain is composed of billions of neurons, each of which do something different, but many of which do something similar. Given that protein expression determines what a cell does, I would define a theoretical "cell type" as a group of cells that express a similar set of proteins and encodes a similar type of information. (I say "similar type of information" since neurons of the same type could encode different versions of the same information. For example, I would say all direction selective retinal ganglion cells are the same cell type whether they encode up or down movement. If different directions express different proteins, then that would change things (which I believe may actually be true)).

While that is a simple theoretical definition of a cell type, there is also the issue of how we deal with cell types in practice. Right now, the most popular way to approach cell-type specificity is transgenic mice (while I poked fun at Nature titles last month, a majority of them were about specific circuits controlling behaviour). To understand how cell-type specificity is done, it would be useful to have an example, so let's look at how Haubensak et. al. in David Anderson's lab did it.

Haubensak et al were interested in the central amygdala, which is involved in conditioned fear. To find genes specific to the central amygdala, they performed a microarray analysis on the transcripts from the brain region (today, people would probably scour the Allen Brain Atlas). This analysis yielded a few candidate genes like PKC-delta, CRH, and enkaphalin. Once they had their candidate genes, they performed immunohistochemistry and patching to characterize the cells' expression patterns and electrical properties. When they were satisfied that PKC-delta neurons were interesting, they injected viruses containing ChR2 into the central amygdala of PKC-delta-Cre mice, and performed their main fear conditioning experiments.

This illustrates the functional definition of neuron types today: Cre-driver lines with viruses injected in a specific location.

This functional definition of a cell type is incomplete a number of ways. It ignores that multiple types of cells in a given region could express the same gene; that is, there could be two types of cells that express PKC-delta. It also ignores that gene expression changes during development (this is partially overcome by the temporal precision of viruses). And it ignores the importance of that gene for the function of the cell.

As a general rule, I would guess that the farther an "identifying gene" gets from the synapse, the less specific it is as a cell marker. The closest genes to the synapse would be neurotransmitters and receptors. For example, in the arcuate hypothalamus, there are cell types identified for the neurotransmitters AgRP, POMC, and KNDy, each of which segregate. Acknowledging that there might be subsets of these cells (e.g. those that project to different places), these cell types are probably well defined. The major exceptions to this neurotransmitter heuristic are probably the neurotransmitters glutamate and GABA, which are expressed so ubiquitously that they are non-specific. (This has not, of course, stopped people from publishing well with VGat or Vglut2 cell lines.)

Examples of cell types in arcuate nucelus as determined by neurotransmitter. The set of neurons expressing GABA includes three different subsets of neurons, an example of how GABA and glutamate are poor cell markers. 

Moving farther from the synapse, many interneurons can be identified by the types of calcium binding proteins they express (e.g. calretinin vs calcitonin). Given how important calcium regulation is, these cells might have different firing properties, and be functionally different. Even farther from the synapse, you can identify cells by secondary messengers like PKC-delta. To be honest, though, without a clear link between secondary messengers and actual function, I am skeptical whether they can truly be cell specific markers.

Besides identifying cell types based on their transcripts, you can subdivide cells based on the projection patterns as well. For example, the AgRP set pictured above is actually composed of multiple groups of cells the project independently to different brain areas.

At the moment, our ability to identify and manipulate cell types is lacking. We use brain atlases and microarrays to identify populations of cells, but don't really know the boundaries of each cell type. So how can we define them better?

Single neuron transcriptomics

My favourite answer is single neuron transcriptomics. The basic idea is that you can isolate the mRNA from a single neuron, amplify it via PCR, and then sequence it; this would give you a quantitative assessment of all the RNA that a single cell expresses. You then repeat this process on dozens or hundreds of cells chosen at random, and use clustering algorithms to see which cells have similar expression profiles, which would constitute a defined cell type.

A paper came out in Nature Biotechnology this year which illustrates the process in developing human cortex. They took tissue, dissociated it, and isolated individual cells using microfluidics. Once single cells were isolated in their little compartments, they were able to perform all the steps outlined above: lysation, amplification, and sequencing.

Clustering of cells based on single-cell transcripts. On the x-axis are individual genes. On the y-axis are individual neurons. The color of each rectangle is the gene expression in the cell. Neurons could be clustered into 4 groups. Some groups, like group 1, could be further subdivided.

From Fig. 3 in Pollen, et. al., 2014

Once they quantified expression for all the mRNA, they performed PCA to identify the most informative genes; then hierarchical clustering of cells using the 500 most PCA-loaded genes. Following clustering, they identified four major groups of neurons that corresponded to cells at different stages of cell division (see above). They were furthermore able to subdivide these clusters, like group 1, into smaller subsets at different stages of differentiation. Using this technique they found that two transcription factors, Egr1, and Fos, were more important in development and differentiation than previously realized.

While cell sorting is a great technique for developing cells, it is less useful for mature neurons which send neurites everywhere; if you were to try to shove a neuron in a cell sorting tube, you would rip off the axons and dendrites, and leak mRNA everywhere. To do single cell transcriptomics in adult neurons, James Eberwine's lab at UPenn developed a method they call TIVA (see figure below).

How TIVA works. "The TIVA tag is composed of several functional groups: biotin, Cy3, poly(A) tail binding 2′-F RNA poly(U) oligo (orange), photocleavable linker (PL), 2′-OMe RNA poly(A) oligo (yellow), Cy5, disulfide bond (S-S) and CPP."

1.) The CPP tag on TIVA allows it to penetrate cells; once TIVA is in the cytosol, CPP is cleaved. You can see which cells have taken up TIVA by monitoring Cy5 fluorescence while exciting Cy3 (Cy3/Cy5 is a FRET pair).

2.) To release the RNA binding sequence, shine light on a single cell, cleaving the PL linker. You can monitor this via the Cy3/Cy5 fluorescence ratio.

3) Once released, TIVA will bind to the polyA tails of mRNA.

4) You lyse the cell, then extract TIVA and associated mRNA via the biotin tag.

6) ???

7) Amplify and sequence

8) Profit!

Fig. 1 from Lovatt et al (2014)

They first verified that TIVA worked in cultures before turning to hippocampal slices. In slices they bath applied TIVA, and monitored its uptake in cells. (panel b, below). They then photocleaved the construct in a single cell (cell ii), and verified the cleavage using fluorescence (panel c, right). Finally, they sucked up the tissue in a pipette, extracted the cell, and sequenced. Using sequencing they were able to verify that they were sequencing neurons, as opposed to the surrounding tissue (panel d).

Single neuron sequencing in hippocampal slices.

TIVA, as presented in the paper, has a couple major downsides. First, it cannot be targeted to specific cell populations. This means that you would have to profile all the cells in a given region, rather than specific subsets your are interested in. If there are many small populations in a given area, identifying them would require large sample sizes. Second, the throughput for the technique appears to be one neuron per hemisphere; in the paper they sequenced on the order of tens of neurons. If you wanted to phenotype one hundred cells in a brain region, you would have to sacrifice fifty mice, which is a hefty burden. On the positive side, once you characterize a brain, you won't have to profile the region again.

Whatever the current caveats, I think single cell transcriptomics is the future. It will allow us to accurately profile cells entire transcriptome, and cluster then to truly identify cells that have similar identities.

The future

So what does this mean for you (or me), dear neuroscientist?

First, I would bet that whatever cell type you (or I) are working on is not actually a single homogenous cell type, but probably composed of multiple cell types. This should not stop you from doing experiments, since the experiments of today are always superseded by tomorrow; but if you encounter ambiguous results, this is likely one source.

Two, intersectional methods are going to be essential going forward. For example, the POMC neurons shown in the first figure can release glutamate, GABA, or neither. If you truly want to understand the subsets of POMC neurons, you would need to use cell lines specific for POMC and glutamate, and viruses that could recognize the combination. Experiments like this will certainly require more breeding, but will give a more precise answer. There is probably a nice faculty position for whoever can come up with an easy intersectional transgenic model.

Third, the future is nigh.There are grants on RePORTER using single-cell transcriptomics to characterize prefrontal cortex, retina and zebrafish, and the dorsomedial hypothalamus. Just in the couple of weeks I was mulling this blog post, a paper in Nature Neuroscience came out characterizing dorsal root ganglion cells. If you want to Better Know a Cell Type, now is probably the time to act.

Tl;dr: If you're a systems neuroscientist who hasn't kept abreast of genomics, just remember: single-cell transcriptomics is getting cheaper by the month, and is going to redefine our understanding of cell types in the next few years.

(If you would like more informed opinion, James Eberwine, who's lab developed TIVA, wrote a commentary in Nature Neuroscience earlier this year.)

References

Haubensak, W., Kunwar, P., ..., Anderson, D. (2010). Genetic dissection of an amygdala microcircuit that gates conditioned fear Nature, 468 (7321), 270-276 DOI: 10.1038/nature09553

Lovatt, D., Ruble, B., ..., Eberwine, J. (2014). Transcriptome in vivo analysis (TIVA) of spatially defined single cells in live tissue Nature Methods, 11 (2), 190-196 DOI: 10.1038/nmeth.2804

Pollen, A., Nowakowski, T., ..., West, J. (2014). Low-coverage single-cell mRNA sequencing reveals cellular heterogeneity and activated signaling pathways in developing cerebral cortex Nature Biotechnology, 32 (10), 1053-1058 DOI: 10.1038/nbt.2967