Single cell 3'end transcriptome variable poly-adenosine acidification analysis method

Abstract

The invention discloses a single cell 3'end transcriptome variable polyadenylation (APA) analysis method, and belongs to the field of bioinformatics. The analysis process provided by the invention is based on 3'end transcriptome sequencing data of a single cell, combines multiple processes of UMI deduplication, false comparison removal, site merging, internal amplification removal, polyA site removal with few supporting cells, site verification and the like, and can realize efficient detection and quantification of APA at the single cell level. The method is especially suitable for single cell APA analysis of single cell 3'end transcriptome sequencing data, including 3 'end single cell sequencing data after 10x scRNA-Seq cell subset combination, 3' end single cell sequencing data based on a C1 chip, and even full-length transcriptome data.

Description

Technical Field

[0001] This invention relates to the field of bioinformatics, specifically to an analysis workflow for single-cell 3' end sequencing results based on next-generation sequencing (NGS) technology, particularly a method for the detection and analysis of alternative polyadenylation (APA) sites. This method can be widely applied in genomics research, especially for the precise detection of APA sites at the single-cell level. Background Technology

[0002] Alternative polyadenylation (APA) refers to the phenomenon during gene transcription where, due to the presence of multiple polyadenylation signal sites, RNA polymerase II can selectively cleave at different sites and add a poly(A) tail, thereby producing mRNA isoforms with different 3'UTR sequences. This selectivity is one of the important mechanisms for gene expression regulation.

[0003] Single-cell omics refers to the analysis of individual cells through genomics, transcriptomics, proteomics, and epigenomics to reveal cellular heterogeneity, function, and regulatory mechanisms. The rise of single-cell omics has provided unprecedented opportunities for studying cell development, cell fate determination, disease mechanisms, and the discovery of therapeutic targets. Current techniques mainly include single-cell isolation technology, single-cell nucleic acid sequencing technology (RNA sequencing (scRNA-seq), whole-genome sequencing, and single-cell exonuclease sequencing), single-cell proteomics analysis technology, and single-cell spatial transcriptomics technology. While various omics analyses at the single-cell level are gradually being perfected, selective polyadenylation (APA) is an important and widespread post-transcriptional regulatory process at the true single-cell level. It affects important processes such as mRNA stability, translation efficiency, and protein sublocalization by selecting different PolyA sites. Many methods exist for analyzing selective polyadenylation, but each has its own limitations. For example, single-cell transcriptome sequencing, exemplified by 10x's scRNA-seq, involves a large cell volume but low sequencing yield per cell, making it unsuitable for calculating the utilization rate of distant polyA sites. Typically, sequencing data from a cell subpopulation are combined before calculating the distant utilization rate, which is a multi-cell level rather than a single-cell level approach. Furthermore, the lower edge of the sequencing read peak map varies in distance from the actual polyA sites, making the polyA sites obtained using this method inaccurate and often significantly different from their true locations, serving only as an approximation. Similarly, methods based on deep learning or even large-scale learning models (LLMs) suffer from poor interpretability. Methods that count the number of sequencing reads falling into each polyA site based on existing polyA database sites are unsuitable for discovering new polyA sites, and some polyA sites in the existing database may be erroneous. Analysis methods based on control samples require a global reliance on a fixed sample, with other samples used relative to it to calculate whether each transcript increases or decreases. Such methods are unsuitable for sparse single-cell data with many undetectable reads for various genes.

[0004] Therefore, it is of great significance to provide an accurate and reliable method for studying polyA sites at the single-cell level, and to improve the sequencing and analysis methods for identifying, screening, and filtering polyA supporting reads. Summary of the Invention

[0005] This invention proposes a method for analyzing variable polyadenylation of the 3' end transcriptome in single cells.

