Ensemblex: an accuracy-weighted ensemble genetic demultiplexing framework for population-scale scRNAseq sample pooling

24 Multiplexing samples from distinct individuals prior to sequencing is a promising step toward 25 achieving population-scale single-cell RNA sequencing by reducing the restrictive costs of the 26 technology. Individual genetic demultiplexing tools resolve the donor-of-origin identity of pooled 27 cells using natural genetic variation but present diminished accuracy on highly multiplexed 28 experiments, impeding the analytic potential of the dataset. In response, we introduce Ensemblex: 29 an accuracy-weighted, ensemble genetic demultiplexing framework that integrates four distinct 30 algorithms to identify the most probable subject labels. Using computationally and experimentally 31 pooled samples, we demonstrate Ensemblex’s superior accuracy and illustrate the implications of 32 robust demultiplexing on biological analyses. 33 34

single-cell RNA sequencing, multiplexing, sample pooling, genetic demultiplexing, 35 induced pluripotent stem cells , differential gene expression, dopaminergic neurons, doublet 36 detection, accuracy-weighted probability, high-throughput sequencing 37 38 39 40 41 42 43 44 45 46 47 48 49 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 3

50 Single-cell RNA sequencing (scRNAseq) continues to revolutionize our molecular understanding 51 of biology by providing unprecedented insight into the transcriptional landscape of individual 52 cells. Unlike bulk RNAseq, where the RNA from all cells within a tissue is sequenced to produce 53 total expressional profiles across all cells, scRNAseq captures transcriptional signatures at a single-54 cell resolution, elucidating the diverse gene expression across distinct cell types and subtypes. 55 Differential gene expression (DGE) can then be calculated between subgroups of cells to reveal 56 cell type-specific expression changes between patient or treatment groups. However, scRNAseq 57 has come at the expense of increased costs, hindering its application for population-scale analyses, 58 which are critical for deriving clinico-pathological associations and characterizing the genetic 59 heterogeneity of complex diseases in biomedical sciences (1, 2). 60 61 In addition to the expense of separately capturing and sequencing cells from individual donors, the 62 costs of scRNAseq are exacerbated for cell cultures , such as those derived from induced 63 pluripotent stem cells (iPSC) (1). In particular, neurological diseases are difficult to study in human 64 tissue because access to post-mortem brains is limit ed and experimental manipulations are not 65 possible; in contrast, iPSC-derived cultures of neurons and other brain cells grown from 66 reprogrammed skin or blood cells of human donors are an excellent model of the brain (3). 67 However, iPSCs from each donor must be individually plated and differentiated in parallel , 68 presenting prohibitively high consumable and labour costs that render the methodology unfeasible 69 for population-scale analyses. Multiplexing cultures by pooling cells from multiple donors prior 70 to growth and differentiation , droplet capture , and sequencing, is one solution to address this 71

as it reduces costs by a factor of the number of samples multiplexed (4). Similarly, 72 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 4 samples such as tumor biopsies can be pooled at acquisition to realize the same benefits. In turn, 73 genetic demultiplexing tools are cost-effective, statistical frameworks that use the natural genetic 74 variation at sites of single-nucleotide polymorphisms (SNP) observed in the transcriptome to 75 cluster cells on the basis of their donor’s genotype. Importantly, genetic demultiplexing can be 76 informed by prior genotype information of the donors to improve demultiplexing accuracy and 77 facilitate the assignment of each cell back to its specific donor -of-origin, which is critical for 78 downstream analyses aiming to investigate discrepancies between subjects. At present, six genetic 79 demultiplexing tools have been developed for scRNAseq: Demuxalot (5) and Demuxlet (6) both 80 require prior genotype information as input ; Freemuxlet (6) relies entirely on the de novo 81 transcriptome and does not incorporate prior genotype information; and ScSplit (7), Souporcell 82 (8), and Vireo (9) provide versions of the algorithm that can work with and without prior genotype 83 information (Table 1). 84 85 A robust genetic demultiplexing tool is tasked with mitigating the addition of technical artifacts 86 into scRNAseq datasets by correctly classifying each pooled cell to its donor -of-origin, correctly 87 identifying heterogenic doublets (erroneous barcodes composed of two or more cells from distinct 88 subjects), and quantifying its confidence in the demultiplexed labels so that low -confidence 89 classifications can be eliminated from downstream analyses. While benchmarking analyses on the 90 available genetic demultiplexing tools have shown effectiveness for demultiplexing small sample 91 sizes, limitations emerge as the number of multiplexed samples approach a population scale (6) 92 (7) (8) (9). For example, using computationally multiplexed samples, Neavin et al. evaluated the 93 performance of genetic demultiplexing tools as the number of samples approached a population 94 scale and observed diminished demultiplexing accuracy with increasing numbers of pooled 95 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 5 samples, as well as notable classification discrepancies between tools (10). Furthermore, even at 96 small sample sizes, divergent assignments between genetic demultiplexing tools are common (8) 97 (9) (11). Another feature that has been shown to affect genetic demultiplexing performance is the 98 underrepresentation of samples in a pool , which is especially relevant for cell culture -based 99 multiplexed experiments, as variable growth rates in vitro across cell lines is common (12) (8) (9). 100 Genetic demultiplexing tools have also shown low concordance for identifying heterogenic 101 doublets, which should be removed prior to downstream analyses to avoid technical noise in the 102 data (10). Importantly, benchmarking analyses have repeatedly highlighted ScSplit’s poor 103 performance relative to the remaining tools (9) (10) (8) (11). The sum of these limitations calls to 104 question the robustness of the individual genetic demultiplexing tools for resolving the donor 105 identities of highly multiplexed samples, which represents an important hurdle for feasibly 106 achieving population-scale scRNAseq analysis. 107 108 In response to the divergent assignments commonly observed across tools, a consensus framework, 109 whereby only cells that show matching sample labels across all individual tools are retained for 110 downstream analyses, may appear sufficient to resolve the risk of introducing technical noise into 111 the data from misclassified cells. However, consensus frameworks are restricted to performing 112 only as well as the worst -performing tool , and genetic demultiplexing performance is highly 113 dataset dependent (10); thus, the overall performance of a consensus framework can vary 114 immensely between datasets. To this end, Neavin et al. proposed a majority vote framework for 115 genetic demultiplexing, whereby a cell is assigned to the sample called by the majority of tools 116 (10). However, this approach can be vulnerable to a subset of tools performing poorly on the 117 dataset, does not allocate additional weight to the votes of tools that perform more favourably on 118 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 6 the dataset, cannot account for instances when ties occur amongst tools , and cannot capture cells 119 that are correctly classified by only one tool. The sum of these limitations leads to the unnecessary 120 removal of cells from downstream analyses, reducing statistical power, especially for highly 121 multiplexed pools where each donor, on average, will have a lower representation of cells in the 122 pool. Moreover, the ability to capture the transcriptional profiles of rare cell types with scRNAseq 123 provides a notable advancement over bulk RNAseq and can strongly influence biological 124 interpretations (13); thus, investigators are reluctant to discard valuable cells in order to maximize 125 the analytic potential of their dataset. 126 127 To address the need for a robust genetic demultiplexing framework that can maximize the number 128 of confidently classified cells retained for downstream analyses , achieve high demultiplexing 129 accuracy for population-scale scRNAseq sample pooling, and maintain reliability across different 130 datasets, we developed Ensemblex: an accuracy -weighted ensemble genetic demultiplexing 131 framework designed to identify the most probable sample labels from each of its constituent tools 132 — Demuxalot, Demuxlet/Freemuxlet, Souporcell, and Vireo. Our ensemble method capitalizes on 133 combining distinct statistical frameworks for genetic demultiplexing while adapting to the overall 134 performance of its constituent tools on the respective dataset , making it resilient against a poorly 135 performing tool and facilitating a higher yield of cells for downstream analyses. The Ensemblex 136 workflow is assembled into a three-step pipeline — 1) accuracy-weighted probabilistic ensemble; 137 2) graph -based doublet detection; 3) Ensemble -independent doublet detection — and can 138 demultiplex pools with or without prior genotype information. 139 140 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 7 Here, we showcase Ensemblex’s improved demultiplexing performance across a variety of settings 141 through benchmarking analyses on a total of 141 computationally multiplexed pools with known 142 ground-truth sample labels ranging in size from 4 to 80 samples. We applied the ensemble method 143 to three diverse, experimentally multiplexed datasets: 1) non -small cell lung cancer (NS CLC) 144 dissociated tumor cells from 7 individuals with donor-specific oligonucleotide labels ; 2) iPSC -145 derived dopaminergic neurons (DaN) from 22 healthy individuals; and 3) iPSC -derived neural 146 stem cells (NSC) from 9 individuals with attention deficit hyperactivity disorder (ADHD) and 7 147 healthy controls. We demonstrate Ensemblex’s robustness across distinct datasets, its ability to 148 return a high proportion of confidently classified cells for downstream analysis, and the 149 implications that its improved demultiplexing performance has on biological interpretations of 150 multiplexed experiments. 151 Table 1. Summary of individual genetic demultiplexing tools. 152 Genetic demultiplexing tool Prior genotype information for genetic demultiplexing Included in the Ensemblex framework Demuxalot (5) Required Yes Demuxlet (6) Required Yes Freemuxlet (6) Not supported Yes ScSplit (7) Optional No Souporcell (8) Optional Yes Vireo (9) Optional Yes 153 154

