SMN1 and SMN2 gene haplotype assembly system and method based on PacBio sequencing

Abstract

The present invention relates to a PacBio sequencing-based SMN1 and SMN2 gene haplotype assembly system and method. The present invention extracts sequences of 1 kb regions upstream and downstream of SMN1 and SMN2 by performing interval extraction on a generated bam format file, and converts the sequences into FASTQ format files. The bases of chromosome 5 other than the 3 kb region upstream and downstream of SMN1 are replaced with "N", and an alignment index is constructed for the reference sequence to form an SMN1 alignment reference file. The FASTQ format file is aligned to the SMN1 alignment reference file to obtain a re-aligned bam format file. Key SNP sites obtained from the re-aligned bam format file are preliminarily grouped, and haplotype assembly is performed to obtain haplotype sequences of SMN1 and SMN2. The present invention effectively identifies the haplotype where the pathogenic site is located, and provides more reliable data support for the typing, diagnosis and clinical research of SMA.

Description

Technical Field

[0001] The present invention relates to the field of bioinformatics analysis and genomics technology, and in particular to a system and method for assembling haplotypes of SMN1 and SMN2 genes based on PacBio sequencing. Background Art

[0002] Spinal muscular atrophy (SMA) is a common autosomal recessive genetic disease. The primary causative gene is the SMN1 (survival motor neuron 1) gene, which is highly similar to the SMN2 gene and exists in multiple copies. Due to the high sequence homology between the SMN1 and SMN2 genes and the complex nature of the differing bases or variants, sequencing analysis and mutation detection are challenging.

[0003] With the rapid development of third-generation sequencing technology, high-fidelity (HiFi) long-read data obtained based on the PacBio platform can directly sequence large fragments of complex genomic regions with low error rates, providing a powerful tool for typing the SMN1 / SMN2 gene region. However, due to the limitations of repeated regions in the reference genome and interference from only a few differential sites between genes, accurate typing of SMN1 and SMN2 and assembly of complete haplotypes still require the development of high-precision bioinformatics methods and processes.

[0004] Currently, only the Paraphase software can be used to assemble haplotypes for homologous genes. However, for highly repetitive regions such as SMN1 / SMN2, problems such as insufficient assembly range and missing assembly of some exons still occur. Therefore, a new method and implementation tool for high-precision variant identification and haplotype assembly in SMA-related gene regions based on PacBio HiFi long-read sequencing data is urgently needed to enable more comprehensive and accurate analysis of SMN1 and SMN2. Summary of the Invention

[0005] In response to the above-mentioned deficiencies in the prior art, the present invention aims to provide a PacBio sequencing-based haplotype assembly system and method for SMN1 and SMN2 genes, so as to accurately identify complex copy number variations, insertions and deletions, and single nucleotide variations of the genes.

[0006] In order to solve the above problems, the present invention adopts the following technical solutions:

[0007] On the one hand, the present invention provides a PacBio sequencing-based SMN1 and SMN2 gene haplotype assembly system, comprising an alignment module, a sequence extraction module, a reference genome construction module, a re-alignment module, and an assembly module;

[0008] The alignment module is used to align the PacBio HiFi original sequencing sequence to the human genome reference sequence using minimap2 software based on the map-hifi mode to generate a bam format file;

[0009] The sequence extraction module is used to perform interval extraction on the generated bam format file using samtools, extract the sequences of the 1 kb regions upstream and downstream of SMN1 and SMN2, and convert them into FASTQ format files;

[0010] The reference genome construction module is used to extract the human reference genome chromosome chr5 sequence, retain the reference sequence 3kb upstream and downstream of SMN1, replace the other bases of chromosome 5 except the 3kb region upstream and downstream of SMN1 with "N", and use minimap2 software to construct an alignment index for the reference sequence to form an SMN1 alignment reference file;

[0011] The re-alignment module is used to align the FASTQ format file to the SMN1 alignment reference file to obtain a re-aligned bam format file;

[0012] The assembly module is used to preliminarily group SMN1 and SMN2 according to the key SNP sites obtained by re-aligning the bam format file, and perform haplotype assembly to obtain the haplotype sequences of SMN1 and SMN2.

[0013] As an implementation method, the 1 kb upstream and downstream regions of SMN1 and SMN2 are chr5:70924086-70954015 and chr5:70048668-70078595, respectively.

[0014] As an embodiment, the region 3kb upstream and downstream of SMN1 is chr5:70922086-70956015.