[0006] To achieve the above objectives, the present invention is implemented through the following technical solution: A method for analyzing variable polyadenylation at the 3' end of a single-cell transcriptome includes the following steps: (1) Obtaining high-throughput sequencing data: Single cells are separated by microarray or flow cytometry. The nucleic acids of single cells are amplified by polyT primers, libraries are constructed and sequenced. The sequencing results are split to obtain the sequencing data of the 3' transcriptome of each cell. (2) Screening polyA support reads and aligning and obtaining polyA site coordinates: Screening polyA support reads from the sequencing data in step (1); Aligning the selected polyA support reads to the reference genome using alignment software; After alignment, retaining high-quality alignments, removing false alignments, and then using unique molecular identifiers to remove duplicates caused by polymerase chain reaction in the sequencing data of each cell, obtaining the Bam file of the alignment coordinates of the sequencing reads of each cell, converting the Bam file to a bed file, obtaining the coordinates of the 3' end of the reads through the bed file, i.e., the polyA site, performing polyA site merging, filtering and verification, and using the verified reliable polyA sites as the final polyA sites; False alignments may be introduced during library preparation using polyT sequences. These false alignments can result in false polyA sites, affecting subsequent analysis. Therefore, it's necessary to actively remove these false polyA sites generated by the false alignment sequences. Specifically, the method for removing false alignments is as follows: In the alignment file, filter out alignments dominated by polyA or polyT. That is, in the sequencing reads aligned to the reference genome, the sequence corresponding to M in the CIGAR string is the one that plays a matching role, and these alignments contain more than 50% consecutive polyA or polyT sequences. (3) Quantification of single-cell APA data: After merging all cells obtained in step (2) and verifying reliable polyA site locations and merge records, count the number of polyA support reads at polyA sites in each single-cell bed file, and generate a sequencing read matrix with polyA sites as rows, cells as columns, and the support count value as the content value. (4) Quantification of gene distal usage: The proportion of distal polyA usage for each gene is calculated to obtain the analysis results. Specifically, it is an index weighted by the percentage of the final polyA read in the gene according to its position, reflecting the difference in polyA site usage status of each gene at the single-cell level. The analysis results are used to characterize the differences in variable polyadenylation between different cell types or cell states; (5) Use the snakemake script to automate the execution of the process and improve data processing efficiency.

[0007] Furthermore, the polyT primers in step (1) are primers with cell tags and unique molecular identifiers.

[0008] Furthermore, the process of selecting polyA-supporting reads in step (2) involves selecting sequencing reads containing 10 consecutive A's or ending with 6 consecutive A's, and removing these A's and downstream bases.

[0009] Furthermore, the alignment software mentioned in step (2) is selected from alignment software that supports RNA alignment, such as STAR, HISAT2, TopHat2 / TopHat.

[0010] Furthermore, the high-quality alignment described in step (2) is the only alignment with a high MAPQ value.

[0011] Furthermore, the polyA site filtering in step (2) includes the following steps: s1: Merging, filtering, and annotation of polyA sites: All BED files containing 3' end positions of cells are merged into one BED file, sorted by chromosome classification, and the 3' end positions of each read are obtained, i.e., the original polyA sites; original polyA sites clustered less than 12 nt are merged and the merging process is recorded. Among the polyA sites merged into a cluster, the highest genomic coordinate of the read is used as the coordinate of the polyA site, and the sum of all reads in the cluster is used as the read count of the polyA site; internal amplification is removed, and polyA sites with few supporting cells, such as polyA sites supported by less than 1%-5% of cells, are filtered out to obtain the final polyA sites; the final polyA sites are annotated according to the reference genome. The annotation content includes: gene name, gene region or upstream and downstream base sequences; s2, polyA site verification: The polyA sites obtained in step s1 are verified using at least three different verification methods to determine the reliability of the final polyA sites.

[0012] Furthermore, the verification method described in step s2 is selected from: Plot the base frequencies 100 nt upstream and downstream of the polyA site; Alternatively, draw a polyA signal enrichment map of nucleotides 60 nt upstream of polyA; Or plot the polyA number distribution at each genomic location; Alternatively, compare the number of overlaps with sites in existing publicly available polyA databases. Alternatively, compare the distribution of relative distances with sites in the existing public polyA database.