and Discussion 155 Evaluating the performance of existing individual genetic demultiplexing tools 156 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 8 To evaluate the performance of individual genetic demultiplexing tools , we generated 157 computationally multiplexed pools using scRNAseq of 80 different iPSC lines from Parkinson’s 158 disease patients and healthy controls, which were differentiated towards a DaN state as part of the 159 Foundational Data Initiative for Parkinson’s Disease (FOUNDIN-PD) (14). Processed scRNAseq 160 data from the independent iPSC lines were merged to simulate sample-pooling using a previously 161 described protocol (9), which provided known ground -truth donor and doublet labels. We 162 generated 96 in silico pools ranging in size from 4 to 80 multiplexed samples, where each sample 163 corresponded to a unique donor-of-origin. The in silico pools averaged 17,396 cells per pool with 164 a constant 15% doublet rate. 165 166 Leveraging whole-genome sequencing (WGS) of the 80 donors from which the iPSC lines were 167 derived and the four genetic demultiplexing tools that can utilize prior genotype information — 168 Demuxalot, Demuxlet, Souporcell, and Vireo -GT — we first investigated the proportion of 169 correctly classified cells by the individual tools (Figure 1A). Across the 96 in silico pools, all tools 170 showed decreasing demultiplexing performance as the number of samples within the pool 171 increased. Souporcell demonstrated the largest decrease in the proportion of correctly classified 172 cells as the number of multiplexed samples increased from 4 (mean = 90.60%) to 80 (mean = 173 53.27%). In accordance with previous findings (10, 15) , the individual genetic demultiplexing 174 tools performed better on singlet classification than doublet detection, highlighting an avenue for 175 improved genetic demultiplexing accuracy by increasing the rate of heterogenic doublet 176 identification (Figure 1A). 177 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 9 178 Figure 1. Evaluation of existing individual genetic demultiplexing tools. Evaluation of genetic 179 demultiplexing tools with prior genotype information on 96 in silico pools with known ground -180 truth sample labels ranging in size from 4 to 80 multiplexed induced pluripotent stem cell (iPSC) 181 lines from genetically distinct individuals, averaging 17,396 cells per pool and a 15% doublet rate. 182 A) Line graphs showing the proportion of correctly classified singlets, doublets, and all cells by 183 each individual genetic demultiplexing tool across varying numbers of multiplexed iPSC lines in 184 a single pool (sample number). The large dots show the mean proportion of correct classifications 185 by an individual tool across replicates at a given sample size (n = 9 per pool size). The blue points 186 show the proportion of cells that were correctly classified by at least one individual genetic 187 demultiplexing tool: Demuxalot, Demuxlet, Souporcell, or Vireo -GT. B) Bar chart showing the 188 mean proportion of total cells from an individual pool correctly classified by only one genetic 189 demultiplexing tool. Error bars represent one standard deviation from the mean. (n = 9 per pool 190 size) C) Bar chart showing the proportion of correctly classified singlet cells labelled as 191 “unassigned” (ambiguous singlet assignments) due to assignment probabilities below the 192 recommended threshold of the respective genetic demultiplexing tool. Error bars represent one 193 standard deviation from the mean. (n = 9 per pool size). 194 195 We also investigated the proportion of cells that were correctly classified by at least one genetic 196 demultiplexing tool to designate the best possible performance of an ensemble method that 197 successfully incorporates every correct classification from its constituent tool s (Figure 1A ). 198 Across the 96 in silico pools, an average of 93.64% of cells were correctly classified by at least 199 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 10 one tool. In comparison, Demuxlet, which demonstrated the best overall performance amongst 200 individual tools , correctly classified 86.73% of cells, on average. Demuxalot was consistently 201 responsible for the highest proportion of cells correctly classified by only one tool ; 1.21% of 202 pooled cells , on average, were correctly classified by Demuxalot only, followed by Demuxlet 203 (mean = 0.83%), Vireo-GT (mean = 0.29%), and Souporcell (mean = 0.26%) (Figures 1B; 204 Additional File 1: Figure S1). Conversely, a consensus framework , correctly classified only 205 81.06% of cells, on average (data not shown). Based on these results, we reasoned that an ensemble 206 genetic demultiplexing method that can identify the most probable sample label from its 207 constituent tools, independent of a consensus assignment, would increase the yield of correctly 208 classified cells. 209 210 Next, we explored the frequency at which correctly classified singlets were labelled as unassigned 211 because their assignment probability failed to meet the tool’s recommended probability threshold. 212 Across the 96 in silico pools, Vireo-GT consistently showed the highest proportion of correctly 213 classified singlets with insufficient assignment probabilities (Vireo-GT mean = 7.86%) followed 214 by Demuxalot (mean = 5.91%), Demuxlet (mean = 2.44%) and Souporcell (mean = 2.34%) 215 (Figure 1C ). While a stringent probability threshold is important to prevent erroneous 216 classifications in downstream analyses, w e reasoned that the unnecessary removal of correctly 217 classified cells could be mitigated by a carefully calibrated ensemble method that allocates 218 additional assignment confidence to cells with matching sample labels across constituent tools , 219 despite low internal tool-specific assignment probabilities. 220 221 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 11 We repeated the above analyse s u sing the same 96 computationally multiplexed pools and the 222 genetic demultiplexing tools that do not require prior genotype information : Freemuxlet, 223 Souporcell, and Vireo . Here, we observed the same overarching limitations as when 224 demultiplexing with prior genotype information: 1) decreasing demultiplexing performance as the 225 number of multiplexed samples increased; 2) poor doublet detection performance compared to 226 singlet classification; 3) high rates of cells only correctly classified by a single tool ; and 4) 227 discarded correctly classified cells due to insufficient assignment probabilities (Additional File 1: 228 Figure S2). When we compared demultiplexing with and without prior genotype information, we 229 observed a trend towards a higher proportion of cells being correctly classified when prior 230 genotype information was available , as previously seen in separate benchmarking analyses (9) 231 (Additional File 1: Figure S3). 232 233 Validating the Ensemblex framework on pools with known ground-truth sample labels 234 To mitigate the limitations of the individual genetic demultiplexing tools and maximize the 235 analytic potential of multiplexed scRNAseq datasets, we developed Ensemblex (Figure 2A). The 236 Ensemblex workflow begins by demultiplexing pooled samples with four distinct demultiplexing 237 algorithms, followed by three steps: 1) accuracy-weighted probabilistic ensemble; 2) graph-based 238 doublet detection ; and 3) ensemble-independent doublet detection (Figure 2B) . As output, 239 Ensemblex returns its own cell-specific sample labels and corresponding assignment probabilities, 240 as well as the sample labels and corresponding assignment probabilities for each of its constituent 241 tools. 242 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 12 243 Figure 2. Characterization of the Ensemblex framework. Ensemblex is a probabilistic -244 weighted ensemble genetic demultiplexing framework for s ingle-cell RNA sequencing analysis, 245 which was designed to leverage the most probable sample labels from each of its constituent tools: 246 Demuxalot, Demuxlet, Souporcell, and Vireo when using prior genotype information or 247 Demuxalot, Freemuxlet, Souporcell, and Vireo when prior genotype information is not available. 248 A) The Ensemblex workflow begins with demultiplexing pooled cells from genetically distinct 249 individuals by each of the constituent tools. The outputs from each individual demultiplexing tool 250 are then used as input into the Ensemblex framework. B) The Ensemblex framework comprises 251 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 13 three distinct steps that are assembled into a pipeline: 1) accuracy-weighted probabilistic ensemble, 252 2) graph -based doublet detection, and 3) ensemble -independent doublet detection. C-D) Line 253 graphs showng t he contribution of each step of the Ensemblex framework on 96 in silico pools 254 with known ground -truth sample labels ranging in size from 4 to 80 multiplexed induced 255 pluripotent stem cell (iPSC) lines from genetically distinct individuals, averaging 17,396 cells per 256 pool and a 15% doublet rate. The average proportion of correctly classified C) singlets and D) 257 doublets across replicates at a given pool size is shown after sequentially applying each step of the 258 Ensemblex framework with prior genotype information ( n = 9 per pool size) . The right panels 259 show the average proportion of correct classifications across all 96 pools; error bars represent one 260 standard deviation from the mean. The blue points show the proportion of cells that were correctly 261 classified by at least one individual genetic demultiplexing tool: Demuxalot, Demuxlet, 262 Souporcell, or Vireo-GT. 263 264 In response to our observation that certain cells are correctly classified by only one tool, we 265 implemented the accuracy-weighted probabilistic ensemble component (Step 1) of the Ensemblex 266 framework. In brief, this unsupervised weighting model identifies the most probable sample label 267 for each cell by assigning weights to each tool’s assignment probabilities based on their estimated 268 balanced accuracy for the dataset (see “Methods”) (Figures 2B) (16). Ensemblex then retains the 269 sample label with the highest cumulative probability across its constituents. However, one 270 challenge for this framework is computing the balanced accuracy of the constituent tools for 271 experimentally multiplexed pools that lack ground-truth labels. Therefore, to estimate the balanced 272 accuracy of a particular constituent tool (e.g., Demuxalot) without ground-truth labels, Ensemblex 273 leverages the cells with a consensus assignment across the three remaining tools ( e.g., Demuxlet, 274 Souporcell, and Vireo-GT) as a proxy for ground -truth. To validate this approach, we utilized in 275 silico pools with known ground truth sample labels to compute the Adjusted Rand Index (ARI) 276 between Ensemblex’s sample labels when the balanced accuracy of the constituent tool s was 277 computed using consensus labels or ground -truth labels. Here, we consistently observed a mean 278 ARI > 0.99, independent of the number of multiplexed samples in a pool, suggesting high 279 assignment concordance between the two approaches (Additional File 1: Figure S4). Applying 280 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 14 the accuracy-weighted probabilistic ensemble component to the 96 in silico pools correctly 281 classified 94.98% of singlets, on average, across all pools, approaching the number of singlets that 282 were correctly classified by at least one constituent tool (mean = 96.48%) (Figure 2C). In contrast, 283 only 66.01% of doublets, on average, were correctly identified across all pools after Step 1 , 284 compared to 76.59% of doublets that were correctly identified by at least one constituent tool 285 (Figure 2D). 286 287 Given that previous analyses have demonstrate d strong doublet call discordance across genetic 288 demultiplexing tools (10), it was unsurprising that Step 1 of the Ensemblex framework performed 289 poorly on doublet identification. Therefore, instead of relying on the cell type classifications of the 290 constituent tools (i.e., singlet or doublet), we elected to leverage the doublet-related features (e.g., 291 doublet probability; see “Methods”) returned by the constituent tools to identify the cells with the 292 highest doublet likelihood, independent of the existing classifications. We implemented this 293 approach in the graph-based doublet detection component (Step 2) of the Ensemblex framework , 294 which was specifically designed to increase the rate of true doublet detection . Step 2 begins by 295 identifying the top n most confident doublets in the pool (see “Methods”). Then, based on the 296 Euclidean distance s in principal component analysis ( PCA) space, the cells that appear most 297 frequently amongst the nearest-neighbors of the high confident doublets and exceed the optimized 298 percentile threshold for the nearest-neighbor frequency are labelled as doublets by Ensemblex 299 (Figure 2B ; Additional File 1: Figure S5; see “Methods”). Upon applying the graph-based 300 doublet detection component to the 96 in silico pools following Step 1 , Ensemblex correctly 301 identified 76.00% of doublets, on average: a 9.99% increase in doublet detection from Step 1. In 302 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 15 turn, the average proportion of correctly classified singlets across all pools (94.43%) decreased by 303 only 0.55% (Figure 2D). 304 305 The ensemble-independent doublet detection component (Step 3) of the Ensemblex framework 306 was implemented to further improve doublet detection. Step 3 was motivated by our observation 307 that certain tools, namely Demuxalot and Vireo, showed high doublet detection specificity (mean 308 = 0.99) on in silico pools with known ground-truth sample labels, but that Steps 1 and 2 failed to 309 incorporate a subset of these correct doublet calls (Additional File 1: Figure S6). Therefore, by 310 default, Ensemblex accepts the doublet calls made by Demuxalot and Vireo -GT (Figure 2B ). 311 Applying the ensemble-independent doublet detection component to the 96 in silico pools 312 following Steps 1 and 2 further increased the average proportion of correctly identified doublets 313 across all pools by 1.58% for a total of 77.63% of doublets detected , while only decreasing the 314 average proportion of correctly classified singlets by 0.13% for a total of 94.30% of singlets 315 correctly classified (Figures 2C and 2D ). Notably, owing to the graph -based doublet detection 316 component, the average proportion of doublets identified by Ensemblex exceeded the average 317 proportion of doublets that were correctly classified by at least one constituent tool. 318 319 While the three-step workflow of the Ensemblex pipeline was designed to maximize the balance 320 between singlet classification and doublet identification, we do prioritize the identification of 321 doublets at the expense of a slightly lower singlet yield to minimize technical noise in the data. 322 However, we recognize that different experimental designs will require varying levels of doublet 323 detection stringency ; thus , users can modify the percentile thresholds for graph -based doublet 324 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 16 detection and nominate different tools for ensemble -independent doublet detection (see 325 “Methods”). 326 327 Benchmarking Ensemblex on pools with known ground-truth sample labels 328 To benchmark Ensemblex against Demuxalot, Demuxlet, Souporcell, and Vireo-GT with prior 329 genotype information, we first utilized the 96 in silico pools with known ground -truth sample 330 labels to assess how Ensemblex’s demultiplexing performance varied as the number of multiplexed 331 samples approached a cohort scale (4-80 samples). Unlike doublets, singlets were only considered 332 correctly classified if their assignment probability exceeded the recommended threshold of the 333 respective tool. On average across all pools, Ensemblex showed a higher proportion of correctly 334 classified singlets ( mean = 92.19%), doublets ( mean = 77.63%), and all cells (mean = 90.12%) 335 than the other tools. In comparison, Demuxlet, widely considered the “gold standard” tool, 336 correctly classified 89.72% of singlets, 68.57% of doublet s, and 86.73% of all cells, on average 337 (Figures 3A -3C). Importantly, the discrepancy in the proportion of correctly classified cells 338 between Ensemblex and the next-best tool was amplified as the number of multiplexed samples 339 increased from 4 (2.78%) to 80 ( 3.52%), demonstrating that our ensemble method was able to 340 partially mitigate decreased demultiplexing accuracy as the pools approach a population scale. 341 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 17 342 Figure 3. Ensemblex ground-truth benchmarking on computationally multiplexed pools. The 343 genetic demultiplexing tools with prior genotype information were evaluated on 96 in silico pools 344 with known ground -truth sample labels ranging in size from 4 to 80 multiplexed induced 345 pluripotent stem cell (iPSC) lines from genetically distinct individuals, averaging 17,396 cells per 346 pool and a 15% doublet rate. A singlet was considered correctly classified if the assigned sample 347 label matched the ground-truth sample label and the assignment probability exceeded the 348 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 18 recommended threshold for the respective tool; a doublet was considered correctly classified if the 349 assigned sample label matched the ground -truth sample label, regardless of the assignment 350 probability. A-I) Line graphs showing the performance of Ensemblex and the individual genetic 351 demultiplexing tools across evaluation metrics. The large dots show the mean value for each tool 352 across replicates at a given sample size (n = 9 per pool size). A) Proportion of correctly classified 353 singlets. B) Proportion of correctly classified doublets. C) Proportion of correctly classified cells. 354 D) Adjusted Rand Index between each tool’s sample labels and the ground-truth sample labels. E) 355 Balanced accuracy of each tool. F) Matthew’s Correlation Coefficient of each tool. G) Area under 356 the receiver operating characteristic curve (AUC) of the singlet assignment probability for each 357 tool. H) Proportion of usable cells returned by each tool. Usable cells were defined as cells 358 classified by singlets with an assignment probability exceeding the recommended threshold of the 359 respective tool. I) Error rate amongst the usable cells returned by each tool; erroneous 360 classifications comprised of true doublets labeled as singlets or true singlets assigned to the wrong 361 sample. 362 363 Next, we applied evaluation metrics for classification models to gauge the overall performance of 364 the genetic demultiplexing tools. We first computed the ARI to evaluate the similarity between the 365 demultiplexed sample labels and the ground -truth sample labels . Here, Ensemblex showed the 366 highest ARI with the ground truth sample labels across all pools (mean = 0.76), followed by 367 Demuxalot ( mean = 0.67) and Demuxlet ( mean = 0.66) (Figure 3D ). We then computed the 368 balanced accuracy to evaluate the binary classification performance — singlet or doublet — of 369 each genetic demultiplexing tool as well as the Matthew’s Correlation Coefficient (MCC), which 370 previous work has suggest ed is more reliable and informative for classification cases where 371 positive (singlet) and negative (doublet) cases have the same analyt ic importance (17). Across all 372 pools, Ensemblex showed the highest balanced accuracy (mean = 0.80) and MCC (mean = 0.64), 373 whereas Demuxalot and Demuxlet showed average balanced accuracies of 0.74 and 0.75, 374 respectively, and both tools show ed an average MCC of 0.54 (Figures 3E and 3F). To evaluate 375 how well Ensemblex’s confidence score (see “Methods”) and each constituent tool’s assignment 376 probability corresponded to the accuracy of their singlet classification, we plotted the area under 377 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 19 the receiver operating characteristic curve (AUC). Although Demuxalot (mean = 0.99) and Vireo-378 GT (mean = 0.99) showed the highest AUC across all pools on average, Ensemblex’s AUC was 379 comparable (mean = 0.98) (Figure 3G). 380 381 Finally, we investigated the proportion of usable cells returned by each demultiplexing tool and 382 the error rate amongst usable cells. We define usable cells as singlet classifications exceeding the 383 recommended probability threshold of the respective tool, while the error rate amongst usable cells 384 constituted incorrectly classified singlets to the wrong donor-of-origin or true doublets incorrectly 385 classified as singlets. We observed that, on average, Ensemblex returned the highest proportion of 386 usable cells across all pools (82.66%), followed by Demuxlet (81.66%), Souporcell ( 81.01%), 387 Demuxalot (79.99%), and Vireo-GT (77.53%) (Figure 3H). Importantly, Ensemblex showed the 388 lowest error rate amongst usable cells (4.34%), followed by Demuxalot ( 4.43%), Demuxlet 389 (5.77%), Vireo-GT (6.16%), and Souporcell (21.82%) (Figure 3I). 390 391 Using computationally multiplexed pools comprised of 24 iPSC lines, we further assessed how the 392 performance of Ensemblex varied as a function of the number cells in a pool when prior genotype 393 information was available. Here, we observed that our ensemble method consistently outperformed 394 the individual demultiplexing tools (Additional File 1: Figure S 7). When cells are pooled 395 experimentally, it is reasonable to expect some iPSC lines to be underrepresented in the pool. 396 Therefore, to assess Ensemblex’s demultiplexing performance in the pre sence of an 397 underrepresented iPSC line, we produced computationally multiplexed pools comprising of 24 398 samples, with one sample showing varying degrees of under representation. Again, we observed 399 that Ensemblex consistently outperformed the individual tools (Additional File 1: Figure S8). 400 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 20 Finally, we repeated the above analyses to assess whether the benefits of using Ensemblex to 401 demultiplex with prior genotype information extended to cases where prior genotype information 402 is not available . In doing so , we observed a trend towards better overall performance by 403 Ensemblex; however, the discrepancy between Ensemblex and the top -performing individual 404 tools, namely Freemuxlet and Souporcell, was less pronounced than when demultiplexing with 405 prior genotype information (Additional File 1: Figures S9-S11). 406 407 Taken together, these results indicate that the Ensemblex framework mitigates the limitations of 408 the individual tools, leading to greater overall demultiplexing performance across computationally 409 multiplexed pools with known ground-truth labels . Ultimately, Ensemblex’s improved 410 demultiplexing performance translates to a higher recovery of usable cells for downstream 411 analyses as well as a higher accuracy amongst usable cells, limiting the unnecessary removal of 412 cells from the dataset and mitigating the introduction of technical artifacts into biological analyses. 413 414 Evaluating Ensemblex on experimentally pooled samples with donor-specific oligonucleotide 415 labels 416 To determine whether Ensemblex’s improved performance across the in silico pools is reflected in 417 real-world multiplexed experiments, we applied Ensemblex to an experimentally multiplexed pool 418 composed of NSCLC dissociated tumor cells from 7 donors, hereafter referred to as the NSCLC 419 dataset (18). Importantly, these NSCLC cells were labelled with donor-specific Cell Multiplexing 420 Oligonucleotides (CMOs) , pro viding a proxy for ground -truth sample labels to evaluate the 421 performance of the genetic demultiplexing tools. For this experiment, we used HTOdemux (19) to 422 assign the cells back to their donor -of-origin based on the CMO expression profiles. HTOdemux 423 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 21 confidently assigned 19,695 cells, of which 15,534 (78.87%) were assigned to individual donors 424 and 4,161 (21.13%) were assigned as doublets; 769 cells (3.76%) were unassignable at a positive 425 quantile of 0.99 and were excluded from downstream analyses (Figures 4A). Application of the 426 Ensemblex framework with prior genotype information to the NSCLC dataset achieved a singlet 427 true positive (TP) rate of 96.92% and doublet TP rate of 66.21% ( Figure 4B). To evaluate the 428 benefits of utilizing the entire Ensemblex workflow (Step s 1-3), we investigated the contribution 429 of each step of the Ensemblex framework to the overall demultiplexing accuracy. Applying graph-430 based doublet detection (Step 2) and ensemble -independent doublet detection (Step 3) to the 431 accuracy weighted assignments obtained from Step 1 increased the proportion of correctly 432 identified doublets by 14%, while slightly decreasing the proportion of correctly classified singlets 433 by 0.05% (Additional File 1: Table S1 ). Although users can elect to utilize different step -434 combinations of the Ensemblex pipeline, these results reaffirm that leveraging the entire workflow 435 maximizes the overall demultiplexing accuracy by achieving a meticulous balance between singlet 436 classification and doublet identification. 437 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 22 438 Figure 4. Evaluating Ensemblex on experimentally multiplexed cells using donor -specific 439 oligonucleotide labels as a proxy for ground -truth. Non-small cell lung cancer (NSCLC) 440 dissociated tumor cells from 7 individuals were pooled and labelled with donor-specific 441 oligonucleotide-labels. Cells were demultiplexed according to their expression of donor -specific 442 oligonucleotide labels by HTOdemux ; HTOdemux’s sample labels were used as a proxy for 443 ground truth. True positives (TP) singlets were defined as cells classified as singlets by both 444 HTOdemux and Ensemblex with matching sample labels; false positives (FP) singlets were 445 defined as cells classified as singlets by both HTOdemux and Ensemblex but assigned to different 446 donors. TP doublets were defined as cells classified as doublets by both HTOdemux and 447 Ensemblex; FP doublets were defined as cells classified as singlets by HTOdemux and doublets 448 by Ensemblex; false negatives (FN) doublets were defined as cells classified as doublets by 449 HTOdemux and singlets by Ensemblex. A) T-distributed Stochastic Neighbor Embedding (t-SNE) 450 visualization of HTOdemux’s sample labels. B) T-SNE visualization of Ensemblex’s 451 demultiplexing performance using HTOdemux’s sample labels as ground truth for singlets (left) 452 and doublets (right) . C) Bar plots showing the s inglet TP and FP rate s for each genetic 453 demultiplexing tool using HTOdemux’s sample labels as ground truth . D) Bar plots showing the 454 doublet TP and FP rates for each genetic demultiplexing tool using HTOdemux’s sample labels as 455 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 23 ground truth. E) Scatter plot showing the proportion of usable cells (confidently classified singlets) 456 and the corresponding usable cell error rate for each genetic demultiplexing tool. F) Adjusted Rand 457 Index, balanced accuracy, Matthew’s Correlation Coefficient, and area under the receiver operating 458 characteristic curve (AUC) of the singlet assignment probability for each genetic demultiplexing 459 tool. 460 461 Upon comparing Ensemblex’s demultiplexing performance with prior genotype information on 462 the NSCLC dataset to the individual genetic demultiplexing tools, it emerged that our ensemble 463

