scRNAseq analysis system

Abstract

The present invention provides a system for assessing gene expression levels from scRNA-Seq experiments. The system uses an in-silico cell calling tool to distinguish between partitions containing cells and background partitions that did not fully capture cells. The system also provides a barcode processing tool that tests barcodes from at least the cell and UMI tiers against a whitelist and optionally converts them to a compact index format. Additionally, the system can implement a multi-mapping tool to store information regarding reads that map to multiple locations within a reference.

Description

Technical Field

[0001] The present disclosure relates to expression analysis using methods such as single-cell RNA-seq.

Background Art

[0002] Single-cell RNA sequencing (scRNA-seq) refers to various protocols that involve sequencing RNA from cells and, in most embodiments, providing a measure of gene expression levels from the sequence data. Some approaches to scRNA-Seq rely on isolating cells in droplets that have the potential to assay a large number of cells per experiment. Common droplet-based protocols include Drop-seq (described in Macosko, 2015, Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets, Cell 161(5):1202-14, incorporated by reference) and inDrop (see Klein, 2015, Droplet barcoding for single-cell transcriptomics applied to embryonic stem cells, Cell 161(5):1187-201, incorporated by reference).

[0003] These droplet-based approaches to scRNA-Seq generally use two oligonucleotide barcodes that are added to RNA molecules or their cDNA copies. The sequences of these two barcodes appear in the sequence data obtained from cellular RNA. Typically, one of these barcodes functions as a cell barcode that is the same across an entire droplet but different between droplets, providing a cell identification tag common to all molecules from one cell. Using cell barcodes enables pooling nucleic acids from multiple cells during sequencing and subsequent in silico separation of the sequence data by cell. Another barcode often used in scRNA-Seq is called a unique molecular identifier (UMI). UMIs tag each individual molecule with a unique sequence of nucleobases. When those molecules are amplified, for example, by PCR, amplicons from one molecule share a common UMI sequence. After sequencing, sequence reads containing the same gene sequence data and UMI sequence data can be treated as duplicates. In a process sometimes called read deduplication, all identical reads are treated as one, and after deduplication, unique reads are counted. The count of unique reads from a gene is a count of the number of mRNA transcripts from that gene present in the cell. The number of mRNA transcripts from various genes in a cell provides a measure of the gene expression level for that cell.

[0004] Several factors make it difficult to provide a measure of gene expression levels from sequence data. For example, barcodes are susceptible to errors that occur during sequencing and amplification. In addition, much of the sequence data is derived from amplification of free-floating nucleic acids in droplets that either do not contain cells or contain damaged cells or other debris. SUMMARY OF THE INVENTION MEANS FOR SOLVING THE PROBLEM [[ID=ll]]

[0005] The present invention provides a system for evaluating gene expression levels from scRNA-Seq experiments. The system uses in-silico tools for "cell calling", i.e., for distinguishing partitions (e.g., droplets) containing cells from background partitions that did not fully capture a cell. Additionally, the system provides tools for "barcode processing" where barcodes are optionally converted to a compact format and / or tested against a whitelist in a manner that significantly improves the processing speed during deduplication while preserving information. Further, the system may implement tools for "multi-mapping" that store information regarding reads that map to multiple locations within a reference (such reads were conventionally typically simply discarded). Additionally, the system includes analytical tools for analysis, such as clustering, visualization, classification, and metric extraction, that are useful for indicating which clusters of cells are present in a sample, e.g., based on gene expression patterns.

[0006] These tools for cell calling, barcoding, multimapping, and analysis are particularly suitable for scRNA-Seq using particle-templated instant partition (PIP) that essentially isolates a very large number of cells into emulsion droplets simultaneously. To perform scRNA-Seq using PIP, cells and particles are mixed together in an aqueous mixture. The particles can be hydrogel beads, and these hydrogel beads can carry reagents such as enzymes and / or oligonucleotides. The reaction mixture can be provided in the aqueous phase together with any of various reagents, or on or inside the hydrogel bead particles. These reagents can include, in addition to any enzyme or oligonucleotide, dNTPs, cofactors, dyes, salts, chemicals, others, or combinations thereof. The aqueous mixture will typically be provided with several hydrogel bead particles corresponding to several aqueous partitions known as droplets or "PIPs" to be created. Cells can be introduced by dilution at a concentration such that PIPs with an average number of cells per PIP of 1 or less will be created. Once the aqueous mixture is prepared, oil can be added to the aqueous phase (optionally while providing a surfactant to the mixture). Then, the oil / aqueous mixture can be sheared, for example, by vortexing a tube (such as a 0.5 mL microcentrifuge tube or a larger conical tube). The shearing energy causes each of the hydrogel bead particles to function as a template, almost instantaneously surrounding partitions of the aqueous liquid with oil and forming a plurality of stable monodisperse emulsion droplets, also known as partitions, and thus, particle-templated instant partitions (PIPs). Any number of PIPs can be formed simultaneously and almost instantaneously during vortexing in one reaction volume.

[0007] Almost without exception, each PIP contains one particle, at most one cell, a cell lysis reagent (such as a proteinase enzyme, etc.), and an oligo for RNA hybrid capture, such as a barcoded polyT capture oligo filled within the particle or covalently bound to the particle. The PIP may further contain other reagents such as reverse transcriptase, dNTP, other capture oligos (such as template switching oligos) that are either linked to the particle or free in solution, primers, cofactors, etc. If the partitions or PIPs are formed immediately, each will contain, on average, one or zero cells. The tubes can be transferred to a thermocycler or other temperature control device. The provided reagents can be used to lyse the cells to release the RNA, or the RNA can be captured by the oligos (generally by hybrid capture and optionally by enzyme binding, such as by ligase or transposase). Typically, all of the oligos of one hydrogel bead particle share a copy of a common "cell barcode" that is (essentially) unique to that particle. In addition, each oligo may contain a UMI that is unique (or at least nearly unique, or essentially or effectively unique) among the oligos of that particle. After the RNA molecules are captured by the oligos, typical RNA-Seq proceeds by cDNA synthesis, amplification and ligation of sequencing adapters, and sequencing to generate sequence data containing nucleotide sequences derived from the RNA molecules and from the oligos (such as from the barcodes). Systems for barcode processing, cell calling, multi-mapping, and analysis using such sequence data are described herein.

