Single-cell analysis of population context advances RNAi screening at multiple levels

Abstract

Isogenic cells in culture show strong variability, which arises from dynamic adaptations to the microenvironment of individual cells. Here we study the influence of the cell population context, which determines a single cell's microenvironment, in image-based RNAi screens. We developed a comprehensive computational approach that employs Bayesian and multivariate methods at the single-cell level. We applied these methods to 45 RNA interference screens of various sizes, including 7 druggable genome and 2 genome-wide screens, analysing 17 different mammalian virus infections and four related cell physiological processes. Analysing cell-based screens at this depth reveals widespread RNAi-induced changes in the population context of individual cells leading to indirect RNAi effects, as well as perturbations of cell-to-cell variability regulators. We find that accounting for indirect effects improves the consistency between siRNAs targeted against the same gene, and between replicate RNAi screens performed in different cell lines, in different labs, and with different siRNA libraries. In an era where large-scale RNAi screens are increasingly performed to reach a systems-level understanding of cellular processes, we show that this is often improved by analyses that account for and incorporate the single-cell microenvironment.

Figure 1

Cell-to-cell variability in comparative RNAi screens of virus infection. (A) Experimental and computational overview of the 34 small-scale and 7 druggable genome screens analysed in this study. Individual screens are denoted as a combination of the virus and cell line abbreviation (for example, DEN-1_KY for DEN-1 infection in HeLa Kyoto cells). (B) Bootstrapped hierarchical clustering of the average population context-determined patterns of virus infection for each screen (Supplementary information). Infection efficiencies per quantile bin of the multidimensional population context and cellular state space were calculated, and z-score log₂ transformed. Median values along all dimensions per screen were clustered (Supplementary information). (C) Heat map of virus infection-induced phenotypes (Supplementary information, and Supplementary Figure 4). Grey boxes indicate insufficient data points. (D) Heat map depicting average cell-to-cell variability infection pattern strength and reproducibility. The strength and reproducibility was calculated as the average of bootstrapped cell-to-cell variability pattern correlations observed between pooled populations of cells selected from different sets of perturbations per screen (Supplementary information). (E) Table of viruses with characteristic virus properties, virus names, families, and abbreviations (see for more information Supplementary Table 2). (F) A schematic outline of the cell population context and cellular state space model used throughout this study. Non-mitotic and non-apoptotic cells denote interphase cells. (G) Example scatter plots of the data used to calculate the pattern strength and reproducibility. SV40_CNX (green dots, correlation coefficient of 0.91) shows a very strong and reproducible pattern of cell-to-cell variability, whereas MHV_TDS (red dots, correlation coefficient of 0.05) shows no measurable pattern of infection. Individual dots represent infection indices for 2 bins of equal cellular state and population context and with at least 30 cells in each bin, but sampled from different sets of perturbations.

Figure 2

RNAi causes widespread changes in cell population context leading to predictable changes in virus infection. (A–C) siRNA-induced changes in the cell population context ((A) Local cell density distributions; (B) Cell size distributions; (C) Fraction of cells located in cell islet edges) are shown from an example small-scale RNAi screen. Individual coloured lines and dots represent values from individual populations of an example RNAi screen in HeLa MZ cells, colour coded according to the number of cells in each perturbed population (see colour scale on the right of panel C). Black lines (and transparent grey regions) and dots represent the median values (and s.d.) of non-targeting siRNA controls. (D) Two examples of wells with equal cell number (∼2600 cells) but different cell population contexts. Segmented nuclei are pseudo-coloured according to the local cell density, cell islet edges, and predicted DEN-1 infection probabilities. The upper panels show HeLa Kyoto cells transfected with the second siRNA targeting STK40, resulting in an average local cell density of 4.5 (a.u.), 32% of cells located on cell islet edges, and a predicted infection index of 0.091. The lower panels show HeLa Kyoto cells transfected with the third siRNA targeting DYRK3, which has an average local cell density of 2.2 (a.u.), of cells located on cell islet edges, and a predicted infection index of 0.248. (E) RNAi results for VSV infection expressing GFP in HeLa MZ cells. Normalized z-scores of fold-change in infection (see Supplementary information) are presented per siRNA perturbation. Values are ranked according to the strength in direct perturbations of the infection process inside single cells. Blue circles, direct effects (i.e. population context-independent) for each siRNA; Green circles, direct effects of non-targeting (scrambled) siRNAs; Red circles, direct effects of positive controls (siRNAs silencing GFP); Light blue squares, total siRNA effects (direct and population context-mediated effects); Black circles, population context model-expected infection indices. Median values of triplicates are shown. Error bars depict standard deviations of direct siRNA effects (n=3). (F) Correlation between all measured infection indices per siRNA and their model predicted infection indices for each small-scale screen (see Supplementary information and Supplementary Figure 13). Inset shows an example scatter plot of these values for DEN-1 infection in HeLa Kyoto cells (correlation=0.85). The bars on the right axis show the relative histogram of all correlation values.

