The nature and nurture of cell heterogeneity: accounting for macrophage gene-environment interactions with single-cell RNA-Seq

Identifying different macrophage states (phenotypes)

Exploratory RNA-seq data analysis (Fig. 1b) revealed a clear differentiation between one- and eight-hour cultured cells. It also highlighted a subgroup of cells present at both time points that enriched for a gene set previously demonstrated to be differentially expressed in macrophages treated with the (anti-)inflammatory TGFB cytokine (Additional file 1). This provided initial evidence for the potential biological relevance of the observed heterogeneity among the macrophages. To more formally study this latent/emergent (i.e. previously unknown) heterogeneity and the gene expression differences underlying it, we developed a novel hybrid model-driven and non-parametric clustering method constrained to report only cell clusters that were well represented across replicates. We provide further details in the Methods, with all analysis code and results provided in the Additional files 1 and 2 respectively. Our more formal clustering analysis confirmed the existence of the lower-abundance cell state at both time points (called ‘cluster one’ in Fig. 1c&d) but also a reproducible third state (called ‘cluster three’) emerging after eight hours. A search for associations between proportions of the cell states and culture conditions revealed a lower abundance of cluster three when in standard media with LPS (Fisher’s Exact -log₁₀P =5.14, Fig. 1e&f), suggesting another inflammatory sub-phenotype in the cells. With these results demonstrating heterogeneity that would otherwise be missed with traditional sequencing, we then asked if any of these cell types resembled the tissue-resident phenotype of interest.

Identifying the macrophage phenotype of interest

To better understand these latent cell states, we used our modification on the Heskes Rank Product method [15] to estimate upper and lower p-value bounds on differential and heterogeneous gene expression. As detailed in the Methods, we defined differential expression as a global/overall shift in gene expression, while heterogeneous/context-specific expression was defined as highly variable expression across cell states or culture microenvironments. We selected this non-parametric gene ranking approach not only for its focus on result reproducibility, but for its speed and ease of data fusion that make it well suited to single-cell analysis.

Compared against one-hour cultures, most cells in eight-hour cultures (cluster two cells) progressed towards an anti-inflammatory transcriptional signature consistent with the tissue-resident macrophage phenotype of interest. This is evidenced, for example, by the expression of IL1RN, encoding IL-1 receptor antagonist, which increased 18-fold over time in cluster two (Fig. 2a). Recombinant IL1RN is used to treat severe inflammatory conditions mediated by its ligand, the archetypal pro-inflammatory cytokine IL-1 [16]. In contrast to other cells, cluster one cells did not move towards a transcriptional profile consistent with a tissue-resident phenotype, and tended to maintain up- or down-regulation of cluster-specific genes. This pattern included lower expression of GAPDH and TPI1, genes involved in glycolysis, a pathway known to vary across macrophage activation phenotypes [17]. Compared with the other clusters, cluster one expressed higher levels of FOXP1 across all culture conditions (Fig. 2d, −log₁₀P of 6.03 and 4.85 at one and eight hours, both globally significant at a 5% FDR). FOXP1 is a transcription factor involved in maintaining embryonic stem cell pluripotency that must be turned off for complete monocyte differentiation into macrophages [18]. Concerned that selective expression of FOXP1 might represent technical differentiation heterogeneity, we used the imaging data to search for other contributing factors. While we did not find evidence for cell motility or morphology associations, we noted that cluster one cells were more likely to have come into contact with the culture chamber retention beads used to prevent cells from escaping (Fig. 1a, Fisher’s Exact -log₁₀P = 4.10). We also noted that cells cultured towards the edges of the chips were more likely to touch retention beads (Fisher’s Exact -log₁₀P = 1.79), but that the chip edge positions were enriched for in cluster one, seemingly independent of the bead association (Fisher’s Exact -log₁₀P = 3.19, Additional file 1: Figure S6). As cell imaging was performed hourly, this may in part reflect false negatives for detecting bead contact of cells in a less mature and adherent state. Notwithstanding the evidence for both differentiation and environmental factors underlying cluster one, the low correlation between cluster one and off-chip bulk tissue samples (Additional file 1: Figure S7) led us to conclude that it may represent a potentially interesting but low-abundance phenotype that is unrelated to our main question of SAMHD1 biology in tissue-resident macrophages.