[0013] The criteria for judging the reliability of the final polyA site in the verification method are as follows: (A) The standard for judging based on the base frequency map of 100nt upstream and downstream of the polyA site is: the base frequency of 100nt upstream and downstream of the polyA site is expected to have an A base peak near the upstream (-20nt) and a sudden change in C base frequency at the polyA site (0nt). (B) The criteria for judging the polyA signal enrichment map of nucleotides 60 nt upstream of polyA is: the motif map of 60 nt upstream of polyA site, the expected main motif should be A[A / T]TAAA; (C) The criteria for judging the polyA number distribution map of each genomic location are: pie chart of the number of polyA sites at different genomic locations. It is expected that the number of polyA sites in the 3'UTR region of the gene should be relatively large, for example, the proportion of sites falling in PA, 3'UTR, and 3'UTR extension region exceeds 30%; (D) The criterion for comparing the number of overlaps with existing polyA sites in the public polyA database is: the overlap between the detected polyA sites and polyA sites in the public polyA-DB2 database. A distance of less than 12 nt is considered an overlap, and it is expected that there should be at least 30% overlap. (E) The criterion for comparing the distribution of relative distances between the detected polyA sites and the polyA sites in the publicly available polyA database is: the main peak of the expected distance distribution map of the detected polyA sites and the polyA sites in the publicly available polyA-DB2 database should be mainly near 0.

[0014] Furthermore, the removal of internal amplification in step s1 is as follows: removing amplifications in the A region of the genome that may not be caused by polyA tails. Specifically, the filtering conditions can be: filtering out polyA sites downstream of the polyA site where there are more than 15 A's in 6 consecutive A's and / or 20 bases.

[0015] Furthermore, the cell state described in step (4) includes cell activation, cell differentiation, or disease-related state.

[0016] The single-cell 3' end variable transcriptome polyadenylation (APA) analysis method of this invention uses an analytical device comprising a processor and a memory. The processor includes an acquisition module, a site identification module, a usage calculation module, and a result output module. The coordinated use of these modules achieves the analysis of selective polyadenylation of the single-cell 3' end transcriptome. The processor and memory described in this analytical device are existing conventional technologies and devices, and the acquisition module, site identification module, usage calculation module, and result output module within the processor are also existing publicly disclosed technologies in the field of computer science.

[0017] The present invention provides a method for analyzing variable polyadenylation at the 3' end of the single-cell transcriptome, which has the following advantages: (1) High detection accuracy: Through an optimized analysis process, this invention uses polyA support reads to accurately capture polyA sites at the single-cell level, calculate their usage ratio, and clarify their usage status. (2) The analysis process proposed in this invention is based on the 3' end transcriptome sequencing data of single cells, and combines multiple processes such as UMI deduplication, removal of false alignments, site merging, removal of internal amplification, removal of polyA sites with few supporting cells, and site verification, which can achieve efficient detection and quantification of APA at the single-cell level.

[0018] (3) High scalability: This method is particularly suitable for single-cell APA analysis of single-cell 3' transcriptome sequencing data, including 3' single-cell sequencing data after merging 10x scRNA-Seq cell subsets, 3' single-cell sequencing data based on C1 chip, and even full-length transcriptome data, such as Smart-Seq2 sequencing data.

[0019] (4) Automation and visualization support: Automated data processing is achieved through the snakemake script, and the relative height of each polyA site peak is easily viewed through the visualization tool based on igv.js. Attached Figure Description