Figure 3

Single-cell modelling of causality in population context-determined cell-to-cell variability of infection in RNAi screens. (A–C) Examples of the causal modelling (see Supplementary information) of direct (A), population context-mediated (B) and masked (C) perturbation effects, as well as their corresponding infection and population context measurements in various RNAi screens. Black arrowheads in the causal graphs denote inferred directionality of effects, whereas their absence indicates effects where both directions had the same likelihood. Red and green edges denote the sign of the pairwise single-cell bootstrapped correlations between nodes of the network (red=negative, green=positive). Orange nodes in the causal graphs represent virus infection; Light blue nodes represent siRNAs indicated by the gene name and their siRNA number (grouped by parentheses if effects for different siRNAs are identical). Grey nodes represent population context parameters (Pop. size, population size; Edge, cell islet edges; LCD, Local cell density). Population context and cellular state nodes not causally linked to infection have been omitted for clarity. The bar plots show normalized z-score values (Supplementary information) for infection (orange bars) and 3 × z-score values (to put them on the same scale as infection) for population context measurements (grey bars). Error bars denote s.d. (n=3) around median values. (D) Quantification of fully indirect (blue) and fully masked (red) siRNA effects for the YFV and VACV large-scale screens. siRNA phenotypes were classified based on their direct and total measurements, as indicated in the legend.

Figure 4

Accounting for population context-propagated effects changes hit lists and improves reproducibility, cell line comparability, and siRNA phenotype consistency. (A) Percentage of gene hits (observed with at least two independent siRNAs in triplicate) that changed after accounting for cell population context in: (left) the top-300 genes of the druggable genome screens, which silencing reduced virus infection; (second from left) the average of the top-3 genes in the 34 small-scale screens, which silencing reduced or increased infection; (third from left) the top-300 genes of the druggable genome screens, which silencing reduced interphase nuclei size; (right) the top-150 genes of the medium-scale RNAi screens on early endosome abundance (stained for EEA1). Brackets on top of the bar graphs indicate the improvement of hit-list overlap between screens of the same cellular activity performed in different cell lines and/or laboratories, with % overlap before and after the arrow representing results before and after accounting for population context-mediated effects. (B) Infection indices at intersections of the multidimensional population context and cellular state space (as in Figure 1F) for interphase cells, cells not located on cell islet edges, at intermediate population sizes. Examples for Vaccinia virus (VACV) infection in HeLa MZ cells (upper) and HeLa Kyoto cells (lower) are shown (see Supplementary Figure 20 for more examples). VACV infection preferentially occurs in different subpopulations of cells in the two cell line strains. (C) Scatter plots for measurements of total effects (upper) and direct effects (lower) of siRNAs on VACV infection in HeLa MZ and HeLa Kyoto cells. Black dots depict median values for individual siRNAs (n=3). Values have been z-score log₂ normalized over all siRNA phenotypes. Total siRNA effects do not significantly correlate between cell lines (ρ=0.12, P<0.16), whereas direct siRNA effects significantly correlate between the two RNAi screens (ρ=0.42, P<10⁻⁶). (D) Comparison of siRNA effect correlations (upward bars) and their significance (downward bars) between virus infections performed in different cell lines (see also Supplementary Figure 17). Correlations between direct siRNA effects are indicated in blue; Correlations between total siRNA effects are depicted in orange. Dashed black line in lower bar graph indicates the P<0.001 significance threshold. Additionally, absolute numbers of shared hits (up and down) are shown for total (before arrow) and direct (after arrow) effects. (E, F) Absolute increases of siRNA consistency (for calculation, see Supplementary information) are shown for the large-scale RNAi screens of YFV and VACV infection (E) and of interphase nuclei size in HeLa CNX and HeLa MZ cells (F). X-axes indicate the siRNA rank threshold at which the perturbation is considered significant. The number of genes for which 2/3 or 3/3 siRNAs give consistent phenotypes are shown, minus the number of corresponding cases expected by random (the random model is depicted on the right of F), and absolute numbers are indicated for certain examples. Numbers of genes with increased consistency for 2/3 or 3/3 siRNAs, after accounting for population context-propagated effects, are color coded as red and green, respectively (see legend). See Supplementary Figure 31 for the results of all seven large-scale infection screens.