obtained the highest singlet and doublet TP rates (Figures 4C and 4D). Ensemblex and 464 Demuxlet also showed the lowest singlet false positive (FP) rates (0.25% and 0.21%, respectively), 465 indicating that singlets were least frequently assigned to the wrong donor -of-origin by these two 466

compared to Demuxalot (1.87%), Vireo-GT (3.91%), and Souporcell ( 11.94%). 467 Souporcell and Vireo-GT returned the highest proportion of usable cells (confidently classified 468 singlets; 88.21% and 86.51%, respectively); albeit, at the expense of high usable cell error rates 469 (22.91% and 13.53%, respectively) (Figure 4E). In turn, Ensemblex, Demuxalot, and Demuxlet 470 showed lower error rates across the usable cells (8.75%, 8.91%, and 9.51%, respectively), amongst 471 which Ensemblex returned the highest proportion of usable cells (83.77%) compared to Demuxalot 472 (83.64%) and Demuxlet ( 83.43%). Here, the relatively high error rate amongst usable cells 473 returned by each demultiplexing tool is attributed to true doublets classified as singlets . Finally, 474 we computed the ARI, balanced accuracy, MCC, and AUC for singlet detection for each tool and 475 observed that Ensemblex again outperformed the remaining tools (Figure 4F). We repeated the 476 above analyses without prior genotype information and observed a similar trend towards better 477 overall performance by Ensemblex (Additional File 1: Table S2 and Figure S12). Together, these 478

corroborate that Ensemblex’s improved performance on the in silico pools extends to 479 experimentally multiplexed samples. 480 481 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 24 Application of Ensemblex to experimentally pooled, highly multiplexed subjects 482 To evaluate Ensemblex’s demultiplexing performance on experimentally pooled, highly 483 multiplexed scRNAseq datasets with prior genotype information, we used pools containing iPSC 484 lines from 22 donors that were differentiated towards DaN by Jerber et al., hereafter referred to as 485 the DaN dataset (12) (Figure 5A). To capture the transcriptional changes throughout neurogenesis, 486 Jerber et al. performed scRNAseq of the iPSC lines grown in pooled cultures at days 11, 30, and 487 52 of differentiation (Figure 5A ). Using three technical replicates from each timepoint, we 488 obtained 84,746 cells after performing quality control as previously described (12) (Additional 489 File 1: Table S3). Each technical replicate was demultiplexed independently by Ensemblex and 490 its constituent tools. 491 492 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 25 Figure 5. Application of Ensemblex to highly multiplexed, experimentally pooled cultures of 493 differentiated dopaminergic neurons. A) Time line of iPSC pooling, dopaminergic neuron 494 (DaN) differentiation, and sample collection from the DaN dataset by Jerber et al. (12). Three 495 technical replicates at each time point (days 11, 30 and, 52 of differentiation) from pools containing 496 22 individual iPSC lines were used in the analysis. Across all timepoints and technical replicates, 497 84,746 cells were obtained for analysis. B) Uniform manifold approximation and projection 498 (UMAP) plots showing confidently assigned singlets or predicted doublets (blue) and ambiguous 499 singlets (singlet assignments with insufficient assignment probabilities; red) returned by each 500 demultiplexing tool. C) Stacked bar chart showing the proportion of confidently assigned singlets 501 or predicted doublets (blue) and ambiguous singlets (red) across technical replicates at each time 502 point returned by each demultiplexing tool. D) Boxplot showing the proportion of confidently 503 classified singlets across technical replicates and time points by each demultiplexing tool. 504 Wilcoxon rank-sum tests were used to compare the proportion of confidently classified singlets by 505 Ensemblex to that of its constituents (n = 9 pools). E) Bar chart showing the proportion of 506 overlapping ambiguous singlet assignments amongst demultiplexing tools across technical 507 replicates and time points (n = 9 pools) . F) Boxplot showing the Adjusted Rand Index (ARI) 508 assessing cluster stability across a range of 11 clustering resolutions (n clustering iterations = 25) 509 after removing doublets identified by each demultiplexing tool. Wilcoxon rank -sum tests were 510 used to compare the clustering ARI after removing Ensemblex doublets to the clustering ARI after 511 removing doublets identified by each constituent tool. * Adjusted P-value < 0.05; ** adjusted P -512 value < 0.01; *** adjusted P-value < 0.001 513 514 To characterize the relationship between Ensemblex and its constituent demultiplexing tools, we 515 computed the ARI between Ensemblex’s sample labels and those of its constituent as well as the 516 percent contribution of each tool to Ensemblex’s final sample labels ( Table 2). Notably, we 517 observed that across day 30 technical replicates Demuxlet showed an ARI of 0.063 with 518 Ensemblex and only contributed 29.74% to Ensemblex’s final sample labels. In contrast, across 519 day 11 and 52 technical replicates Demuxlet showed an ARI of 0.928 and 0.884, respectively, and 520 contributed 95.91% and 90.55%, respectively, to Ensemblex’s final sample label s. Importantly, 521 Demuxlet’s variable contribution to Ensemblex’s sample labels across sequencing time points 522 demonstrates our ensemble method’s ability to adapt to the relative performance of its constituent 523 tools and override the classifications of a poorly performing tool on the respective dataset. 524 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 26 Table 2. Application of Ensemblex to pooled cultures of dopaminergic neurons from 22 525 healthy controls. 526 ARI between Ensemblex and constituent tool assignments Percent contribution to Ensemblex assignments n usable cells n doublets Day 11 Day 30 Day 52 Day 11 Day 30 Day 52 Demuxalot 0.987 0.955 0.982 97.29% 94.75% 97.57% 75,962 8,279 Demuxlet 0.928 0.062 0.884 95.91% 29.74% 90.55% 57,567 6,614 Souporcell 0.883 0.876 0.912 91.62% 91.82% 93.84% 76,811 7,740 Vireo-GT 0.961 0.879 0.958 95.95% 88.80% 95.16% 75,933 6,115 Ensemblex NA NA NA NA NA NA 76,222 8,307 DoubletFinder NA NA NA NA NA NA NA 4,597 Pooled cultures of induced pluripotent stem cell (iPSC) lines from 22 healthy donors were 527 differentiated towards a dopaminergic neuron (DaN) fate and sequenced on days 11, 30, and 52 of 528 differentiation by Jerber et al. (12). For the analysis we used three technical replicates for each 529 sequencing timepoint. Each pool was demultiplexed independently by Ensemblex and its 530 constituent tools with prior genotype information. The Adjusted Rand Index (ARI) between 531 Ensemblex’s assignments and those of the constituent tools was computed across technical 532 replicates corresponding to each differentiation timepoint. The percent contribution represents the 533 proportion of assignments from each constituent tool that matched Ensemblex’s assignments. 534 Usable cells were defined as singlet classifications whose assignment probability exceeded the 535 recommended threshold of the respective tool. Abbreviations: NA = Not applicable. 536 537 To elucidate the discrepancy in Demuxlet’s contribution to Ensemblex’s sample labels across 538 sequencing time points, we investigated the proportion of ambiguous singlet assignment s from 539 Ensemblex and its constituents. Ambiguous singlets are defined as singlet classifications whose 540 assignment probabilities failed to meet the recommended threshold of the respective tool, leaving 541 the identity of the pooled cell unresolved. Across 84,746 cells, Souporcell (195 singlets; 0.23% of 542 cells) and Ensemblex (217 singlets; 0.26% of cells) showed the lowest proportion of ambiguous 543 singlet assignments, followed by Demuxalot (505 singlets; 0.60% of cells) and Vireo-GT (2,698 544 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 27 singlets; 3.18% of cells ). Strikingly, Demuxlet showed 20,565 ambiguous singlet assignments 545 (24.27% of cells), with 92.04% derived from day 30 technical replicates, reflecting Demuxlet’s 546 remarkably low contribution to Ensemblex’s sample labels for cells sequenced at this timepoint 547 (Figures 5B and 5C). In accordance with previous analyses (9, 10), Demuxlet was consistently 548 amongst the top performing constituent tools throughout our benchmarking analyses. Yet, its poor 549 performance across day 30 technical replicates illustrates how the accuracy of individual tools can 550 vary greatly between datasets, highlighting the importance of utilizing multiple distinct algorithms 551 for genetic demultiplexing. We compared the mean proportion of confidently classified singlets 552 across technical replicates from each time point (n = 9) between Ensemblex (99.72%) and each 553 constituent demultiplexing tool using a Wilcoxon rank-sum test. After correction for multiple 554 hypothesis testing, w e observed that the mean proportion of confidently classified singlets by 555 Ensemblex was significantly higher than Demuxalot (mean = 99.36%, P -value = 3.55e-3), 556 Demuxlet (mean = 75.82%, P-value = 1.55e-5), and Vireo-GT (mean = 96.71%, P-value = 1.55e-557 5) (Figure 5D). Thus, despite Demuxlet’s unusually poor performance across day 30 technical 558 replicates, Ensemblex still confidently classified 27,520 singlets (99.61% of singlet assignments) 559 from these pools. Indeed, our ensemble method mitigates the consequences of a poorly performing 560 constituent tool by outweighing the erroneous classifications. In contrast , using a consensus 561 framework returned only 7,446 confidently classified singlets from day 30 technical replicates 562 (20,074 fewer cells than Ensemblex), limiting the availability of data for downstream analyses. 563 564 To further evaluate the ambiguity amongst singlet classification, we investigated the intersection 565 of ambiguous singlets across demultiplexing tools, reasoning that cells that are most challenging 566 to demultiplex would be labelled as ambiguous across all tools (Figure 5E). The singlets that were 567 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 28 assigned as ambiguous by Ensemblex showed the highest ambiguous singlet rate across the 568 remaining tools (mean across all constituent tools = 73.04%; mean across Demuxalot, Demuxlet, 569 and Vireo-GT = 92.32%). In contrast, while Souporcell showed the lowest ambiguous singlet rate 570 overall, only 15.90% of its unassigned singlets, on average, were ambiguous across the remaining 571 tools. These results indicate that the cells labelled as ambiguous by Ensemblex represent the cells 572 that are most challenging to classify across the distinct demultiplexing algorithms. Indeed, limiting 573 Ensemblex’s ambiguous singlet assignments to those that are most difficult to classify is critical 574 for maintaining a balance between maximizing the number of usable cells and minimizing the 575

