Methods for annotation of genome structural variation

Abstract

The present disclosure relates to methods for annotating genomic structural variations and predicting the formation of fusion proteins. The present disclosure also provides systems, devices and computer-readable media for annotating genomic structural variations and predicting the formation of fusion proteins.

Description

Technical Field

[0001] The present invention belongs to the field of bioinformatics, and particularly relates to a method and system for annotating structural variations (SVs) in genomes. Technical Background

[0002] Structural variants (SVs) are the most common form of variation in the human genome. Structural variations include deletion, duplication, insertion, inversion, translocation, and more complex forms of variation. Among germline mutations, structural variation polymorphisms are related factors to some common diseases. Among systemic mutations, structural variations that usually appear in more complex forms are related to the development of cancer.

[0003] Structural variation is a complex form of variation. In addition to the above-mentioned deletion, amplification, insertion, inversion, translocation and other diverse variation types, it may occur within genes or between genes. It may span tens of bp to tens of megabp on the same chromosome, or it may be a translocation between chromosomes. The rich variation forms make it very difficult to annotate and interpret structural variation. At present, the interpretation of structural variation mostly stays on the annotation of genes or their structures. It is impossible to reasonably and effectively annotate and predict clinical information such as its mutation type, speculated impact, whether it is speculated to form a fusion protein, the structural composition of the fusion protein, the prediction of functional impact, HGVS and pathogenicity rating. This has greatly restricted the application of structural variation in the field of precision medicine.

[0004] In summary, the field currently lacks methods and systems that can effectively annotate and interpret structural variations, especially structural variations in the human genome. Summary of the invention

[0005] Listed below are some of the terms and their definitions used in this disclosure.

[0006] NGS: High-Throughput Sequencing, also known as Next Generation Sequencing (NGS), is relative to traditional Sanger Sequencing.

[0007] FASTQ format file: a text format that saves biological sequences (usually nucleic acid sequences) and their sequencing quality score information.

[0008] SAM format file: SAM (Sequence Alignment / Map format) file is a sequence alignment format standard developed by Sanger. It is a text format with TAB as the delimiter.

[0009] BAM format file: A binary file of SAM format file.

[0010] BED format file: BED (Browser Extensible Data) file is a format that defines track feature information such as comments. It is a text format with TAB as the delimiter.

[0011] VCF format file: VCF (Variant Call Format) file is a text format for saving gene sequence variations.

[0012] HGVS: A standard for describing sequence variations developed by the Human Genome Variation Society.

[0013] SV: Structural Variation.

[0014] This method provides a complete method and system for annotation of human genome structural variation and structure and function prediction. In addition to annotation of gene name and structure, this method can further include annotation of mutation type, putative impact, prediction of whether fusion protein is formed, structural composition of fusion protein, functional impact prediction, HGVS annotation and pathogenicity rating and other clinical information to help relevant professionals evaluate the impact of structural variation.

[0015] In one aspect, the present disclosure provides a method for annotating a genomic structural variation (SV), the method comprising:

[0016] a. Providing SV information, the information including the type, orientation and two breakpoint coordinates of the SV, wherein the SV type includes deletion, duplication, inversion and translocation;

[0017] b. annotating the gene, transcript, and position of each breakpoint in the transcript, as well as the direction information of the transcript, the minimum exon number and the maximum exon number of the CDS according to the breakpoint coordinates, or annotating the breakpoint as an intergenic region when the breakpoint is located in the region between genes; and

[0018] c. Annotate the mutation type produced by the SV based on the SV type and the annotation of each breakpoint in step b.

[0019] In some embodiments of the above method, step a includes reading the analysis results of the upstream SV detection software to extract SV information. In a preferred embodiment, the SV variation detection software is selected from Delly, Manta, LUMPY, GRIDSS and markSV.

[0020] In some embodiments, step b comprises:

[0021] b1. annotating the breakpoints as genes or intergenic regions at the positions according to the breakpoint coordinates;

[0022] b2. Determine all transcripts of the genes annotated in b1 and the exons, introns, CDS and UTR information of the transcripts;

[0023] b3. annotating the breakpoint as being located in the UTR region, exon M or intron N according to the position of the breakpoint in the transcript, wherein M is the exon number of the CDS region where the breakpoint is located, and N is the exon number of the CDS with the smaller exon number when the breakpoint is between two CDS regions; and

[0024] b4. Obtain the direction information (+ / - strand information), the minimum exon number and the maximum exon number of the CDS of the transcript of the gene.

[0025] If there are multiple transcripts of a gene, a specified transcript can be selected for annotation. If no transcript is specified, the longest transcript of the gene can be selected for annotation.

[0026] In some embodiments of the above method, annotating the mutation type includes annotating an SV whose two breakpoints are not in the same gene, including:

[0027] If both breakpoints are in the intergenic region, the mutation type is annotated as an intergenic variant;

[0028] If at least one breakpoint was in a gene region, the mutation type was annotated as gene fusion.

[0029] In some embodiments of the above method, annotating the mutation type comprises annotating an SV in which two breakpoints are located in the same gene and the SV type is a deletion, comprising:

[0030] (1) If the two breakpoints are on the same exon, the mutation type is annotated as a splicing region variation or a protein-altering variation based on the distance between at least one breakpoint and the edge of the exon; preferably, the protein-altering variation is annotated as a reading frame deletion or a frameshift variation based on whether the number of missing nucleotides is a multiple of 3 nucleotides;

[0031] (2) If the two breakpoints are on the same intron, the mutation type is annotated as a splice site variation, a splice region variation or an intron variation according to the distance between at least one breakpoint and the upstream and downstream edges of the intron; preferably, the splice site variation is annotated as a splice acceptor variation or a splice donor variation according to the distance between at least one breakpoint and the upstream edge of the intron or the downstream edge of the intron;

[0032] (3) If at least one breakpoint is in the UTR region and the other breakpoint is in the upstream and downstream intergenic regions of the same UTR, the mutation type is annotated as 5'UTR variation or 3'UTR variation according to the UTR region;

