Sequencing analysis method of single cell transcriptome based on ONT sequencing

Abstract

The invention discloses a sequencing analysis method of a single cell transcriptome based on ONT sequencing, and relates to the technical field of gene sequencing analysis. The method comprises the following steps: preparing and constructing a cDNA library of a sample to be detected, and sequencing to obtain third-generation original data; performing quality control on the third-generation original data to obtain third-generation effective data; performing function annotation on the reference genome; software ONT-SC-pipeline is used for carrying out data splitting, comparison and quantification, and quantitative data is obtained; sequentially carrying out basic analysis, Marker gene enrichment analysis and advanced analysis on quantitative data; the basic analysis comprises top gene distribution analysis, hypervariant gene analysis, clustering and grouping analysis, correlation analysis of each cell group, dimension reduction analysis and Marker gene identification; the advanced analysis comprises cell annotation analysis, cell communication analysis, cell trajectory analysis, variable shear analysis and fusion gene analysis.

Description

Technical Field

[0001] The present invention relates to the technical field of gene sequencing analysis, and in particular to a sequencing analysis method for a single cell transcriptome. Background Art

[0002] In biomedical research and clinical practice, single-cell analysis technology is playing an increasingly critical role in revealing cellular heterogeneity, analyzing complex biological processes and disease pathogenesis. Traditional population-level analysis methods mask subtle differences between cells, which are crucial for understanding disease development, cellular decision-making processes, and responses to environmental changes. Single-cell RNA sequencing (scRNA-seq) technology has greatly promoted the detection and quantification of thousands of RNA molecules from a single cell, providing a unique perspective for in-depth exploration of the internal mechanisms of biological systems.

[0003] Despite significant progress in scRNA-seq technology, its limitations in capturing and analyzing full-length RNA molecules have weakened the complete transcript information obtained from single cells and restricted our comprehensive understanding of post-transcriptional modifications, alternative splicing events, and gene expression regulation mechanisms.

[0004] In this context, Oxford Nanopore Technologies (ONT)'s single-cell full-length sequencing technology came into being, which provides an unprecedented complete view of transcripts by directly sequencing full-length RNA molecules in single cells. This technology not only reveals 5' and 3' ends, coding regions, introns and alternative splicing events, but also greatly broadens the understanding of the complexity of single-cell gene expression, providing a new path for biological research and clinical diagnosis.

[0005] However, ONT single-cell full-length sequencing technology faces challenges in sample preparation, data processing, and cost-effectiveness. These challenges involve maintaining the integrity of RNA molecules, improving the quality of sequencing data, and reducing operating costs. Therefore, developing new optimization methods to address these challenges is of vital importance to expanding the application of ONT single-cell full-length sequencing technology in scientific research and diagnosis. Summary of the invention

[0006] The present invention provides a sequencing analysis method for a single-cell transcriptome based on ONT sequencing. The method utilizes 10×Genomics single-cell capture, labeling technology and ONT sequencing technology to more accurately and comprehensively analyze transcripts and genes in cells, and deeply mine transcription information in the life process, thereby providing more valuable data and applications for fields such as life science research, medical diagnosis and drug development. The sequencing analysis method for a single-cell transcriptome provided by the present invention is specifically implemented by the following technologies.

[0007] A single-cell transcriptome sequencing analysis method based on ONT sequencing comprises the following steps:

[0008] Preparing and constructing a cDNA library of the sample to be tested, sequencing it, and obtaining the third-generation original data of the sample to be tested; performing quality control on the third-generation original data to obtain the third-generation valid data of the sample to be tested;

[0009] Functional annotation of genes and transcripts of the reference genome;

[0010] Combining the three generations of effective data with the gene data of the reference genome, the software ONT-SC-pipeline was used to perform data filtering, data splitting, comparison and quantification to obtain quantitative data;

[0011] Performing double cell identification and removal on the quantitative data, and further filtering the data;

[0012] The filtered data were subjected to basic analysis, marker gene enrichment analysis and advanced analysis in turn; the basic analysis included top gene distribution analysis, highly variable gene analysis, clustering analysis, correlation analysis of various cell groups, dimensionality reduction analysis and marker gene identification; the advanced analysis included cell annotation analysis, cell communication analysis, cell trajectory analysis, variable splicing analysis and fusion gene analysis.

[0013] Furthermore, the quality control processing method is: screening according to the sequencing quality value, and eliminating sequences with a length of ≤30bp and a quality value of ≤5.

[0014] Furthermore, the method of using the software ONT-SC-pipeline for data segmentation, comparison and quantification is as follows: using vsearch for data segmentation, and using minimap2 for data comparison, extracting the information of the segmentation result and the comparison result and merging them to obtain the quantitative data.