[0015] As an implementable embodiment, the assembly module includes a variation identification and sequence grouping unit, a filtering unit, a haplotype assembly unit and an output unit;

[0016] The variant identification and sequence grouping unit is used to traverse each sequence, filter low-quality alignments and secondary alignments of the re-aligned bam format files, identify SNVs, insertions, deletions and complex variants, record site depth and site variation frequency, and preliminarily group reads according to the characteristic differential variations of SMN1 and SMN2;

[0017] The filtering unit is used to record the site depth and site variation frequency, and filter the noise in the error variation-frequent region according to the interference site set, and filter the low-frequency noise variation;

[0018] The haplotype assembly unit is used to classify unassigned reads into corresponding haplotype branches or create new haplotypes based on variation frequency, reads coverage and iterative assembly strategy;

[0019] The output unit is used to generate a SAM output file, add an HP tag, and annotate each read as SMN1 / 2_hap_ID or Unassigned; convert the annotated SAM output file to bam and create an index, and view the typing results in a visualization tool.

[0020] As an implementation method, the method of classifying unassigned reads into corresponding haplotype branches or creating new haplotypes based on variation frequency, read coverage, and iterative assembly strategy includes:

[0021] Using the preliminary grouping of SMN1 and SMN2 as the initial cluster, the sequences were clustered and expanded based on all variant sites, and small or noisy groups were eliminated using variant frequency and read coverage.

[0022] For sequences that have not yet been assigned to any haplotype, an iterative assembly and variation alignment strategy is used to assign them to an existing haplotype or a new haplotype branch;

[0023] After the assembly is completed, the variation frequency and sequence consistency within each haplotype are filtered and alleles are completed to obtain the final haplotype set.

[0024] In another aspect, the present invention provides a method for assembling haplotypes of SMN1 and SMN2 genes based on PacBio sequencing, comprising:

[0025] Use minimap2 software based on the map-hifi mode to align the PacBio HiFi raw sequencing sequence to the human genome reference sequence and generate a bam format file;

[0026] Samtools was used to extract the intervals of the generated bam format files, extract the sequences of the 1 kb upstream and downstream regions of SMN1 and SMN2, and convert them into FASTQ format files;

[0027] Extract the human reference genome chromosome chr5 sequence, retain the reference sequence of 3 kb upstream and downstream of SMN1, replace all other bases of chromosome 5 except the 3 kb region upstream and downstream of SMN1 with "N", and use minimap2 software to construct an alignment index for the reference sequence to form an SMN1 alignment reference file;

[0028] Align the FASTQ format file to the SMN1 alignment reference file to obtain a re-aligned bam format file;

[0029] SMN1 and SMN2 were preliminarily grouped according to the key SNP sites obtained from the re-alignment bam format file, and haplotype assembly was performed to obtain the haplotype sequences of SMN1 and SMN2.

[0030] As an implementation method, the 1 kb upstream and downstream regions of SMN1 and SMN2 are chr5:70924086-70954015 and chr5:70048668-70078595, respectively.

[0031] As an embodiment, the region 3kb upstream and downstream of SMN1 is chr5:70922086-70956015.

[0032] As an implementation method, the key SNP sites obtained by re-aligning the bam format file are preliminarily grouped for SMN1 and SMN2, and haplotype assembly is performed to obtain the haplotype sequences of SMN1 and SMN2, including:

[0033] Traverse each sequence, filter low-quality and secondary alignments in the re-aligned bam format file, identify SNVs, insertions, deletions, and complex variants, and record site depth and site variation frequency; preliminarily group reads based on the characteristic differential variation of SMN1 and SMN2;

[0034] It is used to record site depth, site variation frequency, and filter noise in areas with frequent erroneous variations based on the interference site set, as well as to filter low-frequency noise variations.

[0035] Based on the variation frequency, read coverage and iterative assembly strategy, unassigned reads are classified into corresponding haplotype branches or newly created haplotypes;

[0036] Generate a SAM output file and add the HP tag, annotating each read as SMN1 / 2_hap_ID or Unassigned; convert the annotated SAM output file to bam and create an index, and view the typing results in a visualization tool.

[0037] As an implementation method, the method of classifying unassigned reads into corresponding haplotype branches or creating new haplotypes based on variation frequency, read coverage, and iterative assembly strategy includes:

