SCRNA SEQ analysis system

Abstract

The invention provides a system for evaluating the expression level of a gene from an scRNA-Seq experiment. The system uses a computer-simulated cell detection tool to distinguish between partitions containing cells and background partitions in which cells are not fully captured. The system also provides a barcode processing tool that tests barcodes from at least the cellular level and the UMI level for a whitelist and optionally converts these barcodes into a compact index format. In addition, the system may implement a multiple mapping tool to preserve information of reads mapped to multiple locations in the reference.

Description

Technical Field

[0001] The present disclosure relates to expression analysis using methods such as single cell RNA sequencing. Background Art

[0002] Single-cell RNA sequencing (scRNA-seq) refers to a variety of protocols involving sequencing RNA from cells, and in most embodiments, provides a measure of gene expression levels from the sequence data. Some methods of scRNA-Seq rely on isolating cells into droplets, with the potential to assay a large number of cells per experiment. Popular droplet-based protocols include Drop-seq (described in Macosko, 2015, "Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets", Cell, Vol. 161, No. 5: pp. 1202-1214, incorporated by reference) and inDrop (see Klein, 2015, "Droplet barcoding for single-cell transcriptomics applied to embryonic stem cells", Cell, Vol. 161, No. 5: pp. 1187-1201, incorporated by reference).

[0003] Those droplet-based methods of scRNA-Seq typically use two oligonucleotide barcodes added to an RNA molecule or its cDNA copy. The sequences of these two barcodes appear in the sequence data obtained from the cell RNA. Typically, one of these barcodes is the same throughout the droplet, but different between droplets, and acts as a cell barcode, providing a common cell identification tag for all molecules from one cell. The use of cell barcodes allows nucleic acids from multiple cells to be pooled during sequencing, and then the sequence data is separated by computer simulation of cells. Another barcode commonly used in scRNA-Seq is called a unique molecular identifier (UMI). UMI tags each individual molecule with a unique sequence of nucleic acid bases. When those molecules are amplified, for example, by PCR, the 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 considered as duplicates. In a process sometimes referred to as deduplication reads, all identical reads are considered as one read, 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 for each gene in a cell provides a measure of the gene expression level of that cell.

[0004] Many factors make it difficult to provide a measure of gene expression levels from sequence data. For example, barcodes are prone to errors during sequencing and amplification. In addition, much sequence data comes from amplification of stray nucleic acids in droplets that contain no cells or damaged cells or other debris. Summary of the invention

[0005] The present invention provides a system for evaluating gene expression levels from scRNA-Seq experiments. The system uses computer simulation tools for "cell detection", that is, distinguishing partitions containing cells (e.g., droplets) from background partitions that do not fully capture cells. In addition, the system provides tools for "barcode processing", in which barcodes are optionally converted to a compact format and / or tested against a whitelist in a way that greatly increases processing speed while retaining information during deduplication. In addition, the system can implement tools for "multiple mappings", which retain information about reads mapped to multiple locations in a reference (these reads were previously usually just discarded). In addition, the system includes analytical tools for analysis, such as clustering, display, classification, and metric extraction, such as for displaying which cell clusters are present in a sample based on gene expression patterns.

[0006] Those tools for cell detection, barcoding, multiple mapping and analysis are particularly suitable for scRNA-Seq using particle-templated instant partitions (PIPs), which essentially separate very large numbers of cells into emulsion droplets at the same time. To perform scRNA-Seq using PIPs, cells and particles are mixed together in an aqueous mixture. The particles can be hydrogel beads, and those hydrogel beads can carry reagents such as enzymes and / or oligonucleotides. The reaction mixture can be provided with any of a variety of reagents, either in the aqueous phase or on or within the hydrogel bead particles. In addition to any enzymes or oligonucleotides, those reagents may include dNTPs, cofactors, dyes, salts, chemicals, other reagents, or combinations thereof. The aqueous mixture will typically have a number of hydrogel bead particles that corresponds to the number of aqueous partitions (also called droplets or "PIPs") to be prepared. The cells can be introduced at a concentration (i.e., dilution) such that the PIPs are made in which the average number of cells per PIP is no greater than 1. After the aqueous mixture is prepared, the oil can be added to the aqueous phase (optionally, a surfactant is also provided to the mixture). The oil / aqueous mixture can then be sheared, for example, by vortexing a test tube (e.g., a 0.5 mL microcentrifuge tube or larger conical tube). The shear energy causes each hydrogel bead particle to act as a template, almost instantly partitioning the aqueous liquid surrounded by the oil to form multiple stable monodisperse emulsion droplets (also called partitions), hence the term particle template instant partition (PIP). Arbitrarily large numbers of PIPs can be formed simultaneously in one reaction volume, and they form almost instantly during vortexing.

[0007] With very few exceptions, each PIP includes one particle, at most one cell, a cell lysis reagent (e.g., a protease), an oligonucleotide for RNA hybrid capture, such as a barcoded poly-T capture oligonucleotide packaged within the particle or covalently attached to the particle. The PIP may also include other reagents, such as reverse transcriptase, dNTPs, other capture oligonucleotides attached to the particle or free in solution (e.g., template switching oligonucleotides), primers, cofactors, etc. When the partitions or PIPs are formed immediately, each partition will contain an average of 1 or 0 cells. The test tube can be transferred to a thermal cycler or other temperature control device. Using the provided reagents, the cells can be lysed to release RNA, and the RNA can be captured onto the oligonucleotides (usually by hybrid capture, and optionally by enzymatic ligation, such as by a ligase or transposase). Typically, all oligonucleotides of a hydrogel bead particle share a copy of a common "cell barcode" that is (substantially) unique to the particle. In addition, each oligonucleotide may include a UMI that is unique (or at least nearly unique or substantially or effectively unique to the oligonucleotide) at least among the oligonucleotides of the particle. After the RNA molecules are captured onto the oligonucleotides, a typical RNA-Seq will proceed through cDNA synthesis, amplification, ligation of sequencing adapters, and sequencing to generate sequence data that includes nucleotide sequences from the RNA molecules and from the oligonucleotides (e.g., from the barcodes). Systems for barcoding, cell calling, multiplex mapping, and analysis using the sequences are described herein.

