WO2020041946A1 - 基于高通量测序检测同源序列的方法和装置 - Google Patents

基于高通量测序检测同源序列的方法和装置 Download PDF

Info

Publication number
WO2020041946A1
WO2020041946A1 PCT/CN2018/102546 CN2018102546W WO2020041946A1 WO 2020041946 A1 WO2020041946 A1 WO 2020041946A1 CN 2018102546 W CN2018102546 W CN 2018102546W WO 2020041946 A1 WO2020041946 A1 WO 2020041946A1
Authority
WO
WIPO (PCT)
Prior art keywords
sequence
homologous
throughput sequencing
reads
homologous sequence
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Ceased
Application number
PCT/CN2018/102546
Other languages
English (en)
French (fr)
Inventor
张海萍
杨林
黄国栋
曾鹏
高雅
陈芳
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
BGI Shenzhen Co Ltd
Original Assignee
BGI Shenzhen Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by BGI Shenzhen Co Ltd filed Critical BGI Shenzhen Co Ltd
Priority to PCT/CN2018/102546 priority Critical patent/WO2020041946A1/zh
Priority to CN201880096241.7A priority patent/CN112513292B/zh
Publication of WO2020041946A1 publication Critical patent/WO2020041946A1/zh
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

Links

Images

Classifications

    • C—CHEMISTRY; METALLURGY
    • C12—BIOCHEMISTRY; BEER; SPIRITS; WINE; VINEGAR; MICROBIOLOGY; ENZYMOLOGY; MUTATION OR GENETIC ENGINEERING
    • C12Q—MEASURING OR TESTING PROCESSES INVOLVING ENZYMES, NUCLEIC ACIDS OR MICROORGANISMS; COMPOSITIONS OR TEST PAPERS THEREFOR; PROCESSES OF PREPARING SUCH COMPOSITIONS; CONDITION-RESPONSIVE CONTROL IN MICROBIOLOGICAL OR ENZYMOLOGICAL PROCESSES
    • C12Q1/00—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions
    • C12Q1/68—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions involving nucleic acids

