An analysis method for hbv expression level and fragment distribution based on single cell rna-seq

Abstract

The application relates to the technical field of biological industry, and discloses an analysis method for HBV expression level and fragment distribution based on single-cell RNA-seq, constructs position information of each gene of HBV in a linear transcriptome and a circular genome, and obtains a mapping relationship of the transcriptome and the genome; compares original sequencing data obtained by scRNA-seq to an HBV transcriptome reference sequence, identifies human sequences and HBV sequences, removes the human sequences, and obtains HBV source transcriptome sequences; converts the position of the sequencing sequence on the transcriptome into a comparison position on the genome through the mapping relationship of the transcriptome and the genome; calculates the number of sequencing sequences covered on each base of the genome, draws the coverage depth of each base of the HBV genome, and obtains the expression level and fragment distribution of single-cell RNA-seq sequencing data on the HBV genome. The application fills the blank of the HBV analysis method based on single-cell resolution.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of biological industry, and particularly relates to an analysis method of HBV expression level and fragment distribution based on single-cell RNA-seq. BACKGROUND

[0002] Gene expression refers to the process of transcription of deoxyribonucleic acid (DNA) into messenger ribonucleic acid (mRNA) and translation of RNA into protein. This process is a key step in realizing genetic instructions in living organisms. Digital gene expression profiling (DGE) and transcriptome analysis (RNA-Seq) are new methods that use next-generation high-throughput sequencing technology and high-performance computing analysis technology to sequence and accurately analyze the gene expression of a specific tissue and state of a species.

[0003] Single-cell RNA-seq (scRNAseq) technology is a high-throughput sequencing analysis technology for the transcriptome level of a single cell, which can detect the gene expression level in each cell. This technology mainly captures the 3' region of RNA for library sequencing, and the data after sequencing mainly covers the 3' or 5' region of the gene.

[0004] Single-cell RNA sequencing (scRNA-seq) technology provides unprecedented resolution and insight for biological research. Here are some main advantages of scRNA-seq technology:

[0005] 1. High resolution: scRNA-seq can analyze gene expression at the level of individual cells, which allows researchers to identify subtle differences between cells and even discover new cell types or subpopulations.

[0006] 2. Revealing heterogeneity: Traditional bulk cell population sequencing methods only provide averaged data, masking significant differences that may exist between cells. scRNA-seq can detail the heterogeneity within a cell population, which is particularly important for understanding complex tissue structure and function.

[0007] 3. Exploring dynamic changes: By sampling the same type of cells at different time points, scRNA-seq can track changes in gene expression, helping to understand the dynamic changes in cell development, differentiation, and response to external stimuli.

[0008] 4. Low input requirement: scRNA-seq can be used to analyze even a small number of cells isolated from a small amount of sample, making it an ideal tool for studying rare cell types or samples that are difficult to obtain.

[0009] 5. Multi-dimensional information: In addition to gene expression levels, some advanced scRNA-seq methods can also provide other types of molecular information, such as variable splicing, allele-specific expression, and the function of non-coding RNA.

[0010] 6. Disease research and personalized medicine: scRNA-seq helps to understand the changes at the cellular level in disease states, providing evidence for the development of new diagnostic markers and treatment strategies, and promoting the development of personalized medicine.

[0011] 7. Drug response evaluation: This technology can be used to evaluate the response of individual cells to drugs, helping to predict drug efficacy and reduce side effects, thereby optimizing treatment plans.

[0012] However, there are significant limitations in using single-cell RNA sequencing (scRNA-seq) technology to analyze the expression level and fragment distribution of Hepatitis B virus (HBV). The main reasons are as follows:

[0013] 1. Integrated expression mechanism: After HBV infects cells, its genome can be integrated into the host human genome and expressed along with the transcription of host genes. This integrated expression mode means that there may be significant differences in the expression level and fragment distribution of HBV between different individuals and different cell types. In order to better understand the mechanism of action of HBV, it is necessary to accurately measure these differences. However, scRNA-seq technology is mainly used to detect mRNA expression in a non-integrated state, and is not good at identifying and quantifying specific changes caused by viral genome integration.

[0014] 2. Low abundance and background noise: Due to the relatively low content of HBV in infected cells, combined with the background noise generated by a large number of mRNAs in host cells, it is difficult to effectively capture and quantify HBV RNA using scRNA-seq. This makes it difficult to accurately assess the expression level of HBV.

[0015] 3. Lack of 3' polyA tail: HBV mRNA lacks a 3' polyA tail structure, which is a key problem because most scRNA-seq library construction methods rely on polyA tails to select mRNAs for sequencing. This means that HBV mRNA is not easily captured, resulting in a very low content of HBV-related sequences in the final sequencing data, further limiting accurate measurement of its expression level.

[0016] 4. Circular genome alignment challenge: The genome of HBV is circular, while existing bioinformatics tools are mainly designed for linear genomes for aligning and analyzing sequencing data. This leads to inaccurate or incomplete matches when trying to map data from HBV back to the reference genome, which can affect the understanding of HBV genome fragment distribution.

[0017] In view of this, there is currently no HBV analysis method directly derived from single-cell sequencing data. Existing HBV analysis techniques such as Northern Blotting, RT-PCR, RNA-seq (bulk), Reporter Assays, etc. can provide important information about HBV expression levels and the expression abundance of each region of the genome, but they are not based on single-cell resolution data. Although traditional bulk RNA-seq and other molecular biology techniques can provide an overall view of HBV expression, in some cases, especially when faced with a heterogeneous cell population, these methods may not be sufficient to reveal subtle differences in complex situations. SUMMARY

[0018] The present application aims to provide an analysis method of HBV expression level and fragment distribution based on single-cell RNA-seq to solve the problem of the lack of methods to calculate HBV expression level derived from single-cell sequencing data and locate the expression abundance of each region on the HBV genome, filling the gap in HBV analysis methods based on single-cell resolution.

[0019] To achieve the above-mentioned purpose, the technical scheme adopted by the present application is as follows:

[0020] An analysis method of HBV expression level and fragment distribution based on single-cell RNA-seq, comprising the following steps:

