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

Abstract

The invention relates to an SMN1 and SMN2 gene haplotype assembly system and an SMN1 and SMN2 gene haplotype assembly method based on PacBio sequencing. The method comprises the following steps: performing interval extraction on a generated bam-format file, extracting sequences of upstream and downstream 1kb regions of SMN1 and SMN2, converting the sequences into FASTQ-format files, replacing other bases of chromosomes 5 except the upstream and downstream 3kb regions of SMN1 with 'N', constructing a comparison index for a reference sequence, forming an SMN1 comparison reference file, comparing the FASTQ-format file to the SMN1 comparison reference file, and comparing the FASTQ-format file to the SMN1 comparison reference file. The method comprises the following steps of: obtaining a re-comparison bam-format file, carrying out preliminary grouping according to key SNP sites obtained by the re-comparison bam-format file, and carrying out haplotype assembly to obtain haplotype sequences of SMN1 and SMN2; according to the method, the haplotype of the pathogenic site is effectively identified, and more reliable data support is provided for SMA typing, diagnosis and clinical research.

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 main pathogenic gene is the SMN1 (survival motor neuron 1) gene, which is highly similar to the SMN2 gene and has a multi-copy structure. Since the SMN1 and SMN2 gene sequences are highly homologous and the different bases or variant types are complex, sequencing analysis and mutation detection are very difficult.

[0003] With the rapid development of third-generation sequencing technology, high-fidelity (HiFi) long-read data obtained based on the PacBio platform can directly measure large fragment sequences in complex genomic regions with a low error rate, providing a powerful means for typing detection of SMN1 / SMN2 gene regions. However, due to the interference of repeated regions in the reference genome and only a few differential sites between genes, it is still necessary to design high-precision bioinformatics methods and processes to accurately type SMN1 and SMN2 and assemble complete haplotypes.

[0004] Currently, only the software Paraphase can be used for haplotype assembly of homologous genes, but for highly repetitive regions such as SMN1 / SMN2, there are still problems such as insufficient assembly range and missing assembly of some exons. Therefore, a new method and implementation tool for high-precision variant identification and haplotype assembly for SMA-related gene regions based on PacBio HiFi long-read sequencing data is urgently needed to analyze SMN1 and SMN2 more comprehensively and accurately. Summary of the invention

[0005] In view of the above-mentioned deficiencies in the prior art, the present invention aims to provide a PacBio sequencing-based SMN1 and SMN2 gene haplotype assembly system and method to achieve accurate identification of complex copy number variations, insertions and deletions, and single nucleotide variations of the gene.

