Association mapping on single-cell RNA sequencing data

Abstract

In accordance with aspects of the present disclosure, sequencing systems and methodologies described herein may be utilized in various contexts for determining associations between features of a single-cell data set (e.g., perturbations, accessibility, genotype) and gene expression. As discussed herein, various embodiments of the present technique comprise a computational pipeline that accurately maps such associations in a computationally rapid and efficient manner. By way of example, a single-cell data sample comprising hundreds of thousands of cells or more may be processed in less than a day.

Description

ASSOCIATION MAPPING ON SINGLE-CELL RNA SEQUENCING DATAFIELD OF THE TECHNOLOGY DISCLOSED

[0001] The technology disclosed relates to the use of techniques that are implemented on computers and digital data processing systems for the purpose of analyzing single-cell RNA sequencing data to determine associations between gene expression and features present in the data set (e.g., perturbations, accessibility, genotype, and so forth).BACKGROUND

[0002] The subject matter discussed in this section should not be assumed to be prior art merely as a result of its mention in this section. Similarly, a problem mentioned in this section or associated with the subject matter provided as background should not be assumed to have been previously recognized in the prior art. The subject matter in this section merely represents different approaches, which in and of themselves can also correspond to implementations of the claimed technology.

[0003] Aspects of the present disclosure relate generally to devices, systems, and methods providing biological or chemical analysis. Various protocols in biological or chemical research involve performing a large number of controlled reactions on local support surfaces or within predefined reaction chambers. The designated reactions may then be observed or detected, and subsequent analysis may help identify or reveal properties of chemicals involved in the reaction. For example, in some multiplex assays, an unknown analyte having an identifiable label (e.g., fluorescent label) may be exposed to thousands of known probes under controlled conditions. Each known probe may be deposited into a corresponding well of a flow cell channel. Observing any chemical reactions that occur between the known probes and the unknown analyte within the wells may help identify or reveal properties of the analyte. Other examples of such protocols include known DNA sequencing processes, such as sequencing-by-synthesis (SBS) or cyclic-array sequencing.

[0004] While a variety of devices, systems, and methods have been made and used to performbiological or chemical analysis, it is believed that no one prior to the inventor(s) has made or used the devices and techniques described herein.

[0005] With the preceding in mind, the various processes in the human body are complex in nature, having both genetic and epigenetic components. By way of example, genetic variations can help explain many diseases, particularly complex or metabolic diseases that may present many and varied observable characteristics or symptoms, all of which may vary in severity between individuals. With respect to such complex diseases, every human has a unique genetic code and there are many genetic variations within a group of individuals, some of which may be linked to such complex diseases. Correspondingly, it may be difficult to identify which genes are likely to be of clinical interest in the context of a given genetic disease. This may make identifying an association between a given gene or gene variant and phenotypes of interest difficult. This in turn has implications when studying genetic-based diseases as well as when developing pharmaceutical and other treatments for such disorders.

[0006] Conventional methods for association mapping with single-cell data, which may be of interest in such contexts due to the granularity of the data and the ability to target specific spatial regions, tissue, organs, and so forth, include correlation-based and model-based methods. Correlation-based methods use a simple measure of correlation (such as Pearson or Spearman correlation) and are fast to compute but do not account for any confounding factors. Conversely, model-based methods can account for confounders, but are slow and are difficult to scale.SUMMARY

[0007] A summary of certain embodiments disclosed herein is set forth below. It should be understood that these aspects are presented merely to provide the reader with a brief summary of these certain embodiments and that these aspects are not intended to limit the scope of this disclosure. Indeed, this disclosure may encompass a variety of aspects that may not be set forth below.

[0008] The techniques described herein may be utilized in various contexts for the generation of single-cell RNA sequencing data that may be useful in determining associations between features of the data set (e.g., perturbations, accessibility, genotype, or, more generally, any characteristic of value that is part of the variation between the cells) and gene expression. By wayof context, it may be useful to associate signals associated with a feature of a dataset with the level of gene expression measured from single-cell sequencing data (scRNA-seq). Mapping such signals onto scRNA-seq data is useful with respect to a variety of assays including, but not limited to: (1) perturbation sequencing (Perturb-seq) in which a library of perturbations is introduced into cells and the association between each perturbation and single-cell gene expression is determined; (2) expression quantitative trait loci (eQTL) mapping, in which expression levels are associated with genotypes; and (3) peak-to-gene mapping, in which chromatin accessibility signals and gene expression measurements from the same cells (i.e., multiomic sequencing) are associated. In such examples, the scRNA-seq readout is associated with differentiating features associated with the cells, such as their perturbation, their genotype, or their accessibility signals. As discussed herein, various embodiments of the present technique comprise a computational pipeline that accurately maps such associations in a computationally rapid and efficient manner. In practice, such a computational pipeline may be implemented in or for a specialized computational device, e.g., a nucleic acid sequencing machine, such as a next generation sequencing (NGS) device. Such a sequencing device may be employed in generating and / or processing the data in accordance with the techniques described herein. As may be appreciated such sequencing devices and systems may be based on or otherwise employ specialized, non-conventional circuitry, including specialized processing circuitry, specialized memory circuitry or structures, application-specific integrated circuitry, specialized bus or communication structures, and so forth, that are beneficial in addressing the particular issues and problems arising from both the acquisition of and processing of vast amounts of nucleic acid sequence data, which is particularly true in the context of singlecell nucleic acid data as described herein. In certain such implementations one or more graphical processing units (GPUs) are employed to facilitate fast and efficient performance while maintaining accuracy.

[0009] The presently described technique may be particularly useful in the context of large data sets (e.g., data sets over 50,000 cells or 100,000 cells, with each cell constituting a single-cell data set). In particular, with the advent of next generation sequencing technology and corresponding developments in singe-cell isolation and techniques, large-scale scRNA-seq techniques have become viable for obtaining massive scRNA-seq data sets (e.g., data sets encompassing single-cell readout data for hundreds of thousands to millions of cells). Such large-scale scRNA-seq data sets cannot be processed for association mapping using conventionalapproaches based on correlation and / or model refitting for each feature set. By way of example, in Perturb-seq operations thousands of perturbations are associated with the expression of each gene. Such Perturb-seq datasets may scale up to millions of cells, which is an untenable extension of conventional techniques. However, the outputs of the techniques discussed herein may be understood to transform, using a generic or specialized computing environment, the distinctive input data (e.g., both single-cell expression and feature data) into valuable and easily accessible association insights between gene expression and feature permutations for vast numbers of individual cell in a sample from a subject. In this manner, the presently described techniques implement a specific, specialized software solution to address a specific problem (i.e., association of gene expression data and feature permutations for single cells in a sample comprising millions of cells or more) that may be implemented in a generic or specialized computing environment.

[0010] With this in mind, and as discussed in greater detail herein, the presently described techniques incorporate various features and benefits that facilitate accurate analysis at scale. First, in certain implementations the technique employs a model that accounts for over-dispersion and for structure present in the single-cell sequencing data, thereby enabling accurate mapping of association signals. Second, in many scRNA-seq association tasks, each gene is associated with thousands of predictors. As discussed herein, in certain implementations the analysis is streamlined so that a model for the expression of each gene is only fitted once, while each of the predictors is tested by estimating the effect of that predictor on the model without explicitly refitting the model (an analysis referred to as a “score test” herein in which a corresponding score based on the derivative of a likelihood function is determined). Third, to further facilitate large-scale analysis, in certain implementations the score test analysis is done in parallel (such as via multi-threaded processing) on CPUs or GPUs. In context where GPUs are employed, such a parallel implementation may leverage the ability of the GPUs to provide large-scale parallelized matrix multiplication to achieve improvement over traditional approaches, such as in terms of speed of implementation, which may be in a matter of hours for a sample comprising hundreds of thousands of cells.

[0011] In accordance with certain embodiments, a method is provided fortesting a plurality of features in single-cell data. In accordance with such a method, a model for gene expression is fit for single-cell data. The model corresponds to one model per gene. The effect of each feature onthe model is estimated without refitting the model.

[0012] In accordance with further embodiments, a computational pipeline for estimating perturbation effects in single-cell data is provided. In accordance with such an embodiment of a computational pipeline, single-cell RNA expression data is accessed or acquired. Perturbation data for each cell is determined. Preprocessing and clustering of the single-cell RNA expression data is performed. A model based on the perturbation and RNA expression data for each cell is calculated. The effect of each perturbation on the model is estimated based on respective derivatives calculated for each perturbation.

[0013] In accordance with additional embodiments, a nucleic acid sequencing system is provided. In accordance with such an embodiment, the nucleic acid sequencing system comprises one or more processors configured to execute code and process data and one or more memory devices configured to store processor executable code and nucleic acid sequence data. The processor executable code, when executed by the one or more processors, causes acts to be performed comprising: fitting a model for gene expression for single-cell data, wherein the model corresponds to one model per gene; and estimating the effect of each feature on the model without refitting the model.

[0014] Various refinements of the features noted above may exist in relation to various aspects of the present disclosure. Further features may also be incorporated in these various aspects as well. These refinements and additional features may exist individually or in any combination. For instance, various features discussed below in relation to one or more of the illustrated embodiments may be incorporated into any of the above-described aspects of the present disclosure alone or in any combination. The brief summary presented above is intended only to familiarize the reader with certain aspects and contexts of embodiments of the present disclosure without limitation to the claimed subject matter.BRIEF DESCRIPTION OF THE DRAWINGS

[0015] These and other features, aspects, and advantages of the present invention will become better understood when the following detailed description is read with reference to the accompanying drawings in which like characters represent like parts throughout the drawings,wherein:

[0016] FIG. 1 depicts a schematic view of an example of a system that may be used to provide biological or chemical analysis, in accordance with aspects of the present disclosure.

[0017] FIG. 2 depicts a cross-sectional view of an example of a flow cell that may be used in the system of FIG. 1, in accordance with aspects of the present disclosure.

[0018] FIG. 3 depicts a schematic view of another example of a system that may be used to provide biological or chemical analysis, in accordance with aspects of the present disclosure.

[0019] FIG. 4 depicts a schematic view of an example of a networked system in which the system of FIG. 1 or 3 may be incorporated, in accordance with aspects of the present disclosure.

[0020] FIG. 5 depicts a schematic view of an example of a base calling arrangement that may be carried out using the system of FIG. 1, 3, or 4, in accordance with aspects of the present disclosure.

[0021] FIG. 6 depicts a schematic view of an example of a base caller training technique that may be implemented using the system of FIG. 1, 3, or 4, in accordance with aspects of the present disclosure.

[0022] FIG. 7 depicts a feature matrix and an expression matrix for a respective cell, in accordance with aspects of the present disclosure.

[0023] FIG. 8 depicts a plot illustrating an example of over-dispersion in single-cell count data, in accordance with aspects of the present disclosure.

[0024] FIG. 9 depicts plots illustrating structure present within single-cell count data, in accordance with aspects of the present disclosure.

[0025] FIG. 10 depicts a plot of a likelihood function along with a corresponding derivative corresponding to a score test at the point corresponding to the null hypothesis, in accordance with aspects of the present disclosure.

[0026] FIGS. 11 A-l 1C depict plots exhibiting a strong perturbation effect in accordance with aspects of the present disclosure.

[0027] FIGS. 12A-12C depict plots exhibiting a weak perturbation effect in accordance with aspects of the present disclosure.

[0028] FIG. 13 depicts an example of a process flow for calculating a score test, in accordance with aspects of the present disclosure.

[0029] FIG. 14 depicts a plot illustrating improvement in the area under the receiver operating characteristic (AUROC) curve, in accordance with aspects of the present disclosure.

[0030] FIG. 15 depicts a plot illustrating improvement in the area under the precision recall curve (AUPRC), in accordance with aspects of the present disclosure.