[0021] Step 1: Construct the position information of each gene of HBV in the linear transcriptome and circular genome to obtain the mapping relationship between the transcriptome and the genome;

[0022] Step 2: Align the raw sequencing data obtained by scRNAseq to the HBV transcriptome reference sequence to determine the breakpoint position when HBV integrates into the human genome, and identify human and HBV sequences;

[0023] Step 3: Remove human sequences, retain the remaining HBV sequences and the corresponding transcriptome coordinate positions starting from the potential breakpoint position of HBV to obtain HBV-derived transcriptome sequences;

[0024] Step four, the HBV source transcriptome sequence obtained in step three is mapped to the transcriptome and genome mapping relationship constructed in step one, the position of the sequencing sequence on the transcriptome is converted to the position of the sequencing sequence on the genome.

[0025] Step five, the number of sequencing sequences covering each base of the genome is calculated, the coverage depth of each base of the HBV genome is plotted, and the gene position is marked and displayed, to obtain the expression level and fragment distribution of single-cell RNA-seq sequencing data on the HBV genome.

[0026] Preferably, in step one, the HBV genome reference sequence and the HBV transcriptome reference sequence are obtained respectively, the annotation information of the HBV genome is obtained, the position information of each gene of the HBV in the transcriptome and the genome is constructed, and the mapping relationship between the transcriptome and the genome is obtained; the mapping relationship between the transcriptome and the genome is that the starting position 1376 and the ending position 1840 of NC_003977.2 of the HBV genome are X gene, which is X protein; the starting position 1816 and the ending position 2454 of NC_003977.2 are C gene, which is pre-capsid protein; the starting position 1903 and the ending position 2454 of NC_003977.2 are C gene, which is capsid protein; the starting position 2309 and the ending position 3182 of NC_003977.2 are P gene, which is polymerase; the starting position 1 and the ending position 1625 of NC_003977.2 are P gene, which is polymerase; the starting position 2850 and the ending position 3182 of NC_003977.2 are S gene, which is large envelope protein; the starting position 1 and the ending position 837 of NC_003977.2 are S gene, which is large envelope protein; the starting position 3174 and the ending position 3182 of NC_003977.2 are S gene, which is middle membrane protein; the starting position 1 and the ending position 837 of NC_003977.2 are S gene, which is middle membrane protein; the starting position 157 and the ending position 837 of NC_003977.2 are S gene, which is small envelope protein.

[0027] Preferably, in step two, the sequencing raw data FASTQ file obtained by scRNAseq is aligned to the transcriptome reference sequence of HBV by using bwa mem software, only the sequencing sequence information that can be aligned to the HBV genome is retained, and the BAM format alignment record file aligned to the HBV genome is generated as the HBV source sequence and the rest is human source sequence.

[0028] Preferably, in step three, according to the CIGAR record information of each sequencing sequence in the BAM file, the sequencing sequence with soft cut structure is extracted, and the soft cut part is cut off, and the remaining sequence and the corresponding transcriptome coordinate position are retained; the sequencing sequence without soft cut structure is retained.

[0029] Preferably, in the process of determining the breakpoint position when HBV integrates into human genome, the sequence alignment, soft cutting structure distribution and transcriptome coordinate position are comprehensively considered by the following ways:

[0030] For sequence alignment, the alignment results of sequencing sequences and HBV transcriptome reference sequences are analyzed in detail. If the base matching rate of sequencing sequences and reference sequences in a certain region is lower than the preset threshold, and the low matching rate shows continuous distribution, which is significantly different from the matching mode of normal alignment region, the region is marked as a potential breakpoint position region; the preset threshold is 70%-80%;

[0031] For the distribution of soft cutting structure, the sequencing sequences with soft cutting structure are statistically analyzed, the starting position and ending position of each soft cutting structure in the sequence are counted, and these positions are mapped to the transcriptome coordinates; when the appearance frequency of soft cutting structure in a specific transcriptome coordinate interval exceeds the preset frequency threshold, the coordinate interval is marked as a potential breakpoint position interval; the preset frequency threshold is that the number of soft cutting structure appearing in every 1000 sequencing sequences exceeds 10 times;

[0032] In combination with the transcriptome coordinate position, the transcriptome coordinate continuity of sequencing sequences is analyzed in the human-HBV genome chimeric region; if the transcriptome coordinate jumps or interrupts, or does not match the known human genome and HBV genome transcriptome coordinate range, the corresponding position is marked as a potential breakpoint position region;

[0033] Finally, the potential breakpoint position regions and potential breakpoint position intervals marked by the above three ways are comprehensively considered. When multiple marked regions overlap or densely distribute in a smaller range, the common region or dense region is determined as the breakpoint position when HBV integrates into human genome.

[0034] Preferably, the dense distribution in a smaller range refers to that, in a range with a center at a certain position and a radius not more than 50 bp, there are 3 or more potential breakpoint position regions marked by different ways.

[0035] Preferably, in the process of analyzing sequence alignment to determine potential breakpoint position region, the continuous distribution of low matching rate refers to that the number of continuous bases in the low matching rate region exceeds 10.

[0036] Preferably, in the process of counting soft cutting structure distribution to determine potential breakpoint position interval, the length distribution of soft cutting structure is also considered; when the soft cutting structure not only appears in a certain transcriptome coordinate interval with a frequency exceeding the preset frequency threshold, but also has an average length at least 10 bp longer than that in normal region and a maximum length at least 20 bp longer than that in normal region.

[0037] Preferably, in step four, the location of the sequencing sequence alignment on the transcriptome is converted into the alignment position on the genome as follows: gene NC_003977.2, start position 1712, end position 1819, nucleotide sequence AAAGACTGTGTGTTTAATGAGTGGGAGGAGTT.

[0038] Preferably, in step four, the location of the sequencing sequence alignment on the transcriptome is converted into the alignment position on the genome as follows: gene NC 003977.2, start position 1511, end position 1633, nucleotide sequence C0GACCGACCACGGGGCGCACCTCTCTTTACG.

[0039] The advantages of the present application are:

[0040] The present application solves the problem of the lack of methods for calculating the expression level of HBV from single-cell sequencing data and positioning the expression abundance of each region on the HBV genome, filling the gap in HBV analysis methods based on single-cell resolution.

[0041] 1. Improved accuracy: Since the HBV genome is a circular structure, and existing sequence alignment tools are mainly designed for linear genomes, this leads to potential quantitative bias when aligning the S and P genes that span the start and end positions. This solution first aligns the sequencing data to a linear transcriptome, and then indirectly aligns it back to the HBV genome based on the mapping relationship between the transcriptome and the genome, thereby enabling more accurate quantification of the expression levels of these key genes.

[0042] 2. Ability to distinguish overlapping genes and multiple transcripts: There are many overlapping genes and multiple transcripts of a gene in the HBV genome, which increases the difficulty of correctly analyzing their respective expression patterns. This solution uses the method of first aligning the data to the transcriptome, and then mapping it back to the genome position, which helps to more accurately distinguish the specific expression of different overlapping genes and different transcripts of the same gene.

[0043] 3. Provide additional information: Using sequencing sequences that are chimeric with human-HBV genome, not only can the expression of HBV itself be evaluated, but also can be used to identify the exact breakpoint position of HBV integrated into the human genome. This method is not limited to monitoring expression levels, but also provides important clues about the integration mechanism of the virus, which is of great significance for in-depth understanding of HBV infection, replication process and its carcinogenic potential.

[0044] The HBV expression level and fragment distribution analysis method based on single-cell RNA-seq of this invention has unique advantages. For example, it innovatively transforms the circular HBV genome into linear transcriptomic RNA for alignment, and then utilizes the mapping relationship between the transcriptome and genome to achieve accurate analysis, overcoming the inaccuracy problem caused by the circular structure of HBV in traditional methods. Because of this, it can be implemented through a single software program. The software integrates a complete workflow from sequence acquisition, data alignment, sequence identification, transcriptomic sequence processing, position transformation to expression level and fragment distribution calculation and display. Its optimized algorithms and multi-threading technology ensure efficient and accurate data processing. It also provides built-in visualization and auxiliary interpretation functions to facilitate researchers' understanding of the results. At the same time, the robust system architecture ensures data security and traceability, providing a one-stop, efficient, accurate, and reliable analysis platform for HBV research and promoting the development of related fields. Attached Figure Description

[0045] Figure 1 This is a schematic diagram illustrating the location information of each HBV gene in the transcriptome and genome in an embodiment of the present invention.

[0046] Figure 2 This is a schematic diagram illustrating the sequence and transcriptome alignment information extracted from HBV in scRNAseq sequencing data according to an embodiment of the present invention.

[0047] Figure 3 This is a schematic diagram showing the alignment positions of sequencing data on the genome in an embodiment of the present invention.

[0048] Figure 4 This is a distribution map of the expression level of sequencing data on the HBV genome in an embodiment of the present invention.

[0049] Figure 5 This is a complete flowchart of the HBV expression level and fragment distribution analysis method based on scRNA-seq according to an embodiment of the present invention. Detailed Implementation

[0050] The following detailed description illustrates the specific implementation method:

[0051] Example 1

[0052] This embodiment is basically as shown in the appendix. Figure 1 To be continued Figure 5 As shown: The method for analyzing HBV expression levels and fragment distribution based on single-cell RNA-seq in this embodiment includes the following steps:

[0053] The first step is to build the mapping relationship between the transcriptome and the genome: obtain the HBV genome reference sequence, the HBV transcriptome reference sequence, obtain the annotation information of the HBV genome, build the position information of each HBV gene in the transcriptome and the genome, and obtain the mapping relationship between the transcriptome and the genome; In this embodiment, the HBV genome sequence (https: / / www.ncbi.nlm.nih.gov / nuccore / NC_003977.2?report=fasta), the HBV transcriptome sequence (https: / / ftp.ncbi.nlm.nih.gov / genomes / all / GCF / 000 / 861 / 825 / GCF_000861825.2_ViralProj15428 / GCF_000861825.2_ViralProj15428_cds_from_genomic.fna.gz) and the HBV genome annotation file (https: / / ncbi.nlm.nih.gov / datasets / gene / GCF_000861825.2 / ) are downloaded from the National Center for Biotechnology Information (NCBI) of the United States, and the mapping relationship between the transcriptome and the genome is built according to the attached Figure 1 The position information of each HBV gene in the transcriptome and the genome is built, the X gene is from the start position 1376 to the end position 1840 of the HBV genome NC_003977.2, which is the X protein; the C gene is from the start position 1816 to the end position 2454 of the HBV genome NC_003977.2, which is the pre-capsid protein; the C gene is from the start position 1903 to the end position 2454 of the HBV genome NC_003977.2, which is the capsid protein; the P gene is from the start position 2309 to the end position 3182 of the HBV genome NC_003977.2, which is the polymerase; the P gene is from the start position 1 to the end position 1625 of the HBV genome NC_003977.2, which is the polymerase; the S gene is from the start position 2850 to the end position 3182 of the HBV genome NC_003977.2, which is the large envelope protein; the S gene is from the start position 1 to the end position 837 of the HBV genome NC_003977.2, which is the large envelope protein; the S gene is from the start position 3174 to the end position 3182 of the HBV genome NC_003977.2, which is the middle membrane protein; the S gene is from the start position 1 to the end position 837 of the HBV genome NC_003977.2, which is the middle membrane protein; the S gene is from the start position 157 to the end position 837 of the HBV genome NC_003977.2, which is the small envelope protein.