[0008] Barcode processing by the system of the present invention may use a hierarchical barcode system, and optionally, all levels include one of a specified list of possible barcodes. The levels may refer to levels of information specific to the barcode, such as cell barcodes and then UMIs. There may also be intermediate-level barcodes, such as high-level barcodes for different patients or samples, and / or barcodes specific to genes (e.g., when sequence-specific primers can be used). Typically, scRNA-Seq includes cell and UMI barcodes, but other levels may be included for any category of information, such as date or location. Including a barcode in a reagent oligonucleotide leads to sequence data that includes the sequence from the barcode. A computer system may match the barcode to a list or database to assign sequence reads to cells or individual molecules. Barcode matching can be performed very efficiently by matching each level in isolation. Some embodiments of the computer system of the present invention allow for a Hamming distance of 1. Limiting the Hamming distance to 1 means that the bases in a sequence read, e.g., "read 1" ("R1") in a FASTQ file corresponding to the position of each level, can differ from the barcode by only one base and still match that barcode. Hierarchical barcodes can be replaced by a new set of generated barcodes to simplify the barcode list sent to a reference mapping module such as STAR and reduce its size. For example, each barcode can be replaced by a number, such as its numerical position on the list of barcodes. In some embodiments, new barcodes are generated only for barcodes that have at least one matching read. That is, if trillions of barcodes are calculated to form the input to an experiment and the system processes the resulting FASTQ file and identifies only billions of barcodes in the reads, a new restricted list of barcodes can be assigned new barcodes with a significant savings in downstream computational resources.The barcode processing method helps keep downstream steps computationally tractable while maintaining the ability to provide a measure of expression level when mapping reads to references and deduplicating reads using UMI counts.

[0009] Compared to a microfluidic platform (where a limited number of cells are successively isolated into partitions), the use of PIP (where multiple cells are simultaneously isolated into aqueous partitions) reveals several issues associated with the throughput and amount of data initially achieved by PIP. For example, using PIP reveals that it is beneficial to have an automated and accurate tool, aka cell calling, to correctly call whether a partition contains cells as opposed to background partitions that do not fully capture cells. The present disclosure recognizes that, as a practical matter, in the amount of data and throughput associated with PIP, it can be automatically done to distinguish partitions containing cells from background partitions (that did not fully capture cells). For example, from sequence data, all of the identified barcodes can be ranked in order by the number of UMIs associated with each barcode. Barcodes associated with a relatively large number of UMIs are more likely to originate from cells than background partitions, while barcodes associated with a relatively small number of UMIs are more likely to originate from background partitions. The system of the present invention can present the user with a sensitivity / specificity selection (e.g., select one of 5 preset options). Behind the scenes, the computer system calculates the number of UMIs (e.g., on the y-axis) across the ranks of the ranked barcodes (x-axis), fits a function to the curve, finds the derivative to find the inflection points, and inputs the selections from those points along the curve into the user interface. The user can select a sensitivity level for cell calling and the analysis proceeds at the selected sensitivity. For example, a user attempting to classify cell types abundant in a tissue sample may select low sensitivity (preferring to exclude substantially all background partitions and discard some cell data to some extent). A user interested in very rare cells (e.g., cancer cells) in a blood cell sample may select very high sensitivity. The computer system can provide guidance to the user when selecting sensitivity. After any cell calling step, the system performs barcode processing.

[0010] The present invention uses multimapping to store information when an array read maps to two or more locations within a reference. The sequencing step of RNA-Seq typically, often in FASTQ format, provides multiple sequence reads. Analysis of those reads typically involves mapping each read to a reference such as an mRNA or gene list specifically used for expression analysis, or even simply the publicly available human genome. Mapping reads to a reference generally involves some version of an alignment algorithm such as pairwise alignment, alignment by the Burrows-Wheeler transform, comparing k-mers of the string to hash k-mers of the reference, or using a suffix or prefix tree constructed from the reference and / or the read, or any combination thereof. Various tools are available in the art for performing alignment according to such methods (see, for example, Dobin, 2013, STAR: ultrafast universal RNA-seq aligner, Bioinformatics 29(1):15-21, which is incorporated by reference). The alignment finds the best matching location within the reference for each read and, implicitly, the gene (represented within the reference) from which the RNA transcript (represented within the read) was produced. There are several reasons, rooted in biology and informatics, as to why a read may map to a reference at two or more locations. For example, the reference may contain homologous genes (e.g., from duplications) or copy number variations, or the read may be short enough to allow for two or more distinct mappings to the reference. Prior art methods have simply discarded such ambiguous mapping reads. In contrast, the system of the present invention processes such results as multimapping of the reads and stores the reads as multimapped reads. Storing information regarding multimapped reads enables analysis of scRNA-Seq sequencing data and discovery of biological information regarding the expression or copy number variation of homologous genes in a sample.

[0011] After barcode processing, cell calling, and read mapping, the system of the present invention provides tools for the analysis of the results of scRNA-Seq data. The analysis according to the present invention generally includes cell clustering, differential expression analysis, and in particular, the extraction of various metrics for providing measures of the expression levels of various genes in each cell. The system of the present invention provides various tools for assisting in display, review, and interpretation. For example, the clustering method can provide a classification of cells. In particular, the system of the present invention can be used to output basic metrics, barcode rank plots, clustering maps, differential expression tables, or combinations thereof. Preferred embodiments of the system operate in a computer system including at least one processor coupled to a memory subsystem. The functions provided by the system can be implemented in a local computer, a server system, or a cloud-based computing system. Software modules for implementing the described functions can be developed in any suitable environment, such as, for example, python, Ruby on rails, c++, others, or combinations thereof. Versions of the system of the present invention are executable on mac osx®, Linux®, and windows® operating systems. Preferred embodiments operate in a server or cloud environment and are operable to receive sequence data, for example, in FASTQ or FASTA format, from a sequencing resource such as a next-generation sequencing (NGS) instrument or a genomics facility. The system can align reads against a reference that generates a sequence alignment map (SAM) or a binary alignment map (BAM) and provide an output to a user computer. Details and variations of embodiments of the functions of the system are described herein.

[0012] In certain aspects, the present invention provides a system for expression analysis using cell calling. The system includes a processor coupled to a memory subsystem that includes instructions executable by the processor, the instructions causing the system to receive sequence data generated by sequencing RNA from a plurality of partitions, for each partition, correlate the count of barcode sequences in the sequence data to the probability that the partition contained one fully isolated cell, receive a user selection for the sensitivity to the probability that the partition contained one fully isolated cell, and analyze the mRNA levels among the partitions that meet the user selection sensitivity.