[0038] Using the preliminary grouping of SMN1 and SMN2 as the initial cluster, the sequences were clustered and expanded based on all variant sites, and small or noisy groups were eliminated using variant frequency and read coverage.

[0039] For sequences that have not yet been assigned to any haplotype, an iterative assembly and variation alignment strategy is used to assign them to an existing haplotype or a new haplotype branch;

[0040] After the assembly is completed, the variation frequency and sequence consistency within each haplotype are filtered and alleles are completed to obtain the final haplotype set.

[0041] The present invention achieves more comprehensive SMN1 / SMN2 haplotype coverage by utilizing interval extraction from the reference genome, constructing local micro-range reference sequences, performing preliminary grouping of reads at key SNP sites, and combining iterative haplotype assembly. Compared to the existing software Paraphase, this method achieves full-genomic coverage of SMN1 / SMN2 haplotypes while maintaining a consistent number of haplotypes. This allows for the precise localization of pathogenic variants in the entire SMN1 and SMN2 genomes, as well as the identification of cis-trans relationships among multiple pathogenic variants. This significantly enhances the typing, diagnosis, and treatment of SMA. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] Figure 1 Schematic diagram of a PacBio sequencing-based haplotype assembly system for SMN1 and SMN2 genes in Example 1 of the present invention.

[0043] Figure 2 This is a flow chart of a method for assembling SMN1 and SMN2 gene haplotypes based on PacBio sequencing in Example 1 of the present invention.

[0044] Figure 3 This is the Paraphase haplotype assembly diagram of the T1 sample in Example 2 of the present invention.

[0045] Figure 4 This is the haplotype assembly diagram of the present invention for the T1 sample in Example 2 of the present invention.

[0046] Figure 5 This is the Paraphase haplotype assembly diagram of the T2 sample in Example 2 of the present invention.

[0047] Figure 6 This is the haplotype assembly diagram of the present invention for the T2 sample in Example 2 of the present invention.

[0048] Figure 7 This is the Paraphase haplotype assembly diagram of the T3 sample in Example 2 of the present invention.

[0049] Figure 8 This is the haplotype assembly diagram of the T3 sample in Example 2 of the present invention. DETAILED DESCRIPTION

[0050] The present invention will be further described in detail below with reference to specific embodiments.

[0051] It should be noted that these embodiments are only used to illustrate the present invention, rather than to limit the present invention. Simple improvements to the method based on the concept of the present invention fall within the scope of protection claimed by the present invention. Example

[0052] See also Figure 1 , which is a PacBio sequencing-based SMN1 and SMN2 gene haplotype assembly system, including an alignment module 100, a sequence extraction module 200, a reference genome construction module 300, a re-alignment module 400 and an assembly module 500.

[0053] The alignment module 100 is used to align the PacBio HiFi original sequencing sequence to the human genome reference sequence using the minimap2 software based on the map-hifi mode to generate a bam format file.

[0054] This embodiment provides the specific processing process of the above modules:

[0055] Using the long sequence alignment software minimap2 in map-hifi mode, we aligned the PacBio high-fidelity HIFI sequences to the human hg38 reference genome, generating a binary alignment bam file. Using the bam processing software samtools, we sorted the sequences in the bam file by position, generating the sorted.bam file.

[0056] minimap2 -t 8 -ax map-hifi hg38.mmi pacbio_hifi_reads.fq > hifi_map2hg38.sam

[0057] samtools view -Sb hifi_map2hg38.sam > hifi_map2hg38.bam

[0058] samtools sort hifi_map2hg38.bam > sorted.bam

[0059] samtools index sorted.bam.

[0060] The sequence extraction module 200 is used to perform interval extraction on the generated bam format file using samtools, extract the sequences of the 1 kb regions upstream and downstream of SMN1 and SMN2, and convert them into FASTQ format files.

[0061] This embodiment provides the specific processing process of the above modules:

[0062] Using the bam processing software samtools, we extracted the generated sorted.bam files and aligned them to the 1 kb upstream and downstream regions of SMN1 and SMN2 (SMN1 chr5:70924086-70954015, SMN2 chr5:70048668-70078595), generating the files SMN1-SMN2UD1KB.sorted.bam. We then used bamToFastq software to convert SMN1-SMN2UD1KB.sorted.bam to the FASTQ format files SMN1-SMN2UD1KB.sorted.bam.

[0063] samtools view -b sorted.bam chr5:70924086-70954015 > SMN1_UD1KB.bam