[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 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 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:70048668-70078595 and chr5:70924086-70954015, respectively.

[0014] As an implementation method, the 3kb upstream and downstream region of SMN1 is chr5:70046668-70080595.

[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 variation identification and sequence grouping unit is used to traverse each sequence, filter low-quality alignments and secondary alignments of re-aligned bam format files, identify SNVs, insertions, deletions and complex variations, and record site depth and site variation frequency, and preliminarily group reads according to the characteristic difference variation of SMN1 and SMN2;

[0017] The filtering unit is used to record the site depth, 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 implementable method, the method of classifying unassigned reads into corresponding haplotype branches or creating new haplotypes based on variation frequency, reads coverage and iterative assembly strategy includes:

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

[0022] For sequences that have not yet been assigned to any haplotype, they are assigned to an existing haplotype or a newly created haplotype branch through iterative assembly and variation alignment strategies;

[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 original sequencing sequence to the human genome reference sequence and generate a bam format file;

[0026] Use samtools to extract intervals from the generated bam format files, extract the sequences of the 1 kb regions upstream and downstream 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 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 build 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:70048668-70078595 and chr5:70924086-70954015, respectively.

[0031] As an implementation method, the 3kb upstream and downstream region of SMN1 is chr5:70046668-70080595.

[0032] As an implementable 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 out 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 variations 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 filter low-frequency noise variations.

[0035] Based on the variation frequency, reads 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 the visualization tool.

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

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

[0039] For sequences that have not yet been assigned to any haplotype, they are assigned to an existing haplotype or a newly created haplotype branch through iterative assembly and variation alignment strategies;

[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 beneficial effects of the present invention are as follows: the present invention achieves more comprehensive SMN1 / SMN2 haplotype coverage by utilizing interval extraction of reference genome, constructing local micro-range reference sequences, performing preliminary grouping of reads at key SNP sites and combining iterative haplotype assembly. Compared with the existing software Paraphase, the present invention can achieve full gene region coverage of SMN1 / SMN2 haplotypes under the premise of consistent haplotype numbers, can achieve accurate positioning of pathogenic variants in the full gene region of SMN1 and SMN2, can identify the cis-trans relationship of multiple pathogenic variants, and can play a greater role in the typing, diagnosis and treatment research 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 haplotypes of SMN1 and SMN2 genes 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 present invention for the T3 sample in Example 2 of the present invention. DETAILED DESCRIPTION

[0050] The present invention is further described in detail below in conjunction with 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 all fall within the scope of protection claimed by the present invention.

[0052] Example 1

[0053] See also Figure 1 , which is a SMN1 and SMN2 gene haplotype assembly system based on PacBio sequencing, 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.

[0054] 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.

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

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

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

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

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

[0060] samtools index sorted.bam.

[0061] 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.

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

[0063] The generated sorted.bam was extracted using the bam processing software samtools and aligned to the sequences of the upstream and downstream 1kb regions of SMN1 and SMN2 (SMN1 chr5:70048668-70078595, SMN2 chr5:70924086-70954015) to generate the file SMN1-SMN2UD1KB.sorted.bam. The bamToFastq software was used to convert SMN1-SMN2UD1KB.sorted.bam into the sequence FASTQ format file SMN1-SMN2UD1KB.sorted.bamtofastq.fastq.

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

[0065] samtools view -b sorted.bam chr5:70924086-70954015 > SMN2_UD1KB.bam

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

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

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

[0069] 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.

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

[0071] Extract the human hg38 version of the reference genome chr5 chromosome sequence, retain the 3kb region upstream and downstream of SMN1 (SMN1chr5:70046668-70080595), and replace the bases in other regions (A / T / C / G / a / t / c / g) with N. Use the long sequence alignment software minimap2 to construct the alignment reference file hg38_chr5_SMN1_UD3KB.mmi.

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

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

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

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

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

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

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

[0079] The assembly module 500 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.

[0080] 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 .

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

[0082] 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.

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

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

[0085] The preliminary grouping of SMN1 and SMN2 was used as the initial clustering, and the sequences were clustered and expanded based on all variant sites. The small or noisy groups were eliminated using the variant frequency and read coverage.

[0086] For sequences that have not yet been assigned to any haplotype, they are assigned to an existing haplotype or a newly created haplotype branch through iterative assembly and variation alignment strategies;

[0087] 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.

[0088] The output unit 540 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 into bam and create an index, and view the typing results in a visualization tool.

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

[0090] (1) Variant identification and sequence grouping

[0091] (a) Initialization / reading

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

[0093] -defPos = 70951946

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

[0095] - 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}.

[0096] # Interference site set / hash

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

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

[0099] # Define data structure

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

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

[0102] - depth[pos] # site coverage

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

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

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

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

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

[0108] - Get:

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

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

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

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

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

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

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

[0117] if (defPos covers this reads) then

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

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

[0120] else

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

[0122] else

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

[0124] (2) Variation filtering (varfilter)

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

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

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

[0128] Arrange the remaining variation to varfilter.

[0129] (3) Haplotype assembly

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

[0131] - refGenomeBase[pos] = referenceBase.

[0132] (b) Generate ALLELES

[0133] - Based on the remaining variants in varfilter + reference alleles, the candidate alleles of each locus are counted.

[0134] (c) Establish seqAvailablePos[queryID]

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

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

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

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

[0139] - Define a series of functions:

[0140] Function PreAssemble(gene):

[0141] - Perform preliminary assembly of reads of specified genes (such as SMN1 / SMN2)

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

[0143] - Eliminate too few or noisy clusters

[0144] - Remove assembled reads from seqGroup

[0145] Function AlleleImputation():

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

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

[0148] Function Assemble(reads, constraint):

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

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

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

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

[0153] Function Assemble_RA(matches, group):

[0154] - Reassign unassigned reads based on multiple existing haps

[0155] - If still cannot be assigned, it will remain in unassignedReadsGroup

[0156] Function GetRefAllele(pos):

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

[0158] - Otherwise get from refGenomeBase[pos]

[0159] Function arrays_cmp(array1, array2):

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

[0161] Function arrays_intersection(array1, array2):

[0162] - Returns the intersection of the two

[0163] Function lowNumsThreshold(array):

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

[0165] (f) Call in sequence:

[0166] - PreAssemble('SMN1')

[0167] - PreAssemble('SMN2')

[0168] -AlleleImputation()

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

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

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

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

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

[0174] (4) Output results:

[0175] - 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 labels YC:Z:<R,G,B>

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

[0177] - Delete temporary sam files.

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

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

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

[0181] S300, 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 build an alignment index for the reference sequence to form an SMN1 alignment reference file.

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

[0183] 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.

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

[0185] Among them, the 3kb upstream and downstream region of SMN1 is chr5:70046668-70080595.

[0186] 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 to obtain the haplotype sequences of SMN1 and SMN2 including:

[0187] Traverse each sequence, filter out 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 variations of SMN1 and SMN2;

[0188] Record the site depth and site variation frequency, and filter the noise in the error variation-frequent region based on the interference site set, and filter the low-frequency noise variation;

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

[0190] 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 the visualization tool.

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

[0192] The preliminary grouping of SMN1 and SMN2 was used as the initial clustering, and the sequences were clustered and expanded based on all variant sites. The small or noisy groups were eliminated using the variant frequency and read coverage.

[0193] For sequences that have not yet been assigned to any haplotype, they are assigned to an existing haplotype or a newly created haplotype branch through iterative assembly and variation alignment strategies;

[0194] 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.

[0195] Example 2

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

[0197] The published Paraphase assembly process for SMN1 / SMN2 haplotypes is as follows:

[0198] Step 1: Same as step S100 in Example 1.

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

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

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

[0202] Table 1 Comparison of SMN1 / SMN2 haplotype results of Paraphase assembly

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

[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 base quality greater than Q30;

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

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

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

[0213] As can be seen from the above table, 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 quantity 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 rather than to limit it. 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 may 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 human genome reference sequence and 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 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; 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.

2. The SMN1 and SMN2 gene haplotype assembly system based on PacBio sequencing according to claim 1, characterized in that: The 1 kb upstream and downstream regions of SMN1 and SMN2 are chr5:70048668-70078595 and chr5:70924086-70954015 respectively.

3. The SMN1 and SMN2 gene haplotype assembly system based on PacBio sequencing according to claim 1, characterized in that: The 3kb upstream and downstream region of SMN1 is chr5:70046668-70080595.

4. The SMN1 and SMN2 gene haplotype assembly system based on PacBio sequencing according to claim 1, characterized in that: The assembly module includes a variation identification and sequence grouping unit, a filtering unit, a haplotype assembly unit and an output unit; The variation identification and sequence grouping unit is used to traverse each sequence, filter low-quality alignments and secondary alignments of re-aligned bam format files, identify SNVs, insertions, deletions and complex variations, and record site depth and site variation frequency, and preliminarily group reads according to the characteristic difference variation of SMN1 and SMN2; The filtering unit is used to record the site depth, 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, 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.

5. The SMN1 and SMN2 gene haplotype assembly system based on PacBio sequencing according to claim 4, characterized in that: The method of classifying the unassigned reads into corresponding haplotype branches or creating new haplotypes based on the variation frequency, reads coverage and iterative assembly strategy includes: The preliminary grouping of SMN1 and SMN2 was used as the initial clustering, and the sequences were clustered and expanded based on all variant sites. The small or noisy groups were eliminated using the variant frequency and read coverage. For sequences that have not yet been assigned to any haplotype, they are assigned to an existing haplotype or a newly created haplotype branch through iterative assembly and variation alignment strategies; 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.

6. A method for assembling haplotypes of SMN1 and SMN2 genes based on PacBio sequencing, characterized in that: include: Align the PacBio HiFi original sequencing sequence to the human genome reference sequence and generate a bam format file; Use samtools to extract intervals from the generated bam format files, extract the sequences of the 1 kb regions upstream and downstream of SMN1 and SMN2, and convert them into FASTQ format files; 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 build 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.

7. The method for assembling SMN1 and SMN2 gene haplotypes based on PacBio sequencing according to claim 6, characterized in that: The 1 kb upstream and downstream regions of SMN1 and SMN2 are chr5:70048668-70078595 and chr5:70924086-70954015 respectively.

8. The method for assembling haplotypes of SMN1 and SMN2 genes based on PacBio sequencing according to claim 6, characterized in that: The 3kb upstream and downstream region of SMN1 is chr5:70046668-70080595.

9. The method for assembling haplotypes of SMN1 and SMN2 genes based on PacBio sequencing according to claim 6, characterized in that: 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: Traverse each sequence, filter out 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 variations of SMN1 and SMN2; Record the site depth and site variation frequency, and filter the noise in the error variation-frequent region based on the interference site set, and filter the low-frequency noise variation; Based on the variation frequency, reads 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 the visualization tool.

10. The method for assembling SMN1 and SMN2 gene haplotypes based on PacBio sequencing according to claim 9, characterized in that: The method of classifying the unassigned reads into corresponding haplotype branches or creating new haplotypes based on the variation frequency, reads coverage and iterative assembly strategy includes: The preliminary grouping of SMN1 and SMN2 was used as the initial clustering, and the sequences were clustered and expanded based on all variant sites. The small or noisy groups were eliminated using the variant frequency and read coverage. For sequences that have not yet been assigned to any haplotype, they are assigned to an existing haplotype or a newly created haplotype branch through iterative assembly and variation alignment strategies; 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