A pairwise cytokine code explains the organism-wide response to sepsis.

Sepsis is a systemic response to infection with life-threatening consequences. Our understanding of the molecular and cellular impact of sepsis across organs remains rudimentary. Here, we characterize the pathogenesis of sepsis by measuring dynamic changes in gene expression across organs. To pinpoint molecules controlling organ states in sepsis, we compare the effects of sepsis on organ gene expression to those of 6 singles and 15 pairs of recombinant cytokines. Strikingly, we find that the pairwise effects of tumor necrosis factor plus interleukin (IL)-18, interferon-gamma or IL-1ß suffice to mirror the impact of sepsis across tissues. Mechanistically, we map the cellular effects of sepsis and cytokines by computing changes in the abundance of 195 cell types across 9 organs, which we validate by whole-mouse spatial profiling. Our work decodes the cytokine cacophony in sepsis into a pairwise cytokine message capturing the gene, cell and tissue responses of the host to the disease.

A flexible empirical Bayes approach to multivariate multiple regression, and its improved accuracy in predicting multi-tissue gene expression from genotypes.

Predicting phenotypes from genotypes is a fundamental task in quantitative genetics. With technological advances, it is now possible to measure multiple phenotypes in large samples. Multiple phenotypes can share their genetic component; therefore, modeling these phenotypes jointly may improve prediction accuracy by leveraging effects that are shared across phenotypes. However, effects can be shared across phenotypes in a variety of ways, so computationally efficient statistical methods are needed that can accurately and flexibly capture patterns of effect sharing. Here, we describe new Bayesian multivariate, multiple regression methods that, by using flexible priors, are able to model and adapt to different patterns of effect sharing and specificity across phenotypes. Simulation results show that these new methods are fast and improve prediction accuracy compared with existing methods in a wide range of settings where effects are shared. Further, in settings where effects are not shared, our methods still perform competitively with state-of-the-art methods. In real data analyses of expression data in the Genotype Tissue Expression (GTEx) project, our methods improve prediction performance on average for all tissues, with the greatest gains in tissues where effects are strongly shared, and in the tissues with smaller sample sizes. While we use gene expression prediction to illustrate our methods, the methods are generally applicable to any multi-phenotype applications, including prediction of polygenic scores and breeding values. Thus, our methods have the potential to provide improvements across fields and organisms.

Fine-mapping from summary data with the "Sum of Single Effects" model.

In recent work, Wang et al introduced the "Sum of Single Effects" (SuSiE) model, and showed that it provides a simple and efficient approach to fine-mapping genetic variants from individual-level data. Here we present new methods for fitting the SuSiE model to summary data, for example to single-SNP z-scores from an association study and linkage disequilibrium (LD) values estimated from a suitable reference panel. To develop these new methods, we first describe a simple, generic strategy for extending any individual-level data method to deal with summary data. The key idea is to replace the usual regression likelihood with an analogous likelihood based on summary data. We show that existing fine-mapping methods such as FINEMAP and CAVIAR also (implicitly) use this strategy, but in different ways, and so this provides a common framework for understanding different methods for fine-mapping. We investigate other common practical issues in fine-mapping with summary data, including problems caused by inconsistencies between the z-scores and LD estimates, and we develop diagnostics to identify these inconsistencies. We also present a new refinement procedure that improves model fits in some data sets, and hence improves overall reliability of the SuSiE fine-mapping results. Detailed evaluations of fine-mapping methods in a range of simulated data sets show that SuSiE applied to summary data is competitive, in both speed and accuracy, with the best available fine-mapping methods for summary data.

Lymph-derived chemokines direct early neutrophil infiltration in the lymph nodes upon Staphylococcus aureus skin infection.