[0033] (4) If one of the two breakpoints is on the 5'UTR and the other is on the 3'UTR, the mutation type is annotated as copy number deletion variation;

[0034] (5) For other cases, the mutation type was annotated as large genomic rearrangement.

[0035] In some embodiments, the above step (1) includes: if the two breakpoints are on the same exon, then the protein alteration variation is annotated as a frame deletion (the number of missing nucleotides is a multiple of 3 nucleotides) or a frameshift variation (the number of missing nucleotides is not a multiple of 3 nucleotides) according to whether the number of missing nucleotides is a multiple of 3 nucleotides; wherein if at least one breakpoint is within 2bp of the inner edge of the exon, the mutation type is covered as a splicing region variation.

[0036] In some embodiments, the above step (2) includes: if the two breakpoints are on the same intron, the mutation type is annotated as an intron variation; if at least one breakpoint is within 2bp of the upstream outer edge of the exon, the mutation type is covered as a splice acceptor variation, if at least one breakpoint is within 2bp of the downstream outer edge of the exon, the mutation type is covered as a splice donor variation, and if at least one breakpoint is between 2-20bp of the outer edge of the exon, the mutation type is covered as a splice region variation.

[0037] In some embodiments of the above method, annotating the mutation type comprises annotating an SV in which two breakpoints are located in the same gene and the SV type is a duplicate, comprising:

[0038] (1) If the two breakpoints are on the same exon or on adjacent exons, the mutation type is annotated as a protein-altering variation; preferably, the protein-altering variation is annotated as a frame duplication or a frameshift variation based on whether the number of repeated nucleotides is a multiple of 3 nucleotides;

[0039] (2) If one of the two breakpoints is in an intron and the other is in the same intron or in an exon or UTR adjacent to the intron, the mutation type is annotated as an intronic variant;

[0040] (3) If at least one breakpoint was in the UTR region and the other breakpoint was in the upstream and downstream intergenic regions of the same UTR region, the mutation type was annotated as 5'UTR variation or 3'UTR variation according to the UTR region;

[0041] (4) If one of the two breakpoints is on the 5'UTR and the other is on the 3'UTR, the mutation type is annotated as copy number duplication variation;

[0042] (5) For other cases, the mutation type was annotated as large genomic rearrangement.

[0043] In some embodiments of the above method, annotating the mutation type comprises annotating an SV in which two breakpoints are located in the same gene and the SV type is an inversion, comprising:

[0044] (1) If the two breakpoints are on the same exon, the mutation type is annotated as a splice region variation or a protein-altering variation based on the distance between at least one breakpoint and the edge of the exon; preferably, the protein-altering variation is annotated as a gain-of-terminus or a missense variation based on whether there is a terminator in the coding frame after the mutation;

[0045] (2) If the two breakpoints are on the same intron, the mutation type is annotated as a splice site variation, a splice region variation or an intron variation according to the distance between at least one breakpoint and the upstream or downstream edge of the exon; preferably, the splice site variation is annotated as a splice acceptor variation or a splice donor variation according to the distance between at least one breakpoint and the upstream edge of the exon or the downstream edge of the exon;

[0046] (3) If at least one breakpoint was in the UTR region and the other breakpoint was in the upstream and downstream intergenic regions of the same UTR region, the mutation type was annotated as 5'UTR variation or 3'UTR variation according to the UTR region;

[0047] (4) For other cases, the mutation type was annotated as large genomic rearrangement.

[0048] In some embodiments, the above step (1) includes: if the two breakpoints are on the same exon, the mutation type is annotated as termination gain (with terminator) or missense variation (without terminator) according to whether there is a terminator in the coding frame after the mutation; if at least one breakpoint is within 2bp of the inner edge of the exon, the mutation type is covered as splice region variation.

[0049] In some embodiments, the above step (2) includes: if the two breakpoints are on the same intron, the mutation type is annotated as an intron variation; if at least one breakpoint is within 2bp of the upstream outer edge of the exon, the mutation type is covered as a splice acceptor variation, if at least one breakpoint is within 2bp of the downstream outer edge of the exon, the mutation type is covered as a splice donor variation, and if at least one breakpoint is between 2-20bp of the outer edge of the exon, the mutation type is covered as a splice region variation.

[0050] In some embodiments, the methods of the present disclosure further comprise annotating the putative impact of the mutation type.

[0051] In some embodiments, copy number duplication variation is annotated as having a medium predicted impact, and copy number deletion variation and large genomic rearrangement are annotated as having a high predicted impact. In some embodiments, in addition to the above, the predicted impact of other mutation types is annotated according to the prediction criteria of the snpEff software.

[0052] In some embodiments, the prediction criteria of the snpEff software are as follows:

[0053]

[0054] In some embodiments, the method of the present disclosure further comprises annotating a new fusion sequence generated by an SV whose mutation type is gene fusion, comprising:

[0055] a. Construct fusion sequences according to SV direction, including:

[0056] For SVs with FR (forward-reverse) orientation, the sequences on the two genome replies are connected in order; for SVs with RF (reverse-forward) orientation, the sequences on the two gene replies are connected in reverse order; for SVs with FF (forward-forward) orientation, the sequences on the downstream gene replies are reverse-complemented and connected to the upstream gene in order; for SVs with RR (reverse-reverse) orientation, the sequences on the upstream gene replies are reverse-complemented and connected to the downstream gene in order;

[0057] b. According to the composition relationship and transcript direction of the two genes in the fusion sequence, the fusion sequence structure is determined as follows:

[0058] For genes with fusion sequences upstream of the genome, if the transcript chain direction is forward, the fusion sequence retains the exon sequence upstream of the exon where the breakpoint is located; if the transcript chain direction is reverse, the fusion sequence retains the exon sequence downstream of the exon where the breakpoint is located;

[0059] For genes with fusion sequences downstream of the genome, if the transcript chain direction is forward, the fusion sequence retains the exon sequence downstream of the exon where the breakpoint is located; if the transcript chain direction is reverse, the fusion sequence retains the exon sequence upstream of the exon where the breakpoint is located;