[0015] Furthermore, scDblFinder software was used to identify and remove double cells and low-activity cells, and cells with a high proportion of ribosomal RNA and a high proportion of mitochondrial genes were eliminated.

[0016] Furthermore, SeuratV5 software was used to analyze the distribution of top genes; the function FindvariableFeatures was used to analyze highly variable genes; and SeuratV5 software was used for clustering analysis, dimensionality reduction analysis, and marker gene identification.

[0017] Furthermore, the dimensionality reduction analysis is specifically performed using at least one of the PCA, T-SNE, and UMAP analysis methods in SeuratV5 software.

[0018] Specifically, if the sample to be tested is a multiple sample, the harmony software is used instead of the original PCA analysis method to correct and integrate the filtered multiple sample single cell data.

[0019] Furthermore, the hypergeometric distribution algorithm was used to perform enrichment analysis of marker genes.

[0020] Furthermore, known marker genes and scina software were combined for cell annotation analysis.

[0021] Furthermore, CellPhoneDB V5 was used for cell communication analysis.

[0022] Furthermore, the software Monocle3 was used for cell trajectory analysis.

[0023] Furthermore, isoswitch software was used for variable splicing analysis.

[0024] Furthermore, the software JAFFA was used for fusion gene analysis.

[0025] Furthermore, monocle3 software was used for cell trajectory analysis.

[0026] Furthermore, isoswitch software was used for alternative splicing analysis.

[0027] Furthermore, fusion gene analysis was performed using JAFFA software.

[0028] Furthermore, based on the CellPhoneDB V5 results, ggplot2 and igraph methods were used to visualize the communication results of a single pathway between various cell subpopulations; the contribution of each subunit to the overall ligand-receptor pair and the subunits of a single ligand-receptor pair in the pathway were visualized.

[0029] Compared with the prior art, the present invention is beneficial in that:

[0030] 1. Compared with the second-generation sequencing method, ONT sequencing does not need to interrupt the fragments, can detect the complete transcripts, and can quantify the transcripts.

[0031] 2. Compared with the existing Pacbio third-generation sequencing, it has more cost advantages and has extremely high cost-effectiveness while taking into account multiple sequencing results. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] Figure 1 A flow chart of the sequencing analysis method of single-cell transcriptome based on ONT sequencing provided by the present invention.

[0033] Figure 2Flowchart for data filtering, data splitting, alignment, and quantification.

[0034] Figure 3 This is a correlation diagram of the number of UMIs before filtering in a specific implementation method: the X-axis is the number of UMIs, the Y-axis is the number of genes, and the numbers are Pearson correlation coefficients. The ordinate of the left figure is the mitochondrial ratio, and the ordinate of the right figure is the number of genes.

[0035] Figure 4 It is a statistical diagram before and after filtering in a specific implementation method.

[0036] Figure 5 Box plot of TOP genes.

[0037] Figure 6 This is a map of highly variable gene points.

[0038] Figure 7 This is the UMAP analysis result diagram.

[0039] Figure 8 This is the proportion of each Cluster.

[0040] Fig. 9 This is the correlation diagram for each category.

[0041] Fig.10 It is the heat map of Marker genes.

[0042] Fig.11 It is the UMAP diagram of Marker gene expression.

[0043] Fig.12 This is the bubble diagram of GO enrichment of Marker genes.

[0044] Fig.13 Annotate statistical graphs for cells.

[0045] Fig.14 Annotate UMAP images for cells.

[0046] Fig.15 This is a statistical chart of the number of cell communications.

[0047] Fig.16 This is the trajectory diagram after selecting root.

[0048] Fig.17 This is the top gene trend chart.

[0049] Fig.18 This is an example of variable shearing.

[0050] Fig.19 This is an example diagram of a fusion gene. DETAILED DESCRIPTION

[0051] The technical solution of the present invention will be described clearly and completely below. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of them. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0052] The sequencing analysis method of single-cell transcriptome based on ONT sequencing provided by the present invention has the following overall sequencing analysis process: Figure 1 The overall steps of this method are as follows:

[0053] Preparing and constructing a cDNA library of the sample to be tested, sequencing it, and obtaining the third-generation original data of the sample to be tested; performing quality control on the third-generation original data to obtain the third-generation valid data of the sample to be tested;

[0054] Functional annotation of genes and transcripts of the reference genome;

[0055] Combining the three generations of effective data with the gene data of the reference genome, the software ONT-SC-pipeline was used to perform data filtering, data splitting, comparison and quantification to obtain quantitative data;