[0013] The system may be operable to present at least three (e.g., five) distinct computed sensitivity levels to a user (via a graphical user interface displayed on a computing device having an input / output interface) and receive a user selection via the computing device. In some embodiments, the system calculates a function of the barcode count relative to the barcode rank, divides the function by a predetermined value, and selects a cutoff level for the barcode count such that when the cutoff level is exceeded, the partition is considered to contain one fully isolated cell, thereby correlating the barcode sequences to probabilities. The system may be operable to re-analyze the mRNA levels using a new user selection for sensitivity after analyzing the mRNA levels. In certain embodiments, analyzing the mRNA levels includes assigning sequence reads to cells using cell barcodes, deduplicating sequence reads using universal molecular identifiers, mapping the deduplicated reads to a reference, and counting the deduplicated reads that map to genes in the reference as a measure of the expression level of those genes. The system may be operable to downsample the sequence data by mapping fewer reads to the reference than all of the deduplicated reads. In some embodiments, the number of deduplicated reads that are mapped is selected such that the mRNA levels from the sequence data will be normalized to levels calculated from at least one other experiment.

[0014] Aspects of the present invention provide a system for expression analysis using multimapping. Such a system includes a processor coupled to a memory subsystem that includes instructions executable by the processor, the instructions causing the system to receive sequence data (e.g., Illumina output) generated by sequencing RNA from a single cell, map at least one sequence read from the sequence data to a reference that includes reference gene information, identify at least a first location and a second location within the reference to which the sequence read maps with at least a threshold matching score, and store the sequence read in the memory subsystem with markup identifying the read as a mapping to at least the first or second location within the reference. Preferably, the system is operable to map a plurality of reads to the reference and select a first or second location within the reference for a sequence read based on where a greater number of the plurality of reads map. The system may be operable to store the sequence read with markup identifying the read as a mapping to both the first and second locations. In some embodiments, the system assigns a first weight to the mapping to the first location and a second weight to the mapping to the second location. The first and second weights may be at least partially based on respective first and second alignment scores between the sequence read and the reference. The system may be operable to use the markup and the reference gene information to provide a report describing gene duplication or copy number variation in a single cell.

[0015] In some aspects, the present invention provides a system for expression analysis using barcode processing. Such a system includes a processor coupled to a memory subsystem containing instructions executable by the processor, the instructions causing the system to receive sequence data generated by sequencing RNA from a single cell, compare barcodes from sequence reads within the sequence data to a barcode whitelist, and, if a barcode within one sequence read does not match a whitelist barcode at a predetermined Hamming distance value, perform further analysis to omit one sequence read from further analysis in which barcodes from a first tier are compared to a cell barcode whitelist and barcodes from a second tier are compared to a UMI whitelist for each sequence read, deduplicate reads in which the first and second tier barcodes match the whitelist within a predetermined Hamming distance value, and provide a measure of mRNA level from the deduplicated reads. The system may be operable to replace barcodes within sequence reads with index values that occupy less space within the memory subsystem prior to the deduplication step.

[0016] For each of the above aspects, the present invention provides a corresponding method that includes performing the recited functions using the described system. Thus, the present invention provides a system and method for assessing gene expression levels from scRNA-Seq experiments. The system and method use an in-silico cell calling tool to distinguish partitions containing cells from background partitions that did not fully capture a cell. The system and method also provide a barcode processing tool that tests barcodes from at least the cell and UMI tiers against a whitelist and optionally converts them to a compact index format. Further, the system and method may implement a multimapping tool to save information regarding reads that map to multiple locations within a reference. BRIEF DESCRIPTION OF THE DRAWINGS

[0017]

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

Figure 8

Figure 9

Mode for Carrying Out the Invention

[0018] The present disclosure provides a software platform for analyzing single-cell RNA data that can be obtained using particle template instant partitioning (PIP). The system of the present invention proposes a comprehensive analysis solution that provides users with detailed metrics, gene expression profiles, and basic cell quality and clustering indicators. The output of the system can also be used for subsequent specialized tertiary analysis streams.

[0019] System Requirements ● Preferred embodiments of the system can operate on Windows, macOS, and Linux operating systems. ● Preferably, the system includes at least about 64 GB of RAM and 1 TB of hard disk space (to accommodate multiple results). ● In some embodiments, the version of the system operating on Windows uses a local installation of Docker or the like, which may enable the system to seamlessly interoperate with modules that are not themselves compatible with Windows.

[0020] Overview of the Pipeline The system of the present invention may perform any combination of the following steps. 1) FASTQ processing, including barcode matching, QC, and (optionally) downsampling 2) Mapping 3) UMI counting 4) Cell calling 5) Clustering 6) Differential expression 7) Metric extraction 8) Report generation

[0021] Barcode Processing The sequencing process of scRNA-Seq typically generates multiple sequence reads in one or more FASTQ files. One advantage of PIP (over legacy microfluidics) is that multiple cells can be isolated simultaneously within a partition and then queried to generate a single FASTQ file (without the need to concatenate FASTQ files from different runs / days). The process can start by processing the FASTQ file to extract reads associated with barcodes, removing unwanted technical sequences from the data before mapping, and optionally downsampling the data. For example, a module of a system developed in Python, C++, or other such environment can optionally read from the FASTQ file and process each entry to read the barcode sequence and the RNA sequence. Typically, the barcode sequence will be compared to a list of barcodes and the RNA sequence data will be mapped to a reference to identify genes or transcripts.

[0022] In certain embodiments, barcode processing by the systems of the present invention may use a hierarchical barcode system, and optionally, all levels may include one of a specified list of possible barcodes. Levels may refer to levels of information specific to a barcode, such as cell barcodes and then UMIs. There may also be high-level barcodes, such as from different patients or samples, and / or intermediate-level barcodes, such as barcodes specific to a gene (e.g., where sequence-specific primers may be used). Typically, scRNA-Seq includes cell and UMI barcodes, but other levels may be included for any category of information, such as date or location. Including barcodes in reagent oligonucleotides leads to sequence data that includes the sequences from the barcodes. A computer system may match the barcodes against a list or database to assign sequence reads to cells or individual molecules. Barcode matching can be done very efficiently by matching each level in isolation. Some embodiments of the systems of the present invention allow for a Hamming distance of 1. Limiting the Hamming distance to 1 means that the bases in a sequence read, e.g., “read 1” (“R1”) in a FASTQ file corresponding to the position of each level, can differ from the barcode by only one base and still match that barcode.

[0023] To simplify the barcode list sent to a reference mapping module, such as STAR, and reduce its size, hierarchical barcodes may be replaced by a new set of generated barcodes. For example, each barcode may be replaced by a number, such as its numerical position on the list of barcodes. In some embodiments, the new barcodes are generated only for barcodes that have at least one matching read.