A large number of neutrophils infiltrate the lymph node (LN) within 4 h after Staphylococcus aureus skin infection (4 h postinfection [hpi]) and prevent systemic S. aureus dissemination. It is not clear how infection in the skin can remotely and effectively recruit neutrophils to the LN. Here, we found that lymphatic vessel occlusion substantially reduced neutrophil recruitment to the LN. Lymphatic vessels effectively transported bacteria and proinflammatory chemokines (i.e., Chemokine [C-X-C motif] motif 1 [CXCL1] and CXCL2) to the LN. However, in the absence of lymph flow, S. aureus alone in the LN was insufficient to recruit neutrophils to the LN at 4 hpi. Instead, lymph flow facilitated the earliest neutrophil recruitment to the LN by delivering chemokines (i.e., CXCL1, CXCL2) from the site of infection. Lymphatic dysfunction is often found during inflammation. During oxazolone (OX)-induced skin inflammation, CXCL1/2 in the LN was reduced after infection. The interrupted LN conduits further disrupted the flow of lymph and impeded its communication with high endothelial venules (HEVs), resulting in impaired neutrophil migration. The impaired neutrophil interaction with bacteria contributed to persistent infection in the LN. Our studies showed that both the flow of lymph from lymphatic vessels to the LN and the distribution of lymph in the LN are critical to ensure optimal neutrophil migration and timely innate immune protection in S. aureus infection.

Not just fibrotic: endothelial-derived TGFß maintains contractile function and lymphatic muscle phenotype during homeostasis.

Cell-cell communication within the lymphatic vasculature during homeostasis is incompletely detailed. Although many discoveries highlight the pathological roles of transforming growth factor-beta (TGFß) in chronic vascular inflammation and associated fibrosis, only a small amount is known surrounding the role of TGFß-signaling in homeostatic lymphatic function. Here, we discovered that pharmacological blockade of TGFß receptor 1 (TGFßR1) negatively impacts rat mesenteric lymphatic vessel pumping, significantly reducing vessel contractility and surrounding lymphatic muscle coverage. We have identified mesenteric lymphatic endothelial cells themselves as a source of endogenous vascular TGFß and that TGFß production is significantly increased in these cells via activation of a number of functional pattern recognition receptors they express. We show that a continuous supply of TGFß is essential to maintain the contractile phenotype of neighboring lymphatic muscle cells and support this conclusion through in vitro analysis of primary isolated lymphatic muscle cells that undergo synthetic differentiation during 2-D cell culture, a phenomenon that could be effectively rescued by supplementation with recombinant TGFß. Finally, we demonstrate that lymphatic endothelial production of TGFß is regulated, in part, by nitric oxide in a manner we propose is essential to counteract the pathological over-production of TGFß. Taken together, these data highlight the essential role of homeostatic TGFß signaling in the maintenance of lymphatic vascular function and highlight possible deleterious consequences of its inhibition.NEW & NOTEWORTHY The growth factor TGFß is commonly associated with its pathological overproduction during tissue fibrosis rather than its homeostatic functions. We expose the lymphatic endothelium as a source of endogenous TGFß, the impact of its production on the maintenance of surrounding lymphatic muscle cell phenotype, and internally regulated mechanisms of its production. Overall, these results highlight the intricate balance of TGFß-signaling as an essential component of maintaining lymphatic contractile function.

Flow-dependent regulation of rat mesenteric lymphatic vessel contractile response requires activation of endothelial TRPV4 channels.

OBJECTIVES: The objective of our study is to evaluate the involvement of the transient receptor potential vanilloid 4 (TRPV4) in the alteration of lymphatic pumping in response to flow and determine the signaling pathways involved. METHODS: We used immunofluorescence imaging and western blotting to assess TRPV4 expression in rat mesenteric lymphatic vessels. We examined inhibition of TRPV4 with HC067047, nitric oxide synthase (NOS) with L-NNA and cyclooxygenases (COXs) with indomethacin on the contractile response of pressurized lymphatic vessels to flow changes induced by a stepwise increase in pressure gradients, and the functionality of endothelial TRPV4 channels by measuring the intracellular Ca2+ response of primary lymphatic endothelial cell cultures to the selective agonist GSK1016790A. RESULTS: TRPV4 protein was expressed in both the endothelial and the smooth muscle layer of rat mesenteric lymphatics with high endothelial expression around the valve sites. When maintained under constant transmural pressure, most lymphatic vessels displayed a decrease in contraction frequency under conditions of flow and this effect was ablated through inhibition of NOS, COX or TRPV4. CONCLUSIONS: Our findings demonstrate a critical role for TRPV4 in the decrease in contraction frequency induced in lymphatic vessels by increases in flow rate via the production and action of nitric oxide and dilatory prostanoids.

Characterizing and inferring quantitative cell cycle phase in single-cell RNA-seq data analysis.