[0056] Performing double cell identification and removal on the quantitative data, and further filtering the data;

[0057] The filtered data were subjected to basic analysis, marker gene enrichment analysis and advanced analysis in turn; the basic analysis included top gene distribution analysis, highly variable gene analysis, clustering analysis, correlation analysis of various cell groups, dimensionality reduction analysis and marker gene identification; the advanced analysis included cell annotation analysis, cell communication analysis, cell trajectory analysis, variable splicing analysis and fusion gene analysis.

[0058] Specifically, the steps of this method are as follows.

[0059] 1. Test Preparation

[0060] Experimental Materials: Experimental Materials: Mouse liver tissue.

[0061] Related reagents, kits and consumables: all purchased from currently commercially available products.

[0062] 2. Preparation of single cell suspension

[0063] 1. Dissect the mouse, cut off part of the liver tissue, wash the blood in the tissue with PBS or saline, put the washed tissue into a 1mL preservation tube filled with tissue preservation solution, and transport it to the laboratory at 4℃;

[0064] 2. Add 1× PBS solution to a 6-well plate, remove the tissue from the tissue preservation solution and wash it in PBS to remove residual blood and mucus; take an appropriate amount of the washed sample and add it to the dissociation solution, and cut the tissue into 1 mm pieces with surgical scissors. 3 left and right tissue blocks;

[0065] 3. Place the dissociation solution with tissue in a constant temperature metal bath and incubate at 37°C for 15-40 minutes. When the tissue is completely dissociated, mix the dissociation solution with a wide-mouthed pipette tip every 10 minutes. If there are large pieces of tissue, cut them into pieces again with scissors.

[0066] 4. After the dissociation is complete, gently mix the dissociation solution and filter it using a 40-70μm sterile cell sieve;

[0067] 5. Rinse the cell sieve with 5 mL of 1640 medium + 2% FBS;

[0068] 6. Mix the filtered cell suspension and centrifuge at 300g for 5 min at 4°C; discard the supernatant;

[0069] 7. Take 1mL 1640 medium + 2% FBS to resuspend the cell pellet; after mixing, take 10μL of the suspension for AO / PI staining and count and monitor the cells;

[0070] 8. Store the cells on ice or at 4°C until the next step of the experiment.

[0071] 3. Obtain and amplify cell cDNA

[0072] 1. Prepare 10×Genomics capture and labeling cell reaction master mix;

[0073] 2. Place PCR tubes on ice and add Master Mix and other required components to each reaction.

[0074] 3. Process the chip according to the official instructions of 10×Genomics;

[0075] 4. Load the single cell suspension onto the 10×Genomics chip;

[0076] 5. Run 10×Genomics Chromium Controller to prepare Gel Beads-in-emulsion (GEMs);

[0077] 6. Transfer the prepared GEMs to a PCR tube and place it in a PCR instrument for reverse transcription according to the reaction program shown in Table 1 below. Set the temperature of the heated cover to 53°C.

[0078] Table 1

[0079] program temperature time 1 53℃ 45min 2 85℃ 5min 3 4℃ hold

[0080] 7. After reverse transcription is completed, add Dynabeads Cleanup Mix to purify the reverse transcription product, and the first chain of cDNA is obtained after purification;

[0081] 8. Amplify the first strand of cDNA, prepare the amplification reaction system shown in Table 2 below, and amplify according to the reaction procedure shown in Table 3 below.

[0082] Table 2

[0083] Components Volume (μl) cDNA 35 Amp Mix 50 cDNA Primers 15 <![CDATA[H2O]]> 0

[0084] Table 3

[0085]

[0086] 9. Use the SPRIselect reagent kit to purify the amplified product to obtain a double-stranded cDNA amplification product.

[0087] 4. Construction of cDNA library and sequencing based on ONT sequencing platform

[0088] 1. Take 500 ng of amplified cDNA for end repair, and the amplification reaction system is shown in Table 4. Incubation conditions: 20°C for 30 min, 65°C for 20 min, and store at 4°C.

[0089] Table 4

[0090] Components volume DNA 48μl NEBNext FFPE buffer (NEB) 3.5μl NEBNext End-prep buffer (NEB) 2μl FFPE DNA repair Mix (NEB) 2μl End-prep enzyme Mix(NEB) 3.5μl Total 60μl

[0091] 2. Use 1 volume of magnetic beads to purify the repair product, 61 μl of elution buffer to elute the DNA, and take 1 μl of the purified product for quantification using Qubit;

