{'fraction_non_alphanumeric': 0.03769314593869625, 'fraction_numerical': 0.02120763363044511, 'mean_word_length': 4.881676253081348, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 11, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Cardelino: Integrating whole exomes and single-cell transcriptomes to reveal phenotypic impact of somatic variants

Decoding the clonal substructures of somatic tissues sheds light on cell growth, development and differentiation in health, ageing and disease. DNA-sequencing, either using bulk or using single-cell assays, has enabled the reconstruction of clonal trees from frequency and co-occurrence patterns of somatic variants. However, approaches to systematically characterize phenotypic and functional variations between individual clones are not established. Here we present cardelino (https://github.com/PMBio/cardelino), a computational method for inferring the clone of origin of individual cells that have been assayed using single-cell RNA-seq (scRNA-seq). After validating our model using simulations, we apply cardelino to matched scRNA-seq and exome sequencing data from 32 human dermal fibroblast lines, identifying hundreds of differentially expressed genes between cells from different somatic clones. These genes are frequently enriched for cell cycle and proliferation pathways, indicating a key role for cell division genes in non-neutral somatic evolution. Key findings A novel approach for integrating DNA-seq and single-cell RNA-seq data to reconstruct clonal substructure for single-cell transcriptomes. Evidence for non-neutral evolution of clonal populations in human fibroblasts. Proliferation and cell cycle pathways are commonly distorted in mutated clonal populations.


Introduction
Ageing, environment and genetic factors can impact mutational processes, thereby shaping the acquisition of somatic mutations across the life span (Burnet, 1974;Martincorena and Campbell, 2015;Stransky et al. , 2011;Hodis et al. , 2012;Huang et al. , 2018) . The maintenance and evolution of somatic mutations in different sub-populations of cells can result in clonal structure, both within healthy and disease tissues. Targeted, whole-genome and whole-exome DNA sequencing of bulk cell populations has been utilized to reconstruct the mutational processes that underlie somatic mutagenesis (Nik-Zainal et al. , 2012;Alexandrov et al. , 2013;Forbes et al. , 2017;Bailey et al. , 2018;Ding et al. , 2018) as well as clonal trees (Roth et al. , 2014;Deshwar et al. , 2015;Jiang et al. , 2016) .
Availability of single-cell DNA sequencing methods (scDNA-seq; Navin et al. , 2011;Wang et al. , 2014;Navin, 2015) combined with new computational approaches have helped to improve the reconstruction of clonal populations (Kim and Simon, 2014;Navin and Chen, 2016;Jahn et al. , 2016;Kuipers et al. , 2017;Roth et al. , 2016;Salehi et al. , 2017;Malikic et al. , 2017) . However, the functional differences between clones and their molecular phenotypes remain largely unknown. Systematic characterisation of the phenotypic properties of clones could reveal mechanisms underpinning healthy tissue growth and the transition from normal to malignant behaviour.
An important step towards such functional insights would be access to genome-wide expression profiles of individual clones, yielding genotype-phenotype connections for clonal architectures in tissues. Recent studies have explored mapping scRNA-seq profiles to clones with distinct copy number states in cancer, thus providing a first glimpse at clone-to-clone gene expression differences in disease (Müller et al. , 2016;Tirosh et al. , 2016;Fan et al. , 2018) . However, generally-applicable methods for inferring the clone of origin of single cells to study genotype-transcriptome relationships are not yet established.
To address this, we have developed cardelino: a computational method that exploits variant information in scRNA-seq reads to map cells to their clone of origin. After validating our model using simulations, we demonstrate that cardelino allows for accurate assignment of full-length single-cell transcriptomes to the clonal substructure in 32 normal dermal fibroblast lines. With linked somatic variants, clone and gene expression information, we investigate gene expression differences between clones at the level of individual genes and in pathways, which provides new insights into the dynamics of clones. These findings also extend recent studies using bulk DNA-seq data, predominantly in epithelial cells, that have revealed oncogenic mutations and evidence of selective clonal dynamics in normal tissue samples (Behjati et al. , 2014;Martincorena et al. , 2015;Simons, 2016b;Martincorena et al. , 2016;Simons, 2016a) . Our approach can be applied to a broad range of somatic substructure analyses in population or disease settings to reveal previously inaccessible differences in molecular phenotypes between cells from the same individual.

Mapping single-cell transcriptomes to somatic clones with cardelino
We present a two-step strategy for integrating somatic clonal substructure and transcriptional heterogeneity within a population of cells. First, a clonal tree is inferred using variant allele frequencies from bulk or single-cell DNA sequencing data (Roth et al. , 2014;Deshwar et al. , 2015;Jiang et al. , 2016;Kim and Simon, 2014;Navin and Chen, 2016;Jahn et al. , 2016;Kuipers et al. , 2017;Roth et al. , 2016;Salehi et al. , 2017;Malikic et al. , 2017) . Subsequently, cardelino performs probabilistic assignment of individual single-cell transcriptomes to inferred clones, using variant information extracted from single-cell RNA-seq reads ( Fig. 1a ).
Briefly, cardelino is based on a Bayesian mixture model, where individual mixture components correspond to somatic clones in a given population of cells. To enable accurate cell assignment, we employ a beta-binomial error model that accounts for dropout and allelic imbalance, parameters that affect the sensitivity for detecting true somatic variants from sparse scRNA-seq data. Initially, we assess the accuracy of our approach using simulated data that mimic typical clonal structures and properties of scRNA-seq as observed in real data (4 clones, 10 variants per branch, 20% of variants with read coverage, 200 cells, 20 repeat experiments; Methods ). Cardelino achieves high overall performance (Precision-Recall AUC=0.947; Fig. 1b ); for example, at a threshold of P=0.5 for the cell assignment confidence (posterior probability of cell assignment), the model assigns 79% of all cells with 89.0% accuracy.
We explore the effect of key dataset characteristics on cell assignment, including the number of mutations per clonal branch ( Fig. 1c ), the expected number of variants with non-zero scRNA-seq coverage per cell ( Fig. 1d ), the number of clones and the overall extent and variability of allelic imbalance across genes ( Methods; Supp. Fig. S1 ). We also assess individual modelling assumptions in cardelino, confirming the expected benefits from a full Bayesian treatment and the explicit error model for allelic dropout as used in cardelino ( Fig. S1 ).
Taken together, these results indicate that cardelino is broadly applicable to robustly assign individual single-cell transcriptomes to clones, thereby reconstructing clone-specific transcriptome profiles.

Cardelino assigns single cell transcriptomes to clones in human dermal fibroblasts
Next, we apply cardelino to 32 human dermal fibroblast lines derived from healthy donors that are part of the UK human induced pluripotent stem cell initiative (HipSci; Kilpinen et al. , 2017). For each line, we generated deep whole exome sequencing data (WES; median read coverage: 254), and matched Smart-seq2 scRNA-seq profiles using pools of three lines in each processing batch ( Methods ). We assayed between 30 and 107 cells per line (median 61 cells after QC; median coverage: 484k reads; median genes observed: 11,108; Supp. Table S2 ).
Initially, we consider high-confidence somatic single nucleotide variants (SNVs) identified based on WES data ( Methods ) to explore the mutational landscape across lines. This reveals considerable variation in the total number of somatic SNVs, with 41-612 variants per line ( Fig. 2a ). The majority of SNVs can be attributed to the well-documented UV signature, COSMIC Signature 7 (primarily C to T mutations; Forbes et al. , 2017) , agreeing with expected mutational patterns from UV exposure of skin tissues ( Fig. 2a; Supp. Fig. S2; Methods ).
To understand whether the somatic SNVs confer any selective advantage in skin fibroblasts, we used SubConalSelection to identify neutral and selective dynamics at a per-donor level (Williams et al. , 2018) . Other established methods such as dN/dS (Martincorena et al. , 2018) and alternative methods using the SNV frequency distribution (Simons, 2016a;Williams et al. , 2016) are not conclusive in the context of this dataset due to lack of statistical power resulting from the low number of mutations Methods ). In other words, a third of the samples from this cohort of healthy donors contain clones evolving adaptively, which we can investigate in more detail in terms of transcriptome phenotype.
Next, we reconstruct the clonal trees in each line using WES-derived estimates of the variant allele frequency of somatic variants that are also covered by scRNA-seq reads ( Methods ). Canopy (Jiang et al. , 2016) identifies two to four clones per line ( Fig. 2a ). Following this, we use cardelino to map scRNA-seq profiles from a total of 1,728 cells to clones from the corresponding lines ( Methods ), with an assignment rate of 33-100% (at posterior probability P>0.5; corresponding to 4 to 107 cells; median 49.5 cells).
To assess the confidence of these cell assignments, we consider simulated cells drawn from a clonal structure that matches the corresponding line, observing high concordance ( R 2 = 0.92) between the expected and empirical cell-assignment rates ( Fig. 2b ). Lines with clones that harbour fewer distinguishing variants are associated with lower assignment rates ( Supp. Fig. S4 ), at consistently high cell assignment accuracy (median 0.976, mean 0.950; Supp. Fig. S5 ), indicating that the posterior probability of assignment is calibrated across different settings. We also consider the impact of technical features of scRNA-seq data on cell assignment, finding no evidence of biased cell assignments ( Supp. Fig. 6-10 ). Finally, clone prevalences estimated from Canopy and the fractions of cells assigned to the corresponding clones are markedly concordant (overall R 2 = 0.61), providing additional confidence in the cardelino cell assignments ( Fig. 2c ).

Differences in gene expression between clones suggest phenotypic impact of somatic mutations
Initially, we focus on the fibroblast line with the largest number of somatic SNVs ( joxm; white female aged 45-49; Fig. 2a ), with 612 somatic SNVs (112 detected both in WES and scRNA-seq) and 79 QC-passing cells, 97% of which could be assigned to one of three clones ( Fig. 3a ). Principal component analysis of the scRNA-seq profiles of these cells reveals global transcriptome substructure that is aligned with the somatic clonal structure in this population of cells ( Fig. 3b ). Additionally, we observe differences in the fraction of cells in different cell cycle stages, where clone1 has the fewest cells in G1, and the largest fraction in S and G2M ( Fig. 3b, inset plot; PC1 in Supp. Fig. 12-13 ). This suggests that clone1 is proliferating most rapidly. Next, we consider differential expression analysis of individual genes between the two largest clones ( Taken together, the results suggest that somatic substructure in this cell population results in clones that exhibit measurably different expression phenotypes across the transcriptome, with significant differential expression in cell cycle and growth pathways.

Cell cycle and proliferation pathways frequently vary between clones
To quantify the overall effect of somatic substructure on gene expression variation across the entire dataset, we fit a linear mixed model to individual genes ( Methods ), partitioning gene expression variation into a line (likely donor) component, a clone component, technical batch ( i.e . processing plate), cellular detection rate (proportion of genes with non-zero expression per cell) and residual noise. As expected, the line component typically explains a substantially larger fraction of the expression variance than clone (median 5.5% for line, 0.5% for clone), but there are 194 genes with a substantial clone component (>5% variance explained by clone; Fig. 4a ). Even larger clone effects are observed when estimating the clone component in each line separately, which identifies between 331 and 2,162 genes with a substantial clone component (>5% variance explained by clone; median 825 genes; Fig. 4b) . This indicates that there are line-specific differences in the set of genes that vary with clonal structure.
Next, we carry out a systematic differential expression (DE) analysis to assess transcriptomics differences between any pair of clones for each line ( . We also identify a small number of genes that contain somatic variants in a subset of clones, resulting in differential expression between wild-type and mutated clones ( Supp. Fig. S14 ).
To investigate the transcriptomic changes between cells in more detail, we use gene set enrichment

Discussion
Here, we develop and apply a computational approach for integrating somatic clonal structure with single-cell RNA-seq data This allows us to identify molecular signatures that differ between clonal cell populations. Our approach is based on first inferring clonal structure in a population of cells using WES data, followed by the assignment of individual single-cell transcriptomes using a computational approach called cardelino. Our method enables the efficient reconstruction of clone-specific transcriptome profiles from high-throughput assays. Our integrative analysis of bulk WES and scRNA-seq from 32 human fibroblast cell lines reveals substantial phenotypic effects of somatic variation, including in healthy tissue.
Central to our approach is cardelino, a robust model for the probabilistic assignment of cells to clones based on variants contained in scRNA-seq reads. Our approach is conceptually related to de-multiplexing approaches for single-cell transcriptomes from multiple genetically distinct individuals (Kang et al. , 2018) . However, cardelino addresses a substantially more challenging problem: to distinguish cells from the same individual based on the typically small number of somatic variants ( e.g. dozens) that segregate between clones in a population of cells.
Harnessing transcriptomic phenotypic information for cells assigned to clones in fibroblast lines, we identify substantial and convergent gene expression differences between clones across lines, which are notably enriched for pathways related to proliferation and the cell cycle. Analysis of clonal evolutionary dynamics using somatic variant allele frequency distributions from WES data reveals strong evidence for positive selection of clones in ten of 32 lines. These results support previous observations of clonal populations undergoing positive selection in normal human eyelid epidermis assayed by targeted DNA sequencing Simons, 2016b;Martincorena et al. , 2016;Simons, 2016a) . We shed light on the phenotypic origin of this adaptive evolution, as we identify expression of gene sets implicated in proliferation and cancer such as the E2F and MYC pathways.
This surprising result in healthy tissue, suggests pervasive inter-clonal phenotypic variation with important functional consequences.
While we apply cardelino to clonal trees from bulk WES data in this study, our method is general and can be used to assign cells to clonal trees inferred from either bulk or single-cell DNA-seq data. The methods presented here can be applied to any system in which somatic variants tag clonal populations of cells and can be accessed with scRNA-seq assays. Thus, cardelino will be applicable to high-resolution studies of clonal gene expression in both healthy and malignant cell populations as well as in vitro models.
Assignment of cells to clones relies on coverage of somatic variants in scRNA-seq reads, so cell populations with relatively fewer somatic variants may require full-length transcriptome sequencing at higher coverage per cell to enable confident assignments. Increasing both the number of genetically distinct individuals and the numbers of cells assayed per donor would improve the resolution at which differentially expressed genes and pathways can be identified. Larger sample sizes would also motivate further advances in the downstream analysis, for example to account for uncertainty in clonal tree inference and cell-clone assignment when assessing differentially expressed genes. Larger cell counts in particular would also enable joint modeling approaches combining DNA-seq and scRNA-seq data to infer clonal trees and simultaneously assign cells to clones. Our inference methods in cardelino are computationally efficient, so will comfortably scale to multi-site samples and many thousands of cells.
Taken together, our results highlight the utility of cardelino to study gene expression variability in clonal cell populations and suggest that even in nominally healthy human fibroblast cell lines there are clonal populations with growth advantages, opening new avenues to study cell behaviour in clonal populations.

The cardelino model
The input to cardelino is a clonal tree inferred using methods such as Canopy, yielding clonal fractions Fundamentally, we model the alternate read count using a binomial model. For a given site in a given cell, there are two possibilities: the variant is "absent" in the clone the cell is assigned to or the variant is "present" , as encoded in the configuration matrix C . Thus, when considering the "success probability" for the binomial model, where success is defined as observing an alternate read in the scRNA-seq reads, we consider two (sets of) parameters: 0 for homozygous reference alleles (variant absent), and 1 for heterozygous variants (variant present).
The prior probability that cell j belongs to clone k could be assigned as the clonal fraction f k , but to avoid biasing cell assignment towards highly prevalent clones for cells with little read information (where an informative prior is particularly influential) we use a uniform prior such that P(I j = k) = 1 /K for all k . Given this prior distribution, the posterior probability of cell j belonging to clone k can be expressed as: where I j is the identity of the specific clone cell j is assigned to, and a j and d j are the observed alternate and total read count vectors, respectively, for variants 1 to N in cell j . The parameter vector is a set of the unknown binomial success parameters of binomial distributions for modelling the allelic read counts as described above. Specifically, 0 denotes the binomial success rate for the alternative allele when c i,k = 0 (variant absent), thereby accounting for sequencing errors or errors in the clonal tree configuration, and 1 ={ 1 , 2 ..., N } denotes a vector of binomial parameters, one for each variant, for c i,k = 1 . The latter binomial rates model the effect of allelic imbalance, which means the probability of observing alternate reads at frequencies that differ from 0.5 for true heterozygous sites (see Supp. Methods for details). The likelihood for cell j given an assignment to clone k follows then as a product of binomial distributions as follows, .
To capture the uncertainty in the binomial success probabilities, we introduce beta prior distributions on 0 and 1 . To ensure sensible prior distributions, we estimate the beta parameters from the scRNA-seq at known germline heterozygous variants for highly expressed genes ( Supp. Fig. S18 ).
For example, in the fibroblast dataset considered here, this approach yielded prior parameters of beta(0.3, 29.7) for 0 and beta(2.25, 2.65) for i , i >0 .
The posterior probability of cell assignment and the full posterior distribution of these parameters are inferred by a Gibbs sampler, and the details of the algorithm can be found in Supp. Methods , where we also present alternative simpler inference schemes using a constant Binomial success rates or an EM algorithm (see also Supp. Fig. S1 ). Despite the full Bayesian approach, cardelino is computationally efficient, enabling the assignment of hundreds of cells within minutes using a single compute node. These methods will comfortably scale to datasets with many thousands of cells.  (Picelli et al. , 2014) , and library preparation was performed using an Illumina Nextera kit. Samples were sequenced using paired-end 75bp reads on an Illumina HiSeq 2500 machine.

Bulk whole-exome sequencing data and somatic variant calling
We obtained bulk whole-exome sequencing data from HipSci fibroblast (median read coverage: 254) and derived iPS cell lines (median read coverage: 79) released by the HipSci project (Streeter et al. , 2016;Kilpinen et al. , 2017) . Sequenced reads were aligned to the GRCh37 build of the human reference genome (Church et al. , 2011) using bwa (Li, 2013) . To identify single-variant somatic mutation sites in the fibroblast lines, for each exome sample, we searched for sites with a non-reference base in the read pileup using bcftools/mpileup (Li et al. , 2009) . In the initial pre-filtering we retained sites with a per-sample coverage of at least 20 reads, at least three alternate reads in either fibroblast or iPS samples and an allele frequency less than 5% in the ExAC browser (Karczewski et al. , 2017) and 1000 Genomes data (The 1000 Genomes Project Consortium, 2015) . A Fisher exact test (Fisher, 1922) implemented in bcftools/ad-bias was then used to identify sites with significantly different variant allele frequency (VAF) in the exome data between fibroblast and iPS samples for a given line (Benjamini-Hochberg FDR < 10%). Sites were removed if any of the following conditions held: VAF < 1% or VAF > 45% in high-coverage fibroblast exome data; fewer than two reads supporting the alternative allele in the fibroblast sample; VAF > 80% in iPSC data (to filter out potential homozygous alternative mutations); neither the iPS VAF or fibroblast VAF was below 45% (to filter out variants with a "significant" difference in VAF but are more likely to be germline than somatic variants). We further filtered sites to require uniqueness of sites across donors as it is highly unlikely to observe the same point mutation in more than one individual, so such sites almost certainly represent technical artefacts.

Identification of mutational signatures
Signature exposures were estimated using the sigfit package (Gori and Baez-Ortega, 2018) , providing the COSMIC 30 signatures as reference (Forbes et al. , 2017) , and with a highest posterior density (HPD) threshold of 0.9. Signatures were determined to be significant when the HPD did not overlap zero. Two signatures (7 and 11) were significant in two or more donors.

Identification of selection dynamics
Several methods have been developed to detect deviations from neutral growth in cell populations (Simons, 2016a;Williams et al. , 2016Williams et al. , , 2018Martincorena et al. , 2018) . Methods such as dN/dS or models assessing the fit of neutral models to the data need a high number of mutations to determine selection/neutrality. Given the relatively low number of mutations found in the donors in this study, these models are not applicable. We used the package SubClonalSelection (https://github.com/marcjwilliams1/SubClonalSelection.jl) in Julia 0.6.2 which works with a low number of mutations (> 100 mutations; Williams et al. , 2018) . The package simulates the fit of a neutral and a selection model to the allele frequency distribution, and returns a probability for the selection model to fit the data best.
At small allele frequencies the resolution of the allele frequency distribution is limited by the sequencing depth. We chose a conservative lower resolution limit of (Shin et al. , 2017) . At .05 f min = 0 the upper end of the allele frequency distribution we chose a cut-off at to account for ploidy .45 ). For the classification of the donors, we introduced cut-offs on the resulting selection probability = 2 of the algorithm. Donors with a selection probability below 0.3 are classified as 'neutral', above 0.7 as 'selected'. Donors which are neither 'selected' nor 'neutral' remain 'undetermined'. See figure 2A and S3 for the results of the classification and fit of the models to the data.

Single-cell gene expression quantification and quality control
Raw scRNA-seq data in CRAM format was converted to FASTQ format with samtools (v1.5), before reads were adapter-and quality-trimmed with TrimGalore! (github.com/FelixKrueger/TrimGalore) (Martin, 2011) . We quantified transcript-level expression using Ensembl v75 transcripts (Flicek et al. , 2014) by supplying trimmed reads to Salmon v0.8.2 and using the "--seqBias", "--gcBias" and "VBOpt" options (Patro et al. , 2017) . Transcript-level expression values were summarised at gene level (estimated counts) and quality control of scRNA-seq data was done with the scater package (McCarthy et al. , 2017) and normalisation with the scran package (Lun, Bach, et al. , 2016;Lun, McCarthy, et al. , 2016) . Cells were retained for downstream analyses if they had at least 50,000 counts from endogenous genes, at least 5,000 genes with non-zero expression, less than 90% of counts from the 100 most-expressed genes in the cell, less than 20% of counts from ERCC spike-in sequences and a Salmon mapping rate of at least 40% ( Supp . Table S2 ). This filtering approach retains 63.7% of assayed cells.

Deconvolution of donors from pools
To increase experimental throughput in processing cells from multiple distinct donor individuals ( i.e. lines), and to ensure an experimental design robust to batch effects, we pooled cells from three donors in each processing batch, as described above. As such, we do not know the donor identity of each cell at the time of sequencing and cell-donor identity must be inferred computationally. Thus, for both donor and, later, clone identity inference it is necessary to obtain the count of reads supporting the reference and alternative allele at informative germline and somatic variant sites. Trimmed FASTQ reads (described above) were aligned to the GRCh37 p13 genome with ERCC spike-in sequences with STAR in basic two-pass mode (Dobin et al. , 2012) using the GENCODE v19 annotation with ERCC spike-in sequences (Searle et al. , 2010) . We further use picard (Broad Institute, 2015) and GATK version 3.8 (McKenna et al. , 2010) to mark duplicate reads ( MarkDuplicates ), split cigar reads ( SplitNCigarReads ), realign indels ( IndelRealigner ), and recalibrate base scores ( BaseRecalibrator ).
For cell-donor assignment we used the GATK HaplotypeCaller to call variants from the processed single-cell BAM files at 304,405 biallelic SNP sites from dbSNP (Sherry et al. , 2001) build 138 that are genotyped on the Illumina HumanCoreExome-12 chip, have MAF > 0.01, Hardy-Weinberg equilibrium P < 1e-03 and overlap protein-coding regions of the 1,000 most highly expressed genes in HipSci iPS cells (as determined from HipSci bulk RNA-seq data). We merged the per-cell VCF output from GATK HaplotypeCaller across all cells using bcftools version 1.7 (Danecek et al. , 2011(Danecek et al. , , 2016 and filtered variants to retain those with MAF > 0.01, quality score > 20 and read coverage in at least 3% of cells.
We further filtered the variants to retain only those that featured in the set of variants in the high-quality, imputed, phased HipSci genotypes and filtered the HipSci donor genotype file to include the same set of variants.
We used the donor_id function in the cardelino package to assign cells to donors. This function assigns cells to donors by modelling alternative allele read counts with given genotypes of input donors. For a single germline variant, the three base genotypes (as minor allele counts) can be 0, 1 and 2. For doublet genotype profiles generated by combining pairs of donor genotypes, two additional combinatory genotypes, 0.5 and 1.5 are allowed. We assume that each genotype has a unique binomial distribution whose parameters are estimated by an EM algorithm in a framework similar to clone assignment (described above). When we enable doublet detection, the posterior probabilities probability. If the provisionally assigned donor is one of the three donors known to have been pooled together for the specific plate, then we deem the cell to be confidently assigned to that donor, otherwise we deem the cell to have "unassigned" donor. With this approach, 97.4% of cells passing QC (see above) are confidently assigned to a donor ( Supp. Fig. S19 ). Of the cells that are not confidently assigned to a donor, 2.1% are identified as doublets by cardelino and 0.5% remain "unassigned" due to low variant coverage or low posterior probability . Thus, we have 2,338 QC-passing, donor-assigned cells for clonal analysis.

Clonal inference
We inferred the clonal structure of the fibroblast cell population for each of the 32 lines (donors) using Canopy (Jiang et al. , 2016) . We used read counts for the variant allele and total read counts at filtered somatic mutation sites from high-coverage whole-exome sequencing data from the fibroblast samples as input to Canopy. In addition to the variant filtering described above, input sites were further filtered for tree inference to those that had non-zero read coverage in at least one cell assigned to the corresponding line. We used the BIC model selection method in Canopy to choose the optimal number of clones per donor. Here, for each of the 32 lines, we considered the highest-likelihood clonal tree produced by Canopy, along with the estimated prevalence of each clone and the set of somatic variants tagging each clone as the given clonal tree for cell-clone assignment.

Cell-clone assignment
For cell-clone assignment we required read the read counts supporting reference and alternative alleles at somatic variant sites. We used the bcftools version 1.7 mpileup and call methods to call variants at somatic variant sites derived from bulk whole-exome data, as described above, for all confidently assigned cells for each given line. Variant sites were filtered to retain variants with more than three reads observed across all cells for the line and quality greater than 20. We retained cells with at least two somatic variants with non-zero read coverage (2,044 cells across 32 lines). From the filtered VCF output of bcftools we obtained the number of reads supporting the alternative allele and the total read coverage for each somatic variant site with more than three reads covering the site, in total, across all the donor's cells. In general, read coverage of somatic mutation sites in scRNA-seq is sparse, with over 80% of sites for a given cell have no overlapping reads. We used the scRNA-seq read counts at the line's somatic variant sites to assign QC-passing cells from the line to clones using the clone_id function in the cardelino R package.

Simulations to benchmark cell to clone assignment
We simulated data to test the performance of cardelino as follows. First, given a clonal tree configuration C ( N -by-K binary matrix), a given number of cells are generated (e.g. 500, see below), whose genotypes are sampled from K clones following a multinomial distribution parameterised by

Differential expression, pathway and variance component analysis
Expression analyses between clones required further filtering of cells for each line. Analyses were conducted using cells that passed the following filtering procedure for each line: (1) clones identified in the line were retained if at least three cells were confidently assigned to the clone; (2) cells were retained if they were confidently assigned to a retained clone. Lines were retained for DE testing if they had at least 15 cells assigned to retained clones, allowing us to conduct expression analyses for 31 out of the 32 lines (all except vils ).
Differential gene expression (DE) testing was conducted using the quasi-likelihood F-test method (Lund et al. , 2012) in the edgeR package (Robinson et al. , 2010;McCarthy et al. , 2012) as recommended by Soneson and Robinson (Soneson and Robinson, 2018) . To test for differences in expression between cells assigned to different clones in a donor, we fit a linear model for single-cell gene expression with cellular detection rate (proportion of QC-filtered genes expressed in a cell; numeric value), plate on which the cell was processed (a factor) and assigned clone (a factor) as predictor variables. The quasi-likelihood F test was used to identify genes with: (1) any difference in average expression level between clones (analogous to analysis of variance), and (2) differences in average expression between all pairs of clones ("pairwise contrasts"). We considered 10,876 genes that were sufficiently expressed (an average count >1 across cells in all donors) to test for differential expression.
To test for significance of overlap of DE genes across donors, we sampled sets of genes without replacement the same size as the number of DE genes (FDR < 10%) for each line. For each permutation set, we then computed the number of sampled genes shared between between donors.
We repeated this procedure 1,000 times to obtain distributions for the number of DE genes shared by multiple donors if shared genes were obtained purely by chance.
Gene set enrichment (pathway) analyses were conducted using the camera (Wu and Smyth, 2012) method in the limma package (Smyth, 2004;Ritchie et al. , 2015) . Using log 2 -fold-change test statistics for 10,876 genes for pairwise contrasts between clones from the edgeR models above as input, we applied camera to test for enrichment for the 50 Hallmark gene sets from MSigDB, the Molecular Signatures Database (Liberzon et al. , 2011) . For all differential expression and pathway analyses we adjusted for multiple testing by estimating the false discovery rate (FDR) using independent hypothesis weighting (Ignatiadis et al. , 2016) , as implemented in the IHW package, with average gene expression supplied as the independent covariate.
Expression variance across cells is decomposed into multiple components in a linear mixed model, including cellular detection rate (proportion of genes with non-zero expression per cell) as a fixed effect and plate (i.e. experimental batch), donor (i.e. line; only when combining cells across all donors) and clone (nested within donor for combined-donor analysis) as random effects. We fit the linear mixed model on a per-gene basis using the variancePartition R package (Hoffman and Schadt, 2016) .

Code availability
The cardelino methods are implemented in an open-source, publicly available R package ( github.com/PMBio/cardelino ). The code used to process and analyse the data is available (github.com/davismcc/fibroblast-clonality), with a reproducible workflow implemented in Snakemake (Köster and Rahmann, 2012) . Descriptions of how to reproduce the data processing and analysis workflows, with html output showing code and all figures presented in this paper, are available at davismcc.github.io/fibroblast-clonality. Docker images providing the computing environment and software used for data processing ( hub.docker.com/r/davismcc/fibroblast-clonality/ ) and data analyses in R (hub.docker.com/r/davismcc/r-singlecell-img/) are publicly available.

Data availability
Single-cell RNA-seq data have been deposited in the ArrayExpress database at EMBL-EBI