Cellular heterogeneity in gene expression is driven by cellular processes, such as cell cycle and cell-type identity, and cellular environment such as spatial location. The cell cycle, in particular, is thought to be a key driver of cell-to-cell heterogeneity in gene expression, even in otherwise homogeneous cell populations. Recent advances in single-cell RNA-sequencing (scRNA-seq) facilitate detailed characterization of gene expression heterogeneity and can thus shed new light on the processes driving heterogeneity. Here, we combined fluorescence imaging with scRNA-seq to measure cell cycle phase and gene expression levels in human induced pluripotent stem cells (iPSCs). By using these data, we developed a novel approach to characterize cell cycle progression. Although standard methods assign cells to discrete cell cycle stages, our method goes beyond this and quantifies cell cycle progression on a continuum. We found that, on average, scRNA-seq data from only five genes predicted a cell's position on the cell cycle continuum to within 14% of the entire cycle and that using more genes did not improve this accuracy. Our data and predictor of cell cycle phase can directly help future studies to account for cell cycle-related heterogeneity in iPSCs. Our results and methods also provide a foundation for future work to characterize the effects of the cell cycle on expression heterogeneity in other cell types.

The Bayesian lens and Bayesian blinkers.

I discuss the benefits of looking through the 'Bayesian lens' (seeking a Bayesian interpretation of ostensibly non-Bayesian methods), and the dangers of wearing 'Bayesian blinkers' (eschewing non-Bayesian methods as a matter of philosophical principle). I hope that the ideas may be useful to scientists trying to understand widely used statistical methods (including confidence intervals and [Formula: see text]-values), as well as teachers of statistics and practitioners who wish to avoid the mistake of overemphasizing philosophy at the expense of practical matters. This article is part of the theme issue 'Bayesian inference: challenges, perspectives, and prospects'.

Bayesian multivariate reanalysis of large genetic studies identifies many new associations.

Genome-wide association studies (GWAS) have now been conducted for hundreds of phenotypes of relevance to human health. Many such GWAS involve multiple closely-related phenotypes collected on the same samples. However, the vast majority of these GWAS have been analyzed using simple univariate analyses, which consider one phenotype at a time. This is despite the fact that, at least in simulation experiments, multivariate analyses have been shown to be more powerful at detecting associations. Here, we conduct multivariate association analyses on 13 different publicly-available GWAS datasets that involve multiple closely-related phenotypes. These data include large studies of anthropometric traits (GIANT), plasma lipid traits (GlobalLipids), and red blood cell traits (HaemgenRBC). Our analyses identify many new associations (433 in total across the 13 studies), many of which replicate when follow-up samples are available. Overall, our results demonstrate that multivariate analyses can help make more effective use of data from both existing and future GWAS.

Estimating recent migration and population-size surfaces.

In many species a fundamental feature of genetic diversity is that genetic similarity decays with geographic distance; however, this relationship is often complex, and may vary across space and time. Methods to uncover and visualize such relationships have widespread use for analyses in molecular ecology, conservation genetics, evolutionary genetics, and human genetics. While several frameworks exist, a promising approach is to infer maps of how migration rates vary across geographic space. Such maps could, in principle, be estimated across time to reveal the full complexity of population histories. Here, we take a step in this direction: we present a method to infer maps of population sizes and migration rates associated with different time periods from a matrix of genetic similarity between every pair of individuals. Specifically, genetic similarity is measured by counting the number of long segments of haplotype sharing (also known as identity-by-descent tracts). By varying the length of these segments we obtain parameter estimates associated with different time periods. Using simulations, we show that the method can reveal time-varying migration rates and population sizes, including changes that are not detectable when using a similar method that ignores haplotypic structure. We apply the method to a dataset of contemporary European individuals (POPRES), and provide an integrated analysis of recent population structure and growth over the last ∼3,000 years in Europe.

Discovery and characterization of variance QTLs in human induced pluripotent stem cells.

Quantification of gene expression levels at the single cell level has revealed that gene expression can vary substantially even across a population of homogeneous cells. However, it is currently unclear what genomic features control variation in gene expression levels, and whether common genetic variants may impact gene expression variation. Here, we take a genome-wide approach to identify expression variance quantitative trait loci (vQTLs). To this end, we generated single cell RNA-seq (scRNA-seq) data from induced pluripotent stem cells (iPSCs) derived from 53 Yoruba individuals. We collected data for a median of 95 cells per individual and a total of 5,447 single cells, and identified 235 mean expression QTLs (eQTLs) at 10% FDR, of which 79% replicate in bulk RNA-seq data from the same individuals. We further identified 5 vQTLs at 10% FDR, but demonstrate that these can also be explained as effects on mean expression. Our study suggests that dispersion QTLs (dQTLs) which could alter the variance of expression independently of the mean can have larger fold changes, but explain less phenotypic variance than eQTLs. We estimate 4,015 individuals as a lower bound to achieve 80% power to detect the strongest dQTLs in iPSCs. These results will guide the design of future studies on understanding the genetic control of gene expression variance.