[0092] 3. The sequencing adapter ligation system is shown in Table 5 below; incubation conditions: 25°C, 60 min.

[0093] Table 5

[0094] Components volume Repaired DNA 60μl LNB(SQK-LSK114) 25μl Ligation Adapter Mix(SQK-LSK114) 5μl Quick T4 DNA ligase (NEB) 10μl Total 100μl

[0095] 4. Use 0.4 times volume of magnetic beads to purify DNA, elute DNA with 25ul Elution buffer (SQK-LSK114), take 1ul library, and use Qubit for quantification;

[0096] 5. The library system was configured as shown in Table 6 below, and a cDNA library was constructed.

[0097] Table 6

[0098] Components volume library 24μl Library Buffer (SQK-LSK114) 51μl Squencing Buffer(SQK-LSK114) 75μl Total 150μl

[0099] 6. Load the library into the R10.4.1 sequencing chip and sequence it on the PromethION sequencer (Oxford Nanopore Technologies, Oxford, UK) for 72 hours to obtain the third-generation raw data.

[0100] 5. Quality Control Processing

[0101] In the above step 4, the third generation raw data of the Nanopore sequencing data is used, and the format is pod5 format containing all the raw sequencing signals. After base calling by Dorado software (version: 0.7.2), the pod5 format data is converted into fastq format to obtain the third generation valid data.

[0102] 6. Data Segmentation, Comparison and Quantification

[0103] This embodiment uses the software ONT-SC-pipeline for data processing, including data splitting, comparison and quantification. The flow chart is as follows: Figure 2 Specifically:

[0104] Use the vsearch function to extract the adapter and obtain the original data, filter the low-quality reads (i.e., remove the reads that do not contain the adapter), complete the data splitting, and obtain the clean data (splitting result).

[0105] Combined with the reference genome, minimap2 was used to align the clean data, and sringtie was used to identify new transcripts, align transcripts, and extract the relationships between reads ID, transcripts, and genes. At the same time, vsearch was used to identify the Barcode and UMI of the clean data, cluster Barcode and UMI, reallocate Barcode and UMI, and extract the relationships between reads ID, Barcode, and UMI.

[0106] The relationships among reads ID, gene and transcript, Barcode, and UNI were combined and calculated to obtain quantitative data, i.e., gene and transcript expression levels.

[0107] The software ONT-SC-pipeline specifically uses the nextflow workflow management system, integrating vsearch, seqkit, minimap and other methods to make them work together. Compared with traditional shell or python series scripts, this implementation solves the problems of multi-threading difficulties and slow running speed, and realizes the rapid realization from fastq raw data to gene expression. At the same time, it has good scalability and can be well linked with other processes.

[0108] The basic information of the full-length sequence data after quality control filtering is shown in Table 7. The alignment results and quantitative results are shown in Tables 8 and 9, respectively.

[0109] Table 7 Statistics of full-length sequences after quality control filtering (excerpt)

[0110] sample n_reads rl_mean n_fl fl_rate(%) A1 508909 982 465435 91.45

[0111] In Table 7, sample: sample name; n_reads: total number of reads; rl_mean: average length; n_fl: number of full-length reads; fl_rate: ratio of full-length reads.

[0112] Table 8 Comparison results (excerpt)

[0113] Sample Total reads mapped reads map_rate unmapped reads unmapped_rate A2 114,234 114,234 100.00％ 0 0.00％

[0114] In Table 8, sample: sample name; Total reads: total reads; mapped reads: mapped reads; map_rate: mapping rate; unmapped reads: unmapped reads; unmapped_rate: unmapped rate.

[0115] Table 9 Quantitative results table (excerpt)

[0116]

[0117] In Table 9, the first column is the gene name and the other columns are the cell numbers.

[0118] Among the captured cells, some cells have low activity or are even dead cells. Generally, cell filtering is performed based on the number of genes expressed by the cells. The number of genes expressed by different sample types is inconsistent, and different thresholds are used for specific projects, with 200 or 500 thresholds being common. The mitochondrial gene filtering threshold is generally set at 0.1 and 0.25. Figure 3 As shown, the X-axis is the number of UMIs, the Y-axis is the number of genes, and the numbers are the Pearson correlation coefficients. The statistical results of the filtered cells are shown in Table 10 below.

[0119] Table 10 Cell statistics after filtration (excerpt)

[0120] cell_ID orig.ident nCount_RNA nFeature_RNA scDblFinder.class A1_AAACCCAAGACTTCAC-1 A1 1274 881 singlet A1_AAACCCAAGATACAGT-1 A1 4771 2373 singlet A1_AAACCCAAGGGCAACT-1 A1 2418 1571 singlet A1_AAACCCAAGTCGCCCA-1 A1 6593 3228 singlet A1_AAACCCAGTCATGCAT-1 A1 1726 980 singlet