[0008] Barcode processing according to the system of the present invention can use a hierarchical barcode system, optionally with each level including one of a specified list of possible barcodes. A level can refer to a level of information specific to a barcode, such as a cell barcode followed by a UMI. There can also be higher-level barcodes, such as for different patients or samples, and / or intermediate-level barcodes, such as gene-specific barcodes (e.g., if sequence-specific primers can be used). Typically scRNA-Seq will include cell and UMI barcodes, but may include other levels of any category of information, such as date or location. Including barcodes in reagent oligonucleotides results in sequence data including sequences from barcodes. The computer system can match the barcodes to a list or database to assign sequence reads to cells or individual molecules. By matching each level individually, barcode matching can be done efficiently. Some embodiments of the computer system of the present invention allow a Hamming distance of 1. Limiting the Hamming distance to 1 means that the base in the sequence read corresponding to each hierarchical position (e.g., "Read 1" ("R1") within a FASTQ file) can differ from a barcode by one base and still match that barcode. In order to simplify and reduce the size of the barcode list sent to a reference mapping module such as STAR, the hierarchical barcodes can be replaced with a set of newly generated barcodes. For example, each barcode can be replaced with a number (such as its numerical position on the barcode list). In some embodiments, new barcodes are only generated for barcodes that have at least one matching read. That is, if trillions of barcodes are calculated to have formed the input to an experiment, and the system processes the resulting FASTQ file and identifies only billions of barcodes in the reads, then this new, limited list of barcodes can be assigned new barcodes, thereby greatly saving downstream computing resources. The method of barcode processing helps keep downstream steps computationally tractable while remaining able to provide a measure of expression levels when reads are mapped to a reference and UMI counts are used to de-duplicate reads.

[0009] The use of PIPs (where a large number of cells are isolated into aqueous partitions simultaneously) reveals many challenges first encountered by PIPs related to throughput and data volume, as compared to microfluidic platforms (where a limited number of cells are isolated into partitions serially). For example, the use of PIPs reveals that it would be beneficial to have automated and accurate tools for correctly calling out whether a partition contains cells, as opposed to background partitions (also referred to as cell calls) that are incompletely captured cells. The present disclosure recognizes that, as a practical matter, partitions containing cells and background partitions (incompletely captured cells) can be automatically distinguished at the data volume and throughput associated with PIPs. For example, from the sequence data, all barcodes identified can be sorted by the number of UMIs associated with each barcode. Barcodes associated with a relatively high number of UMIs are more likely to be generated from cells rather than from background partitions, while barcodes associated with a relatively low number of UMIs are more likely to be generated from background partitions. The system of the present invention can provide the user with a choice of sensitivity / specificity (e.g., selecting one of five preset options). Following these scenarios, the computer system may compute the number of UMIs (e.g., on the y-axis) versus the rank of the sorted barcodes (x-axis), fit a function to a curve, take derivatives to find inflection points, and populate the user interface with points selected from those along the curve. The user may select a sensitivity level for cell detection, and the analysis proceeds at the selected sensitivity. For example, a user trying to classify a cell type that is abundant in a tissue sample may select a low sensitivity (willing to discard some cell data in favor of excluding substantially all background partitions). 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 may provide guidance to the user on selecting a sensitivity. After any cell detection steps, the system performs barcode processing.

[0010] When sequence reads are mapped to more than one position in a reference, the present invention uses multiple mappings to preserve information. The sequencing step of RNA-Seq typically provides multiple sequence reads, typically in FASTQ format. The analysis of those reads typically involves mapping each read to a reference, such as a list of mRNA or gene sequences specifically for expression analysis, or even just a published human genome. Mapping reads to a reference typically includes some version of an alignment algorithm, such as a pairwise alignment, an alignment of a Burrows-Wheeler transformation, a comparison of a string's k-mer with a hashed k-mer of a reference, or an alignment using a suffix or prefix tree constructed from a reference and / or reads, or any combination thereof. There are a variety of tools in the art that can be used to align according to such methods (see, e.g., Dobin, 2013, "STAR: ultrafast universal RNA-seq aligner", Bioinformatics, Vol. 29, No. 1: pp. 15-21, incorporated by reference). The alignment finds the best matching position within the reference for each read and, by implication, the gene (represented in the reference) that produces the RNA transcript (represented in the read). There are many reasons in biology and informatics why a read may be mapped to a reference in more than one position. For example, a reference may include homologous genes (e.g., from repeats) or copy number variations, or a read may be short enough to allow two or more different mappings to a reference. Prior art methods simply discard such ambiguous mapping reads. In contrast, the system of the present invention treats such results as multiple mappings of the read and saves the read as a multiple mapping read. Saving information about multiple mapping reads allows analysis of scRNA-Seq sequencing data to discover biological information about homologous gene expression or copy number variation in a sample.

[0011] After barcode processing, cell detection and read mapping, the system of the present invention provides tools for analyzing the results of scRNA-Seq data. The analysis according to the present invention generally includes clustering of cells, differential expression analysis and extraction of various metrics, especially to provide, for example, a measure of the expression level of various genes in each cell. The system of the present invention provides a variety of tools to help display, view and interpret. For example, clustering methods can provide classification of cells. In addition, the system of the present invention can also be used to output basic metrics, barcode sorting diagrams, clustering diagrams, differential expression tables or combinations thereof. The preferred embodiment of the system operates in a computer system, which includes 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. The software module that performs the functions can be developed in any suitable environment, such as python, Ruby on rails, c++, other or a combination thereof. The system version of the present invention can be executed in mac osx, Linux and windows operating systems. The preferred embodiment operates in a server or cloud environment and is operable to receive sequence data from sequencing resources such as next generation sequencing (NGS) instruments or genome facilities, for example in FASTQ or FASTA format. The system can align the reads to a reference, generate a sequence alignment map (SAM) or a binary alignment map (BAM), and provide output to a user computer.Details and variations of embodiments of the system's functionality 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, the memory subsystem including instructions executable by the processor to cause the system to perform the following operations: receiving sequence data generated by sequencing RNA from a plurality of partitions; for each partition, associating counts of barcode sequences in the sequence data with the probability that the partition contains a complete isolated cell; receiving a user's selection of a sensitivity for the probability that the partition contains a complete isolated cell; and analyzing the mRNA level in the partition that meets the user's selected sensitivity.