Nr4A1 modulates inflammation-associated intestinal fibrosis and dampens fibrogenic signaling in myofibroblasts.

Intestinal fibrosis is a common complication of the inflammatory bowel diseases (IBDs), contributing to tissue stiffening and luminal narrowing. Human nuclear receptor 4A 1 (NR4A1) was previously reported to regulate mesenchymal cell function and dampen fibrogenic signaling. NR4A1 gene variants are associated with IBD risk, and it has been shown to regulate intestinal inflammation. Here, we tested the hypothesis that NR4A1 acts as a negative regulator of intestinal fibrosis through regulating myofibroblast function. Using the SAMP1/YitFc mouse, we tested whether two pharmacological agents known to enhance NR4A1 signaling, cytosporone B (Csn-B) or 6-mercaptopurine (6-MP), could reduce fibrosis. We also used the dextran sulfate sodium (DSS) model of colitis and assessed the magnitude of colonic fibrosis in mouse nuclear receptor 4A 1 (Nr4a1-/-) and their wild-type littermates (Nr4a1+/+). Lastly, intestinal myofibroblasts isolated from Nr4a1-/- and Nr4a1+/+ mice or primary human intestinal myofibroblasts were stimulated with transforming growth factor-ß1 (TGF-ß1), in the presence or absence of Csn-B or 6-MP, and proliferation and ECM gene expression assessed. Csn-B or 6-MP treatment significantly reduced ileal thickness, collagen, and overall ECM content in SAMP1/YitFc mice. This was associated with a reduction in proliferative markers within the mesenchymal compartment. Nr4a1-/- mice exposed to DSS exhibited increased colonic thickening and ECM content. Nr4a1-/- myofibroblasts displayed enhanced TGF-ß1-induced proliferation. Furthermore, Csn-B or 6-MP treatment was antiproliferative in Nr4a1+/+ but not Nr4a1-/- cells. Lastly, activating NR4A1 in human myofibroblasts reduced TGF-ß1-induced collagen deposition and fibrosis-related gene expression. Our data suggest that NR4A1 can attenuate fibrotic processes in intestinal myofibroblasts and could provide a valuable clinical target to treat inflammation-associated intestinal fibrosis.NEW & NOTEWORTHY Fibrosis and increased muscle thickening contribute to stricture formation and intestinal obstruction, a complication that occurs in 30%-50% of patients with CD within 10 yr of disease onset. More than 50% of those who undergo surgery to remove the obstructed bowel will experience stricture recurrence. To date, there are no drug-based approaches approved to treat intestinal strictures. In the current submission, we identify NR4A1 as a novel target to treat inflammation-associated intestinal fibrosis.

Empirical Bayes shrinkage and false discovery rate estimation, allowing for unwanted variation.

We combine two important ideas in the analysis of large-scale genomics experiments (e.g. experiments that aim to identify genes that are differentially expressed between two conditions). The first is use of Empirical Bayes (EB) methods to handle the large number of potentially-sparse effects, and estimate false discovery rates and related quantities. The second is use of factor analysis methods to deal with sources of unwanted variation such as batch effects and unmeasured confounders. We describe a simple modular fitting procedure that combines key ideas from both these lines of research. This yields new, powerful EB methods for analyzing genomics experiments that account for both sparse effects and unwanted variation. In realistic simulations, these new methods provide significant gains in power and calibration over competing methods. In real data analysis, we find that different methods, while often conceptually similar, can vary widely in their assessments of statistical significance. This highlights the need for care in both choice of methods and interpretation of results.

Dispersal syndromes drive the formation of biogeographical regions, illustrated by the case of Wallace's Line.