[0020] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0021] Figure 1 This is a schematic flowchart of the sequencing and analysis method of the present invention; Figure 2This is a single-cell level polyA site verification diagram of T cells during antibody activation in Example 1; where A is a base frequency diagram of 100 nt upstream and downstream of the polyA site, it is expected that there should be an A base peak near the upstream (-20 nt) and a sudden change in C base frequency at the polyA site (0 nt); B is a motif diagram of 60 nt upstream of the polyA site, it is expected that the main motif should be A[A / T]TAAA; C is a pie chart of the number of polyA sites at different genomic locations, it is expected that the number of polyA sites in the PA, 3'UTR, and 3'UTR extension regions of the gene should be relatively large, accounting for at least 30% cumulatively; D is the overlap between the polyA sites detected by this method (scPolyA-seq2) and the polyA sites in the public database polyA-DB2, the distance is set to less than 12 nt as overlap, it is expected that there should be at least 30% overlap; E is the nearest distance distribution diagram between the polyA sites detected by this method (scPolyA-seq2) and the polyA sites in the public database polyA-DB2, the expected distance should be mainly near 0.

[0022] Figure 3 This is a schematic diagram of spurious alignments in Example 1. A represents a spurious polyA site obtained from a spurious alignment on the positive strand, which is a polyT region on the genome; the red arrow indicates the gene direction. B represents the dominant matching sequence in the positive strand spurious alignment, which is the polyT sequence in the bam file. C represents a spurious polyA site obtained from a spurious alignment on the negative strand, which is a polyA region on the genome; the red arrow indicates the gene direction. D represents the dominant matching sequence in the negative strand spurious alignment, which is the polyA sequence in the bam file.

[0023] Figure 4 This is a graph showing the polyA sites of the CD47 gene sample in Example 1. In Figure A, the graph shows the polyA sites and polyA supporting read peaks on the CD47 gene, with one row representing one sample. The bottom row, CD4 PBMC, shows the polyA sites detected by the method of this invention, and polyA DB2 (hg38) shows the polyA sites in the publicly available database polyA_DB2. Row GRCh38 represents the 3'UTR (Untranslated Regions) region of the CD47 gene in the reference genome hg38. (B) is a graph showing the polyA sites and polyA supporting read peaks of a single-cell sample, with one row representing one cell. The bottom row, CD4 PBMC, shows the polyA sites detected by the method of this invention, and polyA DB2 (hg38) shows the polyA sites in the publicly available database polyA_DB2. Row GRCh38 represents the 3'UTR (Untranslated Regions) region of a single-cell sample in the reference genome. Detailed Implementation

[0024] The present invention will be further described in detail below with reference to specific embodiments and accompanying drawings.

[0025] A method for analyzing variable polyadenylation (APA) in the 3' end transcriptome of a single cell includes the following steps: (1) Obtaining high-throughput sequencing data: Single cells are separated by microarray or flow cytometry, and amplified using poly(T) primers with cell barcode (CB) and UMI (Unique Molecular Identifier, UMI). Library construction and sequencing are performed, and the sequencing results are split according to CB to obtain sequencing data for each cell and obtain 3' transcriptome data of the cells. The number of cells can be appropriately reduced to increase the sequencing volume of each cell.

[0026] The sequencing results were split according to CB to obtain the sequencing data of the 3' transcriptome of each cell.

[0027] (2) Screening polyA support reads and aligning and obtaining polyA site coordinates: Screen polyA support reads from the 3' transcriptome sequencing data of each cell in step (1), such as reads containing 10 consecutive A bases or 6 consecutive A bases at the end. Remove these A bases and the bases after them to obtain polyA support reads. Align the selected polyA support reads to the reference genome using alignment software. After alignment, use parameters such as setting high MAPQ and unique alignment to screen and retain high-quality alignments, remove false alignments, and then use unique molecular identifiers to deduplicate the polymerase chain reaction in the sequencing data of each cell to obtain the Bam file of the sequencing read coordinates of each cell. Convert the Bam file to a bed file, merge the bed files of all cells, and obtain the 3' coordinates of each read through the bed file, which is the original polyA coordinates. Perform polyA site merging, filtering, and verification to obtain the final polyA site. If the polyA site verification step fails, return to the alignment step and continue to fine-tune the parameters.