[0060] c. According to the composition relationship and transcription direction of the two genes in the fusion sequence, there are four situations: (1) ++, i.e., the 5' end of the fusion sequence is the left side of the fusion sequence, and the 3' end is the right side of the fusion sequence; (2) --, i.e., the 3' end of the fusion sequence is the left side of the fusion sequence, and the 5' end is the right side of the fusion sequence; (3) +-, i.e., the transcription directions of the two genes in the fusion sequence conflict, and the 5' end and 3' end cannot be determined; and (4) -+, i.e., the transcription directions of the two genes in the fusion sequence conflict, and the 5' end and 3' end cannot be determined.

[0061] Among them, situations (1) and (2) are predicted to be able to form fusion proteins, and situations (3) and (4) are predicted to be unable to form fusion proteins.

[0062] In some embodiments, the methods of the present disclosure further comprise annotating the functional impact of the SV, comprising:

[0063] For intergenic SVs, if one of the two genes is a proto-oncogene and predicted to produce a fusion protein, the functional impact was annotated as loss of function: if one gene is a tumor suppressor gene, the functional impact was annotated as loss of function;

[0064] For SVs within genes, if the gene is a tumor suppressor gene, the functional impact is annotated as loss of function; for SVs within genes known to cause gain of function, the functional impact is annotated as gain of function;

[0065] For other cases the functional impact is annotated as None.

[0066] In some embodiments, the intragenic SV known to cause gain of function may be selected from BRAF exon11-18 amp, EGFR 18-25 amp, MET exon15-21 amp and BRAF exon5-9 del. In addition to the above, other intragenic SVs that cause gain of function are known in the art.

[0067] In some embodiments, the methods of the present disclosure further comprise annotating the pathogenicity level of the SV, wherein:

[0068] For intergenic SVs, if they were annotated as loss-of-function, the evidence annotation was PVS1 and the rating annotation was potentially pathogenic;

[0069] For SVs within genes, the level of evidence and rating were annotated with reference to the ACMG CNV guidelines (e.g., 2019 ACMG CNV guidelines).

[0070] In some embodiments, the method of the present disclosure further comprises generating annotations for the HGVS data structure based on the type and breakpoint coordinates of the SV. In some embodiments, the HGVS data structure may comprise a g.HGVS, c.HGVS and / or p.HGVS data structure.

[0071] In another aspect, the present disclosure relates to a system for annotating genomic structural variations (SVs), the system comprising:

[0072] SV information reading module; and

[0073] SV Annotation Module,

[0074] The SV information reading module is configured to read SV information, wherein the information includes the type, orientation and two breakpoint coordinates of the SV, wherein the SV type includes deletion, duplication, inversion and translocation; and the SV annotation module is configured to annotate the SV using the method described in the first aspect of the present disclosure.

[0075] In one aspect, the present disclosure relates to an apparatus for annotating genomic structural variations (SVs), comprising:

[0076] memory for storing computer program instructions; and

[0077] a processor for executing computer program instructions,

[0078] When the computer program instructions are executed by the processor, the device performs the method described in the first aspect of the present disclosure.

[0079] In another aspect, the present disclosure relates to a computer-readable medium storing computer program instructions, wherein when the computer program instructions are executed by a processor, the method described in the first aspect of the present disclosure is implemented. DETAILED DESCRIPTION

[0080] Example 1. Exemplary Structural Variation (SV) Detection Methods

[0081] The following uses the markSV method as an example to illustrate a method for obtaining information related to structural variations (SVs) in a genomic region.

[0082] 1. Generation of sequence alignment files: After the library preparation, the samples to be tested are sequenced on the Illumina sequencing platform to generate FASTQ files. After quality control of the FASTQ files, the alignment software BWA-MEM is used to align the FASTQ files with the human reference genome (hg19 / b37) and generate SAM files. The SAM files are converted to BAM files using the samtools software, and the BAM files are used as input files for subsequent testing.

[0083] 2. Calculation of insert length outliers: Read the reads in the BAM file to estimate the parameters of the insert length distribution and the threshold of the outlier. If the BAM file contains more than 1 million reads, 1 million are randomly selected for calculation to reduce the time of calculating the insert length distribution parameters. The insert length distribution defaults to log-normal distribution. The parameters μ and σ of the insert length distribution are calculated and the distribution outlier threshold is obtained according to the 3sigma principle. If the insert length of a pair of reads exceeds μ+3σ, it is considered an abnormal insert length.

[0084] 3. Read the BAM file, filter out reads with MapQ less than 10 and reads located in long repeat regions (homopolymer sequences or microsatellite sequences longer than 45 bp), and retain the remaining valid reads to proceed to the next step.

[0085] 4. Extract SV-related signals from valid reads: (1) Reads containing soft-clipping are treated as SR signals and enter the SR signal analysis process; (2) Reads whose reply direction is not forward-reverse are treated as PE signals and enter the PE signal analysis process; (3) Reads whose reply direction is forward-reverse but whose insert length exceeds the insert length distribution outlier threshold calculated in step 2 above are treated as PE signals and enter the PE signal analysis process.

[0086] 5.SR signal analysis process:

[0087] 5.1 The reads entering the SR signal analysis process are first grouped according to QNAME, and in each group, the primary and secondary replies are paired according to the FLAG value, and the reads that cannot be paired are filtered out;

[0088] 5.2 The paired SR signal read length determines the SV type according to the position and direction information of the main reply and the suboptimal reply, and forms the SR signal cluster after clustering according to the SV type and breakpoint coordinate information;

[0089] 5.3 Count the number (unique) and total number (total) of reads containing different templates in the SR signal cluster according to the start and end coordinates of the read template. Record the QNAME of the reads in the SR signal cluster with unique = 1 and total greater than 3, which will be filtered out in the PE signal analysis process; and

[0090] 5.4 Cluster and merge the two SR signals whose two coordinate fluctuations of the breakpoints are less than or equal to 10 bp.

[0091] 6.PE signal analysis process:

[0092] 6.1 The read pairs entering the PE signal analysis process are first grouped according to QNAME, and R1 and R2 are paired according to the FLAG value in each group to filter out the reads that cannot be paired;

