es and mef single - cell measurements ( 96 cells ) from islam et al were used .
the mouse embryo data was taken from deng et al , using the read alignments described in the manuscript . to identify a subset of genes that can be used to fit error models for a particular single - cell measurements , all pairs of individual cells belonging to a given subpopulation ( e.g. all mef cells )
to do so , the observed abundance a given transcript in each cell was modeled as a mixture of the  drop - out 
this way , the expression of a gene with observed rpm levels of r1 and r2 in cells c1 and c2 respectively was modeled as : { r1poisson(0)dropoutinc1{r1nb(r2)r2nb(r1)amplified}r2poisson(0)dropoutinc2 } the background read frequency for the dropout components was set at 00.1 .
the mixing between the three components was determined by a multinomial logistic regression on a mixing parameter m = log(r1 ) + log(r2 ) .
alternatively , the initial three - component segmentation can be determined based on a user - defined background threshold , which is a lot less computationally intensive .
the genes that were assigned to the  amplified  components were noted , and a set of genes appearing in the  amplified  components in at least 20% of all pair - wise comparisons of cells of the same subpopulation ( excluding the cell for which the model is being fit ) was used to fit the individual error models , as described below .
the expected expression magnitude of these genes was estimated as a median observed magnitude between all the cell measurements in which a gene was classified to be in the  amplified  component .
the aim of the 20% threshold is to have a sufficiently large number of measurements for a given gene so that the median expression magnitude estimate would be reliable , and the model parameters resulting from the fitting procedure correlate well for a range of values corresponding to 6 - 12 cells ( supplementary figure 3d ) . to fit an individual error model c for a measurement of a single cell c ,
the observed rpm values were modeled as a function of an expected expression magnitude , using the set of estimates for a subset of genes described in the previous paragraph . the rpm level rc observed for a gene in cell c
was modeled as a mixture of a  drop - out  and  amplified  components , as a function of an expected expression magnitude e : { rcnb(e)amplifiedrcpoisson(0)dropout } with the mixing parameter m = log(e ) , 0 = 0.1 .
for each cell the model c was fit using em algorithm based on the set of genes for which expected expression magnitudes have been obtained .
the resulting estimates of parameters for the negative binomial and concomitant ( mixing ) regression were used as a description of an error model c in the subsequent analysis .
following bayesian approach , the posterior probability of a gene being expressed at an average level x in a subpopulation of cells s , was determined as an expected value ( e ) : ps(x)=e[cbp(xrc,c ) ] where b is a bootstrap sample of s , and p(x|rc , c)is the posterior probability for a given cell c : p(x|rc , c ) = pd(x)ppoisson ( x + ( 1  pd ( x)pnb(x|rc ) , where pd is the probability of observing a drop - out event in a cell c for a gene expressed at an average level x in s , ppoisson 
 x and pnb(x|rc ) are the probabilities of observing expression magnitude of rc in case of a drop - out ( poisson ) or successful amplification ( nb ) of a gene expressed at a level x in a cell c , with the parameters of the distributions determined by the c fit . for the differential expression analysis , the posterior probability that the gene shows a fold expression difference of f between subpopulations s and g was evaluated as p(f)=xxps(x)pg(fx ) , where x is the valid range of expression levels .
the posterior distributions were renormalized to unity , and an empirical p value was determined to test for significance of expression difference . the results
if scde , deseq , cuffdiff2 and singlecellassay ( sca ) were benchmarked against an expression dataset by moliner et al . that measured bulk mef and es cells grown using the same suspension growth protocol as used by islam et al .. the ability to recover top 1000 genes showing highest expression difference in moliner et al .
was assessed using roc / auc ( figure 2c , supplementary figure 4 ) ranking genes by significance of differential expression as determined by different methods .
standard measure of the genome - wide similarity between two single - cell measurements was determined as a pearson linear coefficient on log - transformed rpm values .
genes that did not show expression signals in any of the cells were excluded from the analysis .
the bray - curtis similarity measure was also calculated on log - transformed values ( linear - based values showed lower performance ) .
direct drop - out  similarity measure aims to estimate pearson linear correlation excluding likely  drop - out  events in any given cell . to achieve that we evaluate average correlation across 1000 sampling rounds , in each round probabilistically excluding likely drop - out observations .
specifically , in each round , an observation of a given gene at an expression level x in a particular cell was substituted with a missing value with probability pd(x)k , where pd(x ) is the probability of a drop - out event in the current cell at an expression magnitude level x , and k=0.9 is additional factor ( to stabilize similarity measure in cases when drop out rates are very high in a given cell ) . the overall similarity between any two cells
was then calculated as an average ( across 1000 sampling rounds ) pearson linear correlation between log - transformed values of observations that are valid ( not missing ) in both cells .
the  reciprocal drop - out  similarity measure aims to reduce the impact of drop - out events on the pearson linear correlation measure by weighting down the contribution of genes that are not likely to be reliably measured in both cells .
for instance , if a gene was observed at a level x1 in the first cell , we will weigh its contribution by the likelihood that such level of expression can be reliably detected ( i.e. without drop - out ) in the second cell .
this kind of reciprocal weighting minimizes the contribution of discrepant ( i.e. amplified vs. drop - out ) measurements to the overall similarity .
similarity was calculated as a weighted pearson linear correlation on log - transformed rpm values , weighting the contribution of each gene by k(1pd1(x2))(1pd2(x2))+(1k ) , where pd1(x2 ) is a probability of observing a dropout event in cell 1 for an expression magnitude x2 at which the gene was observed in the cell 2 .
we find that both direct and reciprocal similarity measures show robust improvements in classification performance for a range of k values between above 0.85 ( see supplementary figure 3e ) .
all similarity measures do well when all 90 + cells and a complete gene set are considered . to provide a meaningful comparison we measured performance on more challenging classification problems based on partial data .
specifically , a subset of 20 random es and 20 mef single - cell measurements was sampled in each iteration .
furthermore , increasing fraction of top differentially - expressed genes was excluded from the analysis ( figure 2d , x - axis ) to pose a more challenging classification problem .
the fraction of correctly classified cells was determined based on the top - level split of the resulting clustering .
es and mef single - cell measurements ( 96 cells ) from islam et al were used .
the mouse embryo data was taken from deng et al , using the read alignments described in the manuscript .
to identify a subset of genes that can be used to fit error models for a particular single - cell measurements , all pairs of individual cells belonging to a given subpopulation ( e.g. all mef cells ) were analyzed using three - component mixture model .
to do so , the observed abundance a given transcript in each cell was modeled as a mixture of the  drop - out  ( poisson ) and  amplification 
this way , the expression of a gene with observed rpm levels of r1 and r2 in cells c1 and c2 respectively was modeled as : { r1poisson(0)dropoutinc1{r1nb(r2)r2nb(r1)amplified}r2poisson(0)dropoutinc2 } the background read frequency for the dropout components was set at 00.1 .
the mixing between the three components was determined by a multinomial logistic regression on a mixing parameter m = log(r1 ) + log(r2 ) .
alternatively , the initial three - component segmentation can be determined based on a user - defined background threshold , which is a lot less computationally intensive .
the genes that were assigned to the  amplified  components were noted , and a set of genes appearing in the  amplified  components in at least 20% of all pair - wise comparisons of cells of the same subpopulation ( excluding the cell for which the model is being fit ) was used to fit the individual error models , as described below .
the expected expression magnitude of these genes was estimated as a median observed magnitude between all the cell measurements in which a gene was classified to be in the  amplified  component .
the aim of the 20% threshold is to have a sufficiently large number of measurements for a given gene so that the median expression magnitude estimate would be reliable , and the model parameters resulting from the fitting procedure correlate well for a range of values corresponding to 6 - 12 cells ( supplementary figure 3d ) . to fit an individual error model c for a measurement of a single cell c ,
the observed rpm values were modeled as a function of an expected expression magnitude , using the set of estimates for a subset of genes described in the previous paragraph . the rpm level rc observed for a gene in cell c
was modeled as a mixture of a  drop - out  and  amplified  components , as a function of an expected expression magnitude e : { rcnb(e)amplifiedrcpoisson(0)dropout } with the mixing parameter m = log(e ) , 0 = 0.1 .
for each cell the model c was fit using em algorithm based on the set of genes for which expected expression magnitudes have been obtained .
the resulting estimates of parameters for the negative binomial and concomitant ( mixing ) regression were used as a description of an error model c in the subsequent analysis .
following bayesian approach , the posterior probability of a gene being expressed at an average level x in a subpopulation of cells s , was determined as an expected value ( e ) : ps(x)=e[cbp(xrc,c ) ] where b is a bootstrap sample of s , and p(x|rc , c)is the posterior probability for a given cell c : p(x|rc , c ) = pd(x)ppoisson ( x + ( 1  pd ( x)pnb(x|rc ) , where pd is the probability of observing a drop - out event in a cell c for a gene expressed at an average level x in s , ppoisson 
 x and pnb(x|rc ) are the probabilities of observing expression magnitude of rc in case of a drop - out ( poisson ) or successful amplification ( nb ) of a gene expressed at a level x in a cell c , with the parameters of the distributions determined by the c fit . for the differential expression analysis , the posterior probability that the gene shows a fold expression difference of f between subpopulations s and g was evaluated as p(f)=xxps(x)pg(fx ) , where x is the valid range of expression levels .
the posterior distributions were renormalized to unity , and an empirical p value was determined to test for significance of expression difference .
the results if scde , deseq , cuffdiff2 and singlecellassay ( sca ) were benchmarked against an expression dataset by moliner et al . that measured bulk mef and es cells grown using the same suspension growth protocol as used by islam et al .. the ability to recover top 1000 genes showing highest expression difference in moliner et al .
was assessed using roc / auc ( figure 2c , supplementary figure 4 ) ranking genes by significance of differential expression as determined by different methods .
standard measure of the genome - wide similarity between two single - cell measurements was determined as a pearson linear coefficient on log - transformed rpm values .
genes that did not show expression signals in any of the cells were excluded from the analysis .
the bray - curtis similarity measure was also calculated on log - transformed values ( linear - based values showed lower performance ) .
the  direct drop - out  similarity measure aims to estimate pearson linear correlation excluding likely  drop - out  events in any given cell . to achieve that we evaluate average correlation across 1000 sampling rounds , in each round probabilistically excluding likely drop - out observations .
specifically , in each round , an observation of a given gene at an expression level x in a particular cell was substituted with a missing value with probability pd(x)k , where pd(x ) is the probability of a drop - out event in the current cell at an expression magnitude level x , and k=0.9 is additional factor ( to stabilize similarity measure in cases when drop out rates are very high in a given cell ) . the overall similarity between any two cells
was then calculated as an average ( across 1000 sampling rounds ) pearson linear correlation between log - transformed values of observations that are valid ( not missing ) in both cells .
reciprocal drop - out  similarity measure aims to reduce the impact of drop - out events on the pearson linear correlation measure by weighting down the contribution of genes that are not likely to be reliably measured in both cells .
for instance , if a gene was observed at a level x1 in the first cell , we will weigh its contribution by the likelihood that such level of expression can be reliably detected ( i.e. without drop - out ) in the second cell .
this kind of reciprocal weighting minimizes the contribution of discrepant ( i.e. amplified vs. drop - out ) measurements to the overall similarity .
similarity was calculated as a weighted pearson linear correlation on log - transformed rpm values , weighting the contribution of each gene by k(1pd1(x2))(1pd2(x2))+(1k ) , where pd1(x2 ) is a probability of observing a dropout event in cell 1 for an expression magnitude x2 at which the gene was observed in the cell 2 .
we find that both direct and reciprocal similarity measures show robust improvements in classification performance for a range of k values between above 0.85 ( see supplementary figure 3e ) .
all similarity measures do well when all 90 + cells and a complete gene set are considered . to provide a meaningful comparison we measured performance on more challenging classification problems based on partial data .
specifically , a subset of 20 random es and 20 mef single - cell measurements was sampled in each iteration .
furthermore , increasing fraction of top differentially - expressed genes was excluded from the analysis ( figure 2d , x - axis ) to pose a more challenging classification problem .
the fraction of correctly classified cells was determined based on the top - level split of the resulting clustering .
 single - cell data provides means to dissect the composition of complex tissues and specialized cellular environments . 
 however , the analysis of such measurements is complicated by high levels of technical noise and intrinsic biological variability . 
 we describe a probabilistic model of expression magnitude distortions typical of single - cell rna sequencing measurements , which enables detection of differential expression signatures and identification of subpopulations of cells in a way that is more tolerant of noise . 
Online Methods
Datasets and initial abundance estimates
Fitting individual error models
Differential expression analysis
Comparison of differential expression performance
Similarity measures and subpopulation analysis
Supplementary Material