[0013] The system can be used to provide at least three (e.g., 5) different calculated 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 associates a barcode sequence with a probability by calculating a function of the barcode count to the barcode sort, dividing the function by a predetermined value and selecting a cutoff level of the barcode count, above which the partition is considered to contain a complete isolated cell. The system can be used to reanalyze the mRNA level with a new user selection of sensitivity after analyzing the mRNA level. In certain embodiments, analyzing the mRNA level includes assigning sequence reads to cells using cell barcodes, de-repeating sequence reads using universal molecular identifiers, mapping de-repeating reads to references, and using the de-repeating read counts of genes mapped to the reference as a measure of the expression level of the gene. The system can be used to downsample the sequence data by mapping less than all de-repeating reads to a reference. In some embodiments, the number of mapped de-repeating reads is selected so that the mRNA level from the sequence data will be normalized to the level calculated from at least one other experiment.

[0014] Aspects of the present invention provide a system for expression analysis with multiple mappings. Such a system includes a processor coupled to a memory subsystem, the memory subsystem including instructions executable by the processor to cause the system to perform the following operations: receiving sequence data (e.g., Illumina output) generated by sequencing RNA from a single cell; mapping at least one sequence read from the sequence data to a reference comprising reference gene information; identifying at least a first position and a second position in the reference, wherein the sequence read is mapped with at least a threshold match score; storing the sequence read in the memory subsystem with a tag that identifies the read as being mapped to at least the first position or the second position in the reference. Preferably, the system is used to map multiple reads to the reference, and select the first position or the second position of the sequence read in the reference based on the position to which a larger number of the multiple reads are mapped. The system can be used to store the sequence read with a tag that identifies the read as being mapped to both the first position and the second position. In some embodiments, the system assigns a first weight to the mapping to the first position and a second weight to the mapping to the second position. The first weight and the second weight can each be based at least in part on a respective first alignment score and a second alignment score between the sequence read and the reference.The system can be used to provide a report describing gene duplication or copy number variation in a single cell using marker and reference gene information.

[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, the memory subsystem including instructions executable by the processor to cause the system to perform the following operations: receiving sequence data generated by sequencing RNA from a single cell; comparing barcodes from sequence reads in the sequence data to a barcode whitelist; omitting further analysis of a sequence read when a barcode in the sequence read fails to match a whitelist barcode with a predetermined Hamming distance value, wherein for each sequence read, the barcode from the first level is compared to the cell barcode whitelist, and the barcode from the second level is compared to the UMI whitelist; deduplicating reads whose first and second level barcodes match the whitelist within the predetermined Hamming distance value; and providing a measure of mRNA levels from the deduplicated reads. The system can be used to replace the barcodes in the sequence reads with index values ​​that occupy less space in the memory subsystem before the deduplication step.

[0016] For each of the foregoing aspects, the present invention provides a corresponding method, including using the system to perform the functions. Thus, the present invention provides systems and methods for evaluating gene expression levels from scRNA-Seq experiments. The system and method use a computer simulated cell detection tool to distinguish between partitions containing cells and background partitions that do not fully capture cells. The system and method also provide a barcode processing tool that tests barcodes from at least the cell level and the UMI level for whitelisting and optionally converts these barcodes to a compact index format. In addition, the system and method can implement a multi-mapping tool to save information about reads mapped to multiple locations in a reference. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 Barcode sequencing maps are shown.

[0018] Figure 2 Barcode sequencing images of HEK / 3T3 cells are shown.

[0019] Figure 3 Barcode sequencing plots of peripheral blood mononuclear cells (PBMCs) are shown.

[0020] Figure 4 Barcode sequencing maps of breast tissue samples are shown.

[0021] Figure 5 Cell detection at five sensitivity levels for PBMC samples is shown.

[0022] Figure 6 Clustering results for HEK / 3T3 samples are shown.

[0023] Figure 7 Clustering results for PBMC samples are shown.

[0024] Figure 8 is a barnyard plot showing the number of transcripts by organism type.

[0025] Fig. 9 A system of the present invention is shown. DETAILED DESCRIPTION

[0026] 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 provides a comprehensive analysis solution, providing 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 flows.

[0027] System requirements

[0028] The preferred embodiment of the system runs on Windows, macOS, and Linux operating systems.

[0029] • Preferably, the system includes at least about 64 GB of RAM and 1 TB of hard disk space (to accommodate multiple results).

[0030] - In some embodiments, the version of the system that runs on Windows uses a native installation of Docker or a similar program, which can allow the system to interoperate seamlessly with modules that are not natively compatible with Windows.

[0031] Pipeline Overview

[0032] The system of the present invention may perform any combination of the following steps:

[0033] 1) FASTQ processing, including barcode matching, QC, and (optionally) downsampling

[0034] 2) Mapping

[0035] 3) UMI Count

[0036] 4) Cell detection

[0037] 5) Clustering

[0038] 6) Differential expression

[0039] 7) Metric Extraction

[0040] 8) Report Generation

[0041] Barcode processing

[0042] The sequencing process of scRNA-Seq produces multiple sequence reads, typically in one or more FASTQ files. One benefit of PIP (over traditional microfluidics) is that a large number of cells can be isolated in a partition at the same time and then interrogated to generate a single FASTQ file (without concatenating FASTQ files from different runs / days). Processing can begin by processing the FASTQ file to: extract reads associated with the barcodes; remove unwanted technical sequences from the data before mapping; and optionally downsample the data. A module of the system (e.g., optionally developed in python, C++, or other such environments) can read in from the FASTQ file and process each entry to read the barcode sequence and the RNA sequence. Typically, the barcode sequence is compared to a list of barcodes, and the RNA sequence data is mapped to a reference to identify a gene or transcript.

[0043] In certain embodiments, barcode processing according to the system of the present invention may use a hierarchical barcode system, wherein each level optionally includes one of a specified list of possible barcodes. A level may refer to a level of information specific to a barcode, such as a cell barcode followed by a UMI. There may also be higher-level barcodes, such as for different patients or samples, and / or intermediate-level barcodes, such as gene-specific barcodes (e.g., if sequence-specific primers can be used). Typically scRNA-Seq will include cell and UMI barcodes, but may include other levels of any type of information, such as date or location. Including barcodes in reagent oligonucleotides results in sequence data that includes sequences from barcodes. A computer system may match barcodes to a list or database to assign sequence reads to cells or individual molecules. Barcode matching may be done efficiently by matching each level individually. Some embodiments of the system of the present invention allow a Hamming distance of 1. Constraining the Hamming distance to 1 means that the base in the sequence read corresponding to each hierarchical position (eg, "read 1" ("R1") within a FASTQ file) can differ from the barcode by one base and still match that barcode.

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