[0031] FIG. 16 depicts the enrichment of STRING data plotted against the top N% of gene pairs with the highest correlation for the presently described technique and two reference techniques.

[0032] FIG. 17 depicts the enrichment of STRING data plotted against the top N gene pairs with the highest correlation for the presently described technique and two reference techniques.

[0033] FIG. 18 depicts enrichment of CORUM pair data plotted against the top 7V% of gene pairs with the highest correlation for the presently described technique and two reference techniques.

[0034] FIG. 19 depicts a plot of target coverage by cell line along with a table of median coverage by cell line.

[0035] FIG. 20 depicts the enrichment of STRING data plotted against the top N% of gene pairs with the highest correlation for various cell lines using the presently described technique.

[0036] FIG. 21 depicts the enrichment of STRING data plotted against the top Argene pairs with the highest correlation for various cell lines using the presently described technique.

[0037] FIG. 22 depicts the enrichment of CORUM pair data plotted against the top N% of gene pairs with the highest correlation for various cell lines using the presently described technique.

[0038] FIG. 23 depicts the enrichment of CORUM pair data plotted against the top N gene pairs with the highest correlation for various cell lines using the presently described technique.

[0039] FIG. 24 depicts CPU usage of the presently described techniques on 8.2 million cells.

[0040] FIG. 25 depicts memory utilization usage of the presently described techniques on 8.2 million cells.

[0041] FIG. 26 depicts job duration of the presently described techniques on 8.2 million cells.DETAILED DESCRIPTION

[0042] The following discussion is presented to enable any person skilled in the art to make and use the technology disclosed and is provided in the context of a particular application and its requirements. Various modifications to the disclosed implementations will be readily apparent to those skilled in the art, and the general principles defined herein may be applied to other implementations and applications without departing from the spirit and scope of the technology disclosed. Thus, the technology disclosed is not intended to be limited to the implementations shown, but is to be accorded the widest scope consistent with the principles and features disclosed herein.

[0043] The following detailed description of certain examples will be better understood when read in conjunction with the appended drawings. To the extent that the figures illustrate diagrams of the functional blocks of various examples, the functional blocks are not necessarily indicative of the division between hardware components. Thus, for example, one or more of the functional blocks (e.g., processors (central processing units (CPUs) or graphics processing units (GPUs)) or memories) may be implemented in a single piece of hardware (e.g., a general-purpose signal processor or random-access memory, hard disk, or the like). Similarly, the programs may be standalone programs, may be incorporated as subroutines in an operating system, may be functions in an installed software package, and the like. It should be understood that the various examples are not limited to the arrangements and instrumentality shown in the drawings.

[0044] I Overview of System for Biological or Chemical Analysis

[0045] Examples described herein may be used in various biological or chemical processes and systems for academic analysis, commercial analysis, or other analysis. More specifically, examples described herein may be used in various processes and systems where it is desired todetect an event, property, quality, or characteristic that is indicative of a designated reaction (e.g., methylation). Bioassay systems such as those described herein may be configured to perform a plurality of designated reactions that may be detected individually or collectively. For example, bioassay systems may be used to sequence a dense array of nucleic acid features through iterative cycles of enzymatic manipulation and image acquisition. In some examples, nucleic acids can be attached to a surface and amplified. Examples of such amplification are described in U.S. Pat. No.7,741,463, entitled “Method of Preparing Libraries of Template Polynucleotides,” issued June 22, 2010, the disclosure of which is incorporated by reference herein, in its entirety; and / or U.S. Pat. No. 7,270,981, entitled “Recombinase Polymerase Amplification,” issued September 18, 2007, the disclosure of which is incorporated by reference herein, in its entirety.

[0046] Components that are used in the bioassay systems may include one or more microfluidic channels that deliver reagents or other reaction components to a reaction site. The reaction sites may be randomly distributed across a substantially planar surface; or may be patterned across a substantially planar surface. Each of the reaction sites may be imaged to detect light from the reaction site. The signals indicating photons emitted from the reaction sites and detected by image sensors may provide illumination values. These illumination values may be combined into an image indicating photons as detected from the reaction sites. These images may be further analyzed to identify compositions, reactions, conditions, etc., at each reaction site.

[0047] II. Examples of Fluidics Devices and Fluid Flow Paths - Example of System with Higher Volume Throughput

[0048] FIG. 1 illustrates a schematic diagram of an example of a system (100) that may be used to perform an analysis on one or more samples of interest. In some implementations, the sample may include one or more clusters of nucleotides (e.g., DNA) that have been linearized to form a single stranded DNA (sstDNA). In the implementation shown, system (100) is configured to receive a flow cell cartridge assembly (102) including a flow cell assembly (103) and a sample cartridge (104). System (100) includes a flow cell receptacle (122) that receives flow cell cartridge assembly (102), a vacuum chuck (124) that supports flow cell assembly (103), and a flow cell interface (126) that is used to establish a fluidic coupling between system (100) and flow cell assembly (103). Flow cell interface (126) may include one or more manifolds. System (100) further includes a sipper manifold assembly (106), a sample loading manifold assembly (108), anda pump manifold assembly (110). System (100) also includes a drive assembly (112), a controller (114), an imaging system (116), and a waste reservoir (118). Controller (114) is electrically and / or communicatively coupled to drive assembly (112) and to imaging system (116); and is configured to cause drive assembly (112) and / or the imaging system (116) to perform various functions as disclosed herein.

[0049] In the present example, flow cell assembly (103) includes a flow cell (128) having a channel (130) and defining a plurality of first openings (132), which are fluidically coupled to the channel (130) and arranged on a first side (134) of the channel (130). Flow cell (128) further includes a plurality of second openings (136) fluidically coupled to the channel (130) and arranged on a second side (138) of the channel (130). Fluid may thus flow through flow cell (128) via channel. While the flow cell (128) is shown including one channel (130), flow cell (128) may include two or more channels (130). Flow cell assembly (103) also includes a flow cell manifold assembly (140) coupled to flow cell (128) and having a first manifold fluidic line (142) and a second manifold fluidic line (144). Flow cell manifold assembly (140) may be in the form of a laminate including a plurality of layers as discussed in more detail below.

[0050] In the implementation shown, first manifold fluidic line (142) has a first fluidic line opening (146) and is fluidically coupled to each of the first openings (132) of flow cell (128); and second manifold fluidic line (144) has a second fluidic line opening (148) and is fluidically coupled to each of the second openings (136). As shown, flow cell assembly (103) includes gaskets (150) coupled to flow cell manifold assembly (140) and fluidically coupled to fluidic line openings (146, 148). In some implementations where flow cell (128) includes a plurality of channels (130), flow cell manifold assembly (140) may include additional fluidic lines (152) that couple first fluidic line openings (146) to a single manifold port (154). In such implementations, a single gasket (150) may be coupled to flow cell manifold assembly (140) that surrounds the manifold port (154) and is in fluidic communication with a plurality of channels (130). In operation, flow cell interface (126) engages with corresponding gaskets (150) to establish a fluidic coupling between system (100) and flow cell (128). The engagement between flow cell interface (126) and gaskets (150) reduces or eliminates fluid leakage between flow cell interface (126) and flow cell (128).

[0051] In the implementation shown, first manifold fluidic line (142) has a portion (156) that is substantially parallel to a longitudinal axis (158) of channel (130); and second manifold fluidicline (144) has a portion (160) that is substantially parallel to longitudinal axis (158) of channel (130). Additionally, first manifold fluidic line (142) is shown being at least partially adjacent a first end (162) of flow cell (128) and spaced from a second end (164) of flow cell (128); and second manifold fluidic line (144) is shown being at least partially adjacent second end (164) of flow cell (128) and spaced from first end (162). Other arrangements of manifold fluidic lines (142, 144) may prove suitable, however.

[0052] In the implementation shown, system (100) includes a sample cartridge receptacle (166) that receives sample cartridge (104) that carries one or more samples of interest (e g., an analyte). System (100) also includes a sample cartridge interface (168) that establishes a fluidic connection with sample cartridge (104). Sample loading manifold assembly (108) includes one or more sample valves (170). Pump manifold assembly (110) includes one or more pumps (172), one or more pump valves (174), and a cache (176). Valves (170, 174) and pumps (172) may take any suitable form. Cache (176) may include a serpentine cache and may temporarily store one or more reaction components during, for example, bypass manipulations of the system (100). While cache (176) is shown being included in pump manifold assembly (110), cache (176) may alternatively be located elsewhere (e.g., in sipper manifold assembly (106) or in another manifold downstream of a bypass fluidic line (178), etc.).

[0053] Sample loading manifold assembly (108) and pump manifold assembly (110) flow one or more samples of interest from sample cartridge (104) through a fluidic line (180) toward flow cell cartridge assembly (102). In some implementations, sample loading manifold assembly (108) may individually load or address each channel (130) of flow cell (128) with a respective sample of interest. The process of loading channel (130) with a sample of interest may occur automatically using system (100). As shown in FIG. 1, sample cartridge (104) and sample loading manifold assembly (108) are positioned downstream of flow cell cartridge assembly (102). In the implementation shown, sample loading manifold assembly (108) is coupled between flow cell cartridge assembly (102) and pump manifold assembly (110). To draw a sample of interest from sample cartridge (104) and toward pump manifold assembly (110), sample valves (170), pump valves (174), and / or pumps (172) may be selectively actuated to urge the sample of interest toward pump manifold assembly (110). Sample cartridge (104) may include a plurality of sample reservoirs that are selectively fluidically accessible via the corresponding sample valves (170). Toindividually flow the sample of interest toward channel (130) of flow cell (128) and away from pump manifold assembly (110), sample valves (170), pump valves (174), and / or pumps (172) may be selectively actuated to urge the sample of interest toward flow cell cartridge assembly (102) and into respective channels (130) of flow cell (128).

[0054] Drive assembly (112) interfaces with sipper manifold assembly (106) and pump manifold assembly (110) to flow one or more reagents that interact with the sample within flow cell (128). In some scenarios, a reversible terminator is attached to the reagent to allow a single nucleotide to be incorporated onto a growing DNA strand. In some such implementations, one or more of the nucleotides has a unique fluorescent label that emits a color when excited. The color (or absence thereof) is used to detect the corresponding nucleotide. In the implementation shown, imaging system (116) excites one or more of the identifiable labels (e.g., a fluorescent label) and thereafter obtains image data for the identifiable labels. The labels may be excited by incident light and / or a laser and the image data may include one or more colors emitted by the respective labels in response to the excitation. The image data (e.g., detection data) may be analyzed by system (100). Examples of features and functionalities that may be incorporated into imaging system (116) will be described in greater detail below.

[0055] After the image data is obtained, drive assembly (112) interfaces with sipper manifold assembly (106) and pump manifold assembly (110) to flow another reaction component (e.g., a reagent) through flow cell (128) that is thereafter received by waste reservoir (118) via a primary waste fluidic line (182) and / or otherwise exhausted by system (100). Some reaction components may perform a flushing operation that chemically cleaves the fluorescent label and the reversible terminator from the sstDNA. The sstDNA may then be ready for another cycle.

[0056] The primary waste fluidic line (182) is coupled between pump manifold assembly (110) and waste reservoir (118). In some implementations, pumps (172) and / or pump valves (174) of pump manifold assembly (110) selectively flow the reaction components from flow cell cartridge assembly (102), through fluidic line (180) and sample loading manifold assembly (108) to primary waste fluidic line (182). Flow cell cartridge assembly (102) is coupled to a central valve (184) via flow cell interface (126). Central valve (184) is coupled with flow cell interface (126) via a fluidic line (185). An auxiliary waste fluidic line (186) is coupled to central valve (184) and to waste reservoir (118). In some implementations, auxiliary waste fluidic line (186) receives excess fluidof a sample of interest from flow cell cartridge assembly (102), via central valve (184), and flows the excess fluid of the sample of interest to waste reservoir (118) when back loading the sample of interest into flow cell (128), as described herein.