[0024] Quality Control In some embodiments, the R2 (Read 2) fastq file contains cDNA constructs and also contains technical sequences related to the PIP chemistry, such as template switch oligo (TSO) and polyA sequences. The occurrence rate of these sequences may be inversely correlated with the fragment size. Since reads containing extra sequences are less likely to be mapped, it may be preferable to remove the TSO sequence from the 5' end and the polyA sequence from the 3' end in the R2 fastq file. Additionally, depending on the type of chemistry, a fixed number of non-cDNA bases may be removed from the 5' end. Reads that are less than 20 bases long after trimming may be discarded from further analysis.

[0025] After all reads have been matched to barcodes, inevitably, some barcodes will be associated with a small number of reads. Since those barcodes are very unlikely to be classified as cells later, it may be economical to exclude them and the associated reads from further analysis. The user can define the minimum number of reads per barcode. If such a threshold is specified, reads associated with low-count barcodes will not be passed to STAR for mapping. See "Cell calling", which is discussed in more detail elsewhere in this specification.

[0026] Downsampling In some experiments, it may be desirable to standardize the sequencing depth so that experimental conditions can be reliably compared. Accordingly, the system of the present invention is operable to downsample data to a specified number of reads. The reads can be randomly selected, and the probability of selection is determined by the total number of reads (either provided by the user or determined by counting the lines in the fastq input).

[0027] Output In certain embodiments, the fastq processing step results in the following output. - <output-root> / metrics / barcode_stats.json: General statistics regarding the quality of barcode matching. - <output-root> / metrics / barcodes / Generated_barcode_read_info_table: List of all generated barcodes passed to STAR for mapping and counting. The list includes the hierarchy associated with the original barcode and the number of exact and error-corrected matches for each barcode. - <output-root> / metrics / barcodes / barcode_whitelist: List of all generated barcodes provided to STAR as "whitelist" for mapping and counting. - <output-root> / barcodes_fastqs / *: Intermediate fastq files used by STAR for mapping. Unless otherwise indicated, these files are deleted at the end of the analysis.

[0028] Mapping The sequencing step of RNA-Seq typically, and often, provides multiple sequence reads in FASTQ format. The analysis of those reads typically involves mapping each read to a reference such as an mRNA or gene list specifically used for expression analysis, or even simply the publicly available human genome. Mapping the reads to a reference generally involves some version of an alignment algorithm such as pairwise alignment, Burrows-Wheeler transform alignment, comparing the k-mers of the string to the hash k-mers of the reference, or using a suffix or prefix tree constructed from the reference and / or the reads. Various tools for performing alignment according to such methods are available in the art. A particular preferred embodiment uses a read mapping module called STAR. See Dobin, 2013, STAR: ultrafast universal RNA-seq aligner, Bioinformatics 29(1):15-21, which is incorporated by reference. In particular, the system of the present invention uses STAR (Spliced Transcripts Alignment to a Reference) and related STARsolo to map reads to the transcriptome and quantify transcript abundance. STAR uses a special index constructed from the original transcriptome of interest. Some common indexes such as for human and mouse can be downloaded from sources such as Fluent BioSciences. However, for some unique situations, the system of the present invention allows the user to create a unique or experiment-specific index.

[0029] The documentation of the STAR product is available from the issuer. In a preferred embodiment, the system of the present invention is, for example, as described in the STAR documentation, <output-root>Store the output of read mapping (e.g., by STAR) in a directory such as / starsolo.

[0030] Note that if STAR is not yet supported on Windows and Mac, the system will use a Docker image to run STAR. When run on those operating systems (or explicitly requested on Linux), the system of the present invention will launch STAR from the Docker image.

[0031] Read mapping generally refers to detecting the location within a reference where an array read maps best. Those systems include programming logic for using information that may potentially be revealed when an array read maps well to two or more locations on the reference.

[0032] The present invention uses multi-mapping to store information when an array read maps to two or more locations within a reference. There are several reasons rooted in biology and informatics as to why a read may map to a reference at two or more locations. For example, the reference may contain homologous genes (e.g., from duplications) or copy number variations, or the read may be short enough to allow for two or more distinct mappings to the reference. Prior art methods have simply discarded such ambiguous mapping reads. In contrast, the system of the present invention processes such results as multi-mapping of reads and stores the reads as multi-mapped reads. Storing information about multi-mapped reads enables the analysis of scRNA-Seq sequencing data and the discovery of biological information regarding the expression or copy number variation of homologous genes in a sample. Importantly, reads that multi-map are preferably not discarded (as compared to prior art systems that discard such reads). In some embodiments, reads are assigned to one of the locations according to a selection logic such as randomly, or (when there are two locations) by alternating the assignment between the two locations when multiple reads map to two, or by mapping each read to the location where most of the reads map, or by giving "partial credit" (e.g., proportional to the number of different locations or weighted by an alignment score), or other, or combinations thereof.

[0033] UMI count After STAR maps the reads in the fastq file to a reference, it creates a Binary Alignment Map (BAM), which is the binary version of the Sequence Alignment Map or SAM file. The system then parses the BAM file to create a dataset with all the molecules (barcode + UMI combinations) in the sample mapped to the reference. The molecular information is saved in an HDF5 format file with a suitable name, e.g., molecule_info.h5. The file can be parsed by the system to create a raw count (feature barcode) matrix. In some embodiments, the matrix is included in the files: matrix.mtx.gz, barcodes.tsv.gz, and features.tsv.gz, and is included, for example, in a directory such as raw_matrix. These files form a sparse matrix that contains the full count of unique UMIs for each barcode and gene. The format of the matrix is preferably consistent with most downstream analysis tools, including third-party analysis tools.

[0034] Cell calling After obtaining the UMI count matrix, the barcodes are separated into cell-containing PIPs and background PIPs, i.e., PIPs that did not fully capture the cells. This separation is based on a barcode rank plot, which orders the barcodes based on the number of unique UMIs associated with each barcode. Compared to microfluidic platforms where a limited number of cells are successively isolated into partitions, the use of PIPs (where multiple cells are simultaneously isolated into aqueous partitions) reveals several issues associated with the throughput and amount of data initially achieved by the PIPs. For example, using PIPs reveals that it is beneficial to have an automated and accurate tool, also known as cell calling, to correctly call whether a partition contains cells as opposed to being a background partition that does not fully capture the cells. The present disclosure recognizes that, as a practical matter, it is beneficial for the partition containing cells to be automatically, and even "on the fly", distinguished from background partitions (that did not fully capture the cells) in terms of the amount of data and throughput associated with the PIPs. For example, from the sequence data, all of the identified barcodes can be ranked in order by the number of UMIs associated with each barcode. Barcodes associated with a larger number of UMIs are more likely to originate from cells, while barcodes associated with a smaller number of UMIs are more likely to originate from background partitions.