[0045] Quality Control

[0046] In some embodiments, the R2 (read 2) fastq file includes the cDNA construct and also includes technical sequences related to the PIP chemistry, such as template switching oligonucleotides (TSO) and poly-A sequences. The occurrence of these sequences may be negatively correlated with the fragment size. Because the read is less likely to be mapped if it contains foreign sequences, it may be preferred to remove the TSO sequence from the 5' end and the poly-A sequence from the 3' end in the R2 fastq file. In addition, depending on the chemical type, a fixed number of non-cDNA bases can be removed from the 5' end. Reads that are less than 20 bases long after trimming can be discarded from further analysis.

[0047] After all reads are matched to barcodes, inevitably some barcodes will be associated with a small number of reads. Those barcodes are highly unlikely to be classified as cells later, so it may be economical to exclude them and the associated reads from further analysis. The user can define the minimum number of reads for each barcode. If such a threshold is specified, reads associated with low count barcodes will not be passed to STAR for mapping. See "Cell Detection" discussed in more detail elsewhere in this article.

[0048] Downsampling

[0049] In some experiments, it may be desirable to normalize the sequencing depth so that experimental conditions can be reliably compared. Therefore, the system of the present invention is used to downsample the data to a specified number of reads. The reads can be randomly selected, with the probability of selection being determined by the total number of reads (provided by the user or determined by counting the number of lines in the fastq input).

[0050] Output

[0051] In certain embodiments, the fastq processing step produces the following output:

[0052] - <output-root> / metrics / barcode_stats.json: General statistics about the quality of barcode matching.

[0053] - <output-root> / metrics / barcodes / Generated_barcode_read_info_table:

[0054] A list of all generated barcodes that are passed to STAR for mapping and counting. The list includes the level associated with the original barcode, and the number of perfect matches and error correction matches for each barcode.

[0055] - <output-root> / metrics / barcodes / barcode_whitelist: A list of all generated barcodes that are provided to STAR as a "whitelist" for mapping and counting.

[0056] - <output-root> / barcodes_fastqs / *: Intermediate fastq files used by STAR for mapping. These files are deleted at the end of the analysis unless otherwise instructed.

[0057] Mapping

[0058] The sequencing step of RNA-Seq generally provides multiple sequence reads, generally in FASTQ format. The analysis of those reads generally involves mapping each read to a reference, such as a list of mRNA or gene sequences specifically for expression analysis, or even just a published human genome. Mapping the reads to a reference generally includes a certain version of an alignment algorithm, such as a comparison of pairs, a Burrows-Wheeler transformed comparison, a comparison of a string's k-mer with a hashed k-mer of a reference, or one or any combination of comparisons using a suffix or prefix tree constructed from a reference and / or read. Various tools are available in the art for alignment according to such methods. Certain preferred embodiments use a read mapping module known as STAR. See Dobin, 2013, "STAR: ultrafast universal RNA-seq aligner", Bioinformatics, Vol. 29, No. 1: pp. 15-21, incorporated by reference. Specifically, the system of the present invention uses STAR (comparison of spliced ​​transcripts with references) and the related STARsolo to map reads to transcriptomes and quantify transcript abundance. STAR uses a special index constructed from the original transcriptome of interest. Some common indexes such as humans and mice can be downloaded from sources such as Fluent BioSciences. However, for some special cases, the system of the present invention allows users to create unique or experimentally specific indexes.

[0059] Documentation for STAR products is available from the publisher. In a preferred embodiment, the system of the present invention stores the output of read mapping (e.g., via STAR) in a repository such as <output-root> / starsolo directory, for example, as described in the STAR documentation.

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

[0061] Read mapping generally refers to detecting the position within a reference to which a sequence read maps best. The system includes programmed logic for using information potentially revealed when a sequence read maps well to more than one position on a reference.

[0062] The present invention uses multiple mappings to preserve information when a sequence read is mapped to more than one position in a reference. There are many reasons in biology and informatics to explain why a read can be mapped to a reference in more than one position. For example, a reference may include homologous genes (e.g., from duplications) or copy number variations, or the read may be short enough to allow two or more different mappings to a reference. Prior art methods simply discard such ambiguous mapping reads. In contrast, the system of the present invention treats such results as multiple mappings of the read and saves the read as a multiple mapping read. Preserving information about multiple mapping reads allows analysis of scRNA-Seq sequencing data to discover biological information about homologous gene expression or copy number variation in a sample. Importantly, it is preferred that multiple mapped reads are not discarded (compared to prior art systems that discard such reads). In some embodiments, reads are assigned to one of the positions according to a selection logic, such as by random assignment, or (in the case of two positions) alternating between the two positions when multiple reads map to the two positions, or by mapping each read to the position to which the majority of reads map, or by giving "partial credit" (e.g., proportional to the number of different positions, or weighted by alignment score), or otherwise, or a combination thereof.

[0063] UMI Count