[0064] samtools view -b sorted.bam chr5:70048668-70078595 > SMN2_UD1KB.bam

[0065] samtools merge SMN1-SMN2UD1KB.sorted.bam SMN1_UD1KB.bam SMN2_UD1KB.bam

[0066] samtools index SMN1-SMN2UD1KB.sorted.bam

[0067] bamToFastq -i SMN1-SMN2UD1KB.sorted.bam -fq SMN1-SMN2UD1KB.sorted.bamtofastq.fastq.

[0068] The reference genome construction module 300 is used to extract the human reference genome chromosome chr5 sequence, retain the reference sequence of 3kb upstream and downstream of SMN1, replace the other bases of chromosome 5 except the 3kb region upstream and downstream of SMN1 with "N", and use minimap2 software to construct an alignment index for the reference sequence to form an SMN1 alignment reference file.

[0069] This embodiment provides the specific processing process of the above modules:

[0070] The human hg38 reference genome chr5 sequence was extracted, retaining the 3 kb region upstream and downstream of SMN1 (SMN1chr5:70922086-70956015) and replacing the bases in other regions (A / T / C / G / a / t / c / g) with N. The long sequence alignment software minimap2 was used to construct the alignment reference file hg38_chr5_SMN1_UD3KB.mmi.

[0071] minimap2 -d hg38_chr5_SMN1_UD3KB.mmi hg38_chr5_SMN1_UD3KB.fa.

[0072] The realignment module 400 is used to align the FASTQ format file to the SMN1 alignment reference file to obtain a realigned bam format file.

[0073] This embodiment provides the specific processing process of the above modules:

[0074] Using the long sequence alignment software minimap2, the generated SMN1-SMN2UD1KB.sorted.bamtofastq.fastq was aligned to the generated hg38_chr5_SMN1_UD3KB.mmi file. This generated the SMN1 3kb upstream and downstream region realignment file hg38SMN1UD3KB.SMN1-SMN2UD1KB.sorted.realigned.bam.

[0075] minimap2 -t 8 -ax map-hifi hg38_chr5_SMN1_UD3KB.mmi SMN1-SMN2UD1KB.sorted.bamtofastq.fastq > hg38SMN1UD3KB.SMN1-SMN2UD1KB.sorted.realigned.sam

[0076] samtools view -Sb hg38SMN1UD3KB.SMN1-SMN2UD1KB.sorted.realigned.sam >hg38SMN1UD3KB.SMN1-SMN2UD1KB.sorted.realigned.bam

[0077] samtools index hg38SMN1UD3KB.SMN1-SMN2UD1KB.sorted.realigned.bam.

[0078] The assembly module 500 is used to preliminarily group SMN1 and SMN2 according to the key SNP sites obtained from the re-alignment bam format file, and perform haplotype assembly to obtain the haplotype sequences of SMN1 and SMN2.

[0079] The assembly module 500 includes a variation identification and sequence grouping unit 510 , a filtering unit 520 , a haplotype assembly unit 530 and an output unit 540 .

[0080] The variation identification and sequence grouping unit 510 is used to traverse each sequence, filter out low-quality alignments and secondary alignments in the re-aligned bam format file, identify SNVs, insertions, deletions, and complex variations; and preliminarily group the reads (SMN1 group, SMN2 group, others group) based on the characteristic differential variation between SMN1 and SMN2 (chr5:70951946 C>T).

[0081] The filtering unit 520 is used to record the site depth and site variation frequency, and filter the noise in the error variation-frequent region according to the interference site set, and filter the low-frequency noise variation.

[0082] The haplotype assembly unit 530 is used to classify unassigned reads into corresponding haplotype branches or create new haplotypes based on variation frequency, read coverage, and iterative assembly strategy.

[0083] Among them, based on the variation frequency, reads coverage and iterative assembly strategy, unassigned reads are classified into corresponding haplotype branches or newly created haplotypes, including:

[0084] Using the preliminary grouping of SMN1 and SMN2 as the initial cluster, the sequences were clustered and expanded based on all variant sites, and small or noisy groups were eliminated using variant frequency and read coverage.

[0085] For sequences that have not yet been assigned to any haplotype, an iterative assembly and variation alignment strategy is used to assign them to an existing haplotype or a new haplotype branch;

[0086] After the assembly is completed, the variation frequency and sequence consistency within each haplotype are filtered and alleles are completed to obtain the final haplotype set.