[0121] In Table 10, ell_ID: cell ID; orig.ident: sample name; nCount_RNA: total number of RNAs; nFeature_RNA: number of RNA types; scDblFinder.class: number of droplet cells; percent.mito: mitochondrial proportion; percent.ribo: rRNA proportion; dropouts.

[0122] like Figure 4 As shown, left panel: number of genes detected in single cells. Middle panel: number of UMIs detected in single cells. Right panel: percentage of mitochondrial gene expression.

[0123] 7. Fundamental Analysis

[0124] 1. Top gene distribution

[0125] Extract and plot the top N genes with the highest gene expression count in the cells. Figure 5 As shown in the figure, the Y axis is the gene name and the X axis is the ratio of gene counts. The larger the ratio, the higher the gene expression.

[0126] 2. Highly variable genes

[0127] The expression matrix of the remaining cells after quality control was normalized; then, the function FindvariableFeatures was used to identify genes with high expression changes.

[0128] like Figure 6 As shown, red represents some of the genes with the most obvious changes, and the gene name is marked on the right.

[0129] 3. Clustering and grouping analysis

[0130] In the face of numerous cells, this embodiment clusters the cells according to similarity, and then analyzes each cluster after clustering. Figure 7 As shown in Table 11, the horizontal axis is the category (Cluster) obtained by clustering, and the vertical axis is the proportion of the Cluster cells in each sample. The statistical results of cluster information are shown in Table 11.

[0131] Table 11 Cluster information statistics (excerpt)

[0132] Cluster number stat(%) 0 765 17.14 1 693 15.53 2 616 13.80 3 445 9.97 4 390 8.74

[0133] 4. Correlation analysis of various cell populations

[0134] The first 1000 cells with the largest standard deviation were selected for correlation analysis of each Cluster. Figure 8 As shown, different colors indicate the size of the correlation.

[0135] 5. Dimensionality reduction analysis

[0136] A notable feature of high-throughput single-cell sequencing data is that the data volume is large and a large number of cells can be reflected at one time. Therefore, it is very important to display the characteristics of cell data through dimensionality reduction and visualization.

[0137] In this embodiment, the single-cell dimensionality reduction graphs generally use t-SNE (t-distributed neighborhood embedding algorithm) and UMAP (Uniform Manifold Approximation and Projection).

[0138] (1) PCA

[0139] In transcriptome research, the relationship between samples is determined by many variables (gene expression). Principal component analysis is an attempt to regroup many original indicators (gene expression) with certain correlation into a new set of independent comprehensive indicators, thereby reducing the complexity of the problem and studying the relationship between gene expression and samples.

[0140] In the two-dimensional PCA analysis results, a scatter plot will be displayed with principal component 1 (PC1) and principal component 2 (PC2) as the X-axis and Y-axis respectively, and each point represents a cell. In the figure, if the distance between two cells is farther, it means that the expression pattern of the genes of the two cells is more different. Conversely, it means that the overall expression pattern of the corresponding cells is closer. Therefore, PCA analysis is often used to evaluate the quality of cell reproducibility.

[0141] If the samples are the result of processing from different batches, there will be large differences even for the same sample. In this case, we will use the harmony software package to perform batch correction. In this case, the displayed graph is the corrected result.

[0142] (2) T-SNE

[0143] t-SNE (t-distributed Stochastic Neighbor Embedding) is a nonlinear dimensionality reduction technique that calculates the distances between all nodes and maps high-dimensional data to low-dimensional space through probability distribution.

[0144] (3)UMAP

[0145] UMAP analysis, full name: Uniform Manifold Approximation and Projection. UMAP analysis is a new dimensionality reduction algorithm. When calculating spatial distribution, it calculates the distance between K adjacent points. It can ensure local information while keeping global information undistorted. Therefore, it has been increasingly used in the literature in recent years.

[0146] like Fig. 9 As shown, one point represents one cell, and different colors represent different Clusters.

[0147] 6. Marker gene analysis

[0148] (1) Marker gene identification analysis

[0149] The identification of marker genes is still based on the R package Seurat, and the function used is FindAllMarkers. In this embodiment, marker genes are identified for each cell cluster. The identification is to perform marker analysis on the target cell cluster and all remaining cells.

[0150] Marker analysis row number FindAllMarkers uses the wilcox statistical method. The marker gene standard is: logfc>0.5, pvalue≤0.05. Marker genes are generally selected from genes with log2fc in the cell cluster. Based on these marker genes, downstream visualization, cell annotation and other analyses can be performed.