[0035] Figure 1 shows a barcode rank plot. Typically, cell barcodes are concentrated in the top mode of the rank plot, while background barcodes are concentrated in the low UMI count mode. Thus, the objective of cell calling is to find the points within the first "knee" area that separates the two modes. One approach to finding that point involves dividing the number of UMIs at the top of the rank plot by a constant denominator (e.g., 10). The method for finding the inflection point can be very sensitive to the shape of the rank plot, and the shape of the rank plot varies greatly between different sample types. The system of the present invention can automatically find the points on the barcode rank plot that divide the cell partitions from the background partitions, or the system can present the user with one or more options that assist in enabling the user to set or select that point, or the system can perform a combination of automatic selection or proposal in combination with user selection.

[0036] The system of the present invention can present the user with a sensitivity / specificity selection (e.g., select one of five pre-set options). Here, the plot is provided to show the results of calling cell barcodes up to one-tenth of the top-ranked barcodes for HEK / 3T3, PBMC, and breast tissue samples.

[0037] Figure 2 shows a barcode rank plot for HEK / 3T3 cells.

[0038] Figure 3 shows a barcode rank plot for peripheral blood mononuclear cells (PBMC).

[0039] Figure 4 shows a barcode rank plot for breast tissue samples.

[0040] Each cell part is emphasized in bold. Cell calls that require HEK / 3T3 appear to match the visual inspection of the barcode rank plot: the entire high UMI mode is captured. However, in the case of PBMCs and breast cells, the top mode has a more negative slope and is less clearly defined, so the number of cell barcodes is clearly underestimated. This could be a deliberately made choice by the user. For example, a user studying PBMCs may need to optimize the assay for data volume (largely ignoring or excluding background PIP) for the quality of the information required due to underlying biological reasons compared to HEK / 3T3.

[0041] In some embodiments, to adapt to differences between samples and facilitate accurate cell calling, the system of the present invention provides cell calling at multiple (e.g., five) different sensitivity levels using the following steps. 1) Select a starting point S near the top of the rank plot. To avoid selecting outliers, the system uses the point at which the slope of the plot first drops below 10% as the starting point. 2) Let the number of UMIs of the barcode at the starting point be M. For each cell calling sensitivity level L, find the number of UMIs U as follows.

[0042]

Equation

[0043] The system can use any suitable approach to determine spots that include knees or inflection points. For example, behind the scenes (or, for example, on a screen graphically presented to the user), the computer system can calculate the number of UMIs (e.g., on the y-axis) across the rank (x-axis) of the ranked barcodes, fit a function to the curve, find the inflection point by taking the derivative, and input options for selected points along the curve into the user interface. The user can select a sensitivity level for cell calling, and the analysis proceeds at the selected sensitivity. For example, a user attempting to classify cell types that are abundant in a tissue sample may select a low sensitivity (preferring to exclude substantially all background partitions and discard some cell data to some extent). A user interested in very rare cells (e.g., cancer cells) in a blood cell sample may select a very high sensitivity. The computer system can provide guidance to the user when selecting sensitivity.

[0044] Figure 5 shows cell calls at five sensitivity levels for a PBMC sample.

[0045] Cell calling inherently involves a trade-off between the number of cells and their quality. The number of cells can be arbitrarily increased by calling the rank plot deeper, but this increases the likelihood of selecting cells with low RNA content or otherwise damaged cells. The optimal number of cells to recover depends on the sample type and the purpose of each experiment. For example, if some cell types in the sample (such as neutrophils or red blood cells) have low RNA expression, they may be present in large numbers in the lower UMI regions of the plot, or the UMI count distribution may even be unimodal. In that case, high-sensitivity cell calling may be guaranteed. On the other hand, when high clustering accuracy is required for accurate cell type characterization, it is preferable to use low-sensitivity cell calling. Cell calls at different sensitivity levels allow the user to observe a wide range of possible results and determine the correct sensitivity for their data.

[0046] Manual selection of the number of cells In some cases, the user may desire more precise control of the number of cells called than that provided at the five sensitivity levels. In some embodiments, the system enables the user to specify the desired number of cells, in which case the exact number of barcodes is selected from the top of the rank plot. The outputs described herein will be generated for the selected set of barcodes. For some embodiments of the system of the present invention, this may be run from the command line, for example, using the prefix force_N, where N is the number of cells input by the user.

[0047] Output After the cells are called, the directory can be created as follows for each sensitivity level within the parent directory. A text file listing the indices (starting from 1) of the barcodes selected as cells, and an image with a barcode rank plot <output-root> / cell_calling. Additionally, filtered count matrices can be generated for each sensitivity level and used for downstream analysis. The filtered matrix is <output-root>It is stored in / cr_filtered_quant. A person skilled in the art will understand how to describe the functions for storing these outputs in a suitable programming or development environment such as Python.

[0048] Analysis In addition to barcode processing, cell calling, and read mapping, the system of the present invention provides tools for the analysis of the results of scRNA-Seq data. The analysis according to the present invention generally includes cell clustering, differential expression analysis, and in particular, the extraction of various metrics for providing a measure of the expression levels of various genes in each cell. The system of the present invention provides various tools for assisting in display, review, and interpretation. For example, the clustering method may provide a useful tool for cell classification.

[0049] Clustering After selecting the cell barcodes, it may be desirable to perform clustering analysis to quantify and visualize the heterogeneity within the cell population and identify different cell types. For all sets of called cells (e.g., different sensitivities or numbers forced by the user), the system of the present invention can create a clustering map using approaches such as K-means clustering with different values of K, as well as Leiden (graph-based) clustering that automatically determines the number of clusters based on nearest neighbors.

[0050] Preprocessing step The starting point for clustering is preferably the sparse row-column representation of the filtered count matrix, where each row represents a cell and each column represents a gene. Several processing steps can be performed before the actual clustering algorithm is executed. 1) Cell normalization: Normalize the data of each cell to unit norm. This is performed so that cells having similar expression profiles but different relative RNA abundances can still be clustered together. 2) Logarithmic transformation: The matrix is transformed using the ln(x + 1) function to bring the matrix closer to a normal distribution. 3) High-variance gene selection: To maximize the efficiency of clustering, a subset of genes with high variance is selected. For all genes, the variance of expression across cells is calculated. The genes are then ranked by variance, and the top N percentile of genes are selected, where N is user input. 4) Scaling: As is common in clustering analysis, it may be preferable to scale the feature columns (genes) so that each column has a mean of 0 and a standard deviation of 1. 5) PCA: The system of the present invention can use Principal Component Analysis (PCA) to reduce the dimensionality of the data from the scaled set of high-variance genes to a small user-controlled number of components that best explain the variance in the data.