[0087] The output unit 540 is used to generate a SAM output file and add an HP tag to annotate each read as SMN1 / 2_hap_ID or Unassigned; the annotated SAM output file is converted to bam and indexed, and the typing results are viewed in a visualization tool.

[0088] This embodiment provides the specific processing process of the above modules:

[0089] (1) Variant identification and sequence grouping

[0090] (a) Initialization / Reading

[0091] - refgenome = "path / to / reference-genome.fa"

[0092] -defPos=70951946

[0093] # Determine the key sites of SMN1 / SMN2

[0094] - errorPos = {70920535,70920820,70920829,70920886..70920899,70921060..70921085,70921810,70922100,70922929,70923056,70923087,70923097,70923099,70923210..70923228,70923424,70923747,70924129,70924594,70924610,70924859..70924892,70926342..70926357,70926274,70926534,70924686,70927608..70927628,70927961,70928141,70928300..70928325,70929146..70929167,70929398..70929427,70929579,70929880,70930594,70931216..70931243,70932550,70932884..70932897,70933345,70933630,70934232,70934675,70934978,70935162,70935773,70935834,70935982,70936151,70936599,70936612..70936646,70936804,70938110,70938227,70934970,70940003,70941839,70942568,70942834,70942983,70943024,70943030,70943289,70943467,70943791,70944690,70945126,70945571,70945633,70945931,70946513,70946926..70946941,70947107,70947528,70947686,70947980,70948155,70948975,70948993,70949682,70949903,70949924,70950068,70950417,70950439,70950501,70950644,70951281,70951909,70952289,70952951,70953000,70954485,70954725,70954855..70954907,70955214,70955275,70955625,70956174}.

[0095] # Interference site set / hash

[0096] - Input: bamFile (PacBio sequencing alignment result file)

[0097] - Open bamFile and iterate over the aligned records using Bio::DB::Sam or other methods

[0098] # Define data structure

[0099] - seqvar[queryID][pos] # Save the variation information of each read at each variation site

[0100] - var[pos][variant] # Record multiple variants that may appear at a certain site and their corresponding reads

[0101] - depth[pos] # site coverage

[0102] - seqinfo[queryID] # Record the alignment interval (start, end) of each read

[0103] - seqGroup = { 'SMN1': [], 'SMN2': [], 'others': []} # Preliminary grouping of reads ID

[0104] - ALLELES[pos] # All alleles and counts at this site.

[0105] (b) Traverse the BAM alignment records to identify variants and sequence groups (for each alignment):

[0106] - If it is a secondary alignment / supplemental alignment (flag & 256 or flag & 2048), skip it

[0107] - Get:

[0108] queryID, alignment start and end positions (start, end)

[0109] query sequence (queryseq), reference sequence (ref), mismatch information

[0110] - Statistical coverage: depth[pos]++ for each site in the start..end interval

[0111] - Based on mismatch, distinguish SNV / insertion / deletion / complex variation:

[0112] var[pos][ "ref>alt" ].seq ← push( queryID )

[0113] seqvar[queryID][pos].base = "altBases" ...

[0114] - Determine whether defPos falls within [start, end]:

[0115] if (defPos covers this reads) then

[0116] if (the reads base at defPos == 'T')

[0117] seqGroup['SMN2'].push(queryID)

[0118] else

[0119] seqGroup['SMN1'].push(queryID)

[0120] else

[0121] seqGroup['others'].push(queryID).

[0122] (2) Variation filtering (varfilter)

[0123] - For each variation in var[pos]:

[0124] Calculate the number of mutations / depth[pos] => rate

[0125] If the rate is too low or too high or in the interference position (errorPos) => delete

[0126] Arrange the remaining variation to varfilter.

[0127] (3) Haplotype assembly

[0128] (a) Read the reference genome into memory:

[0129] - refGenomeBase[pos] = referenceBase.

[0130] (b) Generate ALLELES

[0131] - Count the candidate alleles for each site based on the remaining variants in varfilter + reference alleles.

[0132] (c) Establish seqAvailablePos[queryID]

[0133] - List the alleles covered by the reads (including insertions, deletions, SNVs, etc.)

[0134] - Recorded as: seqAvailablePos[queryID][pos] = base.

[0135] (d) Initialize the haplotype structure hap = {}

[0136] - RA_reads is used to mark reads for subsequent reallocation

[0137] - Define a series of functions:

[0138] Function PreAssemble(gene):