[0057] Sipper manifold assembly (106) includes a shared line valve (188) and a bypass valve (190). Shared line valve (188) may be referred to as a reagent selector valve. Central valve (184) and the valves (188, 190) of sipper manifold assembly (106) may be selectively actuated to control the flow of fluid through fluidic lines (192, 194, 196). Sipper manifold assembly (106) may be coupled to a corresponding number of reagent reservoirs (198) via reagent sippers (200). Reagent reservoirs (198) may contain fluid (e g., reagent and / or another reaction component). In some implementations, sipper manifold assembly (106) includes a plurality of ports. Each port of sipper manifold assembly (106) may receive one of the reagent sippers (200). Reagent sippers (200) may be referred to as fluidic lines. Some forms of reagent sippers (200) may include an array of sipper tubes extending downwardly along the z-dimension from ports in the body of sipper manifold assembly (106). Reagent reservoirs (198) may be provided in a cartridge, and the tubes of reagent sippers (200) may be configured to be inserted into corresponding reagent reservoirs (198) in the reagent cartridge so that liquid reagent may be drawn from each reagent reservoir (198) into the sipper manifold assembly (106).

[0058] Shared line valve (188) of sipper manifold assembly (106) is coupled to central valve (184) via shared reagent fluidic line (196). Different reagents may flow through shared reagent fluidic line (196) at different times. In some versions, when performing a flushing operation before changing between one reagent and another, pump manifold assembly (110) may draw wash buffer through shared reagent fluidic line (196), central valve (184), and flow cell cartridge assembly (102).

[0059] Bypass valve (190) of sipper manifold assembly (106) is coupled to central valve (184) via dedicated reagent fluidic lines (194, 196). Each of the dedicated reagent fluidic lines (194, 196) may be associated with a single reagent. The fluids that may flow through dedicated reagent fluidic lines (194, 196) may be used during sequencing operations and may include a cleave reagent, an incorporation reagent, a scan reagent, a cleave wash, and / or a wash buffer.

[0060] Bypass valve (190) is also coupled to cache (176) of pump manifold assembly (110)via bypass fluidic line (178). One or more reagent priming operations, hydration operations, mixing operations, and / or transfer operations may be performed using bypass fluidic line (178). The priming operations, the hydration operations, the mixing operations, and / or the transfer operations may be performed independent of flow cell cartridge assembly (102). Thus, the operations using bypass fluidic line (178) may occur during, for example, incubation of one or more samples of interest within flow cell cartridge assembly (102). That is, shared line valve (188) may be utilized independently of bypass valve (190) such that bypass valve (190) may utilize bypass fluidic line (178) and / or cache (176) to perform one or more operations while shared line valve (188) and / or central valve (184) simultaneously, substantially simultaneously, or offset synchronously perform other operations.

[0061] Drive assembly (112) includes a pump drive assembly (202) and a valve drive assembly (204). Pump drive assembly (202) may be adapted to interface with one or more pumps (172) to pump fluid through flow cell (128) and / or to load one or more samples of interest into flow cell (128). Valve drive assembly (204) may be adapted to interface with one or more of the valves (170, 174, 184, 188, 190) to control the position of the corresponding valves (170, 174, 184, 188, 190).

[0062] Controller (114) of the present example includes a user interface (206), a communication interface (208), one or more processors (210) (e.g., central processing units (CPUs) or graphics processing units (GPUs)), and a memory (212) storing instructions executable by the one or more processors (210) to perform various functions including the disclosed implementations. User interface (206), communication interface (133), and memory (212) are electrically and / or communicatively coupled to the one or more processors (210). User interface (206) may be adapted to receive input from a user and to provide information to the user associated with the operation of system (100) and / or an analysis taking place. User interface (206) may include a touch screen, a display, a keyboard, a speaker(s), a mouse, a track ball, and / or a voice recognition system.

[0063] Communication interface (208) is adapted to enable communication between system (100) and a remote system(s) (e.g., computers) via a network(s) (e.g., the Internet, an intranet, a local-area network (LAN), a wide-area network (WAN), a coaxial-cable network, a wireless network, a wired network, a satellite network, a digital subscriber line (DSL) network, a cellularnetwork, a Bluetooth connection, a near field communication (NFC) connection, etc.). Some of the communications provided to the remote system may be associated with analysis results, imaging data, etc. generated or otherwise obtained by system (100). Some of the communications provided to system (100) may be associated with a fluidics analysis operation, patient records, and / or a protocol(s) to be executed by system (100).

[0064] The one or more processors (210) and / or system (100) may include one or more of a processor-based system(s) or a microprocessor-based system(s). In some implementations, the one or more processors (210) and / or system (100) includes one or more of a programmable processor, a programmable controller, a microprocessor, a microcontroller, a graphics processing unit (GPU), a digital signal processor (DSP), a reduced-instruction set computer (RISC), an application specific integrated circuit (ASIC), a field programmable gate array (FPGA), a field programmable logic device (FPLD), a logic circuit, and / or another logic-based device executing various functions including the ones described herein.

[0065] Memory (212) may include one or more of a semiconductor memory, a magnetically readable memory, an optical memory, a hard disk drive (HDD), an optical storage drive, a solid-state storage device, a solid-state drive (SSD), a flash memory, a read-only memory (ROM), erasable programmable read-only memory (EPROM), electrically erasable programmable readonly memory (EEPROM), a random-access memory (RAM), a non-volatile RAM (NVRAM) memory, a compact disc (CD), a compact disc read-only memory (CD-ROM), a digital versatile disk (DVD), a Blu-ray disk, a redundant array of independent disks (RAID) system, a cache and / or any other storage device or storage disk in which information is stored for any duration (e.g., permanently, temporarily, for extended periods of time, for buffering, for caching).

[0066] III. Examples of Flow Cell Structures

[0067] As noted above, a system (100) may execute reactions in a flow cell (128) and / or perform analysis on one or more samples of interest in a flow cell (128). The following describes examples of forms that such flow cells (128) may take, it being understood that flow cells (128) may take various other forms and have various other features in addition to or in lieu of the features described below.

[0068] Example of Single-Surface Patterned Flow Cell

[0069] FIG. 2 shows an example of a flow cell (400) that includes a patterned substrate (402), which includes depressions (404) separated by interstitial regions (406), and surface chemistry (410, 412) positioned in the depressions (404). Depressions (404) may be in the form of microwells or nanowells. Depressions (404) may be configured to contain nucleic acid strands or other oligonucleotides and thereby provide a reaction site for SBS and / or for other kinds of processes. In some versions, each depression (404) has a cylindraceous configuration, with a generally circular cross-sectional profile. In some other versions, each depression (404) has a polygonal (e.g., hexagonal, octagonal, square, rectangular, elliptical, etc.) cross-sectional profile. Alternatively, depressions (404) may have any other suitable configuration. It should also be understood that depressions (404) may be arranged in any suitable pattern, including but not limited to a grid pattern.

[0070] Surface chemistry (410, 412) of the present example includes functionalized coating layer (410) and primers (412). While not shown, it is to be understood that the depressions (404) may also have surface preparation or treatment chemistry (e.g., silane or a silane derivative) positioned between the substrate (402) and the functionalized coating layer (410). This same surface preparation or treatment chemistry may also be positioned on the interstitial regions (406). In the present example, a hydrogel (440) is applied before lid (420) is bonded to substrate (402). Hydrogel (440) covers surface chemistry (410, 412) in depressions (404), and at least a portion of the patterned substrate (402) (e.g., those interstitial regions (406) that are not also bonding regions (422)). By way of example only, hydrogel (440) may comprise PAZAM, crosslinked polyacrylamide, agarose gel, etc.

[0071] Flow cell (400) of this example further includes a lid (420) bonded to bonding region(s) (422) of patterned substrate (402). In the example shown in FIG. 2, lid (420) includes a top portion (424) that is connected to several sidewalls (426), and these components (424, 426) define a portion of each of the six flow channels (430A, 430B, 430C, 430D, 430E, 430F). The respective sidewalls (426) isolate one flow channel (430A, 430B, 430C, 430D, 430E, 430F) from each adjacent flow channel (430A, 430B, 430C, 430D, 430E, 430F). Each flow channel (430A, 430B, 430C, 430D, 430E, 430F) is in selective fluid communication with a respective set of depressions (404).

[0072] Lid (420) may be bonded to bonding region (422) of substrate (402) using any suitabletechnique, such as laser bonding, diffusion bonding, anodic bonding, eutectic bonding, plasma activation bonding, glass frit bonding, or other methods known in the art. In some versions, a spacer layer (428) may be used to bond lid (420) to bonding region (422). Spacer layer (428) may comprise any material that will seal at least some of interstitial regions (404) (e.g., bonding region (422)) of substrate (402) and lid (420) together. While not shown, lid (420) or the patterned substrate (402) may include inlet and outlet ports that are to fluidically engage other ports (not shown), such as those of sample cartridge interface (168), for directing fluid(s) into the respective flow channels (430A, 430B, 430C, 430D, 430E, 430F) (e g., from a reagent cartridge or other fluid storage system) and out of the flow channel (e.g., to waste reservoir (118) or another waste removal system). Flow channels (430A, 430B, 430C, 430D, 430E, 430F) may serve to, for example, selectively introduce reaction components or reactants to hydrogel (440) and the underlying surface chemistry (410, 412) in order initiate designated reactions in / at depressions (404).

[0073] While flow cell (400) includes a pattern of depressions (404) to provide an array of reaction sites, other variations may provide reaction sites on or at various other kinds of structural features, including but not limited to continuously planer surfaces and / or protruding surfaces, etc. By way of further example only, flow cell (400) may be constructed and operable in accordance with at least some of the teachings of U.S. Pat. No. 10,919,033, entitled “Flow Cells with Hydrogel Coating,” issued February 16, 2021, the disclosure of which is incorporated by reference herein, in its entirety.

[0074] IV. Examples of Imaging System Features

[0075] As noted above, system (100) includes an imaging system (116) that excites one or more identifiable labels (e.g., a fluorescent label) in samples in reaction sites provided by depressions (404, 462, 464) of a flow cell (128, 400); and thereafter obtains image data for the identifiable labels. This image data is used to identify nucleotides as part of a nucleic acid sequencing process. Alternatively, the image data may be used for various other purposes. The following description provides details on how some versions of imaging system (116) may be configured and operable.

[0076] FIG. 3 illustrates a schematic diagram of another example of a system (500) that maybe used to perform an analysis on one or more samples of interest. Except as otherwise described below, system (500) of this example may be configured and operable like systems (100) described above. System (500) is configured to perform a large number of parallel reactions within a flow cell (510). Flow cell (510) may be configured and operable like flow cells (400) described above or may have any other suitable configuration. Flow cell (510) may thus include one or more flow channels that receive a solution from system (500) and direct the solution toward reaction sites of flow cell (510).

[0077] System (500) includes a system controller (520) that may communicate with the various components, assemblies, and sub-systems of the system (500). Controller (520) may be configured and operable like controllers (114) described above. An imaging assembly (522) of system (500) includes a light emitting assembly (550) that emits light that reaches reaction sites on flow cell (510). Light emitting assembly (550) may include an incoherent light emitter (e.g., emit light beams output by one or more excitation diodes), or a coherent light emitter such as emitter of light output by one or more lasers or laser diodes. In some implementations, light emitting assembly (550) may include a plurality of different light sources (not shown), each light source emitting light of a different wavelength range. Some versions of light emitting assembly (550) may also include one or more collimating lenses (not shown), a light structuring optical assembly (not shown), a projection lens (not shown) that is operable to adjust a structured beam shape and path, epifluorescence microscopy components, and / or other components. Although system (500) is illustrated as having a single light emitting assembly (550), multiple light emitting assemblies (550) may be included in some other implementations.