[0051] Following the preprocessing steps, clustering can be performed in the new PCA space.

[0052] In some embodiments, the system performs clustering using the Leiden algorithm, which is a standard in the art for graph-based clustering in scRNA-seq datasets. A k-nearest-neighbor (KNN) graph is constructed from the PCA matrix, converted to a directed node-edge graph, and passed to Leiden. The Leiden algorithm uses modularity optimization of the graph community to determine the ideal cluster membership. Leiden has the additional advantage of minimal user input (e.g., number of clusters, cut-off, etc.) compared to other methods. See, for example, Franzen, 2020, Alona: a web server for single-cell RNA-seq analysis, Bioinformatics 36(12):3910-3912; Wu, 2020, Tools for the analysis of high-dimensional single-cell RNA sequencing data, Nat Rev Neph 16:408-421, and Zhu, 2020, Single-cell clustering based on shared nearest neighbor and graph partitioning, Interdiscip Sci 12(2):117-130, all of which are incorporated by reference.

[0053] In certain embodiments, the system also uses a common K-means algorithm to divide the cells into discrete clusters. This involves running the algorithm for different values of K, which represent a predetermined number of clusters. The minimum and maximum of the range of K values are controlled by the user and should be adjusted according to the number of cell types expected to be present in the sample.

[0054] Differential expression Once clustering is complete, it may be important to determine the primary genes associated with each cluster. For this step, the system of the present invention may use a Wilcoxon rank sum test (Mann-Whitney U) for independent variables. This non-parametric test is used to determine whether two samples are from the same distribution by ranking the values associated with two groups together and calculating the sum of those ranks within each group. Those ranks are used to calculate a test statistic (U), based on which the clustered cell population and the non-clustered cell population can be compared. For each cluster, the system calculates the differential expression for the genes associated with the cells of that cluster across all other cells. The top genes sorted by descending Z-score for each cluster are exported to a csv file. This process is performed for different values of K in K-means clustering, as well as for graph-based Leiden clustering.

[0055] Output In particular, the system of the present invention may be used to output basic metrics, barcode rank plots, clustering maps, differential expression tables, or combinations thereof. The directory may include, for example, the following for each or any set of the cells called: <output-root>It can be created in a directory such as / clustering. ● A csv file listing the graph-based clustering and the cluster assignment of each cell from all K-means runs. ● A csv file containing the silhouette scores for the graph-based clustering and all K-means runs. The silhouette score is a measure of clustering quality that takes into account both the within-cluster distance and the between-cluster distance. Note that this is calculated in the PCA space used for clustering and not in the visualized UMAP space. ● UMAP plots for the graph-based clustering and all K-means runs. UMAP (Uniform Manifold Approximation and Projection) is a common method for projecting multi-dimensional data into a 2D space. ● A csv file for the graph-based clustering and all K-means runs containing the top-ranked genes for each cluster. ● A csv file similar to the top gene file but containing additional statistics, Z-scores, p-values, adjusted p-values, and log2 fold changes for each gene.

[0056] The following plots show the clustering results for the HEK / 3T3 and PBMC samples. The number of clusters shown in each case is the one that produced the highest silhouette score.

[0057] Figure 6 shows the clustering results for the HEK / 3T3 sample.

[0058] Figure 7 shows the clustering results for the PBMC sample.

[0059] Metric Extraction In the final step of the analysis, various metrics are calculated and reported for each sensitivity level (or for a single set if the number of cells is specified by the user).

[0060] In some embodiments, the system of the present invention may include one or any combination of several basic metrics. For example, one or more of the following metrics may be included in all analyses. ● Total number of reads in the input fastq file ● Percentage of mapped reads (note that the denominator does not include reads excluded from low-count barcodes and reads that are too short after trimming). ● Number of called cells ● Average reads per cell (this is the total number of input reads divided by the number of cells). ● Percentage of reads in cells (this is the number of reads associated with cell-containing barcodes divided by the total number of input reads). ● UMI duplication rate in cells (this is the number of mapped reads in a cell divided by the number of unique reads (duplicate-excluded reads from the same transcript) in the cell. For example, if all UMIs had 3 reads containing it, i.e., 2 duplicate-excluded reads per UMI, the UMI duplication rate would be 3). ● Total number of UMIs in cells ● Median UMI per cell (this is based on the distribution of counts of unique UMIs (transcripts) associated with each cell barcode). ● Number of genes expressed in cells (this number includes all genes with non-zero counts in the filtered matrix, giving equal weight to all genes regardless of expression level). ● Median gene per cell (this number includes all genes with non-zero counts in the filtered matrix, giving equal weight to all genes regardless of expression level). ● Percentage of cells where more than 50% of the transcripts are from mitochondrial genes ● Percentage of cells where less than 5% of the transcripts are from mitochondrial genes

[0061] Estimation of Barnard metrics and multiplexing rate Embodiments of the system provide tools useful for assessing the quality of a particular chemistry or reaction setup in an scRNA-Seq protocol such as Drop-seq or inDrop, or one of those protocols. One way to assess the quality of single cell RNA technology is typically a complex species (a "barnyard") experiment involving a mixture of human and mouse cells. Such experiments can provide information about the incidence of multiplets, or beads that capture more than a single cell. It can also assess the degree of interference from background RNA by measuring cross-species contamination. More human genes are found in mouse cells and vice versa, and the higher the likelihood that more background contamination is present in the system.

[0062] When the reference transcriptome is a composite transcriptome of human and mouse, the system of the present invention is operable to provide an additional set of metrics specific to barnyard experiments, as follows. ● The number of human cells (this is the number of cell barcodes from which at least 85% of the transcripts are derived from the human reference). ● The number of mouse cells (this is the number of cell barcodes from which at least 85% of the transcripts are derived from the mouse reference). ● The number of multiplets (this is the number of cell barcodes from which less than 85% of the transcripts are derived from the human or mouse reference, and since a multiplet can consist of two human cells or two mouse cells, or (rarely) contain more than two cells, it should be noted that this probably underestimates the true number of multiplets). ● The total number of UMIs in human cells ● The median human UMI per human cell ● The number of human genes expressed in human cells ● The median human gene per human cell ● The percentage of mouse UMIs in human cells ● The total number of UMIs in mouse cells ● The median mouse UMI per mouse cell ● The number of mouse genes expressed in mouse cells ● Median mouse gene per mouse cell ● Percentage of human UMIs in mouse cells

[0063] An important output of the Barnyard analysis is called the Barnyard plot.