Definitions

  • the invention relates to the technical field of bioinformatics, in particular to a method and a device for detecting homologous sequences based on high-throughput sequencing.
  • Homologous genes are two or more genes with sequence similarity greater than 80%. Based on the results of high-throughput sequencing data, the reads of homologous gene regions cannot be correctly aligned when they are aligned, resulting in multiple alignments. In most cases, such reads cannot accurately reflect the target position. Base case. Therefore, the correct comparison of homologous genes will encounter certain difficulties, which makes it impossible to use the existing analysis process to directly perform mutation analysis on the offline data. For example, high-throughput methods are needed to identify RHD blood types in clinical practice. Because the RHD gene has a highly homologous RHCE gene (96% similarity), the offline data cannot be correctly compared to the RHD gene. This makes some genetic diseases that contain homologous genes impossible to detect accurately.
  • Rh blood group D antigen (D antigen is defined as Rh-positive or Rh-negative if it is expressed on the erythrocyte membrane) is the main erythrocyte antigen that causes severe neonatal hemolytic disease.
  • the gene encoding the D antigen is the RHD gene.
  • an Rh-positive individual has one RHD gene [RHD heterozygote, RHD (+) / RHD (-)] or two RHD genes [RHD homozygote, RHD (+) / RHD (+)], Rh negative individuals lack the RHD gene [RHD deletion homozygote, RHD (-) / RHD (-)], some complex Rh-negative individuals often have gene fusion or some exons are missing.
  • RHD genes or RHD zygote were mainly based on the phenotype of the Rh small factor, or through indirect methods such as complex family surveys, or based on the amount of RhD antigens.
  • Direct measurement technology is a restriction fragment length polymorphism (RFLP) method. This method uses a pair of PCR primers to amplify the downstream and fusion Rh boxes simultaneously, and then uses restriction endonucleases to perform digestion. It is possible to determine the three RHD zygote types in a single experiment, but the method is complicated, time-consuming and targeted at the RHD-negative type peculiar to Caucasians.
  • RFLP restriction fragment length polymorphism
  • Hearing and speech disabilities rank first among all types of disabilities. About 1 to 3 deaf children per 1,000 newborns each year, 60% of which are related to genetic factors. There are 27.8 million hearing-impaired people in China, about 78 million people with deafness mutation carriers, and about 4 million people with drug-sensitive mutations (high-risk groups). Nearly half of the new deaf children in China each year are drug-induced deafness. . Drug-induced deafness is also inherited in maternal family populations. Therefore, early detection of drug-induced deafness genes, avoiding damage to hearing by medication, and preventing family members and children from medication-induced deafness are urgently needed problems. The gene related to drug deafness is CYP2D6.
  • CYP2D6 has a gene CYP2D7 with a similarity of 94%, it is difficult to determine using existing information analysis and comparison algorithms. Whether sequencing reads originate from CYP2D6 or CYP2D7, so a method is needed to distinguish between CYP2D6 and CYP2D7 sequencing reads, so as to accurately detect mutations that occur on CYP2D6.
  • the invention provides a method and a device for detecting homologous sequences based on high-throughput sequencing, which can solve the problem of accurate positioning of the source of homologous sequences and achieve the purpose of accurately detecting mutations.
  • an embodiment provides a method for detecting homologous sequences based on high-throughput sequencing, including:
  • the number of reads of the reference sequence from which another homologous sequence has been removed is calculated and compared, and the homologous sequence information of the sample is determined according to the number of reads.
  • the specific amplification product is a product obtained by specific amplification using multiple pairs of primers targeting multiple regions of a pair of homologous sequences.
  • the aforementioned homologous sequence is a homologous gene.
  • the aforementioned homologous genes are RHD and RHCE genes, or CYP2D6 and CYP2D7 genes.
  • the specific amplification product is a product obtained by specific amplification using multiple pairs of primers targeting multiple exon regions of a homologous gene.
  • sequence-specific site is a single nucleotide variation (SNV) site.
  • SNV single nucleotide variation
  • the above reference sequence from which another homologous sequence has been removed is a reference sequence in which another homologous sequence is completely replaced with an N base sequence.
  • the determining the homologous sequence information of the sample according to the number of reads is specifically: determining the mutation information of the homologous sequence of the sample according to the number of reads.
  • the determining the homologous sequence information of the sample according to the number of reads is specifically: determining that the homologous sequence is normal and occurs in the genome according to the difference in the number of reads of two different sources compared to the homologous sequence. Missing or duplicate conditions.
  • the alignment result is corrected according to the CIGAR value to accurately distinguish the high-throughput sequencing result according to the sequence-specific site.
  • the results of the alignment are corrected according to the CIGAR value.
  • the method further includes: removing the linker sequences at both ends of the high-throughput sequencing result.
  • the method further includes:
  • an embodiment provides a device for detecting homologous sequences based on high-throughput sequencing, including:
  • An obtaining unit configured to obtain a high-throughput sequencing result of specific amplification products of a pair of homologous sequences of a sample, where the specific amplification products include at least one sequence-specific site for distinguishing homologous sequences;
  • the first alignment unit is configured to align the above-mentioned high-throughput sequencing results with a reference sequence, and divide the above-mentioned high-throughput sequencing results into two groups according to the sequence-specific sites, each group belonging to a homologous sequence ;
  • a second alignment unit for aligning each group of high-throughput sequencing results belonging to one homologous sequence with a reference sequence from which another homologous sequence has been removed;
  • the statistics unit is configured to count the number of reads of the reference sequence from which another homologous sequence has been removed, and determine the homologous sequence information of the sample according to the number of reads.
  • an embodiment provides a computer-readable storage medium including a program that can be executed by a processor to implement the method as in the first aspect.
  • homologous sequences are distinguished by sequence-specific sites, and the high-throughput sequencing results of specific amplification products derived from the homologous sequences are divided into two groups, and then the data of each group is compared and removed separately.
  • a reference sequence of a homologous sequence so as to obtain the number of sequencing reads that belong to each type of sequence in the homologous sequence, so as to accurately distinguish the source of the sequencing sequence and then accurately detect the mutation.
  • FIG. 1 is a flowchart of a method for detecting a homologous sequence based on high-throughput sequencing according to an embodiment of the present invention
  • FIG. 2 is a structural block diagram of a device for detecting homologous sequences based on high-throughput sequencing according to an embodiment of the present invention
  • FIG. 3 is a flowchart of RHD blood group identification and analysis according to an embodiment of the present invention.
  • FIG. 4 is a diagram of a result of detecting primer uniformity in RHD blood group identification according to an embodiment of the present invention
  • FIG. 5 is a statistical diagram of the number of sequencing reads that distinguish specific sites to different genes in RHD blood group identification according to an embodiment of the present invention
  • FIG. 6 is a flowchart of CYP2D6 gene detection of drug-induced deafness according to an embodiment of the present invention
  • FIG. 7 is a diagram showing a result of detecting primer homogeneity in CYP2D6 gene detection of drug-induced deafness according to an embodiment of the present invention.
  • FIG. 8 is a statistical diagram of the number of sequencing reads that distinguish specific sites to different genes in the detection of CYP2D6 gene for drug-induced deafness in the embodiment of the present invention.
  • the present invention provides a method for detecting homologous sequences based on high-throughput sequencing, which method comprises using a specific primer to amplify a pair (2) of homologous sequences, and the primer amplification region contains At least one sequence-specific site is used to distinguish a pair of homologous sequences.
  • the first alignment is used to find the possible location of the sequencing reads.
  • the sequence-specific sites are used to distinguish different homologous sequences to distinguish the good sequencing reads.
  • Re-comparison is performed to accurately locate the mutation, and the mutation is detected through the comparison result.
  • an embodiment provides a method for detecting homologous sequences based on high-throughput sequencing, including:
  • S101 Obtain a high-throughput sequencing result of a specific amplification product of a pair of homologous sequences of a sample.
  • the specific amplification product includes at least one sequence-specific site for distinguishing the homologous sequences.
  • the “sample” refers to a sample targeted by the detection method of the present invention, which may be a clinical sample, including a healthy person sample and a patient sample, such as blood and cerebrospinal fluid samples derived from a healthy person or patient. These samples are subjected to nucleic acid (such as DNA) extraction by techniques known in the art, to obtain nucleic acid sequence fragments in the sample, and using primers that specifically target homologous sequences to amplify target fragments to obtain specific amplification products. Throughput sequencing results.
  • the sequencing platform is not limited, and may be any second-generation high-throughput sequencing platform, including but not limited to Illumina, Ion Torrent, BGISEQ or MGISEQ sequencing platform, and the like.
  • the "homologous sequence” generally refers to a sequence having a sequence similarity greater than 80%, but this is only an exemplary way of defining a homologous sequence.
  • the type of the homologous sequence is not limited, and may be a homologous gene sequence (for example, a sequence containing a readable coding frame), or a non-genetic type homologous sequence. Examples of typical but non-limiting homologous sequences are the RHD and RHCE genes, and the CYP2D6 and CYP2D7 genes.
  • the “specific amplification product” is obtained by specifically amplifying a corresponding position of a homologous sequence by a specific amplification primer.
  • the specific amplification product is the result of multiplex amplification, that is, the product obtained by specific amplification using multiple pairs of primers targeting multiple regions of a pair of homologous sequences.
  • multiplex amplification that is, the product obtained by specific amplification using multiple pairs of primers targeting multiple regions of a pair of homologous sequences.
  • specific amplification is performed using multiple pairs of primers targeting multiple exon regions of the homologous gene to obtain a specific amplification product.
  • sequence-specific site refers to a site that is different from each other at corresponding positions of a pair of homologous sequences, and the sequence-specific site contains at least 1 bp specific base, such as a single nucleotide Variation (SNV) loci.
  • sequence-specific site may be a base insertion, deletion, or copy number variation.
  • the high-throughput sequencing result refers to a result of a certain preprocessing of the offline data, for example, removing the adapter sequence at both ends of the sequencing reads of the offline data can improve the accuracy of the comparison. Rate and data validity.
  • S102 Align the high-throughput sequencing results with a reference sequence, and divide the high-throughput sequencing results into two groups based on sequence-specific sites, each group belonging to a homologous sequence.
  • the "reference sequence” generally refers to the genomic sequence of the species corresponding to the homologous sequence, such as the human reference genome sequence, and especially the human reference genome hg19.
  • the base type at the sequence-specific sites can be used as a marker to distinguish high-throughput sequencing results.
  • multiple (possibly thousands of) sequencing reads derived from a pair of homologous sequences are classified into each homologous sequence type.
  • the alignment result is corrected based on the CIGAR value to accurately distinguish the high-throughput sequencing result according to the sequence-specific site.
  • the CIGAR value correction in the present invention can improve the accuracy of sequencing reads distinguished into their respective groups, and is therefore a preferred embodiment.
  • the purpose of removing another homologous sequence from the reference sequence is to obtain the absolute position of a certain homologous sequence (such as the RHD gene) for comparison to the reference sequence, so as to avoid sequence alignment to another homology sequence.
  • the source sequence (such as the RHCE gene) causes sequencing reads to give the wrong alignment.
  • removing another homologous sequence from a reference sequence by replacing the entire homologous sequence with an N base sequence.
  • step S103 the alignment result is corrected according to the CIGAR value, which can improve the alignment accuracy of sequencing reads, and is therefore a preferred embodiment.
  • step S103 considering the influence of abnormal sequencing results on subsequent statistical accuracy, after the comparison in step S103, it also includes any one or more of the following: (a) filtering out the length of the inserted fragment greater than a preset value ( (E.g., 500) or the results of aligning the sequences at different ends to different chromosomes; (b) removing the base sites whose sequencing quality value is lower than a preset value (for example, 10). Performing the statistics in step S104 after this process can improve the accuracy of the results.
  • a preset value (E.g., 500) or the results of aligning the sequences at different ends to different chromosomes
  • a preset value for example, 10
  • S104 Count the number of reads of the reference sequence from which another homologous sequence has been removed, and determine the homologous sequence information of the sample according to the number of reads.
  • the number of reads aligned on the reference sequence reflects the existence of a homologous sequence genotype in the reference genome, such as the dose. Therefore, by counting the number of reads of the reference sequence from which another homologous sequence has been removed, it is possible to determine the homologous sequence information of the sample, such as the mutation information of the homologous sequence, such as the homologous sequence is normal in the genome and has a deletion. Or repeated situations.
  • an embodiment of the present invention provides a device for detecting homologous sequences based on high-throughput sequencing, including: an obtaining unit 201 High-throughput sequencing results of specific amplification products of a pair of homologous sequences used to obtain a sample, the specific amplification products containing at least one sequence-specific site for distinguishing homologous sequences; first alignment Unit 202 is configured to compare the above-mentioned high-throughput sequencing results with a reference sequence, and divide the above-mentioned high-throughput sequencing results into two groups according to the sequence-specific sites, each group belonging to a homologous sequence; the second An alignment unit 203 is used to compare each group of high-throughput sequencing results belonging to one homologous sequence with a reference sequence from which another homologous sequence has been removed; and a statistics unit 204 is used to perform statistical comparison to removal The number of reads of the reference sequence of another homologous sequence.
  • an embodiment of the present invention provides a computer-readable storage medium including a program, which can be executed by a processor to implement the method for detecting a homologous sequence based on the high-throughput sequencing of the present invention.
  • the program may be stored in a computer-readable storage medium.
  • the storage medium may include: a read-only memory, a random access memory, a magnetic disk, an optical disk, a hard disk, etc.
  • the computer executes the program to realize the above functions.
  • the program is stored in the memory of the device, and when the processor executes the program in the memory, all or part of the functions described above can be implemented.
  • the program may also be stored in a storage medium such as a server, another computer, a magnetic disk, an optical disk, a flash disk, or a mobile hard disk, and saved by downloading or copying.
  • a storage medium such as a server, another computer, a magnetic disk, an optical disk, a flash disk, or a mobile hard disk, and saved by downloading or copying.
  • Design 10 pairs of primers for the 10 exons of the RHD gene and the RHCE gene are designed 10 pairs of primers for the 10 exons of the RHD gene and the RHCE gene. Each pair of primers amplifies one exon region of the RHD gene and the RHCE gene, respectively. 20 products are obtained by amplification, and the 20 products are compared. For classification, accurately distinguish the source of sequencing reads according to specific sites, and then perform comparison again, and finally compare the results of the comparison to calculate the coverage depth of sequencing reads on all exons. By comparing the coverage of a certain exon of the RHD gene and the RHCE gene, to determine whether the exon of the RHD gene is deleted or duplicated.
  • This embodiment includes an experimental part and a biological information analysis part.
  • the experimental part includes: designing specific primers for homologous genes and performing multiplex PCR amplification to complete the preparation of high-throughput sequencing libraries.
  • the homologous gene regions obtained by primer amplification contain at least 1 bp differential sequences to distinguish homologous gene sequences.
  • the biological information analysis section includes: comparing the sequences amplified by the primers, finding the possible positions of the sequences through the first alignment, using differentiating sites to distinguish the sequences of different homologous genes, and performing the distinguished sequences. Re-alignment, through the detection of mutations in the results of the comparison, and comparing the coverage depth of the sequencing reads of a certain homologous region to determine whether there are deletions or duplications in that region. Mutation detection is performed through precisely located sequencing reads to accurately detect genetic disease mutations.
  • the RHD blood group identification analysis process includes:
  • the RHD gene and the RHCE gene each contain 10 exon regions, and 10 pairs of primers are designed to amplify the exon sequences corresponding to the RHD gene and the RHCE gene, respectively.
  • Each exon sequence contains at least 1 bp of specific base. Base to accurately distinguish the RHD gene from the RHCE gene.
  • a second-generation sequencing library was obtained by PCR amplification, and sequencing reads were obtained by high-throughput sequencing for each of the 10 exons.
  • the primers are processed according to the starting position of the primer, the starting position of the sequence, and the length of the primer to ensure that the primer sequence is removed most accurately, thereby retaining the most accurate true sequence information.
  • Target region primer design Primer design for the RHD gene. 10 pairs of primers cover 10 exon regions of the RHD gene. Each pair of primers simultaneously amplifies the same homologous exon of RHD and RHCE. The obtained amplification product contains at least 1bp specific sequence, used in subsequent sequencing to distinguish the source of the sequence. The primer sequence is shown in Table 1.
  • the specific primer pool 1 is obtained by mixing the above-mentioned primers in equal molar numbers.
  • the PCR amplification enzyme uses the KAPA2G Fast Multiplex PCR Kit product (Cat. No. KK5801) from the American kapa company:
  • the amplification system is shown in Table 3:
  • step 1 98 °C, 2min Step 2 98 °C, 10s Step 3 62 °C, 2min Step 4 72 °C, 30s Step 5
  • steps 2-4 15 cycles Step 6 72 °C, 5min
  • Agencourt AMPure XP magnetic beads (American Beckman Coulter Co., Ltd.) were added in 1 volume, and purified according to the instructions. After purification, the DNA was dissolved with 20 ⁇ l of distilled water.
  • the PCR amplifying enzyme used KAPA2G Fast Multiplex PCR Kit product (Cat. No. KK5801) from Kapa Company, USA.
  • step 1 98 °C, 2min Step 2 98 °C, 10s Step 3 62 °C, 2min Step 4 72 °C, 30s Step 5
  • steps 2-4 15 cycles Step 6 72 °C, 5min
  • Agencourt AMPure XP magnetic beads (American Beckman Coulter Co., Ltd.) were added in 1 volume, and purified according to the instructions. After purification, the DNA was dissolved with 20 ⁇ l of distilled water.
  • the BGISEQ-500 platform was used for sequencing, and the sequencing type was 50bp at both ends.
  • the target sequence was aligned to the human reference genome hg19 by the BWA-ALN algorithm.
  • the bwa version was 0.7.15 and the samtools version was 0.1.18.
  • the purpose is to convert the numerical expression form of the FLAG value in the second column to the letter expression form, which can be used to distinguish R1. Or R2, which is used to calculate the degree of primer enrichment for each target region.
  • the sequence is considered to be the data of the target region.
  • 3M1I46M Compared with the reference genome, the first 3 bases of this sequence can be compared to the reference genome, the fourth base is the extra base, and the fifth base can be compared to the reference genome. , So the fourth base needs to be deleted; 3M1D47M: Compared with the reference genome, the first 3 bases of this sequence can be compared to the reference genome, one base is missing at the fourth position, and the fourth base can be started from Align to the reference genome, so you need to add the letter D in the fourth position; 48M2S: Compared with the reference genome, the first 48 bases of this sequence can be compared to the reference genome, and the last 2 bases cannot be compared to the reference Genome, so the last 2 bases need to be deleted; 3S47M: Compared with the reference genome, the first 3 bases of this sequence cannot be compared to the reference genome, and the reference genome can be compared from the 4th base, so Need to delete the first 3 bases.
  • a sequence covers the position chr1: 25599086. If the sequence reads are aligned to the base G at this position, the sequence read is considered to belong to the RHD gene. If the sequence base is C, the sequence read is considered to belong to RHCE. gene.
  • the quality value corresponding to each base is ASCII converted to the corresponding decimal value, and then the corresponding quality value is obtained by subtracting 33. If the value is less than 10, the base is replaced with an *.
  • the target region sequence was counted and the difference in quantity between the two was compared (Figure 5).
  • the results showed that there was no significant difference in the depth coverage of the 10 exons of RHD and RHCE in RHD homozygous positive individual 1.
  • the 10 exons of RHD were normal, so it can be judged that the RHD gene is a homozygous RHD ( +) / RHD (+);
  • RHD exon coverage in RHD heterozygous positive individuals 2 is about half that of RHCE, and half of RHD is missing compared to RHCE, so it can be judged that this RHD gene is heterozygous RHD (+) / RHD (-);
  • the coverage of RHD exons in RHD-negative individuals 3 is almost absent, and there is almost no sequencing reads coverage compared with RHCE, so it can be judged that the 10 exons of the RHD gene are deleted and are homozygous.
  • homologous genes are captured based on multiplex PCR, and the next-generation sequencing and information analysis methods are used to accurately perform dose analysis on a homologous gene region.
  • the differences in homologous sequences can be used to determine whether genes are missing or duplicated .
  • CYP2D6 There are several SNP sites in CYP2D6 gene related to drug-induced deafness. Different SNP base information is related to drug metabolism. However, CYP2D6 has a homology gene CYP2D7 with 94% similarity. It is difficult to avoid amplification based on PCR. To CYP2D7.
  • sample 1 Five pairs of specific primers were used to detect two samples (sample 1, sample 2) of the base information of the known drug site, the amplified products were sequenced on the machine, the sequencing results were analyzed, and the specific sites were used to distinguish CYP2D6 and CYP2D7 sequencing reads, and then based on the discrimination results to accurately detect the base information of CYP2D6 and drug-related metabolic sites.
  • the detection process is shown in Figure 6.
  • Target region primer design Primer design for CYP2D6 gene, 5 pairs of primers cover 5 drug-related SNP sites (see Table 11) of CYP2D6 gene, each pair of primers simultaneously amplify the same site of CYP2D6, CYP2D7, the resulting extension
  • the amplification product contains at least a 1 bp specific sequence (Table 12) and is used for subsequent sequencing to distinguish the source of the sequence.
  • Table 12 The absolute positions and corresponding base types used to distinguish CYP2D6 and CYP2D7
  • the PCR amplifying enzyme used KAPA2G Fast Multiplex PCR Kit product (Cat. No. KK5801) from Kapa Company, USA.
  • the specific primer pool 2 is shown in Table 14:
  • the specific primer pool 2 is composed of an equal number of moles of the above primers.
  • the amplification system is shown in Table 15 below:
  • step 1 98 °C, 2min Step 2 98 °C, 10s Step 3 62 °C, 2min Step 4 72 °C, 30s Step 5
  • steps 2-4 15 cycles Step 6 72 °C, 5min
  • the PCR amplifying enzyme used KAPA2G Fast Multiplex PCR Kit product (Cat. No. KK5801) from Kapa Company, USA.
  • the general primers are shown in Table 5.
  • the amplification system is shown in Table 17 below:
  • step 1 98 °C, 2min Step 2 98 °C, 10s Step 3 62 °C, 2min Step 4 72 °C, 30s Step 5
  • steps 2-4 15 cycles Step 6 72 °C, 5min
  • Agencourt AMPure XP magnetic beads (American Beckman Coulter Co., Ltd.) were added in 1 volume, and purified according to the instructions. After purification, the DNA was dissolved with 20 ⁇ l of distilled water.
  • the BGISEQ-500 platform was used for sequencing, and the sequencing type was 50bp at both ends.
  • the target sequence was aligned to the human reference genome hg19 by the BWA-ALN algorithm.
  • the bwa version was 0.7.15 and the samtools version was 0.1.18.
  • the sequence is considered to be the data of the target region.
  • the CIGAR value is used to correct the result of the comparison and restore the reads to the original state.
  • the quality value corresponding to each base is ASCII converted to the corresponding decimal value, and then the corresponding quality value is obtained by subtracting 33. If the value is less than 10, the base is replaced with an *.
  • the target region sequence was counted and the difference in quantity between the two was compared ( Figure 8).
  • the results show that in one sample, the number of sequencing reads that are distinguished into CYP2D6 and CYP2D7 is approximately the same.
  • Table 19 shows two sample locus result detections.
  • the results show that: first, the different sequencing reads sources are distinguished by gene-specific sites, and then the target position base information is obtained based on the distinguished sequencing reads. In these two samples, the base information of the target site can be correctly identified.

Landscapes

  • Chemical & Material Sciences (AREA)
  • Organic Chemistry (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Zoology (AREA)
  • Wood Science & Technology (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Microbiology (AREA)
  • Immunology (AREA)
  • Physics & Mathematics (AREA)
  • Molecular Biology (AREA)
  • Biotechnology (AREA)
  • Biophysics (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • General Engineering & Computer Science (AREA)
  • General Health & Medical Sciences (AREA)
  • Genetics & Genomics (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

一种基于高通量测序检测同源序列的方法和装置，所述方法包括：获取样本的一对同源序列的特异性扩增产物的高通量测序结果，特异性扩增产物包含至少一个用于区分同源序列的序列特异性位点；将高通量测序结果与参考序列进行比对，并根据序列特异性位点将高通量测序结果区分为两组，每组归属于一个同源序列；将每组归属于一个同源序列的高通量测序结果分别与去除了另一个同源序列的参考序列进行比对；和统计比对到去除了另一个同源序列的参考序列的reads数目，根据reads数确定样本的同源序列信息。本发明能够解决同源序列来源的精确定位问题，实现准确检测突变的目的。

Description

基于高通量测序检测同源序列的方法和装置 技术领域
本发明涉及生物信息学技术领域，具体涉及一种基于高通量测序检测同源序列的方法和装置。
背景技术
同源基因是指序列相似度大于80％的2个或多个基因。基于高通量测序的数据结果，同源基因区域的reads在比对时，无法正确比对到合适位置而导致出现多重比对的情况，大多数情况下这样的reads没法正确反映目标位置的碱基情况。因此对于同源基因的正确比对会遇到一定困难，这使得无法使用现有的分析流程对下机数据直接进行突变分析。例如，临床上需要基于高通量的方法对RHD血型进行鉴定，由于RHD基因存在高度同源的RHCE基因(96％相似性)，下机数据无法正确比对到RHD基因上。这导致一些包含同源基因的遗传病无法进行准确的检测。
Rh血型D抗原(D抗原在红细胞膜表达与否定义为Rh阳性或Rh阴性)是引起严重新生儿溶血病的主要红细胞抗原，编码D抗原的基因为RHD基因，一般一名Rh阳性个体拥有一条RHD基因[RHD杂合子，RHD(+)/RHD(-)]或二条RHD基因[RHD纯合子，RHD(+)/RHD(+)]，Rh阴性个体则缺失RHD基因[RHD缺失纯合体，RHD(-)/RHD(-)]，有些复杂的Rh阴性个体往往出现基因融合或者出现某几个外显子缺失的现象。以往RHD基因数目或RHD合子型测定主要是根据Rh小因子表型估计，或通过复杂的家系调查，或根据RhD抗原的量等间接方法鉴别，2000年国际上建立了第一个RHD基因合子型直接测定技术，为限制性片段长度多态性(RFLP)方法，该方法采用一对PCR引物同时扩增下游和融合Rh盒子，然后采用限制性内切酶进行酶切，再电泳分析片段大小多态性，一次实验结果可以判断三种RHD合子型，但是该方法操作复杂、耗时较长且针对的是白种人特有的RHD阴性型别。2001年和2002年，中国香港学者和奥地利学者分别独立报道了采用扩增突变阻滞检测系统和实时聚合酶链式反应(Real-time PCR)技术检测已知Rh阳性表型个体的RHD基因数目，该技术同样较为复杂且要求特殊的仪器，而且不能区分RHD(+)/RHD(-)和RHD(-)/RHD(-)合子型以及缺失外显子个数及类别的详细信息。
听力语言残疾居各类残疾之首，每年1000个新生儿中约有1～3个聋儿，其中60％的聋病与遗传因素相关。中国听力障碍者2780万人，耳聋基因突变携带者约为7800万人，药物敏感性突变携带者(高危人群)约400万人，在我国每年新增的聋儿中有近半数为药物性耳聋。药物性耳聋还会在母系家族人群中遗传。因此，尽早发现药物致聋基因，避免用药损伤听力，防止家人和孩子用药致聋是急需解决的问题。与药物致聋相关的基因为CYP2D6，如果对其检测可以预防绝大部分药物导致的耳聋，但由于CYP2D6有一个相似度高达94％的基因CYP2D7，使用现有的信息分析比对算法很难确定测序reads究竟是来源于CYP2D6还是CYP2D7，因此需要一种方法来区分CYP2D6和CYP2D7的测序reads，从而准确地检测发生 在CYP2D6上的突变。
对于药物致聋基因CYP2D6用药位点检测，现有的方法可以采用质谱或sanger测序来检测，但这两种方法都遇到通量低等问题，而高通量测序能够一次对多个位点多个基因进行检测，但是高通量测序遇到的问题是无法准确区分来源于CYP2D6和CYP2D7的reads，因此检测的准确性遇到挑战。
发明内容
本发明提供一种基于高通量测序检测同源序列的方法和装置，能够解决同源序列来源的精确定位问题，实现准确检测突变的目的。
根据第一方面，一种实施例中提供一种基于高通量测序检测同源序列的方法，包括：
获取样本的一对同源序列的特异性扩增产物的高通量测序结果，上述特异性扩增产物包含至少一个用于区分同源序列的序列特异性位点；
将上述高通量测序结果与参考序列进行比对，并根据上述序列特异性位点将上述高通量测序结果区分为两组，每组归属于一个同源序列；
将每组归属于一个同源序列的高通量测序结果分别与去除了另一个同源序列的参考序列进行比对；和
统计比对到去除了另一个同源序列的参考序列的reads数目，根据上述reads数目确定上述样本的同源序列信息。
优选地，上述特异性扩增产物是使用靶向一对同源序列的多个区域的多对引物进行特异性扩增得到的产物。
优选地，上述同源序列是同源基因。
优选地，上述同源基因是RHD和RHCE基因，或CYP2D6和CYP2D7基因。
优选地，上述特异性扩增产物是使用靶向同源基因的多个外显子区域的多对引物进行特异性扩增得到的产物。
优选地，上述序列特异性位点是单核苷酸变异(SNV)位点。
优选地，上述去除了另一个同源序列的参考序列是另一个同源序列全部替换为N碱基序列的参考序列。
优选地，上述根据上述reads数目确定上述样本的同源序列信息，具体为：根据上述reads数目确定上述样本的同源序列的突变信息。
优选地，上述根据上述reads数目确定上述样本的同源序列信息，具体为：根据对比到上述同源序列的两个不同来源的reads数目的差异，确定上述同源序列在基因组上属于正常、发生缺失或重复的情况。
优选地，上述将上述高通量测序结果与参考序列进行比对后，根据CIGAR值矫正比对后的结果，以准确地根据上述序列特异性位点将上述高通量测序结果进行区分。
优选地，上述将每组归属于一个同源序列的高通量测序结果分别与去除了另一个同源序列的参考序列进行比对后，根据CIGAR值矫正比对后的结果。
优选地，上述将上述高通量测序结果与参考序列进行比对之前，还包括：去除上述高通量测序结果两端的接头序列。
优选地，上述将每组归属于一个同源序列的高通量测序结果分别与去除了另一个同源序列的参考序列进行比对之后，还包括：
过滤掉插入片段长度大于预设值或两端序列比对到不同染色体的结果；和/或
去除测序质量值低于预设值的碱基位点。
根据第二方面，一种实施例中提供一种基于高通量测序检测同源序列的装置，包括：
获取单元，用于获取样本的一对同源序列的特异性扩增产物的高通量测序结果，上述特异性扩增产物包含至少一个用于区分同源序列的序列特异性位点；
第一比对单元，用于将上述高通量测序结果与参考序列进行比对，并根据上述序列特异性位点将上述高通量测序结果区分为两组，每组归属于一个同源序列；
第二比对单元，用于将每组归属于一个同源序列的高通量测序结果分别与去除了另一个同源序列的参考序列进行比对；和
统计单元，用于统计比对到去除了另一个同源序列的参考序列的reads数目，根据上述reads数目确定上述样本的同源序列信息。
根据第三方面，一种实施例中提供一种计算机可读存储介质，包括程序，该程序能够被处理器执行以实现如第一方面的方法。
本发明的方法，通过序列特异性位点区分同源序列，将来源于同源序列的特异性扩增产物的高通量测序结果区分为两组，然后将每组数据分别比对去除了另一个同源序列的参考序列，从而得到分别属于同源序列中每一种类型序列的测序reads的数量，实现了准确区分测序序列来源，进而实现准确检测突变的目的。
附图说明
图1为本发明实施例的基于高通量测序检测同源序列的方法流程图；
图2为本发明实施例的基于高通量测序检测同源序列的装置结构框图；
图3为本发明实施例的RHD血型鉴定分析流程图；
图4为本发明实施例的RHD血型鉴定中引物均一性检测结果图；
图5为本发明实施例的RHD血型鉴定中特异性位点区分到不同基因的测序reads数目统计图；
图6为本发明实施例的药物致聋CYP2D6基因检测流程图；
图7为本发明实施例的药物致聋CYP2D6基因检测中引物均一性检测结果图；
图8为本发明实施例的药物致聋CYP2D6基因检测中特异性位点区分到不同基因的测序reads数目统计图。
具体实施方式
下面通过具体实施方式结合附图对本发明作进一步详细说明。在以下的实施方式中，很 多细节描述是为了使得本发明能被更好的理解。然而，本领域技术人员可以毫不费力的认识到，其中部分特征在不同情况下是可以省略的，或者可以由其他元件、材料、方法所替代。
另外，说明书中所描述的特点、操作或者特征可以以任意适当的方式结合形成各种实施方式。同时，方法描述中的各步骤或者动作也可以按照本领域技术人员所能显而易见的方式进行顺序调换或调整。因此，说明书和附图中的各种顺序只是为了清楚描述某一个实施例，并不意味着是必须的顺序，除非另有说明其中某个顺序是必须遵循的。
现有技术中，基于高通量测序的数据结果，同源基因区域的测序reads在比对时，无法正确比对到合适位置而导致出现多重比对的情况，大多数情况下这样的reads没法正确反映目标位置的碱基情况。因此对于同源基因的正确比对会遇到一定困难，这使得无法使用现有的分析流程对下机数据直接进行突变分析，造成对于存在同源基因的相关遗传病无法进行准确检测。
针对现有技术中的问题，本发明提供一种基于高通量测序检测同源序列的方法，该方法包括使用特异性引物扩增一对(2个)同源序列，该引物扩增区域包含至少一个序列特异性位点，用于区分一对同源序列，通过第一次比对找到测序reads可能存在的位置，通过序列特异性位点区分不同的同源序列，将区分好的测序reads进行重新比对，从而对突变进行精确的定位，通过比对结果进行突变的检测。
如图1所示，一种实施例中提供一种基于高通量测序检测同源序列的方法，包括：
S101：获取样本的一对同源序列的特异性扩增产物的高通量测序结果，上述特异性扩增产物包含至少一个用于区分同源序列的序列特异性位点。
本发明实施例中，“样本”即本发明的检测方法所针对的样本，可以是临床上的样本，包括健康人样本和病人样本，例如来源于健康人或病人的血液、脑脊液样本等。这些样本经本领域公知的技术进行核酸(如DNA)提取，获取样本中的核酸序列片段，使用特异性靶向同源序列的引物扩增目标片段，得到特异性扩增产物，测序后得到高通量测序结果。测序平台不限，可以是任何第二代高通量测序平台，包括但不限于Illumina、Ion Torrent、BGISEQ或MGISEQ测序平台等。
本发明实施例中，“同源序列”一般是指序列相似度大于80％的序列，但这仅是一种示例性的界定同源序列的方式。本发明实施例中，同源序列的类型没有限制，可以是同源基因序列(例如，含有可读编码框的序列)，也可以是非基因类型的同源序列。典型但非限定性的同源序列的例子是RHD和RHCE基因，以及CYP2D6和CYP2D7基因。
本发明实施例中，“特异性扩增产物”是通过特异性扩增引物对同源序列的相应位置进行靶向扩增得到的。针对同源序列设计特异性扩增引物，每对引物同时扩增一对同源序列的相应位置，该位置至少包含一个序列特异性位点，该位点用于准确区分序列来源，准确定位测序reads来源，达到准确检测突变的目的。
在一些实施例中，特异性扩增产物是多重扩增的结果，即使用靶向一对同源序列的多个区域的多对引物进行特异性扩增得到的产物。例如，在同源序列是同源基因的情况下，使用靶向同源基因的多个外显子区域的多对引物进行特异性扩增得到特异性扩增产物。
本发明实施例中，“序列特异性位点”是指在一对同源序列的对应位置上彼此不同的位点，该序列特异性位点至少包含1bp特异性碱基，例如单核苷酸变异(SNV)位点。在其它实施例中，序列特异性位点还可以是碱基插入、缺失或拷贝数变异等。
在一些实施例中，作为本发明的方法的输入数据，高通量测序结果是指下机数据经过一定预处理的结果，例如对下机数据去除测序reads两端的接头序列，能够提高比对准确率和数据有效性。
S102：将高通量测序结果与参考序列进行比对，并根据序列特异性位点将高通量测序结果区分为两组，每组归属于一个同源序列。
本发明实施例中，“参考序列”一般是指同源序列对应的物种的基因组序列等，例如人类参考基因组序列等，尤其是人类参考基因组hg19等。
由于同源序列在序列特异性位点上碱基不同，因此序列特异性位点上的碱基型可以作为区分高通量测序结果的标志。经过该步骤，来源于一对同源序列的多条(可能成千上万条)测序reads就被分到了每种同源序列类型中。
作为本发明的优选方案，将高通量测序结果与参考序列进行比对后，根据CIGAR值矫正比对后的结果，以准确地根据序列特异性位点将高通量测序结果进行区分。CIGAR值矫正在本发明中能够提高测序reads区分到各自所属组别中的准确度，因此是一种优选的实施方式。
S103：将每组归属于一个同源序列的高通量测序结果分别与去除了另一个同源序列的参考序列进行比对。
本发明实施例中，从参考序列中去除另一个同源序列的目的是，为了得到某一同源序列(例如RHD基因)比对到参考序列上的绝对位置，避免序列比对到另一同源序列(例如RHCE基因)而导致测序reads产生错误的比对位置。
在一些实施例中，通过将另一个同源序列全部替换为N碱基序列，来实现从参考序列中去除另一个同源序列。
类似于步骤S102，在步骤S103中，根据CIGAR值矫正比对后的结果，能够提高测序reads的比对准确度，因此是一种优选的实施方式。
在一些实施例中，考虑到非正常测序结果对后续统计准确性的影响，在步骤S103的比对之后还包括如下任一项或多项：(a)过滤掉插入片段长度大于预设值(例如500)或两端序列比对到不同染色体的结果；(b)去除测序质量值低于预设值(例如10)的碱基位点。经过该处理再进行步骤S104的统计能够提高结果准确性。
S104：统计比对到去除了另一个同源序列的参考序列的reads数目，根据reads数目确定样本的同源序列信息。
本发明实施例中，比对到参考序列上的reads数目，反映了一个同源序列基因型在参考基因组中的存在情况，例如剂量。因此，统计比对到去除了另一个同源序列的参考序列的reads数目，就能够确定样本的同源序列信息，例如同源序列的突变信息，诸如同源序列在基因组上属于正常、发生缺失或重复的情况。
在一些实施例中，所谓“根据reads数目确定样本的同源序列信息”，具体为：根据reads 数目确定样本的同源序列的突变信息。在一些实施例中，所谓“根据reads数目确定样本的同源序列信息”，具体为：根据对比到同源序列的两个不同来源的reads数目的差异，确定同源序列在基因组上属于正常、发生缺失或重复的情况。
如图2所示，对应于本发明的基于高通量测序检测同源序列的方法，本发明一种实施例中提供一种基于高通量测序检测同源序列的装置，包括：获取单元201，用于获取样本的一对同源序列的特异性扩增产物的高通量测序结果，上述特异性扩增产物包含至少一个用于区分同源序列的序列特异性位点；第一比对单元202，用于将上述高通量测序结果与参考序列进行比对，并根据上述序列特异性位点将上述高通量测序结果区分为两组，每组归属于一个同源序列；第二比对单元203，用于将每组归属于一个同源序列的高通量测序结果分别与去除了另一个同源序列的参考序列进行比对；和统计单元204，用于统计比对到去除了另一个同源序列的参考序列的reads数目，根据上述reads数目确定上述样本的同源序列信息。
相应地，本发明一种实施例中提供一种计算机可读存储介质，包括程序，该程序能够被处理器执行以实现如本发明的基于高通量测序检测同源序列的方法。
本领域技术人员可以理解，上述实施方式中各种方法的全部或部分功能可以通过硬件的方式实现，也可以通过计算机程序的方式实现。当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘、光盘、硬盘等，通过计算机执行该程序以实现上述功能。例如，将程序存储在设备的存储器中，当通过处理器执行存储器中程序，即可实现上述全部或部分功能。另外，当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序也可以存储在服务器、另一计算机、磁盘、光盘、闪存盘或移动硬盘等存储介质中，通过下载或复制保存到本地设备的存储器中，或对本地设备的系统进行版本更新，当通过处理器执行存储器中的程序时，即可实现上述实施方式中全部或部分功能。
以下通过实施例详细说明本发明的技术方案，应当理解，实施例仅是示例性的，不能理解为对本发明保护范围的限制。
实施例1：针对RHD样本进行检测
针对RHD基因和RHCE基因的10个外显子设计10对引物，每对引物分别扩增RHD基因和RHCE基因的一个外显子区域，通过扩增分别得到20个产物，对20个产物进行比对分类，根据特异性位点精确区分测序reads的来源，然后再一次进行比对，最后对比对后的结果进行统计，统计所有外显子上的测序reads覆盖深度。通过对比RHD基因和RHCE基因某一外显子上的覆盖度，来判断RHD基因的这个外显子是否发生缺失和重复，相同地，将所有外显子缺失和重复的情况进行统计，得到所有外显子的缺失情况，最后得到检测个体的合子类型RHD(+)/RHD(+)、RHD(+)/RHD(-)或RHD(-)RHD(-)，达到准确进行RHD血型鉴定的目的。
本实施例，包括实验部分和生物信息分析部分。实验部分包括：针对同源基因设计特异性引物并进行多重PCR扩增，完成高通量测序文库的制备，引物扩增得到的同源基因区域至少包含1bp差异序列，用于区分同源基因序列。生物信息分析部分包括：针对引物扩增得到的序列进行比对，通过第一次比对找到序列可能存在的位置，利用差异性位点区分不同的同源 基因的序列，将区分好的序列进行重新比对，通过比对的结果进行突变的检测，通过对比2个同源基因某个区域的测序reads覆盖深度，来判断该区域是否存在缺失或重复。通过精准定位的测序reads来进行突变检测，从而对遗传病突变进行准确检测。
具体而言，如图3所示，RHD血型鉴定分析流程包括：
(1)RHD基因和RHCE基因分别包含10个外显子区域，设计10对引物分别同时扩增RHD基因和RHCE基因对应的外显子序列，每个外显子序列包含至少1bp的特异性碱基，用于准确区分RHD基因和RHCE基因。
(2)通过PCR扩增得到二代测序上机文库，通过高通量测序以进行测序，分别得到10个外显子的测序reads信息。
(3)将测序得到的reads比对参考基因组序列，根据指定位置的特异性位点(RHD和RHCE分别对应的碱基)，精准区分RHD和RHCE的序列，将含有RHD特异性位点的reads归类为RHDreads集合，将含有RHCE特异性位点的reads归类为RHCEreads集合。
(4)将分别得到的RHD集合的reads再一次比对参考基因组序列，此时RHCE基因在人类参考基因组上的碱基序列替换为N，这样可以保证RHD的序列无法比对到RHCE上，而只能比对到RHD上，从而得到最准确的比对结果，并得到每条read在参考基因组上正确的位置信息，由于后续SNV的检测，避免RHD序列比对到RHCE基因的现象；同理得到RHCE集合的reads也再一次比对参考基因组序列，此时RHD基因上在人类参考基因组上的碱基序列全部替换为N。
(5)根据引物的起始位置和序列的起始位置以及引物的长度，来对引物进行处理，保证最准确的对引物序列进行去除，从而保留最准确的真实序列信息。
(6)通过去除低质量的位点，可以降低错误率或者其他噪音对四碱基统计结果的影响。
以下是本实施例的具体实验部分和生物信息分析部分。
采用10对特异性引物分别对已知基因型别的RHD纯合阳性个体1(2条染色体10个外显子都正常)、RHD杂合阳性个体2(一条染色体10个外显子完全缺失)、RHD阴性个体3(2条染色体完全缺失)进行PCR扩增，扩增得到的产物上机测序，对测序结果分析来判断RHD型别。
目标区域引物设计：针对RHD基因进行引物设计，10对引物覆盖RHD基因10个外显子区域，每对引物同时扩增RHD、RHCE的同一个同源外显子，得到的扩增产物至少包含1bp特异性序列，用于后续测序区分序列的来源。引物序列如表1。
表1 特异性引物池1
Figure PCTCN2018102546-appb-000001
Figure PCTCN2018102546-appb-000002
特异性引物池1由上述引物等摩尔数混合得到。
实验部分：
1.第一轮PCR扩增
PCR扩增酶采用美国kapa公司的KAPA2G Fast Multiplex PCR Kit产品(货号KK5801)：
PCR反应体系如下表2：
表2
Figure PCTCN2018102546-appb-000003
扩增体系如下表3：
表3
步骤1 98℃，2min
步骤2 98℃，10s
步骤3 62℃，2min
步骤4 72℃，30s
步骤5 重复步骤2-4，15个循环
步骤6 72℃，5min
加入1倍体积的Agencourt AMPure XP磁珠(美国贝克曼库尔特有限公司)50μl，按照说明书进行纯化，纯化后用20μl蒸馏水溶解DNA。
2.第二轮PCR扩增
PCR扩增酶采用美国kapa公司的KAPA2G Fast Multiplex PCR Kit产品(货号KK5801)。
PCR反应体系如下表4：
表4
Figure PCTCN2018102546-appb-000004
通用引物如表5。
表5 通用引物
Figure PCTCN2018102546-appb-000005
*通用引物1的5’端进行了磷酸化修饰，用于后续BGISEQ-500平台上的单链环化。
扩增体系如下表6：
表6
步骤1 98℃，2min
步骤2 98℃，10s
步骤3 62℃，2min
步骤4 72℃，30s
步骤5 重复步骤2-4，15个循环
步骤6 72℃，5min
加入1倍体积的Agencourt AMPure XP磁珠(美国贝克曼库尔特有限公司)50μl，按照说明书进行纯化，纯化后用20μl蒸馏水溶解DNA。
3.上机测序
文库质检合格后采用华大基因BGISEQ-500平台进行测序，测序类型双端50bp。
数据分析部分：
1.使用cutadapt去除两端含有接头的序列。
2.通过bwtsw算法对人类参考基因组构建比对索引，bwa版本为0.7.15。
3.通过BWA-ALN算法将目标序列比对到人类参考基因组hg19，bwa版本为0.7.15，samtools版本为0.1.18。
将比对到hg19的数据进行提取同时将生成的*.map.bam转化为*.map.flag，其目的是将第二列的FLAG值数值表达形式转变为字母表达形式，可以用于区分R1或R2，从而用于计算各目标区域的引物富集程度。
4.对每个样本的原始测序reads数、去除测序接头后的reads数目、比对率、目标区域数据所占比例等基本信息进行统计(表7)。同时对每个目标区域的富集程度进行统计，从而评估引物的均一性，图4结果显示：10对引物测序得到的深度差别不大，所有区域的深度都在 平均深度的0.4X以上，最低深度和最高深度的差值在5倍以内。
表7 下机数据统计
Figure PCTCN2018102546-appb-000006
显示：下机数据的数据利用率达到97％，比对率达到99％，目标区域比例98％。
在计算引物富集程度时，如果当前这条序列有任一1bp在目标区域内，就认为此条序列就是这个目标区域的数据。
5.通过sam文件的第六列内容，即CIGAR值将比对后的结果进行矫正，将reads还原成最原始的状态，具体矫正实例如下，结果如表8所示。
矫正实例：3M1I46M：与参考基因组相比，此条序列前3个碱基可以比对到参考基因组，第4个碱基为多出的碱基，第5个碱基开始可以比对到参考基因组，因此需要将第4个碱基删除；3M1D47M：与参考基因组相比，此条序列前3个碱基可以比对到参考基因组，第四个位置缺失一个碱基，第4个碱基开始可以比对到参考基因组，因此需要在第四个位置添加字母D；48M2S：与参考基因组相比，此条序列前48个碱基可以比对到参考基因组，最后2个碱基无法比对到参考基因组，因此需要将最后2个碱基删除；3S47M：与参考基因组相比，此条序列前3个碱基无法比对到参考基因组，从第4个碱基开始可以比对到参考基因组，因此需要将开始3个碱基删除。
表8 矫正前和矫正后结果
Figure PCTCN2018102546-appb-000007
6.用blast对RHD和RHCE两个基因的10个外显子区域进行比对，利用比对后的测序reads差异位点的序列区分RHD/RHCE，并对应生成只包含RHD数据和只包含RHCE数据的两个文件。
用于区分的绝对位置以及RHD和RHCE对应的碱基型别如下表9所示。
表9 区分RHD和RHCE来源的特异性位点及碱基
Figure PCTCN2018102546-appb-000008
例如，一条序列覆盖chr1：25599086这个位置，如果测序reads比对到这个位置出现的碱基为G，认为这条测序read属于RHD基因，如果出现的碱基是C，认为这条测序read属于RHCE基因。
7.对参考基因组序列进行处理，将RHCE基因处的碱基替换为N，记为HG19 RHCE -参考集；对参考基因组序列进行处理，将RHD基因处的碱基替换为N，记为HG19 RHD -参考集。
将步骤6中得到的RHD文件和RHCE文件分别比对HG19 RHCE -参考集和HG19 RHD -参考集。这样可以保证RHD的序列无法比对到RHCE上，而只能比对到RHD上，同理RHCE的序列 也无法比对到RHD上，从而得到最准确的比对结果，并得到每条测序read在参考基因组上正确的位置信息，由于后续SNV的检测，避免RHD序列比对到RHCE基因的现象。
8.过滤掉插入片段长度大于500或者PE测序reads比对到不同染色体的序列。
9.通过重新比对后生成的sam文件的第六列内容，即CIGAR值将比对后的结果进行矫正，将测序reads还原成最原始的状态。
10.去除低质量(测序质量值小于10)的碱基位点，同时进行标记，表10示出了一个实例。
将每个碱基对应的质量值ASCII转化为对应的十进制数值，然后减去33即可得到对应的质量值，如果这个值小于10，则将此碱基用*号代替。在python中可使用公式：碱基质量值＝ord(ASCII)–33。
表10 去低质量点实例
Figure PCTCN2018102546-appb-000009
11.对目标位点或区域进行统计。
对目标区域序列进行计数，并比较两者的数量差异(图5)。结果显示：RHD纯合阳性个体1中RHD和RHCE10个外显子的深度覆盖差别不大，与RHCE相比RHD的10个外显子都正常，因此可以判断该RHD基因是纯合的RHD(+)/RHD(+)；RHD杂合阳性个体2中RHD外显子的覆盖度是RHCE的一半左右，与RHCE相比RHD缺失一半，因此可以判断该RHD基因是杂合的RHD(+)/RHD(-)；RHD阴性个体3中RHD外显子的覆盖度基本没有，与RHCE相比几乎不存在测序reads覆盖，因此可以判断该RHD基因10个外显子缺失，且是纯合缺失RHD(-)/RHD(-)。
可见，基于多重PCR捕同源基因的特异性序列，结合二代测序和信息分析方法，来准确对某一同源基因区域进行剂量分析，通过同源序列的差异可以得到基因是否发生缺失和重复。
实施例2：针对耳聋用药位点进行检测
CYP2D6基因存在几个SNP位点与药物致聋相关，不同的SNP碱基信息和药物代谢能力相关，但CYP2D6有一个相似度为94％的同源基因CYP2D7，基于PCR的方法很难避免扩增到CYP2D7。
采用5对特异性引物分别对2个已知用药位点碱基信息的样本(样本1、样本2)进行检测，扩增得到的产物上机测序，对测序结果分析，用特异性位点区分CYP2D6和CYP2D7的测序reads，然后根据区分结果准确检测CYP2D6和药物相关代谢位点的碱基信息。检测流程如图6所示。
目标区域引物设计：针对CYP2D6基因进行引物设计，5对引物覆盖CYP2D6基因5个用药相关的SNP位点(如表11)，每对引物同时扩增CYP2D6、CYP2D7的同一个位点，得到的扩增产物至少包含1bp特异性序列(表12)，用于后续测序区分序列的来源。
表11 CYP2D6基因5个用药相关的SNP位点
rs号 碱基
rs35742686 -/A
rs3892097 A/G
rs5030865 A/T/C/G
rs1065852 C/T
rs28624811 A/T
表12 用于区分CYP2D6和CYP2D7的绝对位置及对应的碱基型别
Figure PCTCN2018102546-appb-000010
实验部分：
1.第一轮PCR扩增
PCR扩增酶采用美国kapa公司的KAPA2G Fast Multiplex PCR Kit产品(货号KK5801)。
PCR反应体系如下表13：
表13
Figure PCTCN2018102546-appb-000011
特异性引物池2如表14：
表14
Figure PCTCN2018102546-appb-000012
特异性引物池2由上述引物等摩尔数混合组成。
扩增体系如下表15：
表15
步骤1 98℃，2min
步骤2 98℃，10s
步骤3 62℃，2min
步骤4 72℃，30s
步骤5 重复步骤2-4，15个循环
步骤6 72℃，5min
加入1倍体积的Agencourt AMPure XP磁珠(美国贝克曼库尔特有限公司)50μl，按照 说明书进行纯化，纯化后用20μl蒸馏水溶解DNA。
2.第二轮PCR扩增
PCR扩增酶采用美国kapa公司的KAPA2G Fast Multiplex PCR Kit产品(货号KK5801)。
反应体系如下表16所示：
表16
Figure PCTCN2018102546-appb-000013
通用引物如表5所示。
扩增体系如下表17所示：
表17
步骤1 98℃，2min
步骤2 98℃，10s
步骤3 62℃，2min
步骤4 72℃，30s
步骤5 重复步骤2-4，15个循环
步骤6 72℃，5min
加入1倍体积的Agencourt AMPure XP磁珠(美国贝克曼库尔特有限公司)50μl，按照说明书进行纯化，纯化后用20μl蒸馏水溶解DNA。
3.上机测序
文库质检合格后采用华大基因BGISEQ-500平台进行测序，测序类型双端50bp。
数据分析部分：
1.使用cutadapt去除两端含有接头的序列。
2.通过bwtsw算法对人类参考基因组构建比对索引，bwa版本为0.7.15。
3.通过BWA-ALN算法将目标序列比对到人类参考基因组hg19，bwa版本为0.7.15，samtools版本为0.1.18。
将比对到hg19的数据进行提取；同时将生成的*.map.bam转化为*.map.flag，其目的是为了将第二列的FLAG值数值表达形式转变为字母表达形式，可以用于区分R1或R2，从而用 于计算各目标区域的引物富集程度。
4.对每个样本的原始测序reads数、去除测序接头后的reads数目、比对率、目标区域数据所占比例等基本信息进行统计(表18)。同时对每个目标区域的富集程度进行统计，从而评估引物的均一性(图7)，结果显示：5对引物测序得到的深度差别不大，所有区域的深度都在平均深度的0.5X以上，最低深度和最高深度的差值在3倍以内。
表18 CYP2D6基因检测下机数据基本信息统计
Figure PCTCN2018102546-appb-000014
显示：下机数据的数据利用率达到98％，比对率达到99％，目标区域比例99％。
在计算引物富集程度时，如果当前这条序列有任一1bp在目标区域内，就认为此条序列就是这个目标区域的数据。
5.通过sam文件的第六列内容，即CIGAR值将比对后的结果进行矫正，将reads还原成最原始的状态。
6.通过比对后的reads差异位点的序列区分CYP2D6和CYP2D7，并对应生成只包含CYP2D6数据和只包含CYP2D7数据的两个文件。
对参考基因组序列进行处理，将CYP2D6基因处的碱基替换为N，记为HG19 CYP2D6 -参考集；对参考基因组序列进行处理，将CYP2D7基因处的碱基替换为N，记为HG19 CYP2D7 -参考集。
7.将步骤6中得到的CYP2D6文件和CYP2D7文件分别比对HG19 CYP2D7 -参考集和HG19 CYP2D6 -参考集。
8.过滤掉插入片段长度大于500或者PE测序reads比对到不同染色体的序列。
9.通过重新比对后生成的sam文件的第六列内容，即CIGAR值将比对后的结果进行矫正，将测序reads还原成最原始的状态。
10.去除低质量(测序质量值小于10)的碱基位点，同时进行标记。
将每个碱基对应的质量值ASCII转化为对应的十进制数值，然后减去33即可得到对应的质量值，如果这个值小于10，则将此碱基用*号代替。在python中可使用公式：碱基质量值＝ord(ASCII)–33。
11.对目标位点或区域进行统计。
对目标区域序列进行计数，并比较两者的数量差异(图8)。结果显示：在一个样本中，区分到CYP2D6和CYP2D7的测序reads数目大致相当。
表19示出了两个样本位点结果检测。
表19 位点结果检测(CYP2D6基因检测)
Figure PCTCN2018102546-appb-000015
Figure PCTCN2018102546-appb-000016
结果表明：先通过基因特异性位点区分不同的测序reads来源后，根据区分后的测序reads来得到目标位置碱基信息，在这两个样本中，可以正确鉴定目标位点的碱基信息。
以上应用了具体个例对本发明进行阐述，只是用于帮助理解本发明，并不用以限制本发明。对于本发明所属技术领域的技术人员，依据本发明的思想，还可以做出若干简单推演、变形或替换。

Claims (15)

  1. 一种基于高通量测序检测同源序列的方法，其特征在于，所述方法包括：
    获取样本的一对同源序列的特异性扩增产物的高通量测序结果，所述特异性扩增产物包含至少一个用于区分同源序列的序列特异性位点；
    将所述高通量测序结果与参考序列进行比对，并根据所述序列特异性位点将所述高通量测序结果区分为两组，每组归属于一个同源序列；
    将每组归属于一个同源序列的高通量测序结果分别与去除了另一个同源序列的参考序列进行比对；和
    统计比对到去除了另一个同源序列的参考序列的reads数目，根据所述reads数目确定所述样本的同源序列信息。
  2. 根据权利要求1所述的方法，其特征在于，所述特异性扩增产物是使用靶向一对同源序列的多个区域的多对引物进行特异性扩增得到的产物。
  3. 根据权利要求1所述的方法，其特征在于，所述同源序列是同源基因。
  4. 根据权利要求1所述的方法，其特征在于，所述同源基因是RHD和RHCE基因，或CYP2D6和CYP2D7基因。
  5. 根据权利要求3或4所述的方法，其特征在于，所述特异性扩增产物是使用靶向所述同源基因的多个外显子区域的多对引物进行特异性扩增得到的产物。
  6. 根据权利要求1所述的方法，其特征在于，所述序列特异性位点是单核苷酸变异(SNV)位点。
  7. 根据权利要求1所述的方法，其特征在于，所述去除了另一个同源序列的参考序列是另一个同源序列全部替换为N碱基序列的参考序列。
  8. 根据权利要求1所述的方法，其特征在于，所述根据所述reads数目确定所述样本的同源序列信息，具体为：根据所述reads数目确定所述样本的同源序列的突变信息。
  9. 根据权利要求1所述的方法，其特征在于，所述根据所述reads数目确定所述样本的同源序列信息，具体为：根据对比到所述同源序列的两个不同来源的reads数的差异，确定所述同源序列在基因组上属于正常、发生缺失或重复的情况。
  10. 根据权利要求1所述的方法，其特征在于，所述将所述高通量测序结果与参考序列进行比对后，根据CIGAR值矫正比对后的结果，以准确地根据所述序列特异性位点将所述高通量测序结果进行区分。
  11. 根据权利要求1所述的方法，其特征在于，所述将每组归属于一个同源序列的高通量测序结果分别与去除了另一个同源序列的参考序列进行比对后，根据CIGAR值矫正比对后的结果。
  12. 根据权利要求1所述的方法，其特征在于，所述将所述高通量测序结果与参考序列进行比对之前，还包括：去除所述高通量测序结果两端的接头序列。
  13. 根据权利要求1所述的方法，其特征在于，所述将每组归属于一个同源序列的高通量测序结果分别与去除了另一个同源序列的参考序列进行比对之后，还包括：
    过滤掉插入片段长度大于预设值或两端序列比对到不同染色体的结果；和/或
    去除测序质量值低于预设值的碱基位点。
  14. 一种基于高通量测序检测同源序列的装置，其特征在于，所述装置包括：
    获取单元，用于获取样本的一对同源序列的特异性扩增产物的高通量测序结果，所述特异性扩增产物包含至少一个用于区分同源序列的序列特异性位点；
    第一比对单元，用于将所述高通量测序结果与参考序列进行比对，并根据所述序列特异性位点将所述高通量测序结果区分为两组，每组归属于一个同源序列；
    第二比对单元，用于将每组归属于一个同源序列的高通量测序结果分别与去除了另一个同源序列的参考序列进行比对；和
    统计单元，用于统计比对到去除了另一个同源序列的参考序列的reads数目，根据所述reads数目确定所述样本的同源序列信息。
  15. 一种计算机可读存储介质，其特征在于，包括程序，所述程序能够被处理器执行以实现如权利要求1-13中任一项所述的方法。
PCT/CN2018/102546 2018-08-27 2018-08-27 基于高通量测序检测同源序列的方法和装置 Ceased WO2020041946A1 (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
PCT/CN2018/102546 WO2020041946A1 (zh) 2018-08-27 2018-08-27 基于高通量测序检测同源序列的方法和装置
CN201880096241.7A CN112513292B (zh) 2018-08-27 2018-08-27 基于高通量测序检测同源序列的方法和装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/CN2018/102546 WO2020041946A1 (zh) 2018-08-27 2018-08-27 基于高通量测序检测同源序列的方法和装置

Publications (1)

Publication Number Publication Date
WO2020041946A1 true WO2020041946A1 (zh) 2020-03-05

Family

ID=69642718

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2018/102546 Ceased WO2020041946A1 (zh) 2018-08-27 2018-08-27 基于高通量测序检测同源序列的方法和装置

Country Status (2)

Country Link
CN (1) CN112513292B (zh)
WO (1) WO2020041946A1 (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115851895A (zh) * 2022-12-23 2023-03-28 中国医学科学院输血研究所 基于三代测序的RHD、RHCE mRNA全长的测定方法及试剂盒
CN120340601A (zh) * 2025-06-20 2025-07-18 杭州华大序风科技有限公司 基因变异检测方法及装置、电子设备及存储介质

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114457169A (zh) * 2022-03-07 2022-05-10 南京鼓楼医院 一种基于高通量测序的RhD基因型检测方法
CN114944188B (zh) * 2022-05-19 2026-02-06 广州微远基因科技有限公司 样本同源性判定模型及其建立方法和应用
CN116959580B (zh) * 2023-08-30 2025-09-19 予果生物科技（北京）有限公司 一种基于靶向高通量测序序列的比对方法及其应用
CN118703609A (zh) * 2024-04-08 2024-09-27 上海荷谱诊断技术有限公司 一种同时测定cyp2d6基因多态性和拷贝数的方法、试剂盒及系统
CN119920328B (zh) * 2025-03-31 2025-09-26 中国科学院微生物研究所 一种病毒物种鉴定方法、鉴定系统、设备和介质

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102965367A (zh) * 2012-12-04 2013-03-13 中国农业科学院棉花研究所 一种获得植物候选抗病基因序列的方法
CN105112569A (zh) * 2015-09-14 2015-12-02 中国医学科学院病原生物学研究所 基于宏基因组学的病毒感染检测及鉴定方法
WO2016023962A1 (en) * 2014-08-13 2016-02-18 Progenika Biopharma S.A. Consensus-based allele detection

Family Cites Families (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2003056030A2 (en) * 2001-11-08 2003-07-10 The Johns Hopkins University Methods and systems of nucleic acid sequencing
CN104531883B (zh) * 2015-01-14 2018-02-02 北京圣谷同创科技发展有限公司 Pkd1基因突变的检测试剂盒及检测方法
CN106367475B (zh) * 2015-07-23 2019-08-30 上海生物信息技术研究中心 一种mmr基因突变检测试剂盒
CN105567830A (zh) * 2016-01-29 2016-05-11 江汉大学 一种植物转基因成分的检测方法
WO2017156290A1 (en) * 2016-03-09 2017-09-14 Baylor College Of Medicine A novel algorithm for smn1 and smn2 copy number analysis using coverage depth data from next generation sequencing
CN107653299A (zh) * 2016-07-23 2018-02-02 成都十洲科技有限公司 一种基于高通量测序的基因芯片探针序列的获取方法
CN107688727B (zh) * 2016-08-05 2020-07-14 深圳华大基因股份有限公司 生物序列聚类和全长转录组中转录本亚型识别方法和装置
CN106372459B (zh) * 2016-08-30 2019-03-15 天津诺禾致源生物信息科技有限公司 一种基于扩增子二代测序拷贝数变异检测的方法及装置
CN106282356B (zh) * 2016-08-30 2019-11-26 天津诺禾医学检验所有限公司 一种基于扩增子二代测序点突变检测的方法及装置
WO2018112249A1 (en) * 2016-12-15 2018-06-21 Illumina, Inc. Methods and systems for determining paralogs
CN108103204B (zh) * 2017-12-15 2018-10-02 东莞博奥木华基因科技有限公司 基于多重PCR及二代测序的Rh血型分型方法及装置
CN108424907B (zh) * 2018-05-09 2021-10-15 北京大学 一种高通量dna多位点精确碱基突变方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102965367A (zh) * 2012-12-04 2013-03-13 中国农业科学院棉花研究所 一种获得植物候选抗病基因序列的方法
WO2016023962A1 (en) * 2014-08-13 2016-02-18 Progenika Biopharma S.A. Consensus-based allele detection
CN105112569A (zh) * 2015-09-14 2015-12-02 中国医学科学院病原生物学研究所 基于宏基因组学的病毒感染检测及鉴定方法

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115851895A (zh) * 2022-12-23 2023-03-28 中国医学科学院输血研究所 基于三代测序的RHD、RHCE mRNA全长的测定方法及试剂盒
CN115851895B (zh) * 2022-12-23 2025-03-21 中国医学科学院输血研究所 基于三代测序的RHD、RHCE mRNA全长的测定方法及试剂盒
CN120340601A (zh) * 2025-06-20 2025-07-18 杭州华大序风科技有限公司 基因变异检测方法及装置、电子设备及存储介质
CN120340601B (zh) * 2025-06-20 2025-09-26 杭州华大序风科技有限公司 基因变异检测方法及装置、电子设备及存储介质

Also Published As

Publication number Publication date
CN112513292B (zh) 2023-12-26
CN112513292A (zh) 2021-03-16

Similar Documents

Publication Publication Date Title
WO2020041946A1 (zh) 基于高通量测序检测同源序列的方法和装置
JP7637139B2 (ja) がん予測パイプラインにおけるｒｎａ発現コールを自動化するためのシステムおよび方法
Aune et al. Expression of long non-coding RNAs in autoimmunity and linkage to enhancer function and autoimmune disease risk genetic variants
Sheng et al. Multi-perspective quality control of Illumina RNA sequencing data analysis
CN106715711B (zh) 确定探针序列的方法和基因组结构变异的检测方法
TWI636255B (zh) 癌症檢測之血漿ｄｎａ突變分析
Guo et al. Three-stage quality control strategies for DNA re-sequencing data
US12106825B2 (en) Computational modeling of loss of function based on allelic frequency
WO2015149034A9 (en) Gene fusions and gene variants associated with cancer
WO2015042980A1 (zh) 确定染色体预定区域中snp信息的方法、系统和计算机可读介质
Fang et al. DNA methylation entropy is associated with DNA sequence features and developmental epigenetic divergence
US10787708B2 (en) Method of identifying a gene associated with a disease or pathological condition of the disease
Kubiritova et al. On the critical evaluation and confirmation of germline sequence variants identified using massively parallel sequencing
WO2020047694A1 (zh) 确定新发突变在胚胎中的遗传状态的方法和装置
CN111508561A (zh) 同源序列和同源序列中串联重复序列的检测方法、计算机可读介质和应用
Yamamoto et al. Functional landscape of genome-wide postzygotic somatic mutations between monozygotic twins
WO2024137407A1 (en) Methods and targets of dna methylation entropy
CN116287309A (zh) 鉴定鹅繁殖障碍的分子标记、引物、pcr方法和应用
JP2020520679A (ja) 無細胞核酸から得られた配列分析データに係わる背景対立因子の頻度分布を生成する方法、及びそれを利用して無細胞核酸から変異を検出する方法
Borràs et al. The use of transcriptomics in clinical applications
Meng Ethics statement
Lobon Garcia et al. Somatic mutations detected in Parkinson disease could affect genes with a role in synaptic and neuronal processes
Devine et al. Xuefang Zhao, Ryan L. Collins, Wan-Ping Lee, 5 Alexandra M. Weber, 6, 7 Yukyung Jun, 5 Qihui Zhu, 5 Ben Weisburd, 2 Yongqing Huang, 8 Peter A. Audano, 9 Harold Wang, Mark Walker, 2, 3 Chelsea Lowther, Jack Fu, Human Genome Structural Variation Consortium, Mark B. Gerstein, 10
Vermeulen Improving and estimating Y chromosome loss in blood and brain tissues using high-throughput sequencing
Culibrk Copy number variation in metastatic cancer: methods and analysis of somatic copy number variation in advanced human cancers

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 18931705

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 18931705

Country of ref document: EP

Kind code of ref document: A1