A comparison of gene expression in cluster three versus other cells (Fig. 2b) was consistent with a shift towards macrophages with a tissue remodelling phenotype. This was evidenced, for example, by greater expression of TIMP3, an inhibitor of extracellular matrix degradation, and CHI3L1, a secreted glycoprotein also thought to drive tissue remodelling and a known genetic risk factor for asthma severity [19]. Since cluster three’s transcriptional phenotype most closely resembles that of cluster two (Additional file 1: Figure S18), it seems plausible that this represents an emerging secondary phenotype that requires macrophage signaling when cells are exposed to LPS (Fig. 1e&f). Future studies over longer culture periods will help clarify the microenvironmental dependence and persistence of this cluster.

We focus the remainder of our observations on the main body of cells (cluster two), resembling the desired tissue-resident phenotype. While no cluster two enrichment for particular motility or morphology imaging features were found over these short culture periods, we briefly note an imaging subtype of cells observed to be phagocytosing the retention beads. Phagocytosing macrophages did not assume a distinct transcriptomic cluster 一 likely due to short culture times 一 nevertheless they did enrich for oxidative stress and mitochondrial genes, examples being PRDX1 and MT-CO3 (−log₁₀P of 4.60 and 4.43 respectively, both globally significant at 5% FDR). This unexpected result suggests an avenue for single-cell studies to explore the temporal dynamics of phagocytosis [20].

Changes in macrophage behavior (cluster two) with SAMHD1 knockout

After filtering our data to focus on the cell subtype of interest (cluster two), we tested for varying knockout and wild-type differential and heterogeneous expression over time. The most striking feature of the globally significant knockout effects was that they were predominantly microenvironment specific (Fig. 3a). Not only does this stress the importance of studying gene-environment interactions in cellular genetics models with known phenotypic plasticity, but in this work allows us to comment on macrophage signaling contributions. Specifically, we note a highly significant association with SOD1 expression, a gene that encodes the cytosolic isozyme of superoxide dismutase. SOD1 is a superoxide radical scavenger that may confer some protection against HIV-1 neuropathy as part of the oxidative stress response [21], but which has no described associations with HIV-1 restriction factors. As with SAMHD1 loss of function, gain-of-function mutations in SOD1 are associated with neuroinflammation and degeneration, clinically manifesting as amyotrophic lateral sclerosis, which is likely a result of toxic protein aggregates [22]. Overall, SOD1 transcript levels in this study were similar in wild-type and SAMHD1 knockout cells. However, whereas SOD1 expression in wild-type cells stimulated with LPS increased over time as expected, expression in knockout cells increased over time in the absence of LPS activation. Not only does this point to a qualitative difference in expression but 一 as can be seen from the barely detectable expression changes in standard media (Fig. 3a) 一 the magnitude of the difference is much greater in conditioned media. This SOD1 oxidative stress response thus has both a notable genetic and signaling component, which may be dampened by identifying and blocking the augmenting signaling factor(s).

Further exploring this oxidative stress response, we noticed a change in SOD1’s coexpression (i.e. genes correlating with it) in knockouts versus wild-types (Fig. 3c). Analysing which genes tend to coexpress with SOD1 allows for biological contextualisation of this effect, as genes are only expressed when required, and so expression correlation implies coregulation. A key difference here being the increased coexpression resolution offered by single-cell analysis to generate gene networks in both perturbed and unperturbed cells, compared with conventional sequencing that is limited to only studying average gene correlations with perturbation. We asked whether there are genes with which SOD1 alters its rank (Spearman) correlation across cells when SAMHD1 is knocked out. Fig. 3b lists the globally significant signatures from MSigDB v5.1 [23], consisting of gene sets more correlated with SOD1 in knockout versus wild-type cells. In other words, these highlight genes and pathways coregulating with the altered oxidative stress response in the knockout macrophages. In addition to the expected enrichment for reactive oxygen species (ROS) modulators in this list, we note a striking enrichment for proteasome genes. The proteosome is known to be activated and dysregulated with SOD1 mutation, leading to aberrant cellular proteostasis [24]. However, it is plausible in these cells that the altered proteostasis is simply a result of the macrophages’ stress response. In particular we note an enrichment for proteasome genes involved with the activation of NF-kappaB (NF-kB), a pro-inflammatory transcription factor on which multiple macrophage signaling pathways converge [25]. A similar proteostasis theme is observed for the regulation of ornithine decarboxylase, a macrophage anti-inflammatory enzyme [26] that metabolises a well described marker of macrophage pro- versus anti-inflammatory phenotype polarisation. With these pathway associations validating this approach by supporting known biology of increased inflammation with SAMHD1 knockout - albeit within the context of oxidative stress and altered proteostasis - we asked if other pathways in Fig. 3b point to unexpected biological themes.