of technical artifacts into downstream analyses from misclassified cells. 576 577 Next, we compared the doublet predictions made by each genetic demultiplexing tool and 578 DoubletFinder, a doublet detection tool that predicts doublets by estimating the similarity of the 579 transcriptional profile of a pooled cell to artificial doublets generated by combining the 580 transcriptional profiles of randomly selected cell pairs (20). Although the average number of 581 unique molecular identifiers (UMI) per cell across doublets identified by each tool was 582 significantly higher than the consensus singlets ( Additional File 1: Figure S13), we observed a 583 notable discrepancy in the number of doublets identified by each tool ; DoubletFinder identified 584 the fewest doublets (n = 4,597), while Ensemblex identified the most doublets (n = 8,307) (Table 585 1). Accordingly, all tools identified doublets that every other tool assigned as singlets (Additional 586 File 1: Figure S13). While Ensemblex identified the highest number of doublets, it still returned 587 a higher number of confidently classified singlets (n = 76,222) than Demuxalot ( n = 75,962), 588 Demuxlet (n = 57,567), and Vireo-GT (n = 75,933). Thus, even though the Ensemblex framework 589 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 29 prioritizes the identification of doublets at the expense of a slight ly lower singlet classification 590 rate, our ensemble method still returns a high proportion of usable cells for downstream analyses. 591 592 To evaluate the impact of doublet removal on the stability of clusters in the DaN dataset, we 593 performed 25 different random start iterations of the Louvain network detection at various 594 clustering resolutions after removing the doublets identified by each tool (21). Removing the 595 doublets identified by Ensemblex resulted in the highest ARI (mean ARI = 0.942), on average, 596 across clustering resolutions (Figure 5F ), suggesting the greatest cluster stability . However, 597 Wilcoxon rank -sum tests only revealed a statistically significant difference in the cluster 598 assignment ARI between Ensemblex and Souporcell (mean ARI = 0.922, P-value = 1.08e-2) after 599 correction for multiple hypothesis testing. Nonetheless, the highest cluster stability after removal 600 of Ensemblex’s putative doublets illustrates how improved doublet detection can translate to 601 improved biological analyses and is reflective of its superior doublet identification performance 602 on the benchmarking analyses. 603 604 Evaluating the impact of demultiplexing tools on differential gene expression analysis 605 To evaluate the impact of genetic demultiplexing tool s on scRNAseq DGE analysis, we 606 multiplexed iPSC-derived NSCs from individuals with ADHD and controls (Figure 6A). NSCs 607 were pooled and cultured until 100% confluence was reached. Two multiplexing experiments were 608 performed: Experiment 1 (n ADHD = 7; n control = 6) and Experiment 2 (n ADHD = 9; n control 609 = 7). After filtering cells for > 500 total and unique RNA transcripts , we obtained 30,433 cells 610 across both pools. Louvain clustering on the integrated scRNAseq dataset identified 12 clusters, 611 which were annotated as eight putative cell types (Figure 6B). 612 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 30 613 Figure 6. Evaluating the impact of discordant assignments between genetic demultiplexing 614 tools on d ifferential gene expression analysis. A) Schematic illustrating the workflow for the 615 neural stem cell (NSC) dataset. Pooled induced pluripotent stem cell (iPSC) -derived neural stem 616 cell cultures from individuals with attention deficit hyperactivity disorder (ADHD) and controls 617 were collected in two separate experiments. NSCs were dissociated for single -cell RNA 618 sequencing and prior genotype information of the pooled subjects was obtained through 619 microarray genotyping. The pools were demultiplexed by Ensemblex and its constituents with 620 prior genotype information and differential gene expression (DEG) was computed between ADHD 621 and controls. B) Uniform manifold approximation and projection (UMAP) plot showing the 622 putative cell types. C) Summary of the number of usable cells — singlets above the recommended 623 probability threshold of the respective demultiplexing tool — assigned to ADHD donors and 624 controls and the number of identified doublets by each demultiplexing tool. D) Boxplot showing 625 the Adjusted Rand Index (ARI) assessing cluster stability across a range of 11 clustering 626 resolutions (n clustering iterations = 25) after removing doublets identified by each demultiplexing 627 tool. A one-way Analysis of Variance (ANOV A) test comparing the ARI after removing doublets 628 identified by each tool revealed a significant difference between tools (n = 11 clustering 629 resolutions; P-value = 1.18e-3). E) Proportion of ADHD and control cells identified as putative 630 doublets by Ensemblex that were assigned as singlets by the constituent demultiplexing tools. F) 631 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 31 Heatmap showing the number of cell-type specific DEGs between ADHD and controls using the 632 subject labels of each demultiplexing tool. G) Heatmap showing the number of cell-type specific 633 DEGs between ADHD and controls using the subject labels of each demultiplexing tool and 634 removing putative doublets identified by Ensemblex. Cell-types not shown in the heatmaps had no 635 DEGs passing the adjusted P-value = 0.5| threshold across all tools. 636 637 We independently demultiplexed both pools using Ensemblex and its constituents to assign the 638 cells back to their donor -of-origin with prior genotype information (Figure 6C). The number of 639 cells assign ed to ADHD and control donors by each genetic demultiplexing tool is shown in 640 Additional File 1: Table S6. Importantly, the NSC dataset provides a valuable illustration of the 641 consequences of unnecessarily discarding cells from downstream analyses. For example, 642 Ensemblex and Vireo -GT returned 2,387 and 882 confidently assigned GRIA1high neurons, 643 respectively, whereas a consensus approach would have confidently assigned only 563 GRIA1high 644 neurons (Additional File 1: Table S6). 645 646 Each genetic demultiplexing tool predicted the ADHD cells to be vastly underrepresented 647 compared to the control cells; Ensemblex assigned 2,739 cells to individuals with ADHD and 648 19,880 cells to controls, suggesting that the ADHD iPSC lines were lost throughout the culturing 649 and sequencing process (Figure 6C). Additionally, we observed a notable difference in the number 650 of identified doublets across the tools; Vireo-GT identified the fewest doublets (n = 2,707), while 651 Demuxlet identified the most doublets ( n = 8,329) (Figure 6C). We aimed to characterize the 652 change in cluster stability after removing the doublets identified by each tool and observed that 653 removing the doublets identified by Ensemblex resulted in the highest ARI (mean ARI = 0.995), 654 on average, across clustering resolutions (Figure 6D ). A one-way ANOV A test comparing the 655 clustering ARI after removal of doublets identified by each tool revealed a significant difference 656 between tools (P-value = 1.18e-3). Demuxlet (n = 8,329) identified more doublets than Ensemblex 657 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 32 (n = 6,373), but exhibited lower cluster stability (ARI), suggesting that increased cluster stability 658 is not merely representative of the number of doublet s removed but rather the quality of doublet 659 removal. 660 661 Given the underrepresentation of ADHD cells across the dataset, we elected to investigate the cells 662 that were identified as doublets by Ensemblex but assigned as singlets by the constituent tools and 663 how these putative doublets were distributed across samples according to disorder status. 664 Demuxalot (n = 388) and Demuxlet (n = 726) assigned a relatively low number of Ensemblex’s 665 doublets as singlets , which represented 0.66% and 4.58% of ADHD sample assignments, 666 respectively, and 1.97% and 3.58% of control sample assignments, respectively (Figure 6E). In 667 contrast, Souporcell (n = 3,902) and Vireo-GT (n = 1,334) assigned a relatively high number of 668 Ensemblex’s doublets as singlets , which represented 31.97% and 24.88% of ADHD sample 669 assignments, respectively, and 11.65% and 3.97% of control sample assignments, respectively, 670 illustrating how variable doublet detection can impact the assembly of cells assigned to don or 671 categories and which cells are retained for downstream analyses. 672 673 Finally, w e used the model-based analysis of single -cell transcriptomics (MAST) statistical 674 framework to compute cell-type specific DGE between individuals with ADHD and controls using 675 the demultiplexed sample labels from each tool (22). We observed a significant discrepancy in the 676 number of cell type -specific differentially expressed gene s (DEGs; adjusted P-value 0.5) depending on the demultiplexing tool used (Figure 6F). Most 678 notably, for glia cells Souporcell identified 116 DEGs; Vireo-GT identified 98 DEGs; Ensemblex 679 identified 7 DEGs; Demuxalot identified 6 DEGs; and Demuxlet identified 1 DEG. Similar 680 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 33 patterns were observed across SOX2high NSCs, POU5F1high neural progenitor cells (NPC) , 681 S100Bhigh NPCs, and DCXhigh neurons, whereby Souporcell or Vireo-GT’s sample labels resulted 682 in a remarkably high number of DEGs compared to Ensemblex, Demuxalot, and Demuxlet. Given 683 that Souporcell and Vireo-GT made relatively few doublet call s and that 31.97% and 24.88% of 684 ADHD sample assignments made by Souporcell and Vireo -GT, respectively, were putative 685 doublets identified by Ensemblex, we elected to repeat the DGE analysis using the demultiplexed 686 sample labels from each tool but this time we removed all putative doublets identified by 687 Ensemblex. In doing so, we observed a decrease in the number of DEGs identified by Souporcell 688 and Vireo-GT across cell types , suggesting that the putative doublets identified by Ensemblex, 689 which were classified as singlets by Souporcell and Vireo -GT, were driving the initial signals 690 (Figure 6G ). For example, the number of glia-specific DEGs decreased from 116 to 0 with 691 Souporcell’s sample labels, and 98 to 0 with Vireo-GT’s sample labels. Given that the NSC dataset 692 lacked ground -truth sample labels, we could not definitively determine which cells were true 693 doublets; however, the increase in clustering ARI after removal of Ensemblex’s putative doublets 694 (Figure 6D), coupled with Ensemblex’s improved doublet identification performance on pools 695 with known ground -truth sample labels ( Figure 2B ), afforded confidence to assume that our 696 ensemble method performed favorably. Nonetheless, this analysis reveals that the choice of 697 demultiplexing tool can greatly impact biological analyses. 698 699

700 Multiplexing protocols, coupled with the introduction of genetic demultiplexing tools constituted 701 a significant advancement for scRNAseq by providing a feasible means to dramatically increase 702 the throughput of biological replicates . As the demand for population-scale scRNAseq analysis 703 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 34 continues to grow with the maturation of singe-cell technologies, the prospect of multiplexing 704 entire cohorts has emerged. However, the realization of this goal is impeded by the limitations of 705 the current genetic demultiplexing tools. These include decreasing demultiplexing performance as 706 the number of multiplexed samples increase s (9, 10) , relatively poor doublet detection 707 performance (10), relatively high rates of cells that can only be correctly classified by single 708 algorithms, the unnecessary removal of correctly classified cells due to insufficient assignment 709 probabilities, and highly variable demultiplexing performance between datasets (10). In this work 710 we presented Ensemblex, which offers a unique solution to these limitations by meticulously 711 implementing distinct demultiplexing algorithms into a robust, accuracy -weighted ensemble 712 framework that is exceptionally equipped to classify highly multiplexed pools. 713 714 We applied Ensemblex to a diverse array of computationally and experimentally multiplexed 715 scRNAseq datasets. Benchmarking analyses on pools with known ground -truth sample labels 716 revealed Ensemblex’s superior demultiplexing performance across pools reaching 80 multiplexed 717 samples, which translated to a higher proportion of cells retained for downstream analyses and 718 lower error rates amongst classified cells. Ensemblex also demonstrated a notable advancement 719 for identifying heterogenic doublets, which is a well -documented limitation of the genetic 720 demultiplexing tools currently available (9, 10, 15) . While previous analyses indicated that the 721 number of multiplexed samples in a pool directly impacted doublet detection efficiency (15), we 722 showed that Ensemblex’s ability to identify doublets remained relatively constant when >24 723 samples were multiplexed. Our findings suggest that super loading cells prior to sequencing —724 which will result in a higher number of usable cells but a higher a doublet rate (6) — followed by 725 heterogenic doublet detection by Ensemblex, may be a viable approach for implementing 726 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 35 population-scale multiplexing in practice. We also demonstrated that the performance of individual 727 genetic demultiplexing tools can be highly dataset -dependent, reflecting the findings of previous 728 work (10). However, due to its unsupervised weighting model, we showed that Ensemblex is 729 resistant to poorly performing constituent tools, maximizing the consistency of its demultiplexing 730 performance. Nonetheless, if each constituent tool performs poorly on a given dataset, the poor 731 performance will be reflected in Ensemblex’s demultiplexing accuracy. Finally, we illustrated that 732 discordant sample assignments amongst genetic demultiplexing tools can greatly impact DGE 733 analyses, necessitating that investigators carefully consider their choice of genetic demultiplexing 734 tool. Although untested, we anticipate that the impacts of discordant sample assignments amongst 735 genetic demultiplexing tools on biological interpretations would be exacerbated for computational 736 analyses that consider the specific donor identity of the pooled cells, such as expression 737 quantitative trait loci (eQTL) analyses, as opposed to donor groups (i.e., case and control). Due to 738 Ensemblex’s ability to seamlessly integrate multiple algorithms into an adaptable framework, we 739 argue that our ensemble method achieves unmatched reliability for experimentally multiplexed 740 pools that lack ground truth sample labels. 741 742 Undoubtedly, a limitation of utilizing an ensemble method for genetic demultiplexing is the 743 necessity to run each individual demultiplexing algorithm, which can be computationally 744 expensive. Yet, in the absence of comparing demultiplexed sample labels across tools, poor 745 performance by a given individual algorithm on experimentally multiplexed pools is undetectable, 746 and the risk of introducing technical artifacts and losing usable cells for downstream analyses is 747 prominent. As such, we believe that the relatively high computational cost of Ensemblex is a 748 worthwhile investment to maximize the biological insight obtained from multiplexed scRNAseq 749 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 36 datasets. To mitigate the burden of genetic demultiplexing by multiple individual tools, we provide 750 a coherent pipeline that runs each constituent demultiplexing tool in parallel and seamlessly 751 processes the respective output files with the Ensemblex algorithm. 752 753 Compared to when demultiplexing was informed by prior genetic data of the pooled samples, the 754 improvement of Ensemblex over its constituent tools was far less pronounced for genotype-free 755 demultiplexing cases. A ll demultiplexing tools, including Ensemblex, showed drops in 756 demultiplexing performance when >16 samples were multiplexed in a pool without prior genotype 757 information. Nonetheless, Ensemblex still constitutes an advancement over the individual tools for 758 genotype-free demultiplexing cases due to the robustness achieved by incorporating distinct 759 demultiplexing algorithms, which protects against the prospect of poorly performing individual 760 tools on the respective dataset. Furthermore, an intrinsic limitation of demultiplexing without prior 761 genotype information is that samples cannot be directly linked to metadata, leaving the sample 762 identity of the inferred clusters unresolved (9). Although challenging, this limitation can be 763 mitigated by identifying a small subset of discriminatory variants from the reconstructed genotypes 764 of the constituent demultiplexing tools, which could be used to manually assign the computed 765 clusters to samples if such discriminatory variants are known by the investigator. While the 766 Ensemblex pipeline provides users the option to demultiplex pools with or without prior genotype 767 information, we assert that users take caution when electing to perform population -scale 768 multiplexing experiments without using prior genetic data. 769 770 Genetic demultiplexing tools have been used extensively for scRNAseq analysis across many 771 disciplines in the biological sciences, including microbiology (8), model organisms (15), cancer 772 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 37 biology (23), and neurodegenerative disease (12). Recent work has also evaluated the utility of 773 genetic demultiplexing tools for different single-cell, read-based modalities such as single-nuclei 774 RNA sequencing (snRNAseq) and single -nuclei assay for transposase -accessible chromatin 775 sequencing (scATACseq) (24). Although untested, we expect Ensemblex to prove beneficial in 776 demultiplexing for these assays, but comprehensive benchmarking with the appropriate datasets is 777 required and was not explored here. 778 779 We expect numerous biological fields to exploit the benefits of Ensemblex through its application 780 to highly multiplexed pools comprising cells from many genetically distinct individuals. 781 Specifically for biomedical sciences, the preparation and labour costs of scRNAseq remains 782 prohibitively expensive for analyzing entire cohorts of patients, which is critical for characterizing 783 the genetic heterogeneity and etiological diversity of disease, and for maint aining sufficient 784 statistical power for detecting associations between transcriptional changes and clinical or 785 pathological observations (1). By increasing the throughput of biological replicates, multiplexing 786 has rendered the prospect of analyzing entire patient cohort s with single -cell transcriptomics 787 feasible. Highly-multiplexed scRNAseq experiments have already been presented in the literature 788 and, to the best of our knowledge, have pooled up to 24 samples in a single dish (12). However, 789 we demonstrated that Ensemblex’s demultiplexing accuracy remains relatively constant when >24 790 samples are multiplexed at concentrations that abide by the current limitations of experimental 791 protocols, suggesting that Ensemblex equips the research community with the necessary 792 computational framework to expand the upper limits of the number of genetically distinct 793 individuals in a single pool. 794 795 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 38 While multiplexing mitigates the labour and consumable costs of scRNAseq analysis, the cost of 796 sequencing remains expensive and the increasing number of genetically distinct individuals in a 797 single pool necessitates that a greater number of cells must be sequenced to ensure adequate 798 representation. Accordingly, Ensemblex is equipped to demultiplex pools comprising cells from 799 more genetically distinct individuals than is feasible with the current laboratory technologies. 800 However, we expect that the cost of sequencing will continue to decrease with the maturation of 801 the technology, and our tool will be in place for when the anticipated wet lab advancements are 802 realized. Overall, we conclude that Ensemblex constitutes a notable advancement towards the 803 pressing demand for population-scale single-cell transcriptomics. 804 805