[0064] After STAR maps the reads in the fastq file to the reference, it creates a .bam file (a binary alignment map (BAM), which is a binary version of a 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 that are mapped to the reference. The molecular information is saved to a file in HDF5 format, for example, with an appropriate name such as molecule_info.h5. The system can analyze this file to create a raw count (feature barcode) matrix. In some embodiments, the matrix is ​​contained in the following files: matrix.mtx.gz, barcodes.tsv.gz, and features.tsv.gz, for example, all in a directory such as raw_matrix. These files form a sparse matrix that contains the complete counts 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.

[0065] Cell detection

[0066] After obtaining the UMI count matrix, the barcodes are separated into PIPs containing cells and background PIPs (i.e., PIPs that do not fully capture cells). This separation is based on a barcode ranking graph that ranks the barcodes based on the number of unique UMIs associated with the barcodes. The use of PIPs (where a large number of cells are isolated simultaneously into aqueous partitions) as compared to microfluidic platforms (where a limited number of cells are isolated sequentially into partitions) reveals many challenges first encountered by PIPs related to throughput and data volume. For example, the use of PIPs reveals that it would be beneficial to have automated and accurate tools for correctly calling out whether a partition contains cells, as opposed to background partitions (also known as cell calls) that are incompletely captured cells. The present disclosure recognizes that, as a practical matter, at the data volume and throughput associated with PIPs, distinguishing between partitions containing cells and background partitions (incompletely captured cells) would benefit from being done automatically, or even "on the fly." For example, from the sequence data, all barcodes identified can be sorted by the number of UMIs associated with each barcode. Barcodes associated with a higher number of UMIs are more likely to arise from cells, while barcodes associated with a lower number of UMIs are more likely to arise from background partitions.

[0067] Figure 1 A barcode sorting plot is shown. Typically, cellular barcodes are concentrated in the top pattern of the sorting plot, while background barcodes are concentrated in the low UMI count pattern. Therefore, the goal of cell detection is to find a point in the first "inflection" region that separates the two patterns. One method of finding this point involves dividing the number of UMIs at the top of the sorting plot by a constant denominator (e.g., 10). Methods for finding the inflection point may be highly sensitive to the shape of the sorting plot, which varies greatly between different sample types. The system of the present invention can automatically find the point on the barcode sorting plot that separates the cellular partition from the background partition, or the system can provide one or more options to the user to help the user set or select this point, or the system can perform a combination of automatic selection or suggestions in conjunction with user selections.

[0068] The system of the present invention can provide the user with a sensitivity / specificity selection (e.g., one of five preset options). The graphs presented here show the cell barcode detection results for HEK / 3T3, PBMC, and breast tissue samples, up to 1 / 10 of the upper limit of the number of barcodes.

[0069] Figure 2 Barcode sequencing images of HEK / 3T3 cells are shown.

[0070] Figure 3 Barcode sequencing plots of peripheral blood mononuclear cells (PBMCs) are shown.

[0071] Figure 4 Barcode sequencing maps of breast tissue samples are shown.

[0072] The cell fraction in each plot is highlighted in bold. Cell detection for HEK / 3T3 appears consistent with visual inspection of the barcode sorting plot: the entire high UMI pattern is captured. However, in the case of PBMCs and breast cells, the number of cell barcodes is clearly underestimated, as the top pattern has a more negative slope and is less well defined. This may be a choice made intentionally by the user. For example, a user studying PBMCs compared to HEK / 3T3 may need to optimize the data volume determination (largely ignoring or excluding background PIPs) depending on the quality of information sought due to underlying biological reasons.

[0073] In some embodiments, to accommodate differences between samples and facilitate accurate cell detection, the system of the present invention provides cell detection at multiple (e.g., five) different sensitivity levels using the following steps:

[0074] 1) Select a starting point S near the top of the ordination graph. To avoid selecting outliers, the system uses the point where the slope of the graph first falls below 10% as the starting point.

[0075] 2) Let M be the number of UMIs of the barcode at the starting point. For each cell detection sensitivity level L, find the number of UMIs U such that:

[0076]

[0077] 3) For each cell detection sensitivity level Li, all barcodes with UMI counts higher than or equal to Ui are selected as barcodes containing cells.

[0078] The system may use any suitable method to determine the point comprising the turn or inflection point. For example, behind the scene (or on a screen, e.g., graphically displayed to a user), the computer system may compute the number of UMIs (e.g., on the y-axis) versus the order of the ordered barcodes (x-axis), fit a function to a curve, take the derivative to find the inflection point, and populate a user interface with options for selecting a point along the curve. The user may select a sensitivity level for cell detection, and the analysis proceeds at the selected sensitivity. For example, a user trying to classify a cell type that is present in large quantities in a tissue sample may select a low sensitivity (willing to discard some cell data in favor of excluding substantially all background partitions). 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 may provide guidance to the user regarding selecting a sensitivity.

[0079] Figure 5 Cell detection at five sensitivity levels for PBMC samples is shown.

[0080] By its very nature, cell detection involves a trade-off between the number of cells and their quality. One can arbitrarily increase the number of cells by detecting the sorted graph more deeply, but this increases the possibility of selecting low RNA content or otherwise damaged cells. The optimal number of cells recovered depends on the sample type and purpose of each experiment. For example, if some cell types in the sample (such as neutrophils or erythrocytes) have low RNA expression, they may be present in high numbers in the lower UMI areas of the graph, or the UMI count distribution may even be unimodal. In this case, high sensitivity cell detection is guaranteed. On the other hand, if high clustering accuracy is required to accurately characterize cell types, it is better to use low sensitivity cell detection. Cell detection at different sensitivity levels allows users to observe a variety of possible results and determine the correct sensitivity for their data.

[0081] Manual selection of cell number

[0082] In some cases, the user may want more precise control over the number of cells detected, beyond the range provided by the five sensitivity levels. In some embodiments, the system allows the user to specify a desired number of cells, in which case that exact number of barcodes is selected from the top of the sorting graph. The output described herein will be generated for the selected set of barcodes. For some embodiments of the system of the present invention, this can be performed from the command line, for example, using the prefix force_N, where N is the number of cells entered by the user.

[0083] Output

[0084] After cells are detected, directories can be created for each sensitivity level within a parent directory, such as:

[0085] <output-root> / cell_calling, contains a text file listing the indices (starting from 1) of the barcodes selected as cells and an image with a barcode ranking plot. In addition, a filtered count matrix can be generated for each sensitivity level and can be used for downstream analysis. The filtered matrix is ​​stored in <output-root>Those skilled in the art will understand how to write functions in a suitable programming or development environment such as python to save those outputs.

[0086] analyze

[0087] In addition to barcode processing, cell detection and read mapping, the system of the present invention also provides tools for analyzing the results of scRNA-Seq data. Analysis according to the present invention generally includes clustering of cells, differential expression analysis and extraction of various metrics, especially to provide, for example, a measure of the expression level of various genes in each cell. The system of the present invention provides a variety of tools to help display, view and interpret. For example, clustering methods can provide tools useful for cell classification.

[0088] Clustering

[0089] After selecting cell barcodes, cluster analysis may be required to quantify and visualize heterogeneity within the cell population and identify different cell types. For each group of cells detected (e.g., different sensitivities or user-imposed numbers), the system of the present invention can create a cluster graph using methods such as K-means clustering with different K values ​​and using Leiden (graph-based) clustering, which automatically determines the number of clusters based on nearest neighbors.