[0139] - Perform preliminary assembly of reads for a specified gene (such as SMN1 / SMN2)

[0140] - Clustering to form hap[gene][n].allele

[0141] - Eliminate too few or noisy clusters

[0142] - Remove assembled reads from seqGroup

[0143] Function AlleleImputation():

[0144] - Calculate the coverage interval of each haplotype at start..end

[0145] - Perform allele imputation or consistency correction (noise filtering)

[0146] Function Assemble(reads, constraint):

[0147] - Try to match the specified set of reads to an existing hap

[0148] - If the match is unique, assign it to the corresponding hap

[0149] - If there are multiple matches or no matches, create a new hap or mark it as unassigned

[0150] - Stored in hap.unassignedReadsGroup[...] for subsequent redistribution

[0151] Function Assemble_RA(matches, group):

[0152] - Perform secondary redistribution of unassigned reads based on multiple existing hap

[0153] - If it still cannot be allocated, it will remain in the unassignedReadsGroup

[0154] Function GetRefAllele(pos):

[0155] - If the varfilter exists for this site, the most reasonable reference allele is selected

[0156] - Otherwise get from refGenomeBase[pos]

[0157] Function arrays_cmp(array1, array2):

[0158] - Determine whether the contents of two arrays are the same

[0159] Function arrays_intersection(array1, array2):

[0160] - Returns the intersection of the two

[0161] Function lowNumsThreshold(array):

[0162] - Dynamically calculate the threshold for removing low-frequency groups (such as max * 0.2 or fixed threshold, etc.)

[0163] (f) Call in sequence:

[0164] - PreAssemble('SMN1')

[0165] - PreAssemble('SMN2')

[0166] - AlleleImputation()

[0167] - Get the unassembled reads of seqGroup['SMN1'], seqGroup['SMN2'], seqGroup['others'] => unassembleSeqs

[0168] - Loop Assemble(unassembleSeqs, constraint=1), then AlleleImputation()

[0169] * If there is no new allocation or the loop ends => exit

[0170] - Repeat the above steps (constraint=0 no longer limits the overlapping region) to assemble the remaining reads

[0171] - After the end, the final haplotype information is stored in hap

[0172] (4) Output:

[0173] - According to each haplotype in hap => add label HP:Z to each read:<hap_id> If it is the main haplotype, you can add color label YC:Z:<R,G,B>

[0174] - Output to new .sam file => Convert to .bam => Index

[0175] - Delete temporary sam files.

[0176] See also Figure 2 , a PacBio sequencing-based method for assembling haplotypes of the SMN1 and SMN2 genes (not intended for disease diagnosis), including:

[0177] S100. Use minimap2 software based on the map-hifi mode to align the PacBio HiFi raw sequencing sequence to the human genome reference sequence and generate a bam format file.

[0178] S200. Use samtools to perform interval extraction on the generated bam format file, extract the sequences of the 1 kb upstream and downstream regions of SMN1 and SMN2, and convert them into FASTQ format files.

[0179] S300. Extract the human reference genome chromosome chr5 sequence, retain the reference sequence of 3 kb upstream and downstream of SMN1, replace the other bases of chromosome 5 except the 3 kb region upstream and downstream of SMN1 with "N", and use minimap2 software to build an alignment index for the reference sequence to form an SMN1 alignment reference file.

[0180] S400 , aligning the FASTQ format file to the SMN1 alignment reference file to obtain a re-aligned bam format file.

[0181] S500: Preliminary grouping of SMN1 and SMN2 is performed based on the key SNP sites obtained from the re-alignment bam format file, and haplotype assembly is performed to obtain the haplotype sequences of SMN1 and SMN2.

[0182] Among them, the 1kb upstream and downstream regions of SMN1 and SMN2 are chr5:70924086-70954015 and chr5:70048668-70078595 respectively.

[0183] Among them, the region 3kb upstream and downstream of SMN1 is chr5:70922086-70956015.

[0184] Among them, SMN1 and SMN2 were preliminarily grouped according to the key SNP sites obtained from the re-alignment bam format file, and haplotype assembly was performed. The haplotype sequences of SMN1 and SMN2 obtained include:

[0185] Traverse each sequence, filter low-quality and secondary alignments in the re-aligned bam format file, identify SNVs, insertions, deletions, and complex variants, and record site depth and site variation frequency; preliminarily group reads based on the characteristic differential variation of SMN1 and SMN2;