[0028] Spurious alignments are removed by filtering out polyA or poly(T) dominated alignments in the alignment file. These spurious alignments can create false polyA sites, misleading subsequent analyses. The method for removing spurious alignments is to ensure that the sequence corresponding to M in the CIGAR string is the one that performs the matching function in the sequencing reads aligned to the reference genome, where consecutive polyA or polyT sequences account for more than 50% of the alignments. Specifically, the verification of polyA sites includes the following steps: s1: Merging, filtering, and annotation of polyA sites: All BED files containing 3' end positions of cells are merged into one BED file, sorted by chromosome classification, and the 3' end positions of each read are obtained, i.e., the original polyA sites; original polyA sites clustered less than 12 nt are merged and the merging process is recorded. Among the polyA sites merged into a cluster, the genomic coordinates of the highest read are used as the coordinates of that polyA site, and the sum of all reads in the cluster is used as the read count of that polyA site; internal amplification is removed, i.e., amplification of A regions on the genome that may not be caused by polyA tails is removed. For example, polyA sites with 6 consecutive A bases or more than 15 A bases in 20 bases downstream of the polyA site can be filtered out; polyA sites with few supporting cells, such as polyA sites supported by less than 1% of cells, are filtered out to obtain the final polyA sites; the final polyA sites are annotated according to the reference genome. The annotation includes basic information such as the gene name, gene region, or upstream and downstream base sequences; s2, polyA site verification: The polyA sites obtained in step s1 are verified using at least three different verification methods to determine the reliability of the final polyA sites.

[0029] The main verification methods include: Plot the base frequencies 100 nt upstream and downstream of the polyA site; Plot a polyA signal enrichment map of nucleotides 60 nt upstream of polyA; Plot the polyA count distribution at each genomic location; Alternatively, compare the number of overlaps with sites in existing publicly available polyA databases; Alternatively, compare the distribution of relative distances with sites in the existing public polyA database.

[0030] The criteria for judging the reliability of the final polyA site in the verification method are as follows: (A) The judgment criteria based on the base frequency map 100nt upstream and downstream of the polyA site are: the base frequency 100nt upstream and downstream of the polyA site is expected to have an A base peak near the upstream (-20nt) and a sudden change in C base frequency at the polyA site (0nt). (B) The criteria for judging the polyA signal enrichment map of nucleotides 60 nt upstream of polyA is: the motif map of 60 nt upstream of polyA site, the expected main motif should be A[A / T]TAAA; (C) The criteria for judging the polyA number distribution map of each genomic location are: pie chart of the number of polyA sites at different genomic locations. It is expected that the number of polyA sites in the 3'UTR region of the gene should be relatively large, for example, the proportion of sites falling in PA, 3'UTR, and 3'UTR extension region exceeds 30%; (D) The criterion for comparing the number of overlaps with existing polyA sites in the public polyA database is: the overlap between the detected polyA sites and polyA sites in the public polyA-DB2 database. A distance of less than 12 nt is considered an overlap, and it is expected that there should be at least 30% overlap. (E) The criterion for comparing the distribution of relative distances between the detected polyA sites and the polyA sites in the publicly available polyA database is: the main peak of the expected distance distribution map of the detected polyA sites and the polyA sites in the publicly available polyA-DB2 database should be mainly near 0.

[0031] (3) Quantification of single-cell APA data: After merging all cells obtained in step (2) and verifying reliable polyA site locations and merge records, count the number of polyA support reads at polyA sites in each single-cell bed file, and generate a sequencing read matrix with polyA sites as rows, cells as columns, and the support count value as the content value. (4) Quantification of gene distal usage: The proportion of distal polyA usage for each gene is calculated to obtain the analysis results. Specifically, it is an index weighted by position to the percentage of the final polyA read in the gene, reflecting the difference in polyA site usage status at the single-cell level for each gene. This analysis result is used to characterize the selective polyadenylation differences between different cell types or cell states. Cell states include cell activation, cell differentiation, or disease-related states.

[0032] (5) The snakemake script is used to automate the process and improve data processing efficiency. The BAM and BED files generated during the analysis are loaded, and a custom igv.js local server is used to visualize the polyA site peaks at the single-cell level.

[0033] The detailed flowchart is as follows: Figure 1 As shown.