[0064] Figure 8 is a Barnyard plot showing the number of mouse transcripts versus the number of human transcripts, with different shades indicating human cells, mouse cells, and multiplets.

[0065] Output After all metrics are calculated, the results for each set of cells called can be reported in files such as the following. ● Basic metrics are in csv and json formats <output-root>It is stored in / metrics. ● When applicable, the burn yard metrics are in csv and json formats <output-root> / Stored in the barnyard.

[0066] Analysis report The operating system of the present invention for analyzing sequence data from scRNA-Seq experiments optionally generated a report such as a document in a format such as XML or HTML. Such reports are <output-root>It can be stored in a suitable directory such as / report and may include basic metrics, barcode rank plots, clustering maps, and differential expression tables. The user can toggle between different clustering modes, between graph-based and K-means at different K values, which will change the content of the maps and tables, and when the "bone yard" mode is used, their statistics and bone yard plots will also be included accordingly.

[0067] In a standard analysis, a report can be generated for each sensitivity level. An additional composite report (e.g., "combined.html") can include all sensitivity levels. The user can be given the option to browse through different sensitivities in the composite report and select the optimal sensitivity level or determine the number of cells to force subsequent reanalysis.

[0068] In the case of an analysis using a manually selected (forced) number of cells, a report will be generated for this specific number and a new tab can be added to the composite report.

[0069] Execution of the pipeline Preferred embodiments of the system operate in a computer system including at least one processor coupled to a memory subsystem. The functions provided by the system can be implemented in a local computer, a server system, or a cloud-based computing system. Software modules for implementing the described functions can be developed in any suitable environment, such as, for example, python, Ruby on rails, C++, others, or combinations thereof. Versions of the system of the present invention are executable on mac osx, Linux, and windows operating systems. Preferred embodiments operate in a server or cloud environment and are operable to receive sequence data, for example, in FASTQ or FASTA format, from a sequencing resource such as a next-generation sequencing (NGS) instrument or a genomics facility. The system can align reads against a reference to generate a sequence alignment map (SAM) or a binary alignment map (BAM) and provide an output to a user computer. Details and variations of embodiments of the functions of the system are described herein.

[0070] Figure 9 shows a system 901 including at least one local computer 905. System 901 can also include a server or cloud computer 909. Preferably, any local computer 905 or server or cloud computer 909 included in system 901 includes at least one processor coupled to a memory subsystem. Preferably, any local computer 905 or server or cloud computer 909 included in system 901 is operable to receive sequence data from a sequencing device 921 or a genomics facility via a communication network 917, whereby the machines of system 901 communicate and interoperate with each other. The memory subsystem within local computer 905 or server or cloud computer 909 of system 901 preferably includes program instructions executable by a processor to cause system 901 to perform the analysis functions described herein.

[0071] Thus, the system of the present invention provides a complete analysis solution for a single sample. The system of the present invention can take in a pair-end set of fastq files and return a set of outputs including a full barcode x gene count matrix, a set of filtered count matrices for different cell calling sensitivity levels, and various sample metrics. For mixed human / mouse samples ("yardburner" experiments), the results may include specific metrics related to species separation. The composite results are summarized in an HTML report. In certain embodiments, the system of the present invention is executed from a command line prompt on Linux, Windows, and Mac. Note that file names and directory names containing spaces may be enclosed in quotes.

[0072] Please refer to the following description of all required and optional arguments for a particular embodiment of the system of the present disclosure. To see the full list of arguments from the command line, execute the following. $ pipseeker count --help.

[0073] Required arguments $ pipseeker count --input-path <path to fastq files> --output-root <destination for all outputs> --STAR-index-path <path to mapping reference> --input-path (string)

[0074] Path to the source fastq files. If this is a directory name, all files ending with.fastq or fastq.gz will be included as input. Alternatively, if there are files from multiple samples within a single directory, a prefix can be specified at the end of the path, which will limit the files included to only those whose names begin with that prefix.

[0075] For example, consider a situation where the input directory contains the results of 2 lanes from 2 different samples. $ ls my_output_dir sample_1_L001_R1.fastq.gz sample_1_L001_R2.fastq.gz sample_1_L002_R1.fastq.gz sample_1_L002_R2.fastq.gz sample_2_L001_R1.fastq.gz sample_2_L001_R2.fastq.gz sample_2_L002_R1.fastq.gz sample_2_L002_R2.fastq.gz

[0076] Using --input-path my_output_dir would result in an undesirable situation where all files are included in the analysis as components of the same sample. Instead, in two separate runs, use --input-path my_output_dir / sample_1 and --input-path my_output_dir / sample_2 to include only the relevant files for each sample. --output-root(string)

[0077] The directory where all PIPseeker outputs will be stored. It may be the same as the input path or a different directory. --star-index-path(string)

[0078] Path to the directory containing the indexed reference to be used for mapping.

[0079] Optional arguments Typically: --chemistry(string,default:mark6)

[0080] Enable the user to input the chemical substances used to generate the sample. The accepted values are specific to the embodiments (e.g., mark3 and mark6). --downsample-to(integer))

[0081] An integer representing the target number of reads to be mapped when downsampling is desired. If specified, the reads will be randomly selected and only the selected subset will be supplied to STAR for mapping. --input-reads(integer)

[0082] An integer representing the number of input reads when --downsample-to is specified. The total number of reads must be known in advance to calculate the downsampling ratio.

[0083] --input-reads should include the total number of reads combined across all input files. If not specified, the reads will be counted manually, so using this argument can shorten the execution time. --verbosity(integer,default:1)

[0084] Verbosity level: 0 (silent), 1 (concise), or 2 (verbose).

[0085] Barcoding --min-reads-per-barcode(integer,default:1) The minimum number of reads associated with a barcode. If a barcode is associated with fewer reads than the specified number, all reads from this barcode will be discarded from further analysis. Increasing this number can shorten the mapping time by not processing barcodes that are less likely to represent cells. Note that the low-count part of the barcode rank plot will be incomplete and the mapping statistics may change slightly.

[0086] --retain-barcoded-fastqs Flag for retaining intermediate barcoded fastqs. In most cases, it is not necessary. If not specified, those fastqs will be discarded after analysis is complete. Note that the input fastqs are always retained regardless of this flag.

[0087] Mapping --star-threads(integer,default:8) Number of CPU threads to use for STAR mapping (default: 8).

[0088] Cell calling --force-cells(integer) Force a specific number of barcodes to be considered cells. If specified, this number of barcodes will be selected from the top of the barcode rank plot. This will replace the default sensitivity-based cell calling mode.