A striking pathway result was the repeated enrichment for proteasome genes involved in the degradation of cell cycle proteins. While the direct causal relationship between cell cycle and SAMHD1 is considered to be via its cyclin-dependent kinase phosphorylation and inhibition [27], these results point to an additional reverse relationship. A tumor suppressor role for SAMHD1 has been proposed based on its maintenance of genome integrity and cancer associations with downregulation [12, 13]. However, the pathways in Fig. 3b counter-intuitively point to p53 dependent and independent G₁ arrest as part of the DNA damage response. These results highlight a necessary fine balance in SAMHD1 activity in terms of cell cycle control, with too little SAMHD1 arresting cell cycle progression via other tumour suppressors in response to DNA damage. This relationship between SAMHD1 and G₁ arrest has been demonstrated in dividing fibroblasts that turn senescent [28], though an important difference here is that it suggests the wild-type G₀ (post-mitotic) macrophages shift to a G₁ block with SAMHD1 knockout. Finally, and perhaps most interesting in terms of the HIV-1 host-pathogen relationship, is the highly significant pathway enrichment for the degradation of APOBEC3G as part of the altered stress response in SAMHD1 knockouts. This HIV-1 restriction factor causes proviral DNA hypermutation via cytidine deamination, and so shares no known overlap in antiviral activity with SAMHD1 [29].

Stem cells, generation of macrophages and experimental media

The human embryonic stem cell line HUES-2 was obtained from the HUES Facility, University of Harvard [36]. Feeder-free PSC cells were cultured in mTeSR^TM-1 medium (Stem Cell Technologies) on Matrigel (Corning)-coated tissue culture dishes, passaged with TrypLE (Invitrogen) with the addition of 10 μmol/L Rho-kinase inhibitor Y-27632 (Abcam). A double-nicking CRISPR-Cas9 approach was used to generate SAMHD1-knockout stem cell lines [37]. Plasmid pX462 (gift from Feng Zhang; Addgene plasmids cat. 48141, [38]), expressing the guide RNA, D10A-mutated Cas9 and a puromycin-selection cassette was adapted to target SAMHD1 at exon 4 (GTGTATCAATGATTCGGACGAGG and CGATACATCAAACAGCTGGGAGG; PAM underlined) or exon 5 target sites (CGTTCACTTATCTGCAGCTCTGG and GGATGTCTAGTTCACGCACTGGG; PAM underlined) using protocols previously described [38, 39]. PSCs were transfected with all four plasmids targeting SAMHD1 using the Neon® Transfection system (Invitrogen) according to manufacturer’s guidelines (2 × 10⁶ cells electroporated with 15 μg DNA using a 100 μL tip at 1000 V, 40 ms pulse width, one pulse), cultured without antibiotics for 48 h and then for 48 h with selection in 0.4 μg/mL puromycin (Sigma). Single-cell clones were generated by plating transfected PSCs at low density onto mitotically-inactivated mouse embryonic feeder (MEF) cells [40, 41] on gelatin-coated tissue culture plates in stem cell medium (KO-DMEM, 2 mmol/L L-Glutamine, 100 mmol/L nonessential amino acids, 20% serum replacement, and 8 ng/mL FGF2; Invitrogen). Clones with modifications at the SAMHD1 locus were identified by high resolution melt analysis. Sequencing confirmed that clone E2 had an out-of-frame insertion (29 and 49 bp) into each allele of SAMHD1 exon 4 and G9 had an out-of-frame deletion (43 bp) in one allele and an in-frame deletion (39 bp) that deleted the essential allosteric GTP binding site (amino acids 135 to 147) and would alter the catalytic site. Macrophages derived from these clones were screened for SAMHD1 expression by western blotting using a mouse anti-SAMHD1 antibody (clone 2D7, Insight Biotechnology Ltd) and a rabbit anti-GAPDH antibody (Sigma) (Additional file 1: Figure S1).