[0186] Record site depth and site variation frequency, and filter out noise in areas with frequent erroneous variations based on the interference site set, as well as low-frequency noise variations;

[0187] Based on the variation frequency, read coverage and iterative assembly strategy, unassigned reads are classified into corresponding haplotype branches or newly created haplotypes;

[0188] Generate a SAM output file and add the HP tag, annotating each read as SMN1 / 2_hap_ID or Unassigned; convert the annotated SAM output file to bam and create an index, and view the typing results in a visualization tool.

[0189] Among them, based on the variation frequency, reads coverage and iterative assembly strategy, unassigned reads are classified into corresponding haplotype branches or newly created haplotypes, including:

[0190] Using the preliminary grouping of SMN1 and SMN2 as the initial cluster, the sequences were clustered and expanded based on all variant sites, and small or noisy groups were eliminated using variant frequency and read coverage.

[0191] For sequences that have not yet been assigned to any haplotype, an iterative assembly and variation alignment strategy is used to assign them to an existing haplotype or a new haplotype branch;

[0192] After the assembly is completed, the variation frequency and sequence consistency within each haplotype are filtered and alleles are completed to obtain the final haplotype set.

[0193] Example 2

[0194] Comparison of SMN1 / SMN2 haplotypes of three samples using this method and Paraphase assembly

[0195] Use the published Paraphase method to assemble SMN1 / SMN2 haplotypes:

[0196] Step 1: Same as step S100 in embodiment 1.

[0197] Step 2: Use the homologous gene haplotype assembly software Paraphase to assemble haplotypes of the aligned sequences of the SMA-related genes SMN1 and SMN2 regions based on sorted.bam.

[0198] Paraphrase -b sorted.bam -r genome.fa --samtools samtools --minimap2minimap2 -o . / -g smn1 --genome 38.

[0199] 3 cases of SMN1 / SMN2 long fragment capture Pacbio sequencing generated HiFi sequence data, using Example 1 and Paraphase method process to assemble haplotypes of SMN1 / SMN2, and the results were compared, see Table 1 and Figure 3-8 The basic statistics of Pacbio sequencing data of the three test samples are shown in Table 2.

[0200] Table 1 Comparison of SMN1 / SMN2 haplotypes assembled by Paraphase

[0201]

[0202] Table 2 Basic statistics of Pacbio sequencing data for three test samples

[0203]

[0204] Sample: sample ID;

[0205] TotalReads: total number of sequences;

[0206] TotalBases: total number of bases;

[0207] MidLength: median sequence length;

[0208] > Q20 (%): the percentage of base quality greater than Q20;

[0209] > Q30 (%): the percentage of bases with a quality greater than Q30;

[0210] MappedRatio (%): ratio of human genome sequences that can be aligned;

[0211] AveDepth (X): average depth of the target area;

[0212] Coverage (%): coverage of the target area.

[0213] As can be seen from the table above, the present invention can significantly improve the assembly level and accuracy of the full-length region of the SMN1 / SMN2 gene (especially exon 1) while maintaining or exceeding the haplotype quantitative resolution of existing methods, providing more complete data support for subsequent diagnostic analysis and clinical research.

[0214] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not limiting. Although the present invention has been described with reference to the preferred embodiments of the present invention, it should be understood by those skilled in the art that various changes can be made in form and details without departing from the spirit and scope of the present invention as defined in the appended claims.

Claims