[0078] In the present example, the light from light emitting assembly (550) is directed by dichroic mirror assembly (546) through an objective lens assembly (542) onto a sample of a flow cell (510), which is positioned on a motion stage (570). In the case of fluorescent microscopy of a sample, a fluorescent element associated with the sample of interest fluoresces in response to the excitation light, and the resultant light is collected by objective lens assembly (542) and is directed to an image sensor of camera system (540) to detect the emitted fluorescence. In some implementations, a tube lens assembly may be positioned between the objective lens assembly (542) and the dichroic mirror assembly (546) or between the dichroic mirror (546) and the image sensor of the camera system (540). A moveable lens element may be translatable along alongitudinal axis of the tube lens assembly to account for focusing on an upper interior surface or lower interior surface of the flow cell (510) and / or spherical aberration introduced by movement of the objective lens assembly (542).

[0079] In the present example, a filter switching assembly (544) is interposed between dichroic mirror assembly (546) and camera system (540). Filter switching assembly (544) includes one or more emission filters that may be used to pass through particular ranges of emission wavelengths and block (or reflect) other ranges of emission wavelengths. For example, emission filters may be used to direct different wavelength ranges of emitted light to different image sensors of the camera system (540) of imaging assembly (522). For instance, the emission filters may be implemented as dichroic mirrors that direct emission light of different wavelengths from flow cell (510) to different image sensors of camera system (540). In some variations, a projection lens is interposed between filter switching assembly (544) and camera system (540). Filter switching assembly (544) may be omitted in some versions.

[0080] System (500) further includes a fluid delivery assembly (590) that may direct the flow of reagents (e.g., fluorescently labeled nucleotides, buffers, enzymes, cleavage reagents, etc.) to (and through) flow cell (510) and waste valve (580). Fluid delivery assembly (590) may be configured and operable like the various fluid delivery components described herein. System (500) of the present example also includes a temperature station actuator (530) and heater / cooler (532) that may optionally regulate the temperature of conditions of the fluids within the flow cell (510). In some implementations, the heater / cooler (532) may be fixed to sample stage (570), upon which the flow cell (510) is placed, and / or may be integrated into sample stage (570).

[0081] Flow cell (510) may be removably mounted on sample stage (570), which may provide movement and alignment of flow cell (510) relative to objective lens assembly (542). Sample stage (570) may have one or more actuators to allow sample stage (570) to move in any of three dimensions. For example, actuators may be provided to allow sample stage (570) to move in the x, y, and z directions relative to objective lens assembly (542), tilt relative to objective lens assembly (542), and / or otherwise move relative to objective lens assembly (542). Movement of sample stage (570) may allow one or more sample locations on flow cell (510) to be positioned in optical alignment with objective lens assembly (542). Movement of sample stage (570) relative to objective lens assembly (542) may be achieved by moving sample stage (570) itself, by movingobjective lens assembly (542), by moving some other component of imaging assembly (522), by moving some other component of system (500), or any combination of the foregoing. For instance, in some implementations, the sample stage (570) may be actuatable in the x and y directions relative to the objective lens assembly (542) while a focus component (562) or z-stage may move the objective lens assembly (542) along the z direction relative to the sample stage (570).

[0082] In some implementations, a focus component (562) may be included to control positioning of one or more elements of objective lens assembly (542) relative to the flow cell (510) in the focus direction (e.g., along the z-axis or z-dimension). Focus component (562) may include one or more actuators physically coupled to the objective lens assembly (542), the optical stage, the sample stage (570), or a combination thereof, to move flow cell (510) on sample stage (570) relative to the objective lens assembly (542) to provide proper focusing for the imaging operation. In the present example, the focus component (562) utilizes a focus tracking module (560) that is configured to detect a displacement of the objective lens assembly (542) relative to a portion of the flow cell (510) and output data indicative of an in-focus position to the focus component (562) or a component thereof or operable to control the focus component (562), such as controller (520), to move the objective lens assembly (542) to position the corresponding portion of the flow cell (510) in focus of the objective lens assembly (542).

[0083] In some implementations, an actuator of focus component (562) or for sample stage (570) may be physically coupled to objective lens assembly (542), the optical stage, sample stage (570), or a combination thereof, such as, for example, by mechanical, magnetic, fluidic, or other attachment or contact directly or indirectly to or with the stage or a component thereof. The actuator of focus component (562) may be configured to move objective lens assembly (542) in the z-direction while maintaining sample stage (570) in the same plane (e.g., maintaining a level or horizontal attitude, perpendicular to the optical axis). In some implementations, sample stage (570) includes an x direction actuator and a y direction actuator to form an x-y stage. Sample stage (570) may also be configured to include one or more tip or tilt actuators to tip or tilt sample stage (570) and / or a portion thereof, to account for any slope in its surfaces.

[0084] Camera system (540) may include one or more image sensors to monitor and track the imaging (e.g., sequencing) of flow cell (510). Camera system (540) may be implemented, for example, as a CCD or CMOS image sensor camera, but other image sensor technologies (e.g.,active pixel sensor) may be used. By way of further example only, camera system (540) may include a dual-sensor time delay integration (TDI) camera, a single-sensor camera, a camera with one or more two-dimensional image sensors, and / or other kinds of camera technologies. While camera system (540) and associated optical components are shown as being positioned above flow cell (510) in FIG. 3, one or more image sensors or other camera components may be incorporated into system (500) in numerous other ways as will be apparent to those skilled in the art in view of the teachings herein. For instance, one or more image sensors may be positioned under flow cell (510), such as within the sample stage (570) or below the sample stage (570); or may even be integrated into flow cell (510).

[0085] V. Examples of Data Processing Features

[0086] A. Example of Networked Data Processing Arrangement

[0087] As noted above, a system (100, 500) may include a controller (114, 520) that is configured to process data, execute algorithms, etc., as needed to perform a sequencing operation or other kind of operation. In some scenarios, system (100, 500) may be coupled with other devices via a network to perform further data processing, data storage, execution of algorithms, etc. FIG.4 shows an example of such an arrangement. In particular, FIG. 4 shows a networked system (800) that includes a sequencing device (810), a server device (820), a client device (830), and a local device (840), with all devices (810, 820, 830, 840) being coupled together via a network (850). Network (850) may take any suitable form as will be apparent to those skilled in the art in view of the teachings herein.

[0088] As shown in FIG. 4, sequencing device (810) comprises a computing device and a sequencing device system (812) for sequencing a genomic sample or other nucleic-acid polymer. In some versions, by executing sequencing device system (812) using a processor, sequencing device (810) analyzes nucleotide fragments or oligonucleotides extracted from genomic samples to generate nucleotide reads or other data utilizing computer implemented methods and systems either directly or indirectly on sequencing device (810). More particularly, sequencing device (810) receives nucleotide- sample slides (e.g., flow cells (128, 400, 510)) comprising nucleotide fragments extracted from samples and further copies and determines the nucleobase sequence of such extracted nucleotide fragments. It should be understood that sequencing device (810) mayrepresent a version of systems (100, 500) described above.

[0089] In some versions, the sequencing device (810) utilizes SBS to sequence nucleotide fragments into nucleotide reads and determine nucleobase calls for the nucleotide reads. By executing sequencing device system (812), sequencing device (810) may further store the nucleobase calls as part of base-call data that is formatted as a binary base call (BCL) fde and send the BCL fde to the local device (840) and / or the server device(s) (820). Sequencing device (810) may communicate the BCL fde and / or other data to local device (840) and / or client device (830) via network (850) or directly (i.e., bypassing network (850)).

[0090] In some scenarios, local device (840) is located at or near a same physical location of sequencing device (810). For instance, local device (840) and sequencing device (810) may be integrated into a single computing device. Local device (840) may run sequencing system (814) to generate, receive, analyze, store, and transmit digital data, such as by receiving base-call data or determining variant calls based on analyzing such base-call data. By executing software in the form of sequencing system (814), local device (840) may align nucleotide reads with a structural variation graph genome (824) and determine genetic variants based on the aligned nucleotide reads. Local device (840) may also send data to client device (830), including a variant call file (VCF) or other information indicating nucleobase calls, sequencing metrics, error data, or other metrics.

[0091] Server device(s) (820) may be located remotely from the local device (840) and sequencing device (810). Server device(s) (820) may comprise a distributed collection of servers, where server device(s) (820) include a number of server devices distributed across network (850) and located in the same or different physical locations. Similar to local device (840), server device(s) (820) may include a version of sequencing system (814). Accordingly, server device(s) (820) may generate, receive, analyze, store, and transmit digital data, such as by receiving basecall data or determining variant calls based on analyzing such base-call data. As indicated above, sequencing device (810) may send (and server device(s) (820) may receive) base-call data from sequencing device (810). Server device(s) (820) may also send data to client device (830), including VCFs or other sequencing related information.

[0092] As indicated above, as part of server device(s) (820) or local device (840), sequencingsystem (814) may generate or implement a structural variation graph genome with alternate contiguous sequences representing structural variant haplotypes. For instance, system (814) may identify candidate structural variants of a threshold frequency (or that otherwise satisfy another occurrence threshold) within a genomic sample database. From among the candidate structural variants, sequencing system (814) selects structural variant haplotypes based on one or both of satisfying another occurrence threshold and finding flanking variants adjacent to particular structural variant haplotypes. Sequencing system (814) may likewise select reference haplotypes of genomic regions corresponding to the selected structural variant haplotypes from a reference genome. Based on the selected haplotypes, sequencing system (814) generates a structural variation graph genome comprising both alternate contiguous sequences representing the structural variant haplotypes and reference sequences representing the reference haplotypes. Based on comparing nucleotide reads of a genomic sample with alternate contiguous sequences representing structural variant haplotypes, sequencing system (814) can determine nucleobase calls for the genomic sample.

[0093] By executing a sequencing application (832), client device (830) may generate, store, receive, and send digital data. In particular, client device (830) may receive sequencing data from local device (840) or receive call files (e.g., BCL) and sequencing metrics from sequencing device (810). Furthermore, client device (830) may communicate with local device (840) or server device(s) (820) to receive a VCF comprising nucleobase calls and / or other metrics, such as a basecall-quality metrics or pass-filter metrics. Client device (830) may accordingly present or display information pertaining to variant calls or other nucleobase calls within a graphical user interface of sequencing application (832) to a user associated with client device (830). For example, client device (830) may present structural variant calls and / or sequencing metrics for a sequenced genomic sample within a graphical user interface of sequencing application (832).

[0094] As shown in FIG. 4, sequencing application (832) is included in client device (830). Sequencing application (832) may include a web application or a native application stored and executed on client device (830) (e.g., a mobile application, desktop application). Sequencing application (832) may include instructions that (when executed) cause client device (830) to receive data from sequencing system (814) and present, for display at client device (830), basecall data or data from a VCF. Furthermore, sequencing application (832) may instruct client device(830) to display summaries for multiple sequencing runs.

[0095] As further illustrated in FIG. 4, a version of sequencing system (814) may be located and implemented (e.g., entirely or in part) on client device (830) or sequencing device (810). In some versions, sequencing system (814) is implemented by one or more other components of networked system (800), such as local device (840). In particular, sequencing system (814) may be implemented in a variety of different ways across sequencing device (810), local device (840), server device(s) (820), and client device (830). For example, sequencing system (814) may be downloaded from server device(s) (820) to sequencing system (814) and / or local device (840) where all or part of the functionality of sequencing system (814) is performed at each respective device within networked system (800).

[0096] B. Examples of Base Calling Schemes