[0090] Preprocessing steps

[0091] The starting point for clustering is preferably a sparse matrix representation of the filtered count matrix, where each row represents a cell and each column represents a gene. Before running the actual clustering algorithm, several processing steps can be performed:

[0092] 1) Cell normalization: Normalize the data of each cell to a unit standard. The purpose of this is to allow cells with similar expression profiles but different relative RNA abundances to still be clustered together.

[0093] 2) Logarithmic transformation: Use the ln(x+1) function to transform the matrix so that it is closer to a normal distribution.

[0094] 3) Highly variable gene selection: A subset of genes with high variance is selected to maximize the efficiency of clustering. For each gene, the expression variance between cells is calculated. The genes are then ranked by variance, and the top ranked genes at the Nth percentile are selected, where N is a user input.

[0095] 4) Scaling: According to the convention of cluster analysis, it is preferred to scale the feature columns (genes) so that the mean of each column is 0 and the standard deviation is 1.

[0096] 5) PCA: The system of the present invention can use principal component analysis (PCA) to reduce the dimensionality of the data from a scaled set of highly variable genes to a small, user-controlled number of components that best explain the variance in the data.

[0097] After the preprocessing step, clustering can be performed in the new PCA space.

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

[0099] In certain embodiments, the system also uses the popular K-means algorithm to divide cells into distributed clusters. The system runs the algorithm for different K values ​​representing a predetermined number of clusters. The minimum and maximum values ​​of the K value range are controlled by the user and should be customized according to the number of cell types expected to be present in the sample.

[0100] Differential expression

[0101] Once clustering is complete, it may be of interest to determine the main genes associated with each cluster. For this step, the system of the present invention can use the Wilcoxon rank sum test (Mann-Whitney U) for the independent variable. This nonparametric test is used to determine whether two samples are from the same distribution, which is achieved by arranging the values ​​associated with the two groups together and calculating the sum of those arrangements in each group. These arrangements are used to calculate the test statistic (U), which can compare clusters and non-cluster cell populations. For each cluster, the system calculates the differential expression of the genes associated with the cells in the cluster relative to the genes associated with all other cells. The top ranked genes arranged in descending order of the Z scores of each cluster are exported to a csv file. This process is performed for different K values ​​in K-means clustering and Leiden clustering based on a graph.

[0102] Output

[0103] In addition, the system of the present invention can also be used to output basic metrics, barcode ranking diagrams, cluster diagrams, differential expression tables, or combinations thereof. <output-root>A directory is created in the / clustering directory for each or any group of cells detected, which contains the following:

[0104] A csv file listing the cluster assignments for each cell from graph-based clustering and from all K-means runs.

[0105] A csv file containing the graph-based clustering and silhouette scores for all K-means runs. The silhouette score is a measure of clustering quality that takes into account both intra-cluster and inter-cluster distances. Note that it is calculated in the PCA space used for clustering, not in the UMAP space used for visualization.

[0106] Graph-based clustering and UMAP plots for all K-means runs. UMAP (Uniform Manifold Approximation and Projection) is a popular method for projecting multidimensional data onto a two-dimensional space.

[0107] CSV files of graph-based clustering and all K-means runs, containing the top genes for each cluster.

[0108] A csv file similar to the top genes file, but includes additional statistics for each gene: Z-score, p-value, adjusted p-value, and log2 fold change.

[0109] The figure below shows the clustering results for HEK / 3T3 samples and PBMC samples. The number of clusters shown in each case is the number of clusters that produced the highest silhouette score.

[0110] Figure 6 Clustering results for HEK / 3T3 samples are shown.

[0111] Figure 7 Clustering results for PBMC samples are shown.

[0112] Metrics Extraction

[0113] In the final step of the analysis, various metrics are calculated and reported for each sensitivity level (or for a single set when the user specifies the number of cells).

[0114] In some embodiments, the system of the present invention may include one or any combination of multiple basic metrics. For example, one or more of the following metrics may be included in each analysis:

[0115] Total number of reads in the input fastq file

[0116] • 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).

[0117] The number of cells detected

[0118] The average number of reads per cell (this is the total number of input reads divided by the number of cells).

[0119] • The percentage of reads in cells (this is the number of reads associated with barcodes containing cells divided by the total number of input reads).

[0120] The UMI duplication rate in the cell (this is the number of mapped reads in the cell divided by the number of unique reads (duplicate reads from the same transcript) in the cell. For example, if each UMI has three reads containing it, i.e., two deduplicated reads per UMI, the UMI duplication rate would be 3).

[0121] Total number of UMIs in the cell

[0122] The median UMI per cell (this is based on the distribution of unique UMI (transcript) counts associated with each cell barcode).

[0123] • The number of genes expressed in the cell (this number includes all genes with non-zero counts in the filtered matrix, and gives each gene the same weight regardless of expression level).

[0124] - Median genes per cell (this number includes all genes with non-zero counts in the filtered matrix, and gives equal weight to each gene regardless of expression level).

[0125] Percentage of cells in which more than 50% of transcripts are derived from mitochondrial genes

[0126] Percentage of cells in which less than 5% of transcripts are derived from mitochondrial genes

[0127] Farm metrics and estimation of multicellularity rates

[0128] Embodiments of the system provide tools for evaluating the quality of scRNA-Seq protocols (such as Drop-seq or inDrop) or a specific chemistry or reaction setup in one of these protocols. One way to evaluate the quality of single-cell RNA technology is a mixed species ("farm") experiment, typically involving a mixture of human and mouse cells. Such an experiment can provide information about the occurrence rate of multiple cells 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. The more human genes are found in mouse cells (and vice versa), the more background contamination there may be in the system.

[0129] When the reference transcriptome is a mixed human and mouse transcriptome, the system of the invention is used to provide an additional set of metrics specific to farm experiments:

[0130] • Number of human cells (this is the number of cell barcodes with at least 85% of transcripts from a human reference).

[0131] • Number of mouse cells (this is the number of cell barcodes with at least 85% of transcripts from the mouse reference).

[0132] Number of multiplex cells (this is the number of cell barcodes that have less than 85% of their transcripts from the human or mouse reference; note that a multiplex may consist of two human cells or two mouse cells, or (rarely) include more than two cells, so this may underestimate the true number of multiplex cells).

