|Home | About | Journals | Submit | Contact Us | Français|
Single-cell gene expression analysis has contributed to a better understanding of the transcriptional heterogeneity in a variety of model systems, including those used in research in developmental, cancer, and stem cell biology. Nowadays, technological advances facilitate the generation of large gene expression datasets in high-throughput format. Strategies are needed to pertinently visualize this information in a tissue–structure related context, so as to improve data analysis and aid the drawing of meaningful conclusions. Here we describe an approach that utilizes spatial properties of the tissue source to enable the reconstruction of hollow sphere–shaped tissues and organs from single-cell gene expression data in three-dimensional space. To demonstrate our method, we used cells of the mouse otocyst and the renal vesicle as examples. This protocol presents a straightforward computational expression analysis workflow and is implemented on the MATLAB and R statistical computing and graphics software platforms. Hands-on time for typical experiments can be less than 1 h using a standard desktop PC or Mac.
Measurement of gene expression levels at single-cell resolution has proven to be an accurate tool for in-depth analyses of cellular differentiation1 and cancer development2, and in describing regulatory mechanisms of cell fate decision processes3. Typically, multiplex quantitative RT-PCR (qRT-PCR) or whole-transcriptome analysis (RNA-Seq) of individual cells can be conducted for many genes simultaneously, resulting in the accumulation of data that is complex and information-rich. Existing analysis strategies generally utilize basic differential gene expression calculations4, different clustering algorithms5, and dimension-reduction procedures5. Here, we describe a technique that amalgamates single-cell expression data with geometric modeling in the three-dimensional (3D) space. We utilize two expression datasets, one derived from qRTPCR studies of 267 individual cells from the mouse otocyst, the precursor of the vertebrate inner ear6 and the second from RNA-Seq analysis of 57 cells from the renal vesicle, the first polarized epithelial precursor of the nephron7. The otocyst is a transient structure in early inner ear development that has the morphology of a hollow sphere, similar to the renal vesicle. An envelope-like arrangement serves as a template for the overall organization of individual cells in both tissues. Genes with previously known spatial expression domains are used to calculate a blueprint that recapitulates the position of each cell in the context of an orthogonal coordinate system representing the major organ axes. This platform provides a novel way to efficiently communicate biological data, and it enables the analysis of multidimensional data in an accessible and informative format.
This protocol employs either qRT-PCR or RNA-Seq data collected from single cells and provides instructions on how to analyze and visualize the data in a representative 3D model of the analyzed spherical organ. The study of gene expression in individual cells affords a more accurate representation of cell-to-cell variations, in contrast to bulk measurements that only reflect the stochastic average of such expression8. A number of multivariate analysis techniques have been proposed to scrutinize large sets of data. Most of the existing approaches deploy various clustering algorithms to describe potential subpopulations and apply dimension-reduction protocols, such as principal component analysis (PCA), to resolve patterns of shared transcriptional identities9,10. Although these methods are fundamental tools in generally organizing the biological heterogeneity of tissue-derived cells, they fail to take into account the original spatial organization of the tissue or organ. Gene expression data that reflect spatially encrypted properties can provide much deeper insight in the effort to understand cellular processes than gene expression information alone. The control and function of many cellular processes are tightly linked to and affected by the cells' spatial distribution11. Single-cell gene expression analysis methods generally require the complete dissociation of tissues and organs. As a consequence, once individual cells are dissociated, spatial information is lost. The method described in this protocol enables the recovery of spatial information. Computational reconstruction of multicellular structures in which individual cells have defined positional parameters and assignment of gene expression values on a cell-to-cell basis presents an invaluable tool to characterize cellular identities in the context of their (micro)-environment. Our protocol provides geometric modeling of shell-like, spherically shaped organs and tissues by integrating gene expression data derived from single cells with computational dimension-reduction methods, such as PCA.
This protocol can generally be applied to a variety of tissues that fit the morphological requirements, including the otocyst, renal vesicle, optic vesicle12, seminal vesicle13, or Kupffer's vesicle14 to name only a few. We note that the native cellular organization can vary between species and developmental stages. Particularly, spherical organs with more than one cellular layer, such as the blastocyst, probably fall short of being adequately represented with this protocol. Nevertheless, we envision that — given the appropriate numerical equations — this method can be expanded and used widely to comprehensively reconstruct sphere-shaped tissues and organs based on quantitative mRNA expression data of individual cells, in high-throughput and high-resolution.
In-situ hybridization and immunohistochemistry enable expression analysis of mRNA and proteins. 3D reconstructions of expression domains identified with these methods can be conducted using stacks of microtome or optical sections15. However, the throughput of these 3D reconstructions is low and only a few genes can be tested in parallel. On the other hand, microarray platforms as well as population-based RNA deep sequencing techniques enable the simultaneous measurements of thousands of genes16, yet these measurements are generally performed in bulk cell populations. Our method combines the benefits of high-throughput gene expression data acquisition with restored spatial information at single-cell level (Figure 1). The resulting technology enables the analysis of intricate gene expression data within the 3D context of multicellular systems.
To our knowledge, no comparable techniques exist that utilize single-cell gene-expression data to generate a comprehensive, spatially delineated expression atlas in a quantitative manner. SINGuLAR, a computational platform based on the statistical programming language R, developed by Fluidigm, enables analysis of large-scale quantitative expression data using techniques such as PCA, hierarchical clustering, and violin plot diagrams; yet none of these techniques acknowledge the structural context of the cell-derived anatomical configuration, and they present the data in 2D format only (http://www.fluidigm.com/singular-analysis-toolset.html). Other approaches have been developed that rely on various image acquisition and computer-based reconstruction protocols and do not directly measure RNA levels in individual cells17. More recently, techniques have been developed to measure the RNA complement of the genome within cells of intact tissues with subcellular resolution18.
Structures that are morphologically more sophisticated than those covered in the Procedure can similarly be computed and their geometric modeling is only limited by the availability of appropriate mathematical calculations and equations to describe the object in the 3D space; these approaches are, however, not covered in the present protocol. Because the anatomical characterization by mathematical equations constitutes an integral component of the protocol, we note that the geometric formulation should always be carefully examined and revised on an experiment-to-experiment basis. Furthermore, the protocol highlights the visualization aspect of single cell data and provides simple subsequent quantitation measures. We anticipate that the Procedure will serve as a basis for the research community to further develop our approach, so that the practicality of it will be improved and refined. Such improvements and refinements include the implementation of alternative mathematical equations that enable the description of non-spherically-shaped and hollow organs in 3D space, adjustment of analytical parameters that are tissue-specific and data-collection-platform dependent, as well as the incorporation of functional features that aid in the quantitation process, such as the calculation of mean gene expression per anatomical domain. We envision that computational reconstruction of multicellular structures using single-cell gene expression data, coupled with tissue-contextual, automated quantitation features will transform the field of cell and developmental biology.
The chances of success of this protocol are highly dependent on both the gene selection (in case of qRT-PCR) and feasibility of mathematically describing the cellular structure of a tissue or organ in 3D space (Figure 2). The fact that a priori knowledge about a minimum set of genes is necessary to implement this protocol represents a challenge in some instances, as satisfactory expression data is not always available (Figure 3). Therefore, the Procedure may not always be readily applicable and traditional initial validation experiments (e.g., in situ hybridization) are recommended to address such deficits. We also note that the present Protocol is not readily applicable to specimens derived from tumor samples. In these specimens, transcriptional heterogeneity is often accompanied by major morphological differences that vary from case to case. Additionally, the lack of a structured organization, which can be typically found in developing organs, exacerbates the difficulty of establishing axis-confined domains and makes this approach rather impracticable.
The methodology described herein exploits PCA as the underlying core technique. PCA is a mathematical approach to reduce the dimensionality of a multivariable dataset by concurrently retaining as much information as possible19 (Figure 3b). The result is a transformed coordinate system in which newly identified variables (the principal components) are arranged as linear combinations of the original variable vectors (the genes). However, expression profiles do not always conform to a linear relationship between genes in which case conclusive data visualization and interpretation may be difficult. Alternative techniques, such as non-linear PCA or other non-linear dimension reduction algorithms may help supplementing this guide, if more refined transformations are needed20.
As with all applications that measure the abundance of RNA species, degradation effects can critically hamper the downstream analysis. We therefore remind the researcher to comply with appropriate RNA handling precautions. The present protocol assumes that proper and adequate quality control measures of raw data have been applied.
Figure 2 summarizes the workflow of the present protocol, beginning with the determination of the tissue structure that one seeks to analyze, followed by the design of assays (required for qRT-PCR approaches only), and concluding with the in silico representation in 3D space and subsequent quantitation. Commonly single cell qRT-PCR approaches target the transcriptome in a sequence-specific manner unlike RNA-Seq approaches where quantitative data of the whole transcriptome is measured. Therefore, experimental setups using qRT-PCR platforms require the design of specific primer pairs – called the assays – to amplify transcripts of interests. Naturally, this step is not necessary when using single cell RNA-Seq data. In experimental setups where qRT-PCR data is employed, the module `Assay Design' consequently presents a vital stage in successfully applying the here-described methodology.
A careful and thorough selection of genes is crucial in effectively visualizing the data in 3D, as the expression data directly affects the degree of information that is preserved in the first principal components (Figure 3a,b). It is important to include genes with known expression domains that provide spatial information, permitting the delineation of putative organ or body axes, such as dorsal and/or ventral or the distal and/or proximal axes (Figures 4, ,5).5). Therefore, prior knowledge of a finite number of genes, which we call `anchor genes', with well-defined gene expression motifs is essential to exploit the full potential of this integrative technique. The inner ear dataset, utilized in the example application of the approach, relies on two anchor genes, with well-described expression domains in the otocyst, to determine all three axes of the organ structure in 3D space (Figure 6). The accuracy of this model was validated by confirmation or extension of the spatial expression of 35 additional genes with known expression in the mouse otocyst6. In the case of the nephron precursors, the analysis of which is also covered in the Procedure, we initially utilized twelve anchor genes to distinguish between proximal and distal domains. Of these twelve markers, eight qualified as conclusive anchor genes, similar to the conclusions reported by Brunskill et al.7
In addition to genes expected to be detectable in only a subset of cells, the list of assays should contain positive control genes with ubiquitous expression, such as Actb and Gapdh, as well as markers with expected absent expression (indicating off-target cells) for general quality control purposes. Depending on the collection method used for cell enrichment, the addition of representative negative markers is of particular importance. It will allow the researcher to exclude potential contaminating cells that are not part of the organ or tissue of interest and could interfere with the downstream analysis. Finally, any number of transcripts can be added for which expression patterns are unknown. With regard to the 267 otocyst-derived cells, the present protocol is based on the analysis of 96 genes of which about 40 had prior defined expression information for the organ6. Examples for genes with known expression patterns are Oc90, Gbx2, and Lfng that were predicted to be present in only a portion of all analyzed cells21–23. Actb, Gapdh and enhanced green fluorescent protein (Egfp) reporter gene served as universally expressed control assays and Pax6, which is expressed in the hindbrain but not the otocyst, was used as a negative marker. For the present analysis of the 57 renal-vesicle-derived cells aimed at showing general applicability of the protocol, we included the following distal anchor genes as specified by the authors, Brunskill et al.7: Pou3f3, Dll1, Sox9, Dkk1, Papss2, Greb1, Pcsk9, Lhx1, Bmp2, as well as the proximal genes Cdh6, Wt1, Tmem100.
A detailed description of experimental steps prior to ready-to-use data accumulation that include cell isolation and raw data acquisition is not provided in this protocol and can be found elsewhere24,25. We employed fluorescence-activated cell sorting (FACS) to collect individual cells from the target population (e.g., Pax2-Cre fate-labeled otocyst and neuroblasts of E10.5 mouse embryos)6. Similarly, nephron progenitor cells of P4 mice were assembled using FACS7. Alternative methods can be employed to enrich for the desired cells, such as magnetic cell sorting26, laser capture microdissection27 or micropipette aspiration28. We also note that cell collection needs to be carried out in a timely manner. Pooling of cells that originate from different batches (i.e., embryos) has to be treated with caution to ensure that variation in gene expression indicates mainly spatial and not temporal differences6,29. Regardless of the procedure applied, the total number of cells collected should correlate with measures of organ or tissue complexity, such as the expected number of different cell types of the tissue one seeks to geometrically model, as well as the tissue or organ's overall volume. The larger and more complex the multicellular structure of the tissue or organ, the more cells are necessary to adequately represent it computationally. It is reasonable to expect that an increase in the number of cells included in the analysis corresponds to an increased accuracy of the mathematical image and improved degree of gene expression map resolution. For RNA-Seq this increase in cell numbers used for analysis can be achieved cost-effectively by barcoding nucleic acids of individual cells to facilitate multiplexing of hundreds of cells for sequencing29.
The present Protocol uses qRT-PCR data from 382 early inner ear cells generated on the Fluidigm Biomark HD platform6 from otocysts and neuroblast cells of mouse embryos. For 3D reconstruction, we focus on the 267 otocyst-derived cells identified by cluster analysis in Durruthy-Durruthy, R. et al.6. The RNA-Seq data from 57 early kidney cells was prepared using Fluidigm's C1 system30 for automated reverse transcription and cDNA amplification as described in Brunskill E.W. et al.7. We note that 57 cells may not be an adequate number to representatively approximate the renal vesicle in its entirety.
This encompasses procedures like normalization and data outlier removal. To date, no widely agreed measure has been proposed to normalize single cell gene expression data. In approaches based on qRT-PCR, normalization techniques successfully applied are cell-based2, gene-based1, or both3. It should be the researchers' decision to evaluate different normalization techniques to determine which one is the most appropriate for their experimental paradigm. In the present protocol, in which qRT-PCR-derived data are utilized, we implemented cell-specific normalization factors that account for all genes across all cells, an approach we utilized in the study that is the basis of this Protocol6. In short, in Durruthy-Durruthy, R. et al.6, for every cell, the median Log2Ex value across all genes was calculated. The difference between this cell-characteristic value and the mean of all median Log2Ex values was subtracted from all Log2Ex values.
We note that it is the researcher's decision to define outlier cells, as this decision relies on various parameters (e.g., tissue type and collection method) and no universally valid approach exists. We assume that proper quality control measures are applied to assess the integrity of the raw data.
We include this module to highlight the possible necessity to first identify putative distinct subpopulations that greatly differ from each other with respect to their global transcriptional identity, suggesting dissimilar tissue site origins. In the case of the inner ear precursor cells, we filtered out a population of neuroblast cells that were not part of the otocyst cell population6.
In the Procedure, we provide instructions on how to implement the algorithm with two different computational environments, MATLAB (in the Procedure) and R (in Supplementary Methods). Alternative statistical analysis and visualization platforms can be used as well, such as Gnu Octave (https://www.gnu.org/software/octave/). We strongly recommend consulting documentations and manuals of MATLAB and R, as this protocol does not provide a comprehensive overview of how to use these software packages. More exhaustive and detailed resources for MATLAB can be found at http://www.mathworks.com/help/matlab/ and http://www.mathworks.com/academia/student_center/tutorials/launchpad.html that include tutorials, videos, and Q&A sections. Tutorials for R are available at http://cran.r-project.org/manuals.html. Table 1 summarizes the MATLAB and R toolboxes used in this protocol. The present protocol contains updated algorithms than were used in the primary research publication6. We also want to remind the reader that suggested parameters depend on the tissue sample examined.
The step-by-step guide in the Procedure is designed for both the novice and expert in MATLAB or R. The goal is to provide a transparent and straightforward workflow that enables the researcher to address questions concerning single cell-associated gene expression data in the context of spatially construed features.
A general understanding of the morphology of the tissue from which the cells are isolated is required. Also, expression data for some genes (e.g., anchor genes) must exist that are necessary to (1) map putative body axes onto the model and (2) to validate the distribution of the cells using independent information. The herein mathematically described spherical characterization can be modified and adapted to different geometric systems such as cylindrical, cubical, and potentially more complex multi-segmental structures.
In order to familiarize themselves with the protocol, we suggest readers implement the Procedure using the example data provided and referred to below, using either of both available options (Step 2 and 3 option A versus Step 2 and 3 option B, see Procedure), before proceeding to analyze their own data.
In the Procedure, to demonstrate the workflow, we employed the data set from our recent study6, in which we analyzed 382 single cells from the E10.5 mouse inner ear anlage consisting of 267 otocyst cells and 115 otic neuroblasts. In this protocol, we focus on the 267 cells of the otocyst. The numerical data is presented on a log-scale as Log2Ex values6. In addition we include the analysis of single-cell gene expression data that have been acquired using an RNA-Seq approach7. Here, cells are derived from P4 mouse renal vesicles and RPKM (reads per kilobase per million mapped reads) values were log2-transformed. Linear-scaled zero-values were set to – 16, the next negative integer of the most negative log2-transformed value of the dataset. We aim to visualize the data in a way that allows data extraction in a spatial context. The computational basis to achieve this objective is PCA-based 3D data projection.
Multiple file formats can be used to import the data into the MATLAB or R platforms. The raw data format created by the Fluidigm Software is a comma-separated value file (CSV), which can be opened and processed with any text editor, as well as Microsoft Excel. Processed RNA-Seq data is commonly available in .TXT or .CSV file format.
Supplementary Table S2 of Durruthy-Durruthy, R. et al.6 contains the 382-cell (otocyst and neuroblast) data table in a pre-processed and summarized format in an Excel file (mmc2.xlsx, available from the journal's webserver)6. It is not the intent of this protocol to discuss the methods (e.g., measures to assess quality of single cell data, normalization) for obtaining ready-to-use data from recorded raw data, which should be determined at the researcher's discretion. The example RNA-Seq data7 can be downloaded from the Gene Expression Omnibus (GEO) database under accession number GSE59130. Alternatively, both files are included as Supplementary Data in this protocol. Within this zip file are Otocyst_data.csv and renalvesicle_data.csv, which include only the relevant data, comprised of expression values of the 267 otocyst cells and 57 renal vesicle cells, respectively. Download and save either of these files to the working directory.
Download and install the latest version of MATLAB from http://www.mathworks.com/products/matlab/. Please refer to the product manual for information regarding the download and installation procedure. Another useful resource for MATLAB can be found at: http://www.mathworks.com/videos/getting-started-with-matlab-68985.html?s_tid=main_tutorial_ML_rp. Please verify that the required MATLAB toolboxes are installed (Table 1).
Download and unpack the file Supplementary Software 1. Save all function files in the MATLAB working directory (the same directory where the data files are saved). Each program file is a function that is used throughout the protocol. For additional information about each file please type `help program file' at the MATLAB command line.
Download and install the latest release version of R from http://cran.r-project.org/. Consult the R installation and administration manual if necessary. A useful reference for R commands can be found at http://cran.r-project.org/doc/contrib/Short-refcard.pdf. On Unix-based operating systems (e.g., Linux and OSX), the R package that plots data in three dimensions requires that X11 be installed. Additional required packages are specified in the Procedure section and in Table 1.
Download and unpack the file Supplementary Software 2. Save all script files in the R working directory (the same directory where the data files are saved). For additional information about each file please type `?program file' at the command.
CRITICAL: Below is described the application of the protocol using MATLAB. Two alternatives are given: one that makes use of the otocyst dataset (Step 2 and Step 3, Options A), the other that makes use of the renal vesicle dataset (Step 2 and Step 3, Options B). An analogous step-by-step procedure that uses the software environment R is provided in Supplementary Methods. Users may choose either one based on their preferences.
CRITICAL: In the following, MATLAB commands will be preceded by a `>>' sign, which is the typical prompt in MATLAB, not part of the executed command. Each command has to be executed in the sequence specified in the protocol. A semicolon at the end of a command tells MATLAB not to display any output from the command. An individual command preceded by a `>>' sign may stretch continuously over multiple lines.
CRITICAL: Accurate approximation of spatial relationships between single cells in the three-dimensional space requires the inclusion of features that demonstrate anatomical asymmetries in expression that are axially manifested. Conclusive representation of the tissue or organ involves, prior to PCA, identifying gene markers with spatially associated information and dismissing genes with cell-to-cell expression level variability that does not correlate with cellular position. The following Procedure step 3 describes the selection of marker genes with known asymmetric expression distribution that selectively and broadly label cells of the particular domain one seeks to model (e.g., Oc90 for the dorsal otocyst domain). We refer to these markers as `anchor genes' (see Experimental Design section of Introduction for further information). In the script file DetermineConclusiveAG.m within the Supplementary Software 1 we list the conceptual steps of this strategy in detail in the programmer comments.
CRITICAL: This subsection of the Procedure is Optional. It can be implemented to examine the data variance that is retained in the first few principal components (Figure 3b). Cases in which the variance in the first three components is not proportionally higher than in the subsequent components, may result in less accurate representations of the tissue or organ. This is crucial as data projection onto a lower-dimensional subspace, which is defined by principal components, is always accompanied with some degree of information loss. A relatively large total variance value that is captured by the first three components in relation to the following components suggests that information loss is attenuated which in turn improves pattern recognition of the data. This section also enables the visual inspection of the data when projected onto the first three components. Identification of subgroups based on how the data spreads in this subspace may lead to non-homogenous distributions of cells when projected onto a spherical model. In both cases (low variance in first three components, subgroup recognition) an alternative selection of anchor genes is advised.
CRITICAL: Boundary formation between groups of cells is a crucial event in development, as it enables cells with distinct functions to be kept physically separated from each other. In addition, the shape of a boundary and its orientation with respect to the body axes can influence important patterning events of an organ as development progresses and control its morphogenic roadmap. Here, boundaries are computed on the basis of co-expression of all conclusive anchor genes for each domain. Specific details of how the borders between neighboring domains are established can be found in the programmer comments in the script file ExtractCellsForAxis.m within the Supplementary Software 1.
CRITICAL STEP: Depending on the number of axes determined previously (step 3) the sphere is partitioned into two, four or eight domains. The PlotCellsIn3D function outputs the assignment of cells to domains (hemispheres, quadrants or octants) in the variable cellAssignments, which is a vector of length equal to the number of cells containing the numeric assignments of cells to the different hemispheres (values are 1 or 2), quadrants (i.e., Dorsal-Ventral and Lateral-Medial axes), or octants on for instance the Dorsal-Ventral, Lateral-Medial and Anterior-Posterior axes. The assignments can be inferred from Figure 6. For brevity, we do not specify the entire code of this function. Its main functional principles are documented in detail in the file PlotCellsIn3D.m within Supplementary Software 1. We also postulate that all three major body axes are aligned perpendicular to each other. This will facilitate visualization, characterization and quantification of the mathematical projection.
Steps 1 – 2: Importing data into MATLAB workspace: 10 min.
Steps 3 – 6: Determine anchor gene correlated markers for principal component analysis: 5 min.
Steps 7 – 9 (Optional): Perform principal component analysis and examine variance distribution across first principal components: 10 min.
Steps 10 – 13: Establish boundaries and project the data onto 3D space: 10 min.
Step 14: Quantitation of the data: 5 min.
Step 15: One-line command to plot expression data in 3D space: 2 min.
Principal component analysis is a mathematical procedure that aims to reduce the number of informative variables into a smaller set. As a measure of biological variability, the proportion of the total variance retained within each principle component is maximized and PCs are ranked in descending order (PC1 with highest variance, PC2 with second highest variance, etc.). In our setting, the first three components comprise 49% (otocyst) and 28% (renal vesicle) of the total variance when PCA is performed with all genes whose expression differences across cells are presumably based on spatial information (Figure 3a,b). We cannot provide a definitive figure that would enable the researcher to estimate whether the subspace that is spanned by the first three components retains `enough' data variability. Additional studies in the future that will produce single-cell gene expression data of spherical organs will likely contribute to a more determined assessment of this measure. Therefore, it is up to the researcher to decide whether a sufficient portion of the total variance is accounted by the first three components. For the two example data sets used in the protocol, we quantitated the expression distribution of various genes and found that their distribution aligned well with previously defined expression domains in the developing organ as described in the literature for the ear6 and renal vesicle7 (Figure 5). Based on these findings, we concluded that the three-component-subspace generated by PCA captures an adequate amount of data variability. If the first three PCs account for an insufficient portion of the total variance, it may exacerbate the difficulty of recognizing patterns and of performing quantitative downstream evaluations. To overcome such shortcomings, a revised gene collection may be needed.
The graphical output obtained at the end of the Procedure should architecturally resemble the crude geometry of a shell-like spherical object (Fig. 4). Morphologically more complex organs and tissues will demand alternative mathematical approaches, as outlined in the Limitations section.
Simple quantification methods based on different segmentation measures (e.g., other than Apical-Posterior, Dorsal-Ventral, Lateral-Medial, Proximal-Distal distinction) of the object can also be applied to support different approaches to characterize gene expression patterns.
We thank all members of the Heller laboratory for comments on the manuscript. This work was supported by NIH grants DC006167 and DC012250 to S.H., by P30 core support (DC010363), by the Stanford Initiative to Cure Hearing Loss, and in part by FP7-Health-2013-Innovation a cooperative grant by the European Commission.
Contributions R.D.D. and S.H. designed the experiments. R.D.D. performed the experiments. R.D.D. and A.G. performed data analysis. R.D.D., A.G. and S.H. wrote the manuscript.
Competing financial interests S.H. is affiliated with Inception 3, Inc. The rest of the authors declare no competing financial interests.