[0097] FIG. 5 illustrates a system (900) that employs two or more base callers for base calling operations on the raw images (i.e., sensor data) output by image sensors in a sequencing machine sequencing machine (910). Sequencing machine (910) of this example includes a flow cell (912), which includes a plurality of tiles (914). Each tile (914) includes a plurality of clusters (916). Sequencing machine (910) may be understood to represent a version of systems (100, 500) or sequencing device (810) described above; while flow cell (912) may be understood to represent a version of flow cells (128, 400, 510) described above. Sequencing machine (910) thus outputs sensor data (920) comprising raw images from the tiles (914) of flow cell (912).

[0098] In the present example, system (900) comprises a first base caller (922) and a second base caller (926), though some variations may include more than two base callers (922, 926). Each base caller (922, 926) of this example outputs corresponding base call classification information. For example, first base caller (922) outputs first base call classification information (924); and second base caller (926) outputs second base call classification information (928). A base calling combining module (930) generates final base calls (932), based on one or both first base call classification information (924) and / or second base call classification information (928). In some versions, first base caller (922) is a neural -network based base-caller; while second base caller (926) is a non-neural network based base-caller. For example, first base caller (922) may include a non-linear system employing one or more neural network models for base calling. The first basecaller (922) may also be referred to as a DeepRTA (Deep Real Time Analysis) base caller or Deep Neural Network base caller.

[0099] By way of further example only, second base caller (926) may include, at least in part, a linear system used for base calling. For example, some versions of second base caller (926) do not employ a neural network for base calling (or use a smaller neural network model for base calling, compared to a larger neural network model used by first base caller (922)). Second base caller (926) may also be referred to as an RTA (Real Time Analysis) base caller. An RTA base caller may use linear intensity extractors to extract features from sequencing images for base calling. In some such versions, RTA performs a template generation step to produce a template image that identifies locations of clusters (916) on a tile (914) using sequencing images from some number of initial sequencing cycles called template cycles. The template image is used as a reference for subsequent registration and intensity extraction steps. The template image is generated by detecting and merging bright spots in each sequencing image of the template cycles, which in turn involves sharpening a sequencing image (e.g., using the Laplacian convolution), determining an “on” threshold by a spatially segregated Otsu approach, and subsequent five-pixel local maximum detection with subpixel location interpolation.

[0100] In another example, locations of clusters (916) on a tile (914) are identified using fiducial markers. A solid support upon which a biological specimen is imaged may include such fiducial markers, to facilitate determination of the orientation of the specimen or the image thereof in relation to probes that are attached to the solid support. Examples of fiducials include, but are not limited to, beads (with or without fluorescent moieties or moieties such as nucleic acids to which labeled probes can be bound), fluorescent molecules attached at known or determinable features, or structures that combine morphological shapes with fluorescent moieties.

[0101] RTA then registers a current sequencing image against the template image. This is achieved by using image correlation to align the current sequencing image to the template image on a sub-region, or by using non-linear transformations (e.g., a full six-parameter linear affine transformation). RTA generates a color matrix to correct cross-talk between color channels of the sequencing images. RTA implements empirical phasing correction to compensate noise in the sequencing images caused by phase errors. After different corrections are applied to the sequencing images, RTA extracts signal intensities for each spot location in the sequencingimages. For example, for a given spot location, signal intensity may be extracted by determining a weighted average of the intensity of the pixels in a spot location. For example, a weighted average of the center pixel and neighboring pixels may be performed using bilinear or bicubic interpolation. In some implementations, each spot location in the image may comprise a few pixels (e.g., 1-5 pixels). RTA then spatially normalizes the extracted signal intensities to account for variation in illumination across the sampled imaged. For example, intensity values may be normalized such that a 5th and 95th percentiles have values of 0 and 1, respectively. The normalized signal intensities for the image (e.g., normalized intensities for each channel) may be used to calculate mean chastity for the plurality of spots in the image.

[0102] In some implementations, RTA uses an equalizer to maximize the signal-to-noise ratio of the extracted signal intensities. The equalizer may be trained (e.g., using least square estimation, adaptive equalization algorithm) to maximize the signal-to-noise ratio of cluster intensity data in sequencing images. In some implementations, the equalizer is a lookup table (LUT) bank with a plurality of LUTs with subpixel resolution, also referred to as “equalizer filters” or “convolution kernels.” By way of example only, the number of LUTs in the equalizer may depend on the number of subpixels into which pixels of the sequencing images can be divided. For example, if the pixels are divisible into w-by-n subpixels (e.g., 5 x 5 subpixels), then the equalizer generates ir LUTs (e.g., 25 LUTs).

[0103] In some implementations of training the equalizer, data from the sequencing images is binned by well subpixel location. It should be understood that a “well” may include depressions (404) of a flow cell (400) or any other kind of reaction site (e g., in a flow cell or otherwise). In an example of sequencing images being binned by well subpixel location, for a 5 x 5 LUT, l / 25th of the wells have a center that is in bin (1,1) (e.g., the upper left comer of a sensor pixel), l / 25th of the wells are in bin (1,2), and so on. The equalizer coefficients for each bin may be determined using least squares estimation on the subset of data from the wells corresponding to the respective bins. This way, the resulting estimated equalizer coefficients are different for each bin. Each LUT / equalizer filter / convolution kernel has a plurality of coefficients that are learned from the training. The number of coefficients in a LUT may correspond to the number of pixels that are used for base calling a cluster. For example, if a local grid of pixels (image or pixel patch) that is used to base call a cluster is of size p xp (e.g., 9x 9 pixel patch), then each LUT has p1coefficients(e g., 81 coefficients). The training may produce equalizer coefficients that are configured to mix / combine intensity values of pixels that depict intensity emissions from a target cluster being base called and intensity emissions from one or more adjacent clusters in a manner that maximizes the signal-to-noise ratio. The signal maximized in the signal -to-noise ratio is the intensity emissions from the target cluster, and the noise minimized in the signal-to-noise ratio is the intensity emissions from the adjacent clusters, i.e., spatial crosstalk, plus some random noise (e.g., to account for background intensity emissions). The equalizer coefficients are used as weights and the mixing / combining includes executing element-wise multiplication between the equalizer coefficients and the intensity values of the pixels to calculate a weighted sum of the intensity values of the pixels, i.e., a convolution operation.

[0104] RTA then performs base calling by fitting a mathematical model to the optimized intensity data. Suitable mathematical models that can be used include, for example, a k-means clustering algorithm, a k-means-like clustering algorithm, expectation maximization clustering algorithm, a histogram-based method, and the like. Four Gaussian distributions may be fit to the set of two-channel intensity data such that one distribution is applied for each of the four nucleotides represented in the data set. In some implementations, an expectation maximization (EM) algorithm may be applied. As a result of the EM algorithm, for each X, Y value (referring to each of the two channel intensities respectively) a value may be generated which represents the likelihood that a certain X, Y intensity value belongs to one of four Gaussian distributions to which the data is fitted. Where four bases give four separate distributions, each X, Y intensity value will also have four associated likelihood values, one for each of the four bases. The maximum of the four likelihood values indicates the base call. For example, if a cluster is “off’ in both channels, the base call is G. If the cluster is “off’ in one channel and “on” in another channel the base call is either C or T (depending on which channel is on), and if the cluster is “on” in both channels the base call is A.

[0105] In some implementations of RTA, the base calling errors get averaged out across many training examples. In some other implementations, the ground truth may be sourced using aligned genomic data, which may provide better quality because aligned genomic data may use reference genome and truth information that incorporate the knowledge gained from multiple sequencing platforms and sequencing runs to average out the noise. The ground truth may include base-specific intensity values (or feature values) that reliably represent intensity profiles of bases A, C, G, and T, respectively. A base caller like RTA base caller (926) base calls clusters by processing the sequencing images and producing, for each base call, color-wise intensity values / outputs. The color-wise intensity values may be considered base-wise intensity values because, depending on the type of chemistry (e g., 2-color chemistry or 4-color chemistry), the colors map to each of the bases A, C, G, and T. The base with the closest matching intensity profile is called.

[0106] A trainer may train base caller (926) and generate the trained coefficients of the sharpening masks using various training techniques. FIG. 6 shows one implementation of an adaptive technique that may be used to train base caller (926), e.g., using an offline or online mode. Here, the logic is y = x.h + d, where x is the input pixel intensities, h is the sharpening mask coefficients, d is the DC offset. In some implementations, x and h are row and column vectors respectively, with length 81. This vector model is equivalent to a dot product of 9 x 9 matrices representing input pixels and coefficients. The cost is the expected value of error squared. The gradient update moves each coefficient in a direction that reduces the expected value of error squared. Applying this update generates a new estimate of the coefficients that moves them in a direction that (on average) reduces the mean squared error (MSE). In some implementations, Mu is a small constant used to change the adaptation rate / convergence speed. A DC term update can be calculated in a similar way. A gain term update also can be calculated in a similar way.

[0107] In some implementations, since linear interpolation is applied on the coefficient sets, the updates are applied slightly differently in the following manner:

[0108] h(q, n+1) = h(q, n) + lambda q. mu.x(n). e(n)

[0109] In the equation above, h(q, n) is weight q at cycle n, lambda q is the linear interpolation weight for a particular set of coefficients and can include four updates per output due to linear interpolation in two dimensions. The recursive least-squares technique extends the least squares technique to a recursive algorithm.

[0110] C. Examples of Structural Variation Graph Genome Generation