1. A PacBio sequencing-based SMN1 and SMN2 gene haplotype assembly system, characterized in that: It includes alignment module, sequence extraction module, reference genome construction module, re-alignment module and assembly module; The alignment module is used to align the PacBio HiFi original sequencing sequence to the hg38 version of the human genome reference sequence to generate a bam format file; The sequence extraction module is used to perform interval extraction on the generated bam format file using samtools, extract the sequences of the 1 kb regions upstream and downstream of SMN1 and SMN2, and convert them into FASTQ format files; The reference genome construction module is used to extract the hg38 version of the human reference genome chromosome chr5 sequence, retain the reference sequence 3kb upstream and downstream of SMN1, replace the other bases of chromosome chr5 except the 3kb region upstream and downstream of SMN1 with "N", and use minimap2 software to construct an alignment index for the reference sequence to form an SMN1 alignment reference file; The re-alignment module is used to align the FASTQ format file to the SMN1 alignment reference file to obtain a re-aligned bam format file; The assembly module is used to preliminarily group SMN1 and SMN2 according to the key SNP sites obtained by re-aligning the bam format file, and perform haplotype assembly to obtain the haplotype sequences of SMN1 and SMN2; The assembly module includes a variation identification and sequence grouping unit, a filtering unit, a haplotype assembly unit and an output unit; The variant identification and sequence grouping unit is used to traverse each sequence, filter low-quality and secondary alignments of the re-aligned bam format file, identify SNVs, insertions, deletions, and complex variants, record site depth and site variation frequency, and preliminarily group reads based on the characteristic differential variation chr5:70951946 C>T between SMN1 and SMN2; The filtering unit is used to record the site depth and site variation frequency, and filter the noise in the error variation-frequent region according to the interference site set, and filter the low-frequency noise variation; The haplotype assembly unit is used to classify unassigned reads into corresponding haplotype branches or create new haplotypes based on variation frequency, reads coverage and iterative assembly strategy; The output unit is used to generate a SAM output file and add an HP tag, annotating each read as SMN1 / 2_hap_ID or Unassigned; Convert the annotated SAM output file to bam and create an index, and view the typing results in a visualization tool; The classification of unassigned reads into corresponding haplotype branches or newly created haplotypes based on variation frequency, read coverage, and iterative assembly strategy includes: Using the preliminary grouping of SMN1 and SMN2 as the initial cluster, the sequences were clustered and expanded based on all variant sites, and small or noisy groups were eliminated using variant frequency and read coverage. For sequences that have not yet been assigned to any haplotype, an iterative assembly and variation alignment strategy is used to assign them to an existing haplotype or a new haplotype branch; After the assembly is completed, the variation frequency and sequence consistency within each haplotype are filtered and alleles are completed to obtain the final haplotype set.

2. A method for assembling SMN1 and SMN2 gene haplotypes for non-disease diagnosis purposes based on PacBio sequencing, characterized in that: include: Align the PacBio HiFi raw sequencing sequence to the hg38 version of the human genome reference sequence and generate a bam format file; Samtools was used to extract the intervals of the generated bam format files, extract the sequences of the 1 kb upstream and downstream regions of SMN1 and SMN2, and convert them into FASTQ format files; Extract the hg38 version of the human reference genome chromosome chr5 sequence, retaining the reference sequence 3kb upstream and downstream of SMN1, and replace all other bases in chromosome chr5 except the 3kb region upstream and downstream of SMN1 with "N". Use minimap2 software to construct an alignment index for the reference sequence to form an SMN1 alignment reference file; Align the FASTQ format file to the SMN1 alignment reference file to obtain a re-aligned bam format file; SMN1 and SMN2 were preliminarily grouped according to the key SNP sites obtained from the re-alignment bam format file, and haplotype assembly was performed to obtain the haplotype sequences of SMN1 and SMN2; The key SNP sites obtained by re-aligning the bam format file are used to preliminarily group SMN1 and SMN2, and perform haplotype assembly to obtain the haplotype sequences of SMN1 and SMN2, including: Each sequence was traversed, and low-quality and secondary alignments in the re-aligned bam format files were filtered out to identify SNVs, insertions, deletions, and complex variants, and the site depth and site variation frequency were recorded. Reads were preliminarily grouped based on the characteristic differential variation chr5:70951946 C>T between SMN1 and SMN2. Record site depth and site variation frequency, and filter out noise in areas with frequent erroneous variations based on the interference site set, as well as low-frequency noise variations; Based on the variation frequency, read coverage and iterative assembly strategy, unassigned reads are classified into corresponding haplotype branches or newly created haplotypes; Generate a SAM output file and add the HP tag, annotating each read as SMN1 / 2_hap_ID or Unassigned; convert the annotated SAM output file to bam and create an index, and view the typing results in a visualization tool; The classification of unassigned reads into corresponding haplotype branches or newly created haplotypes based on variation frequency, read coverage, and iterative assembly strategy includes: Using the preliminary grouping of SMN1 and SMN2 as the initial cluster, the sequences were clustered and expanded based on all variant sites, and small or noisy groups were eliminated using variant frequency and read coverage. For sequences that have not yet been assigned to any haplotype, an iterative assembly and variation alignment strategy is used to assign them to an existing haplotype or a new haplotype branch; After the assembly is completed, the variation frequency and sequence consistency within each haplotype are filtered and alleles are completed to obtain the final haplotype set.

Images