[0093] 6.2 The paired PE signal reads are used to determine the SV type based on the position and direction information of R1 and R2; the reads with the QNAME recorded in the above step 5.4 are filtered out, and the PE signal clusters are formed after clustering based on the SV type and breakpoint coordinate information;

[0094] 6.3 Cluster and merge the two PE signals whose two coordinate fluctuations of the breakpoints are both less than or equal to 10bp.

[0095] 7 Signal integration: SR signal clusters and PE signal clusters are merged according to the breakpoint coordinate information and the QNAME of the read length to form a comprehensive signal cluster, which includes the merged SR signal clusters and PE signal clusters, or the individual SR signal clusters or PE signal clusters that failed to be merged.

[0096] 8. Filters:

[0097] 8.1 Evaluate the similarity of the integrated signal clusters and filter out the integrated signal clusters with high similarity, as follows: extract the reference genome sequence of 220 bp upstream and downstream of the two coordinates of the breakpoint, and use BLASTn to compare them. If the identity exceeds 80% and the comparison length exceeds 50 bp, it is considered to be highly similar;

[0098] 8.2 Evaluate the repetitive regions for the integrated signal clusters and filter out the integrated signal clusters located in the repetitive regions, as follows: extract the read sequences at the two coordinates of the breakpoints of the integrated signal clusters, and identify whether the sequences contain homopolymers (1 base sequence repeated more than 10 times), microsatellites (2-6 consecutive base sequences repeated more than 4 times), or a combination of the two; if the read sequence at at least one coordinate contains homopolymer sequences and / or microsatellite sequences, it is considered to be located in the repetitive region;

[0099] 8.3 Evaluate the chimeras for the integrated signal clusters and filter out the integrated signal clusters with chimeras, as follows: extract the soft-clipping sequence and the starting coordinates of the soft-clipping reads at each breakpoint coordinate, and group the reads with both the soft-clipping sequence and the starting coordinates of the soft-clipping into one group; take the group with the starting coordinates of the soft-clipping consistent with the breakpoint coordinates as the primary group, merge the other groups into the secondary group, and calculate the entropy values ​​of the two groups; if the entropy value on at least one breakpoint coordinate exceeds 0.5, the integrated signal cluster is considered to have chimeras.

[0100] 8.4 Evaluate the simple templates for the comprehensive signal clusters and filter out the comprehensive signal clusters that are simple templates, as follows: extract the fluctuation range of the start and end coordinates of the read template in the comprehensive signal clusters. If the fluctuation range of the start and end coordinates is less than or equal to 3 bp, it is considered to be a simple template;

[0101] 8.5 Evaluate the average mapq of the integrated signal clusters and filter out the integrated signal clusters with low average mapq, specifically as follows: extract the mapq values ​​of the reads at the two coordinates of the breakpoints of the integrated signal clusters, and calculate the average mapq; if the average mapq of the reads at at least one coordinate is less than or equal to the filtering threshold of 25, the integrated signal cluster is considered to have a low average mapq;

[0102] 8.6 Evaluate whether the integrated signal cluster is in the BED file definition area and filter out the integrated signal clusters that are not in the BED file definition area, specifically as follows: extract two coordinates of the integrated signal cluster breakpoints, and if both coordinates are not in the BED file definition area, it is considered that the integrated signal cluster is not in the BED file definition area; if no BED file is specified, this step is not performed;

[0103] 8.7 Filter out the integrated signal clusters with SV size less than 500;

[0104] 8.8 According to the starting coordinates and ending coordinates of the reads in the integrated signal clusters, determine the number of read pairs (unique) and the total number of read pairs (total) containing different templates in each integrated signal cluster, and filter out integrated signal clusters with unique values ​​less than 4.

[0105] 9. Output report: Annotate the comprehensive signal clustering and output the result report.

[0106] Example 2. Structural variation (SV) annotation process

[0107] 1. Read the structural variation result file

[0108] Read and parse the analysis results of upstream structural variation detection software (such as markSV software) and extract the information required for annotation, such as SV type, orientation, breakpoint coordinates (including chromosome number and genome coordinates), etc. Supported formats include custom format tsv files and VCF format files.

[0109] 2. Gene and structure information annotation

[0110] Structural variation usually includes two breakpoints, and it is necessary to annotate gene and structural information for each breakpoint. The annotation of gene and structural information is divided into three steps: gene annotation, transcript annotation, and gene structure annotation.

[0111] 2.1 Gene Annotation

[0112] According to the chromosome number and genome coordinates of the breakpoint, the genes and related information at the location are first retrieved in the GENCODE database. If there is no gene at the breakpoint, it is regarded as an intergenic region. At the same time, the genes upstream and downstream of the breakpoint are retrieved, and the point is annotated as an intergenic region between the upstream gene and the downstream gene in the form of intergenic (upstream gene, downstream gene), and the annotation of gene and structural information is ended at the same time. If there are multiple genes at the breakpoint, genes with protein function are preferred.

[0113] 2.2 Transcript Annotation

[0114] According to the gene ID of the previous gene annotation, all transcripts of the gene ID and their UTR, Exon, CDS and other structural information are retrieved from the GENCODE database. If there is a specified transcript for a gene according to the configuration information, the transcript is directly selected to enter the gene structure annotation. Otherwise, the longest transcript is selected as the default transcript to enter the gene structure annotation.

[0115] 2.3 Gene structure annotation

[0116] According to the transcript selected in the previous step, the UTR, Exon, CDS and other structural information contained in the transcript are traversed. If the breakpoint is in the UTR region, it is considered that the breakpoint is in the UTR region of this gene and this transcript. If the breakpoint is in the CDS region, the exon number of the CDS region is obtained, and the breakpoint is considered to be in the exon N region of this gene and this transcript (N is the exon number). If the breakpoint is between two adjacent CDS regions, the one with the smaller exon number of the two CDS regions is obtained, and the breakpoint is considered to be in the intron N region of this gene and this transcript (N is the smaller exon number). At the same time, the direction information (+ / - chain) of this gene and this transcript, the CDS minimum exon number and maximum exon number information are collected.