Feeder-free PSC cells were cultured in mTeSR^TM-1 medium (Stem Cell Technologies) on Matrigel (Corning)-coated tissue culture dishes, passaged with TrypLE (Invitrogen) with the addition of 10 μmol/L Rho-kinase inhibitor Y-27632 (Abcam). A protocol devised in our laboratory was used to generate macrophages from PSC cultures. Briefly, embryoid bodies were formed using the spin method in AggreWells™800 (Stemcell Technologies) plates, each of which was split into two monocyte factories in T175 tissue culture flasks containing ~150 embryoid bodies. The monocytes released into the supernatant were harvested regularly and plated into 96-well plates at 5 × 10⁴ cells per well in macrophage differentiation medium consisting of XVIVO™15 (Lonza) supplemented with 100 ng/mL M-CSF (Invitrogen), 2 mM glutamax (Invitrogen), 100 U/mL penicillin and 100 μg/mL streptomycin (Invitrogen). Four days after plating the media was replaced with fresh macrophage differentiation media with additional 10% fetal bovine serum (FBS; Invitrogen) and the cells were used on day 7 of differentiation. We used three batches (A, B, G, each corresponding to one AggreWells™800 plate) of wild-type cells and two batches (G9: C, D and E2: E, F) each of the two SAMHD1-knockout clones for the Polaris macrophage stimulation experiments.

Four kinds of media (‘conditioned’ and ‘standard’, with or without LPS) were used in the microfluidic chips. To generate conditioned media, stem cell-derived monocytes were plated at 5 × 10⁵ cells/well in a 12-well tissue culture plate and differentiated for 4 days in macrophage differentiation medium followed by 3 days in the presence of 10% FBS. The medium was then replaced with fresh macrophage medium plus 10% FBS with or without 100 ng/mL LPS (Sigma). “Standard media” were generated by incubation of macrophage differentiation medium plus 10% FBS either with or without 100 ng/ml LPS at 37 °C to simulate incubation of conditioned media with cells. For each kind of medium, after 24 h with LPS/mock stimulation, supernatants were recovered by centrifugation, 0.45 μm-filtered and stored at −80 °C then thawed and clarified by centrifugation at 14,100 rcf for 10 min before use.

Polaris™ protocols

Polaris runs followed the protocol ‘Using Polaris to Generate Single-Cell cDNA Libraries for mRNA Sequencing’ (PN 101–0082 A1, Fluidigm). Set-up parameters are shown in Table S1. Polaris integrated fluidic circuits (IFCs) were prepared for cell capture by a priming step, during which the capture chambers were also coated with the extracellular matrix compound fibronectin (25 ng/μl, cat. F4759, Sigma-Aldrich) for handling adherent cells, and capture beads (prepared to Fluidigm specifications) were loaded to prevent the release of captured single cells (Fig. 1). Priming was arranged to finish as the cell mix became ready for loading (Additional file 1: Figure S2A).

For each experimental run, one well of cells from a single batch each of wild-type and knockout cells was used. One sample was stained with CellTracker™ Orange CMRA Dye and the other with both CellTracker™ Orange CMRA and CellTracker™ Green CMFDA Dye (cat. C34551 and C7025, Thermo Fisher Scientific) in Wash Buffer (Fluidigm). Dye concentrations were adjusted upwards between runs during the experiment to improve sensitivity and resolution: 1 μM for runs 1 and 2; 2 μM for runs 3 to 12; and 3 μM for runs 13 to 20. Wild-type and knockout cells were dual- or single-stained in approximately equal numbers of experiments (Additional file 1: Figure SB). The culture medium was removed and the cells were washed twice with 250 μl Wash Buffer (Fluidigm) then stained with 100 μl of staining solution at 37 °C for 15 min before a further wash with 150 μl of Wash Buffer. The Wash Buffer was removed, 100 μl of 0.5 mM EDTA (cat. 15575, Gibco) in PBS was added and the cells were incubated at 37 °C for 15 min before the EDTA was removed and the cells were resuspended gently in 50 μl of Feed Media (X-VIVO™ 15 cat. BE02-061Q (Lonza), 1% Penicillin-Streptomycin (10,000 U/mL), cat. 15140–122, (Gibco)). Cells were counted on a TC-20TM Automated Cell Counter (Bio-Rad) and samples were adjusted to a final concentration of 350–400 cells/μl.

