Linking single-cell measurements of mass, growth rate, and gene expression

Mass and growth rate are highly integrative measures of cell physiology not discernable via genomic measurements. Here, we introduce a microfluidic platform enabling direct measurement of single-cell mass and growth rate upstream of highly multiplexed single-cell profiling such as single-cell RNA sequencing. We resolve transcriptional signatures associated with single-cell mass and growth rate in L1210 and FL5.12 cell lines and activated CD8+ T cells. Further, we demonstrate a framework using these linked measurements to characterize biophysical heterogeneity in a patient-derived glioblastoma cell line with and without drug treatment. Our results highlight the value of coupled phenotypic metrics in guiding single-cell genomics. Electronic supplementary material The online version of this article (10.1186/s13059-018-1576-0) contains supplementary material, which is available to authorized users.


Figure S1
Figure S1 | Fluidic regimes for maintaining cell spacing in the sSMR (a) Schema�c of sSMR presented in Figure 1 deno�ng the array entrance region used for fluidic simula�on presented in (b) (dashed box outline). (b) COMSOL fluidic simula�ons demonstra�ng the loading (le�) and flushing (right) fluidic regimes described in Addi�onal File 1: Note S1. The imaging region used to trigger between each fluidic state is outlined (solid box). a. b.
c.       , and (f) BT159 cells treated with RG7388 (n = 14,138 genes). Genes are plo�ed in rank order where genes with highest posi�ve and nega�ve correla�ons with biophysical parameters are found at the le�-most and right-most por�on of the x axis, respec�vely. For each data set, a null distribu�on of correla�on coefficients was determined by finding the correla�on between gene expression and mass for randomly permuted data. A�er 10 itera�ons, we determined the average standard devia�on of these distribu�ons of correla�on coefficients. Any individual gene that had a correla�on coefficient with an absolute value greater than twice the standard devia�on (P<0.05, denoted by the dashed lines in the plots) was considered significant (red bars), all genes presented as blue bars fell below this threshold. The number of genes showing a significant posi�ve or nega�ve correla�on with the biophysical parameter of interest are shown in each plot.     Plot of mass versus mass-normalized growth rate (growth efficiency) for a subset of the FL5.12 cells depicted in Figure 2 that were captured downstream for scRNA-seq (n = 124). Points are colored by G1/S score. The "cell cycle G1/S phase transi�on" gene ontology term was found to be significantly enriched amongst genes ranked by correla�on with growth efficiency. To determine the G1/S transi�on scores for single FL5.12 cells, we found the average of mean-centered, z-score scaled expression values for the genes in the "cell cycle G1/S phase transi�on" gene ontology term found to correlate significantly with normalized growth rate (n = 13 genes, Addi�onal File 1: Figure S4, Addi�onal File 5:        Table S10, Methods).

Supplementary Notes Note S1 | Maintaining minimum cell spacing in mass sensor array
Loading single cells into the mass sensor array at a fixed, minimum spacing requires the implementation of active switching between two distinct fluidic states. Initially, equivalent pressures are applied to the upstream and downstream ports on the bypass channel leading in to the array (Additional File 1: Figure  S1, ports P1 and P3). In this "loading" configuration, all streamlines are directed into the array and therefore cells in the bypass channel will enter the array. An imaging region at the entrance to the mass sensor array (outlined in Additional File 1: Figure S1) is used as an indication of when a cell has been successfully loaded. Real-time optical peak detection within this region is used to switch from this loading fluidic state to a "flushing" regime wherein the upstream pressures (P1) is increased and the downstream pressure (P3) is decreased such that a vast majority of streamlines continue along the bypass channel with a small fraction entering the array. Because cells are of finite size and occupy several streamlines, they are directed along the bypass channel and not drawn in to the array. Importantly, during this process the pressure at the entrance to the mass sensor array is maintained at a fixed value, therefore any cells that have entered the array continue to flow at a constant speed. Therefore, although the volumetric flow rate is maintained across the array while flushing, no additional cells are loaded. After a desired amount of time has elapsed the system is automatically returned to the loading configuration to obtain the next cell for measurement.

Note S2 | Determining reproducibility of gene signatures related to mass and MAR
In order to determine the reliability and reproducibility of the linked biophysical and gene expression profiles, it was important to compare these signatures with additional results collected from independent experiments. For L1210 cells, single-cell gene expression profiles had previously been collected for cells with known times since division (TSD), a proxy for cell cycle progression [25]. We therefore hypothesized that the list of genes with expression levels that correlated significantly with single-cell mass (an alternative proxy for cell cycle progression) would show significant overlap with genes that correlated strongly with TSD. To determine the extent of this similarity, we constructed two test gene sets for gene set enrichment analyses: one which included genes with a significant positive correlation with cell mass and another which included genes with a significant negative correlation with cell mass (Additional File 1: Figure S6a, Additional File 2: Table S1). These gene subsets were compared to the full L1210 gene list measured previously, with genes ranked by how strongly their expression levels correlated with TSD. Genes with a significant positive correlation with mass were significantly over-represented amongst genes that showed a positive correlation with TSD in prior measurements (FDR<0.05). Similarly, genes with a significant negative correlation with mass were significantly over-represented amongst genes that showed a negative correlation with TSD (FDR<0.05). These results indicate that similar sets of genes are correlated with both TSD and single-cell mass, suggesting consistency between the measurements collected here and those collected previously.
Next, we sought to perform a similar comparison for FL5.12 cells. However, in contrast to L1210 cells, no single-cell gene expression measurements had been collected for these cells previously. We therefore conducted a second, independent experiment where single-cell mass and MAR measurements were collected upstream of scRNA-seq for FL5.12 cells (Additional File 1: Figure S6b,c). Using this independent data set, we generated full gene lists that were ranked by correlation strength with either mass or massnormalized MAR. Then we once again constructed test gene sets, this time containing genes from the original FL5.12 data set with significant correlations (both positive and negative) with either mass or massnormalized MAR (P<0.05). Following the same analysis described above, we found that gene sets correlating with both mass and mass-normalized MAR showed significant overlap between both replicate experiments (FDR<0.05). This once again demonstrates the reproducibility of the gene expression signatures that correlate with single-cell biophysical properties.

Note S3 | Additional cell cycle gene expression analysis
To further validate the cell cycle-related gene expression signatures that correlated with cell mass, we performed an additional set of analyses relying on single-cell gene expression data alone modelled after the cell cycle interpretation approaches presented by Macosko et al. and Kowalczyk et al. [2,24]. Briefly, we utilized gene lists found to be associated with the G1/S and G2/M stages of the cell cycle reported by Whitfield et al. [53]. Rough phase-specific scores were determined by calculating the average expression value (ln(TPM+1)) of genes from these lists that were detected in the data set of interest (L1210, FL5.12, or CD8+ T cells). The lists used for each phase were then filtered to only include genes that correlated strongly (R>0.3) with these rough scores. The average expression values across these remaining genes were then mean centered and divided by their standard deviation to yield the final G1/S and G2/M phasespecific scores. Each individual cell was assigned to either the G1/S or G2/M phases based on which gene list yielded the maximal score. For all cell types, the mass of cells assigned to the G2/M phase were significantly greater than those assigned to the G1/S phase (P<0.001, Mann-Whitney U Test, boxplots in Additional File 1: Figure S5). Furthermore, cells ranked by mass showed clear negative and positive relationships with G1/S and G2/M scores, respectively (heatmaps in Additional File 1: Figure S5). These results offer further evidence of a coordination between single-cell mass and cell cycle gene expression in L1210, FL5.12, and CD8+ T cells at various stages of activation.