[0117] In this way, the gene and structure information annotation obtains the information of the gene, transcript, gene structure, and the direction information of this gene and transcript (+ / - chain), CDS minimum exon number and maximum exon number for use in subsequent processes.

[0118] 3. Mutation type and inferred impact annotation

[0119] The naming of mutation types is mainly implemented according to the framework of SO (Sequence Ontology), and some adjustments are made according to actual conditions.

[0120] 3.1 Structural variation among genes

[0121] For structural variants where the two breakpoints are not in the same gene. If both breakpoints are in the intergenic region, the mutation type is annotated as intergenic variant (SO:0001628). If only one breakpoint is in the intergenic region, the mutation type is annotated as gene fusion (SO:0001565). If both breakpoints are in the gene region, the mutation type is annotated as gene fusion (SO:0001565).

[0122] 3.2 Structural changes within genes different

[0123] For structural variations with two breakpoints in the same gene, specific judgments are made based on the variation type.

[0124] 3.2.1 For the mutation type of Deletion, make the following judgment:

[0125] If the two breakpoints are in the same exon, if the length of the deletion is a multiple of 3 nucleotides, the mutation type is annotated as a frame deletion (SO:0001822), otherwise the mutation type is annotated as a frameshift variant (SO:0001589). Among them, if the breakpoint is 2bp within the edge of the exon (affecting splicing), the mutation type is covered as a splice region variant (SO:0001630).

[0126] If the two breakpoints are on the same intron, the mutation type is annotated as an intron variation (SO:0001627). Among them, if the breakpoint is 2bp outside the upstream edge of the exon (affecting splicing), the mutation type is covered as a splice acceptor variation (SO:0001574); if the breakpoint is 2bp outside the downstream edge of the exon (affecting splicing), the mutation type is covered as a splice donor variation (SO:0001575); if the breakpoint is 2-20bp outside the exon (affecting splicing), the mutation type is covered as a splice region variation (SO:0001630).

[0127] If the two breakpoints are in the same UTR, if it is 5'UTR, the mutation type is annotated as 5'UTR variation (SO:0001623). If it is 3'UTR, the mutation type is annotated as 3'UTR variation (SO:0001624).

[0128] If the two breakpoints are in different exons or introns, or one in an exon and the other in an intron, the mutation type is annotated as large genomic rearrangement.