[0034] Example 1

[0035] A method for analyzing variable polyadenylation (APA) in the 3' end transcriptome of a single cell, the specific steps of which are as follows: Taking the CD47 gene levels in single-cell samples during antibody activation as an example (1) Obtaining high-throughput sequencing data. To obtain high-quality data, flow cytometry was used to separate single cells from the T cell samples into the wells of a 384-well plate. The cells were lysed, and cDNA was reverse transcribed. The cDNA was then amplified using polyT primers with cell barcodes (CB) and unique molecular identifiers (UMI) to construct a library for sequencing, obtaining the 3' transcriptome data of each cell. The number of cells can be appropriately reduced to increase the sequencing throughput per cell. The sequencing results were split according to CB to obtain the 3' transcriptome sequencing data of each cell.

[0036] Using the CD47 gene as an example, the sequencing read counts of different polyA sites of this gene are visualized, reflecting the differences in the utilization rate of different polyA sites of this gene at different antibody activation time points.

[0037] (2) Screening for polyA-supporting reads. Screen for polyA-supporting reads from the sequencing data, such as reads containing 10 consecutive A bases or 6 consecutive A bases at the end, and remove these A bases and the bases following them.

[0038] (3) Align and obtain polyA site coordinates. Align the selected polyA support reads to the reference genome, for example, using STAR alignment software. Only retain high-quality alignments, such as setting a high MAPQ and unique alignment conditions. For STAR alignment results, you can use samtools to filter the data by setting the parameter "-q 255 -F 2820".

[0039] Remove spurious alignments: In the alignment file, filter out alignments dominated by polyA or polyT, as these spurious alignments can create false polyA sites and mislead subsequent analysis. Specifically, the matching indicator is the sequence corresponding to M in the CIGAR string within the aligned genome sequencing reads, where the proportion of consecutive polyA or polyT sites is high (e.g., exceeding 50%). Examples of spurious polyA site reads and the reference genome are shown below. Figure 3 As shown; the polyA peak diagram of the CD47 gene at the sample and single-cell levels is shown in the figure. Figure 4 As shown.

[0040] PCR duplicates (polymerase chain reaction, PCR) were removed from the sequencing data of each cell using UMI. The resulting cell reads were then aligned to a Bam file and converted to a BED file.

[0041] (4) Merging and filtering of polyA sites. Merge the 3' end positions of all cells into one file, separate them by chromosome, and sort them. Obtain the 3' end position of each read, i.e., the original polyA site.

[0042] Merge closely spaced polyA sites, such as those clustered within 12 nt. Record the merging process; among the polyA sites merged into a cluster, the genomic coordinate of the highest read is used as the coordinate of that polyA site, and the sum of all reads in that cluster is used as the read count of that polyA site.

[0043] Remove internal amplification: Remove amplification in the A region of the genome that may not be caused by polyA tails. That is, filter out polyA sites downstream of the polyA site that have a lot of A bases, such as sequences with 6 consecutive A bases downstream of the polyA site or more than 15 A bases in the 20 bases downstream of the polyA site. The parameters can be adjusted according to the subsequent polyA site verification results and the number of polyA bases, and the parameters can be set to be more rigorous or more lenient.

[0044] Sites with low supporting cell counts, such as polyA sites supported by less than 1% of cells, are filtered out to obtain the final polyA sites.

[0045] (5) PolyA site annotation. Based on the reference genome, annotate the obtained polyA sites with information such as the gene name, gene region, and upstream and downstream base sequences.

[0046] (6) PolyA site verification. Multiple methods were used to verify the obtained polyA sites, such as: Plot the base frequencies 100 nt upstream and downstream of the polyA site; Plot a polyA signal (also called motif) enrichment map of the 60 nt upstream of polyA (nucleotide, nt); Plot the polyA count distribution at each genomic location; The number of overlaps between sites in the existing public polyA database and sites detected by this method is compared. The distribution of relative distances between sites in the existing public polyA database and sites detected by this method is compared.