[0089] Clustering --clustering-percent-genes(integer,default:10) Percentage of genes to use for clustering. Genes with the highest variability in expression up to this percentile rank will be included.

[0090] --principal-components(integer,default:15)

[0091] Number of components for PCA dimensionality reduction. --percent-genes(float,default:10)

[0092] Percentage of genes to select for clustering based on the highest variance in expression. --min-clusters-kmeans(integer, default: 2)

[0093] The minimum number of clusters for K-means clustering. The value of K will start from this number and be incremented. --max-clusters-kmeans(integer, default: 12)

[0094] The maximum number of clusters for K-means clustering. The value of K will be incremented until it reaches this number. --nearest-neighbors(integer, default: 10)

[0095] The number of nearest neighbors for graph-based clustering. This defines the minimum number of cells per cluster. --diff-exp-genes(integer, default: 50)

[0096] The number of top differentially expressed genes included in each cluster.

[0097] Metrics --run-barnyard A flag to enable the generation of the "Barnyard" metric and include it in the analysis report.

[0098] Report --id The sample ID to include in the output report.

[0099] --description The description of the sample to include in the output report.

[0100] Final stage re-run: Re-analysis After observing the results of a complete analysis, the user may wish to change parameters related to later stages, such as the range of K values for K-means clustering. The user may also wish to force a particular number of barcodes to represent cells. This should not typically require repeating the most time-consuming barcode generation and mapping steps.

[0101] When using the reanalysis feature, in a preferred embodiment, barcode generation and mapping are skipped and analysis is started directly upon cell calling. If this mode is used, the input path should point to a directory containing the complete analysis results from a previous execution. Arguments related to --output-root, --chemistry, barcode generation, and mapping will be ignored.

[0102] For example, assume that after running a complete pipeline and visually inspecting the results, none of the standard cell calling sensitivity levels appear to capture the correct region of the barcode rank plot. The user can force a particular number of cells without re-running the entire pipeline.

[0103] $ pipseeker reanalyze --input-path my_analysis_result_dir --force-cells 5000 Using the commands, hardware, and functions described, the present invention provides a system comprising a processor coupled to a memory subsystem including instructions executable by the processor to cause the system to evaluate gene expression levels from scRNA-Seq experiments on a processor. The system uses an in-silico cell calling tool to distinguish between partitions containing cells and background partitions that did not fully capture cells. The system also provides a barcode processing tool that tests barcodes from at least the cell and UMI tiers against a whitelist and optionally converts them to a compact index format. Further, the system may implement a multimapping tool to store information regarding reads that map to multiple locations within a reference. The system also includes analysis tools for providing reports including, for example, clustering, cell classification, and measurement of expression levels. The functionality of the system addresses problems revealed by scRNA-Seq using pre-templated instant partitions (PIPs) in which a very large number of cells are simultaneously isolated into droplets with a speed not achievable by legacy microfluidics.

Explanation of Signs

[0104] 901 System 905 Local computer 909 Server or cloud computer 917 Communication network 921 Sequencing device

Claims

1. A system for expression analysis using cell calling, the system comprising a processor coupled to a memory subsystem, the memory subsystem comprising instructions executable by the processor, the instructions causing the system to receive array data generated by sequencing RNA from a plurality of partitions, for each partition, correlate the count of barcode arrays in the array data with the probability that the partition contained one fully isolated cell, receive a user selection for the sensitivity to the probability that the partition contained one fully isolated cell, analyze the mRNA levels between the partitions that meet the user selection sensitivity.

2. The system of claim 1, further operable to present at least three distinct calculated sensitivity levels to a user via a graphical user interface displayed on a computing device having an input / output interface and to receive a user selection via the computing device.

3. The system of claim 1, wherein the system correlates the barcode array with the probability by calculating a function of the barcode count relative to the barcode rank, dividing the function by a predetermined value, and selecting a cutoff level of the barcode count such that if the cutoff level is exceeded, the partition is considered to contain one fully isolated cell.

4. The system of claim 1, further operable to re-analyze the mRNA levels using a new user selection for the sensitivity after analyzing the mRNA levels.

5. Analyzing the mRNA levels includes assigning sequence reads to cells using cell barcodes, deduplicating sequence reads using universal molecular identifiers, mapping deduplicated reads to a reference, and counting deduplicated reads that map to genes in the reference as a measure of the expression level of those genes.

6. The system of claim 5, further operable to downsample the array data by mapping fewer reads than all of the deduplicated reads to the reference.

7. The system according to claim 6, wherein the number of the deduplicated reads to be mapped is selected such that the mRNA level from the array data is normalized to a level calculated from at least one other experiment.

8. A system for expression analysis using multi-mapping, the system comprising a processor coupled to a memory subsystem, the memory subsystem comprising instructions executable by the processor, the instructions causing the system to receive array data generated by sequencing RNA from a single cell, map at least one array read from the array data to a reference including reference gene information, identify at least a first location and a second location within the reference where the array read maps with at least a threshold matching score, store the array read in the memory subsystem, with a markup identifying the read as a mapping to at least the first location or the second location within the reference.

9. The system according to claim 8, further operable to map a plurality of reads to the reference and select the first location or the second location within the reference for the array read based on the location where a greater number of the plurality of reads map.

10. The system according to claim 8, further operable to store the array read with a markup identifying the read as a mapping to both the first location and the second location.

11. The system according to claim 10, wherein the system assigns a first weight to the mapping to the first location and a second weight to the mapping to the second location.

12. The system according to claim 11, wherein the first weight and the second weight are each at least partially based on a respective first alignment score and a second alignment score between the array read and the reference.

13. The system according to claim 8, further operable to provide a report describing gene duplication or copy number variation in the single cell using the markup and the reference gene information. [[ID=|18]]

14. A system for expression analysis using barcode processing, the system comprising a processor coupled to a memory subsystem, the memory subsystem comprising instructions executable by the processor, the instructions causing the system to receive sequence data generated by sequencing RNA from a single cell, compare barcodes from sequence reads within the sequence data to a barcode whitelist, when a barcode within one sequence read does not match a whitelist barcode at a predetermined Hamming distance value, omit the one sequence read from further analysis, the further analysis being such that for each sequence read, barcodes from a first hierarchy are compared to a cell barcode whitelist and barcodes from a second hierarchy are compared to a UMI whitelist, deduplicate reads where the first hierarchy barcode and the second hierarchy barcode match the whitelist within the predetermined Hamming distance value, provide a measure of mRNA level from the deduplicated reads. A system. **Claim 15** The system of claim 14, further operable to replace the barcodes within the sequence reads with index values that occupy less space within the memory subsystem, prior to the deduplication step.

Images