[0133] The total number of UMIs in human cells

[0134] Median human UMI per cell

[0135] The number of human genes expressed in human cells

[0136] Median human genes per human cell

[0137] Percentage of mouse UMIs in human cells

[0138] Total number of UMIs in mouse cells

[0139] Median mouse UMI per mouse cell

[0140] The number of mouse genes expressed in mouse cells

[0141] Median mouse genes per mouse cell

[0142] Percentage of human UMIs in mouse cells

[0143] An important output of farm analysis is called a farm map.

[0144] Figure 8 Here is a farm plot showing the number of mouse transcripts versus the number of human transcripts, with different shades representing human cells, mouse cells, and multiplex cells:

[0145] Output

[0146] After all metrics have been calculated, the results for each set of detected cells can be reported in a file such as:

[0147] Basic metrics are stored in csv and json formats <output-root> / metrics.

[0148] Farm metrics are stored in csv format and json format if applicable <output-

[0149] root> / barnyard.

[0150] Analytical Report

[0151] The operating system of the present invention for analyzing sequence data from scRNA-Seq experiments optionally generates reports, such as documents in formats such as XML or HTML. Such reports can be stored in appropriate directories, such as:

[0152] <output-root> / report and can include basic metrics, barcode ordination plots, clustering plots, and differential expression tables. The user can switch between different clustering modes (graph-based clustering mode and K-means clustering mode with different K values), which will change the contents of the mapping and tables.

[0153] If the "farm" mode is used, those statistics and farm plots will also be included.

[0154] In standard analysis, a report can be generated for each sensitivity level. An additional combined report (e.g., "combined.html") can include all sensitivity levels. The user can choose to browse the different sensitivities in the combined report and select the best sensitivity level or decide on the mandatory number of cells to be used in subsequent re-analysis.

[0155] If a cell number is manually selected (forced) during analysis, a report will be generated for that specific number and a new tab may be added to the combined report.

[0156] Running the pipeline

[0157] The preferred embodiment of the system operates in a computer system, and the computer system includes at least one processor coupled to a memory subsystem. The function provided by the system can be implemented in a local computer, a server system, or a cloud-based computing system. The software module that performs the function can be developed in any suitable environment, such as python, Ruby onrails, C++, other or their combination. The system version of the present invention can be executed in mac osx, Linux and windows operating systems. The preferred embodiment operates in a server or cloud environment, and can operate to receive sequence data from sequencing resources such as next generation sequencing (NGS) instruments or genome facilities, for example, in FASTQ or FASTA format. The system can read segments and reference comparisons, generate sequence alignment maps (SAM) or binary alignment maps (BAM), and output is provided to user computers. The details and variations of the embodiments of the system functions are described herein.

[0158] Fig. 9 A system 901 is shown that includes at least one local computer 905. The system 901 may also include a server or cloud computer 909. Preferably, any local computer 905 or server or cloud computer 909 included in the 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 the system 901 is used to receive sequence data from a sequencing instrument 921 or a genome facility via a communication network 917, through which the machines of the system 901 communicate and interoperate with each other. The memory subsystem in the local computer 905 or server or cloud computer 909 of the system 901 preferably includes program instructions that can be executed by a processor to cause the system 901 to perform the analysis functions described herein.

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

[0160] See below for a description of all required and optional parameters for certain embodiments of the disclosed system. To view a complete list of parameters from the command line, run:

[0161] $pipseeker count --help.

[0162] Required parameters

[0163] $ pipseeker count --input-path<path to fastq files> --output-root<destination for all outputs> --STAR-index-path<path to mapping reference>

[0164] --input-path(string)

[0165] Path to the source fastq files. If this is a directory name, all files with names ending in .fastq or fastq.gz will be included as input. Alternatively, if there are files from multiple samples in one directory, a prefix may be specified at the end of the path, which will restrict the included files to only those with names starting with that prefix.

[0166] For example, consider the case where the input directory contains two-channel results from two different samples:

[0167] $ls my_output_dir

[0168] sample_1_L001_R1.fastq.gz sample_1_L001_R2.fastq.gz

[0169] sample_1_L002_R1.fastq.gz sample_1_L002_R2.fastq.gz

[0170] sample_2_L001_R1.fastq.gz sample_2_L001_R2.fastq.gz

[0171] sample_2_L002_R1.fastq.gz sample_2_L002_R2.fastq.gz

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

[0173] --output-root(string)

[0174] The directory where all PIPseeker output will be stored. This can be the same as the input path, or a different directory.

[0175] --star-index-path(string)

[0176] The path to a directory containing index references for mappings.

[0177] Optional parameters

[0178] generally:

[0179] --chemistry(string,default:mark6)

[0180] Allows the user to enter the chemical used to generate the sample. Accepted values ​​are implementation specific (e.g., Marker 3 and Marker 6).

[0181] --downsample-to(integer)

[0182] An integer representing the target number of reads to be mapped if downsampling is desired. If specified, reads will be selected randomly and only the selected subset will be sent to STAR for mapping.

[0183] --input-reads(integer)

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

[0185] --input-reads should include the total number of reads across all input files combined. If not specified, reads will be counted manually, so using this parameter can reduce run time.

[0186] --verbosity(integer,default:1)

[0187] Verbosity level: 0 (silent), 1 (terse), or 2 (verbose).

[0188] Barcoding

[0189] --min-reads-per-barcode(integer,default:1)

[0190] Minimum number of reads associated with a barcode. If a barcode is associated with less than the specified number of reads, all reads from that barcode will be discarded without further analysis. Increase this number to reduce mapping time by not processing barcodes that are unlikely to be cells. Note that low count portions of the barcode ranking graph will be incomplete and mapping statistics may change slightly.

[0191] --retain-barcoded-fastqs

[0192] Flag to keep intermediate barcode fastqs. In most cases it is not needed. If not specified, those fastqs will be discarded after finishing the analysis. Note that the input fastqs are always kept regardless of this flag.

[0193] Mapping

[0194] --star-threads(integer,default:8)

[0195] Number of CPU threads to use for STAR mapping (default: 8).

[0196] Cell detection

[0197] --force-cells(integer)

[0198] Force a specific number of barcodes to be considered cells. If specified, that number of barcodes will be selected from the top of the barcode sorting graph. This will override the default, sensitivity-based cell calling mode.