Fluidigm guidelines (Fluidigm Single-Cell Preparation Guide, PN 100–7697) were used to establish optimal buoyancy at an 4:1 (cells:cell suspension reagent) ratio for the Polaris experiment. For each run, 15 μl each of differentially stained wild type and knockout cells were mixed with 7.5 μl of Cell Suspension Reagent (Fluidigm, PN 101–0434) 25 μl of the resulting suspension, containing an estimated 7000–8000 cells, was loaded on the Polaris IFC. On-board imaging settings to control automatic cell selection were set for each run, with a threshold that varied between 4000 and 6000 for a constant 1.0 s exposure. After completion of the cell selection step, the IFC was removed from the Polaris system and placed in an incubator at 37 °C and 5% CO₂ for 2 h to allow the cells to settle prior to dosing. Dosing culture media were centrifuged in small volumes (500 μl) at 14,100 rcf for 10 min before 27 μl of medium and 200 μl of Feed Media were loaded into the appropriate wells of the IFC (Additional file 1: Figure S2C). For 1-h dosing runs, the dosing step was stopped manually and for 8-h dosing runs, cells were also dosed at 4 h. Successive runs were scheduled to alternate between one and eight hour dosing.

The Polaris acquired images of all 48 chambers in all available fluorescence channels during cell capture, at the start of dosing and every 1 h thereafter until the end of the run. IFCs were removed from the Polaris for additional, high-resolution imaging on a Leica TCS SP8 confocal microscope (Leica Microsystems) at 3 stages of each run: after cell capture, after the 2-h incubation and at the end of the dosing step, just before lysis. (Additional file 1: Figure S3A). The imaging protocol was designed to acquire bright-field and fluorescence images, for CellTracker Orange and CellTracker Green: 488 nm and 561 nm excitation and hybrid detectors for emission at 500–550 nm and 571–630 nm respectively, were used with a 20× objective and a template to automatically locate the 48 IFC cell isolation chambers, in an image acquisition process lasting ~5 min during which the microscope chamber temperature was maintained at 37 °C. The black vinyl film on the lower surface of the Polaris chip was removed temporarily for each imaging stage. For each time-point and channel, the Polaris stored a single image of the IFC’s 48 cell chambers, 5200 × 1000 pixels in size, in which each pixel had dimensions of approximately 5.5 × 5.5 μm. On the Leica, one image of 512 × 512 pixels, each 1.39 × 1.39 pixels, was captured per cell chamber in bright field, orange and green channels. The initial aims of the image analysis were to differentiate double- and single-stained cells (wild-type and knockout or vice-versa, depending on the run) and to assess cell shape, motility and phagocytosis behaviour. The R package EBImage [42] was used with a set of bespoke R functions to automatically identify individual cell chambers and the cells they carried. Formally, within the large Polaris image(s), cells were detected as clusters of 4 or more of the brightest 30 pixels in the trimmed image corresponding to each cell chamber, for each fluorescent channel. Motility was assessed by measuring changes in the cluster positions between imaging time-points. To assess cell morphology (circularity, the proportion of pixels within the smallest circle encompassing the cluster that belong to the cluster), the higher-resolution Leica images were analysed using a re-scaled version of the same algorithm, requiring clusters of 60 or more of the brightest 411 pixels. Automated imaging assessments were checked extensively by eye, leading to the identification of irregularly shaped or unevenly stained cells. Macrophage phagocytosis of capture beads was assessed manually.