AIM: Biogeographical regions (realms) reflect patterns of co-distributed species (biotas) across space. Their boundaries are set by dispersal barriers and difficulties of establishment in new locations. We extend new methods to assess these two contributions by quantifying the degree to which realms intergrade across geographical space and the contributions of individual species to the delineation of those realms. As our example, we focus on Wallace's Line, the most enigmatic partitioning of the world's faunas, where climate is thought to have little effect and the majority of dispersal barriers are short water gaps. LOCATION: Indo-Pacific. TIME PERIOD: Present day. MAJOR TAXA STUDIED: Birds and mammals. METHODS: Terrestrial bird and mammal assemblages were established in 1-degree map cells using range maps. Assemblage structure was modelled using latent Dirichlet allocation, a continuous clustering method that simultaneously establishes the likely partitioning of species into biotas and the contribution of biotas to each map cell. Phylogenetic trees were used to assess the contribution of deep historical processes. Spatial segregation between biotas was evaluated across time and space in comparison with numerous hard realm boundaries drawn by various workers. RESULTS: We demonstrate that the strong turnover between biotas coincides with the north-western extent of the region not connected to the mainland during the Pleistocene, although the Philippines contains mixed contributions. At deeper taxonomic levels, Sulawesi and the Philippines shift to primarily Asian affinities, resulting from transgressions of a few Asian-derived lineages across the line. The partitioning of biotas sometimes produces fragmented regions that reflect habitat. Differences in partitions between birds and mammals reflect differences in dispersal ability. MAIN CONCLUSIONS: Permanent water barriers have selected for a dispersive archipelago fauna, excluded by an incumbent continental fauna on the Sunda shelf. Deep history, such as plate movements, is relatively unimportant in setting boundaries. The analysis implies a temporally dynamic interaction between a species' intrinsic dispersal ability, physiographic barriers, and recent climate change in the genesis of Earth's biotas.

UNIFYING AND GENERALIZING METHODS FOR REMOVING UNWANTED VARIATION BASED ON NEGATIVE CONTROLS.

Unwanted variation, including hidden confounding, is a well-known problem in many fields, but particularly in large-scale gene expression studies. Recent proposals to use control genes, genes assumed to be unassociated with the covariates of interest, have led to new methods to deal with this problem. Several versions of these removing unwanted variation (RUV) methods have been proposed, including RUV1, RUV2, RUV4, RUVinv, RUVrinv, and RUVfun. Here, we introduce a general framework, RUV*, that both unites and generalizes these approaches. This unifying framework helps clarify the connections between existing methods. In particular, we provide conditions under which RUV2 and RUV4 are equivalent. The RUV* framework preserves an advantage of the RUV approaches, namely, their modularity, which facilitates the development of novel methods based on existing matrix imputation algorithms. We illustrate this by implementing RUVB, a version of RUV* based on Bayesian factor analysis. In realistic simulations based on real data, we found RUVB to be competitive with existing methods in terms of both power and calibration. However, providing a consistently reliable calibration among the data sets remains challenging.

Off-Target Effect of Lovastatin Disrupts Dietary Lipid Uptake and Dissemination through Pro-Drug Inhibition of the Mesenteric Lymphatic Smooth Muscle Cell Contractile Apparatus.

Previously published, off-target effects of statins on skeletal smooth muscle function have linked structural characteristics within this drug class to myopathic effects. However, the effect of these drugs on lymphatic vascular smooth muscle cell function, and by proxy dietary cholesterol uptake, by the intestinal lymphatic network has not been investigated. Several of the most widely prescribed statins (Atorvastatin, Pravastatin, Lovastatin, and Simvastatin) were tested for their in-situ effects on smooth muscle contractility in rat mesenteric collecting lymphatic vessels. Lovastatin and Simvastatin had a concentration-dependent effect of initially increasing vessel contraction frequency before flatlining the vessel, a phenomenon which was found to be a lactone-ring dependent phenomenon and could be ameliorated through use of Lovastatin- or Simvastatin-hydroxyacid (HA). Simvastatin treatment further resulted in mitochondrial depolymerization within primary-isolated rat lymphatic smooth muscle cells (LMCs) while Lovastatin was found to be acting in a mitochondrial-independent manner, increasing the function of RhoKinase. Lovastatin's effect on RhoKinase was investigated through pharmacological testing and in vitro analysis of increased MLC and MYPT1 phosphorylation within primary isolated LMCs. Finally, acute in vivo treatment of rats with Lovastatin, but not Lovastatin-HA, resulted in a significantly decreased dietary lipid absorption in vivo through induced disfunction of mesenteric lymph uptake and trafficking.