[0151] like Fig.10 As shown: the horizontal axis is the gene with the highest avg_log2FC in each Cluster, the vertical axis is different Clusters, and the color of the point represents the normalized expression level.

[0152] like Fig.11 The color of the point indicates the expression level, and the bottom map is the UMAP map.

[0153] (2) Marker gene enrichment analysis

[0154] Enrichment analysis refers to the use of statistical methods to study whether the differentially expressed gene functions are concentrated in certain specific functional categories. Commonly used functional categories include Gene Ontology and KEGG Pathway. In this embodiment, enrichment analysis uses a hypergeometric distribution algorithm to analyze and find GOTerm or KEGG Pathway that is significantly enriched in differentially expressed genes relative to all annotated genes.

[0155] The differentially expressed gene enrichment analysis software is ClusterProfiler (version: 3.14.3). p.adjust(qvalue) is the pvalue after multiple hypothesis testing correction. The value range of p.adjust(qvalue) is [0,1]. The closer it is to zero, the more significant the enrichment.

[0156] like Fig.12 As shown, it is a GO enrichment bubble chart, the horizontal axis is RichFactor (the proportion of Marker genes in this GO to all Marker genes), the vertical axis is Gene Ontology, and the color of the point represents p.adjust (qvalue).

[0157] The selected excerpts of GO enrichment examples are shown in Table 12.

[0158] Table 12 GO enrichment example table (excerpt)

[0159]

[0160] ID: Gene Ontology number; Description: Description of Gene Ontology function; GeneRatio: The ratio of the marker gene annotated to this GO Term to all the marker genes with GO annotation; BgRatio: The ratio of all genes annotated to this GO Term to all genes with GO functional annotation; pvalue: Statistical significance level indicator.

[0161] Generally, pvalue < 0.05 indicates significant enrichment; p.adjust: adjusted pvalue; qvalue: adjusted pvalue, which is different from the p.adjust correction method.

[0162] Each term in the GO database has a clear inclusion and substitution relationship. Based on these relationships, we can draw a directed acyclic graph to view the relationship between each term more clearly.

[0163] However, since there are a large number of terms, it would be very complicated to draw a complete directed acyclic graph, so this embodiment only displays part of the terms. For example, an ellipse or a square represents a GO term, and arrows represent the relationship between different GO terms.

[0164] 7. Cell Annotation

[0165] Single-cell transcriptome sequencing (scRNA-seq) data analysis can identify different cell types through differences in gene expression levels between different cells, thereby studying cell heterogeneity in complex tissues, and ultimately helping us understand the pathogenesis of diseases, cell lineages and differentiation trajectories, or intercellular communication mechanisms. In the actual scRNA-seq data processing process, cell type annotation is a key step in subsequent analysis and research.

[0166] In this embodiment, the annotation of cell types is divided into cell-based annotation and cluster-based annotation according to the annotation object, and can be divided into supervised annotation and unsupervised annotation according to whether label data is required. The statistical results of cell annotation are shown in Table 13 below.

[0167] Table 13 Cell annotation statistics (excerpt)

[0168] type count lateral_root 3330 xylem 1782 branch 1142 lateral_root_cap 976 phloem 921

[0169] like Fig.13 As shown, each subpopulation is stacked in a bar chart, with the horizontal axis representing different samples, the vertical axis representing the number of cells, and different colors representing different cell types.

[0170] like Fig.14 As shown, the figure is a UMAP map, and different colors represent different cell types.

[0171] 8. Cell communication analysis

[0172] Cell communication analysis is receptor-ligand analysis. Cell interactions mediated by biological cell ligands-receptors regulate various bioinformatics processes such as development, differentiation, and inflammation. The mutual communication between cells constitutes a cell regulatory network, which affects the progress of various physiological processes. The tissue ecosystem (microenvironment), especially the tumor microenvironment, is composed of various types of cells, including immune cells, stromal cells, hematopoietic cells, etc. These different types of cells can play an important role in tumor development, drug resistance, immune infiltration, and inflammation through interactions between ligands and receptors (for example, with immune checkpoint inhibitors).

[0173] This embodiment uses the software CellPhoneDB V5 for communication analysis. The software uses network analysis and pairing mode to identify the main input and output signals between cells and determine the communication relationship between them. It also provides a variety of visualization methods for customer convenience.

[0174] In this embodiment, the software CellPhoneDB V5 uses the gene expression data of cells as input and combines the communication between ligands, receptors and their auxiliary factors to simulate intercellular communication. In order to establish intercellular communication, "cell_type" is used as a label for communication analysis.