[0047] The criteria for judging the reliability of the final polyA site using several verification methods are as follows: (A) The judgment criteria based on the base frequency map 100nt upstream and downstream of the polyA site are: the base frequency 100nt upstream and downstream of the polyA site is expected to have an A base peak near the upstream (-20nt) and a sudden change in C base frequency at the polyA site (0nt). (B) The criteria for judging the polyA signal enrichment map of nucleotides 60 nt upstream of polyA is: the motif map of 60 nt upstream of polyA site, the expected main motif should be A[A / T]TAAA; (C) The criteria for judging the polyA number distribution map of each genomic location are: pie chart of the number of polyA sites at different genomic locations. It is expected that the number of polyA sites in the 3'UTR region of the gene should be relatively large, for example, the proportion of sites falling in PA, 3'UTR, and 3'UTR extension region exceeds 30%; (D) The criterion for comparing the number of overlaps with existing polyA sites in the public polyA database is: the overlap between the detected polyA sites and polyA sites in the public polyA-DB2 database. A distance of less than 12 nt is considered an overlap, and it is expected that there should be at least 30% overlap. (E) The criterion for comparing the distribution of relative distances between the detected polyA sites and the polyA sites in the publicly available polyA database is: the main peak of the expected distance distribution map of the detected polyA sites and the polyA sites in the publicly available polyA-DB2 database should be mainly near 0.

[0048] If the indicators do not meet expectations, the parameters need to be adjusted and the analysis process repeated. The verification indicator graph is shown below. Figure 2 As shown. Generally, the first two conditions must be met, that is, the judgment criteria (A) and (B) must be met.

[0049] (7) Quantification of single-cell APA data. Using the BED file of each single cell, for the final polyA site location and merged records, the number of polyA-supported reads for each polyA site is counted, and a sequencing read matrix is ​​generated with polyA sites as rows, cells as columns, and the supported count value as the content value.

[0050] (8) Quantification of distal polyA usage. The distal polyA usage ratio of each gene is calculated to obtain the analysis results. Specifically, it is an index weighted by the percentage of polyA reads in the gene according to position, reflecting the difference in polyA site usage status of each gene at the single-cell level. This analysis result is used to characterize the selective polyadenylation differences between different cell types or cell states, including cell activation, cell differentiation, or disease-related states.

[0051] (9) Process automation and visualization. The snakemake script is used to automate the execution of the process and improve data processing efficiency.

[0052] Load the BAM and BED files generated during the analysis process, and use a custom igv.js local server to visualize the polyA site peaks at the single-cell level.

[0053] The technical solutions provided by the embodiments of the present invention have been described in detail above. Specific examples have been used to illustrate the principles and implementation methods of the embodiments of the present invention. The descriptions of the embodiments above are only for helping to understand the principles of the embodiments of the present invention. At the same time, for those skilled in the art, there will be changes in the specific implementation methods and application scope based on the embodiments of the present invention. Therefore, the content of this specification should not be construed as a limitation of the present invention.

Claims

1. A method for analyzing variable polyadenylation of the 3' end transcriptome in single cells, characterized in that: Includes the following steps: (1) Obtaining high-throughput sequencing data: Single cells are separated by microarray or flow cytometry. The nucleic acids of single cells are amplified by polyT primers, libraries are constructed and sequenced. The sequencing results are split to obtain the sequencing data of the 3' transcriptome of each cell. (2) Screening polyA support reads and aligning and obtaining polyA site coordinates: Screening polyA support reads from the sequencing data in step (1); Aligning the selected polyA support reads to the reference genome using alignment software; After alignment, retaining high-quality alignments, removing false alignments, and then using unique molecular identifiers to remove duplicates caused by polymerase chain reaction in the sequencing data of each cell, obtaining the Bam file of the alignment coordinates of the sequencing reads of each cell, converting the Bam file to a bed file, obtaining the coordinates of the 3' end of the sequencing reads, i.e., the polyA site, through the bed file, merging, filtering and verifying polyA sites, and using the verified reliable polyA sites as the final polyA sites; The removal of false alignments is as follows: in the alignment file, alignments dominated by polyA or polyT are filtered out. That is, in the sequencing reads aligned to the reference genome, the sequence corresponding to M in the CIGAR string is used for matching, and the alignments in which consecutive polyA or polyT accounts for more than 50% are excluded. (3) Quantification of single-cell APA data: After merging all cells obtained in step (2) and verifying reliable polyA site locations and merge records, count the number of polyA support reads at polyA sites in each single-cell bed file, and generate a sequencing read matrix with polyA sites as rows, cells as columns, and the support count value as the content value. (4) Quantification of gene distal usage: The proportion of distal polyA sites used for each gene is calculated to obtain the analysis results. Specifically, it is an index weighted by the percentage of the final polyA reads in the gene according to the position, reflecting the difference in polyA site usage status of each gene at the single-cell level. The analysis results are used to characterize the differences in polyadenylation between different cell types or cell states; (5) Use the snakemake script to automate the execution of the process and improve data processing efficiency.