[0129] If the two breakpoints are located at positions that only affect the UTR region, for example, if one breakpoint is in the 5'UTR region and the other is in an intron or exon on the 5'UTR region (i.e., the upstream and downstream gene regions of the same 5'UTR region), the mutation type is annotated as 5'UTR variation (SO:0001623). If one breakpoint is in the 3'UTR region and the other is in an intron or exon on the 3'UTR region (i.e., the upstream and downstream gene regions of the same 3'UTR region), the mutation type is annotated as 3'UTR variation (SO:0001624).

[0130] If the two breakpoints are one in the 5'UTR and the other in the 3'UTR, the mutation type is annotated as copy number deletion variation.

[0131] For other cases, the mutation type was annotated as large genomic rearrangement.

[0132] 3.2.2 For the variation type of duplication, make the following judgment:

[0133] When the two breakpoints are in the same exon or in two adjacent exons, if the length of the duplication is a multiple of 3 nucleotides, the mutation type is annotated as an in-frame insertion (SO:0001821), otherwise the mutation type is annotated as a frameshift variant (SO:0001589).

[0134] If the two breakpoints are in the same intron, the mutation type is annotated as intronic variation (SO:0001627).

[0135] If the two breakpoints are in the same UTR, if it is 5'UTR, the mutation type is annotated as 5'UTR variation (SO:0001623). If it is 3'UTR, the mutation type is annotated as 3'UTR variation (SO:0001624).

[0136] If the two breakpoints are in two non-adjacent exons, the mutation type is annotated as large genomic rearrangement.

[0137] If the two breakpoints are in two adjacent exons and introns, the mutation type is annotated as intronic variant (SO:0001627).

[0138] If the two breakpoints are located in two non-adjacent exons and introns, the mutation type is annotated as large genomic rearrangement.

[0139] If the two breakpoints are in adjacent UTR and intron, the mutation type is annotated as intronic variation.

[0140] If the two breakpoints are located at positions that only affect the UTR region, for example, if one breakpoint is in the 5'UTR region and the other is in an intron or exon on the 5'UTR region (i.e., the upstream and downstream gene regions of the same 5'UTR region), the mutation type is annotated as 5'UTR variation (SO:0001623). If one breakpoint is in the 3'UTR region and the other is in an intron or exon on the 3'UTR region (i.e., the upstream and downstream gene regions of the same 3'UTR region), the mutation type is annotated as 3'UTR variation (SO:0001624).

[0141] If the two breakpoints are one in the 5'UTR and the other in the 3'UTR, the mutation type is annotated as copy number duplication variation.

[0142] For other cases, the mutation type was annotated as large genomic rearrangement.

[0143] 3.2.3 For the mutation type of inversion, make the following judgment:

[0144] If the two breakpoints are on the same exon, if there is a terminator in the coding frame after the mutation, the mutation type is annotated as termination gain (SO:0001587). Otherwise, the mutation type is annotated as missense mutation (SO:0001583). Among them, if the breakpoint is 2bp within the edge of the exon (affecting splice), the mutation type is covered as splice region mutation (SO:0001630).

[0145] If the two breakpoints are on the same intron, the mutation type is annotated as an intron variation (SO:0001627). Among them, if the breakpoint is 2bp outside the upstream edge of the exon (affecting splicing), the mutation type is covered as a splice acceptor variation (SO:0001574); if the breakpoint is 2bp outside the downstream edge of the exon (affecting splicing), the mutation type is covered as a splice donor variation (SO:0001575); if the breakpoint is 2-20bp outside the exon (affecting splicing), the mutation type is covered as a splice region variation (SO:0001630).

[0146] If the two breakpoints are located at positions that only affect the UTR region, for example, if one breakpoint is in the 5'UTR region and the other is in an intron or exon on the 5'UTR region (i.e., the upstream and downstream gene regions of the same 5'UTR region), the mutation type is annotated as 5'UTR variation (SO:0001623). If one breakpoint is in the 3'UTR region and the other is in an intron or exon on the 3'UTR region (i.e., the upstream and downstream gene regions of the same 3'UTR region), the mutation type is annotated as 3'UTR variation (SO:0001624).

[0147] If the two breakpoints are in different exons or introns, or one in an exon and one in an intron, the mutation type is annotated as large genomic rearrangement (SO:0001572).

[0148] If the two breakpoints were one in the 5'UTR and the other in the 3'UTR, the mutation type was annotated as large genomic rearrangement.

[0149] For other cases, the mutation type was annotated as large genomic rearrangement (SO:0001572).

[0150] For copy number duplication variants, the annotation inferred impact is medium. For copy number deletion variants, the annotation inferred impact is high.

[0151] For mutations of large genomic rearrangement type, the annotation inferred impact is high.

[0152] The inferred effects of the remaining mutation types are annotated according to the prediction criteria of the snpEff software, as follows.

[0153]

[0154] 4. Fusion sequence structure annotation and fusion protein prediction

[0155] For new fusion sequences generated by structural variations between genes, the structure of the fusion sequence is annotated and whether a fusion protein will be formed is predicted to help determine the impact of the structural variation on protein function.

[0156] First, construct the fusion sequence according to the variation type and direction. For intrachromosomal variation types, there are deletions, duplications, inversions, etc. Their corresponding directions are FR (forward-reverse) corresponding to "deletions", RF (reverse-forward) corresponding to "duplications", and FF and RR corresponding to "inversions" (among which, inversions have two types, left-spanning and right-spanning, left-spanning corresponds to the FF direction, and right-spanning corresponds to the RR direction). The directions corresponding to interchromosomal variation types are similar to those within chromosomes, which are also FR, RF, FF, and RR. For the FR type, the sequences on the replies of the two genes are connected in order. For the RF type, the sequences on the replies of the two genes are connected in reverse order. For the FF type, the sequences on the replies of the downstream gene (genome) are reversely complemented and then connected in order to the upstream gene (genome). For the RR type, the sequences on the replies of the upstream gene (genome) are reversely complemented and then connected in order to the downstream gene (genome).

[0157] Then, the structure of the fusion sequence is determined according to the composition relationship and transcription direction of the two genes in the fusion sequence. For the gene on the left side of the fusion sequence, if the gene transcript is +, the sequence at the 5' end of the gene is retained after fusion. If the gene transcript is -, the sequence at the 3' end of the gene is retained after fusion. For the gene on the right side of the fusion sequence, if the gene transcript is +, the sequence at the 3' end of the gene is retained after fusion. If the gene transcript is -, the sequence at the 5' end of the gene is retained after fusion. For the sequence that retains the 5' end, if the breakpoint is on exon N, it is predicted that the exon 1-(N-1) part of the gene is retained; if the breakpoint is on intron N, it is predicted that the intron 1-N part of the gene is retained. For the sequence that retains the 3' end, if the breakpoint is on exon N, it is predicted that the exon (N+1) and the subsequent part of the gene are retained; if the breakpoint is on intron N, it is also predicted that the exon (N+1) and the subsequent part of the gene are retained.

[0158] Finally, determine the transcription direction of the new sequence after fusion. According to the composition relationship and transcription direction of the two genes in the fusion sequence, there are four possibilities: (1). ++, that is, the 5' end of the fusion sequence is the left side of the fusion sequence, and the 3' end is the right side of the fusion sequence. (2). - -, that is, the 3' end of the fusion sequence is the left side of the fusion sequence, and the 5' end is the right side of the fusion sequence, which is opposite to (1). Both of the above situations predict that fusion proteins can be formed. (3). +-, that is, the transcription directions of the two genes in the fusion sequence conflict, and the 5' end and 3' end cannot be determined. (4). - +, that is, the transcription directions of the two genes in the fusion sequence conflict, and the 5' end and 3' end cannot be determined. Both of the above situations predict that fusion proteins cannot be formed.

[0159] 5. HGVS Annotation

[0160] HGVS helps describe and understand structural variation within genes.

[0161] 5.1 Generate g.HGVS

[0162] According to the coordinates of the two breakpoints and the mutation type, g.HGVS can be generated in the form of {refseq}:g.{start}_{end}{var_type}. Among them, refseq is the refseq corresponding to the chromosome number, or the chromosome number can be used directly; start is the coordinate of the smaller breakpoint; end is the coordinate of the larger breakpoint; if it is a deletion, var_type is del, if it is a duplication, var_type is dup, if it is an inversion, var_type is inv.

[0163] 5.2 Generate c.HGVS

[0164] Based on the g.HGVS generated in the previous step, use the position-converter page on mutalyzer to annotate c.HGVS. For a g.HGVS, multiple c.HGVS may be generated due to different genes and transcripts. First select the relevant entries of the gene name annotated in the previous step. If a transcript is specified, use the transcript and its c.HGVS. Otherwise, LRG transcripts and their c.HGVS are preferred. If there is no LRG transcript, transcripts starting with NM or MP and their c.HGVS are preferred.

[0165] 5.3 Generate p.HGVS

[0166] Based on the transcripts and c.HGVS generated in the previous step, use the name-checker page on mutalyzer to annotate p.HGVS. If p.HGVS cannot be annotated, annotate p.HGVS as None. The default p.HGVS is in the three-letter form. At the same time, based on the conversion relationship between the three-letter and single-letter amino acids, a single-letter p.HGVS form is generated for use in subsequent processes.

[0167] The above process of interacting with the mutalyzer web page / interface can be automated using a web crawler to improve efficiency and stability.

[0168] 6. Functional Impact Prediction

[0169] For structural variants between genes, if the partner gene (one of the two genes) is a proto-oncogene and the upstream pipeline predicts the production of a fusion protein, the functional impact is annotated as loss-of-function. If the partner gene is a tumor suppressor gene, the functional impact is also annotated as loss-of-function.

[0170] For intragenic variants, if the gene (the two genes are the same) is a tumor suppressor gene, the functional impact is annotated as loss of function. For intragenic structural variants known to cause gain of function, including BRAF exon11-18 amp, EGFR 18-25amp, MET exon15-21 amp, and BRAF exon5-9 del, the functional impact is annotated as gain of function.

[0171] For other cases the feature impact annotation is None.

[0172] 7.Classification Annotation

[0173] For structural variants between genes, if the upstream pipeline annotated them as loss-of-function, the level of evidence was annotated as PVS1 and the rating was annotated as likely pathogenic.

[0174] For structural variations within genes, the 2019 ACMG CNV Guidelines were mainly used to annotate and judge the level of evidence and rating.

[0175] Example 3. Structural variation annotation

[0176] After a series of steps of library preparation, the samples were sequenced on the Illumina sequencing platform to generate FASTQ files. After quality control of the FASTQ files, the alignment software BWA-MEM was used to align the FASTQ files with the human reference genome (hg19 / b37) and generate SAM files. The SAM files were then converted to BAM files using the samtools software, and the BAM files were used as input files, and the structural variation detection software was used to analyze the structural variation results. The structural variation results were used as input files, and the analysis results were obtained using this method.

[0177] The analysis steps are as described in Example 2. The analysis results are taken as an example of a lung cancer FFPE sample. After the sample was analyzed by the structural variation detection software, there was an interchromosomal translocation (CTX) in the result, the direction was FR, one breakpoint was at chr6:117647625, and the other breakpoint was at chr5:149782845. After analysis by this method, it was concluded that the structural variation was a mutation type of gene fusion, and occurred in the fusion of ROS1 intron 32 (chr6:117647625) and CD74 intron 6 (chr5:149782845). The fusion form is CD74-ROS1 (C6:R33), which indicates that the fusion is the 5'UTR-intron 6 of the CD74 gene connected to the exon 33-3'UTR of the ROS1 gene. And this method predicts that the fusion will form a fusion protein.

[0178] The sample was verified by FISH method and was CD74-ROS1 positive, which showed that the sample had fusion of CD74 and ROS1 at the DNA level. The sample was also verified by IHC method and was CD74-ROS1 positive, which showed that the sample also had fusion of CD74 and ROS1 at the protein level. Both of these verification results are consistent with the prediction results of this method. In addition, this method also provides other annotations and prediction information: the inferred impact is high, the functional impact is predicted to be loss of function, the evidence level is PVS1, and the pathogenicity rating is possibly pathogenic.

[0179] This method not only annotates the gene and its structure for structural variation, but also annotates or predicts its mutation type, inferred impact, whether it is predicted to form a fusion protein, the structural composition of the fusion protein, the prediction of functional impact, HGVS, pathogenicity rating and other clinically concerned information. Thus, the basic information of structural variation and clinically relevant information are annotated or predicted from multiple perspectives, providing a more comprehensive and in-depth interpretation for clinical applications.

Claims

1. A method for annotating a genomic structural variation (SV), the method comprising: a. Providing SV information, the information including the type, direction and two breakpoint coordinates of the SV, wherein the SV type includes deletion, duplication, inversion and translocation; b. annotating the gene, transcript, position of the breakpoint in the transcript, direction information of the transcript, minimum exon number and maximum exon number of the CDS according to the breakpoint coordinates, or annotating the breakpoint as an intergenic region when the breakpoint is located in a region between genes; step b comprises: b1. annotating the breakpoints as genes or intergenic regions at the positions according to the breakpoint coordinates; b2. Determine all transcripts of the genes annotated in b1 and the exons, introns, CDS and UTR information of the transcripts, wherein for genes with specified transcripts, the transcripts are selected for annotation, and for genes without specified transcripts, the longest transcripts are selected for annotation; b3. annotating the breakpoint as being located in the UTR region, exon M or intron N according to the position of the breakpoint in the transcript, wherein M is the exon number of the CDS region where the breakpoint is located, and N is the exon number of the CDS with the smaller exon number when the breakpoint is between two CDS regions; and b4. Obtaining the direction information of the transcript of the gene, the minimum exon number and the maximum exon number of the CDS; and c. Annotate the mutation type produced by the SV based on the SV type and the annotation of each breakpoint in step b. The annotation of mutation types includes annotation of SVs whose two breakpoints are not in the same gene, including: If both breakpoints are in the intergenic region, the mutation type is annotated as an intergenic variant; If at least one breakpoint was in a gene region, the mutation type was annotated as gene fusion; The annotation of mutation types includes annotation of SVs whose two breakpoints are located in the same gene and whose SV type is deletion, including: (1) If the two breakpoints are on the same exon, the mutation type is annotated as a splicing region variation or a protein-altering variation based on the distance between at least one breakpoint and the edge of the exon; the protein-altering variation is annotated as a reading frame deletion or a frameshift variation based on whether the number of missing nucleotides is a multiple of 3 nucleotides; (2) If the two breakpoints are on the same intron, the mutation type is annotated as a splice site variation, a splice region variation, or an intron variation according to the distance between at least one breakpoint and the upstream and downstream edges of the intron; the splice site variation is annotated as a splice acceptor variation or a splice donor variation according to the distance between at least one breakpoint and the upstream edge of the intron or the downstream edge of the intron; (3) If at least one breakpoint is in the UTR region and the other breakpoint is in the upstream and downstream intergenic regions of the same UTR, the mutation type is annotated as 5'UTR variation or 3'UTR variation according to the UTR region; (4) If one of the two breakpoints is on the 5'UTR and the other is on the 3'UTR, the mutation type is annotated as copy number deletion variation; (5) For other cases, the mutation type was annotated as large genomic rearrangement; The annotation of mutation types includes annotation of SVs with two breakpoints located in the same gene and repeated SV types, including: (1) If the two breakpoints are on the same exon or on adjacent exons, the mutation type is annotated as a protein-altering variation; depending on whether the number of repeated nucleotides is a multiple of 3 nucleotides, the protein-altering variation is annotated as a reading frame duplication or a frameshift variation; (2) If one of the two breakpoints is in an intron and the other is in the same intron or in an exon or UTR adjacent to the intron, the mutation type is annotated as an intronic variant; (3) If at least one breakpoint was in the UTR region and the other breakpoint was in the upstream and downstream intergenic regions of the same UTR region, the mutation type was annotated as 5'UTR variation or 3'UTR variation according to the UTR region; (4) If one of the two breakpoints is on the 5'UTR and the other is on the 3'UTR, the mutation type is annotated as copy number duplication variation; (5) For other cases, the mutation type was annotated as large genomic rearrangement; The annotation of mutation types includes annotation of SVs whose two breakpoints are located in the same gene and whose SV type is inversion, including: (1) If the two breakpoints are on the same exon, the mutation type is annotated as a splice region variant or a protein-altering variant based on the distance between at least one breakpoint and the edge of the exon; the protein-altering variant is annotated as a gain-of-terminus variant or a missense variant based on whether there is a terminator in the coding frame after the mutation; (2) If the two breakpoints are on the same intron, the mutation type is annotated as a splice site variation, a splice region variation, or an intron variation according to the distance between at least one breakpoint and the upstream or downstream edge of the exon; the splice site variation is annotated as a splice acceptor variation or a splice donor variation according to the distance between at least one breakpoint and the upstream edge of the exon or the downstream edge of the exon; (3) If at least one breakpoint was in the UTR region and the other breakpoint was in the upstream and downstream intergenic regions of the same UTR region, the mutation type was annotated as 5'UTR variation or 3'UTR variation according to the UTR region; (4) For other cases, the mutation type was annotated as large genomic rearrangement; The method includes annotating a new fusion sequence generated by an SV whose mutation type is gene fusion, including: a. Construct fusion sequences according to SV direction, including: For SVs with FR (forward-reverse) orientation, the sequences on the two genome replies are connected in order; for SVs with RF (reverse-forward) orientation, the sequences on the two gene replies are connected in reverse order; for SVs with FF (forward-forward) orientation, the sequences on the downstream gene replies are reverse-complemented and connected to the upstream gene in order; for SVs with RR (reverse-reverse) orientation, the sequences on the upstream gene replies are reverse-complemented and connected to the downstream gene in order; b. According to the composition relationship and transcript direction of the two genes in the fusion sequence, the fusion sequence structure is determined as follows: For genes with fusion sequences upstream of the genome, if the transcript chain direction is forward, the fusion sequence retains the exon sequence upstream of the exon where the breakpoint is located; if the transcript chain direction is reverse, the fusion sequence retains the exon sequence downstream of the exon where the breakpoint is located; For genes with fusion sequences downstream of the genome, if the transcript chain direction is forward, the fusion sequence retains the exon sequence downstream of the exon where the breakpoint is located; if the transcript chain direction is reverse, the fusion sequence retains the exon sequence upstream of the exon where the breakpoint is located; c. According to the composition relationship and transcription direction of the two genes in the fusion sequence, there are four situations: (1) ++, i.e., the 5' end of the fusion sequence is the left side of the fusion sequence, and the 3' end is the right side of the fusion sequence; (2) --, i.e., the 3' end of the fusion sequence is the left side of the fusion sequence, and the 5' end is the right side of the fusion sequence; (3) +-, i.e., the transcription directions of the two genes in the fusion sequence conflict, and the 5' end and 3' end cannot be determined; and (4) -+, i.e., the transcription directions of the two genes in the fusion sequence conflict, and the 5' end and 3' end cannot be determined. Among them, situations (1) and (2) are predicted to be able to form fusion proteins, and situations (3) and (4) are predicted to be unable to form fusion proteins.

2. The method of claim 1, wherein step a comprises reading analysis results of upstream SV detection software to extract SV information.

3. The method according to claim 2, wherein the SV variation detection software is selected from Delly, Manta, LUMPY, GRIDSS and markSV.

4. The method of claim 1, further comprising annotating the inferred impact of the mutation type; wherein the copy number duplication variation is annotated as having a medium inferred impact, the copy number deletion variation and the large genomic rearrangement are annotated as having a high inferred impact, and the inferred impact of the remaining mutation types is annotated according to the prediction criteria of the snpEff software.

5. The method of claim 1, further comprising annotating the functional impact of the SV, comprising: For intergenic SVs, if one of the two genes is a proto-oncogene and predicted to produce a fusion protein, the functional impact was annotated as loss of function: if one gene is a tumor suppressor gene, the functional impact was annotated as loss of function; For SVs within genes, if the gene is a tumor suppressor gene, the functional impact is annotated as loss of function; for SVs within genes known to cause gain of function, the functional impact is annotated as gain of function; For other cases the functional impact is annotated as None.

6. The method of claim 1, wherein the method further comprises annotating the pathogenicity level of the SV, wherein: For intergenic SVs, if they were annotated as loss-of-function, the evidence annotation was PVS1 and the rating annotation was potentially pathogenic; For SVs within genes, the level of evidence and rating were annotated with reference to the ACMG CNV guidelines.

7. The method of claim 1, wherein the method further comprises generating annotations of the HGVS data structure based on the type and breakpoint coordinates of the SV.

8. The method of claim 7, wherein the HGVS data structure comprises a g.HGVS, c.HGVS and / or p.HGVS data structure.

9. A system for annotating genomic structural variations (SVs), the system comprising: SV information reading module; and SV Annotation Module, The SV information reading module is configured to read SV information, the information including the type, direction and two breakpoint coordinates of the SV, wherein the SV type includes deletion, duplication, inversion and translocation; and the SV annotation module is configured to annotate the SV by the method described in any one of claims 1-8.

10. An apparatus for annotating genomic structural variations (SVs), comprising: memory for storing computer program instructions; and a processor for executing computer program instructions, Wherein when the computer program instructions are executed by the processor, the device performs the method of any one of claims 1-8.

11. A computer-readable medium storing computer program instructions, wherein when the computer program instructions are executed by a processor, the method of any one of claims 1 to 8 is implemented.

Images