[0175] (1) (Using the communication analysis R package CellPhoneDB V5) Statistical analysis of cell communication quantity

[0176] The statistical results (visualization) of all cell subpopulations are: the overall statistical results of the number of receptor-ligand pairs and communication strength between different cell subpopulations. The size of the circles of various colors on the periphery indicates the number of cells. The larger the circle, the more cells there are. The cells emitting arrows express ligands, and the cells pointed by the arrows express receptors. The more ligand-receptor pairs there are, the thicker the line.

[0177] (2) Statistical analysis of cell communication intensity

[0178] like Fig.15 As shown in the visualization statistical analysis results, the size of the peripheral circles of various colors indicates the number of cells. The larger the circle, the more cells there are. The cells emitting arrows express ligands, and the cells pointed by the arrows express receptors. The greater the probability of ligand-receptor pair communication, the thicker the line.

[0179] (3) Single-channel communication visualization

[0180] Multiple ligand-receptor pairs may be on the same pathway. The communication between the ligand-receptor pairs involved in the pathway is analyzed in units of pathways. The communication between each cell subpopulation in a single pathway is visualized. Heat maps, chord diagrams, and network diagrams are just different forms of expression, but the meaning is the same. In the visualization results, the ordinate represents the ligand, the abscissa represents the receptor, and the color represents the communication probability (Communication Prob). The two ends of the line represent the ligand-receptor, and the width of the line represents the communication probability (Communication Prob).

[0181] (4) Ligand-receptor visualization

[0182] There are many ligand receptors in the pathway, and the ligand receptors may also be composed of multiple subunits, and the influence of each subunit may be large or small. Therefore, the contribution of each subunit to the overall ligand receptor and the subunits of a single ligand receptor pair in the pathway can be visualized, with the horizontal axis being the relative contribution and the vertical axis being the different subunits; the influence of each subunit on different cell types can be displayed; multiple ligand-receptor-mediated cell communication bubble diagrams can also be displayed, with the horizontal axis being the communication between cell subpopulations and the vertical axis being the ligand receptor. The point size is the P value, and the color indicates the probability of communication.

[0183] (5) Related gene expression dot plots and violin plots

[0184] The expression levels of pathway ligand receptor-related genes in different cell types can also be visualized.

[0185] 9. Cell trajectory analysis (quasi-time series analysis)

[0186] Pseudotime is a method for measuring the degree of progress of a cell in a biological process such as cell differentiation. It is an abstraction of a biological process and describes the shortest path distance between the current state of a cell and the starting state of the time trajectory. This embodiment uses the software monocle3 to analyze cell trajectories. Fig.16 As shown, the trajectory diagram after selecting the root node.

[0187] (1) Differential genes in pseudo-time trajectory

[0188] Moran's index (morans_I) can represent spatial correlation, and its value is between -1 and 1, where 0 means that the gene has no spatial correlation. In a single cell, it can represent the spatial co-expression effect, and 1 means that the expression values ​​of this gene in cells with similar spatial distances are highly similar.

[0189] like Fig.17 As shown in the figure, the time series diagram of some genes in different Clusters, the horizontal axis is the different time series, and the vertical axis is the expression level.

[0190] (2) Analysis of co-expression modules

[0191] Once you have a set of genes that vary in some way across clusters, you can group them into modules, and the software monocle3 provides exactly this method.

[0192] 10. Variable shear analysis (advanced analysis)

[0193] Thanks to the full-length ONT being able to measure complete transcripts, we were able to analyze the variable splicing of genes through different transcripts of the same gene. We used the software isoswitch for variable splicing analysis and used the plot_assay_stats() function to draw a summary graph, including the number of genes, the number of transcripts, the distribution of isoforms, and the number of genes in each cell type, which can concisely describe the distribution of isoforms in the dataset.

[0194] like Fig.18 The following is an example of variable splicing. Top left: Gene (single transcript, multiple transcripts) and transcript number statistical bar graph; Middle left: Gene corresponding transcript number statistical bar graph; Bottom left: Gene transcript number vs. main transcript number scatter plot; Right: Different Cluster expression gene box plot.

[0195] 11. Fusion gene analysis

[0196] Full-length transcript sequencing technology provides a powerful tool for revealing structural information at the transcript level, and also opens the door to detecting gene fusion events at the single-cell level. Gene fusion, that is, the partial or complete sequence recombination of two or more genes to form a new hybrid gene, can be caused by a variety of mechanisms such as chromosomal translocation, deletion or inversion, and plays a vital role in the occurrence and development of diseases.