Figure 5

Accounting for cell population context leads to changes in the enrichments of functional annotation classes of hits in SV40 infection. Network view of functional annotation classes (based on DAVID (Dennis Jr et al, 2003) annotation clustering), which become depleted (blue) or enriched (green) after accounting for population context-propagated effects. Values between brackets behind annotation class names indicate the number of genes in these classes before (r) and after accounting for cell population context (c). Classes are linked in the network if they share more than 90% of the genes of the smaller class, for either total or direct effects. Two groups of classes, related to membrane trafficking and to kinase-mediated signalling, are highlighted with yellow-filled ellipses. See Supplementary Figure 28 for the complete overview.

Figure 6

RNAi perturbations of population context-determined infection patterns. (A) Single-cell bootstrapped correlations between virus infection and location on cell islet edges are shown for the 147 siRNAs in each of the 34 small-scale screens. Bootstrapped correlations were calculated combining single-cell measurements from triplicate experiments. The grey region around zero indicates weak correlations. Dot colour and dot size depict the median direct effects of siRNAs on infection and the population size, respectively. Insets on the right show the observed patterns for SFV infection in HeLa Kyoto cells, in which TRIO (upper), or ABL1 (lower), is silenced (Blue=DAPI, green=SFV infection) (see also Supplementary Figure 25). Additional RNAi phenotypes of STK40 and QSK, which have multiple siRNAs (number in parentheses) targeting the same gene that cluster together in the dot plots (clustering indicated with black bold lines), are highlighted. (B) Measured indices of SFV infection in HeLa Kyoto cells located or not on cell islet edges. Median and s.d. of triplicate measurements are shown for TRIO-silenced cells (red), ABL1-silenced cells (blue) and scrambled (non-targeting) siRNA-silenced cells (grey). (C) Protein levels of Trio (red) are higher in sparsely growing cells, whereas protein levels of ABL1 (blue) are higher in cells in densely growing cells. Solid lines and dots indicate median values of mean intensity per cell, for varying local cell densities. Light-coloured regions indicate the inter-quartile range of intensities of individual cells. (D) Silencing of TRIO increases CME (as measured by transferrin uptake relative to that of non-perturbed cells) more strongly in single cells at lower local cell densities (red), whereas ABL1 silencing increases CME more strongly in single cells at high local cell densities. The solid lines represent median values, and transparent region 0.5 × inter-quartile ranges of single-cell values.

Figure 7

Towards a systems-level analysis of host factors regulating infection across mammalian viruses. (A) Bootstrapped hierarchical clustering of all 147 direct siRNA effects of the 34 small-scale RNAi screens. Per branch-point, three bootstrapped empirical P-values were calculated: Red, approximately unbiased P-value; Green: sub-tree P-value; Blue: leaf-set (i.e. clade or cluster) P-value. The tree with the highest summed edge P-value from 10⁶ bootstrapped trees is shown (see Supplementary information). (B) Direct effect of DYRK3 (orange) and FRAP1 (blue) silencing for each assay. Per screen, the median value of all three siRNAs, and median of triplicates is shown (i.e. similar to a two out of three criterion). Asterisks indicate gene-silencing phenotypes, which were independently validated (see Supplementary Figure 32). Distribution below the x-axis depicts typical data density of the normal distribution. Note that for visual purposes the x-axis is inverted, with negative values (reduced infection phenotypes) on the right side. The correlation of DYRK3/FRAP1 phenotypes over all assays is −0.48 (P<0.0038). (C) Gene scores for selected genes are shown for all small-scale RNAi screens of virus infection. Bar colours correspond to different genes (see legend). (D) Functional annotation enrichment P-values are visualized for six functional annotation categories on the clustering of the 34 small-scale RNAi screens. Enrichment was calculated using a rank-based Kolmogorov–Smirnov method (see Supplementary information). Node colour indicates enrichment in down-hits (red) or up-hits (green), and node size indicates P-value of enrichment (see legend).