2. The method for analyzing variable polyadenylation of the 3' end transcriptome of a single cell according to claim 1, characterized in that: The polyT primers in step (1) are primers with cell tags and unique molecular identifiers.

3. The method for analyzing variable polyadenylation (APA) in the 3' end transcriptome of a single cell according to claim 1, characterized in that: The process of screening polyA-supporting reads in step (2) involves screening for sequencing reads containing 10 consecutive A's or ending with 6 consecutive A's, and removing these A's and downstream bases.

4. The method for analyzing variable polyadenylation of the 3' end transcriptome in a single cell according to claim 1, characterized in that: The comparison software used in step (2) is selected from STAR, HISAT2, TopHat2 and TopHat.

5. The method for analyzing variable polyadenylation of the 3' end transcriptome of a single cell according to claim 1, characterized in that: The high-quality alignment described in step (2) is the only alignment with a high MAPQ value.

6. The method for analyzing variable polyadenylation of the 3' end transcriptome of a single cell according to claim 1, characterized in that: Step (2) of polyA site filtering includes the following steps: s1: Merging, filtering, and annotation of polyA sites: All BED files containing 3' end positions of cells are merged into one BED file, sorted by chromosome classification, and the 3' end positions of each read are obtained, i.e., the original polyA sites; original polyA sites clustered less than 12 nt are merged and the merging process is recorded. Among the polyA sites merged into a cluster, the genomic coordinate of the highest read is used as the coordinate of the polyA site, and the sum of all reads in the cluster is used as the read count of the polyA site; internal amplification is removed; polyA sites with less than 5% supporting cell count are filtered out; the final polyA sites are obtained; the final polyA sites are annotated according to the reference genome. The annotation content includes: gene name, gene region or upstream and downstream base sequences; s2, polyA site verification: The polyA sites obtained in step s1 are verified using at least three different verification methods to determine the reliability of the final polyA sites.

7. The method for analyzing variable polyadenylation of the 3' end transcriptome of a single cell according to claim 6, characterized in that: The verification method described in step s2 is selected from: Plot the base frequencies 100 nt upstream and downstream of the polyA site; Alternatively, a polyA signal enrichment map of nucleotides 60 nt upstream of polyA can be plotted; Or plot the polyA number distribution at each genomic location; Alternatively, compare the number of overlaps with sites in the existing public polyA database; Alternatively, compare the distribution of relative distances with sites in the existing public polyA database.

8. The method for analyzing variable polyadenylation of the 3' end transcriptome of a single cell according to claim 6, characterized in that: Step s1, which removes internal amplification, means removing A regions on the genome that may not be caused by polyA tails. Specific filtering conditions can be: filtering out polyA sites downstream of polyA sites where there are more than 15 A bases in 6 consecutive A bases and / or 20 bases.

9. The method for analyzing variable polyadenylation of the 3' end transcriptome of a single cell according to claim 1, characterized in that: The cell states described in step (4) include cell activation, cell differentiation, or disease-related states.

Images