[0197] This embodiment uses the software JAFFA for fusion gene analysis. By aligning the clean reads obtained by sequencing to the reference genome and combining the gene annotation information, we can extract sequence fragments spanning two genes and then identify potential fusion genes. To ensure the reliability of the results, we set a coverage depth threshold (default: cutoff ≥ 5) to filter low-quality data, and finally obtained the following fusion gene information, as shown in Table 14 and Fig.19 shown.

[0198] Table 14 (excerpt)

[0199]

[0200] The above specific embodiments describe the implementation of the present invention in detail, but the present invention is not limited to the specific details in the above embodiments. Within the scope of the claims and technical concept of the present invention, the technical solution of the present invention can be modified and changed in many simple ways, and these simple modifications all belong to the protection scope of the present invention.

Claims

1. A single-cell transcriptome sequencing analysis method based on ONT sequencing, characterized in that: The following steps are involved: Preparing and constructing a cDNA library of the sample to be tested, sequencing it, and obtaining the third-generation original data of the sample to be tested; Performing quality control on the three generations of raw data to obtain three generations of valid data of the sample to be tested; Functional annotation of genes and transcripts of the reference genome; Combining the three generations of effective data with the gene data of the reference genome, the software ONT-SC-pipeline was used to perform data splitting, alignment and quantification to obtain quantitative data; Performing double cell identification and removal on the quantitative data, and further filtering the data; The filtered data were subjected to basic analysis, marker gene enrichment analysis and advanced analysis in turn; the basic analysis included top gene distribution analysis, highly variable gene analysis, clustering analysis, correlation analysis of various cell groups, dimensionality reduction analysis and marker gene identification; the advanced analysis included cell annotation analysis, cell communication analysis, cell trajectory analysis, variable splicing analysis and fusion gene analysis.

2. The single-cell transcriptome sequencing analysis method based on ONT sequencing according to claim 1, characterized in that: The quality control method is as follows: screening is performed according to the sequencing quality value, and sequences with a length of ≤30 bp and a quality value of ≤5 are eliminated.

3. The single-cell transcriptome sequencing analysis method based on ONT sequencing according to claim 1, characterized in that: The method of using the software ONT-SC-pipeline for data segmentation, comparison and quantification is as follows: using vsearch for data segmentation, and using minimap2 for data comparison, extracting the information of the segmentation result and the comparison result, and merging them to obtain the quantitative data.

4. The single-cell transcriptome sequencing analysis method based on ONT sequencing according to claim 1, characterized in that: scDblFinder software was used to identify and remove double cells and low-activity cells, and cells with a high proportion of ribosomal RNA and a high proportion of mitochondrial genes were eliminated.

5. The single-cell transcriptome sequencing analysis method based on ONT sequencing according to claim 1, characterized in that: SeuratV5 software was used to analyze the distribution of top genes; the FindvariableFeatures function was used to analyze highly variable genes; SeuratV5 software was used for clustering analysis, dimensionality reduction analysis and Marker gene identification.

6. The single-cell transcriptome sequencing analysis method based on ONT sequencing according to claim 5, characterized in that: The dimensionality reduction analysis is specifically performed using at least one of the PCA, T-SNE, and UMAP analysis methods in the SeuratV5 software.

7. The single-cell transcriptome sequencing analysis method based on ONT sequencing according to claim 6, characterized in that: If the sample to be tested is a multiple sample, the harmony software is used instead of the PCA analysis method to correct and integrate the filtered multiple sample single cell data.

8. The single-cell transcriptome sequencing analysis method based on ONT sequencing according to claim 1, characterized in that: The hypergeometric distribution algorithm was used for enrichment analysis of marker genes.

9. The single-cell transcriptome sequencing analysis method based on ONT sequencing according to claim 1, characterized in that: Combine known marker genes and scina software for cell annotation analysis; use CellPhoneDB V5 for cell communication analysis; use monocle3 software for cell trajectory analysis; use isoswitch software for variable splicing analysis; use JAFFA software for fusion gene analysis; use monocle3 software for cell trajectory analysis; use isoswitch software for variable splicing analysis; use JAFFA software for fusion gene analysis.

10. The single-cell transcriptome sequencing analysis method based on ONT sequencing according to claim 9, characterized in that: Based on the results of CellPhoneDB V5, ggplot2 and igraph methods were used to visualize the communication results of a single pathway between various cell subpopulations; the contribution of each subunit to the overall ligand-receptor relationship and the subunits of a single ligand-receptor pair in the pathway were visualized.

Images