{"paper_id":"0f7e683a-8bd8-49bb-91e3-9ffe64a01b07","body_text":"1 \nEnsemblex: a n accuracy -weighted ensemble genetic demultiplexing  framework for 1 \npopulation-scale scRNAseq sample pooling 2 \n 3 \nAuthors: Michael R. Fiorini1,2, michael.fiorini@mail.mcgill.ca; Saeid Amiri2 4 \nsaeid.amiri@mcgill.ca;  Allison A. Dilliott2,3, allison.dilliot@mcgill.ca;  Cristine M. Yde Ohki4, 5 \ncristinemarie.ydeohki@uzh.ch; Lukasz Smigielski4, lukasz.smigielski@kjpd.uzh.ch; Susanne 6 \nWalitza4,5,6, susanne.walitza@pukzh.ch;  Edward A. Fon2,3, ted.fon@mcgill.ca;  Edna 7 \nGrünblatt4,5,6, edna.gruenblatt@kjpd.uzh.ch;  Rhalena A. Thomas2,3*, rhalena.thomas@mcgill.ca;  8 \nSali M.K. Farhan1,2,3*, sali.farhan@mcgill.ca  9 \n 10 \n1. Department of Human Genetics, McGill University, Montreal, Quebec H3A 2B4, Canada  11 \n2. The Montreal Neurological Institute-Hospital, McGill University, Montreal, Quebec H3A 2B4, 12 \nCanada 13 \n3. Department of Neurology and Neurosurgery, McGill University, Montreal, Quebec H3A 2B4, 14 \nCanada 15 \n4. Department of Child and Adolescent Psychiatry and Psychotherapy , Psychiatric University 16 \nHospital Zurich, University of Zurich, Zurich, Switzerland 17 \n5. Neuroscience Center Zurich, University of Zurich and the ETH Zurich, Zurich, Switzerland 18 \n6. Zurich Center for Integrative Human Physiology, University of Zurich, Switzerland 19 \n* Shared co-senior authorship 20 \nCorresponding authors:  Sali M.K. Farhan, Sali.farhan@mcgill.ca; Rhalena A. Thomas, 21 \nrhalena.thomas@mcgill.ca  22 \n  23 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 2 \nAbstract  24 \nMultiplexing samples from distinct individuals  prior to sequencing is a promising step toward  25 \nachieving population-scale single-cell RNA sequencing by reducing the restrictive costs of the 26 \ntechnology. Individual genetic demultiplexing tools resolve the donor-of-origin identity of pooled 27 \ncells using natural genetic variation  but present diminished accuracy on highly multiplexed  28 \nexperiments, impeding the analytic potential of the dataset. In response, we introduce Ensemblex: 29 \nan accuracy-weighted, ensemble genetic demultiplexing framework that integrates four distinct 30 \nalgorithms to identify the most probable subject labels. Using computationally and experimentally 31 \npooled samples, we demonstrate Ensemblex’s superior accuracy and illustrate the implications of 32 \nrobust demultiplexing on biological analyses. 33 \n 34 \nKeywords: single-cell RNA sequencing, multiplexing, sample pooling, genetic demultiplexing, 35 \ninduced pluripotent stem cells , differential gene expression, dopaminergic neurons, doublet 36 \ndetection, accuracy-weighted probability, high-throughput sequencing 37 \n 38 \n 39 \n 40 \n 41 \n 42 \n 43 \n 44 \n 45 \n 46 \n 47 \n 48 \n 49 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 3 \nBackground 50 \nSingle-cell RNA sequencing (scRNAseq) continues to revolutionize our molecular understanding 51 \nof biology by providing unprecedented insight into the transcriptional landscape of individual 52 \ncells. Unlike bulk RNAseq, where the RNA from all cells within a tissue is sequenced to produce 53 \ntotal expressional profiles across all cells, scRNAseq captures transcriptional signatures at a single-54 \ncell resolution, elucidating the diverse gene expression across distinct cell types  and subtypes. 55 \nDifferential gene expression (DGE) can then be calculated between subgroups of cells to reveal 56 \ncell type-specific expression changes between patient or treatment  groups. However, scRNAseq 57 \nhas come at the expense of increased costs, hindering its application for population-scale analyses, 58 \nwhich are critical for  deriving clinico-pathological associations  and characterizing the genetic 59 \nheterogeneity of complex diseases in biomedical sciences (1, 2). 60 \n 61 \nIn addition to the expense of separately capturing and sequencing cells from individual donors, the 62 \ncosts of scRNAseq are exacerbated for  cell cultures , such as those derived from induced 63 \npluripotent stem cells (iPSC) (1). In particular, neurological diseases are difficult to study in human 64 \ntissue because access to post-mortem brains is limit ed and experimental manipulations are not 65 \npossible; in contrast, iPSC-derived cultures of neurons and other brain cells grown from 66 \nreprogrammed skin or blood cells of human donors are an excellent model of the brain  (3). 67 \nHowever, iPSCs from each donor  must be individually plated  and differentiated in parallel , 68 \npresenting prohibitively high consumable and labour costs that render the methodology unfeasible 69 \nfor population-scale analyses. Multiplexing cultures by pooling cells from multiple donors  prior 70 \nto growth and differentiation , droplet capture , and sequencing,  is one solution to address this 71 \nlimitation as it reduces costs by a factor of the number of samples multiplexed (4). Similarly, 72 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 4 \nsamples such as tumor biopsies can be pooled at acquisition to realize the same benefits. In turn, 73 \ngenetic demultiplexing tools are cost-effective, statistical frameworks that use the natural genetic 74 \nvariation at sites of single-nucleotide polymorphisms (SNP) observed in  the transcriptome to 75 \ncluster cells on the basis of their  donor’s genotype. Importantly, genetic demultiplexing can be 76 \ninformed by prior genotype information of the donors  to improve demultiplexing accuracy and 77 \nfacilitate the assignment of  each cell back to its specific donor -of-origin, which is critical for 78 \ndownstream analyses aiming to investigate discrepancies between subjects. At present, six genetic 79 \ndemultiplexing tools have been developed for scRNAseq: Demuxalot (5) and Demuxlet (6) both 80 \nrequire prior genotype information as input ; Freemuxlet (6) relies entirely on the de novo 81 \ntranscriptome and does  not incorporate prior genotype information; and ScSplit (7), Souporcell 82 \n(8), and Vireo (9) provide versions of the algorithm that can work with and without prior genotype 83 \ninformation (Table 1).  84 \n 85 \nA robust genetic demultiplexing tool is tasked with  mitigating the addition of technical artifacts 86 \ninto scRNAseq datasets by correctly classifying each pooled cell to its donor -of-origin, correctly 87 \nidentifying heterogenic doublets (erroneous barcodes composed of two or more cells from distinct 88 \nsubjects), and quantifying its confidence in the demultiplexed labels so that low -confidence 89 \nclassifications can be eliminated from downstream analyses. While benchmarking analyses on the 90 \navailable genetic demultiplexing tools have shown effectiveness for demultiplexing small sample 91 \nsizes, limitations emerge as the number of multiplexed samples approach a population scale  (6) 92 \n(7) (8) (9). For example, using computationally multiplexed samples, Neavin et al. evaluated the 93 \nperformance of genetic demultiplexing tools as the number of samples approached a population 94 \nscale and observed diminished demultiplexing accuracy with increasing numbers of pooled 95 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 5 \nsamples, as well as notable classification discrepancies between tools  (10). Furthermore, even at 96 \nsmall sample sizes, divergent assignments between genetic demultiplexing tools are common (8) 97 \n(9) (11). Another feature that has been shown to affect genetic demultiplexing performance is the 98 \nunderrepresentation of samples in a pool , which is especially relevant for cell culture -based 99 \nmultiplexed experiments, as variable growth rates in vitro across cell lines is common (12) (8) (9). 100 \nGenetic demultiplexing tools have also shown low concordance for identifying heterogenic 101 \ndoublets, which should be removed prior to downstream analyses to avoid technical noise in the 102 \ndata (10). Importantly, benchmarking analyses have repeatedly highlighted ScSplit’s poor 103 \nperformance relative to the remaining tools (9) (10) (8) (11). The sum of these limitations calls to 104 \nquestion the robustness of the individual genetic demultiplexing tools for resolving the donor 105 \nidentities of highly multiplexed samples, which represents  an important hurdle for feasibly 106 \nachieving population-scale scRNAseq analysis. 107 \n 108 \nIn response to the divergent assignments commonly observed across tools, a consensus framework, 109 \nwhereby only cells that show matching sample labels across all individual tools are retained for 110 \ndownstream analyses, may appear sufficient to resolve the risk of introducing technical noise into 111 \nthe data from misclassified cells. However, consensus frameworks are restricted to performing  112 \nonly as well as the worst -performing tool , and genetic demultiplexing performance is highly 113 \ndataset dependent (10); thus, the overall performance of a consensus framework can vary 114 \nimmensely between datasets. To this end, Neavin et al. proposed a majority vote framework for 115 \ngenetic demultiplexing, whereby a cell is assigned to the sample called by the majority of tools 116 \n(10). However, this approach can be vulnerable to a subset of tools performing poorly on the 117 \ndataset, does not allocate additional weight to the votes of tools that perform more favourably on 118 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 6 \nthe dataset, cannot account for instances when ties occur amongst tools , and cannot capture cells 119 \nthat are correctly classified by only one tool. The sum of these limitations leads to the unnecessary 120 \nremoval of cells from downstream analyses, reducing statistical power, especially for  highly 121 \nmultiplexed pools where each donor, on average, will have a lower representation of cells in the 122 \npool. Moreover, the ability to capture the transcriptional profiles of rare cell types with scRNAseq 123 \nprovides a notable advancement over bulk RNAseq and can strongly influence biological 124 \ninterpretations (13); thus, investigators are reluctant to discard valuable cells in order to maximize 125 \nthe analytic potential of their dataset. 126 \n 127 \nTo address the need for a robust genetic demultiplexing framework that can maximize the number 128 \nof confidently classified cells  retained for downstream analyses , achieve high demultiplexing 129 \naccuracy for population-scale scRNAseq sample pooling, and maintain reliability across different 130 \ndatasets, we developed Ensemblex: an accuracy -weighted ensemble genetic demultiplexing 131 \nframework designed to identify the most probable sample labels from each of its constituent tools 132 \n— Demuxalot, Demuxlet/Freemuxlet, Souporcell, and Vireo. Our ensemble method capitalizes on 133 \ncombining distinct statistical frameworks for genetic demultiplexing while adapting to the overall 134 \nperformance of its constituent tools on the respective dataset , making it resilient against a poorly 135 \nperforming tool and facilitating a higher yield of cells for downstream analyses.  The Ensemblex 136 \nworkflow is assembled into a three-step pipeline — 1) accuracy-weighted probabilistic ensemble; 137 \n2) graph -based doublet detection; 3) Ensemble -independent doublet detection  — and can 138 \ndemultiplex pools with or without prior genotype information.  139 \n 140 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 7 \nHere, we showcase Ensemblex’s improved demultiplexing performance across a variety of settings 141 \nthrough benchmarking analyses on a total of 141 computationally multiplexed pools with known 142 \nground-truth sample labels ranging in size from 4 to 80 samples. We applied the ensemble method 143 \nto three diverse, experimentally multiplexed datasets: 1) non -small cell lung cancer (NS CLC) 144 \ndissociated tumor cells from 7 individuals with donor-specific oligonucleotide labels ; 2) iPSC -145 \nderived dopaminergic neurons  (DaN) from 22 healthy individuals; and 3) iPSC -derived neural 146 \nstem cells (NSC) from 9 individuals with attention deficit hyperactivity disorder (ADHD) and 7 147 \nhealthy controls. We demonstrate Ensemblex’s robustness across distinct datasets, its ability to 148 \nreturn a high proportion of confidently classified cells for downstream analysis, and the 149 \nimplications that its improved demultiplexing performance has on biological interpretations of 150 \nmultiplexed experiments.  151 \nTable 1. Summary of individual genetic demultiplexing tools.  152 \nGenetic demultiplexing tool \nPrior genotype information for \ngenetic demultiplexing \nIncluded in the Ensemblex \nframework \nDemuxalot (5) Required Yes \nDemuxlet (6) Required Yes \nFreemuxlet (6) Not supported Yes \nScSplit (7) Optional No \nSouporcell (8) Optional Yes \nVireo (9) Optional Yes \n 153 \n 154 \nResults and Discussion 155 \nEvaluating the performance of existing individual genetic demultiplexing tools 156 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 8 \nTo evaluate the performance of individual genetic demultiplexing tools , we generated 157 \ncomputationally multiplexed pools using scRNAseq of 80 different iPSC lines from Parkinson’s 158 \ndisease patients and healthy controls, which were differentiated towards a DaN state as part of the 159 \nFoundational Data Initiative for Parkinson’s Disease (FOUNDIN-PD) (14). Processed scRNAseq 160 \ndata from the independent iPSC lines were merged to simulate sample-pooling using a previously 161 \ndescribed protocol (9), which provided known ground -truth donor and doublet labels. We 162 \ngenerated 96 in silico pools ranging in size from 4 to 80 multiplexed samples, where each sample 163 \ncorresponded to a unique donor-of-origin. The in silico pools averaged 17,396 cells per pool with 164 \na constant 15% doublet rate.  165 \n 166 \nLeveraging whole-genome sequencing (WGS) of the 80 donors from which the iPSC lines were 167 \nderived and the four genetic demultiplexing tools that can utilize prior genotype information  — 168 \nDemuxalot, Demuxlet, Souporcell, and Vireo -GT — we first investigated the proportion of 169 \ncorrectly classified cells by the individual tools (Figure 1A). Across the 96 in silico pools, all tools 170 \nshowed decreasing demultiplexing performance as  the number of samples within the pool 171 \nincreased. Souporcell demonstrated the largest decrease in the proportion of correctly classified 172 \ncells as the number of multiplexed samples increased from 4 (mean = 90.60%) to 80 (mean =  173 \n53.27%). In accordance with previous findings  (10, 15) , the individual genetic demultiplexing 174 \ntools performed better on singlet classification than doublet detection, highlighting an avenue for 175 \nimproved genetic demultiplexing accuracy by increasing the rate of heterogenic doublet 176 \nidentification (Figure 1A).  177 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 9 \n 178 \nFigure 1. Evaluation of existing individual genetic demultiplexing tools.  Evaluation of genetic 179 \ndemultiplexing tools with prior genotype information  on 96 in silico pools with known ground -180 \ntruth sample labels ranging in size from 4 to 80 multiplexed induced pluripotent stem cell (iPSC) 181 \nlines from genetically distinct individuals, averaging 17,396 cells per pool and a 15% doublet rate. 182 \nA) Line graphs showing the proportion of correctly classified singlets, doublets, and all cells by 183 \neach individual genetic demultiplexing tool across varying numbers of multiplexed iPSC lines in 184 \na single pool (sample number). The large dots show the mean proportion of correct classifications 185 \nby an individual tool across replicates at a given sample size (n = 9 per pool size). The blue points 186 \nshow the proportion of cells that were correctly classified by at least one individual genetic 187 \ndemultiplexing tool: Demuxalot, Demuxlet, Souporcell, or Vireo -GT. B) Bar chart showing the  188 \nmean proportion of total cells from an individual pool  correctly classified by only one genetic 189 \ndemultiplexing tool. Error bars represent one standard deviation from the mean. (n = 9 per pool 190 \nsize) C) Bar chart showing the  proportion of correctly classified singlet cells labelled as 191 \n“unassigned” (ambiguous singlet assignments)  due to assignment probabilities below the 192 \nrecommended threshold of the respective genetic demultiplexing tool. Error bars represent one 193 \nstandard deviation from the mean. (n = 9 per pool size). 194 \n 195 \nWe also investigated the proportion of cells that were correctly classified by at least one  genetic 196 \ndemultiplexing tool to designate the best possible performance of an ensemble method  that 197 \nsuccessfully incorporates every correct classification from its constituent tool s (Figure 1A ). 198 \nAcross the 96 in silico pools, an average of 93.64% of cells were correctly classified by at least 199 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 10 \none tool. In comparison, Demuxlet, which demonstrated the best overall performance  amongst 200 \nindividual tools , correctly classified 86.73% of cells, on average. Demuxalot was consistently 201 \nresponsible for the highest proportion of cells correctly classified by only one tool ; 1.21% of 202 \npooled cells , on average,  were correctly classified by Demuxalot  only, followed by Demuxlet 203 \n(mean =  0.83%), Vireo-GT (mean =  0.29%), and  Souporcell (mean =  0.26%) (Figures 1B; 204 \nAdditional File 1: Figure S1). Conversely, a consensus framework , correctly classified  only 205 \n81.06% of cells, on average (data not shown). Based on these results, we reasoned that an ensemble 206 \ngenetic demultiplexing method that can identify  the most probable sample label  from its 207 \nconstituent tools, independent of a consensus assignment,  would increase the yield of correctly 208 \nclassified cells.  209 \n 210 \nNext, we explored the frequency at which correctly classified singlets were labelled as unassigned 211 \nbecause their assignment probability failed to meet the tool’s recommended probability threshold. 212 \nAcross the 96 in silico pools, Vireo-GT consistently showed the highest proportion of correctly 213 \nclassified singlets with insufficient assignment probabilities  (Vireo-GT mean = 7.86%) followed 214 \nby Demuxalot (mean = 5.91%), Demuxlet (mean = 2.44%) and Souporcell (mean = 2.34%) 215 \n(Figure 1C ). While a stringent probability threshold is important to prevent erroneous 216 \nclassifications in downstream analyses, w e reasoned that the unnecessary removal of correctly 217 \nclassified cells could be mitigated by a carefully calibrated ensemble method that allocates 218 \nadditional assignment confidence to cells with matching sample labels across constituent tools , 219 \ndespite low internal tool-specific assignment probabilities.  220 \n 221 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 11 \nWe repeated the above analyse s u sing the same  96 computationally multiplexed pools and the 222 \ngenetic demultiplexing tools that do not require prior  genotype information : Freemuxlet, 223 \nSouporcell, and Vireo . Here, we observed the same overarching limitations  as when 224 \ndemultiplexing with prior genotype information: 1) decreasing demultiplexing performance as the 225 \nnumber of multiplexed samples increased; 2) poor doublet detection performance compared to 226 \nsinglet classification; 3)  high rates of cells only correctly classified by a single tool ; and 4) 227 \ndiscarded correctly classified cells due to insufficient assignment probabilities (Additional File 1: 228 \nFigure S2). When we compared demultiplexing with and without prior genotype information, we 229 \nobserved a trend  towards a higher proportion of cells being correctly classified  when prior 230 \ngenotype information was available , as previously seen in separate benchmarking analyses (9) 231 \n(Additional File 1: Figure S3).   232 \n 233 \nValidating the Ensemblex framework on pools with known ground-truth sample labels 234 \nTo mitigate the  limitations of the  individual genetic demultiplexing tools  and maximize the 235 \nanalytic potential of multiplexed scRNAseq datasets, we developed Ensemblex (Figure 2A). The 236 \nEnsemblex workflow begins by demultiplexing pooled samples with four distinct demultiplexing 237 \nalgorithms, followed by three steps: 1) accuracy-weighted probabilistic ensemble; 2) graph-based 238 \ndoublet detection ; and 3) ensemble-independent doublet detection  (Figure 2B) . As output, 239 \nEnsemblex returns its own cell-specific sample labels and corresponding assignment probabilities, 240 \nas well as the sample labels and corresponding assignment probabilities for each of its constituent 241 \ntools. 242 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 12 \n 243 \nFigure 2. Characterization of the Ensemblex  framework. Ensemblex is a probabilistic -244 \nweighted ensemble genetic demultiplexing framework for s ingle-cell RNA sequencing analysis, 245 \nwhich was designed to leverage the most probable sample labels from each of its constituent tools: 246 \nDemuxalot, Demuxlet, Souporcell, and Vireo  when using prior genotype information or 247 \nDemuxalot, Freemuxlet, Souporcell, and Vireo when prior genotype information is not available. 248 \nA) The Ensemblex workflow begins with demultiplexing pooled cells from genetically distinct 249 \nindividuals by each of the constituent tools. The outputs from each individual demultiplexing tool 250 \nare then used as input into the Ensemblex framework. B) The Ensemblex framework comprises 251 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 13 \nthree distinct steps that are assembled into a pipeline: 1) accuracy-weighted probabilistic ensemble, 252 \n2) graph -based doublet detection, and 3) ensemble -independent doublet detection. C-D) Line 253 \ngraphs showng t he contribution of each step of the Ensemblex framework on 96 in silico pools 254 \nwith known ground -truth sample labels ranging in size from 4 to 80 multiplexed induced 255 \npluripotent stem cell (iPSC) lines from genetically distinct individuals, averaging 17,396 cells per 256 \npool and a 15% doublet rate. The average proportion of correctly classified C) singlets and D) 257 \ndoublets across replicates at a given pool size is shown after sequentially applying each step of the 258 \nEnsemblex framework with prior genotype information ( n = 9 per pool size) . The right panels 259 \nshow the average proportion of correct classifications across all 96 pools; error bars represent one 260 \nstandard deviation from the mean. The blue points show the proportion of cells that were correctly 261 \nclassified by at least one individual genetic demultiplexing tool: Demuxalot, Demuxlet, 262 \nSouporcell, or Vireo-GT.  263 \n 264 \nIn response to our observation that certain cells are correctly classified by only one tool, we 265 \nimplemented the accuracy-weighted probabilistic ensemble component (Step 1) of the Ensemblex 266 \nframework. In brief, this unsupervised weighting model identifies the most probable sample label 267 \nfor each cell by assigning weights to each tool’s assignment probabilities based on their estimated 268 \nbalanced accuracy for the dataset (see “Methods”) (Figures 2B) (16). Ensemblex then retains the 269 \nsample label with the highest cumulative probability across its constituents.  However, one 270 \nchallenge for this framework is computing the balanced accuracy of the constituent tools for 271 \nexperimentally multiplexed pools that lack ground-truth labels. Therefore, to estimate the balanced 272 \naccuracy of a particular constituent tool (e.g., Demuxalot) without ground-truth labels, Ensemblex 273 \nleverages the cells with a consensus assignment across the three remaining tools ( e.g., Demuxlet, 274 \nSouporcell, and Vireo-GT) as a proxy for ground -truth. To validate this approach, we utilized in 275 \nsilico pools with known ground truth sample labels to compute the Adjusted Rand Index (ARI) 276 \nbetween Ensemblex’s sample labels when the balanced accuracy of the constituent tool s was 277 \ncomputed using consensus labels or ground -truth labels. Here, we consistently observed a mean 278 \nARI > 0.99,  independent of the number of multiplexed samples in a pool,  suggesting high 279 \nassignment concordance between the two approaches  (Additional File 1: Figure  S4). Applying 280 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 14 \nthe accuracy-weighted probabilistic ensemble component to the 96 in silico  pools correctly 281 \nclassified 94.98% of singlets, on average, across all pools, approaching the number of singlets that 282 \nwere correctly classified by at least one constituent tool (mean = 96.48%) (Figure 2C). In contrast, 283 \nonly 66.01% of doublets, on average, were correctly identified across all pools  after Step 1 , 284 \ncompared to 76.59% of doublets that were correctly identified by at least one constituent tool 285 \n(Figure 2D).  286 \n 287 \nGiven that previous analyses have demonstrate d strong doublet call discordance across genetic 288 \ndemultiplexing tools (10), it was unsurprising that Step 1 of the Ensemblex framework performed 289 \npoorly on doublet identification. Therefore, instead of relying on the cell type classifications of the 290 \nconstituent tools (i.e., singlet or doublet), we elected to leverage the doublet-related features (e.g., 291 \ndoublet probability; see “Methods”) returned by the constituent tools to identify the cells with the 292 \nhighest doublet likelihood, independent of the existing classifications.  We implemented this 293 \napproach in the graph-based doublet detection component (Step 2) of the Ensemblex framework , 294 \nwhich was specifically designed to increase the rate of true doublet detection . Step 2 begins by 295 \nidentifying the top n most confident doublets in the pool  (see “Methods”). Then, based on the 296 \nEuclidean distance s in principal component analysis ( PCA) space, the cells that appear most 297 \nfrequently amongst the nearest-neighbors of the high confident doublets and exceed the optimized 298 \npercentile threshold for the nearest-neighbor frequency are labelled as doublets by Ensemblex 299 \n(Figure 2B ; Additional File 1: Figure S5; see “Methods”). Upon applying the graph-based 300 \ndoublet detection  component to the 96 in silico  pools following Step 1 , Ensemblex correctly 301 \nidentified 76.00% of doublets, on average: a 9.99% increase in doublet detection from Step 1. In 302 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 15 \nturn, the average proportion of correctly classified singlets across all pools (94.43%) decreased by 303 \nonly 0.55% (Figure 2D).   304 \n 305 \nThe ensemble-independent doublet detection  component (Step 3) of the Ensemblex framework 306 \nwas implemented to further improve doublet detection. Step 3 was motivated by our observation 307 \nthat certain tools, namely Demuxalot and Vireo, showed high doublet detection specificity (mean 308 \n= 0.99) on in silico pools with known ground-truth sample labels, but that Steps 1 and 2 failed to 309 \nincorporate a subset of these correct doublet calls (Additional File 1: Figure S6). Therefore, by 310 \ndefault, Ensemblex accepts the  doublet calls made by Demuxalot and Vireo -GT (Figure 2B ). 311 \nApplying the ensemble-independent doublet detection component to the 96 in silico  pools 312 \nfollowing Steps 1 and 2 further increased the average proportion of correctly identified doublets 313 \nacross all pools by 1.58% for a total of 77.63% of doublets detected , while only decreasing the 314 \naverage proportion of correctly classified singlets by  0.13% for a total of 94.30% of singlets 315 \ncorrectly classified (Figures 2C and 2D ). Notably, owing to the graph -based doublet detection 316 \ncomponent, the average proportion of doublets identified by Ensemblex exceeded the average 317 \nproportion of doublets that were correctly classified by at least one constituent tool.  318 \n 319 \nWhile the three-step workflow of the Ensemblex pipeline was designed to maximize the balance 320 \nbetween singlet classification and doublet  identification, we do prioritize the identification of 321 \ndoublets at the expense of a slightly lower singlet yield to minimize technical noise in the data.  322 \nHowever, we recognize that different experimental designs will require varying levels of doublet 323 \ndetection stringency ; thus , users can modify the percentile thresholds for graph -based doublet 324 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 16 \ndetection and nominate different tools for ensemble -independent doublet detection  (see 325 \n“Methods”).  326 \n 327 \nBenchmarking Ensemblex on pools with known ground-truth sample labels 328 \nTo benchmark Ensemblex against Demuxalot, Demuxlet, Souporcell, and Vireo-GT with prior 329 \ngenotype information, we first utilized the 96 in silico pools with known ground -truth sample 330 \nlabels to assess how Ensemblex’s demultiplexing performance varied as the number of multiplexed 331 \nsamples approached a cohort scale (4-80 samples). Unlike doublets, singlets were only considered 332 \ncorrectly classified if their assignment probability exceeded the recommended threshold of the 333 \nrespective tool. On average across all pools, Ensemblex showed a higher proportion of correctly 334 \nclassified singlets ( mean = 92.19%), doublets ( mean = 77.63%), and all cells (mean = 90.12%) 335 \nthan the other tools. In comparison,  Demuxlet, widely considered the “gold standard” tool, 336 \ncorrectly classified 89.72% of singlets, 68.57% of doublet s, and 86.73% of all cells, on average 337 \n(Figures 3A -3C). Importantly, the discrepancy in the proportion of correctly classified cells 338 \nbetween Ensemblex and the next-best tool was amplified as the number of multiplexed samples 339 \nincreased from 4 (2.78%) to 80 ( 3.52%), demonstrating that our ensemble method  was able to 340 \npartially mitigate decreased demultiplexing accuracy as the pools approach a population scale. 341 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 17 \n 342 \nFigure 3. Ensemblex ground-truth benchmarking on computationally multiplexed pools. The 343 \ngenetic demultiplexing tools with prior genotype information were evaluated on 96 in silico pools 344 \nwith known ground -truth sample labels ranging in size from 4 to 80 multiplexed induced 345 \npluripotent stem cell (iPSC) lines from genetically distinct individuals, averaging 17,396 cells per 346 \npool and a 15% doublet rate. A singlet was considered correctly classified if the assigned sample 347 \nlabel matched the ground-truth sample label and the assignment probability exceeded the 348 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 18 \nrecommended threshold for the respective tool; a doublet was considered correctly classified if the 349 \nassigned sample label matched the ground -truth sample label, regardless of the assignment 350 \nprobability. A-I) Line graphs showing the performance of Ensemblex and the individual genetic 351 \ndemultiplexing tools across evaluation metrics. The large dots show the mean value for each tool 352 \nacross replicates at a given sample size (n = 9 per pool size). A) Proportion of correctly classified 353 \nsinglets. B) Proportion of correctly classified doublets. C) Proportion of correctly classified cells. 354 \nD) Adjusted Rand Index between each tool’s sample labels and the ground-truth sample labels. E) 355 \nBalanced accuracy of each tool. F) Matthew’s Correlation Coefficient of each tool. G) Area under 356 \nthe receiver operating characteristic curve (AUC) of the singlet assignment probability for each 357 \ntool. H) Proportion of usable cells returned by each tool. Usable cells were defined as  cells 358 \nclassified by singlets with an assignment probability exceeding the recommended threshold of the 359 \nrespective tool. I) Error rate amongst the usable cells returned by each tool; erroneous 360 \nclassifications comprised of true doublets labeled as singlets or true singlets assigned to the wrong 361 \nsample. 362 \n 363 \nNext, we applied evaluation metrics for classification models to gauge the overall performance of 364 \nthe genetic demultiplexing tools. We first computed the ARI to evaluate the similarity between the 365 \ndemultiplexed sample labels and the ground -truth sample labels . Here, Ensemblex showed the 366 \nhighest ARI with the ground truth sample labels  across all pools (mean = 0.76), followed by  367 \nDemuxalot ( mean = 0.67)  and Demuxlet ( mean = 0.66) (Figure 3D ). We then computed the 368 \nbalanced accuracy to evaluate the binary classification performance — singlet or doublet — of 369 \neach genetic demultiplexing tool as well as the Matthew’s Correlation Coefficient (MCC), which 370 \nprevious work has suggest ed is more reliable and informative for classification cases where 371 \npositive (singlet) and negative (doublet) cases have the same analyt ic importance (17). Across all 372 \npools, Ensemblex showed the highest balanced accuracy (mean = 0.80) and MCC (mean = 0.64), 373 \nwhereas Demuxalot and Demuxlet showed average balanced accuracies of 0.74 and 0.75, 374 \nrespectively, and both tools show ed an average MCC of 0.54 (Figures 3E and 3F). To evaluate 375 \nhow well Ensemblex’s confidence score (see “Methods”) and each constituent tool’s assignment 376 \nprobability corresponded to the accuracy of their singlet classification, we plotted the area under 377 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 19 \nthe receiver operating characteristic curve (AUC). Although Demuxalot (mean = 0.99) and Vireo-378 \nGT (mean = 0.99) showed the highest AUC across all pools on average, Ensemblex’s AUC was 379 \ncomparable (mean = 0.98) (Figure 3G).  380 \n 381 \nFinally, we investigated the proportion of usable cells returned by each demultiplexing tool and 382 \nthe error rate amongst usable cells. We define usable cells as singlet classifications exceeding the 383 \nrecommended probability threshold of the respective tool, while the error rate amongst usable cells 384 \nconstituted incorrectly classified singlets to the wrong donor-of-origin or true doublets incorrectly 385 \nclassified as singlets. We observed that, on average, Ensemblex returned the highest proportion of 386 \nusable cells across all pools  (82.66%), followed by Demuxlet (81.66%), Souporcell ( 81.01%), 387 \nDemuxalot (79.99%), and Vireo-GT (77.53%) (Figure 3H). Importantly, Ensemblex showed the 388 \nlowest error  rate amongst usable cells  (4.34%), followed by Demuxalot ( 4.43%), Demuxlet 389 \n(5.77%), Vireo-GT (6.16%), and Souporcell (21.82%) (Figure 3I).  390 \n 391 \nUsing computationally multiplexed pools comprised of 24 iPSC lines, we further assessed how the 392 \nperformance of Ensemblex varied as a function of the number cells in a pool when prior genotype 393 \ninformation was available. Here, we observed that our ensemble method consistently outperformed 394 \nthe individual demultiplexing tools  (Additional File 1: Figure S 7). When cells are pooled 395 \nexperimentally, it is reasonable to expect some iPSC lines to be underrepresented in the pool.  396 \nTherefore, to assess Ensemblex’s demultiplexing performance in the pre sence of an 397 \nunderrepresented iPSC line, we produced computationally multiplexed pools comprising of 24 398 \nsamples, with one sample showing varying degrees of under representation.  Again, we observed 399 \nthat Ensemblex consistently outperformed the individual tools (Additional File 1: Figure S8). 400 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 20 \nFinally, we repeated the above analyses to assess whether the benefits of using Ensemblex to 401 \ndemultiplex with prior genotype information extended to cases where prior genotype information 402 \nis not available . In doing so , we observed a trend towards better overall performance by 403 \nEnsemblex; however, the discrepancy between Ensemblex and the top -performing individual 404 \ntools, namely Freemuxlet and Souporcell, was less pronounced than when demultiplexing with 405 \nprior genotype information (Additional File 1: Figures S9-S11).  406 \n 407 \nTaken together, these results indicate that the Ensemblex framework mitigates the limitations of 408 \nthe individual tools, leading to greater overall demultiplexing performance across computationally 409 \nmultiplexed pools with known ground-truth labels . Ultimately, Ensemblex’s improved 410 \ndemultiplexing performance translates to a higher recovery of usable cells for downstream 411 \nanalyses as well as a higher accuracy amongst usable cells, limiting the unnecessary removal of 412 \ncells from the dataset and mitigating the introduction of technical artifacts into biological analyses. 413 \n 414 \nEvaluating Ensemblex on experimentally pooled samples with donor-specific oligonucleotide 415 \nlabels 416 \nTo determine whether Ensemblex’s improved performance across the in silico pools is reflected in 417 \nreal-world multiplexed experiments,  we applied Ensemblex to an experimentally multiplexed pool 418 \ncomposed of NSCLC dissociated tumor cells from 7 donors, hereafter referred to as the NSCLC 419 \ndataset (18). Importantly, these NSCLC cells were labelled with donor-specific Cell Multiplexing 420 \nOligonucleotides (CMOs) , pro viding a proxy for ground -truth sample labels to evaluate the 421 \nperformance of the genetic demultiplexing tools. For this experiment, we used HTOdemux (19) to 422 \nassign the cells back to their donor -of-origin based on the CMO expression profiles. HTOdemux 423 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 21 \nconfidently assigned 19,695 cells, of which 15,534 (78.87%) were assigned to individual donors 424 \nand 4,161 (21.13%) were assigned as doublets; 769 cells (3.76%) were unassignable at a positive 425 \nquantile of 0.99 and were excluded from downstream analyses (Figures 4A). Application of the 426 \nEnsemblex framework with prior genotype information to the NSCLC dataset achieved a singlet 427 \ntrue positive (TP)  rate of 96.92% and doublet TP rate of 66.21% ( Figure 4B). To evaluate the 428 \nbenefits of utilizing the entire Ensemblex workflow (Step s 1-3), we investigated the contribution 429 \nof each step of the Ensemblex framework to the overall demultiplexing accuracy. Applying graph-430 \nbased doublet detection (Step 2) and ensemble -independent doublet detection (Step 3) to the 431 \naccuracy weighted assignments obtained from Step 1 increased the proportion of correctly 432 \nidentified doublets by 14%, while slightly decreasing the proportion of correctly classified singlets 433 \nby 0.05% (Additional File 1: Table S1 ). Although users can elect to utilize different step -434 \ncombinations of the Ensemblex pipeline, these results reaffirm that leveraging the entire workflow 435 \nmaximizes the overall demultiplexing accuracy by achieving a meticulous balance between singlet 436 \nclassification and doublet identification.  437 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 22 \n 438 \nFigure 4. Evaluating Ensemblex on experimentally multiplexed cells using donor -specific 439 \noligonucleotide labels as a proxy for ground -truth. Non-small cell lung cancer  (NSCLC) 440 \ndissociated tumor cells from 7 individuals were pooled and labelled with donor-specific 441 \noligonucleotide-labels. Cells were demultiplexed according to their expression of donor -specific 442 \noligonucleotide labels by HTOdemux ; HTOdemux’s sample labels were used as a proxy for 443 \nground truth.  True positives (TP)  singlets were defined as cells classified as singlets by both 444 \nHTOdemux and Ensemblex with matching sample labels; false positives (FP)  singlets were 445 \ndefined as cells classified as singlets by both HTOdemux and Ensemblex but assigned to different 446 \ndonors. TP doublets were defined as cells classified as doublets by both HTOdemux and 447 \nEnsemblex; FP doublets were defined as cells classified as singlets by HTOdemux and doublets 448 \nby Ensemblex; false negatives (FN)  doublets were defined as cells classified as doublets by 449 \nHTOdemux and singlets by Ensemblex. A) T-distributed Stochastic Neighbor Embedding (t-SNE) 450 \nvisualization of HTOdemux’s sample labels. B) T-SNE visualization of  Ensemblex’s 451 \ndemultiplexing performance using HTOdemux’s sample labels as ground truth for singlets (left) 452 \nand doublets (right) . C) Bar plots showing  the s inglet TP and FP rate s for each genetic 453 \ndemultiplexing tool using HTOdemux’s sample labels as ground truth . D) Bar plots showing the 454 \ndoublet TP and FP rates for each genetic demultiplexing tool using HTOdemux’s sample labels as 455 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 23 \nground truth. E) Scatter plot showing the proportion of usable cells (confidently classified singlets) 456 \nand the corresponding usable cell error rate for each genetic demultiplexing tool. F) Adjusted Rand 457 \nIndex, balanced accuracy, Matthew’s Correlation Coefficient, and area under the receiver operating 458 \ncharacteristic curve (AUC) of the singlet assignment probability for each genetic demultiplexing 459 \ntool. 460 \n 461 \nUpon comparing Ensemblex’s demultiplexing performance with prior genotype information  on 462 \nthe NSCLC dataset to the individual genetic demultiplexing tools, it emerged that our ensemble 463 \nmethod obtained the highest singlet and doublet TP rates (Figures 4C and 4D). Ensemblex and 464 \nDemuxlet also showed the lowest singlet false positive (FP) rates (0.25% and 0.21%, respectively), 465 \nindicating that singlets were least frequently assigned to the wrong donor -of-origin by these two 466 \nmethods compared to Demuxalot  (1.87%), Vireo-GT (3.91%), and Souporcell ( 11.94%). 467 \nSouporcell and Vireo-GT returned the highest proportion of usable  cells (confidently classified 468 \nsinglets; 88.21% and 86.51%, respectively); albeit, at the expense of high usable cell error rates 469 \n(22.91% and 13.53%, respectively) (Figure 4E). In turn, Ensemblex, Demuxalot, and Demuxlet  470 \nshowed lower error rates across the usable cells (8.75%, 8.91%, and 9.51%, respectively), amongst 471 \nwhich Ensemblex returned the highest proportion of usable cells (83.77%) compared to Demuxalot 472 \n(83.64%) and Demuxlet ( 83.43%). Here, the relatively high error rate amongst usable cells 473 \nreturned by each demultiplexing tool is attributed to true doublets classified as singlets . Finally, 474 \nwe computed the ARI, balanced accuracy, MCC, and AUC for singlet detection for each tool and 475 \nobserved that Ensemblex again outperformed the remaining tools  (Figure 4F). We repeated the 476 \nabove analyses without prior genotype information  and observed a similar trend towards better 477 \noverall performance by Ensemblex (Additional File 1: Table S2 and Figure S12). Together, these 478 \nresults corroborate that Ensemblex’s improved performance on the in silico  pools extends to 479 \nexperimentally multiplexed samples. 480 \n 481 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 24 \nApplication of Ensemblex to experimentally pooled, highly multiplexed subjects 482 \nTo evaluate Ensemblex’s demultiplexing performance on  experimentally pooled, highly 483 \nmultiplexed scRNAseq datasets with prior genotype information, we used pools containing iPSC 484 \nlines from 22 donors that were differentiated towards DaN by Jerber et al., hereafter referred to as 485 \nthe DaN dataset (12) (Figure 5A). To capture the transcriptional changes throughout neurogenesis, 486 \nJerber et al. performed scRNAseq of the iPSC lines grown in pooled cultures at days 11, 30, and 487 \n52 of differentiation  (Figure 5A ). Using three technical replicates from each timepoint, we 488 \nobtained 84,746 cells after performing quality control as previously described (12) (Additional 489 \nFile 1: Table S3). Each technical replicate was demultiplexed independently by Ensemblex and 490 \nits constituent tools.  491 \n 492 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 25 \nFigure 5. Application of Ensemblex to highly multiplexed, experimentally pooled cultures of 493 \ndifferentiated dopaminergic neurons.  A) Time line of iPSC pooling, dopaminergic neuron 494 \n(DaN) differentiation, and sample collection from the DaN dataset by Jerber et al. (12). Three 495 \ntechnical replicates at each time point (days 11, 30 and, 52 of differentiation) from pools containing 496 \n22 individual iPSC lines were used in the analysis. Across all timepoints and technical replicates, 497 \n84,746 cells were obtained for analysis.  B) Uniform manifold approximation and projection 498 \n(UMAP) plots showing confidently assigned singlets or predicted doublets (blue) and ambiguous 499 \nsinglets (singlet assignments with insufficient assignment probabilities; red)  returned by each 500 \ndemultiplexing tool. C) Stacked bar chart showing the proportion of confidently assigned singlets 501 \nor predicted doublets (blue) and ambiguous singlets (red) across technical replicates at each time 502 \npoint returned by each demultiplexing tool.  D) Boxplot showing the proportion  of confidently 503 \nclassified singlets across technical replicates  and time points  by each demultiplexing tool. 504 \nWilcoxon rank-sum tests were used to compare the proportion of confidently classified singlets by 505 \nEnsemblex to that of its constituents  (n = 9 pools). E) Bar chart showing the  proportion of 506 \noverlapping ambiguous singlet assignments  amongst demultiplexing tools across technical 507 \nreplicates and time points (n = 9 pools) . F) Boxplot showing the Adjusted Rand Index (ARI) 508 \nassessing cluster stability across a range of 11 clustering resolutions (n clustering iterations = 25) 509 \nafter removing doublets identified by each demultiplexing tool. Wilcoxon rank -sum tests were 510 \nused to compare the clustering ARI after removing Ensemblex doublets to the clustering ARI after 511 \nremoving doublets identified by each constituent tool. * Adjusted P-value < 0.05; ** adjusted P -512 \nvalue < 0.01; *** adjusted P-value < 0.001 513 \n 514 \nTo characterize the relationship between Ensemblex and its constituent demultiplexing tools, we 515 \ncomputed the ARI between Ensemblex’s sample labels and those of its constituent as well as the 516 \npercent contribution of each tool to Ensemblex’s final sample labels ( Table 2). Notably, we 517 \nobserved that across day 30 technical replicates  Demuxlet showed an ARI of 0.063 with 518 \nEnsemblex and only contributed 29.74% to Ensemblex’s final sample labels. In contrast, across 519 \nday 11 and 52 technical replicates Demuxlet showed an ARI of 0.928 and 0.884, respectively, and 520 \ncontributed 95.91% and 90.55%, respectively,  to Ensemblex’s final sample label s. Importantly, 521 \nDemuxlet’s variable contribution to Ensemblex’s sample labels across sequencing time points 522 \ndemonstrates our ensemble method’s ability to adapt to the relative performance of its constituent 523 \ntools and override the classifications of a poorly performing tool on the respective dataset.   524 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 26 \nTable 2. Application of Ensemblex to pooled cultures of dopaminergic neurons from 22 525 \nhealthy controls.  526 \n ARI between Ensemblex and \nconstituent tool assignments \n Percent contribution to \nEnsemblex assignments \nn \nusable cells \nn \ndoublets  Day 11 Day 30 Day 52  Day 11 Day 30 Day 52 \nDemuxalot 0.987 0.955 0.982  97.29% 94.75% 97.57% 75,962 8,279 \nDemuxlet 0.928 0.062 0.884  95.91% 29.74% 90.55% 57,567 6,614 \nSouporcell 0.883 0.876 0.912  91.62% 91.82% 93.84% 76,811 7,740 \nVireo-GT 0.961 0.879 0.958  95.95% 88.80% 95.16% 75,933 6,115 \nEnsemblex NA NA NA  NA NA NA 76,222 8,307 \nDoubletFinder NA NA NA  NA NA NA NA 4,597 \nPooled cultures of induced pluripotent stem cell (iPSC) lines from 22 healthy donors were 527 \ndifferentiated towards a dopaminergic neuron (DaN) fate and sequenced on days 11, 30, and 52 of 528 \ndifferentiation by Jerber et al. (12). For the analysis we used three technical replicates for each 529 \nsequencing timepoint. Each pool was demultiplexed independently by Ensemblex and its 530 \nconstituent tools with prior genotype information. The Adjusted Rand Index (ARI)  between 531 \nEnsemblex’s assignments and those of the constituent tools was computed across technical 532 \nreplicates corresponding to each differentiation timepoint. The percent contribution represents the 533 \nproportion of assignments from each constituent tool that matched Ensemblex’s assignments. 534 \nUsable cells were defined as singlet classifications whose assignment probability exceeded the 535 \nrecommended threshold of the respective tool. Abbreviations: NA = Not applicable. 536 \n 537 \nTo elucidate the discrepancy in Demuxlet’s contribution to Ensemblex’s sample labels  across 538 \nsequencing time points, we  investigated the proportion of ambiguous singlet assignment s from 539 \nEnsemblex and its constituents. Ambiguous singlets are defined as singlet classifications whose 540 \nassignment probabilities failed to meet the recommended threshold of the respective tool, leaving 541 \nthe identity of the pooled cell unresolved. Across 84,746 cells, Souporcell (195 singlets; 0.23% of 542 \ncells) and Ensemblex (217 singlets; 0.26% of cells) showed the lowest proportion of ambiguous 543 \nsinglet assignments, followed by Demuxalot (505 singlets; 0.60% of cells) and Vireo-GT (2,698 544 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 27 \nsinglets; 3.18% of cells ). Strikingly, Demuxlet showed 20,565 ambiguous singlet assignments 545 \n(24.27% of cells), with 92.04% derived from day 30 technical replicates, reflecting Demuxlet’s 546 \nremarkably low contribution to Ensemblex’s sample labels for cells sequenced at this timepoint  547 \n(Figures 5B and 5C). In accordance with previous analyses  (9, 10), Demuxlet was consistently 548 \namongst the top performing constituent tools throughout our benchmarking analyses. Yet, its poor 549 \nperformance across day 30 technical replicates illustrates how the accuracy of individual tools can 550 \nvary greatly between datasets, highlighting the importance of utilizing multiple distinct algorithms 551 \nfor genetic demultiplexing. We compared the mean proportion of confidently classified singlets 552 \nacross technical replicates from each time point (n = 9) between Ensemblex (99.72%) and each 553 \nconstituent demultiplexing tool  using a Wilcoxon rank-sum test. After correction for multiple 554 \nhypothesis testing, w e observed that the mean proportion of confidently classified singlets by 555 \nEnsemblex was  significantly higher than Demuxalot (mean = 99.36%, P -value = 3.55e-3), 556 \nDemuxlet (mean = 75.82%, P-value = 1.55e-5), and Vireo-GT (mean = 96.71%, P-value = 1.55e-557 \n5) (Figure 5D). Thus, despite Demuxlet’s unusually poor performance across day 30 technical 558 \nreplicates, Ensemblex still confidently classified 27,520 singlets (99.61% of singlet assignments) 559 \nfrom these pools. Indeed, our ensemble method mitigates the consequences of a poorly performing 560 \nconstituent tool by outweighing the erroneous classifications. In contrast , using a consensus 561 \nframework returned only 7,446 confidently classified singlets  from day 30 technical replicates 562 \n(20,074 fewer cells than Ensemblex), limiting the availability of data for downstream analyses. 563 \n 564 \nTo further evaluate the ambiguity amongst singlet classification, we investigated the intersection 565 \nof ambiguous singlets across demultiplexing tools, reasoning that cells that are most challenging 566 \nto demultiplex would be labelled as ambiguous across all tools (Figure 5E). The singlets that were 567 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 28 \nassigned as ambiguous by Ensemblex  showed the highest ambiguous singlet rate across the 568 \nremaining tools (mean across all constituent tools = 73.04%; mean across Demuxalot, Demuxlet, 569 \nand Vireo-GT = 92.32%). In contrast, while Souporcell showed the lowest ambiguous singlet rate 570 \noverall, only 15.90% of its unassigned singlets, on average, were ambiguous across the remaining 571 \ntools. These results indicate that the cells labelled as ambiguous by Ensemblex represent the cells 572 \nthat are most challenging to classify across the distinct demultiplexing algorithms. Indeed, limiting 573 \nEnsemblex’s ambiguous singlet assignments to those that are most difficult to classify is critical 574 \nfor maintaining a balance between maximizing the number of usable cells and minimizing the 575 \nintroduction of technical artifacts into downstream analyses from misclassified cells.   576 \n 577 \nNext, we compared the doublet predictions made by each genetic demultiplexing tool  and 578 \nDoubletFinder, a doublet detection tool that predicts doublets by estimating the similarity of the 579 \ntranscriptional profile of a pooled cell to artificial doublets generated by combining the 580 \ntranscriptional profiles of randomly selected cell pairs  (20). Although the average number of 581 \nunique molecular identifiers (UMI) per cell across doublets identified by each tool was 582 \nsignificantly higher than the consensus singlets ( Additional File 1: Figure S13), we observed a 583 \nnotable discrepancy in the number of doublets identified by each tool ; DoubletFinder identified 584 \nthe fewest doublets (n = 4,597), while Ensemblex identified the most doublets (n = 8,307) (Table 585 \n1). Accordingly, all tools identified doublets that every other tool assigned as singlets (Additional 586 \nFile 1: Figure S13). While Ensemblex identified the highest number of doublets, it still returned 587 \na higher number of confidently classified singlets  (n =  76,222) than Demuxalot ( n =  75,962), 588 \nDemuxlet (n = 57,567), and Vireo-GT (n = 75,933). Thus, even though the Ensemblex framework 589 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 29 \nprioritizes the identification of doublets at the expense of a slight ly lower singlet classification 590 \nrate, our ensemble method still returns a high proportion of usable cells for downstream analyses.   591 \n 592 \nTo evaluate the impact of doublet removal on the stability of clusters  in the DaN dataset, we 593 \nperformed 25 different random start iterations of the  Louvain network detection at various 594 \nclustering resolutions after removing the doublets identified by each tool (21). Removing the 595 \ndoublets identified by Ensemblex resulted in the  highest ARI (mean ARI = 0.942), on average, 596 \nacross clustering resolutions  (Figure 5F ), suggesting the greatest cluster stability . However, 597 \nWilcoxon rank -sum tests only revealed a statistically significant difference in the cluster 598 \nassignment ARI between Ensemblex and Souporcell (mean ARI = 0.922, P-value = 1.08e-2) after 599 \ncorrection for multiple hypothesis testing. Nonetheless, the highest cluster stability after removal 600 \nof Ensemblex’s putative doublets illustrates how improved doublet detection can translate to 601 \nimproved biological analyses and is reflective of its superior doublet identification performance 602 \non the benchmarking analyses. 603 \n 604 \nEvaluating the impact of demultiplexing tools on differential gene expression analysis 605 \nTo evaluate the impact of genetic demultiplexing tool s on scRNAseq DGE analysis, we 606 \nmultiplexed iPSC-derived NSCs from individuals with ADHD and controls (Figure 6A). NSCs 607 \nwere pooled and cultured until 100% confluence was reached. Two multiplexing experiments were 608 \nperformed: Experiment 1 (n ADHD = 7; n control = 6) and Experiment 2 (n ADHD = 9; n control 609 \n= 7). After filtering cells for > 500 total and unique RNA transcripts , we obtained 30,433 cells 610 \nacross both pools. Louvain clustering on the integrated scRNAseq dataset identified 12 clusters, 611 \nwhich were annotated as eight putative cell types (Figure 6B). 612 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 30 \n 613 \nFigure 6. Evaluating the impact of discordant assignments between genetic demultiplexing 614 \ntools on d ifferential gene expression analysis. A) Schematic illustrating the workflow for the 615 \nneural stem cell (NSC) dataset. Pooled induced pluripotent stem cell (iPSC) -derived neural stem 616 \ncell cultures from individuals with attention deficit hyperactivity disorder (ADHD) and controls 617 \nwere collected in two separate experiments. NSCs were dissociated for single -cell RNA 618 \nsequencing and prior genotype information of the pooled subjects was obtained through 619 \nmicroarray genotyping. The pools were demultiplexed by Ensemblex and its constituents with 620 \nprior genotype information and differential gene expression (DEG) was computed between ADHD 621 \nand controls. B) Uniform manifold approximation and projection (UMAP) plot showing the 622 \nputative cell types. C) Summary of the number of usable cells — singlets above the recommended 623 \nprobability threshold of the respective demultiplexing tool — assigned to ADHD donors and 624 \ncontrols and the number of identified doublets by each demultiplexing tool. D) Boxplot showing 625 \nthe Adjusted Rand Index (ARI) assessing cluster stability across a range of 11 clustering 626 \nresolutions (n clustering iterations = 25) after removing doublets identified by each demultiplexing 627 \ntool. A one-way Analysis of Variance (ANOV A) test comparing the ARI after removing doublets 628 \nidentified by each tool revealed a significant difference between tools (n = 11 clustering 629 \nresolutions; P-value = 1.18e-3). E) Proportion of ADHD and control cells identified as putative 630 \ndoublets by Ensemblex that were assigned as singlets by the constituent demultiplexing tools. F) 631 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 31 \nHeatmap showing the number of cell-type specific DEGs between ADHD and controls using the 632 \nsubject labels of each demultiplexing tool. G) Heatmap showing the number of cell-type specific 633 \nDEGs between ADHD and controls  using the subject labels of each demultiplexing tool  and 634 \nremoving putative doublets identified by Ensemblex. Cell-types not shown in the heatmaps had no 635 \nDEGs passing the adjusted P-value < 0.01 and |Log2FC >= 0.5| threshold across all tools. 636 \n 637 \nWe independently demultiplexed both pools using Ensemblex and its constituents to assign the 638 \ncells back to their donor -of-origin with prior genotype information (Figure 6C). The number of 639 \ncells assign ed to ADHD and control donors by each genetic demultiplexing  tool is shown in 640 \nAdditional File 1: Table S6. Importantly, the NSC dataset provides a valuable illustration of the 641 \nconsequences of unnecessarily discarding cells from downstream analyses. For example, 642 \nEnsemblex and Vireo -GT returned 2,387 and 882 confidently assigned GRIA1high neurons, 643 \nrespectively, whereas a consensus approach would have confidently assigned only 563 GRIA1high 644 \nneurons (Additional File 1: Table S6).  645 \n 646 \nEach genetic demultiplexing tool predicted the ADHD cells to be vastly underrepresented 647 \ncompared to the control cells; Ensemblex assigned 2,739 cells to individuals with ADHD and 648 \n19,880 cells to controls, suggesting that the ADHD iPSC lines were lost throughout the culturing 649 \nand sequencing process (Figure 6C). Additionally, we observed a notable difference in the number 650 \nof identified doublets across the tools; Vireo-GT identified the fewest doublets (n = 2,707), while 651 \nDemuxlet identified the most doublets ( n = 8,329) (Figure 6C). We aimed to characterize the 652 \nchange in cluster stability after removing the doublets identified by each tool and observed that 653 \nremoving the doublets identified by Ensemblex resulted in the highest ARI (mean ARI = 0.995), 654 \non average, across clustering resolutions  (Figure 6D ). A one-way ANOV A test comparing the  655 \nclustering ARI after removal of doublets identified by each tool revealed a significant difference  656 \nbetween tools (P-value = 1.18e-3). Demuxlet (n = 8,329) identified more doublets than Ensemblex 657 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 32 \n(n = 6,373), but exhibited lower cluster stability (ARI), suggesting that increased cluster stability 658 \nis not merely representative of the number of doublet s removed but rather the quality of doublet 659 \nremoval.  660 \n 661 \nGiven the underrepresentation of ADHD cells across the dataset, we elected to investigate the cells 662 \nthat were identified as doublets by Ensemblex but assigned as singlets by the constituent tools and 663 \nhow these putative  doublets were distributed across  samples according to disorder status. 664 \nDemuxalot (n = 388) and Demuxlet (n = 726) assigned a relatively low number of Ensemblex’s 665 \ndoublets as singlets , which represented 0.66% and 4.58% of ADHD sample assignments, 666 \nrespectively, and 1.97% and 3.58% of control sample assignments, respectively (Figure 6E). In 667 \ncontrast, Souporcell (n = 3,902) and Vireo-GT (n = 1,334) assigned a relatively high number of 668 \nEnsemblex’s doublets as singlets , which represented 31.97% and 24.88% of ADHD sample 669 \nassignments, respectively, and 11.65% and 3.97% of control sample assignments, respectively, 670 \nillustrating how variable doublet detection can impact the assembly of cells assigned to don or 671 \ncategories and which cells are retained for downstream analyses.  672 \n 673 \nFinally, w e used the  model-based analysis of single -cell transcriptomics  (MAST) statistical 674 \nframework to compute cell-type specific DGE between individuals with ADHD and controls using 675 \nthe demultiplexed sample labels from each tool (22). We observed a significant discrepancy in the 676 \nnumber of cell type -specific differentially expressed gene s (DEGs; adjusted P-value < 0.01; 677 \nabsolute log2 fold change > 0.5)  depending on the demultiplexing tool used (Figure 6F). Most 678 \nnotably, for glia cells Souporcell identified 116 DEGs; Vireo-GT identified 98 DEGs; Ensemblex 679 \nidentified 7 DEGs; Demuxalot identified 6 DEGs; and Demuxlet identified 1 DEG. Similar 680 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 33 \npatterns were observed  across SOX2high NSCs, POU5F1high neural progenitor cells (NPC) , 681 \nS100Bhigh NPCs, and DCXhigh neurons, whereby Souporcell or Vireo-GT’s sample labels resulted 682 \nin a remarkably high number of DEGs compared to Ensemblex, Demuxalot, and Demuxlet. Given 683 \nthat Souporcell and Vireo-GT made relatively few doublet call s and that 31.97% and 24.88% of 684 \nADHD sample assignments made by Souporcell and Vireo -GT, respectively,  were putative 685 \ndoublets identified by Ensemblex, we elected to repeat the DGE analysis using the demultiplexed 686 \nsample labels from each tool but this time we removed all putative doublets identified by 687 \nEnsemblex. In doing so, we observed a decrease in the number of DEGs identified by Souporcell 688 \nand Vireo-GT across cell types , suggesting that the putative doublets  identified by Ensemblex, 689 \nwhich were classified as singlets by Souporcell and Vireo -GT, were driving the initial  signals 690 \n(Figure 6G ). For example, the number of glia-specific DEGs decreased from 116 to 0 with 691 \nSouporcell’s sample labels, and 98 to 0 with Vireo-GT’s sample labels. Given that the NSC dataset 692 \nlacked ground -truth sample labels, we could not definitively determine which cells were true 693 \ndoublets; however, the increase in clustering ARI after removal of Ensemblex’s putative doublets 694 \n(Figure 6D), coupled with Ensemblex’s improved doublet identification performance on pools 695 \nwith known ground -truth sample labels ( Figure 2B ), afforded confidence to assume that our 696 \nensemble method performed favorably. Nonetheless, this analysis reveals that the choice of 697 \ndemultiplexing tool can greatly impact biological analyses. 698 \n 699 \nConclusion 700 \nMultiplexing protocols, coupled with the introduction of genetic demultiplexing tools constituted 701 \na significant advancement for scRNAseq  by providing a feasible means to dramatically increase 702 \nthe throughput of biological replicates . As the demand for population-scale scRNAseq analysis 703 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 34 \ncontinues to grow with the maturation of singe-cell technologies, the prospect of multiplexing 704 \nentire cohorts has emerged. However, the realization of this goal is impeded by the limitations of 705 \nthe current genetic demultiplexing tools. These include decreasing demultiplexing performance as 706 \nthe number of multiplexed samples increase s (9, 10) , relatively poor doublet detection 707 \nperformance (10), relatively high rates of cells that can only be correctly classified by single 708 \nalgorithms, the unnecessary removal of correctly classified cells due to insufficient assignment 709 \nprobabilities, and highly variable demultiplexing performance between datasets (10). In this work 710 \nwe presented Ensemblex, which offers a unique solution to these limitations by meticulously 711 \nimplementing distinct demultiplexing algorithms into a robust, accuracy -weighted ensemble 712 \nframework that is exceptionally equipped to classify highly multiplexed pools. 713 \n 714 \nWe applied Ensemblex to a diverse array of computationally and experimentally multiplexed 715 \nscRNAseq datasets. Benchmarking analyses on pools with known ground -truth sample labels 716 \nrevealed Ensemblex’s superior demultiplexing performance across pools reaching 80 multiplexed 717 \nsamples, which translated to a higher proportion of cells retained for downstream analyses and 718 \nlower error rates amongst classified cells. Ensemblex also demonstrated a notable advancement 719 \nfor identifying heterogenic doublets, which is a well -documented limitation of the genetic 720 \ndemultiplexing tools currently available (9, 10, 15) . While previous analyses indicated that the 721 \nnumber of multiplexed samples in a pool directly impacted doublet detection efficiency (15), we 722 \nshowed that Ensemblex’s ability to identify doublets remained relatively constant when >24 723 \nsamples were multiplexed. Our findings suggest that super loading cells prior to sequencing —724 \nwhich will result in a higher number of usable cells but a higher a doublet rate (6) — followed by 725 \nheterogenic doublet detection by Ensemblex, may be a viable approach for implementing 726 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 35 \npopulation-scale multiplexing in practice. We also demonstrated that the performance of individual 727 \ngenetic demultiplexing tools can be highly dataset -dependent, reflecting the findings of previous 728 \nwork (10). However, due to its unsupervised weighting model, we showed that Ensemblex is 729 \nresistant to poorly performing constituent tools, maximizing the consistency of its demultiplexing 730 \nperformance. Nonetheless, if each constituent tool performs poorly  on a given dataset, the poor 731 \nperformance will be reflected in Ensemblex’s demultiplexing accuracy. Finally, we illustrated that 732 \ndiscordant sample assignments amongst genetic demultiplexing tools can greatly impact DGE 733 \nanalyses, necessitating that investigators carefully consider their choice of genetic demultiplexing 734 \ntool. Although untested, we anticipate that the impacts of discordant sample assignments amongst 735 \ngenetic demultiplexing tools on biological interpretations would be exacerbated for computational 736 \nanalyses that consider the specific donor identity of the  pooled cells, such as expression 737 \nquantitative trait loci (eQTL) analyses, as opposed to donor groups (i.e., case and control). Due to 738 \nEnsemblex’s ability to seamlessly integrate multiple algorithms into an adaptable framework, we 739 \nargue that our ensemble method achieves unmatched reliability for experimentally multiplexed 740 \npools that lack ground truth sample labels. 741 \n 742 \nUndoubtedly, a limitation of utilizing an ensemble method for genetic demultiplexing is the 743 \nnecessity to run  each individual demultiplexing algorithm, which can be  computationally 744 \nexpensive. Yet, in the absence of comparing demultiplexed sample labels across tools, poor 745 \nperformance by a given individual algorithm on experimentally multiplexed pools is undetectable, 746 \nand the risk of introducing technical artifacts and losing usable cells  for downstream analyses is 747 \nprominent. As such, we believe that the relatively high computational cost of Ensemblex is a 748 \nworthwhile investment to maximize the biological insight obtained from multiplexed scRNAseq 749 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 36 \ndatasets. To mitigate the burden of genetic demultiplexing by multiple individual tools, we provide 750 \na coherent pipeline that runs each constituent demultiplexing tool in parallel and seamlessly 751 \nprocesses the respective output files with the Ensemblex algorithm. 752 \n 753 \nCompared to when demultiplexing was informed by prior genetic data of the pooled samples, the 754 \nimprovement of Ensemblex over its constituent tools was far less pronounced for genotype-free 755 \ndemultiplexing cases. A ll demultiplexing tools, including Ensemblex, showed drops in 756 \ndemultiplexing performance when >16 samples were multiplexed in a pool without prior genotype 757 \ninformation. Nonetheless, Ensemblex still constitutes an advancement over the individual tools for 758 \ngenotype-free demultiplexing cases due to the robustness achieved by incorporating distinct 759 \ndemultiplexing algorithms, which protects against the prospect of poorly performing individual 760 \ntools on the respective dataset. Furthermore, an intrinsic limitation of demultiplexing without prior 761 \ngenotype information is that samples cannot be directly linked to metadata, leaving the sample 762 \nidentity of the inferred clusters unresolved (9). Although challenging, this limitation can be 763 \nmitigated by identifying a small subset of discriminatory variants from the reconstructed genotypes 764 \nof the constituent demultiplexing tools, which could be used to manually assign the computed 765 \nclusters to samples if such discriminatory variants are known by the investigator. While the 766 \nEnsemblex pipeline provides users the option to demultiplex pools with or without prior genotype 767 \ninformation, we assert that users take caution when electing to perform population -scale 768 \nmultiplexing experiments without using prior genetic data.  769 \n 770 \nGenetic demultiplexing tools have been used extensively for scRNAseq analysis across many 771 \ndisciplines in the biological sciences, including microbiology (8), model organisms (15), cancer 772 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 37 \nbiology (23), and neurodegenerative disease (12). Recent work has also evaluated the utility of 773 \ngenetic demultiplexing tools for different single-cell, read-based modalities such as single-nuclei 774 \nRNA sequencing (snRNAseq) and single -nuclei assay for transposase -accessible chromatin 775 \nsequencing (scATACseq) (24). Although untested, we expect Ensemblex to prove beneficial in 776 \ndemultiplexing for these assays, but comprehensive benchmarking with the appropriate datasets is 777 \nrequired and was not explored here.  778 \n 779 \nWe expect numerous biological fields to exploit the benefits of Ensemblex through its application 780 \nto highly multiplexed pools comprising cells from  many genetically distinct individuals. 781 \nSpecifically for biomedical sciences, the preparation and labour costs of scRNAseq remains 782 \nprohibitively expensive for analyzing entire cohorts of patients, which is critical for characterizing 783 \nthe genetic heterogeneity and etiological diversity of disease, and for maint aining sufficient 784 \nstatistical power for detecting associations between transcriptional changes and clinical or 785 \npathological observations (1). By increasing the throughput of biological replicates, multiplexing 786 \nhas rendered the prospect of analyzing entire patient cohort s with single -cell transcriptomics 787 \nfeasible. Highly-multiplexed scRNAseq experiments have already been presented in the literature 788 \nand, to the best of our knowledge, have pooled up to 24 samples in a single dish (12). However, 789 \nwe demonstrated that Ensemblex’s demultiplexing accuracy remains relatively constant when >24 790 \nsamples are multiplexed at concentrations that abide by the current limitations of experimental 791 \nprotocols, suggesting that  Ensemblex equips  the research community with the necessary 792 \ncomputational framework to expand the upper limits of the number of genetically distinct 793 \nindividuals in a single pool.  794 \n 795 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 38 \nWhile multiplexing mitigates the labour and consumable costs of scRNAseq analysis, the cost of 796 \nsequencing remains expensive and the increasing number of genetically distinct individuals in a 797 \nsingle pool necessitates that a  greater number of cells must be sequenced to ensure adequate 798 \nrepresentation. Accordingly, Ensemblex is equipped to demultiplex pools comprising cells from 799 \nmore genetically distinct individuals than is feasible with the current laboratory technologies. 800 \nHowever, we expect that the cost of sequencing will continue to decrease with the maturation of 801 \nthe technology, and our tool will be in place for when the  anticipated wet lab advancements are 802 \nrealized. Overall, we conclude that Ensemblex constitutes a notable advancement towards the 803 \npressing demand for population-scale single-cell transcriptomics. 804 \n 805 \nMethods 806 \nEnsemblex framework overview 807 \nEnsemblex is an ensemble genetic demultiplexing framework for scRNAseq sample  pooling that 808 \nwas designed to identify the most probable sample labels from each of its constituent tools : 809 \nDemuxalot (5), Demuxlet (6), Souporcell (8), and Vireo (9) when demultiplexing with prior 810 \ngenotype information or Demuxalot, Freemuxlet (6), Souporcell, and Vireo when demultiplexing 811 \nwithout prior genotype information. After running each constituent demultiplexing tool in parallel, 812 \nEnsemblex merges the output files  containing the sample-cell assignments from each tool  and 813 \nperforms three distinct steps of the Ensemblex pipeline: 814 \n1. Accuracy-weighted probabilistic ensemble; 815 \n2. Graph-based doublet detection;  816 \n3. Ensemble-independent doublet detection. 817 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 39 \nUpon obtaining the final Ensemblex sample labels (donor-of-origin identity of the pooled cells) , 818 \nthe singlet assignment confidence score is computed.  819 \n 820 \nStep 1: Accuracy-weighted probabilistic ensemble 821 \nEnsemblex utilizes an unsupervised weighting model to identify the most probable sample 822 \nlabel for each cell. Ensemblex weighs each constituent tool’s assignment probability 823 \ndistribution by its estimated balanced accuracy for the dataset in a framework adapted from 824 \nthe work of Large et al. (16). To estimate the balanced accuracy of a particular constituent tool 825 \n(e.g., Demuxalot) for experimentally multiplexed  datasets lacking ground -truth labels, 826 \nEnsemblex uses the cells with a consensus assignment across the three remaining tools (e.g. , 827 \nDemuxlet, Souporcell, and Vireo-GT) as a proxy for ground-truth. The balanced accuracy for 828 \neach tool is calculated using equation 1: 829 \n 830 \n(1) 𝐵𝑎𝑙𝑎𝑛𝑐𝑒𝑑 𝑎𝑐𝑐𝑢𝑟𝑎𝑐𝑦=  \n1\n2 ((\n𝑇𝑃\n𝑇𝑃+𝐹𝑁) + (\n𝑇𝑁\n𝑇𝑁+𝐹𝑃)) 831 \n 832 \nWhere TP is the number of correctly classified singlets; true -negative (TN) is the number of 833 \ncorrectly classified doublets; FP is the number of incorrectly classified singlets; false- negative 834 \n(FN) is the number of incorrectly classified doublets . The probability distribution of each 835 \nconstituent tool ( 𝑝𝑗̂) is then weighted by its estimated balanced accuracy ( 𝑤𝑗) to produce an 836 \naccuracy-weighted ensemble probability for each cell: 837 \n 838 \n(2) 𝑝̂(𝑦 = 𝑖|𝐸) ∝  ∑ 𝑤𝑗𝑝̂𝑗(𝑦 = 𝑖|𝑀𝑗)𝑘\n𝑗=1  839 \n 840 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 40 \nWhere 𝑝̂ is the probability that a barcode belongs to class 𝑖; 𝑦 is the class variable with 𝑐 841 \npossible values, 𝑦 ∈ (1, … , 𝑐); 𝑐 is the number of pooled samples plus 1 to account for 842 \ndoublets; 𝐸 is a vector of the results  of 𝑀 classifiers, 𝐸 = (𝑀1, … , 𝑀𝑘); 𝑀is the individual 843 \nconstituent demultiplexing output from each tool. Given 𝑝̂, Ensemblex assigns each barcode’s 844 \nsample identity (𝑦̂) as the class (sample label) with the maximum probability: 845 \n 846 \n(3) 𝑦̂= arg 𝑚𝑎𝑥𝑖∈(1,…,𝑐) 𝑝̂(𝑦 = 𝑖|𝐸) 847 \n  848 \nStep 2: Graph-based doublet detection 849 \nEnsemblex employs a graph-based approach to identify doublets that are incorrectly labeled as 850 \nsinglets by the accuracy-weighted probabilistic ensemble component (Step 1). For graph-based 851 \ndoublet detection, Ensembl ex leverages pre-defined features returned from each constituent 852 \ntool: 853 \n1. Demuxalot: doublet probability; 854 \n2. Demuxlet/Freemuxlet: singlet log likelihood – doublet log likelihood; 855 \n3. Demuxlet/Freemuxlet: number of single nucleotide polymorphisms (SNP) per cell; 856 \n4. Demuxlet/Freemuxlet: number of reads per cell; 857 \n5. Souporcell: doublet log probability; 858 \n6. Vireo: doublet probability; 859 \n7. Vireo: doublet log likelihood ratio. 860 \nFor each feature independently, the pooled cells are ordered from the most to the least probable 861 \ndoublet and are assigned a percentile rank. Beginning with a percentile threshold of 99.99 , 862 \nEnsemblex screens each cell to identify those that exceed the percentile threshold across all 863 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 41 \nfeatures; cells that exceed the percentile threshold across all features are labeled as “confident 864 \ndoublets”. For each iteration, Ensemblex decreases the percentile threshold by 0.01 and repeats 865 \nthe screening process until it has identified n confident doublets (nCD). Ensemblex performs 866 \na parameter sweep to determine the optimal nCD to use for graph-based doublet detection (see 867 \nbelow).    868 \n 869 \nNext, the above features are input into a PCA using the stats (v3.6.2) R package (25) and a 870 \nEuclidean distance matrix is generated from the first two principal components (PC). For each 871 \nconfident doublet independently, the remaining cells in the pool are assigned a percentile rank 872 \nbased on their proximity in Euclidean space to the confident doublet and the cells that exceed 873 \nthe designated  nearest neighbour percentile threshold ( pT) are identified. For all cells that 874 \nexceeded the designated pT  for any confident doublet  (putative doublets) , Ensemblex 875 \ncomputes the number of times the putative doublet was amongst the nearest neighbours of any 876 \nconfident doublet (fNN); an fNN equal to nCD indicates that a putative doublet was amongst 877 \nthe top nearest neighbours for each confident doublet.  878 \n 879 \nTo optimize the nCD and pT parameters for experimentally pooled samples lacking ground-880 \ntruth labels, Ensemblex performs an automated parameter sweep at each pairwise combination 881 \nof nCD and pT values; nCD values range from 50 to 300, in increments of 50, while pT values  882 \ndepend on the expected doublet rate (exDR) and range from 1 − \n𝑒𝑥𝐷𝑅\n6   to 1 −  𝑒𝑥𝐷𝑅, in 883 \nintervals of  \n1−𝑒𝑥𝐷𝑅\n6 . The distribution of fNN values for each combination of nCD and pT 884 \nparameters are plotted and Pearson’s measure of kurtosis (K), is used to predict which 885 \ncombination of pT and nCD  values optimize the identification of true doublets while 886 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 42 \nminimizing the rate of incorrectly labelled true singlets as doublets. Ensemblex screens for 887 \ncombinations of nCD and pT values that result in negatively skewed fNN distributions with 888 \nhigh K, signifying high peakedness and heavy tails. High peakedness indicates that cells 889 \nexceeding the designated pT concentrated around nCD, reflecting their proximity in Euclidean 890 \nspace to all high confident doublets, while heavy tails indicate that even cells with lower fNN 891 \nvalues were identified as nearest neighbour to many confident doublets.  Ensemblex first 892 \nidentifies the pT that returns the highest K, on average, across nCD values tested in the 893 \nparameter sweep using equation 4:   894 \n 895 \n(4) 𝑝𝑇̂= arg 𝑚𝑎𝑥𝑝𝑇∈{1− 𝑒𝑥𝐷𝑅\n6   ,….,1−𝑒𝑥𝐷𝑅) (\n∑ 𝐾(𝑦=𝑝𝑇)𝑛𝐶𝐷∈{50,100,150,200,250,300}\n2 ) 896 \n 897 \nWhere K of the distribution of fNN values of the putative doublets is defined as: 898 \n 899 \n(5) 𝐾(fNN) = 𝐸[(\n𝑋−𝜇\n𝜎 )\n4\n] 900 \n 901 \nWhere 𝜇 is the mean of the distribution and 𝜎 is the standard deviation. Upon identifying the 902 \noptimal pT value ( 𝑝𝑇̂), Ensemblex plots the K corresponding to 𝑝𝑇̂ across all nCD values 903 \ntested in the parameter sweep. If an inflection point is identifiable, Ensemblex identifies 𝑛𝐶𝐷̂ 904 \nas the nCD value corresponding to the point of inflection on the curve. Otherwise, Ensemblex 905 \nidentifies 𝑛𝐶𝐷̂ as the nCD value corresponding to the highest K. Cells flagged as putative 906 \ndoublets identified using 𝑝𝑇̂ and 𝑛𝐶𝐷̂ are labelled as doublets by Ensemblex. 907 \n 908 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 43 \nStep 3: Ensemble-independent doublet detection 909 \nBenchmarking on computationally multiplexed pools with known ground-truth sample labels 910 \nrevealed that certain genetic demultiplexing tools, namely Demuxalot and Vireo, showed high 911 \ndoublet detection specificity, but that Steps 1 and 2 of the Ensemblex workflow failed to 912 \ncorrectly label a subset of doublet calls by these tools. To mitigate this issue and maximize the 913 \nrate of doublet identification, Ensemblex labels the cells that are identified as doublets by Vireo 914 \nor Demuxalot as double ts by default; however, users can nominate different tools for the 915 \nensemble-independent doublet detection component depending on the desired doublet 916 \ndetection stringency. Doublet specificity was computed using equation 6: 917 \n 918 \n(6) 𝐷𝑜𝑢𝑏𝑙𝑒𝑡 𝑠𝑝𝑒𝑐𝑖𝑓𝑖𝑐𝑖𝑡𝑦=  (\n𝑇𝑁\n𝑇𝑁+𝐹𝑃) 919 \n 920 \nWhere TN is the number of correctly classified doublets; FP is the number of true singlets 921 \nincorrectly classified as doublets. 922 \n 923 \nEnsemblex singlet assignment confidence score 924 \nEnsemblex computes a singlet confidence score to inform which cells should be discarded to 925 \navoid misclassification in downstream analyses.  First, Ensemblex evaluates how well an 926 \nindividual constituent tool’s assignment probability (e.g., Demuxalot) corresponded to the 927 \naccuracy of their assignment, using consensus cells across the three remaining tools ( e.g., 928 \nDemuxlet, Souporcell, Vireo) as a proxy for ground-truth, by fitting a binary logistic regression 929 \nmodel to compute the odds that a singlet was correctly classified given its corresponding 930 \nprobability. Using the binary logistic regression models, Ensemblex computes the AUC using 931 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 44 \nthe empirical method implemented in the ROCit (v2.1.1) R package for each tool (26). Then, 932 \nfor each cell, if Ensemblex’s sample label matche s that of a constituent tool, and if the 933 \nassignment probability of the constituent tool supersedes its probability threshold, the tool’s 934 \ncomputed AUC is added to the accuracy-weighted probabilistic ensemble probability produced 935 \nin Step 1 to yield the confidence score. By default, singlet assignments with a confidence score 936 \nless than 1.00 are labelled as unassigned by Ensemblex. Ensemblex’s confidence score and the 937 \ndesignated threshold is a successful predictor of accurately classified singlets because singlets 938 \nwill only achieve a confidence score ≥ 1 if: 939 \n1.  All constituent tools show the same sample label (accuracy -weighted probabilistic 940 \nensemble probability = 1.00); 941 \n2. At least one constituent tool confidently assigns the cell to an individual donor and the 942 \nconstituent tool’s probability assignment adequately corresponds to the overall 943 \naccuracy of their singlet assignment. 944 \n 945 \nApplication of Ensemblex with and without prior genotype information 946 \nGiven the dependencies of certain tools on prior genotype information, there are notable 947 \ndifferences between the Ensemblex workflows for demultiplexing with and without prior 948 \ngenotype information . When demultiplexing with prior genotype information, Ensemblex 949 \nleverages the sample labels from Demuxalot, Demuxlet, and Vireo -GT with prior genotype 950 \ninformation, and Souporcell without prior genotype information.  When demultiplexing 951 \nwithout prior genotype information, Ensemblex leverages the sample labels from Demuxalot, 952 \nFreemuxlet, Souporcell, and Vireo. However, given that Demuxalot requires prior genotype 953 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 45 \ninformation, Ensemblex uses the estimated donor .vcf file generated by Freemuxlet for input 954 \ninto the Demuxalot algorithm as prior genetic data.  955 \n 956 \nRunning the Ensemblex pipeline 957 \nA complete user guide for running the Ensemblex pipeline can be found at the Ensemblex 958 \nGitHub site: https://neurobioinfo.github.io/ensemblex/site/.  We provide two distinct yet highly 959 \ncomparable pipelines depending on the availability of prior genotype information . Both 960 \npipelines can be downloaded as a singularity image and are comprised of four steps: 961 \n1. Establish the pipeline and working directory; 962 \n2. Prepare input files for constituent genetic demultiplexing tools; 963 \n3. Parallel demultiplexing by constituent genetic demultiplexing tools; 964 \n4. Application of the Ensemblex algorithm for ensemble classification. 965 \n 966 \nAs input into the Ensemblex pipeline, users must provide a .tsv file describing the barcodes of 967 \nthe pooled cells, a. bam sequencing file for the pool, a reference genotype .vcf file (e.g., 1000 968 \nGenome Project) (27), a reference genome sequence .fasta file (e.g., 10X Genomics), and, if 969 \ndemultiplexing with prior genotype information, a .vcf  file describing the genetic data of the 970 \npooled samples.  971 \n 972 \nGenetic demultiplexing by constituent tools 973 \nGenetic demultiplexing by the constituent demultiplexing tools was performed following best 974 \npractices as defined by the authors of the respective tools using Python (v3.8.10).  975 \nDemuxalot  976 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 46 \nCellRanger-generated .bam file, filtered barcode .tsv file, and the corresponding donor .vcf file 977 \nwere used as input into the Demuxalot workflow. Candidate variants for scRNAseq genotyping 978 \nwere retained if the minimum coverage was > 200 and minimum alternative coverage was > 979 \n10. The top 100 SNPs per donor were retained to cluster the cells by genotype. Doublet calls 980 \nwere made with a prior strength of 0.25. 981 \n 982 \nDemuxlet  983 \nWe used the popscle  suite ( https://github.com/statgen/popscle) for Demuxlet. CellRanger-984 \ngenerated .bam file, filtered barcode .tsv file, and the corresponding donor .vcf file were used 985 \nas input into the Demuxlet workflow. The dsc-pileup function was first used to pileup candidate 986 \nvariants around known variant sites with the following parameters: --cp-BQ 40 --min-BQ 13 -987 \n-min-MQ 20 --minTD 0 --min-total 0 --min-uniq 0 --min-snp 0. The Demuxlet algorithm was 988 \nthen applied to cluster the cells by genotype with the following parameters: --geno-error-offset 989 \n0.10 --geno-error-coeff 0.00 --min-callrate 0.50 --doublet-prior 0.50 --cap-BQ 40 --min-BQ 13 990 \n--min-MQ 20 --min-TD 0 --min-total 0 --min-uniq 0 --min-snp 0. 991 \n 992 \nFreemuxlet  993 \nWe used the popscle suite ( https://github.com/statgen/popscle) for Freemuxlet. CellRanger-994 \ngenerated .bam file, filtered barcode .tsv file, and reference genotype .vcf file from the 1000 995 \nGenomes Project, phase 3 (27), were used as input into the Freemuxlet workflow.  The dsc-996 \npileup function was first used to pileup candidate variants around known variant sites with the 997 \nfollowing parameters: --cp-BQ 40 --min-BQ 13 --min-MQ 20 --minTD 0 --min-total 0 --min-998 \nuniq 0 --min-snp 0. The Freemuxlet algorithm was then applied to cluster the cells by genotype 999 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 47 \nwith the following parameters: --doublet-prior 0.50 --bf-thres 5.41 --frac-init-clust 0.50 --inter-1000 \ninit 10 --cap-BQ 40 --min-BQ 13 --min-total 0 --min-uniq 0 --min-snp 0. 1001 \n 1002 \nSouporcell  1003 \nCellRanger-generated .bam file, filtered barcode .tsv file, 10X Genomics reference .fasta file, 1004 \nand the corresponding donor .vcf file when demultiplexing with prior genotype information 1005 \nwere used as input into the Souporcell workflow. A FASTQ file was first generated from the 1006 \n.bam file using the renamer.py script. These reads were mapped to the reference genome using 1007 \nminimap2 with the following parameters: --ax splice –t 8 –G50k –k 21 –w 11 –sr --A2 –B8 –1008 \nO12,32 –E2,1 –r200 –p.5 –N20 –f1000,5000 –n2 –m20 –s40 –g200 –2k50m –secondary=no. 1009 \nThe barcodes and UMI were added back to the .sam file using the retag.py script and the 1010 \nresulting .bam file was sorted and indexed with Samtools. Variants were called using Freebayes 1011 \nwith the following parameters: --iXu –C 2 –q 20 –n 3 –E 1 –m 30 –min-coverage 6. Vartix was 1012 \nused to compute the number of alleles for each cell using the following parameters: --umi –1013 \nmapq 30 –scoring-method coverage. The Souporcell algorithm was then applied to cluster the 1014 \ncells by genotype; when demultiplexing with prior genotype information the --1015 \nknown_genotypes and --known_genotypes_sample_names parameters were included. 1016 \nTroublet was used to identify doublets and the consensus.py script was used for genotype and 1017 \nambient RNA co-inference.  1018 \n 1019 \nVireo 1020 \nCellRanger-generated .bam file, filtered barcode .tsv file, reference genotypes from the 1000 1021 \nGenomes Project, phase 3 (27), and the corresponding donor .vcf file when demultiplexing 1022 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 48 \nwith prior genotype information were used as input to the Vireo workflow. CellSNP was used 1023 \nto identify candidate variants for scRNAseq genotyping with the following parameters: --1024 \nminMAF 0.1 and --minCOUNT 100. The Vireo algorithm was then applied to cluster the cells 1025 \nby genotype with the --forceLearnGT parameter; when demultiplexing with prior genotype 1026 \ninformation (Vireo-GT) the --d and --t GT parameters were used. 1027 \n 1028 \nConsensus demultiplexing framework 1029 \nFor the consensus demultiplexing framework, singlets were considered confidently classified 1030 \nif Demuxalot, Demuxlet, Vireo, and Souporcell assigned a cell to the same donor -of-origin. 1031 \nCells classified as “ambiguous” or doublet by at least one tool were discarded. 1032 \n 1033 \nGeneration of computationally pooled samples for ground-truth benchmarking 1034 \nTo benchmark Ensemblex on computationally pooled samples with known ground -truth sample 1035 \nlabels, we leveraged 80 independently sequenced iPSC lines from Parkinson’s disease patients and 1036 \nhealthy controls, which were differentiated towards a dopaminergic neuronal state and sequenced 1037 \nafter 65 days of differentiation as part of the FOUNDIN-PD (14). Controlled access FASTQ files 1038 \nfrom the independently sequenced iPSC lines were  obtained from https://www.ppmi-info.org/ 1039 \n(accessed 09-17-2023) and processed by the CellRanger counts pipeline (v3.1.0)  with default 1040 \nparameters and aligned to GRCh38 reference genome. The CellRanger-generated .bam and filtered 1041 \nbarcode files were used as input into the synth_pool.py script produced by the authors of Vireo to 1042 \nsimulate sample pooling (9). In brief, reads from a subset of cells from the iPSC line-specific .bam 1043 \nfiles were merged and doublets were generated by combining the reads from random cell pairs. 1044 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 49 \nSample identities were added to each cell’s barcode, revealing the ground-truth sample labels for 1045 \nbenchmarking procedures. 1046 \n 1047 \nTo evaluate how genetic demultiplexing performance varied as a function of the number of 1048 \nmultiplexed samples, we generated 96 computationally multiplexed  pools using the 80 1049 \nFOUNDIN-PD lines with sample sizes of 4, 8, 16, 24, 32, 40,  48, 56, 64, 72, and 80. An equal 1050 \nnumber of cells from each line were used in the in silico pool. For the sample size of four we 1051 \ngenerated six replicates; for the sample sizes of 8-80 we generated nine replicates each. Replicates 1052 \nwere produced with different sample and cell combinations. The 96 in silico pools averaged 17,396 1053 \ncells (minimum = 8,696; maximum = 26,087). For this experiment, we maintained a 15% doublet 1054 \nrate as previously described (9).  1055 \n 1056 \nTo evaluate how genetic demultiplexing performance varied as a function of the number of cells 1057 \nin a pool, we generated 18 computationally multiplexed pools using the 80 FOUNDIN-PD lines 1058 \nwith 8,000, 16,000, 24,000, 32,000, 40,000, and 48,0000 pooled cells ; we generated three 1059 \nreplicates per pool size. Twenty-four samples were multiplexed for each pool and an equal number 1060 \nof cells from each sample were used. Replicates were produced with different sample and cell 1061 \ncombinations. For this experiment, we simulated a doublet rate of 6% per 8,000 pooled cells. 1062 \n 1063 \nTo evaluate if the overall demultiplexing performance varied due to the underrepresentation of a 1064 \ncell line, we generated 15 computationally multiplexed pools using the 80 FOUNDIN-PD lines 1065 \ncomprising 23 multiplexed samples with 1,000 cells and one randomly selected sample that 1066 \nshowed various degrees of underrepresentation, including 100 cells (10%), 300 cells (30%), 500 1067 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 50 \ncells (50%), 700 cells (70%), or 900 cells (90%). Three replicates were generated for each degree 1068 \nof underrepresentation. Replicates were produced with different sample and cell combinations. For 1069 \nthis experiment, we maintained a 18% doublet rate. 1070 \n 1071 \nWGS for the 80 donors from which the FOUNDIN-PD lines were derived was performed on whole 1072 \nblood-extracted DNA as previously described by the Parkinson’s Progression Markers Initiative  1073 \n(PPMI) (28). The c ontrolled-access WGS .vcf files were obtained from https://www.ppmi-1074 \ninfo.org/ (accessed 09-17-2023). Genotypes of common variants ( minor allele frequency > 5%) 1075 \nwere used as prior genotype information for the genetic demultiplexing tools in the benchmarking 1076 \nanalyses.   1077 \n 1078 \nPreparation, processing, and analysis of experimentally pooled samples 1079 \nUnless specified otherwise, e xperimentally pooled samples were processed with the CellRanger 1080 \ncounts pipeline (v5.0.1)  and analyzed with  the Seurat (v5.0.0) R package (29), using the 1081 \nscRNAbox analytical pipeline (30).  1082 \n 1083 \nNon-small cell lung cancer dataset 1084 \nNSCLC dissociated tumor cells  from seven donors were labelled with TotalSeq-B Human 1085 \nTBNK Cocktail (18). Multiplexed cells were then sequenced on an Illumina NovaSeq 6000 to 1086 \nan average read depth of approximately 70,000 reads per cell for gene expression and 25,000 1087 \nreads per cell for  CellPlex. Publicly available gene expression .bam and barcode  .tsv files 1088 \nreturned from the CellRanger multi pipeline (v6.1.2) were obtained from the 10X Genomics 1089 \nDatasets portal (10X Genomics Datasets) and used as input into the Ensemblex pipeline . We 1090 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 51 \nused the sample -specific gene expression .bam files and the BCFtools (v1.16) mpielup  1091 \nfunction to generate genotype likelihoods for prior genotype information (31).  1092 \n 1093 \nWe used HTOdemux to assign the cells back to their donor -of-origin based on the CMO 1094 \nexpression profiles as a proxy for ground-truth sample labels (19). Publicly available feature-1095 \nbarcode expression matrices were filtered to only include CMO labels used for multiplexing 1096 \n— CMO301, CMO302, CMO303, CMO304, CMO306, CMO307,  and CMO308  — and 1097 \nbarcodes with a CMO count > 0. The CMO expression profiles were normalized with Seurat’s 1098 \nNormalizeData function using the CLR normalization method and HTOdemux was applied to 1099 \nthe CMO assay using a positive quantile of 0.99.  1100 \n 1101 \nDopaminergic neuron dataset 1102 \nJerber et al. sequenced multiplexed experiments comprising 22 healthy donor iPSC lines from 1103 \nthe HipSci project (32) (http://www.hipsci.org) on days 11, 30, and 52 of DaN differentiation 1104 \nusing Illumina HiSeq 4000 to an average depth of 40,000-60,000 reads per cell (12). We used 1105 \nthree technical replicates for each timepoint, which are comprehensively described in 1106 \nAdditional File 1: Table S3. Publicly available gene expression .fastq files were obtained from 1107 \nthe European Nucleotide Archive (ENA) with accession number  ERP121676 and processed 1108 \nwith the CellRanger counts pipeline (v5.0.1) with default parameters  using the GRCh37  1109 \nreference genome. The CellRanger-generated. bam files, filtered barcode .tsv files, and .vcf 1110 \nfiles describing the pooled samples (see below) were used as input into the Ensemblex pipeline 1111 \nfor each technical replicate independently. Filtering of the scRNAseq data was performed as 1112 \ndescribed by Jerber et al.  (12). Genes with non-zero counts in at least 0.05% of cells  were 1113 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 52 \nretained. DoubletFinder (v2.0.4) was applied independently to each technical replicate. Time-1114 \npoint specific replicates were integrated with Seurat’s integration algorithm (33) and clustered 1115 \nby the Louvain network detection using the top 50 PCs and 10 nearest neighbours.  1116 \n 1117 \nWhole-exome sequencing (WES) .vcf files corresponding to the 22 pooled HipSci lines were 1118 \nobtained from the ENA with accession number PRJEB7243  (34). Genotypes of common 1119 \nvariants (minor allele frequency > 1%) were used as prior genotype information for the genetic 1120 \ndemultiplexing tools (12).  1121 \n 1122 \nNeural stem cell dataset  1123 \nWe performed two multiplexed experiments comprising iPSCs from individuals with ADHD 1124 \nand heathy controls differentiated into NSCs: Experiment 1 (n ADHD = 7; n control = 6) and 1125 \nExperiment 2 (n ADHD = 9; n control = 7). 1126 \n 1127 \nSubject recruitment 1128 \nPatients diagnosed with ADH D and matching healthy controls between  6−18 years old 1129 \nwere recruited by the Department of Child and Adolescent Psychiatry and Psychotherapy 1130 \nof the University of Zurich, as described previously (35). Inclusion and exclusion criteria 1131 \nfor recruitment of these individuals described previously (35). Additional File 1:  Table 1132 \nS4 provides a list of the individual subjects and their derived cell lines included in this 1133 \nstudy. Salivary DNA from ADHD patients and controls was genotyped using the Infinium 1134 \nGlobal Screening Array (Illumina) , as previously described, and used as prior genotype 1135 \ninformation for genetic demultiplexing (35). 1136 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 53 \n 1137 \nNeural stem cell culture 1138 \nThe generation and characterization of iPSC used in this study and the NSCs differentiation 1139 \nprotocols were previously described in (35) (36). NSCs cultures were seeded in two 1140 \nindependent experiments (designated as “1” and “2”), each of them consisting of  NSCs 1141 \npooled together into two culture dishes and maintained as NSCs until 100% confluence, 1142 \nwhen all iPSC lines were combined into one sample for sequencing. For most cell lines 1143 \ndifferent clones for each iPSC line were used in the two experiments  Additional File 1:  1144 \nTable S5. When applicable, the second clones of the same NSCs lines  were cultured 1145 \nseparately (designated as “.1” and “.2”)  in a second experiment. In the first experiment, 1146 \n56,250 cells per cell line were seeded in the pooled dishes.  In the second experiment the 1147 \nproportions of cells seeded we adjusted to their proliferation profile assessed in (36). Upon 1148 \nreaching 100% confluence, cells were dissociated for scRNAseq experiments and 1149 \ncombined to a single sample for sequencing as described below. 1150 \n 1151 \nDissociation of pooled neural stem cell cultures for single-cell RNA sequencing 1152 \nCells were washed in PBS and then incubated with 1 mL of StemPro Accutase (Gibco) for 1153 \n3 minutes at 37°C. After incubation, 2 mL of PBS, stopping the Accutase reaction, and cells 1154 \nwere gently pipetted up and down between 5 to 10 times to break up clumps of cells before 1155 \ntransfer to a 15 mL conical tube. The cells were centrifuged at 300 x g for 5 minutes and 1156 \nthe supernatant was removed. Following, 334 µL of Neural Expansion Media (NEM) was 1157 \nadded to each cell pellet using a 1000 µL pipette tip until cells were completely 1158 \nresuspended. An additional 666 µL of NEM was added to each well and gently pipette 1159 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 54 \nmixed 5 times. A 100 -µm cell strainer was used to filter the cell suspension  before 1160 \ncentrifugation at 300 x g for 4 minutes. The supernatant was carefully removed, and the 1161 \npellet was resuspended in 3 mL of PBS 1x containing 0.04% Bovine Serum Albumin 1162 \n(BSA) by pipetting up and down 5 times using a 5 mL serological pipette.  The cells were 1163 \ncentrifuged at 300 x g for 10 minutes and further submitted to live cell sorting with the 1164 \nMagnetic Dead Cell Removal Kit (Miltenyi Biotec, 130 -090-101), according to the 1165 \nmanufacturer. The resulting flow-through containing live cells was centrifuged for 300 x g 1166 \nfor 5 minutes and the supernatant was removed carefully to not disturb the cell pellet. Cells 1167 \nwere resuspended in 1 mL of PBS 1x containing 0.04% BSA for automated cell counting. 1168 \nFor each experiment, the cells from the two culture dishes were processed in parallel. Equal 1169 \ncounts of cells were combined for the final cell suspension for scRNAseq preparation at 1170 \nthe Functional Genomics Center Zurich at the University of Zurich. 1171 \n 1172 \nLibrary processing and sequencing 1173 \nAll samples were processed using the 10x Genomics  Chromium 3’ Single Cell Protocol  1174 \nand sequenced using NovaSeq 6000 S1 (Illumina).  For the first sample containing NSC 1175 \npools 1.1 and 1.2, 18,000 NSCs were loaded into one single 10x Genomics Lane to target 1176 \n13,000 cells. For the second sample containing NSC pools 2.1 and 2.2, 29,000 NSCs were 1177 \nloaded to target 18,000 cells.  1178 \n 1179 \n Demultiplexing and scRNAseq analysis  1180 \nFASTQ files were processed with the CellRanger counts pipeline ( v5.0.1) with default 1181 \nparameters and aligned to the GRCh37 reference genome. The CellRanger-generated. bam 1182 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 55 \nfiles, filtered barcode .tsv files, and .vcf files describing the pooled samples were used as 1183 \ninput into the Ensemblex pipeline. Genotypes of common variants (minor allele frequency 1184 \n> 1%) were used as prior genotype information for the genetic demultiplexing tools . The 1185 \nfiltered feature-barcode expression matrices  were used to analyz e the pooled cells  1186 \nfollowing a standard scRNAseq analysis workflow  using Seurat (30). Cells were filtered 1187 \nfor > 500 total and unique RNA transcripts. Doublets were removed using DoubletFinder 1188 \n(v2.0.4). The two NSC samples were integrated using Seurat’s integration algorithm (33). 1189 \nThe top 25 PCs were selected for Louvain network detection to identify clusters using 65 1190 \nnearest neighbours. Twelve clusters were identified at a clustering resolution of 0.25, which 1191 \nwere assigned as eight putative cell types using a combination of known markers and gene 1192 \nenrichment analysis. The top marker genes from each cluster were identified using Seurat’s 1193 \nFindAllMarkers with the Wilcoxon rank-sum test. Significant DEGs (log2 fold change > 1194 \n0.25 and P-value < 0.05 ) were input into EnrichR (37) and cell types were predicted with 1195 \nthe Cell Marker Augmented 2021 (38) and Azimuth Cell Types 2021 (39) libraries. Multiple 1196 \nclusters showed expression profiles for similar broad cell types — Neurons, NPC s, and 1197 \nNSCs. We used Seurat’s FindMarkers function to identify differentially expressed marker 1198 \ngenes between the clusters of the same broad cell type and top marker genes were selected 1199 \nto identify the cell subtypes.  1200 \n 1201 \nFor each putative cell type, DGE was calculated between ADHD and controls using the 1202 \nMAST statistical framework (22, 40) . Pooled cells were assigned as ADHD or control 1203 \nbased on the demultiplexed sample labels from each of the individual genetic 1204 \ndemultiplexing tools. Cells labeled as “ambiguous singlets” or doublets by the individual 1205 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 56 \ntools were excluded from their respective DGE analysis.  P-values were corrected for 1206 \nmultiple hypothesis testing using the Bonferroni method. A gene was considered 1207 \ndifferentially expressed if the adjusted P -value was ≤ 0.0 1 and the absolute value of the 1208 \nLog2 fold-change was ≥ 0.5. To compute DGE using the sample labels from the individual 1209 \ntools after the removal of Ensemblex’s putative doublet calls, we repeated the above 1210 \nprocedures but this time all cells labeled as doublets by the respective tool or Ensemblex 1211 \nwere excluded from the DGE analysis.  1212 \n 1213 \nPerformance metrics and statistical analyses 1214 \nWe performed all statistical analyses using the R statistical software (v4.2.2)  (41). We used the 1215 \nggplot2 R package (v3.4.2) for data visualization (42). 1216 \n 1217 \nSinglet classification 1218 \nA singlet was considered correctly classified if the demultiplexed sample label matched the 1219 \nground-truth sample label (i.e., specific sample ID) and the assignment probability exceeded 1220 \nthe recommended threshold for the respective tool. For computationally multiplexed pools, the 1221 \nproportion of correctly classified singlets was computed as: 1222 \n 1223 \n(7) 𝑃𝑟𝑜𝑝𝑜𝑟𝑡𝑖𝑜𝑛 𝑐𝑜𝑟𝑟𝑒𝑐𝑡 𝑠𝑖𝑛𝑔𝑙𝑒𝑡𝑠=  \n𝑇𝑃\n𝑛 𝑡𝑟𝑢𝑒 𝑠𝑖𝑛𝑔𝑙𝑒𝑡𝑠 1224 \n 1225 \nFor the NSCLC dataset , HTOdemux’s sample labels were considered ground-truth, and the 1226 \nsinglet TP and FP rate were computed as: 1227 \n 1228 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 57 \n(8) 𝑆𝑖𝑛𝑔𝑙𝑒𝑡 𝑇𝑃 𝑟𝑎𝑡𝑒 =  \n𝑇𝑃\n𝑛 𝐻𝑇𝑂𝑑𝑒𝑚𝑢𝑥 𝑠𝑖𝑛𝑔𝑙𝑒𝑡𝑠 1229 \n(9) 𝑆𝑖𝑛𝑔𝑙𝑒𝑡 𝐹𝑃 𝑟𝑎𝑡𝑒=  \n𝐹𝑃\n𝑛 𝐻𝑇𝑂𝑑𝑒𝑚𝑢𝑥 𝑠𝑖𝑛𝑔𝑙𝑒𝑡𝑠 1230 \n 1231 \nDoublet identification 1232 \nA doublet was considered correctly classified if the demultiplexed sample label matched the 1233 \nground-truth sample label, independent of the assignment probability. For computationally 1234 \nmultiplexed pools, the proportion of correctly classified doublets was computed as: 1235 \n 1236 \n(10) 𝑃𝑟𝑜𝑝𝑜𝑟𝑡𝑖𝑜𝑛 𝑐𝑜𝑟𝑟𝑒𝑐𝑡 𝑑𝑜𝑢𝑏𝑙𝑒𝑡𝑠=  \n𝑇𝑁\n𝑛 𝑡𝑟𝑢𝑒 𝑑𝑜𝑢𝑏𝑙𝑒𝑡𝑠 1237 \n 1238 \nFor the NSCLC dataset, TP doublets were defined as cells classified as doublets by both 1239 \nHTOdemux and Ensemblex; FP doublets were defined as cells classified as singlets by 1240 \nHTOdemux and doublets by Ensemblex; FN doublets were defined as cells classified as 1241 \ndoublets by HTOdemux and singlets by Ensemblex. The doublet TP, FP, and FN  rates were 1242 \ncomputed as: 1243 \n 1244 \n(11) 𝐷𝑜𝑢𝑏𝑙𝑒𝑡 𝑇𝑃 𝑟𝑎𝑡𝑒 =  \n𝑇𝑃\n𝑛 𝐻𝑇𝑂𝑑𝑒𝑚𝑢𝑥 𝑑𝑜𝑢𝑏𝑙𝑒𝑡𝑠 1245 \n(12) 𝐷𝑜𝑢𝑏𝑙𝑒𝑡 𝐹𝑃 𝑟𝑎𝑡𝑒=  \n𝐹𝑃\n𝑛 𝑝𝑜𝑜𝑙𝑒𝑑 𝑑𝑟𝑜𝑝𝑙𝑒𝑡𝑠 1246 \n(13) 𝐷𝑜𝑢𝑏𝑙𝑒𝑡 𝐹𝑁 𝑟𝑎𝑡𝑒= 1 −  𝐷𝑜𝑢𝑏𝑙𝑒𝑡 𝑇𝑃 𝑟𝑎𝑡𝑒 1247 \n 1248 \nAdjusted Rand Index 1249 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 58 \nTo evaluate the similarity between two distinct sample clusterings we computed the ARI using 1250 \nthe pdfCluster (v1.0.4) R package (43). For the benchmarking analyses, we computed the ARI 1251 \nbetween the demultiplexed sample labels by each genetic demultiplexing tool and the ground-1252 \ntruth sample labels (computationally pooled samples) or HTOdemux’s sample labels (NSCLC 1253 \ndataset). We followed the same procedure when computing the ARI between Ensemblex’s 1254 \nsample labels and those of its constituent tools (DaN and NSC datasets); however, the ground-1255 \ntruth sample labels were replaced by Ensemblex’s sample labels for these analyses. For 1256 \nexperiments evaluating the impact  of doublets on the stability of clusters in gene expression 1257 \nspace, we computed the ARI between clusters at a given  clustering resolution after removing 1258 \ndoublets identified by each genetic demultiplexing tool. Clustering stability was computed at 1259 \nresolutions of 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, and 1.0. For each clustering 1260 \nresolution, 25 iterations of Louvain clustering were performed while shuffling the order of the 1261 \nnodes in the graph. The ARI between clustering pairs at each clustering resolution was then 1262 \ncomputed.  1263 \n 1264 \nBalanced accuracy 1265 \nBalanced accuracies were computed to evaluate the binary classification performance of each 1266 \ngenetic demultiplexing tool  on imbalanced datasets, where doublets represented a minority 1267 \nclass compared to singlets.  The balanced accuracy of each genetic demultiplexing tool was 1268 \ncomputed against the ground -truth sample labels  (computationally pooled samples ) or 1269 \nHTOdemux’s sample labels (NSCLC dataset) using equation 1.  1270 \n 1271 \nMatthew’s correlation coefficient (MCC) 1272 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 59 \nThe MCC was used as a second metric for evaluating the binary classification performance of 1273 \nthe genetic demultiplexing tool. The MCC of each genetic demultiplexing tool was computed 1274 \nagainst the ground -truth sample labels (computationally pooled samples ) or HTOdemux’s 1275 \nsample labels (NSCLC dataset) using equation 14: 1276 \n 1277 \n(14) 𝑀𝐶𝐶 =  \n𝑇𝑁×𝑇𝑃−𝐹𝑁×𝐹𝑃\n√(𝑇𝑃+𝐹𝑃)(𝑇𝑃−𝐹𝑁)(𝑇𝑁+𝐹𝑃)(𝑇𝑁+𝐹𝑁) 1278 \n 1279 \nArea under the receiver operating characteristic curve for singlet detection 1280 \nTo evaluate how well each genetic demultiplexing tool’s assignment probability corresponded 1281 \nto the accuracy of their singlet assignments when ground-truth sample labels were known, we 1282 \nfit a binary logistic regression model to compute the odds that a singlet was correctly classified 1283 \nby a tool given the corresponding confidence score or probability. Correctly and incorrectly 1284 \nclassified singlets were set as the positive and negative references, respectively. We then used 1285 \nthe binary logistic regression model to compute the receiver operating characteristic curve for 1286 \neach tool , which plots the singlet TP and FP rate s across classification thresholds,  and 1287 \ncalculated the AUC using the empirical method implemented in the ROCit (v2.1.1) R package 1288 \n(26).  1289 \n 1290 \nAbbreviations 1291 \nADHD, attention deficit hyperactivity disorder; ANOV A, Analysis of variance; ARI, Adjusted 1292 \nRand Index; AUC, area under the receiver operating characteristic curve; BSA, Bovine Serum 1293 \nAlbumin; CMO, Cell Multiplexing Oligonucleotides; DaN, dopaminergic neurons; DGE, 1294 \ndifferential gene expression;  DEG, differentially expressed genes;  ENA, European Nucleotide 1295 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 60 \nArchive; eQTL, expression quantitative trait loci;  FN, false -negative; fNN, nearest neighbour 1296 \nfrequency; FOUNDIN -PD; Foundational Data Initiative for Parkinson’s Disease; FP, false 1297 \npositive; iPSC, induced pluripotent stem cell; K, kurtosis; MAST, model-based analysis of single-1298 \ncell transcriptomics;  MCC, Matthew’s Correlation Coefficient; nCD, number of confident 1299 \ndoublets; NEM, neural expansion media; NPC, neural progenitor cell; NSC, neural stem cell; 1300 \nNSCLC, non-small cell lung cancer; PC, principal component; PCA principal component analysis; 1301 \nPPMI, Parkinson’s Progression Markers Initiative; pT, nearest neighbour percentile threshold; 1302 \nscATACseq, single-cell assay for  transposase-accessible chromatin sequencing;  scRNAseq, 1303 \nsingle-cell RNA sequencing;  SNP, single nucleotide polymorphism; snRNAseq, single -nuclei 1304 \nRNA sequencing;  TN, true-negative; TP, true-positive; UMI, unique molecular identified; WES, 1305 \nwhole-exome sequencing; WGS, whole-genome sequencing.  1306 \n 1307 \nDeclarations 1308 \nEthics approval and consent to participate 1309 \nThe iPSC lines (ADHD & controls) used in this  project were approved by the Cantonal Ethics 1310 \nCommittee Zurich (BASEC-Nr.-2016-00101 & BASEC -Nr.-201700825) and followed the latest 1311 \nversion of the Declaration of Helsinki, as previously reported  (35). The subjects and/or parents 1312 \nhave voluntarily consented to participate in this study. 1313 \n 1314 \nConsent for publication 1315 \nNot applicable. 1316 \n 1317 \nAvailability of data and materials 1318 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 61 \nTranscriptional data for the 80 independently sequenced iPSC lines and the corresponding WGS 1319 \ndata are available from the PPMI database (www.ppmi-info.org/access-dataspecimens/download-1320 \ndata), RRID:SCR 006431. For up -to-date information on the study, visit www.ppmi-info.org. 1321 \nProcessed transcriptional data for the NSCLC dataset are available from the 10X Genomics 1322 \nDatasets Portal  (https://www.10xgenomics.com/datasets/20k-mixture-of-nsclc-dtcs-from-7-1323 \ndonors-3-v3-1-with-intronic-reads-3-1-standard). Transcriptional data for the DaN datasets are 1324 \navailable from the ENA with accession number ERP121676. WES data for the 22 HipSci lines 1325 \npooled in the DaN datasets are available from the ENA with accession number PRJEB7243. 1326 \nProcessed scRNAseq data for the NSC dataset are available  from the corresponding author upon 1327 \nreasonable request. The code used for the analyses presented in the work are available at  1328 \nhttps://github.com/neurobioinfo/ensemblex. Ensemblex is freely available under an MIT open -1329 \nsource license at https://zenodo.org/records/11639103.  1330 \n 1331 \nCompeting interests 1332 \nThe authors declare that they have no competing interests. 1333 \n 1334 \nFunding 1335 \nThis work was supported by the Michael J. Fox Foundation  [MJFF-021629 to EAF, RAT, and 1336 \nSMKF]. PPMI — a public-private partnership — is funded by the Michael J. Fox Foundation for 1337 \nParkinson’s Research and funding partners, including 4D Pharma, Abbvie, AcureX, Allergan, 1338 \nAmathus Therapeutics, Aligning Science Across Parkinson's, AskBio, Avid Radiopharmaceuticals, 1339 \nBIAL, BioArctic, Biogen, Biohaven, BioLegend, BlueRock Therapeutics, Bristol -Myers Squibb, 1340 \nCalico Labs, Capsida Biotherapeutics, Celgene, Cerevel Therapeutics, Coave Therapeutics, 1341 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 62 \nDaCapo Brainscience, Denali, Edmond J. Safra Foundation, Eli Lilly, Gain Therapeutics, GE 1342 \nHealthCare, Genentech, GSK, Golub Capital, Handl Therapeutics, Insitro, Jazz Pharmaceuticals, 1343 \nJohnson & Johnson Innovative Medicine, Lundbeck, Merck, Meso Scale Discovery, Mission 1344 \nTherapeutics, Neurocrine Biosciences, Neuron23, Neuropore, Pfizer, Piramal, Prevail 1345 \nTherapeutics, Roche, Sanofi, Servier, Sun Pharma Advanced Research Company, Takeda, Teva, 1346 \nUCB, Vanqua Bio, Verily, V oyager v. 25MAR2024 Therapeutics, the Weston Family Foundation 1347 \nand Yumanity Therapeutics. For funding t he ADHD study, we thank the Neuroscience Centre 1348 \nZurich (ZNZ) for the Zurich -McGill University Neurodevelopmental Disorder Research 1349 \nCollaboration and the  Psychiatric University Hospital Zurich  (PUK) Forschungsfonds Nr. 8702 1350 \n“Fonds für wissenschaftliche Zwecke im Interesse der Heilung von psychiatrischen Krankheiten” 1351 \nand the Candoc PhD grant from the University of Zurich [FK-22-044 to CMYO]. 1352 \n 1353 \nAuthors’ contributions 1354 \nMRF, EAF, RAT, and SMKF conceived the study. MRF developed the Ensemblex framework and 1355 \nwrote the corresponding R code. MRF performed the analyses and produced the figures. MRF and 1356 \nSA developed the Ensemblex pipeline and created the GitHub site. MRF  and SA tested the 1357 \nEnsemblex pipeline. MRF wrote the Ensemblex documentation. MRF,  AAD, RAT and SMKF 1358 \ninterpreted the data sets. CMYO performed all cell cultures  and sequencing preparation for the 1359 \nNSC dataset. MRF, CMYO, and RAT performed the cell type annotations for the NSC dataset. LS 1360 \nand EG provided the NSC genetic data. SW recruited the subjects for the NSC dataset. MRF wrote 1361 \nthe manuscript with input from all authors. EG  supervised the NSC data collection. RAT and 1362 \nSMKF supervised the project.  1363 \n 1364 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 63 \nAcknowledgments 1365 \nThe author’s acknowledge Dan Spiegelman for his help in processing the VCF files for individuals 1366 \npooled in the NSC dataset. Schematic illustrations presented in the manuscript were prepared with 1367 \nBioRender (https://www.biorender.com/). 1368 \n 1369 \nAuthors’ information 1370 \nThis work was supported by the Michael J Fox Foundation in a grant award to EAF , RAT, and 1371 \nSMKF [MJFF-021629]. MRF is supported by a CIHR Canada Graduate Scholarships -Master’s 1372 \nAward, a Fonds de Recherche Santé Québec Master’s Award, and a Brain Canada Rising Stars 1373 \nAward. EAF is supported by a Fonds d’Accéleration des Collaborations en Santé (FACS) grant 1374 \nfrom CQDM/MEI and a Canada Research Chair (Tier 1) in Parkinson’s disease.  R.A.T. received 1375 \nfunding through the McGill Healthy Brains for Healthy Lives (HBHL) Postdoctoral Fellowship 1376 \nand Molson NeuroEngineering Fellowship.  SMKF received funding from Brain Canada and the 1377 \nMontreal Neurological Institute-Hospital. CMYO is supported by the Candoc PhD grant from the 1378 \nUniversity of Zurich (UZH) [FK-22-044]. 1379 \n 1380 \nFigure legends   1381 \nFigure 1. Evaluation of existing individual genetic demultiplexing tools.  Evaluation of genetic 1382 \ndemultiplexing tools with prior genotype information  on 96 in silico pools with known ground -1383 \ntruth sample labels ranging in size from 4 to 80 multiplexed induced pluripotent stem cell (iPSC) 1384 \nlines from genetically distinct individuals, averaging 17,396 cells per pool and a 15% doublet rate. 1385 \nA) Line graphs showing the proportion of correctly classified singlets, doublets, and all cells by 1386 \neach individual genetic demultiplexing tool across varying numbers of multiplexed iPSC lines in 1387 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 64 \na single pool (sample number). The large dots show the mean proportion of correct classifications 1388 \nby an individual tool across replicates at a given sample size (n = 9 per pool size). The blue points 1389 \nshow the proportion of cells that were correctly classified by at least one individual genetic 1390 \ndemultiplexing tool: Demuxalot, Demuxlet, Souporcell, or Vireo -GT. B) Bar chart showing the  1391 \nmean proportion of total cells from an individual pool  correctly classified by only one genetic 1392 \ndemultiplexing tool. Error bars represent one standard deviation from the mean. (n = 9 per pool 1393 \nsize) C) Bar chart showing the  proportion of correctly classified singlet cells labelled as 1394 \n“unassigned” (ambiguous singlet assignments)  due to assignment probabilities below the 1395 \nrecommended threshold of the respective genetic demultiplexing tool. Error bars represent one 1396 \nstandard deviation from the mean. (n = 9 per pool size). 1397 \n 1398 \nFigure 2. Characterization of the Ensemblex framework. Ensemblex is a probabilistic -1399 \nweighted ensemble genetic demultiplexing framework for s ingle-cell RNA sequencing analysis, 1400 \nwhich was designed to leverage the most probable sample labels from each of its constituent tools: 1401 \nDemuxalot, Demuxlet, Souporcell, and Vireo  when using prior genotype information or 1402 \nDemuxalot, Freemuxlet, Souporcell, and Vireo when prior genotype information is not available. 1403 \nA) The Ensemblex workflow begins with demultiplexing pooled cells from genetically distinct 1404 \nindividuals by each of the constituent tools. The outputs from each individual demultiplexing tool 1405 \nare then used as input into the Ensemblex framework. B) The Ensemblex framework comprises 1406 \nthree distinct steps that are assembled into a pipeline: 1) accuracy-weighted probabilistic ensemble, 1407 \n2) graph -based doublet detection, and 3) ensemble -independent doublet detection. C-D) Line 1408 \ngraphs showng t he contribution of each step of the Ensemblex framework on 96 in silico pools 1409 \nwith known ground -truth sample labels ranging in size from 4 to 80 multiplexed induced 1410 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 65 \npluripotent stem cell (iPSC) lines from genetically distinct individuals, averaging 17,396 cells per 1411 \npool and a 15% doublet rate. The average proportion of correctly classified C) singlets and D) 1412 \ndoublets across replicates at a given pool size is shown after sequentially applying each step of the 1413 \nEnsemblex framework with prior genotype information (n = 9 per pool size) . The right panels 1414 \nshow the average proportion of correct classifications across all 96 pools; error bars represent one 1415 \nstandard deviation from the mean. The blue points show the proportion of cells that were correctly 1416 \nclassified by at least one individual genetic demultiplexing tool: Demuxalot, Demuxlet, 1417 \nSouporcell, or Vireo-GT.  1418 \n 1419 \nFigure 3. Ensemblex ground-truth benchmarking on computationally multiplexed pools. The 1420 \ngenetic demultiplexing tools with prior genotype information were evaluated on 96 in silico pools 1421 \nwith known ground -truth sample labels ranging in size from 4 to 80 multiplexed induced 1422 \npluripotent stem cell (iPSC) lines from genetically distinct individuals, averaging 17,396 cells per 1423 \npool and a 15% doublet rate. A singlet was considered correctly classified if the assigned sample 1424 \nlabel matched the ground -truth sample label and the assignment probability exceeded the 1425 \nrecommended threshold for the respective tool; a doublet was considered correctly classified if the 1426 \nassigned sample label matched the ground -truth sample label, regardless of the assignment 1427 \nprobability. A-I) Line graphs showing the performance of Ensemblex and the individual genetic 1428 \ndemultiplexing tools across evaluation metrics. The large dots show the mean value for each tool 1429 \nacross replicates at a given sample size (n = 9 per pool size). A) Proportion of correctly classified 1430 \nsinglets. B) Proportion of correctly classified doublets. C) Proportion of correctly classified cells. 1431 \nD) Adjusted Rand Index between each tool’s sample labels and the ground-truth sample labels. E) 1432 \nBalanced accuracy of each tool. F) Matthew’s Correlation Coefficient of each tool. G) Area under 1433 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 66 \nthe receiver operating characteristic curve (AUC) of the singlet assignment probability for each 1434 \ntool. H) Proportion of usable cells returned by each tool. Usable cells were defined as  cells 1435 \nclassified by singlets with an assignment probability exceeding the recommended threshold of the 1436 \nrespective tool. I) Error rate amongst the usable cells returned by each tool; erroneous 1437 \nclassifications comprised of true doublets labeled as singlets or true singlets assigned to the wrong 1438 \nsample. 1439 \n 1440 \nFigure 4. Evaluating Ensemblex on experimentally multiplexed cells  using donor -specific 1441 \noligonucleotide labels as a proxy for ground -truth. Non-small cell lung cancer  (NSCLC) 1442 \ndissociated tumor cells from 7 individuals were pooled and labelled with donor-specific 1443 \noligonucleotide-labels. Cells were demultiplexed according to their expression of donor-specific 1444 \noligonucleotide labels by HTOdemux ; HTOdemux’s sample labels were used as a proxy for 1445 \nground truth.  True positives (TP)  singlets were defined as cells classified as singlets by both 1446 \nHTOdemux and Ensemblex with matching sample labels; false positives (FP)  singlets were 1447 \ndefined as cells classified as singlets by both HTOdemux and Ensemblex but assigned to different 1448 \ndonors. TP doublets were defined as cells classified as doublets by both HTOdemux and 1449 \nEnsemblex; FP doublets were defined as cells classified as singlets by HTOdemux and doublets 1450 \nby Ensemblex; false negatives (FN)  doublets were defined as cells classified as doublets by 1451 \nHTOdemux and singlets by Ensemblex. A) T-distributed Stochastic Neighbor Embedding (t-SNE) 1452 \nvisualization of HTOdemux’s sample labels. B) T-SNE visualization of  Ensemblex’s 1453 \ndemultiplexing performance using HTOdemux’s sample labels as ground truth  for singlets (left) 1454 \nand doublets (right) . C) Bar plots showing  the s inglet TP and FP rate s for each  genetic 1455 \ndemultiplexing tool using HTOdemux’s sample labels as ground truth . D) Bar plots showing the 1456 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 67 \ndoublet TP and FP rates for each genetic demultiplexing tool using HTOdemux’s sample labels as 1457 \nground truth. E) Scatter plot showing the proportion of usable cells (confidently classified singlets) 1458 \nand the corresponding usable cell error rate for each genetic demultiplexing tool. F) Adjusted Rand 1459 \nIndex, balanced accuracy, Matthew’s Correlation Coefficient, and area under the receiver operating 1460 \ncharacteristic curve (AUC) of the singlet assignment probability for each genetic demultiplexing 1461 \ntool. 1462 \n 1463 \nFigure 5. Application of Ensemblex to highly multiplexed, experimentally pooled cultures of 1464 \ndifferentiated dopaminergic neurons.  A) Time line of iPSC pooling, dopaminergic neuron 1465 \n(DaN) differentiation, and sample collection from the DaN dataset by Jerber et al. (12). Three 1466 \ntechnical replicates at each time point (days 11, 30 and, 52 of differentiation) from pools containing 1467 \n22 individual iPSC lines were used in the analysis. Across all timepoints and technical replicates, 1468 \n84,746 cells were obtained for analysis.  B) Uniform manifold approximation and projection 1469 \n(UMAP) plots showing confidently assigned singlets or predicted doublets (blue) and ambiguous 1470 \nsinglets (singlet assignments with insufficient assignment probabilities ; red ) returned by each 1471 \ndemultiplexing tool. C) Stacked bar chart showing the proportion of confidently assigned singlets 1472 \nor predicted doublets (blue) and ambiguous singlets (red) across technical replicates at each time 1473 \npoint returned by each demultiplexing tool.  D) Boxplot showing the p roportion of confidently 1474 \nclassified singlets across technical replicates  and time points  by each demultiplexing tool. 1475 \nWilcoxon rank-sum tests were used to compare the proportion of confidently classified singlets by 1476 \nEnsemblex to that of its constituents  (n = 9  pools). E) Bar chart showing the  proportion of 1477 \noverlapping ambiguous singlet assignments  amongst demultiplexing tools across technical 1478 \nreplicates and time points (n = 9 pools) . F) Boxplot showing the Adjusted Rand Index (ARI) 1479 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 68 \nassessing cluster stability across a range of 11 clustering resolutions (n clustering iterations = 25) 1480 \nafter removing doublets identified by each demultiplexing tool. Wilcoxon rank -sum tests were 1481 \nused to compare the clustering ARI after removing Ensemblex doublets to the clustering ARI after 1482 \nremoving doublets identified by each constituent tool. * Adjusted P-value < 0.05; ** adjusted P-1483 \nvalue < 0.01; *** adjusted P-value < 0.001 1484 \n 1485 \nFigure 6. Evaluating the impact of discordant assignments between genetic demultiplexing 1486 \ntools on differential gene expression  analysis. A) Schematic illustrating the workflow for the 1487 \nneural stem cell (NSC) dataset. Pooled induced pluripotent stem cell (iPSC) -derived neural stem 1488 \ncell cultures from individuals with attention deficit hyperactivity disorder (ADHD) and controls 1489 \nwere collected in two separate experiments . NSCs were dissociated for single-cell RNA 1490 \nsequencing and prior genotype information of the pooled subjects  was obtained through  1491 \nmicroarray genotyping. The pools were demultiplexed by Ensemblex and its constituent s with 1492 \nprior genotype information and differential gene expression (DEG) was computed between ADHD 1493 \nand controls. B) Uniform manifold approximation and projection (UMAP) plot showing the 1494 \nputative cell types. C) Summary of the number of usable cells — singlets above the recommended 1495 \nprobability threshold of the respective demultiplexing tool — assigned to ADHD donors and 1496 \ncontrols and the number of identified doublets by each demultiplexing tool. D) Boxplot showing 1497 \nthe Adjusted Rand Index (ARI) assessing cluster stability across a range of 11 clustering 1498 \nresolutions (n clustering iterations = 25) after removing doublets identified by each demultiplexing 1499 \ntool. A one-way Analysis of Variance (ANOV A) test comparing the ARI after removing doublets 1500 \nidentified by each tool revealed a significant difference between tools (n = 11 clustering 1501 \nresolutions; P-value = 1.18e-3). E) Proportion of ADHD and control cells identified as putative 1502 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 69 \ndoublets by Ensemblex that were assigned as singlets by the constituent demultiplexing tools. F) 1503 \nHeatmap showing the number of cell-type specific DEGs between ADHD and controls using the 1504 \nsubject labels of each demultiplexing tool. G) Heatmap showing the number of cell-type specific 1505 \nDEGs between ADHD and controls  using the subject labels of each demultiplexing tool  and 1506 \nremoving putative doublets identified by Ensemblex. Cell-types not shown in the heatmaps had no 1507 \nDEGs passing the adjusted P-value < 0.01 and |Log2FC >= 0.5| threshold across all tools. 1508 \n 1509 \nTables 1510 \nTable 1. Summary of individual genetic demultiplexing tools.  1511 \nGenetic demultiplexing tool \nPrior genotype information for \ngenetic demultiplexing \nIncluded in the Ensemblex \nframework \nDemuxalot (5) Required Yes \nDemuxlet (6) Required Yes \nFreemuxlet (6) Not supported Yes \nScSplit (7) Optional No \nSouporcell (8) Optional Yes \nVireo (9) Optional Yes \n 1512 \nTable 2. Application of Ensemblex to pooled cultures of dopaminergic neurons from 22 1513 \nhealthy controls.  1514 \n ARI between Ensemblex and \nconstituent tool assignments \n Percent contribution to \nEnsemblex assignments \nn \nusable cells \nn \ndoublets  Day 11 Day 30 Day 52  Day 11 Day 30 Day 52 \nDemuxalot 0.987 0.955 0.982  97.29% 94.75% 97.57% 75,962 8,279 \nDemuxlet 0.928 0.062 0.884  95.91% 29.74% 90.55% 57,567 6,614 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 70 \nSouporcell 0.883 0.876 0.912  91.62% 91.82% 93.84% 76,811 7,740 \nVireo-GT 0.961 0.879 0.958  95.95% 88.80% 95.16% 75,933 6,115 \nEnsemblex NA NA NA  NA NA NA 76,222 8,307 \nDoubletFinder NA NA NA  NA NA NA NA 4,597 \nPooled cultures of induced pluripotent stem cell (iPSC) lines from 22 healthy donors  were 1515 \ndifferentiated towards a dopaminergic neuron (DaN) fate and sequenced on days 11, 30, and 52 of 1516 \ndifferentiation by Jerber et al. (12). For the analysis we used three technical replicates for each 1517 \nsequencing timepoint. Each pool was demultiplexed independently by Ensemblex and its 1518 \nconstituent tools with prior genotype information. The Adjusted Rand Index (ARI)  between 1519 \nEnsemblex’s assignments and th ose of the  constituent tools was computed across technical 1520 \nreplicates corresponding to each differentiation timepoint. The percent contribution represents the 1521 \nproportion of assignments from each constituent tool that matched Ensemblex’s assignments. 1522 \nUsable cells were defined as singlet classifications whose assignment probability exceeded the 1523 \nrecommended threshold of the respective tool. Abbreviations: NA = Not applicable. 1524 \n 1525 \nReferences 1526 \n1. Fiorini MR, Dilliott AA, Thomas RA, Farhan SMK. Transcriptomics of Human Brain Tissue 1527 \nin Parkinson's Disease: a Comparison of Bulk and Single-cell RNA Sequencing. Mol Neurobiol. 1528 \n2024. 1529 \n2. Ringman JM, Goate A, Masters CL, Cairns NJ, Danek A, Graff-Radford N, et al. Genetic 1530 \nheterogeneity in Alzheimer disease and implications for treatment strategies. Curr Neurol 1531 \nNeurosci Rep. 2014;14(11):499. 1532 \n3. McKinney CE. Using induced pluripotent stem cells derived neurons to model brain 1533 \ndiseases. Neural Regen Res. 2017;12(7):1062-7. 1534 \n4. Howitt G, Feng Y , Tobar L, Vassiliadis D, Hickey P , Dawson MA, et al. Benchmarking single-1535 \ncell hashtag oligo demultiplexing methods. NAR Genomics and Bioinformatics. 1536 \n2023;5(4):lqad086. 1537 \n5. Rogozhnikov A, Ramkumar P , Shah K, Bedi R, Kato S, Escola GS. Demuxalot: scaled up 1538 \ngenetic demultiplexing for single-cell sequencing. bioRxiv. 2021:2021.05. 22.443646. 1539 \n6. Kang HM, Subramaniam M, Targ S, Nguyen M, Maliskova L, McCarthy E, et al. 1540 \nMultiplexed droplet single-cell RNA-sequencing using natural genetic variation. Nat Biotechnol. 1541 \n2018;36(1):89-94. 1542 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 71 \n7. Xu J, Falconer C, Nguyen Q, Crawford J, McKinnon BD, Mortlock S, et al. Genotype-free 1543 \ndemultiplexing of pooled single-cell RNA-seq. Genome Biol. 2019;20(1):290. 1544 \n8. Heaton H, Talman AM, Knights A, Imaz M, Gaffney DJ, Durbin R, et al. Souporcell: robust 1545 \nclustering of single-cell RNA-seq data by genotype without reference genotypes. Nat Methods. 1546 \n2020;17(6):615-20. 1547 \n9. Huang Y , McCarthy DJ, Stegle O. Vireo: Bayesian demultiplexing of pooled single-cell 1548 \nRNA-seq data without genotype reference. Genome Biol. 2019;20(1):273. 1549 \n10. Neavin D, Senabouth A, Hang Lee JT , Ripoll A, Consortium s-e, Franke L, et al. Demuxafy: 1550 \nImprovement in droplet assignment by integrating multiple single-cell demultiplexing and 1551 \ndoublet detection methods. BioRxiv. 2022:2022.03. 07.483367. 1552 \n11. Cardiello JF, Joven Araus A, Giatrellis S, Helsens C, Simon A, Leigh ND. Evaluation of 1553 \ngenetic demultiplexing of single-cell sequencing data from model species. Life Sci Alliance. 1554 \n2023;6(8). 1555 \n12. Jerber J, Seaton DD, Cuomo AS, Kumasaka N, Haldane J, Steer J, et al. Population-scale 1556 \nsingle-cell RNA-seq profiling across dopaminergic neuron differentiation. Nature genetics. 1557 \n2021;53(3):304-12. 1558 \n13. Jindal A, Gupta P , Jayadeva, Sengupta D. Discovery of rare cells from voluminous single 1559 \ncell expression data. Nature communications. 2018;9(1):4719. 1560 \n14. Bressan E, Reed X, Bansal V, Hutchins E, Cobb MM, Webb MG, et al. The Foundational 1561 \nData Initiative for Parkinson Disease: Enabling efficient translation from genetic maps to 1562 \nmechanism. Cell Genom. 2023;3(3):100261. 1563 \n15. Cardiello JF, Araus AJ, Giatrellis S, Helsens C, Simon A, Leigh ND. Evaluation of genetic 1564 \ndemultiplexing of single-cell sequencing data from model species. Life Science Alliance. 1565 \n2023;6(8). 1566 \n16. Large J, Lines J, Bagnall A. A probabilistic classifier ensemble weighting scheme based on 1567 \ncross-validated accuracy estimates. Data Min Knowl Discov. 2019;33(6):1674-709. 1568 \n17. Chicco D, Tötsch N, Jurman G. The Matthews correlation coefficient (MCC) is more 1569 \nreliable than balanced accuracy, bookmaker informedness, and markedness in two-class 1570 \nconfusion matrix evaluation. BioData mining. 2021;14:1-22. 1571 \n18. 20k Mixture of NSCLC DTCs from 7 donors, 3' v3.1 (with intronic reads) [Internet]. 10X 1572 \nGenomics. 2022 [cited January 8th, 2024]. Available from: 1573 \nhttps://www.10xgenomics.com/datasets/20k-mixture-of-nsclc-dtcs-from-7-donors-3-v3-1-with-1574 \nintronic-reads-3-1-standard. 1575 \n19. Stoeckius M, Zheng S, Houck-Loomis B, Hao S, Yeung BZ, Mauck WM, 3rd, et al. Cell 1576 \nHashing with barcoded antibodies enables multiplexing and doublet detection for single cell 1577 \ngenomics. Genome Biol. 2018;19(1):224. 1578 \n20. McGinnis CS, Murrow LM, Gartner ZJ. DoubletFinder: Doublet Detection in Single-Cell 1579 \nRNA Sequencing Data Using Artificial Nearest Neighbors. Cell Syst. 2019;8(4):329-37 e4. 1580 \n21. Macosko EZ, Basu A, Satija R, Nemesh J, Shekhar K, Goldman M, et al. Highly Parallel 1581 \nGenome-wide Expression Profiling of Individual Cells Using Nanoliter Droplets. Cell. 1582 \n2015;161(5):1202-14. 1583 \n22. Finak G, McDavid A, Yajima M, Deng J, Gersuk V, Shalek AK, et al. MAST: a flexible 1584 \nstatistical framework for assessing transcriptional changes and characterizing heterogeneity in 1585 \nsingle-cell RNA sequencing data. Genome Biol. 2015;16:278. 1586 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 72 \n23. Weber LM, Hippen AA, Hickey PF, Berrett KC, Gertz J, Doherty JA, et al. Genetic 1587 \ndemultiplexing of pooled single-cell RNA-sequencing samples in cancer facilitates effective 1588 \nexperimental design. Gigascience. 2021;10(9). 1589 \n24. Bose A, Beal MF. Mitochondrial dysfunction in Parkinson's disease. Journal of 1590 \nneurochemistry. 2016;139:216-31. 1591 \n25. Wickham H. ggplot2. Wiley interdisciplinary reviews: computational statistics. 1592 \n2011;3(2):180-5. 1593 \n26. Khan MRAA. Rocit-an r package for performance assessment of binary classifier with 1594 \nvisualization. 2019. 1595 \n27. Genomes Project C, Auton A, Brooks LD, Durbin RM, Garrison EP , Kang HM, et al. A 1596 \nglobal reference for human genetic variation. Nature. 2015;526(7571):68-74. 1597 \n28. Parkinson Progression Marker I. The Parkinson Progression Marker Initiative (PPMI). Prog 1598 \nNeurobiol. 2011;95(4):629-35. 1599 \n29. Hao Y , Stuart T , Kowalski MH, Choudhary S, Hoffman P , Hartman A, et al. Dictionary 1600 \nlearning for integrative, multimodal and scalable single-cell analysis. Nature biotechnology. 1601 \n2023:1-12. 1602 \n30. Thomas RA, Fiorini MR, Amiri S, Fon EA, Farhan SM. ScRNAbox: Empowering Single-Cell 1603 \nRNA Sequencing on High Performance Computing Systems. bioRxiv. 2023:2023.11. 13.566851. 1604 \n31. Li H. A statistical framework for SNP calling, mutation discovery, association mapping 1605 \nand population genetical parameter estimation from sequencing data. Bioinformatics. 1606 \n2011;27(21):2987-93. 1607 \n32. Kilpinen H, Goncalves A, Leha A, Afzal V, Alasoo K, Ashford S, et al. Common genetic 1608 \nvariation drives molecular heterogeneity in human iPSCs. Nature. 2017;546(7658):370-5. 1609 \n33. Stuart T , Butler A, Hoffman P , Hafemeister C, Papalexi E, Mauck WM, 3rd, et al. 1610 \nComprehensive Integration of Single-Cell Data. Cell. 2019;177(7):1888-902 e21. 1611 \n34. Streeter I, Harrison PW, Faulconbridge A, The HipSci C, Flicek P , Parkinson H, et al. The 1612 \nhuman-induced pluripotent stem cell initiative-data resources for cellular genetics. Nucleic 1613 \nAcids Res. 2017;45(D1):D691-D7. 1614 \n35. Yde Ohki CM, Grossmann L, Doring C, Hoffmann P , Herms S, Werling AM, et al. 1615 \nGeneration of integration-free induced pluripotent stem cells from healthy individuals. Stem 1616 \nCell Res. 2021;53:102269. 1617 \n36. Yde Ohki CM, Walter NM, Bender A, Rickli M, Ruhstaller S, Walitza S, et al. Growth rates 1618 \nof human induced pluripotent stem cells and neural stem cells from attention-deficit 1619 \nhyperactivity disorder patients: a preliminary study. J Neural Transm (Vienna). 2023;130(3):243-1620 \n52. 1621 \n37. Chen EY , Tan CM, Kou Y , Duan Q, Wang Z, Meirelles GV, et al. Enrichr: interactive and 1622 \ncollaborative HTML5 gene list enrichment analysis tool. BMC Bioinformatics. 2013;14:128. 1623 \n38. Zhang X, Lan Y , Xu J, Quan F, Zhao E, Deng C, et al. CellMarker: a manually curated 1624 \nresource of cell markers in human and mouse. Nucleic Acids Res. 2019;47(D1):D721-D8. 1625 \n39. Hao Y , Hao S, Andersen-Nissen E, Mauck WM, 3rd, Zheng S, Butler A, et al. Integrated 1626 \nanalysis of multimodal single-cell data. Cell. 2021;184(13):3573-87 e29. 1627 \n40. Macosko EZ, Basu A, Satija R, Nemesh J, Shekhar K, Goldman M, et al. Highly parallel 1628 \ngenome-wide expression profiling of individual cells using nanoliter droplets. Cell. 1629 \n2015;161(5):1202-14. 1630 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint \n\n 73 \n41. Ihaka R, Gentleman R. R: a language for data analysis and graphics. Journal of 1631 \ncomputational and graphical statistics. 1996;5(3):299-314. 1632 \n42. Wickham H, Wickham H. Data analysis: Springer; 2016. 1633 \n43. Azzalini A, Menardi G. Clustering via nonparametric density estimation: The R package 1634 \npdfCluster. arXiv preprint arXiv:13016559. 2013. 1635 \n 1636 \n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint","source_license":"CC-BY-4.0","license_restricted":false}