For each Polaris run, a bulk control sample of the wild-type and of the knockout cells used in the run was prepared for sequencing. Aliquots containing 2000–4000 of the stained cells loaded on the Polaris were kept at room temperature, incubated along with the loaded IFC for the 2-h incubation step, and kept at 37 °C during dosing. Within 30 min of starting the Polaris lysis step, bulk-cell samples were lysed and processed using the RNA extraction kit RNeasy Micro Kit (cat# 74004, Qiagen) and eluted in 14 μl of RNAse-free water. At the end of the dosing protocol and after any additional imaging, cells were lysed for reverse transcription and amplification for cDNA generation in the Polaris. Master mixes for this procedure were prepared using the SMARTer Ultra Low RNA Kit (Clontech), according to the Fluidigm Polaris protocol with minor modifications. The cell lysis mix (28 μl) contained 8.0 μl Polaris Lysis Reagent, 9.6 μl of a 1/200 dilution of Polaris Lysis Plus Reagent (in PCR-grade water, prepared immediately before use), 9.0 μl SMARTer kit 3′ SMART CDS Primer IIA and (for runs 1–11) 1.4 μl of diluted ERCC spike-in RNA, prepared by adding 1.0 μl of a 1/10 dilution of ERCC ExFold RNA Spike-In Mixes (cat. 4456739, Ambion) to 96.5 μl Polaris Loading Reagent and 2.5 μl SMARTer Kit RNase Inhibitor (40 U/μl). The reverse transcription reaction mix (48 μl) contained 15.5 μl 5× First-Strand Buffer (RNase-free), 1.9 μl DTT, 7.7 μl dNTP mix, 7.7 μl IIA Oligonucleotide, 1.9 μl RNAse Inhibitor, 7.7 μl SMARTscribe Reverse Transcriptase (all SMARTer Kit, Clontech), 2.4 μl Polaris Loading Reagent and 3.2 μl Polaris RT Plus Reagent. The PCR mix (90 μl) contained 63.5 μl PCR-grade water, 10.0 μl 10× Advantage 2 PCR Buffer (not SA – Short Amplicon), 4.0 μl 50× dNTP Mix, 4.0 μl 50× Advantage 2 Polymerase Mix (all Advantage 2 Kit, Clontech) and 4.0 μl SMARTer IS PCR primer. The cDNA products were harvested into a 96-well plate in the arrangement shown in Additional file 1: Figure S2D. Harvest wells with atypical volumes (some with no material, others with an excess) were excluded from further analysis. Bulk control samples comprising 1 μl of RNA extracted from the bulk samples above were each mixed with 4.5 μl of cell lysis mix and processed using the following temperature sequence: 37 °C for 5 min, 72 °C for 3 min, 25 °C for 1 min, then 4 °C (hold), then reverse transcribed using 9.0 μl of the reaction mix, at 42 °C for 90 min followed by enzyme inactivation at 70 °C for 10 min, then 4 °C (hold). Bulk-sample PCRs contained 1 μl of cDNA generated in the previous step and 9.3 μl PCR mix. PCR conditions were as follows: 95 °C for 1 min then 5 cycles of {95 °C/20 s, 58 °C/4 min, 68 °C/6 min}, 9 cycles of {95 °C/20 s, 64 °C/30 s, 68 °C/6 min}, 7 cycles of {95 °C/30 s, 64 °C/30 s, 68 °C/7 min}, then 72 °C for 10 min and 4 °C (hold). After Quant-iT™ PicoGreen (Thermo Fisher) normalization of cDNAs to 0.22 ng/μl, Illumina sequencing libraries were prepared together using Nextera XT DNA Library Preparation Kit (Illumina), according to the manufacturer’s specifications, at quarter-scale on a Beckman-Coulter FXp automated liquid handling instrument, and pooled as up to 192 multiplexed libraries using in-house dual-indexing library tags. Some samples that had failed because of empty cell chambers or the presence of more than one cell in a chamber were included, and two bulk-cell samples per run were included in each pool, which was sequenced as a single lane of 75 b paired-end reads on an Illumina HiSeq 4000 instrument in 4 different runs.

RNA-Seq data generation and initial quality control

Samples were prepared and paired-end sequenced using the Illumina HiSeq™4000 Sequencing platform, as described above. RNA-seq reads were trimmed for Nextera/Illumina adapter sequences using skewer-v0.1.125 [43]. Trimmed reads were mapped to a modified reference genome comprising the human genome, Homo sapiens GRCh37, and fasta sequences for ERCC spike-ins (ThermoFisher). Reads in gzipped fastq format were aligned using Hisat2 version-2.0.0-beta [44] with default parameters. Duplicate reads were marked using MarkDuplicates.jar implemented in Picard tools v1.92. BAM alignments were name sorted with Samtools version 1.1. Alignment metrics were calculated using CollectRnaSeqMetrics.jar implemented in Picard tools v1.92 for full BAM files and with potential PCR duplicates marked. RNA-SeQC [45] was used to calculate sequencing bias, as the median estimators in Picard can result in zero estimates. Reads mapping uniquely to genes annotated in ENSEMBL release 76 were counted using featureCounts [46] implemented in subread-v1.5.0 [47]. Read distribution between various features - assigned reads (mapped uniquely to exons), multiple mapping, ambiguous mapping, No features (mapped uniquely to intronic and intergenic regions) – was obtained from featureCounts results. Read counts were normalized to Transcripts per million (TPM), and number of detected genes per sample were calculated by counting genes with at least 1 TPM. Details of the supplied data and meta-data provided are provided in Additional file 1, including further QC using the R scater package [48]. In brief, we excluded culture chambers with visually confirmed doublets (two cells), numbers of detected genes more similar to bulk controls, and cells with very low starting cDNA. FASTQ files for the data have been uploaded to the Gene Expression Omnibus (GSE87849).

Clustering the macrophages and gene expression modeling

All modeling and analyses were performed in the R environment (version 3.2.1, x86_64-pc-linux-gnu 64-bit, Ubuntu 14.04.2 LTS). Multidimensional scaling (MDS) on the cell rank correlations - without gene weighting - was used to reduce one hour and eight hour cells to five dimensions, as detailed in the provided Additional file 1 code wrapped in the cellStates() function. Inclusion of a greater number of dimensions did not influence the clustering. Additional file 1: Figure S16 & S17 demonstrate the cell densities and individual cells for the five dimensions in the eight hour cells, where some clusters of cells can be seen to be unique to individual chips. To ensure reproducibility, the constraint was added that a cluster with fewer than three cells in more than half of the replicates not be considered in the downstream modelling. While such groups of cells may represent illuminating features in the macrophage model’s behaviour, these are not reproducible effects.

Extracting major cell states was performed using a hybrid of model-driven clustering (Gaussian mixtures) and non-parametric clustering (partitioning around medoids), as detailed in the Additional file 1 code wrapped in the consensusCluster() function. Repeated fitting of mixtures to the data - with one chip replicate omitted per iteration - was used to produce a consensus matrix of the proportion of iterations cells share clusters. This was non-parametrically re-clustered, selecting the maximum number of clusters for which experimental replication was strongly represented in each cluster. We note that although this combination of methods was employed, other tested approaches such as hierarchically clustering the consensus matrix produced similar results due to the consensus and reproducibility constraints. The defined clusters are recorded in the cell meta-data (Additional file 2). The major latent cell states identified with exploratory data analysis remain the only two identified major states in the one hour cells, with a third reproducible cluster emerging at eight hours that shares properties with the main group of cells (Additional file 1: Figure S18).

Altered gene expression was modelled as the change in mean conditioned on (i.e. tested per) cell subtype of interest. Overall mean across cell subtypes (μ) was used as a measure of global shift in gene expression, while mean absolute deviation (MAD) of the subtypes was used as a measure of variability and so context/subtype specificity. For example, for d_i = change in mean expression over time for cell subtype i, μ = (d₁ + d₂ + .. + d_n)/n and MAD = (|d₁-μ| + |d₂-μ| + .. + |d_n-μ|)/n. More robust estimators, such as the use of quantiles, provided similar top hits, so here we present μ and MAD estimates, focusing rather on a rank product framework to determine statistical significance of rank reproducibility per sequencing library [15]. A focus on gene ranks has several advantages well suited to single-cell work: it is non-parametric, robust, comments directly on reproducibility, and allows data fusion or meta-analysis without the need for complex data normalisations. For computational speed, the Heskes rank product algorithm was used to assess bounds on the statistical significance. Reported p-values are the geometric means of the upper and lower bounds provided by the rankprodbounds() function, with the qvalue package used to estimate global significance at 5% false discovery rate [49]. Rank (Spearman) correlations were used to estimate gene co-expression, followed by testing for co-expression signature enrichment with Preranked Gene Set Enrichment Analysis (default settings) [23].