[0111] In some scenarios, a secondary analysis may be performed iteratively while sequence reads are generated by a sequencing system such as systems (100, 500, 814) described herein. Secondary analyses may encompass both alignment of sequence reads to a reference sequence(eg., the human reference genome sequence) and utilization of this alignment to detect differences between a sample and the reference. Secondary analyses may enable detection of genetic differences, variant detection and genotyping, identification of single nucleotide polymorphisms (SNPs), small insertions and deletion (indels) and structural changes in the DNA, such as copy number variants (CNVs) and chromosomal rearrangements.[00112J By performing secondary analyses while sequence reads are generated, system (100, 500, 814) may determine preliminary variant calls iteratively in real-time (or with zero or low latency). Final results of variant determinations may be available soon after (or immediately after) the end of a sequencing run. Alternatively, a sequencing run may be terminated early if variant calls are available with sufficient confidence during the run. In some scenarios, only information related to variant determinations (e.g., variant calls) is transferred off the sequencing system (100, 500, 814). This may decrease, or minimize, the data bandwidth required in comparison to performing the variant determinations in a system that is external. In addition, only variant information may be sent to a computing system (e.g., a cloud computing system) for further processing. In this example, sequencing runs may be terminated prior to completion of an entire sequencing process. For example, if the identity of a pathogen of interest is determined after a number of sequencing cycles of a sequencing run, the sequencing run may be terminated. Thus, the time to a particular answer (e.g., pathogen identification) may be decreased. In some implementations, outputs and intermediate results of system (100, 500, 814) may include histograms of duplicates, exact matches, single and double SNPs, and single and double indels.

[0113] VI. Examples of association mapping on single-cell RNA sequencing data

[0114] The described sequencing systems and methodologies described herein may be utilized in various contexts for the generation of single-cell RNA sequencing data that may be useful in determining associations between features of the data set (e.g., perturbations, accessibility, genotype, or, more generally, any characteristic of value that is part of the variation between the cells) and gene expression. By way of context, it may be useful to associate signals associated with a feature of a dataset with the level of gene expression measured from single-cell sequencing data (scRNA-seq). Mapping such signals onto scRNA-seq data is useful with respect to a variety of assays including, but not limited to: (1) perturbation sequencing (Perturb-seq) in which a library of perturbations is introduced into cells and the association between each perturbation and single-cell gene expression is determined; (2) expression quantitative trait loci (eQTL) mapping, in which expression levels are associated with genotypes; and (3) peak-to-gene mapping, in which chromatin accessibility signals and gene expression measurements from the same cells (i.e., multiomic sequencing) are associated. In such examples, the scRNA-seq readout is associated with differentiating features associated with the cells, such as their perturbation, their genotype, or their accessibility signals. As discussed herein, various embodiments of the present technique comprise a computational pipeline that accurately maps such associations in a computationally rapid and efficient manner. By way of example, a single-cell data sample comprising hundreds of thousands of cells or more may be processed in six hour or less (e.g., 1 hour, 2 hours, 3 hours, 4 hours, 5, hours or 6 hours). In certain such implementations one or more graphical processing units (GPUs) are employed to facilitate fast and efficient performance while maintaining accuracy.

[0115] The presently described technique may be particularly useful in the context of large data sets (e g., data sets over 50,000 cells or 100,000 cells, with each cell constituting a single-cell data set). In particular, with the advent of next generation sequencing technology and corresponding developments in singe-cell isolation and techniques, large-scale scRNA-seq techniques have become viable for obtaining massive scRNA-seq data sets (e.g., data sets encompassing single-cell readout data for hundreds of thousands to millions of cells). Such large-scale scRNA-seq data sets cannot be processed in a reasonably practical manner for association mapping using conventional approaches based on correlation and / or model refitting for each feature set. By way of example, in Perturb-seq operations thousands of perturbations are associated with the expression of each gene. Such Perturb-seq datasets may scale up to millions of cells, which is an untenable extension of conventional techniques.

[0116] With this in mind, and as discussed in greater detail below, the presently described technique incorporates various features and benefits that facilitate accurate analysis at scale. First, in certain implementations the technique employs a model that accounts for over-dispersion and for structure present in the single-cell sequencing data, thereby enabling accurate mapping of association signals. Second, in many scRNA-seq association tasks, each gene is associated with thousands of predictors. As discussed herein, in certain implementations the analysis is streamlined so that a model for the expression of each gene is only fitted once, while each of the predictors is tested by estimating the effect of that predictor on the model without explicitly refitting the model(an analysis referred to as a “score test” herein in which a corresponding score based on the derivative of a likelihood function is determined). Third, to further facilitate large-scale analysis, in certain implementations the score test analysis is done in parallel (such as via multi-threaded processing) on CPUs or GPUs. In context where GPUs are employed, such a parallel implementation may leverage the ability of the GPUs to provide large-scale parallelized matrix multiplication to achieve improvement over traditional approaches, such as in terms of speed of implementation, which may be in a matter of hours for a sample comprising hundreds of thousands of cells.

[0117] As noted above and discussed herein, association mapping in the present context refers to the mapping of some feature (e.g., a perturbation from Perturb-seq, genotype, accessibility, and so forth) associated with a respective cell gene expression to gene expression within the given cell. As may be appreciated, in such a context the feature data and the expression data for each respective cell can each be represented in respective matrices. By way of example, and turning to FIG. 7, a feature matrix 1000 and an expression matrix 1004 for a respective cell are shown. In this example, the feature matrix 100 comprises values of the signal (a) measured for each feature ( / ) in each cell (z) of the sample, where the feature may correspond to a perturbation, a genotype, an accessibility, and so forth. By way of example, in a Perturb-seq data set, the feature signal (cz) will be the signal of the perturbation. Similarly, in multiomic sequencing the feature signal (a) will be peak accessibility. Correspondingly, for eQTL mapping, the feature signal (a) will be the genotype in each single nucleotide polymorphism (SNP). With respect to the expression matrix 1004, this matrix comprises values of the RNA expression (e) of each gene (g) in each cell (z).

[0118] Such single-cell data may exhibit various characteristics that may make analysis more difficult. By way of example, one characteristic that may be present in single-cell data is overdispersion of count data. As used herein, over-dispersion may be understood to be a greater than expected variability (i.e., statistical dispersion) in a set of data (e.g., count data) based on the underlying model. Such over-dispersion may distort the estimated standard errors and test statistics of the overall goodness-of-fit measures, which may render fitting of the observed data difficult or infeasible. Such over-dispersion of the count data may render simple correlation and / or linear regression modeling techniques ineffective. An example of over-dispersion in single-cell count data is exhibited in FIG. 8, depicting a plot of expression variance against average expressionfor single-cell count data and exhibiting over-dispersed characteristics.

[0119] Further, single-cell count data may exhibit structure that must be modeled or otherwise accounted for, resulting in complex models including terms to correct for the underlying structure within the count data. By way of example, such underlying structure within the single-cell data may be attributed to the different cell types that may be studied together in a given sample. In practice, such structure within the single-cell count data may hinder casual inference and potentially result in excessive false-positive associations. An example of structure present within such data is provided by FIG. 9, which exhibits two plots. The leftmost plot depicts RNA (i.e., gene expression) data and an example of underlying structure that may be present in such singlecell data drawn from different sources that may be present within a sample. Similarly, the rightmost plot depicts accessibility (i.e., chromatin binding) data and an example of underlying structure that may be present in such single-cell data drawn from different sources that may be present within a sample.

[0120] With respect to a model that may be employed to associate single-cell data for feature and gene expression associations, an example is provided below. In this example, feature data may correspond to gene accessibility data within different cells.

[0121] By way of example, in one embodiment testing for such associations may be performed using a computational pipeline in accordance with equations (1) and (2) shown below:(1) yt-Pois^i)(2)>where Pois is indicative of a Poisson distribution, A is peak accessibility, S is total counts in the cell, AT is the fraction of mitochondrial reads, the termPCn,ifipcn Laccounts for structure within the data, and the Observation-Level Random Effect term (OLRE) addresses over-dispersion within the data. While a Poisson mixed effect model (e.g., a Poisson generalized linear mixedmodel (GLMM)) with OLRE is one possible suitable approach, in general the model may be any suitable generalized linear model, including but not limited to negative binomial models.

[0122] With respect to the principal components (PC) term accounting for structure, in the provided implementation of the model, X expression principal components may be employed. In one embodiment, X equals 10 principal components, though fewer or more PCs may be employed as needed to account for the desired or observed extent of structure within the data. In particular, and as may be appreciated the principal component term provides a mechanism to break down or partition variation that is observed between all of the different cells into axes of variation (corresponding to the number of principal components) that can be used as covariates in the model to account for the global variation that is present. In this manner, the total or global variation may be used to determine what the main (e.g., statistically significant) axes of variation are within the data. These axes (or PCs) may in turn be used to remove signal unrelated to the feature-gene relationships being studies, such as variation attributable to structure within the data.