Inference and visualization of DNA damage patterns using a grade of membership model.

MOTIVATION: Quality control plays a major role in the analysis of ancient DNA (aDNA). One key step in this quality control is assessment of DNA damage: aDNA contains unique signatures of DNA damage that distinguish it from modern DNA, and so analyses of damage patterns can help confirm that DNA sequences obtained are from endogenous aDNA rather than from modern contamination. Predominant signatures of DNA damage include a high frequency of cytosine to thymine substitutions (C-to-T) at the ends of fragments, and elevated rates of purines (A & G) before the 5' strand-breaks. Existing QC procedures help assess damage by simply plotting for each sample, the C-to-T mismatch rate along the read and the composition of bases before the 5' strand-breaks. Here we present a more flexible and comprehensive model-based approach to infer and visualize damage patterns in aDNA, implemented in an R package aRchaic. This approach is based on a 'grade of membership' model (also known as 'admixture' or 'topic' model) in which each sample has an estimated grade of membership in each of K damage profiles that are estimated from the data. RESULTS: We illustrate aRchaic on data from several aDNA studies and modern individuals from 1000 Genomes Project Consortium (2012). Here, aRchaic clearly distinguishes modern from ancient samples irrespective of DNA extraction, lab and sequencing protocols. Additionally, through an in-silico contamination experiment, we show that the aRchaic grades of membership reflect relative levels of exogenous modern contamination. Together, the outputs of aRchaic provide a concise visual summary of DNA damage patterns, as well as other processes generating mismatches in the data. AVAILABILITY AND IMPLEMENTATION: aRchaic is available for download from https://www.github.com/kkdey/aRchaic. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

A simple new approach to variable selection in regression, with application to genetic fine mapping.

We introduce a simple new approach to variable selection in linear regression, with a particular focus on quantifying uncertainty in which variables should be selected. The approach is based on a new model - the "Sum of Single Effects" (SuSiE) model - which comes from writing the sparse vector of regression coefficients as a sum of "single-effect" vectors, each with one non-zero element. We also introduce a corresponding new fitting procedure - Iterative Bayesian Stepwise Selection (IBSS) - which is a Bayesian analogue of stepwise selection methods. IBSS shares the computational simplicity and speed of traditional stepwise methods, but instead of selecting a single variable at each step, IBSS computes a distribution on variables that captures uncertainty in which variable to select. We provide a formal justification of this intuitive algorithm by showing that it optimizes a variational approximation to the posterior distribution under the SuSiE model. Further, this approximate posterior distribution naturally yields convenient novel summaries of uncertainty in variable selection, providing a Credible Set of variables for each selection. Our methods are particularly well-suited to settings where variables are highly correlated and detectable effects are sparse, both of which are characteristics of genetic fine-mapping applications. We demonstrate through numerical experiments that our methods outperform existing methods for this task, and illustrate their application to fine-mapping genetic variants influencing alternative splicing in human cell-lines. We also discuss the potential and challenges for applying these methods to generic variable selection problems.

Correction: Visualizing the structure of RNA-seq expression data using grade of membership models.

[This corrects the article DOI: 10.1371/journal.pgen.1006599.].

Visualizing the structure of RNA-seq expression data using grade of membership models.

Grade of membership models, also known as "admixture models", "topic models" or "Latent Dirichlet Allocation", are a generalization of cluster models that allow each sample to have membership in multiple clusters. These models are widely used in population genetics to model admixed individuals who have ancestry from multiple "populations", and in natural language processing to model documents having words from multiple "topics". Here we illustrate the potential for these models to cluster samples of RNA-seq gene expression data, measured on either bulk samples or single cells. We also provide methods to help interpret the clusters, by identifying genes that are distinctively expressed in each cluster. By applying these methods to several example RNA-seq applications we demonstrate their utility in identifying and summarizing structure and heterogeneity. Applied to data from the GTEx project on 53 human tissues, the approach highlights similarities among biologically-related tissues and identifies distinctively-expressed genes that recapitulate known biology. Applied to single-cell expression data from mouse preimplantation embryos, the approach highlights both discrete and continuous variation through early embryonic development stages, and highlights genes involved in a variety of relevant processes-from germ cell development, through compaction and morula formation, to the formation of inner cell mass and trophoblast at the blastocyst stage. The methods are implemented in the Bioconductor package CountClust.