806 Ensemblex framework overview 807 Ensemblex is an ensemble genetic demultiplexing framework for scRNAseq sample pooling that 808 was designed to identify the most probable sample labels from each of its constituent tools : 809 Demuxalot (5), Demuxlet (6), Souporcell (8), and Vireo (9) when demultiplexing with prior 810 genotype information or Demuxalot, Freemuxlet (6), Souporcell, and Vireo when demultiplexing 811 without prior genotype information. After running each constituent demultiplexing tool in parallel, 812 Ensemblex merges the output files containing the sample-cell assignments from each tool and 813 performs three distinct steps of the Ensemblex pipeline: 814 1. Accuracy-weighted probabilistic ensemble; 815 2. Graph-based doublet detection; 816 3. Ensemble-independent doublet detection. 817 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 39 Upon obtaining the final Ensemblex sample labels (donor-of-origin identity of the pooled cells) , 818 the singlet assignment confidence score is computed. 819 820 Step 1: Accuracy-weighted probabilistic ensemble 821 Ensemblex utilizes an unsupervised weighting model to identify the most probable sample 822 label for each cell. Ensemblex weighs each constituent tool’s assignment probability 823 distribution by its estimated balanced accuracy for the dataset in a framework adapted from 824 the work of Large et al. (16). To estimate the balanced accuracy of a particular constituent tool 825 (e.g., Demuxalot) for experimentally multiplexed datasets lacking ground -truth labels, 826 Ensemblex uses the cells with a consensus assignment across the three remaining tools (e.g. , 827 Demuxlet, Souporcell, and Vireo-GT) as a proxy for ground-truth. The balanced accuracy for 828 each tool is calculated using equation 1: 829 830 (1) 𝐵𝑎𝑙𝑎𝑛𝑐𝑒𝑑 𝑎𝑐𝑐𝑢𝑟𝑎𝑐𝑦= 1 2 (( 𝑇𝑃 𝑇𝑃+𝐹𝑁) + ( 𝑇𝑁 𝑇𝑁+𝐹𝑃)) 831 832 Where TP is the number of correctly classified singlets; true -negative (TN) is the number of 833 correctly classified doublets; FP is the number of incorrectly classified singlets; false- negative 834 (FN) is the number of incorrectly classified doublets . The probability distribution of each 835 constituent tool ( 𝑝𝑗̂) is then weighted by its estimated balanced accuracy ( 𝑤𝑗) to produce an 836 accuracy-weighted ensemble probability for each cell: 837 838 (2) 𝑝̂(𝑦 = 𝑖|𝐸) ∝ ∑ 𝑤𝑗𝑝̂𝑗(𝑦 = 𝑖|𝑀𝑗)𝑘 𝑗=1 839 840 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 40 Where 𝑝̂ is the probability that a barcode belongs to class 𝑖; 𝑦 is the class variable with 𝑐 841 possible values, 𝑦 ∈ (1, … , 𝑐); 𝑐 is the number of pooled samples plus 1 to account for 842 doublets; 𝐸 is a vector of the results of 𝑀 classifiers, 𝐸 = (𝑀1, … , 𝑀𝑘); 𝑀is the individual 843 constituent demultiplexing output from each tool. Given 𝑝̂, Ensemblex assigns each barcode’s 844 sample identity (𝑦̂) as the class (sample label) with the maximum probability: 845 846 (3) 𝑦̂= arg 𝑚𝑎𝑥𝑖∈(1,…,𝑐) 𝑝̂(𝑦 = 𝑖|𝐸) 847 848 Step 2: Graph-based doublet detection 849 Ensemblex employs a graph-based approach to identify doublets that are incorrectly labeled as 850 singlets by the accuracy-weighted probabilistic ensemble component (Step 1). For graph-based 851 doublet detection, Ensembl ex leverages pre-defined features returned from each constituent 852 tool: 853 1. Demuxalot: doublet probability; 854 2. Demuxlet/Freemuxlet: singlet log likelihood – doublet log likelihood; 855 3. Demuxlet/Freemuxlet: number of single nucleotide polymorphisms (SNP) per cell; 856 4. Demuxlet/Freemuxlet: number of reads per cell; 857 5. Souporcell: doublet log probability; 858 6. Vireo: doublet probability; 859 7. Vireo: doublet log likelihood ratio. 860 For each feature independently, the pooled cells are ordered from the most to the least probable 861 doublet and are assigned a percentile rank. Beginning with a percentile threshold of 99.99 , 862 Ensemblex screens each cell to identify those that exceed the percentile threshold across all 863 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 41 features; cells that exceed the percentile threshold across all features are labeled as “confident 864 doublets”. For each iteration, Ensemblex decreases the percentile threshold by 0.01 and repeats 865 the screening process until it has identified n confident doublets (nCD). Ensemblex performs 866 a parameter sweep to determine the optimal nCD to use for graph-based doublet detection (see 867 below). 868 869 Next, the above features are input into a PCA using the stats (v3.6.2) R package (25) and a 870 Euclidean distance matrix is generated from the first two principal components (PC). For each 871 confident doublet independently, the remaining cells in the pool are assigned a percentile rank 872 based on their proximity in Euclidean space to the confident doublet and the cells that exceed 873 the designated nearest neighbour percentile threshold ( pT) are identified. For all cells that 874 exceeded the designated pT for any confident doublet (putative doublets) , Ensemblex 875 computes the number of times the putative doublet was amongst the nearest neighbours of any 876 confident doublet (fNN); an fNN equal to nCD indicates that a putative doublet was amongst 877 the top nearest neighbours for each confident doublet. 878 879 To optimize the nCD and pT parameters for experimentally pooled samples lacking ground-880 truth labels, Ensemblex performs an automated parameter sweep at each pairwise combination 881 of nCD and pT values; nCD values range from 50 to 300, in increments of 50, while pT values 882 depend on the expected doublet rate (exDR) and range from 1 − 𝑒𝑥𝐷𝑅 6 to 1 − 𝑒𝑥𝐷𝑅, in 883 intervals of 1−𝑒𝑥𝐷𝑅 6 . The distribution of fNN values for each combination of nCD and pT 884 parameters are plotted and Pearson’s measure of kurtosis (K), is used to predict which 885 combination of pT and nCD values optimize the identification of true doublets while 886 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 42 minimizing the rate of incorrectly labelled true singlets as doublets. Ensemblex screens for 887 combinations of nCD and pT values that result in negatively skewed fNN distributions with 888 high K, signifying high peakedness and heavy tails. High peakedness indicates that cells 889 exceeding the designated pT concentrated around nCD, reflecting their proximity in Euclidean 890 space to all high confident doublets, while heavy tails indicate that even cells with lower fNN 891 values were identified as nearest neighbour to many confident doublets. Ensemblex first 892 identifies the pT that returns the highest K, on average, across nCD values tested in the 893 parameter sweep using equation 4: 894 895 (4) 𝑝𝑇̂= arg 𝑚𝑎𝑥𝑝𝑇∈{1− 𝑒𝑥𝐷𝑅 6 ,….,1−𝑒𝑥𝐷𝑅) ( ∑ 𝐾(𝑦=𝑝𝑇)𝑛𝐶𝐷∈{50,100,150,200,250,300} 2 ) 896 897 Where K of the distribution of fNN values of the putative doublets is defined as: 898 899 (5) 𝐾(fNN) = 𝐸[( 𝑋−𝜇 𝜎 ) 4 ] 900 901 Where 𝜇 is the mean of the distribution and 𝜎 is the standard deviation. Upon identifying the 902 optimal pT value ( 𝑝𝑇̂), Ensemblex plots the K corresponding to 𝑝𝑇̂ across all nCD values 903 tested in the parameter sweep. If an inflection point is identifiable, Ensemblex identifies 𝑛𝐶𝐷̂ 904 as the nCD value corresponding to the point of inflection on the curve. Otherwise, Ensemblex 905 identifies 𝑛𝐶𝐷̂ as the nCD value corresponding to the highest K. Cells flagged as putative 906 doublets identified using 𝑝𝑇̂ and 𝑛𝐶𝐷̂ are labelled as doublets by Ensemblex. 907 908 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 43 Step 3: Ensemble-independent doublet detection 909 Benchmarking on computationally multiplexed pools with known ground-truth sample labels 910 revealed that certain genetic demultiplexing tools, namely Demuxalot and Vireo, showed high 911 doublet detection specificity, but that Steps 1 and 2 of the Ensemblex workflow failed to 912 correctly label a subset of doublet calls by these tools. To mitigate this issue and maximize the 913 rate of doublet identification, Ensemblex labels the cells that are identified as doublets by Vireo 914 or Demuxalot as double ts by default; however, users can nominate different tools for the 915 ensemble-independent doublet detection component depending on the desired doublet 916 detection stringency. Doublet specificity was computed using equation 6: 917 918 (6) 𝐷𝑜𝑢𝑏𝑙𝑒𝑡 𝑠𝑝𝑒𝑐𝑖𝑓𝑖𝑐𝑖𝑡𝑦= ( 𝑇𝑁 𝑇𝑁+𝐹𝑃) 919 920 Where TN is the number of correctly classified doublets; FP is the number of true singlets 921 incorrectly classified as doublets. 922 923 Ensemblex singlet assignment confidence score 924 Ensemblex computes a singlet confidence score to inform which cells should be discarded to 925 avoid misclassification in downstream analyses. First, Ensemblex evaluates how well an 926 individual constituent tool’s assignment probability (e.g., Demuxalot) corresponded to the 927 accuracy of their assignment, using consensus cells across the three remaining tools ( e.g., 928 Demuxlet, Souporcell, Vireo) as a proxy for ground-truth, by fitting a binary logistic regression 929 model to compute the odds that a singlet was correctly classified given its corresponding 930 probability. Using the binary logistic regression models, Ensemblex computes the AUC using 931 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 44 the empirical method implemented in the ROCit (v2.1.1) R package for each tool (26). Then, 932 for each cell, if Ensemblex’s sample label matche s that of a constituent tool, and if the 933 assignment probability of the constituent tool supersedes its probability threshold, the tool’s 934 computed AUC is added to the accuracy-weighted probabilistic ensemble probability produced 935 in Step 1 to yield the confidence score. By default, singlet assignments with a confidence score 936 less than 1.00 are labelled as unassigned by Ensemblex. Ensemblex’s confidence score and the 937 designated threshold is a successful predictor of accurately classified singlets because singlets 938 will only achieve a confidence score ≥ 1 if: 939 1. All constituent tools show the same sample label (accuracy -weighted probabilistic 940 ensemble probability = 1.00); 941 2. At least one constituent tool confidently assigns the cell to an individual donor and the 942 constituent tool’s probability assignment adequately corresponds to the overall 943 accuracy of their singlet assignment. 944 945 Application of Ensemblex with and without prior genotype information 946 Given the dependencies of certain tools on prior genotype information, there are notable 947 differences between the Ensemblex workflows for demultiplexing with and without prior 948 genotype information . When demultiplexing with prior genotype information, Ensemblex 949 leverages the sample labels from Demuxalot, Demuxlet, and Vireo -GT with prior genotype 950 information, and Souporcell without prior genotype information. When demultiplexing 951 without prior genotype information, Ensemblex leverages the sample labels from Demuxalot, 952 Freemuxlet, Souporcell, and Vireo. However, given that Demuxalot requires prior genotype 953 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 45 information, Ensemblex uses the estimated donor .vcf file generated by Freemuxlet for input 954 into the Demuxalot algorithm as prior genetic data. 955 956 Running the Ensemblex pipeline 957 A complete user guide for running the Ensemblex pipeline can be found at the Ensemblex 958 GitHub site: https://neurobioinfo.github.io/ensemblex/site/. We provide two distinct yet highly 959 comparable pipelines depending on the availability of prior genotype information . Both 960 pipelines can be downloaded as a singularity image and are comprised of four steps: 961 1. Establish the pipeline and working directory; 962 2. Prepare input files for constituent genetic demultiplexing tools; 963 3. Parallel demultiplexing by constituent genetic demultiplexing tools; 964 4. Application of the Ensemblex algorithm for ensemble classification. 965 966 As input into the Ensemblex pipeline, users must provide a .tsv file describing the barcodes of 967 the pooled cells, a. bam sequencing file for the pool, a reference genotype .vcf file (e.g., 1000 968 Genome Project) (27), a reference genome sequence .fasta file (e.g., 10X Genomics), and, if 969 demultiplexing with prior genotype information, a .vcf file describing the genetic data of the 970 pooled samples. 971 972 Genetic demultiplexing by constituent tools 973 Genetic demultiplexing by the constituent demultiplexing tools was performed following best 974 practices as defined by the authors of the respective tools using Python (v3.8.10). 975 Demuxalot 976 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 46 CellRanger-generated .bam file, filtered barcode .tsv file, and the corresponding donor .vcf file 977 were used as input into the Demuxalot workflow. Candidate variants for scRNAseq genotyping 978 were retained if the minimum coverage was > 200 and minimum alternative coverage was > 979 10. The top 100 SNPs per donor were retained to cluster the cells by genotype. Doublet calls 980 were made with a prior strength of 0.25. 981 982 Demuxlet 983 We used the popscle suite ( https://github.com/statgen/popscle) for Demuxlet. CellRanger-984 generated .bam file, filtered barcode .tsv file, and the corresponding donor .vcf file were used 985 as input into the Demuxlet workflow. The dsc-pileup function was first used to pileup candidate 986 variants around known variant sites with the following parameters: --cp-BQ 40 --min-BQ 13 -987 -min-MQ 20 --minTD 0 --min-total 0 --min-uniq 0 --min-snp 0. The Demuxlet algorithm was 988 then applied to cluster the cells by genotype with the following parameters: --geno-error-offset 989 0.10 --geno-error-coeff 0.00 --min-callrate 0.50 --doublet-prior 0.50 --cap-BQ 40 --min-BQ 13 990 --min-MQ 20 --min-TD 0 --min-total 0 --min-uniq 0 --min-snp 0. 991 992 Freemuxlet 993 We used the popscle suite ( https://github.com/statgen/popscle) for Freemuxlet. CellRanger-994 generated .bam file, filtered barcode .tsv file, and reference genotype .vcf file from the 1000 995 Genomes Project, phase 3 (27), were used as input into the Freemuxlet workflow. The dsc-996 pileup function was first used to pileup candidate variants around known variant sites with the 997 following parameters: --cp-BQ 40 --min-BQ 13 --min-MQ 20 --minTD 0 --min-total 0 --min-998 uniq 0 --min-snp 0. The Freemuxlet algorithm was then applied to cluster the cells by genotype 999 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 47 with the following parameters: --doublet-prior 0.50 --bf-thres 5.41 --frac-init-clust 0.50 --inter-1000 init 10 --cap-BQ 40 --min-BQ 13 --min-total 0 --min-uniq 0 --min-snp 0. 1001 1002 Souporcell 1003 CellRanger-generated .bam file, filtered barcode .tsv file, 10X Genomics reference .fasta file, 1004 and the corresponding donor .vcf file when demultiplexing with prior genotype information 1005 were used as input into the Souporcell workflow. A FASTQ file was first generated from the 1006 .bam file using the renamer.py script. These reads were mapped to the reference genome using 1007 minimap2 with the following parameters: --ax splice –t 8 –G50k –k 21 –w 11 –sr --A2 –B8 –1008 O12,32 –E2,1 –r200 –p.5 –N20 –f1000,5000 –n2 –m20 –s40 –g200 –2k50m –secondary=no. 1009 The barcodes and UMI were added back to the .sam file using the retag.py script and the 1010 resulting .bam file was sorted and indexed with Samtools. Variants were called using Freebayes 1011 with the following parameters: --iXu –C 2 –q 20 –n 3 –E 1 –m 30 –min-coverage 6. Vartix was 1012 used to compute the number of alleles for each cell using the following parameters: --umi –1013 mapq 30 –scoring-method coverage. The Souporcell algorithm was then applied to cluster the 1014 cells by genotype; when demultiplexing with prior genotype information the --1015 known_genotypes and --known_genotypes_sample_names parameters were included. 1016 Troublet was used to identify doublets and the consensus.py script was used for genotype and 1017 ambient RNA co-inference. 1018 1019 Vireo 1020 CellRanger-generated .bam file, filtered barcode .tsv file, reference genotypes from the 1000 1021 Genomes Project, phase 3 (27), and the corresponding donor .vcf file when demultiplexing 1022 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 48 with prior genotype information were used as input to the Vireo workflow. CellSNP was used 1023 to identify candidate variants for scRNAseq genotyping with the following parameters: --1024 minMAF 0.1 and --minCOUNT 100. The Vireo algorithm was then applied to cluster the cells 1025 by genotype with the --forceLearnGT parameter; when demultiplexing with prior genotype 1026 information (Vireo-GT) the --d and --t GT parameters were used. 1027 1028 Consensus demultiplexing framework 1029 For the consensus demultiplexing framework, singlets were considered confidently classified 1030 if Demuxalot, Demuxlet, Vireo, and Souporcell assigned a cell to the same donor -of-origin. 1031 Cells classified as “ambiguous” or doublet by at least one tool were discarded. 1032 1033 Generation of computationally pooled samples for ground-truth benchmarking 1034 To benchmark Ensemblex on computationally pooled samples with known ground -truth sample 1035 labels, we leveraged 80 independently sequenced iPSC lines from Parkinson’s disease patients and 1036 healthy controls, which were differentiated towards a dopaminergic neuronal state and sequenced 1037 after 65 days of differentiation as part of the FOUNDIN-PD (14). Controlled access FASTQ files 1038 from the independently sequenced iPSC lines were obtained from https://www.ppmi-info.org/ 1039 (accessed 09-17-2023) and processed by the CellRanger counts pipeline (v3.1.0) with default 1040 parameters and aligned to GRCh38 reference genome. The CellRanger-generated .bam and filtered 1041 barcode files were used as input into the synth_pool.py script produced by the authors of Vireo to 1042 simulate sample pooling (9). In brief, reads from a subset of cells from the iPSC line-specific .bam 1043 files were merged and doublets were generated by combining the reads from random cell pairs. 1044 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 49 Sample identities were added to each cell’s barcode, revealing the ground-truth sample labels for 1045 benchmarking procedures. 1046 1047 To evaluate how genetic demultiplexing performance varied as a function of the number of 1048 multiplexed samples, we generated 96 computationally multiplexed pools using the 80 1049 FOUNDIN-PD lines with sample sizes of 4, 8, 16, 24, 32, 40, 48, 56, 64, 72, and 80. An equal 1050 number of cells from each line were used in the in silico pool. For the sample size of four we 1051 generated six replicates; for the sample sizes of 8-80 we generated nine replicates each. Replicates 1052 were produced with different sample and cell combinations. The 96 in silico pools averaged 17,396 1053 cells (minimum = 8,696; maximum = 26,087). For this experiment, we maintained a 15% doublet 1054 rate as previously described (9). 1055 1056 To evaluate how genetic demultiplexing performance varied as a function of the number of cells 1057 in a pool, we generated 18 computationally multiplexed pools using the 80 FOUNDIN-PD lines 1058 with 8,000, 16,000, 24,000, 32,000, 40,000, and 48,0000 pooled cells ; we generated three 1059 replicates per pool size. Twenty-four samples were multiplexed for each pool and an equal number 1060 of cells from each sample were used. Replicates were produced with different sample and cell 1061 combinations. For this experiment, we simulated a doublet rate of 6% per 8,000 pooled cells. 1062 1063 To evaluate if the overall demultiplexing performance varied due to the underrepresentation of a 1064 cell line, we generated 15 computationally multiplexed pools using the 80 FOUNDIN-PD lines 1065 comprising 23 multiplexed samples with 1,000 cells and one randomly selected sample that 1066 showed various degrees of underrepresentation, including 100 cells (10%), 300 cells (30%), 500 1067 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 50 cells (50%), 700 cells (70%), or 900 cells (90%). Three replicates were generated for each degree 1068 of underrepresentation. Replicates were produced with different sample and cell combinations. For 1069 this experiment, we maintained a 18% doublet rate. 1070 1071 WGS for the 80 donors from which the FOUNDIN-PD lines were derived was performed on whole 1072 blood-extracted DNA as previously described by the Parkinson’s Progression Markers Initiative 1073 (PPMI) (28). The c ontrolled-access WGS .vcf files were obtained from https://www.ppmi-1074 info.org/ (accessed 09-17-2023). Genotypes of common variants ( minor allele frequency > 5%) 1075 were used as prior genotype information for the genetic demultiplexing tools in the benchmarking 1076 analyses. 1077 1078 Preparation, processing, and analysis of experimentally pooled samples 1079 Unless specified otherwise, e xperimentally pooled samples were processed with the CellRanger 1080 counts pipeline (v5.0.1) and analyzed with the Seurat (v5.0.0) R package (29), using the 1081 scRNAbox analytical pipeline (30). 1082 1083 Non-small cell lung cancer dataset 1084 NSCLC dissociated tumor cells from seven donors were labelled with TotalSeq-B Human 1085 TBNK Cocktail (18). Multiplexed cells were then sequenced on an Illumina NovaSeq 6000 to 1086 an average read depth of approximately 70,000 reads per cell for gene expression and 25,000 1087 reads per cell for CellPlex. Publicly available gene expression .bam and barcode .tsv files 1088 returned from the CellRanger multi pipeline (v6.1.2) were obtained from the 10X Genomics 1089 Datasets portal (10X Genomics Datasets) and used as input into the Ensemblex pipeline . We 1090 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 51 used the sample -specific gene expression .bam files and the BCFtools (v1.16) mpielup 1091 function to generate genotype likelihoods for prior genotype information (31). 1092 1093 We used HTOdemux to assign the cells back to their donor -of-origin based on the CMO 1094 expression profiles as a proxy for ground-truth sample labels (19). Publicly available feature-1095 barcode expression matrices were filtered to only include CMO labels used for multiplexing 1096 — CMO301, CMO302, CMO303, CMO304, CMO306, CMO307, and CMO308 — and 1097 barcodes with a CMO count > 0. The CMO expression profiles were normalized with Seurat’s 1098 NormalizeData function using the CLR normalization method and HTOdemux was applied to 1099 the CMO assay using a positive quantile of 0.99. 1100 1101 Dopaminergic neuron dataset 1102 Jerber et al. sequenced multiplexed experiments comprising 22 healthy donor iPSC lines from 1103 the HipSci project (32) (http://www.hipsci.org) on days 11, 30, and 52 of DaN differentiation 1104 using Illumina HiSeq 4000 to an average depth of 40,000-60,000 reads per cell (12). We used 1105 three technical replicates for each timepoint, which are comprehensively described in 1106 Additional File 1: Table S3. Publicly available gene expression .fastq files were obtained from 1107 the European Nucleotide Archive (ENA) with accession number ERP121676 and processed 1108 with the CellRanger counts pipeline (v5.0.1) with default parameters using the GRCh37 1109