[0054] The second step is to align the sequencing data and determine the sequence source and breakpoint position: as shown in the attached Figure 2 The raw sequencing data obtained by scRNAseq (including the attached Figure 2The human-derived sequences shown by dashed lines and the HBV-derived sequences shown by solid lines are aligned to the HBV transcriptome reference sequence, the human-derived sequences are identified as human-derived sequences, and the HBV-derived sequencing data in the sequencing data is identified as HBV-derived sequences. Specifically, the sequencing raw data FASTQ file obtained by scRNAseq is aligned to the transcriptome reference sequence of HBV using the bwa mem software (version v0.7.17), only the sequencing sequence information that can be aligned to the HBV genome is retained, and a BAM format alignment record file aligned to the HBV genome is generated.

[0055] The third step is to extract and process the HBV-derived sequences: extract the HBV-derived sequences in the second step, and cut off the part of the sequence derived from the human transcriptome, i.e. the human-derived sequences that have been identified and separated, and retain the remaining HBV-derived sequences and the corresponding transcriptome coordinate positions starting from the HBV potential breakpoint position to obtain the HBV-derived transcriptome sequences. Specifically, as shown in Figure 2 According to the CIGAR record information of each sequencing sequence in the BAM file, the sequencing sequence with a soft clip (Soft clip, S) structure is extracted, and the soft clip part is cut off, and the remaining sequence and the corresponding transcriptome coordinate position are retained; the sequencing sequence without a soft clip structure is retained.

[0056] The fourth step is to convert the alignment position: the HBV-derived transcriptome sequence obtained in the third step is converted by the transcriptome and genome mapping relationship constructed in the first step to convert the sequencing sequence alignment position on the transcriptome to the alignment position on the genome; specifically, as shown in Figure 3 The sequence derived from HBV in the scRNAseq sequencing data and the transcriptome alignment position information are combined with the results of the first step, and the transcriptome alignment position of the sequencing data is converted to the genome alignment position:

[0057] Example A, gene NC_003977.2, start position 1712, end position 1819, nucleotide sequence AAAGACTGTGTGTTTAATGAGTGGGAGGAGTT, cell barcode ATGTCCCCATAACAGA, unique molecular identifier GTATGATATTCG;

[0058] Example A, gene NC 003977.2, start position 1511, end position 1633, nucleotide sequence C0GACCGACCACGGGGCGCACCTCTCTTTACG, cell barcode CTCGAGGTCTGAGGCC, unique molecular identifier TTGAAACATTTA;

[0059] Example A, gene NC 003977.2, start position 1511, end position 1606, nucleotide sequence C0GACCGACCACGGGGCGCACCTCTCTTTACG, cell barcode GATTCGATCAAACGAA, unique molecular identifier GAAGCAACCCCT;

[0060] Example A, gene NC 003977.2, start position 998, end position 1147, nucleotide sequence AATTGTGGGTCTTTTGGGTTTTGCCGCCCCTT, cell barcode ACCTACCGTCCOGGTA, unique molecular identifier TATTTTCCAAGG;

[0061] Example A, gene NC 003977.2, start position 998, end position 1147, nucleotide sequence AATTGTGGGTCTTTTGGGTTTTGCCGCCCCTT, cell barcode ACCTACCGTCCOGGTA, unique molecular identifier TATTTTCCAAGG;

[0062] Example B, gene NC 003977.2, start position 1594, end position 1743, nucleotide sequence CACCTCIGCACGTOGCATGGAGACCACCGTGA, cell barcode CIGGACCGTAACACAT, unique molecular identifier CTCCATTAAGAT;

[0063] Example B, gene NC 003977.2, start position 1668, end position 1817, nucleotide sequence TIGGACTTTCAGCAATGTCAACGACCGACCTT, cell barcode CCCTCTCICGTTCAGA, unique molecular identifier TAGTCTAAAAGG;

[0064] Example B, gene NC 003977.2, start position 3093, end position 3182, nucleotide sequence CCTCCTGCCTCCACCAATOGGCAGTCAGGAOG, cell barcode GACGCTGGCAGAGACA, unique molecular identifier GTATTCIGAGCA;

[0065] Example B, gene NC 003977.2, start position 1657, end position 1800, nucleotide sequence TAAGAGGACTCTTGGACTTTCAGCAATGICAA, cell barcode TICAGGATCGAACTCA, unique molecular identifier ACATGCGGGCAT;

[0066] Example C, gene NC 003977.2, start position 1694, end position 1819, nucleotide sequence GACCTIGAGGCATACTTCAAAGACIGTGIGTT, cell barcode AGGGTTTICCAAGCAT, unique molecular identifier CACCGACTAGTG;

[0067] Example C, gene NC 003977.2, start position 1673, end position 1819, nucleotide sequence CTTTCAGCAATGTCAACGACCGACCTIGAGGC, cell barcode AAGCGTIICTGICGTC, unique molecular identifier TTAGGTATCAAT;

[0068] Example C, gene NC 003977.2, start position 482, end position 631, nucleotide sequence TAATTCCAGGATCATCAACAACCAGCACCGGA, cell barcode TTIGGTTTCACTCOG, unique molecular identifier ATATCATGACCT;

[0069] Example C, gene NC 003977.2, start position 482, end position 631, nucleotide sequence TAATTCCAGGATCATCAACAACCAGCACOGGA, cell barcode TTTGGTTICACTCCGT, unique molecular identifier ATATCATGACCT.

[0070] In the fifth step, the expression level and fragment distribution are calculated: the number of sequencing sequences covering each base of the genome is calculated, the coverage depth of each base of the HBV genome is plotted, and the gene position is marked and displayed, to obtain the expression level and fragment distribution of the single-cell RNA-seq sequencing data on the HBV genome as shown in FIG. 6. Figure 4 In the fifth step, the expression level and fragment distribution are calculated: the number of sequencing sequences covering each base of the genome is calculated, the coverage depth of each base of the HBV genome is plotted, and the gene position is marked and displayed, to obtain the expression level and fragment distribution of the single-cell RNA-seq sequencing data on the HBV genome as shown in FIG. 6.

[0071] In the identification of human and HBV sequences in the sequencing data, the specific steps are as follows: first, the sequencing raw data FASTQ file obtained by scRNAseq is aligned to the HBV transcriptome reference sequence using bwa mem software (version v0.7.17). In the alignment process, the sequencing sequence can be matched with the HBV transcriptome reference sequence through the algorithm of the software. Only the sequencing sequence information that can be aligned to the HBV genome is retained, and a BAM format alignment record file aligned to the HBV genome is generated. In this BAM file, detailed alignment information of each sequencing sequence is included.

[0072] For determining the breakpoint position of HBV integration into human genome, the CIGAR record information of each sequencing sequence in the BAM file is mainly relied on. The CIGAR record describes the alignment of the sequencing sequence and the reference sequence in detail, and the Softclip (S) structure is the key information. The sequencing sequence with the Softclip structure is extracted, and the Softclip part is usually the part that is "cut off" due to the incomplete matching of the sequence with the reference sequence in the alignment process. Since the integration of HBV into the human genome can cause such abnormal alignment of the sequencing sequence, the possible breakpoint position can be inferred by analyzing the position and sequence characteristics of the Softclip part. Specifically, the Softclip part is cut off, and the remaining sequence and the corresponding transcriptome coordinate position are retained. For the sequencing sequence without the Softclip structure, the sequence and the corresponding transcriptome coordinate position are retained. Through such analysis of a large number of sequencing sequences, the alignment of the sequence, the distribution of the Softclip structure, and the transcriptome coordinate position are comprehensively considered, and the breakpoint position of HBV integration into the human genome can be accurately determined. For example, if multiple sequencing sequences are found to have Softclip structures at similar positions in a certain region, and the region is associated with the known sequence characteristics of the human genome and HBV genome, it can be highly suspected that the region is the breakpoint position of HBV integration, and the breakpoint position can be determined after further confirmation by existing experimental verification methods.

[0073] The present embodiment solves the problem that there is currently a lack of methods for calculating the expression level of HBV from single-cell sequencing data and positioning the expression abundance of each region on the HBV genome, and fills the gap in HBV analysis methods based on single-cell resolution. Compared with existing HBV analysis methods, the present embodiment has obvious advantages:

[0074] Firstly, existing alignment tools are for sequence alignment of linear genomes, and the genome sequence of HBV (NCBI nucleic acid ID number: NC_003977.2) is also recorded in the database in a linear manner. However, the genome of HBV is a circular genome, and the S and P genes of HBV span the start and end positions of the HBV genome. The existing tools will cause deviations in the quantification of the expression levels of the S and P genes when aligning the sequencing data to the HBV genome. The present embodiment achieves more accurate quantification of gene expression levels by aligning the sequencing data to the linear transcriptome and indirectly aligning the sequencing data to the genome according to the mapping relationship between the transcriptome and the genome;

[0075] Secondly, there are many overlapping genes on the HBV genome, and there are multiple transcripts of the same gene. The method of aligning data to the transcriptome and then mapping to the genome position in the present embodiment can better distinguish the expression levels of overlapping genes and different transcripts of a gene;

[0076] Finally, by sequencing the chimeric human-HBV genome, the breakpoint of HBV integration into human genome can be further determined, providing more abundant research information in addition to expression levels.

[0077] Through this embodiment, the intercellular heterogeneity can be revealed: HBV infection can lead to different response patterns in different cell types or within the same type of cells. Through scRNA-seq, it can be identified which cell types are more susceptible to HBV infection and how these cells regulate their gene expression during infection.

[0078] This embodiment overcomes the unique biological characteristics of HBV (such as lack of polyA tail, circular genome structure), directly applies existing scRNA-seq technology and analysis process to study the contradictory problems of HBV, and accurately quantifies and locates the expression activities on the HBV genome by constructing accurate mapping relationship between transcriptome and genome.

[0079] In the process of developing this technical solution, there are many challenges, among which the technical problems brought by the unique biological characteristics of HBV are particularly prominent. The key to breaking the routine lies in the handling of the circular structure of HBV genome and the gene spanning phenomenon.

[0080] Traditional alignment tools are usually used for sequence alignment of linear genomes, while the genome of HBV is a circular genome, and its S and P genes span the start and end positions of the HBV genome. This special structure causes the expression level of S and P genes to be quantitatively biased when using existing tools to align sequencing data to the HBV genome. This is a technical bias that has existed in the field for a long time, that is, it is generally believed that the existing alignment methods for linear genomes can be directly applied to the analysis of HBV genome.

[0081] To overcome this bias, this solution breaks the conventional thinking and no longer directly aligns the sequencing data to the circular genome, but aligns the sequencing data to the linear transcriptome, and indirectly aligns the sequencing data to the genome according to the mapping relationship between the constructed transcriptome and genome. Specifically, in the first step, the HBV genome sequence, transcriptome sequence and genome annotation file are downloaded from NCBI, the position information of each HBV gene in the transcriptome and genome is carefully constructed, and the mapping relationship between the transcriptome and genome is obtained. In the subsequent steps, the raw sequencing data obtained by scRNA-seq is first aligned to the HBV transcriptome reference sequence, and human and HBV sequences are identified and separated, and then the HBV source transcriptome sequence is converted to the alignment position on the genome through the mapping relationship. This innovative method effectively solves the alignment bias problem caused by the special structure of the HBV genome, and realizes more accurate quantification of gene expression levels.

[0082] Before obtaining the final solution, the R&D team tried to directly use existing alignment tools for linear genomes, such as Bowtie2, to align the raw sequencing data from scRNAseq directly to the circular genome of HBV. However, due to the special structure of the HBV genome, especially the S and P genes spanning the start and end positions, the alignment results of these two genes deviate seriously from the actual situation, and the expression level quantification deviates greatly. After in-depth analysis of the results, it was found that this method could not accurately distinguish different transcripts of genes and the expression levels of overlapping genes, and could not give reliable analysis results for the complex gene expression on the HBV genome.

[0083] Another abandoned solution is to try to make simple modifications to existing alignment tools to adapt to the circular genome of HBV. For example, try to adjust the parameters in the alignment algorithm to consider the gene spanning situation to some extent. But after a lot of experimental verification, this method has certain improvement on part of the data, but the overall effect is still not ideal, and cannot meet the requirements of comprehensive and accurate analysis of HBV genome.

[0084] In terms of accuracy of gene expression level quantification, this embodiment takes S and P genes as an example. Compared with the traditional method of directly aligning to the circular genome, the innovative method of aligning the sequencing data to the transcriptome and then mapping to the genome can more accurately quantify the expression levels of these two genes spanning the start and end positions of the genome. Assuming that in the traditional method, the quantification deviation of S gene expression level may reach 30%-50% (this is the deviation range estimated due to alignment errors caused by circular structure), and after using this solution, through the accurate construction and application of the mapping relationship between transcriptome and genome, this deviation can be controlled within 5%-10%, greatly improving the accuracy of quantification.

[0085] For distinguishing the expression levels of overlapping genes and different transcripts of genes, this solution can clearly distinguish the expression signals of different transcripts and overlapping genes through a unique alignment and analysis process. For example, when analyzing multiple overlapping gene regions on the HBV genome, the traditional method may not be able to accurately distinguish the expression contributions of different genes, resulting in chaotic analysis results. This solution can accurately identify the expression of each gene transcript, increase the expression analysis resolution of overlapping gene regions by at least 2-3 times, and can more detailedly study the expression regulation mechanism of HBV genes.

[0086] In the determination of the breakpoint position of HBV integration into the human genome, the present scheme can provide more abundant research information through the analysis of the sequencing sequence of human-HBV genome chimeric. Compared with the traditional method, more potential breakpoint positions can be detected. Assuming that the traditional method can only detect 5-8 breakpoint positions, the present scheme can increase the number of detected breakpoints to 15-20, greatly increasing the depth and breadth of the study of HBV integration mechanism.

[0087] Example Two

[0088] In this embodiment, in the process of determining the breakpoint position of HBV integration into the human genome, the sequence alignment, soft cut structure distribution and transcriptome coordinate position are considered comprehensively as follows:

[0089] For sequence alignment, the alignment results of the sequencing sequence and the HBV transcriptome reference sequence are analyzed in detail. If the base matching rate of the sequencing sequence and the reference sequence in a certain region is lower than the preset threshold, and the low matching rate shows continuous distribution, which is significantly different from the matching mode of the normal alignment region, the region is marked as a potential breakpoint position region. The preset threshold is 70%-80%, which ensures that it can effectively distinguish between normal and abnormal alignment. For example, when the base matching rate of the sequencing sequence and the reference sequence in a certain region is lower than 75%, and the low matching rate shows continuous distribution, i.e. the number of consecutive bases in the low matching rate region exceeds 10 (the length is determined according to the average fragment length of the HBV transcriptome reference sequence and the allowable range of sequencing error), which is significantly different from the matching mode of the normal alignment region, the region is marked as a potential breakpoint position region. The purpose of such setting is to avoid misjudging the low matching caused by accidental sequencing error as a potential breakpoint position. Because in normal sequencing process, due to various factors, a small amount of base mismatch may occur, but this mismatch is usually scattered and does not show continuous low matching. When the matching rate of consecutive multiple bases is lower than 75%, it is likely that the sequence structure has been changed due to factors such as HBV integration, thereby affecting the alignment results.

[0090] For the distribution of soft cutting structure, statistical analysis is performed on the sequencing sequences with soft cutting structure, the starting position and the ending position of each soft cutting structure in the sequence are counted, and these positions are mapped to the transcriptome coordinates; when the occurrence frequency of soft cutting structure in a specific transcriptome coordinate interval exceeds the preset frequency threshold, the coordinate interval is marked as a potential breakpoint position interval. The preset frequency threshold is that the number of soft cutting structure appearing in each 1000 sequencing sequences exceeds 10 times, which ensures that it can accurately screen out the soft cutting structure aggregation area caused by HBV integration. When the occurrence frequency of soft cutting structure in a specific transcriptome coordinate interval exceeds this threshold, for example, the number of soft cutting structure appearing in a certain 100bp transcriptome coordinate interval reaches 15 times or more, and the average length, maximum length and other parameters of soft cutting structure appear significant differences compared with other normal regions (such as the average length is more than 20% longer than the normal region), the coordinate interval is marked as a potential breakpoint position interval. This is because when HBV integrates into the human genome, it may cause abnormal soft cutting of sequencing sequences at certain positions. By counting the frequency and related length parameters of soft cutting structure, these regions that may be related to HBV integration can be effectively screened out.

[0091] In combination with the transcriptome coordinate position, the transcriptome coordinate continuity of the sequencing sequence in the human-HBV genome chimeric region is analyzed. If the transcriptome coordinate is found to jump, interrupt, or not match the known human genome and HBV genome transcriptome coordinate range, the corresponding position is marked as a potential breakpoint position; coordinate changes at positions prone to integration, such as gene boundary regions and repeat sequence regions, are particularly concerned. For example, when the transcriptome coordinate suddenly jumps from the human genome coordinate range to the HBV genome coordinate range at a certain point, and the sequences before and after the jump are highly similar to the specific conserved sequences of the human genome and the HBV genome, respectively, the jump position is most likely the breakpoint position of HBV integration.

[0092] Finally, the potential breakpoint position regions marked by the above three ways are comprehensively considered, and when multiple marked regions overlap or are densely distributed in a small range, the common region or dense region is determined as the breakpoint position of HBV integration into the human genome. Among them, "densely distributed in a small range" specifically refers to the presence of 3 or more potential breakpoint position regions marked by different ways within a range of 50bp radius centered on a certain position. When this situation occurs, the common region or dense region is determined as the breakpoint position of HBV integration into the human genome. This is because multiple independent judgment factors point to the possible existence of breakpoints in the nearby region, which greatly increases the credibility of the region as the true breakpoint position. Through this comprehensive analysis method, the accuracy and reliability of the determination of the breakpoint position can be effectively improved, and the misjudgment caused by a single factor can be avoided.

[0093] The single-cell RNA-seq-based HBV expression level and fragment distribution analysis method of the present application significantly improves the recognition accuracy, especially in the accurate recognition of breakpoint positions. In traditional analysis, since the gene bank stores linear sequences, the HBV genome is a circular structure, and part of the gene spans the beginning and end. When using traditional tools to align the sequencing data with the circular HBV genome, inaccurate alignment problems may occur in special areas such as spanning the beginning and end and breakpoints due to structural differences. This is because conventional alignment algorithms are designed based on the characteristics of linear genomes and are difficult to handle the special sequence characteristics and position relationships of circular genomes.

[0094] The present application innovatively transcribes the circular HBV genome into linear transcriptome RNA, eliminating the alignment complexity brought by the circular structure, and its sequence is more suitable for the linear reference sequence in the gene bank. By aligning the transcriptome RNA with the reference sequence in the gene bank, not only is the alignment deviation caused by the circular structure avoided, but also the foundation for subsequent accurate recognition of breakpoint positions is laid. From the principles of molecular biology and bioinformatics, the transcription process faithfully records the coding information of the HBV genome, and presents it in a linear form, so that the alignment algorithm can more accurately identify the matching region and reduce the mismatch or omission. At the same time, by constructing the mapping relationship between the transcriptome and the genome, the transcriptome alignment result is accurately mapped back to the genome, further ensuring the accuracy of the analysis of the expression level and fragment distribution of each region on the HBV genome, and realizing the accurate determination of the breakpoint position. Therefore, the method of the present application effectively overcomes the limitations of traditional methods and significantly improves the accurate recognition ability of HBV-related sequences and breakpoint positions.

[0095] Example Three

[0096] In the process of determining the breakpoint position of the integrated HBV into the human genome, the sequence alignment, soft cut structure distribution and transcriptome coordinate position are considered comprehensively as follows:

[0097] For sequence alignment, the alignment results of the sequencing sequence and the HBV transcriptome reference sequence are analyzed in detail. If the base matching rate of the sequencing sequence and the reference sequence in a certain region is less than 75%, and the low matching rate shows a continuous distribution, i.e. the number of continuous bases in the low matching rate region exceeds 15 bp, which is significantly different from the matching mode of the normal alignment region, then the region is marked as a potential breakpoint position region. The length of 15 bp is determined according to the average fragment length of the HBV transcriptome reference sequence and the allowable range of sequencing error, in order to avoid misjudging the occasional sequencing error as a potential breakpoint position.

[0098] For the distribution of soft cutting structure, statistical analysis is performed on the sequencing sequences with soft cutting structure, the starting position and the ending position of each soft cutting structure in the sequence are counted, and these positions are mapped to the transcriptome coordinates. When the occurrence frequency of soft cutting structure in a certain transcriptome coordinate interval exceeds the preset frequency threshold of 10 occurrences per 1000 sequencing sequences, and the average length of soft cutting structure in the interval is at least 10 bp longer than that in the normal region (for example, the average length of the normal region is 15 bp, the average length of the interval is 25 bp or more), and the maximum length is at least 20 bp longer than that in the normal region (for example, the maximum length of the normal region is 20 bp, the maximum length of the interval is 40 bp or more), the coordinate interval is marked as a potential breakpoint position interval.

[0099] In combination with the transcriptome coordinate position, the transcriptome coordinate continuity of the sequencing sequence in the chimeric region of human-HBV genome is analyzed. If the transcriptome coordinate is found to jump, interrupt, or not match the known human genome and HBV genome transcriptome coordinate range, the corresponding position is marked as a potential breakpoint position region. Especially for the found transcriptome coordinate jump, further analysis of the sequence characteristics before and after the jump is performed, if the sequence before the jump has a similarity of 95% or more with a certain conserved sequence of the human genome, and the sequence after the jump has a similarity of 90% or more with a certain conserved sequence of the HBV genome, the jump position is taken as a potential breakpoint position for attention.

[0100] Finally, the potential breakpoint position regions and potential breakpoint position intervals marked by the above three ways are comprehensively considered, when multiple marked regions overlap with each other or within a range of 50 bp radius centered at a certain position, 3 or more potential breakpoint position regions marked by different ways appear, the common region or dense region is determined as the breakpoint position when HBV integrates into the human genome.

[0101] In the analysis of sequence alignment to determine the potential breakpoint position region, the number of consecutive bases for the continuity distribution is 15 bp, which is determined by the average fragment length of 200 bp of the HBV transcriptome reference sequence, combined with the allowed range of sequencing errors, through experiments and data analysis, to avoid misjudgment of low matching caused by accidental sequencing errors as potential breakpoint positions, and to ensure that normal alignment and abnormal alignment caused by HBV integration can be effectively distinguished.

[0102] In the statistics of soft cutting structure distribution to determine the potential breakpoint position interval, the preset frequency threshold is set to 10 occurrences of soft cutting structure per 1000 sequencing sequences. At the same time, the difference requirements for the average length and the maximum length of soft cutting structure and the normal region, such as the average length being at least 10 bp longer than the normal region and the maximum length being at least 20 bp longer than the normal region, are used to improve the accuracy of screening potential breakpoint position intervals.

[0103] Wherein, when determining the potential breakpoint position in combination with the transcriptome coordinate position, the requirement for the similarity of the sequence before and after the jump to the corresponding genomic conserved sequence, i.e. the similarity of the sequence before the jump to a certain conserved sequence of the human genome is 95% or above, and the similarity of the sequence after the jump to a certain conserved sequence of the HBV genome is 90% or above, is determined through comparison and analysis of a large amount of known human genome and HBV genome conserved sequence data, combined with actual verification experiments, which ensures that the transcriptome coordinate jump position caused by HBV integration can be accurately identified as a potential breakpoint position.

[0104] The present embodiment determines the HBV integration breakpoint position through multi-dimensional accurate analysis, which significantly improves the accuracy. Not only can sequencing errors be effectively filtered and potential regions be accurately identified, but also the soft cutting structure and transcriptome coordinate characteristics can be comprehensively considered to avoid the omission of single-dimensional analysis. Multi-evidence comprehensive judgment enhances the reliability and comprehensively guarantees the integrity of the breakpoint determination. Based on the accurate breakpoint position determination, the analysis method of HBV expression level and fragment distribution based on single-cell RNA-seq can complete all the calculation and analysis through one software program, which makes the operation more accurate and fast.

[0105] The above is only an embodiment of the present application, and well-known specific technical solutions and / or common knowledge of characteristics in the scheme are not described in detail. It should be noted that for those skilled in the art, without departing from the technical solutions of the present application, a number of modifications and improvements can be made, which should also be considered as the protection scope of the present application, and these will not affect the effect and practicality of the present application. The protection scope of the present application should be subject to the content of its claims, and the specific implementation mode and the like in the specification can be used to explain the content of the claims.

Claims

1. An analysis method based on single-cell RNA-seq of HBV expression level and fragment distribution, characterized in that, The method comprises the following steps: Step 1: Construct the position information of each gene of HBV in the linear transcriptome and the circular genome to obtain the mapping relationship between the transcriptome and the genome; Step 2: Align the raw sequencing data obtained by scRNAseq to the HBV transcriptome reference sequence to determine the breakpoint position when HBV integrates into the human genome, and identify the human sequence and the HBV sequence; Step 3: Remove the human sequence, retain the remaining HBV sequence and the corresponding transcriptome coordinate position from the potential breakpoint position of HBV to obtain the HBV source transcriptome sequence; Step 4: Convert the position of the sequencing sequence aligned on the transcriptome to the aligned position on the genome by the mapping relationship between the transcriptome and the genome constructed in step 1; Step 5: Calculate the number of sequencing sequences covering each base of the genome, draw the coverage depth of each base of the HBV genome, and mark the gene position to obtain the expression level and fragment distribution of single-cell RNA-seq sequencing data on the HBV genome; In the process of determining the breakpoint position when HBV integrates into the human genome, the sequence alignment, soft cut structure distribution and transcriptome coordinate position are comprehensively considered by the following methods: For sequence alignment, the alignment results of the sequencing sequence and the HBV transcriptome reference sequence are analyzed in detail. If the base matching rate of the sequencing sequence and the reference sequence in a certain region is lower than a preset threshold, and the low matching rate shows a continuous distribution and is significantly different from the matching mode of the normal alignment region, the region is marked as a potential breakpoint position region; the preset threshold is 70%-80%; For the distribution of soft cut structure, the starting position and the ending position of each soft cut structure in the sequence are statistically analyzed, and these positions are mapped to the transcriptome coordinates. When the appearance frequency of the soft cut structure in a specific transcriptome coordinate interval exceeds a preset frequency threshold, the coordinate interval is marked as a potential breakpoint position interval; the preset frequency threshold is that the number of soft cut structures appearing in every 1000 sequencing sequences exceeds 10 times; In combination with the transcriptome coordinate position, the transcriptome coordinate continuity of the sequencing sequence is analyzed in the human-HBV genome chimeric region. If the transcriptome coordinates are found to jump, interrupt or be inconsistent with the known transcriptome coordinate range of the human genome and the HBV genome, the corresponding position is marked as a potential breakpoint position region; Finally, the potential breakpoint position regions and potential breakpoint position intervals marked by the above three methods are comprehensively considered. When multiple marked regions overlap or densely distribute in a small range, the common region or the dense region is determined as the breakpoint position when HBV integrates into the human genome.

2. The method of claim 1, wherein the method is based on single-cell RNA-seq analysis of HBV expression levels and fragment distribution. In step one, the HBV genome reference sequence and the HBV transcriptome reference sequence are obtained respectively, the annotation information of the HBV genome is obtained, the position information of the HBV genes in the transcriptome and the genome is constructed, and the mapping relationship between the transcriptome and the genome is obtained; the mapping relationship between the transcriptome and the genome is that the starting position 1376 and the ending position 1840 of NC_003977.2 are the X gene, which is the X protein; the starting position 1816 and the ending position 2454 of NC_003977.2 are the C gene, which is the pre-capsid protein; the starting position 1903 and the ending position 2454 of NC_003977.2 are the C gene, which is the capsid protein; the starting position 2309 and the ending position 3182 of NC_003977.2 are the P gene, which is the polymerase; the starting position 1 and the ending position 1625 of NC_003977.2 are the P gene, which is the polymerase; the starting position 2850 and the ending position 3182 of NC_003977.2 are the S gene, which is the large envelope protein; the starting position 1 and the ending position 837 of NC_003977.2 are the S gene, which is the large envelope protein; the starting position 3174 and the ending position 3182 of NC_003977.2 are the S gene, which is the medium membrane protein; the starting position 1 and the ending position 837 of NC_003977.2 are the S gene, which is the medium membrane protein; the starting position 157 and the ending position 837 of NC_003977.2 are the S gene, which is the small envelope protein.

3. The method of claim 1, wherein the method is based on single-cell RNA-seq analysis of HBV expression levels and fragment distribution. In step two, the sequencing raw data FASTQ file obtained by scRNAseq is aligned to the HBV transcriptome reference sequence by using the bwa mem software, only the sequencing sequence information that can be aligned to the HBV genome is retained, and the BAM format alignment record file aligned to the HBV genome is generated as the HBV source sequence and the rest is the human sequence.

4. The method of claim 3, wherein the method is based on single-cell RNA-seq analysis of HBV expression levels and fragment distribution. In step three, according to the CIGAR record information of each sequencing sequence in the BAM file, the sequencing sequence with soft cut structure is extracted, and the soft cut part is cut off, and the remaining sequence and the corresponding transcriptome coordinate position are retained; the sequencing sequence without soft cut structure is retained, and the sequence and the corresponding transcriptome coordinate position are retained.

5. The method of claim 1, wherein the method is based on single-cell RNA-seq analysis of HBV expression levels and fragment distribution. The dense distribution in a smaller range refers to that, in a range with a radius of not more than 50 bp centered on a certain position, there are 3 or more potential breakpoint position regions marked by different ways.

6. The method of claim 1, wherein the method is based on single-cell RNA-seq analysis of HBV expression levels and fragment distribution. When analyzing the sequence alignment to determine the potential breakpoint position region, the continuous distribution of the low matching rate refers to that the number of continuous bases in the low matching rate region is more than 10.

7. The method of claim 1, wherein the method is based on single-cell RNA-seq analysis of HBV expression levels and fragment distribution. When determining the potential breakpoint position interval by analyzing the distribution of soft cut structure, the length distribution of the soft cut structure is also considered; when the frequency of the soft cut structure in a certain transcriptome coordinate interval is more than the preset frequency threshold, and the average length of the soft cut structure is at least 10 bp longer than that in the normal region and the maximum length is at least 20 bp longer than that in the normal region.

Images