[0199] Clustering

[0200] --clustering-percent-genes(integer,default:10)

[0201] The percentage of genes used for clustering. Genes with the highest expression variability within this percentile ranking will be included.

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

[0203] The number of components for PCA dimensionality reduction.

[0204] --percent-genes(float,default:10)

[0205] The percentage of genes selected for clustering based on the highest variance in expression.

[0206] --min-clusters-kmeans(integer,default:2)

[0207] The minimum number of clusters to use for K-means clustering. The value of K will increase from this number.

[0208] --max-clusters-kmeans(integer,default:12)

[0209] The maximum number of clusters to use for K-means clustering. The value of K will be incremented until this number is reached.

[0210] --nearest-neighbors(integer,default:10)

[0211] The number of nearest neighbors for graph-based clustering. This defines the minimum number of cells per cluster.

[0212] --diff-exp-genes(integer,default:50)

[0213] The number of the most differentially expressed genes included in each cluster.

[0214] measure

[0215] --run-barnyard

[0216] Flag that enables "farm" metrics to be generated and included in analysis reports.

[0217] Report

[0218] --id

[0219] The sample ID to include in the output report.

[0220] --description

[0221] A sample description is included in the output report.

[0222] Re-run the final stage: Reanalysis

[0223] After observing the results of the complete analysis, the user may wish to change parameters related to the later stages, such as the range of K values ​​used for K-means clustering. The user may also wish to force a certain number of barcodes to be assumed to represent the cell. This should not require repeating the barcoding and mapping steps, which are usually the most time-consuming.

[0224] In a preferred embodiment, using the reanalysis feature will skip barcoding and mapping and start analysis directly when cells are detected. If this mode is used, the input path should point to a directory containing the complete analysis results from a previous run. --output-root, --chemistry. Parameters related to barcoding and mapping will be ignored.

[0225] For example, suppose that after running the full pipeline and visually inspecting the results, none of the standard cell call sensitivity levels appear to capture the correct area of ​​the barcode sequencing map. The user can force a specific number of cells without rerunning the entire pipeline:

[0226] $pipseeker reanalyze--input-path my_analysis_result_dir--force-cells5000

[0227] Using the described commands, hardware, and functions, the present invention provides a system comprising a processor coupled to a memory subsystem, the memory subsystem comprising instructions executable by the processor to enable the system to evaluate gene expression levels from a scRNA-Seq experiment. The system uses a computer simulated cell detection tool to distinguish between partitions containing cells and background partitions that do not fully capture cells. The system also provides a barcode processing tool that tests barcodes from at least the cell level and the UMI level for whitelisting and optionally converts these barcodes to a compact index format. In addition, the system can implement a multi-mapping tool to save information about reads mapped to multiple locations in a reference. The system also includes analysis tools for clustering, cell classification, and providing reports, including, for example, measures of expression levels. The functionality of the system solves the problems posed by scRNA-Seq using pre-templated instant partitioning (PIP), in which a large number of cells are simultaneously separated into droplets at a speed that cannot be achieved with traditional microfluidics.

Claims

1. A system for expression analysis using cell detection, the system comprising a processor coupled to a memory subsystem, the memory subsystem comprising instructions executable by the processor to cause the system to perform the following operations: receiving sequence data generated by sequencing RNA from a plurality of partitions; For each partition, correlating the count of the barcode sequence in the sequence data with the probability that the partition contains a complete isolated cell; receiving a user selection of a sensitivity for the probability that the partition contains an intact isolated cell; and The mRNA levels in the partitions that meet the user-selected sensitivity are analyzed.

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

3. The system of claim 1, wherein the system associates the barcode sequence with the probability by calculating a function of barcode counts versus barcode order, dividing the function by a predetermined value, and selecting a cutoff level of barcode counts above which a partition is considered to contain a complete isolated cell.

4. The system of claim 1, further configured to reanalyze the mRNA levels after analyzing the mRNA levels using a new user selection for the sensitivity.

5. The system of claim 1 , wherein analyzing the mRNA levels comprises assigning sequence reads to cells using cellular barcodes, deduplicating sequence reads using universal molecular identifiers, mapping the deduplicated reads to a reference, and using the deduplicated read counts of genes mapped to the reference as a measure of expression levels of the genes.

6. The system of claim 5, further configured to downsample the sequence data by mapping less than all of the deduplicated reads to the reference.

7. The system of claim 6, wherein the number of the deduplicated reads mapped is selected so that the mRNA levels from the sequence data will be normalized to levels calculated from at least one other experiment.

8. A system for expression analysis using multiple mappings, the system comprising a processor coupled to a memory subsystem, the memory subsystem comprising instructions executable by the processor to cause the system to: receiving sequence data generated by sequencing RNA from a single cell; mapping at least one sequence read from the sequence data to a reference comprising reference gene information; identifying at least a first position and a second position in the reference where the sequence read maps with at least a threshold match score; The sequence reads are stored in the memory subsystem with a tag that identifies the read as mapping to at least the first position or the second position in the reference.

9. The system of claim 8, further wherein the system is used to map multiple reads to the reference, and select a first position or a second position of the sequence read in the reference based on the position to which a larger number of the multiple reads are mapped.

10. The system of claim 8, further configured to store the sequence reads with a tag that identifies the reads as mapping to both the first location and the second location.

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

12. The system of claim 11, wherein the first weight and the second weight are each based at least in part on a respective first alignment score and a second alignment score between the sequence read and the reference.

13. The system according to claim 8, further configured to use the marker and the reference gene information to provide a report describing gene duplication or copy number variation in the single cell.

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 to cause the system to: receiving sequence data generated by sequencing RNA from a single cell; comparing barcodes from sequence reads in the sequence data to a barcode whitelist; When a barcode in a sequence read fails to match a whitelist barcode with a predetermined Hamming distance value, further analysis of the sequence read is omitted, wherein for each sequence read, the barcode from the first level is compared to the cellular barcode whitelist and the barcode from the second level is compared to the UMI whitelist; deduplicating reads of first-level and second-level barcodes that match the whitelist within the predetermined Hamming distance value; as well as A measure of mRNA levels is provided from the deduplicated reads.

15. The system of claim 14, further configured to replace the barcodes in the sequence reads with index values ​​that occupy less space in the memory subsystem prior to the deduplication step.

Images