genome. The CellRanger-generated. bam files, filtered barcode .tsv files, and .vcf 1110 files describing the pooled samples (see below) were used as input into the Ensemblex pipeline 1111 for each technical replicate independently. Filtering of the scRNAseq data was performed as 1112 described by Jerber et al. (12). Genes with non-zero counts in at least 0.05% of cells were 1113 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 52 retained. DoubletFinder (v2.0.4) was applied independently to each technical replicate. Time-1114 point specific replicates were integrated with Seurat’s integration algorithm (33) and clustered 1115 by the Louvain network detection using the top 50 PCs and 10 nearest neighbours. 1116 1117 Whole-exome sequencing (WES) .vcf files corresponding to the 22 pooled HipSci lines were 1118 obtained from the ENA with accession number PRJEB7243 (34). Genotypes of common 1119 variants (minor allele frequency > 1%) were used as prior genotype information for the genetic 1120 demultiplexing tools (12). 1121 1122 Neural stem cell dataset 1123 We performed two multiplexed experiments comprising iPSCs from individuals with ADHD 1124 and heathy controls differentiated into NSCs: Experiment 1 (n ADHD = 7; n control = 6) and 1125 Experiment 2 (n ADHD = 9; n control = 7). 1126 1127 Subject recruitment 1128 Patients diagnosed with ADH D and matching healthy controls between 6−18 years old 1129 were recruited by the Department of Child and Adolescent Psychiatry and Psychotherapy 1130 of the University of Zurich, as described previously (35). Inclusion and exclusion criteria 1131 for recruitment of these individuals described previously (35). Additional File 1: Table 1132 S4 provides a list of the individual subjects and their derived cell lines included in this 1133 study. Salivary DNA from ADHD patients and controls was genotyped using the Infinium 1134 Global Screening Array (Illumina) , as previously described, and used as prior genotype 1135 information for genetic demultiplexing (35). 1136 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 53 1137 Neural stem cell culture 1138 The generation and characterization of iPSC used in this study and the NSCs differentiation 1139 protocols were previously described in (35) (36). NSCs cultures were seeded in two 1140 independent experiments (designated as “1” and “2”), each of them consisting of NSCs 1141 pooled together into two culture dishes and maintained as NSCs until 100% confluence, 1142 when all iPSC lines were combined into one sample for sequencing. For most cell lines 1143 different clones for each iPSC line were used in the two experiments Additional File 1: 1144 Table S5. When applicable, the second clones of the same NSCs lines were cultured 1145 separately (designated as “.1” and “.2”) in a second experiment. In the first experiment, 1146 56,250 cells per cell line were seeded in the pooled dishes. In the second experiment the 1147 proportions of cells seeded we adjusted to their proliferation profile assessed in (36). Upon 1148 reaching 100% confluence, cells were dissociated for scRNAseq experiments and 1149 combined to a single sample for sequencing as described below. 1150 1151 Dissociation of pooled neural stem cell cultures for single-cell RNA sequencing 1152 Cells were washed in PBS and then incubated with 1 mL of StemPro Accutase (Gibco) for 1153 3 minutes at 37°C. After incubation, 2 mL of PBS, stopping the Accutase reaction, and cells 1154 were gently pipetted up and down between 5 to 10 times to break up clumps of cells before 1155 transfer to a 15 mL conical tube. The cells were centrifuged at 300 x g for 5 minutes and 1156 the supernatant was removed. Following, 334 µL of Neural Expansion Media (NEM) was 1157 added to each cell pellet using a 1000 µL pipette tip until cells were completely 1158 resuspended. An additional 666 µL of NEM was added to each well and gently pipette 1159 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 54 mixed 5 times. A 100 -µm cell strainer was used to filter the cell suspension before 1160 centrifugation at 300 x g for 4 minutes. The supernatant was carefully removed, and the 1161 pellet was resuspended in 3 mL of PBS 1x containing 0.04% Bovine Serum Albumin 1162 (BSA) by pipetting up and down 5 times using a 5 mL serological pipette. The cells were 1163 centrifuged at 300 x g for 10 minutes and further submitted to live cell sorting with the 1164 Magnetic Dead Cell Removal Kit (Miltenyi Biotec, 130 -090-101), according to the 1165 manufacturer. The resulting flow-through containing live cells was centrifuged for 300 x g 1166 for 5 minutes and the supernatant was removed carefully to not disturb the cell pellet. Cells 1167 were resuspended in 1 mL of PBS 1x containing 0.04% BSA for automated cell counting. 1168 For each experiment, the cells from the two culture dishes were processed in parallel. Equal 1169 counts of cells were combined for the final cell suspension for scRNAseq preparation at 1170 the Functional Genomics Center Zurich at the University of Zurich. 1171 1172 Library processing and sequencing 1173 All samples were processed using the 10x Genomics Chromium 3’ Single Cell Protocol 1174 and sequenced using NovaSeq 6000 S1 (Illumina). For the first sample containing NSC 1175 pools 1.1 and 1.2, 18,000 NSCs were loaded into one single 10x Genomics Lane to target 1176 13,000 cells. For the second sample containing NSC pools 2.1 and 2.2, 29,000 NSCs were 1177 loaded to target 18,000 cells. 1178 1179 Demultiplexing and scRNAseq analysis 1180 FASTQ files were processed with the CellRanger counts pipeline ( v5.0.1) with default 1181 parameters and aligned to the GRCh37 reference genome. The CellRanger-generated. bam 1182 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 55 files, filtered barcode .tsv files, and .vcf files describing the pooled samples were used as 1183 input into the Ensemblex pipeline. Genotypes of common variants (minor allele frequency 1184 > 1%) were used as prior genotype information for the genetic demultiplexing tools . The 1185 filtered feature-barcode expression matrices were used to analyz e the pooled cells 1186 following a standard scRNAseq analysis workflow using Seurat (30). Cells were filtered 1187 for > 500 total and unique RNA transcripts. Doublets were removed using DoubletFinder 1188 (v2.0.4). The two NSC samples were integrated using Seurat’s integration algorithm (33). 1189 The top 25 PCs were selected for Louvain network detection to identify clusters using 65 1190 nearest neighbours. Twelve clusters were identified at a clustering resolution of 0.25, which 1191 were assigned as eight putative cell types using a combination of known markers and gene 1192 enrichment analysis. The top marker genes from each cluster were identified using Seurat’s 1193 FindAllMarkers with the Wilcoxon rank-sum test. Significant DEGs (log2 fold change > 1194 0.25 and P-value < 0.05 ) were input into EnrichR (37) and cell types were predicted with 1195 the Cell Marker Augmented 2021 (38) and Azimuth Cell Types 2021 (39) libraries. Multiple 1196 clusters showed expression profiles for similar broad cell types — Neurons, NPC s, and 1197 NSCs. We used Seurat’s FindMarkers function to identify differentially expressed marker 1198 genes between the clusters of the same broad cell type and top marker genes were selected 1199 to identify the cell subtypes. 1200 1201 For each putative cell type, DGE was calculated between ADHD and controls using the 1202 MAST statistical framework (22, 40) . Pooled cells were assigned as ADHD or control 1203 based on the demultiplexed sample labels from each of the individual genetic 1204 demultiplexing tools. Cells labeled as “ambiguous singlets” or doublets by the individual 1205 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 56 tools were excluded from their respective DGE analysis. P-values were corrected for 1206 multiple hypothesis testing using the Bonferroni method. A gene was considered 1207 differentially expressed if the adjusted P -value was ≤ 0.0 1 and the absolute value of the 1208 Log2 fold-change was ≥ 0.5. To compute DGE using the sample labels from the individual 1209 tools after the removal of Ensemblex’s putative doublet calls, we repeated the above 1210 procedures but this time all cells labeled as doublets by the respective tool or Ensemblex 1211 were excluded from the DGE analysis. 1212 1213 Performance metrics and statistical analyses 1214 We performed all statistical analyses using the R statistical software (v4.2.2) (41). We used the 1215 ggplot2 R package (v3.4.2) for data visualization (42). 1216 1217 Singlet classification 1218 A singlet was considered correctly classified if the demultiplexed sample label matched the 1219 ground-truth sample label (i.e., specific sample ID) and the assignment probability exceeded 1220 the recommended threshold for the respective tool. For computationally multiplexed pools, the 1221 proportion of correctly classified singlets was computed as: 1222 1223 (7) 𝑃𝑟𝑜𝑝𝑜𝑟𝑡𝑖𝑜𝑛 𝑐𝑜𝑟𝑟𝑒𝑐𝑡 𝑠𝑖𝑛𝑔𝑙𝑒𝑡𝑠= 𝑇𝑃 𝑛 𝑡𝑟𝑢𝑒 𝑠𝑖𝑛𝑔𝑙𝑒𝑡𝑠 1224 1225 For the NSCLC dataset , HTOdemux’s sample labels were considered ground-truth, and the 1226 singlet TP and FP rate were computed as: 1227 1228 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 57 (8) 𝑆𝑖𝑛𝑔𝑙𝑒𝑡 𝑇𝑃 𝑟𝑎𝑡𝑒 = 𝑇𝑃 𝑛 𝐻𝑇𝑂𝑑𝑒𝑚𝑢𝑥 𝑠𝑖𝑛𝑔𝑙𝑒𝑡𝑠 1229 (9) 𝑆𝑖𝑛𝑔𝑙𝑒𝑡 𝐹𝑃 𝑟𝑎𝑡𝑒= 𝐹𝑃 𝑛 𝐻𝑇𝑂𝑑𝑒𝑚𝑢𝑥 𝑠𝑖𝑛𝑔𝑙𝑒𝑡𝑠 1230 1231 Doublet identification 1232 A doublet was considered correctly classified if the demultiplexed sample label matched the 1233 ground-truth sample label, independent of the assignment probability. For computationally 1234 multiplexed pools, the proportion of correctly classified doublets was computed as: 1235 1236 (10) 𝑃𝑟𝑜𝑝𝑜𝑟𝑡𝑖𝑜𝑛 𝑐𝑜𝑟𝑟𝑒𝑐𝑡 𝑑𝑜𝑢𝑏𝑙𝑒𝑡𝑠= 𝑇𝑁 𝑛 𝑡𝑟𝑢𝑒 𝑑𝑜𝑢𝑏𝑙𝑒𝑡𝑠 1237 1238 For the NSCLC dataset, TP doublets were defined as cells classified as doublets by both 1239 HTOdemux and Ensemblex; FP doublets were defined as cells classified as singlets by 1240 HTOdemux and doublets by Ensemblex; FN doublets were defined as cells classified as 1241 doublets by HTOdemux and singlets by Ensemblex. The doublet TP, FP, and FN rates were 1242 computed as: 1243 1244 (11) 𝐷𝑜𝑢𝑏𝑙𝑒𝑡 𝑇𝑃 𝑟𝑎𝑡𝑒 = 𝑇𝑃 𝑛 𝐻𝑇𝑂𝑑𝑒𝑚𝑢𝑥 𝑑𝑜𝑢𝑏𝑙𝑒𝑡𝑠 1245 (12) 𝐷𝑜𝑢𝑏𝑙𝑒𝑡 𝐹𝑃 𝑟𝑎𝑡𝑒= 𝐹𝑃 𝑛 𝑝𝑜𝑜𝑙𝑒𝑑 𝑑𝑟𝑜𝑝𝑙𝑒𝑡𝑠 1246 (13) 𝐷𝑜𝑢𝑏𝑙𝑒𝑡 𝐹𝑁 𝑟𝑎𝑡𝑒= 1 − 𝐷𝑜𝑢𝑏𝑙𝑒𝑡 𝑇𝑃 𝑟𝑎𝑡𝑒 1247 1248 Adjusted Rand Index 1249 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 58 To evaluate the similarity between two distinct sample clusterings we computed the ARI using 1250 the pdfCluster (v1.0.4) R package (43). For the benchmarking analyses, we computed the ARI 1251 between the demultiplexed sample labels by each genetic demultiplexing tool and the ground-1252 truth sample labels (computationally pooled samples) or HTOdemux’s sample labels (NSCLC 1253 dataset). We followed the same procedure when computing the ARI between Ensemblex’s 1254 sample labels and those of its constituent tools (DaN and NSC datasets); however, the ground-1255 truth sample labels were replaced by Ensemblex’s sample labels for these analyses. For 1256 experiments evaluating the impact of doublets on the stability of clusters in gene expression 1257 space, we computed the ARI between clusters at a given clustering resolution after removing 1258 doublets identified by each genetic demultiplexing tool. Clustering stability was computed at 1259 resolutions 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 resolution, 25 iterations of Louvain clustering were performed while shuffling the order of the 1261 nodes in the graph. The ARI between clustering pairs at each clustering resolution was then 1262 computed. 1263 1264 Balanced accuracy 1265 Balanced accuracies were computed to evaluate the binary classification performance of each 1266 genetic demultiplexing tool on imbalanced datasets, where doublets represented a minority 1267 class compared to singlets. The balanced accuracy of each genetic demultiplexing tool was 1268 computed against the ground -truth sample labels (computationally pooled samples ) or 1269 HTOdemux’s sample labels (NSCLC dataset) using equation 1. 1270 1271 Matthew’s correlation coefficient (MCC) 1272 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 59 The MCC was used as a second metric for evaluating the binary classification performance of 1273 the genetic demultiplexing tool. The MCC of each genetic demultiplexing tool was computed 1274 against the ground -truth sample labels (computationally pooled samples ) or HTOdemux’s 1275 sample labels (NSCLC dataset) using equation 14: 1276 1277 (14) 𝑀𝐶𝐶 = 𝑇𝑁×𝑇𝑃−𝐹𝑁×𝐹𝑃 √(𝑇𝑃+𝐹𝑃)(𝑇𝑃−𝐹𝑁)(𝑇𝑁+𝐹𝑃)(𝑇𝑁+𝐹𝑁) 1278 1279 Area under the receiver operating characteristic curve for singlet detection 1280 To evaluate how well each genetic demultiplexing tool’s assignment probability corresponded 1281 to the accuracy of their singlet assignments when ground-truth sample labels were known, we 1282 fit a binary logistic regression model to compute the odds that a singlet was correctly classified 1283 by a tool given the corresponding confidence score or probability. Correctly and incorrectly 1284 classified singlets were set as the positive and negative references, respectively. We then used 1285 the binary logistic regression model to compute the receiver operating characteristic curve for 1286 each tool , which plots the singlet TP and FP rate s across classification thresholds, and 1287 calculated the AUC using the empirical method implemented in the ROCit (v2.1.1) R package 1288 (26). 1289 1290 Abbreviations 1291 ADHD, attention deficit hyperactivity disorder; ANOV A, Analysis of variance; ARI, Adjusted 1292 Rand Index; AUC, area under the receiver operating characteristic curve; BSA, Bovine Serum 1293 Albumin; CMO, Cell Multiplexing Oligonucleotides; DaN, dopaminergic neurons; DGE, 1294 differential gene expression; DEG, differentially expressed genes; ENA, European Nucleotide 1295 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 60 Archive; eQTL, expression quantitative trait loci; FN, false -negative; fNN, nearest neighbour 1296 frequency; FOUNDIN -PD; Foundational Data Initiative for Parkinson’s Disease; FP, false 1297 positive; iPSC, induced pluripotent stem cell; K, kurtosis; MAST, model-based analysis of single-1298 cell transcriptomics; MCC, Matthew’s Correlation Coefficient; nCD, number of confident 1299 doublets; NEM, neural expansion media; NPC, neural progenitor cell; NSC, neural stem cell; 1300 NSCLC, non-small cell lung cancer; PC, principal component; PCA principal component analysis; 1301 PPMI, Parkinson’s Progression Markers Initiative; pT, nearest neighbour percentile threshold; 1302 scATACseq, single-cell assay for transposase-accessible chromatin sequencing; scRNAseq, 1303 single-cell RNA sequencing; SNP, single nucleotide polymorphism; snRNAseq, single -nuclei 1304 RNA sequencing; TN, true-negative; TP, true-positive; UMI, unique molecular identified; WES, 1305 whole-exome sequencing; WGS, whole-genome sequencing. 1306 1307 Declarations 1308 Ethics approval and consent to participate 1309 The iPSC lines (ADHD & controls) used in this project were approved by the Cantonal Ethics 1310 Committee Zurich (BASEC-Nr.-2016-00101 & BASEC -Nr.-201700825) and followed the latest 1311 version of the Declaration of Helsinki, as previously reported (35). The subjects and/or parents 1312 have voluntarily consented to participate in this study. 1313 1314 Consent for publication 1315 Not applicable. 1316 1317 Availability of data and materials 1318 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 61 Transcriptional data for the 80 independently sequenced iPSC lines and the corresponding WGS 1319 data are available from the PPMI database (www.ppmi-info.org/access-dataspecimens/download-1320 data), RRID:SCR 006431. For up -to-date information on the study, visit www.ppmi-info.org. 1321 Processed transcriptional data for the NSCLC dataset are available from the 10X Genomics 1322 Datasets Portal (https://www.10xgenomics.com/datasets/20k-mixture-of-nsclc-dtcs-from-7-1323 donors-3-v3-1-with-intronic-reads-3-1-standard). Transcriptional data for the DaN datasets are 1324 available from the ENA with accession number ERP121676. WES data for the 22 HipSci lines 1325 pooled in the DaN datasets are available from the ENA with accession number PRJEB7243. 1326 Processed scRNAseq data for the NSC dataset are available from the corresponding author upon 1327 reasonable request. The code used for the analyses presented in the work are available at 1328 https://github.com/neurobioinfo/ensemblex. Ensemblex is freely available under an MIT open -1329 source license at https://zenodo.org/records/11639103. 1330 1331 Competing interests 1332 The authors declare that they have no competing interests. 1333 1334 Funding 1335 This work was supported by the Michael J. Fox Foundation [MJFF-021629 to EAF, RAT, and 1336 SMKF]. PPMI — a public-private partnership — is funded by the Michael J. Fox Foundation for 1337 Parkinson’s Research and funding partners, including 4D Pharma, Abbvie, AcureX, Allergan, 1338 Amathus Therapeutics, Aligning Science Across Parkinson's, AskBio, Avid Radiopharmaceuticals, 1339 BIAL, BioArctic, Biogen, Biohaven, BioLegend, BlueRock Therapeutics, Bristol -Myers Squibb, 1340 Calico Labs, Capsida Biotherapeutics, Celgene, Cerevel Therapeutics, Coave Therapeutics, 1341 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 62 DaCapo Brainscience, Denali, Edmond J. Safra Foundation, Eli Lilly, Gain Therapeutics, GE 1342 HealthCare, Genentech, GSK, Golub Capital, Handl Therapeutics, Insitro, Jazz Pharmaceuticals, 1343 Johnson & Johnson Innovative Medicine, Lundbeck, Merck, Meso Scale Discovery, Mission 1344 Therapeutics, Neurocrine Biosciences, Neuron23, Neuropore, Pfizer, Piramal, Prevail 1345 Therapeutics, Roche, Sanofi, Servier, Sun Pharma Advanced Research Company, Takeda, Teva, 1346 UCB, Vanqua Bio, Verily, V oyager v. 25MAR2024 Therapeutics, the Weston Family Foundation 1347 and Yumanity Therapeutics. For funding t he ADHD study, we thank the Neuroscience Centre 1348 Zurich (ZNZ) for the Zurich -McGill University Neurodevelopmental Disorder Research 1349 Collaboration and the Psychiatric University Hospital Zurich (PUK) Forschungsfonds Nr. 8702 1350 “Fonds für wissenschaftliche Zwecke im Interesse der Heilung von psychiatrischen Krankheiten” 1351 and the Candoc PhD grant from the University of Zurich [FK-22-044 to CMYO]. 1352 1353 Authors’ contributions 1354 MRF, EAF, RAT, and SMKF conceived the study. MRF developed the Ensemblex framework and 1355 wrote the corresponding R code. MRF performed the analyses and produced the figures. MRF and 1356 SA developed the Ensemblex pipeline and created the GitHub site. MRF and SA tested the 1357 Ensemblex pipeline. MRF wrote the Ensemblex documentation. MRF, AAD, RAT and SMKF 1358 interpreted the data sets. CMYO performed all cell cultures and sequencing preparation for the 1359 NSC dataset. MRF, CMYO, and RAT performed the cell type annotations for the NSC dataset. LS 1360 and EG provided the NSC genetic data. SW recruited the subjects for the NSC dataset. MRF wrote 1361 the manuscript with input from all authors. EG supervised the NSC data collection. RAT and 1362 SMKF supervised the project. 1363 1364 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 63 Acknowledgments 1365 The author’s acknowledge Dan Spiegelman for his help in processing the VCF files for individuals 1366 pooled in the NSC dataset. Schematic illustrations presented in the manuscript were prepared with 1367 BioRender (https://www.biorender.com/). 1368 1369 Authors’ information 1370 This work was supported by the Michael J Fox Foundation in a grant award to EAF , RAT, and 1371 SMKF [MJFF-021629]. MRF is supported by a CIHR Canada Graduate Scholarships -Master’s 1372 Award, a Fonds de Recherche Santé Québec Master’s Award, and a Brain Canada Rising Stars 1373 Award. EAF is supported by a Fonds d’Accéleration des Collaborations en Santé (FACS) grant 1374 from CQDM/MEI and a Canada Research Chair (Tier 1) in Parkinson’s disease. R.A.T. received 1375 funding through the McGill Healthy Brains for Healthy Lives (HBHL) Postdoctoral Fellowship 1376 and Molson NeuroEngineering Fellowship. SMKF received funding from Brain Canada and the 1377 Montreal Neurological Institute-Hospital. CMYO is supported by the Candoc PhD grant from the 1378 University of Zurich (UZH) [FK-22-044]. 1379 1380 Figure legends 1381 Figure 1. Evaluation of existing individual genetic demultiplexing tools. Evaluation of genetic 1382 demultiplexing tools with prior genotype information on 96 in silico pools with known ground -1383 truth sample labels ranging in size from 4 to 80 multiplexed induced pluripotent stem cell (iPSC) 1384 lines from genetically distinct individuals, averaging 17,396 cells per pool and a 15% doublet rate. 1385 A) Line graphs showing the proportion of correctly classified singlets, doublets, and all cells by 1386 each individual genetic demultiplexing tool across varying numbers of multiplexed iPSC lines in 1387 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 64 a single pool (sample number). The large dots show the mean proportion of correct classifications 1388 by an individual tool across replicates at a given sample size (n = 9 per pool size). The blue points 1389 show the proportion of cells that were correctly classified by at least one individual genetic 1390 demultiplexing tool: Demuxalot, Demuxlet, Souporcell, or Vireo -GT. B) Bar chart showing the 1391 mean proportion of total cells from an individual pool correctly classified by only one genetic 1392 demultiplexing tool. Error bars represent one standard deviation from the mean. (n = 9 per pool 1393 size) C) Bar chart showing the proportion of correctly classified singlet cells labelled as 1394 “unassigned” (ambiguous singlet assignments) due to assignment probabilities below the 1395 recommended threshold of the respective genetic demultiplexing tool. Error bars represent one 1396 standard deviation from the mean. (n = 9 per pool size). 1397 1398 Figure 2. Characterization of the Ensemblex framework. Ensemblex is a probabilistic -1399 weighted ensemble genetic demultiplexing framework for s ingle-cell RNA sequencing analysis, 1400 which was designed to leverage the most probable sample labels from each of its constituent tools: 1401 Demuxalot, Demuxlet, Souporcell, and Vireo when using prior genotype information or 1402 Demuxalot, Freemuxlet, Souporcell, and Vireo when prior genotype information is not available. 1403 A) The Ensemblex workflow begins with demultiplexing pooled cells from genetically distinct 1404 individuals by each of the constituent tools. The outputs from each individual demultiplexing tool 1405 are then used as input into the Ensemblex framework. B) The Ensemblex framework comprises 1406 three distinct steps that are assembled into a pipeline: 1) accuracy-weighted probabilistic ensemble, 1407 2) graph -based doublet detection, and 3) ensemble -independent doublet detection. C-D) Line 1408 graphs showng t he contribution of each step of the Ensemblex framework on 96 in silico pools 1409 with known ground -truth sample labels ranging in size from 4 to 80 multiplexed induced 1410 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 65 pluripotent stem cell (iPSC) lines from genetically distinct individuals, averaging 17,396 cells per 1411 pool and a 15% doublet rate. The average proportion of correctly classified C) singlets and D) 1412 doublets across replicates at a given pool size is shown after sequentially applying each step of the 1413 Ensemblex framework with prior genotype information (n = 9 per pool size) . The right panels 1414 show the average proportion of correct classifications across all 96 pools; error bars represent one 1415 standard deviation from the mean. The blue points show the proportion of cells that were correctly 1416 classified by at least one individual genetic demultiplexing tool: Demuxalot, Demuxlet, 1417 Souporcell, or Vireo-GT. 1418 1419 Figure 3. Ensemblex ground-truth benchmarking on computationally multiplexed pools. The 1420 genetic demultiplexing tools with prior genotype information were evaluated on 96 in silico pools 1421 with known ground -truth sample labels ranging in size from 4 to 80 multiplexed induced 1422 pluripotent stem cell (iPSC) lines from genetically distinct individuals, averaging 17,396 cells per 1423 pool and a 15% doublet rate. A singlet was considered correctly classified if the assigned sample 1424 label matched the ground -truth sample label and the assignment probability exceeded the 1425 recommended threshold for the respective tool; a doublet was considered correctly classified if the 1426 assigned sample label matched the ground -truth sample label, regardless of the assignment 1427 probability. A-I) Line graphs showing the performance of Ensemblex and the individual genetic 1428 demultiplexing tools across evaluation metrics. The large dots show the mean value for each tool 1429 across replicates at a given sample size (n = 9 per pool size). A) Proportion of correctly classified 1430 singlets. B) Proportion of correctly classified doublets. C) Proportion of correctly classified cells. 1431 D) Adjusted Rand Index between each tool’s sample labels and the ground-truth sample labels. E) 1432 Balanced accuracy of each tool. F) Matthew’s Correlation Coefficient of each tool. G) Area under 1433 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 66 the receiver operating characteristic curve (AUC) of the singlet assignment probability for each 1434 tool. H) Proportion of usable cells returned by each tool. Usable cells were defined as cells 1435 classified by singlets with an assignment probability exceeding the recommended threshold of the 1436 respective tool. I) Error rate amongst the usable cells returned by each tool; erroneous 1437 classifications comprised of true doublets labeled as singlets or true singlets assigned to the wrong 1438 sample. 1439 1440 Figure 4. Evaluating Ensemblex on experimentally multiplexed cells using donor -specific 1441 oligonucleotide labels as a proxy for ground -truth. Non-small cell lung cancer (NSCLC) 1442 dissociated tumor cells from 7 individuals were pooled and labelled with donor-specific 1443 oligonucleotide-labels. Cells were demultiplexed according to their expression of donor-specific 1444 oligonucleotide labels by HTOdemux ; HTOdemux’s sample labels were used as a proxy for 1445 ground truth. True positives (TP) singlets were defined as cells classified as singlets by both 1446 HTOdemux and Ensemblex with matching sample labels; false positives (FP) singlets were 1447 defined as cells classified as singlets by both HTOdemux and Ensemblex but assigned to different 1448 donors. TP doublets were defined as cells classified as doublets by both HTOdemux and 1449 Ensemblex; FP doublets were defined as cells classified as singlets by HTOdemux and doublets 1450 by Ensemblex; false negatives (FN) doublets were defined as cells classified as doublets by 1451 HTOdemux and singlets by Ensemblex. A) T-distributed Stochastic Neighbor Embedding (t-SNE) 1452 visualization of HTOdemux’s sample labels. B) T-SNE visualization of Ensemblex’s 1453 demultiplexing performance using HTOdemux’s sample labels as ground truth for singlets (left) 1454 and doublets (right) . C) Bar plots showing the s inglet TP and FP rate s for each genetic 1455 demultiplexing tool using HTOdemux’s sample labels as ground truth . D) Bar plots showing the 1456 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 67 doublet TP and FP rates for each genetic demultiplexing tool using HTOdemux’s sample labels as 1457 ground truth. E) Scatter plot showing the proportion of usable cells (confidently classified singlets) 1458 and the corresponding usable cell error rate for each genetic demultiplexing tool. F) Adjusted Rand 1459 Index, balanced accuracy, Matthew’s Correlation Coefficient, and area under the receiver operating 1460 characteristic curve (AUC) of the singlet assignment probability for each genetic demultiplexing 1461 tool. 1462 1463 Figure 5. Application of Ensemblex to highly multiplexed, experimentally pooled cultures of 1464 differentiated dopaminergic neurons. A) Time line of iPSC pooling, dopaminergic neuron 1465 (DaN) differentiation, and sample collection from the DaN dataset by Jerber et al. (12). Three 1466 technical replicates at each time point (days 11, 30 and, 52 of differentiation) from pools containing 1467 22 individual iPSC lines were used in the analysis. Across all timepoints and technical replicates, 1468 84,746 cells were obtained for analysis. B) Uniform manifold approximation and projection 1469 (UMAP) plots showing confidently assigned singlets or predicted doublets (blue) and ambiguous 1470 singlets (singlet assignments with insufficient assignment probabilities ; red ) returned by each 1471 demultiplexing tool. C) Stacked bar chart showing the proportion of confidently assigned singlets 1472 or predicted doublets (blue) and ambiguous singlets (red) across technical replicates at each time 1473 point returned by each demultiplexing tool. D) Boxplot showing the p roportion of confidently 1474 classified singlets across technical replicates and time points by each demultiplexing tool. 1475 Wilcoxon rank-sum tests were used to compare the proportion of confidently classified singlets by 1476 Ensemblex to that of its constituents (n = 9 pools). E) Bar chart showing the proportion of 1477 overlapping ambiguous singlet assignments amongst demultiplexing tools across technical 1478 replicates and time points (n = 9 pools) . F) Boxplot showing the Adjusted Rand Index (ARI) 1479 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 68 assessing cluster stability across a range of 11 clustering resolutions (n clustering iterations = 25) 1480 after removing doublets identified by each demultiplexing tool. Wilcoxon rank -sum tests were 1481 used to compare the clustering ARI after removing Ensemblex doublets to the clustering ARI after 1482 removing doublets identified by each constituent tool. * Adjusted P-value < 0.05; ** adjusted P-1483 value < 0.01; *** adjusted P-value < 0.001 1484 1485 Figure 6. Evaluating the impact of discordant assignments between genetic demultiplexing 1486 tools on differential gene expression analysis. A) Schematic illustrating the workflow for the 1487 neural stem cell (NSC) dataset. Pooled induced pluripotent stem cell (iPSC) -derived neural stem 1488 cell cultures from individuals with attention deficit hyperactivity disorder (ADHD) and controls 1489 were collected in two separate experiments . NSCs were dissociated for single-cell RNA 1490 sequencing and prior genotype information of the pooled subjects was obtained through 1491 microarray genotyping. The pools were demultiplexed by Ensemblex and its constituent s with 1492 prior genotype information and differential gene expression (DEG) was computed between ADHD 1493 and controls. B) Uniform manifold approximation and projection (UMAP) plot showing the 1494 putative cell types. C) Summary of the number of usable cells — singlets above the recommended 1495 probability threshold of the respective demultiplexing tool — assigned to ADHD donors and 1496 controls and the number of identified doublets by each demultiplexing tool. D) Boxplot showing 1497 the Adjusted Rand Index (ARI) assessing cluster stability across a range of 11 clustering 1498 resolutions (n clustering iterations = 25) after removing doublets identified by each demultiplexing 1499 tool. A one-way Analysis of Variance (ANOV A) test comparing the ARI after removing doublets 1500 identified by each tool revealed a significant difference between tools (n = 11 clustering 1501 resolutions; P-value = 1.18e-3). E) Proportion of ADHD and control cells identified as putative 1502 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 69 doublets by Ensemblex that were assigned as singlets by the constituent demultiplexing tools. F) 1503 Heatmap showing the number of cell-type specific DEGs between ADHD and controls using the 1504 subject labels of each demultiplexing tool. G) Heatmap showing the number of cell-type specific 1505 DEGs between ADHD and controls using the subject labels of each demultiplexing tool and 1506 removing putative doublets identified by Ensemblex. Cell-types not shown in the heatmaps had no 1507 DEGs passing the adjusted P-value = 0.5| threshold across all tools. 1508 1509 Tables 1510 Table 1. Summary of individual genetic demultiplexing tools. 1511 Genetic demultiplexing tool Prior genotype information for genetic demultiplexing Included in the Ensemblex framework Demuxalot (5) Required Yes Demuxlet (6) Required Yes Freemuxlet (6) Not supported Yes ScSplit (7) Optional No Souporcell (8) Optional Yes Vireo (9) Optional Yes 1512 Table 2. Application of Ensemblex to pooled cultures of dopaminergic neurons from 22 1513 healthy controls. 1514 ARI between Ensemblex and constituent tool assignments Percent contribution to Ensemblex assignments n usable cells n doublets Day 11 Day 30 Day 52 Day 11 Day 30 Day 52 Demuxalot 0.987 0.955 0.982 97.29% 94.75% 97.57% 75,962 8,279 Demuxlet 0.928 0.062 0.884 95.91% 29.74% 90.55% 57,567 6,614 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 70 Souporcell 0.883 0.876 0.912 91.62% 91.82% 93.84% 76,811 7,740 Vireo-GT 0.961 0.879 0.958 95.95% 88.80% 95.16% 75,933 6,115 Ensemblex NA NA NA NA NA NA 76,222 8,307 DoubletFinder NA NA NA NA NA NA NA 4,597 Pooled cultures of induced pluripotent stem cell (iPSC) lines from 22 healthy donors were 1515 differentiated towards a dopaminergic neuron (DaN) fate and sequenced on days 11, 30, and 52 of 1516 differentiation by Jerber et al. (12). For the analysis we used three technical replicates for each 1517 sequencing timepoint. Each pool was demultiplexed independently by Ensemblex and its 1518 constituent tools with prior genotype information. The Adjusted Rand Index (ARI) between 1519 Ensemblex’s assignments and th ose of the constituent tools was computed across technical 1520 replicates corresponding to each differentiation timepoint. The percent contribution represents the 1521 proportion of assignments from each constituent tool that matched Ensemblex’s assignments. 1522 Usable cells were defined as singlet classifications whose assignment probability exceeded the 1523 recommended threshold of the respective tool. Abbreviations: NA = Not applicable. 1524 1525

1526 1. Fiorini MR, Dilliott AA, Thomas RA, Farhan SMK. Transcriptomics of Human Brain Tissue 1527 in Parkinson's Disease: a Comparison of Bulk and Single-cell RNA Sequencing. Mol Neurobiol. 1528 2024. 1529 2. Ringman JM, Goate A, Masters CL, Cairns NJ, Danek A, Graff-Radford N, et al. Genetic 1530 heterogeneity in Alzheimer disease and implications for treatment strategies. Curr Neurol 1531 Neurosci Rep. 2014;14(11):499. 1532 3. McKinney CE. Using induced pluripotent stem cells derived neurons to model brain 1533 diseases. Neural Regen Res. 2017;12(7):1062-7. 1534 4. Howitt G, Feng Y , Tobar L, Vassiliadis D, Hickey P , Dawson MA, et al. Benchmarking single-1535 cell hashtag oligo demultiplexing methods. NAR Genomics and Bioinformatics. 1536 2023;5(4):lqad086. 1537 5. Rogozhnikov A, Ramkumar P , Shah K, Bedi R, Kato S, Escola GS. Demuxalot: scaled up 1538 genetic demultiplexing for single-cell sequencing. bioRxiv. 2021:2021.05. 22.443646. 1539 6. Kang HM, Subramaniam M, Targ S, Nguyen M, Maliskova L, McCarthy E, et al. 1540 Multiplexed droplet single-cell RNA-sequencing using natural genetic variation. Nat Biotechnol. 1541 2018;36(1):89-94. 1542 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 71 7. Xu J, Falconer C, Nguyen Q, Crawford J, McKinnon BD, Mortlock S, et al. Genotype-free 1543 demultiplexing of pooled single-cell RNA-seq. Genome Biol. 2019;20(1):290. 1544 8. Heaton H, Talman AM, Knights A, Imaz M, Gaffney DJ, Durbin R, et al. Souporcell: robust 1545 clustering of single-cell RNA-seq data by genotype without reference genotypes. Nat Methods. 1546 2020;17(6):615-20. 1547 9. Huang Y , McCarthy DJ, Stegle O. Vireo: Bayesian demultiplexing of pooled single-cell 1548 RNA-seq data without genotype reference. Genome Biol. 2019;20(1):273. 1549 10. Neavin D, Senabouth A, Hang Lee JT , Ripoll A, Consortium s-e, Franke L, et al. Demuxafy: 1550 Improvement in droplet assignment by integrating multiple single-cell demultiplexing and 1551 doublet detection methods. BioRxiv. 2022:2022.03. 07.483367. 1552 11. Cardiello JF, Joven Araus A, Giatrellis S, Helsens C, Simon A, Leigh ND. Evaluation of 1553 genetic demultiplexing of single-cell sequencing data from model species. Life Sci Alliance. 1554 2023;6(8). 1555 12. Jerber J, Seaton DD, Cuomo AS, Kumasaka N, Haldane J, Steer J, et al. Population-scale 1556 single-cell RNA-seq profiling across dopaminergic neuron differentiation. Nature genetics. 1557 2021;53(3):304-12. 1558 13. Jindal A, Gupta P , Jayadeva, Sengupta D. Discovery of rare cells from voluminous single 1559 cell expression data. Nature communications. 2018;9(1):4719. 1560 14. Bressan E, Reed X, Bansal V, Hutchins E, Cobb MM, Webb MG, et al. The Foundational 1561 Data Initiative for Parkinson Disease: Enabling efficient translation from genetic maps to 1562 mechanism. Cell Genom. 2023;3(3):100261. 1563 15. Cardiello JF, Araus AJ, Giatrellis S, Helsens C, Simon A, Leigh ND. Evaluation of genetic 1564 demultiplexing of single-cell sequencing data from model species. Life Science Alliance. 1565 2023;6(8). 1566 16. Large J, Lines J, Bagnall A. A probabilistic classifier ensemble weighting scheme based on 1567 cross-validated accuracy estimates. Data Min Knowl Discov. 2019;33(6):1674-709. 1568 17. Chicco D, Tötsch N, Jurman G. The Matthews correlation coefficient (MCC) is more 1569 reliable than balanced accuracy, bookmaker informedness, and markedness in two-class 1570 confusion matrix evaluation. BioData mining. 2021;14:1-22. 1571 18. 20k Mixture of NSCLC DTCs from 7 donors, 3' v3.1 (with intronic reads) [Internet]. 10X 1572 Genomics. 2022 [cited January 8th, 2024]. Available from: 1573 https://www.10xgenomics.com/datasets/20k-mixture-of-nsclc-dtcs-from-7-donors-3-v3-1-with-1574 intronic-reads-3-1-standard. 1575 19. Stoeckius M, Zheng S, Houck-Loomis B, Hao S, Yeung BZ, Mauck WM, 3rd, et al. Cell 1576 Hashing with barcoded antibodies enables multiplexing and doublet detection for single cell 1577 genomics. Genome Biol. 2018;19(1):224. 1578 20. McGinnis CS, Murrow LM, Gartner ZJ. DoubletFinder: Doublet Detection in Single-Cell 1579 RNA Sequencing Data Using Artificial Nearest Neighbors. Cell Syst. 2019;8(4):329-37 e4. 1580 21. Macosko EZ, Basu A, Satija R, Nemesh J, Shekhar K, Goldman M, et al. Highly Parallel 1581 Genome-wide Expression Profiling of Individual Cells Using Nanoliter Droplets. Cell. 1582 2015;161(5):1202-14. 1583 22. Finak G, McDavid A, Yajima M, Deng J, Gersuk V, Shalek AK, et al. MAST: a flexible 1584 statistical framework for assessing transcriptional changes and characterizing heterogeneity in 1585 single-cell RNA sequencing data. Genome Biol. 2015;16:278. 1586 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 72 23. Weber LM, Hippen AA, Hickey PF, Berrett KC, Gertz J, Doherty JA, et al. Genetic 1587 demultiplexing of pooled single-cell RNA-sequencing samples in cancer facilitates effective 1588 experimental design. Gigascience. 2021;10(9). 1589 24. Bose A, Beal MF. Mitochondrial dysfunction in Parkinson's disease. Journal of 1590 neurochemistry. 2016;139:216-31. 1591 25. Wickham H. ggplot2. Wiley interdisciplinary reviews: computational statistics. 1592 2011;3(2):180-5. 1593 26. Khan MRAA. Rocit-an r package for performance assessment of binary classifier with 1594 visualization. 2019. 1595 27. Genomes Project C, Auton A, Brooks LD, Durbin RM, Garrison EP , Kang HM, et al. A 1596 global reference for human genetic variation. Nature. 2015;526(7571):68-74. 1597 28. Parkinson Progression Marker I. The Parkinson Progression Marker Initiative (PPMI). Prog 1598 Neurobiol. 2011;95(4):629-35. 1599 29. Hao Y , Stuart T , Kowalski MH, Choudhary S, Hoffman P , Hartman A, et al. Dictionary 1600 learning for integrative, multimodal and scalable single-cell analysis. Nature biotechnology. 1601 2023:1-12. 1602 30. Thomas RA, Fiorini MR, Amiri S, Fon EA, Farhan SM. ScRNAbox: Empowering Single-Cell 1603 RNA Sequencing on High Performance Computing Systems. bioRxiv. 2023:2023.11. 13.566851. 1604 31. Li H. A statistical framework for SNP calling, mutation discovery, association mapping 1605 and population genetical parameter estimation from sequencing data. Bioinformatics. 1606 2011;27(21):2987-93. 1607 32. Kilpinen H, Goncalves A, Leha A, Afzal V, Alasoo K, Ashford S, et al. Common genetic 1608 variation drives molecular heterogeneity in human iPSCs. Nature. 2017;546(7658):370-5. 1609 33. Stuart T , Butler A, Hoffman P , Hafemeister C, Papalexi E, Mauck WM, 3rd, et al. 1610 Comprehensive Integration of Single-Cell Data. Cell. 2019;177(7):1888-902 e21. 1611 34. Streeter I, Harrison PW, Faulconbridge A, The HipSci C, Flicek P , Parkinson H, et al. The 1612 human-induced pluripotent stem cell initiative-data resources for cellular genetics. Nucleic 1613 Acids Res. 2017;45(D1):D691-D7. 1614 35. Yde Ohki CM, Grossmann L, Doring C, Hoffmann P , Herms S, Werling AM, et al. 1615 Generation of integration-free induced pluripotent stem cells from healthy individuals. Stem 1616 Cell Res. 2021;53:102269. 1617 36. Yde Ohki CM, Walter NM, Bender A, Rickli M, Ruhstaller S, Walitza S, et al. Growth rates 1618 of human induced pluripotent stem cells and neural stem cells from attention-deficit 1619 hyperactivity disorder patients: a preliminary study. J Neural Transm (Vienna). 2023;130(3):243-1620 52. 1621 37. Chen EY , Tan CM, Kou Y , Duan Q, Wang Z, Meirelles GV, et al. Enrichr: interactive and 1622 collaborative HTML5 gene list enrichment analysis tool. BMC Bioinformatics. 2013;14:128. 1623 38. Zhang X, Lan Y , Xu J, Quan F, Zhao E, Deng C, et al. CellMarker: a manually curated 1624 resource of cell markers in human and mouse. Nucleic Acids Res. 2019;47(D1):D721-D8. 1625 39. Hao Y , Hao S, Andersen-Nissen E, Mauck WM, 3rd, Zheng S, Butler A, et al. Integrated 1626 analysis of multimodal single-cell data. Cell. 2021;184(13):3573-87 e29. 1627 40. Macosko EZ, Basu A, Satija R, Nemesh J, Shekhar K, Goldman M, et al. Highly parallel 1628 genome-wide expression profiling of individual cells using nanoliter droplets. Cell. 1629 2015;161(5):1202-14. 1630 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint 73 41. Ihaka R, Gentleman R. R: a language for data analysis and graphics. Journal of 1631 computational and graphical statistics. 1996;5(3):299-314. 1632 42. Wickham H, Wickham H. Data analysis: Springer; 2016. 1633 43. Azzalini A, Menardi G. Clustering via nonparametric density estimation: The R package 1634 pdfCluster. arXiv preprint arXiv:13016559. 2013. 1635 1636 .CC-BY 4.0 International licenseavailable under a (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 The copyright holder for this preprintthis version posted June 19, 2024. ; https://doi.org/10.1101/2024.06.17.599314doi: bioRxiv preprint