[0123] With respect to the Observation-Level Random Effect term (OLRE) term accounting for over-dispersion within the single-cell data, adding an observation-level random intercept in which every cell is its own ‘group’ may be used as an alternative to negative binomial regression. With respect to this OLRE approach, the following framework may be useful in understanding and implementing such observation-level random intercepts:(3) y1~Pois Xi)(4) In Aj = X(3 -I- uo i<

[0124] In accordance with the above frame-work and in the context of providing an observation -level random effect, the number of observed molecules in the cell may be considered as coming from a Poisson distribution and the of molecules expressed in each cell follows thedescribed lognormal distribution. By including such an OLRE term, over-dispersion in the singlecell data is addressed by the presently described model and processor-implemented pipeline. While over-dispersion and structure are noted specifically in this example, more generally the model may comprise terms to address multi-cell type factors and / or covariates addressing differences between cells.[00125J In practice, the output of a processor implemented pipeline incorporating the model outlined by equation (2) may take the form of a matrix of p-values and corresponding strengths of effect for evaluated genetic associations or relationships being evaluated. As part of configuring the performance of the pipeline and / the output a user may specify thresholds around features being evaluated (e.g., for multiomic analysis with accessibility data, a cut-off for each peak to be evaluated (e.g., 500 kb in both directions for a 1 Mb total zone centered around each peak or, alternatively, 100 kb in both directions, 1 0 kb in both directions, 200 kb in both directions, 250 kb in both directions, and so forth). In such implementations, the output may then include a signal for all of the genes within a corresponding threshold-defined region.

[0126] In certain approaches, significance testing of such model-based techniques may be performed by comparing the model implemented in the pipeline (i.e., a full model) to a null model that does not include the respective feature data by way of a likelihood ratio test (LRT). As may be appreciated, a likelihood function typically is employed to measure how well a given statistical model explains the actual observed data. Such a likelihood function calculate a probability that the data would be observed under different parameter values of the model. In practice, such as likelihood model may be generated using the joint probability distribution of the random variable used to generate the observations. Thus, when assessed on the observed data points, the likelihood function can be used to assess the model parameters. With respect to a likelihood ratio test, such a test involves the comparison of the goodness of fit of the two models (i.e., the full and the null models) for a given feature value (e.g., perturbation, accessibility, genotype, and so forth) and gene based on the ratio of their likelihoods. In this manner the likelihood ratio test allows assessment of whether the ratio its natural logarithm is significantly different from zero.

[0127] With this in mind, in such approaches the model is fit twice, once without the feature data (the null model) and once with the feature data (i.e., the full model). Further, this may be performed for hundreds, thousands, or tens of thousands of genes for a given sample. The twomodels (i.e., the full model and the null model) are then compared based on a likelihood ratio to determine which better fits the data. Based on the tested associations, pooled genome-wide permutations can be determined and used to identify associations between gene expression and features in the single-cell data.

[0128] However, there are applications (e.g., Perturb-seq) where one may be measuring thousands of feature variations (e.g., perturbations in the Perturb-seq context) and the effects of thousands of perturbations on thousands of genes. Such a model fitting approach is not feasible, in a computational sense, for every perturbation. More generally, in various scRNA-seq association tasks, each gene may be associated with thousands of predictors (e.g., feature permutations or variations), each requiring a separate model-fitting step and the number of fitting steps that would be involved is not computationally feasible.

[0129] In accordance with the presently described techniques, these various problems are addressed. In particular, as noted above, refitting the model for each feature permutation is computationally impractical. In accordance with the presently described approach, instead of refitting the model for each feature permutation, a statistical approach is employed in which a score test metric is employed to determine the significance of a feature permutation value. In accordance with such an approach, only the null model (one null model per gene) is fit. That is, instead of the number of model fittings being F x G (i.e., feature permutations (F) X # of genes (G)), the number of model fittings is simply G.

[0130] Each feature (e.g., perturbation, genotype, accessibility) is examined by calculating a score test for the respective feature. In such a context, the score test may be understood to be the derivative of the likelihood function with the contribution of the additional parameter being added (i.e., the feature permutation contribution) or, more generally, the gradient (e.g., the sum or combination of the derivatives) of the likelihood function with respect to each feature set. In this manner, without fitting a full model, one can estimate how much better (or worse) the model performs with the feature added (i.e., fitted). An example of this may be seen in FIG. 10, in which a plot of a likelihood function 1100 is illustrated along with a corresponding derivative 1104 at the point of the null hypothesis, where the null hypothesis is that the covariate for the feature (e.g., perturbation) is zero, which is essentially equivalent to fitting the model without the covariate in question. The example of FIG. 10 illustrates a situation in which the null hypothesis is not trueand the feature (e.g., perturbation) in question is associated with gene expression. Correspondingly, the point associated with the null hypothesis (or null model) 1104 in this example is not the true maximum and the derivative at that point is not 0, thereby demonstrating that there is an association between the feature in question and gene expression. In this example, the slope is substantially different from 0, potentially indicating a significant association.

[0131] Thus, the slope of the derivative 1104 is indicative of the effect of the feature under analysis (e.g., the extent of association between the respective feature or feature permutations and the respective gene in question). For example, a flat or relatively flat slope is indicative of no effect or a weak effect (i .e., weak association). Conversely, the further from zero the slope of the derivative is, the greater the effect or association. That is, the further the score test is from zero, the greater the likelihood that the likelihood function would change if fitted for the feature being analyzed. Correspondingly, one or more thresholds may be defined for assessing the value of a given score test so as to quantitatively or qualitatively characterize the assessment of the effect of that feature on gene expression. That is, a score test value greater than a given threshold may be characterized as having a significant association, while a score test value less than the same threshold (or a separate threshold) may be characterized as having negligible or insignificant association.

[0132] In this manner, such association testing is simplified. For example, for an assessment of 20,000 genes (G) over F feature permutations, the model need only be fitted G times, as opposed to F x G times. That is, a derivative for each feature is calculated rather than a refit of the model for each feature.

[0133] With the preceding in mind, FIGS. 11A-11C and FIGS. 12A-12C provide two examples based on the above-described approach, the first (FIGS. 11A-11C) illustrating a perturbation having a strong effect and the second (FIGS. 12A-12C) illustrating perturbation having a weak effect. Turning to FIG. 11 A, the plot illustrates the model residuals from fitting the null model as a function of the raw counts. Cells receiving the perturbation are shown as lighter filled circles. As may be observed, the lighter filled circles have corresponding negative residuals, in this example highly negative residuals, demonstrating that these cells have lower expression than would be expected from the model. FIG. 1 IB depicts the raw counts plotted with respect to the perturbation state. FIG. 11C depicts the distribution of model residuals plotted with respect tothe perturbation state. In this example where there is a strong perturbation effect, the effect is demonstrated by the downward-shifted (largely below 0) distribution. Correspondingly, the score test value (i.e., score statistic) for this example is very significant, with a value of 1506.90 and a - value of ~0.

[0134] Conversely, turning to FIGS. 12A-12C, which depict corresponding plots for a different sample in which there is little or no perturbation effect, it may be observed that FIG. 12A exhibits that the residuals for the perturbed cells are distributed throughout the plot, with no discernible trend. This is also evident in FIG. 12C which, unlike the downward-shifted distribution of FIG. 11C, is instead largely centered about 0. This effect is not readily discerned without modeling, which is observable by looking at FIGS. 1 IB and 12B which do not utilize modeling. With respect to the weak or no effect example, the score test value (i.e., score statistic) for this example is very small, with a value of 0.02 and a - value of 0.89.

[0135] With the preceding in mind, in certain implementations it may be useful to employ one or more graphics processing units (GPUs) to perform one or both of the null model fitting and the score test value derivation. By way of example, in implementations employing two or more GPUs, processing may be performed in a parallelized data flow. In such implementations, the use of GPUs and / or parallel data handling may increase the speed at which associations can be tested relative to other processors and / or data processing flows. Conversely, in other implementations it may be suitable to employ one or more central processing units (CPUs) to perform one or both of the null model fitting and the score test value derivation. In practice, a Jax backend may be employed in an implementation which is configured to utilize a GPU or GPUs if found to be present and which otherwise defaults to using a CPU or CPUs if no GPU is seen by the backend. Similarly, other comparable deep learning library tools may provide similar functionality of identifying and preferentially using a GPU if present and otherwise defaulting to use of a CPU.

[0136] In a more general context, without respect to the type of processor employed, parallelization of the described techniques may be accomplished via multithreading techniques that may be implemented at the level of the sequencer device itself, at a downstream workstation, or at the server or datacenter level. Such parallelization may be particularly suitable in the context of large-scale analyses. For example, the presently described techniques may be applied in the single-cell analysis of samples comprising a hundred thousand cells to millions of cells (e.g., a halfa million cells to millions of cells). Further, at such scales the effects of interactions otherwise missed or muted by the use of derivatives for feature effect estimation (as opposed to full fitting of each model for feature effects and interactions) may be countered due to the size of the sample involved. That is, the effects of interactions missed by the move to derivative based estimation may be addressed by the move to millions of cells being assessed in the single-cell sample as opposed to 20,000 such cells in conventional approaches.

[0137] Turning to FIG. 13, an example process flow is depicted. In contrast to prior approaches in which both the null model and a feature fitted model were calculated, in this example workflow only the null model is calculated (block 1140). In the depicted example, single-cell RNA data 1144 (e.g., expression data in conjunction with feature data) is accessed or acquired (such as an output of a next generation sequencing system used to implement a single-cell RNA protocol). Further data or sample manipulation steps may be performed in terms of preparing data for analysis, such as the depicted “perturbation / cell assignment” step 1148 in a Perturb-seq type analysis and the depicted “preprocessing and clustering” step 1152 which may be performed as part of addressing structure and other issues present in single-cell data. As depicted in this example, such steps may be performed independent of one another, such as in parallel.

[0138] With respect to the subsequent calculation of the null model at step 1140, in one implementation the null model may be calculated or fitted in accordance with equation (1) above and:(9) y ~ Norm(0,As noted herein, in such an embodiment the null model is calculated only once using all of the data for a given gene.

[0139] Subsequently, at step 1156, a score test is calculated based on derivatives for each feature permutation present for the gene in question. In the depicted example, the score test is calculated using one or more GPUs, though in practice one or more CPUs may be used instead. The score test calculation may be in accordance with:(10) logwherein random effects are approximated using maximum a posteriori probability estimates (MAP). In some implementations confidence intervals may be obtained by performing permutations on this process. However, in other instances the estimation achieved is sufficient and such additional steps are not needed or performed.[00140J With the preceding in mind, a case study is depicted that relates to a Perturb-seq type screen. In such a Perturb-seq type screen, microRNAs (e.g., small RNA sequences that may be regulators) are inserted into cells of a sample and the effect of the inserted on each of the genes in questions are observed. Such a screen may be useful because it is understood how the microRNAs target and, correspondingly, what the expected targets are. As a result, such a Perturb-seq tool may be useful in assessing how the described analytic pipeline performs compared to looking naively at the RNA expression data and performing the respective feature changes. That is, use of the Perturb technique provides targets that are known to be true as well as targets that cannot be true. Therefore, this known targeting can be used to calculate metrics of how well the presently described classifier (i.e., the score test) is performing.

[0141] Further, in contrast to other Perturb-seq type analyses in which a control is instead employed, the presently described techniques allow one to take, for every perturbation, all of the other perturbations as part of the null model. In particular, the changes of interest will not be shared across all different perturbations but will instead be specific to certain of the perturbations. This improves power considerably with respect to the presently described approaches because, instead of the null model being based on only a few hundred cells (which one hopes are representative of a baseline), the null model can instead be based on the entire experiment or sample run as the baseline, which as noted above may be millions of cells.

[0142] As depicted in FIGS. 14 and 15, an improvement was observed with respect to such metrics when using the presently described techniques. In particular, with respect to the area under the receiver operating characteristic (AUROC) curve, a greater than 100 basis point improvement was observed, as seen in FIG. 14. Similarly, turning to FIG. 15, an improvement (greater than 50basis points) in the context of the area under the precision recall curve (AUPRC) was also observed using the presently described techniques. As may be appreciated, AUROC and AUPRC are both known metrics for evaluating the performance of a classifier and take into account both false positives and false negatives, corresponding to sensitivity and specificity.

[0143] With the preceding in mind, the presently described techniques are suitable for use as part of a tool or suite of tools employed in the tertiary analysis of genetic data, in which biological data mining and interpretation tools are used to abstract processed sequencing data and derived analytics into insights related to physiological or phenotypic phenomena, such as may be related to the causes of disease and / or how to treat or prevent such diseases.

[0144] VII. Experimental Results

[0145] A. Comparison with Prior Techniques

[0146] By way of proof of concept, the techniques described herein were compared against existing techniques. Data sets employed for the purpose of comparison included STRING, a public database the includes protein-protein interactions, and CORUM, a database of protein complexes. Such datasets were deemed suitable as it has been shown that perturbations whose effects are correlated are enriched with interacting pairs and pairs that are in the same protein complex. Correspondingly, this information may be used to benchmark differential expression results, where differential expression results that better prioritize interacting proteins indicate better performance.

[0147] In terms of the prior or conventional techniques selected for comparison with the presently described techniques, two were selected. The first reference technique is referred to herein as Deseq2, which is a bulk RNA-seq analysis pipeline that is based on negative binomial regression. In accordance with this technique, for single cell data, cells that received the same perturbation are aggregated (i.e., pseudo-bulked) and differential expression was observed between each perturbation and set of control perturbations.

[0148] The second reference technique is referred to herein as “Replogle” and follows a processing pipeline described in a previously published genome- wide perturb seq study (Replogle JM, et al. (2022) Mapping information-rich genotype-phenotype landscapes with genome-scale Perturb-seq. Cell. 185(14):2559-2575.E28). In accordance with the Replogle pipeline, gene expression is normalized and z-scores are calculated comparing each perturbation to the controlperturbation cells. For the purpose of comparison, since the Replogle study as published selected specific perturbations which were deemed by the authors as “high quality” for publication, performance relative to Replogle herein was evaluated against both their selected data set as well as against their entire data set.

[0149] As discussed herein, the comparison processing pipeline based on the presently described techniques (i.e., the “present techniques”) employs a baseline negative binomial model that is fit for each gene, which is followed by fast association testing using the presently described score test.

[0150] With the preceding in mind, the experimental results presented below describe analysis of genome-wide (GW) Perturb-seq experiments in which a library of CRISPR interference (CRISPRi) perturbations against all protein-coding genes was used to profile the effect of downregulating each protein coding gene.

[0151] It may be noted that the comparison presented in sub-section (ii) below were performed on a 1.9 million cell subset of K562 data. The comparison presented in sub-section (iii) shows performance across the full dataset of 8.2 million K562 cells.

[0152] (i) Enrichments Across Methods

[0153] With the preceding in mind, and turning to FIGS. 16-17, protein-protein interaction (PPI) enrichments are plotted for the different methodologies noted above. Turning to FIG. 16, enrichment of STRING data having a score of > 800 (vertical axis) is plotted against the top N° / o of gene pairs with the highest correlation (horizontal axis) for the three techniques described above for full gene features and subsets of gene features. As may be observed, the presently disclosed techniques outperform the two prior techniques. Turning to FIG. 17, the equivalent number of top correlated perturbations are compared (horizontal axis) instead of taking the top percentile. In this context as well, it may be observed that the presently disclosed techniques outperform the two prior techniques. Similarly, and turning to FIG. 18, protein complex enrichments are instead plotted for the different methodologies noted above. In particular, the enrichment of CORUM protein complex data is depicted as having enrichment in CORUM pairs (vertical axis) plotted against the top N° / o of gene pairs with the highest correlation (horizontal axis) for the three techniques described above for full gene features and subsets of gene features. As in the priorresults, the presently disclosed techniques outperform the two prior techniques.

[0154] (ii) Performance Across Cell Lines

[0155] In a further study performance of the presently disclosed technique relative to the two reference techniques across cell lines was investigated using genome-wide (GW) CRISPR interference (CRISPRi) data. The CRIPSRI data included approximately 42 million cells that were profded and processed in three batches of experiments (GW1-3), broken down as follows.Table 1

[0156] With respect to this subject matter, and turning to FIG. 19, a plot of target coverage by cell line is provided along with a table of median coverage by cell line. In particular, the plot depicts the cumulative distribution of guide target coverage, with fraction of target (cumulative) represented on the vertical axis and number of cells with target perturbation represented on the horizontal axis.

[0157] With this in mind, FIG. 20 depicts STRING PPI enrichments across cell lines using the presently described techniques. In particular, FIG. 20 depicts the enrichment of STRING data having a score of > 800 (vertical axis) plotted against the top A% of gene pairs with the highest correlation (horizontal axis) for the various cell lines. As shown, the present techniques achieve strong enrichments across the different investigated genome-wide Perturb-seq data sets. Turning to FIG. 21, corresponding results are depicted for the top A gene pairs with the highest correlation (horizontal axis) instead of taking the top percentile. In this context as well, it may be observed that the presently disclosed techniques achieve strong enrichments across the investigated genomewide Perturb-seq data sets. Turning to FIGS. 22 and 23, corresponding analyses and results are depicted for CORUM data instead of STRING data. As with the STRING investigation results, the presently disclosed techniques can be observed to achieve strong enrichments across the investigated genome-wide Perturb-seq data sets.

[0158] (iii) Runtime Performance

[0159] As discussed herein, the presently described techniques provide improvements in terms of computational efficiency and performance. In the tested implementation the presently described techniques were implemented so as to split processing across genes using a Nextflow pipeline. With this setup, the top 12,000 expressed genes were studied, which is larger than the expected number of expressed genes in the K562 cell line), and these 12,000 genes were split into 1,200 chunks parallelized across the cluster.

[0160] With this in mind, and turning to FIG. 24, CPU usage of the presently described techniques on 8.2 million cells is depicted. In this example, the processes RunBaselineChunk and RunScoreTestChunk were run on 1,200 total chunks. In this example, RunPCAAndCovariates is 1 chunk and InitializeMatricesR is 2 processes. Turning to FIG. 25, memory usage (i.e., physical memory usage) is depicted for the same run as the CPU usage of FIG. 24. Similarly, FIG. 26 depicts job duration (i.e., task execution time) for the same run as the CPU usage of FIG. 24. Based on these results, computational efficiency and runtime represent improvements over prior methodologies.

[0161] VIII. Examples of Combinations

[0162] The following examples relate to various non-exhaustive ways in which the teachings herein may be combined or applied. The following examples are not intended to restrict the coverage of any claims that may be presented at any time in this application or in subsequent filings of this application. No disclaimer is intended. The following examples are being provided for nothing more than merely illustrative purposes. It is contemplated that the various teachings herein may be arranged and applied in numerous other ways. It is also contemplated that some variations may omit certain features referred to in the below examples. Therefore, none of the aspects or features referred to below should be deemed critical unless otherwise explicitly indicated as such at a later date by the inventors or by a successor in interest to the inventors. If any claims are presented in this application or in subsequent filings related to this application that include additional features beyond those referred to below, those additional features shall not be presumed to have been added for any reason relating to patentability.

[0163] Example 1

[0164] A method for testing a plurality of features in single-cell data, the method comprising: fitting a model for gene expression for single-cell data, wherein the model corresponds to one model per gene; and estimating the effect of each feature on the model without refitting the model.

[0165] Example 2

[0166] The method of Example 1, wherein estimating the effect of each feature on the model comprises determining the derivative of a likelihood function with respect to each feature.

[0167] Example 3

[0168] The method of Example 2, further comprising generating a score based on the derivative.

[0169] Example 4

[0170] The method of Examples 1-3, wherein the features comprise perturbations, accessibilities, or genotypes.

[0171] Example 5

[0172] The method of Examples 1-4, wherein one or both of fitting the model or estimating the effect of each feature are performed in parallel.

[0173] Example 6

[0174] The method of Examples 1-5, wherein one or both of fitting the model or estimating the effect of each feature are performed using graphics processing units (GPUs).

[0175] Example 7

[0176] The method of Examples 1-6, wherein one or both of fitting the model or estimating the effect of each feature are performed in parallel on an array of GPUs.

[0177] Example 8

[0178] The method of Examples 1-7, wherein the single-cell data comprises discrete cell data for 100,000 cells or more in a sample.

[0179] Example 9

[0180] The method of Examples 1-8, wherein the single-cell data comprises discrete cell data for millions of cells or more in a sample.

[0181] Example 10

[0182] The method of Examples 1-8, wherein testing the plurality of features in single-cell data comprises testing hundreds of thousands of cells or more within 6 hours or less.

[0183] Example 11

[0184] The method of Examples 1-10, wherein the model comprises terms to address structure and over-dispersion in the single-cell data.

[0185] Example 12

[0186] The method of Examples 1-11, wherein the model comprises terms to address multicell type factors and / or covariates addressing differences between cells.

[0187] Example 13

[0188] A computational pipeline for estimating perturbation effects in single-cell data, comprising: accessing or acquiring single-cell RNA expression data; determining perturbation data for each cell; performing preprocessing and clustering of the single-cell RNA expression data; calculating a model based on the perturbation and RNA expression data for each cell; and estimating the effect of each perturbation on the model based on respective derivatives calculated for each perturbation.

[0189] Example 14

[0190] The method of Example 13, wherein the acts of determining perturbation data and performing preprocessing and clustering are performed in parallel on an array of graphics processing units (GPUs).

[0191] Example 15

[0192] The method of Examples 13-14, wherein the model is not refitted when estimating the effect of each perturbation.

[0193] Example 16

[0194] The method of Examples 13-15, wherein the single-cell data comprises discrete cell data for 100,000 cells or more.

[0195] Example 17

[0196] The method of Examples 13-16, wherein the single-cell data comprises discrete cell data for millions of cells or more.

[0197] Example 18

[0198] The method of Examples 13-17, wherein the model comprises terms to address structure and over-dispersion in the single-cell data.

[0199] Example 19

[0200] A nucleic acid sequencing system, comprising: one or more processors configured to execute code and process data; and one or more memory devices configured to store processor executable code and nucleic acid sequence data. The processor executable code, when executed by the one or more processors, causes acts to be performed comprising: fitting a model for gene expression for single-cell data, wherein the model corresponds to one model per gene; and estimating the effect of each feature on the model without refitting the model.

[0201] Example 20

[0202] The method of Example 19, wherein the one or more processors comprise at least one graphics processing units (GPUs).

[0203] VIII. Miscellaneous

[0204] While the foregoing examples are provided in the context of a system (100) that may be used in nucleotide sequencing processes, and in particular the sequencing of methylation data, the teachings herein may also be readily applied in other contexts, including in systems that perform other processes (i.e., other than nucleotide sequencing procedures). The teachings herein are thus not necessarily limited to systems that are used to perform nucleotide sequencing processes.

[0205] It is to be understood that the subject matter described herein is not limited in its application to the details of construction and the arrangement of components set forth in the description herein or illustrated in the drawings hereof. The subject matter described herein is capable of other implementations and of being practiced or of being carried out in various ways. Also, it is to be understood that the phraseology and terminology used herein is for the purpose of description and should not be regarded as limiting. As used herein, an element or step recited in the singular and proceeded with the word “a” or “an” should be understood as not excluding plural of said elements or steps, unless such exclusion is explicitly stated. Furthermore, references to “one example” are not intended to be interpreted as excluding the existence of additional examples that also incorporate the recited features. The use of “including,” “comprising,” or “having” and variations thereof herein is meant to encompass the items listed thereafter and equivalents thereof as well as additional items.

[0206] When used in the claims, the term “set” should be understood as one or more things which are grouped together. Similarly, when used in the claims “based on” should be understood as indicating that one thing is determined at least in part by what it is specified as being “based on.” Where one thing is required to be exclusively determined by another thing, then that thing will be referred to as being “exclusively based on” that which it is determined by.

[0207] Unless specified or limited otherwise, the terms “mounted,” “connected,” “supported,” and “coupled” and variations thereof are used broadly and encompass both direct and indirect mountings, connections, supports, and couplings. Further, “connected” and “coupled” are not restricted to physical or mechanical connections or couplings. Also, it is to be understood that phraseology and terminology used herein with reference to device or element orientation (such as, for example, terms like “above,” “below,” “front,” “rear,” “distal,” “proximal,” and the like) are only used to simplify description of one or more examples described herein, and do not alone indicate or imply that the device or element referred to must have a particular orientation. In addition, terms such as “outer” and “inner” are used herein for purposes of description and are not intended to indicate or imply relative importance or significance.

[0208] It is to be understood that the above description is intended to be illustrative, and not restrictive. For example, the above-described examples (and / or aspects thereof) may be used in combination with each other. In addition, many modifications may be made to adapt a particularsituation or material to the teachings of the presently described subject matter without departing from its scope. While the dimensions, types of materials and coatings described herein are intended to define the parameters of the disclosed subject matter, they are by no means limiting and instead illustrations. Many further examples will be apparent to those of skill in the art upon reviewing the above description. The scope of the disclosed subject matter should, therefore, be determined with reference to the appended claims, along with the full scope of equivalents to which such claims are entitled. In the appended claims, the terms “including” and “in which” are used as the plain-English equivalents of the respective terms “comprising” and “wherein.” Moreover, in the following claims, the terms “first,” “second,” and “third,” etc. are used merely as labels, and are not intended to impose numerical requirements on their objects. Further, the limitations of the following claims are not written in means — plus-function format and are not intended to be interpreted based on 35 U.S.C. §112(f) paragraph, unless and until such claim limitations expressly use the phrase “means for” followed by a statement of function void of further structure.

[0209] The following claims recite aspects of certain examples of the disclosed subject matter and are considered to be part of the above disclosure. These aspects may be combined with one another.

Claims

What is claimed is:

1. A method for testing a plurality of features in single-cell data, the method comprising:fitting a model for gene expression for single-cell data, wherein the model corresponds to one model per gene; andestimating the effect of each feature on the model without refitting the model.

2. The method of claim 1, wherein estimating the effect of each feature on the model comprises determining the derivative of a likelihood function with respect to each feature.

3. The method of claim 2, further comprising generating a score based on the derivative.

4. The method of claim 1, wherein the features comprise perturbations, accessibilities, or genotypes.

5. The method of claim 1, wherein one or both of fitting the model or estimating the effect of each feature are performed in parallel.

6. The method of claim 1, wherein one or both of fitting the model or estimating the effect of each feature are performed using graphics processing units (GPUs).

7. The method of claim 1, wherein one or both of fitting the model or estimating the effect of each feature are performed in parallel on an array of GPUs.

8. The method of claim 1, wherein the single-cell data comprises discrete cell data for 100,000 cells or more in a sample.

9. The method of claim 1, wherein the single-cell data comprises discrete cell data for millions of cells or more in a sample.

10. The method of claim 1, wherein testing the plurality of features in single-cell data comprises testing hundreds of thousands of cells or more within 6 hours or less.

11. The method of claim 1, wherein the model comprises terms to address structure and overdispersion in the single-cell data.

12. The method of claim 1, wherein the model comprises terms to address multi-cell type factors and / or covariates addressing differences between cells.

13. A computational pipeline for estimating perturbation effects in single-cell data, comprising:accessing or acquiring single-cell RNA expression data;determining perturbation data for each cell;performing preprocessing and clustering of the single-cell RNA expression data; calculating a model based on the perturbation and RNA expression data for each cell; and estimating the effect of each perturbation on the model based on respective derivatives calculated for each perturbation.

14. The computational pipeline of claim 13, wherein the acts of determining perturbation data and performing preprocessing and clustering are performed in parallel on an array of graphics processing units (GPUs).

15. The computational pipeline of claim 13, wherein the model is not refitted when estimating the effect of each perturbation.

16. The computational pipeline of claim 13, wherein the single-cell data comprises discrete cell data for 100,000 cells or more.

17. The computational pipeline of claim 13, wherein the single-cell data comprises discrete cell data for millions of cells or more.

18. The computational pipeline of claim 13, wherein the model comprises terms to address structure and over-dispersion in the single-cell data.

19. A nucleic acid sequencing system, comprising:one or more processors configured to execute code and process data; andone or more memory devices configured to store processor executable code and nucleic acid sequence data, the processor executable code, when executed by the one or more processors, causes acts to be performed comprising:fitting a model for gene expression for single-cell data, wherein the model corresponds to one model per gene; andestimating the effect of each feature on the model without refitting the model.

20. The nucleic acid sequencing system of claim 19, wherein the one or more processors comprise at least one graphics processing units (GPUs).

Images