Screening method of methylation sequencing fragment

Abstract

The invention relates to the field of methylation sequencing, in particular to a screening method of methylation sequencing fragments. The method comprises: 1) using a healthy sample, establishing a signal baseline by constructing an n-order Markov chain model, the healthy sample containing a methylation sequencing result, and a sequencing fragment of the healthy sample originating from an rDNA region; and 2) filtering a to-be-detected sample by using the established signal baseline, the to-be-detected sample containing a methylation sequencing result, and a sequencing fragment of the to-be-detected sample being derived from an rDNA region. According to the method disclosed by the invention, background noise from normal tissues is effectively filtered, and methylation signals (especially rDNA methylation signals) related to diseases are effectively screened out.

Description

Technical Field

[0001] The present invention relates to the field of methylation sequencing, and in particular to a method for screening methylation sequencing fragments. Background Art

[0002] DNA methylation is one of the most common epigenetic marks on the human genome, and changes in DNA methylation have been widely used as biomarkers for various biological processes or diseases. rDNA methylation is also an important target of DNA methylation. It has been confirmed that there is a potential association between rDNA methylation and cancer and aging.

[0003] Liquid biopsy is known for its minimally invasive nature, improved patient compliance, and enhanced feasibility, and is widely studied for detecting disease-related methylation markers. However, in liquid biopsy, the background noise from normal tissue is relatively large. How to screen out disease-related methylation signals (especially rDNA methylation signals) under such a large background noise is a technical problem that technicians need to solve urgently. Summary of the invention

[0004] The present invention provides a screening method for sequencing fragments to at least partially solve the above-mentioned technical problems.

[0005] At present, in the method of analyzing and aligning rDNA, the rDNA pattern sequence is directly used for the transcription region of rDNA, and the hg1 9 (human genome) + rDNA pattern sequence (U13369.1) is used for the intergenic spacer. This method has obvious defects, including that the direct alignment of the transcription region with the rDNA pattern sequence will result in many false positive alignments with short read length results; and the hg1 9 human genome and the selected pattern sequence have a ~2200-bp sequence missing.

[0006] In order to solve the problems existing in the prior art, the present invention provides an improved and optimized ribosomal DNA (rDNA) sequence alignment method: through human genome (T2T-CHM13v2.0) + r DNA consensus (KY962518.1) combined strategy, masking all the real rDNA regions on the human genome T2T as much as possible, achieving high quality r DNA comparison.

[0007] In one aspect, the present invention provides an optimized rDNA sequence alignment method. The method provided by the present invention solves the alignment difficulties caused by inaccurate sequences and tandem duplications after second-generation sequencing by constructing a new rDNA reference genome, thereby effectively obtaining rDNA sequence information and accurately identifying rDNA methylation signals.

[0008] Furthermore, the method provided by the present invention can be used for ribosomal DNA methylation signal mining and blood testing. Among them, it includes the technical difficulties of extremely low target signal circulating tumor DNA (ctDNA) content in cfDNA samples of blood of early cancer patients and excessive non-tumor-related DNA noise in the background, which is used to build corresponding algorithms and models to accurately remove cancer-related signals.

[0009] Therefore, in one aspect, the present invention provides an optimized method for ribosomal DNA (rDNA) sequence alignment, comprising constructing a new reference genome as an rDNA masking genome.

[0010] In another aspect, the present invention provides a reference genome for aligning ribosomal DNA (rDNA) sequences, which is an rDNA masking genome. Also, the present invention provides a genomic region for masking.

[0011] In another aspect, the present invention also provides a method for analyzing the rDNA methylation level of a sample, comprising using the optimized rDNA sequence alignment method of the present invention and / or a reference genome constructed thereby. And further, a method for identifying rDNA methylation markers associated with a disease is provided, comprising using the optimized rDNA sequence alignment method of the present invention and / or a reference genome constructed thereby to screen for disease-related rDNA methylation markers. r DNA methylation markers.

[0012] In yet another aspect, the present invention provides a computer program product comprising a hardware processor and a non-transitory computer-readable storage medium storing executable instructions which, when executed by the hardware processor, cause the processor to perform operations including the above-described method.

[0013] In one aspect, the present invention provides a method for screening sequencing fragments, the method comprising:

[0014] 1) Establish a signal baseline using healthy samples; and

[0015] 2) Filter the samples to be tested using the established signal baseline.

[0016] In one embodiment, the sequencing is methylation sequencing, for example, whole genome bisulfite methylation sequencing (WGBS), oxidized bisulfite sequencing (oxBS-seq), simplified genomic DNA methylation sequencing (RRBS), double restriction RRBS (dRRBS), extended bisulfite sequencing (XRBS), single cell and micro cell whole genome methylation sequencing (scWGBS), amplicon hydroxymethylation / methylation sequencing, DNA hydroxymethylation / methylation immunoprecipitation sequencing ((h)MeDIP-seq). Preferably, the sequencing fragment is whole genome bisulfite methylation sequencing (WGBS).

[0017] In one embodiment, in step 1), the number of healthy samples used can be at least 5, for example: at least 5, at least 10, at least 20, at least 50, at least 100, at least 200 or at least 1000, or a sub-range consisting of any values ​​in these ranges.

[0018] In one embodiment, in step 1), the healthy sample is derived from a mammal, for example, a human, a cow, a horse, a dog, a sheep, a dog, a mouse, etc. Preferably, the healthy sample used is derived from a human.

[0019] In one embodiment, in step 1), the healthy sample used is derived from various tissues or organs, for example, skin, muscle, blood, brain, heart, liver, spleen, lung, etc. Preferably, the healthy sample used is derived from blood tissue.

[0020] In one embodiment, in step 1), the healthy sample used is derived from methylation sequencing.

[0021] In one embodiment, in step 1), the healthy sample used contains methylation sequencing results.

[0022] In one embodiment, in step 1), the methylation sequencing results of the healthy samples are subjected to sequence alignment, for example, alignment using the rDNA sequence alignment method described in the present invention.

[0023] In one embodiment, in step 1), the methylation sequencing results of the healthy samples are subjected to sequence alignment and further data processing to obtain the fragment and CpG site distribution, for example, by processing with Methyldackel and / or PANDaseq software.

[0024] In one embodiment, step 1) includes: performing sequence alignment on the methylation sequencing results of healthy samples, and obtaining fragment and CpG site distribution through data processing based on the results of the sequence alignment.

[0025] In one embodiment, in step 1), the sequenced fragments of the healthy sample are derived from the ribosomal DNA (rDNA) region.

[0026] In one embodiment, in step 1), a signal baseline is established by constructing an n-th order Markov chain model.

[0027] In one embodiment, in step 1), the step of constructing a Markov chain model includes: mixing methylation sequencing result data from multiple (for example: at least 5 cases, at least 10 cases, at least 20 cases, at least 50 cases, at least 100 cases, at least 200 cases or at least 1000 cases, or a sub-range consisting of any values ​​in these ranges) healthy samples.

[0028] In one embodiment, in step 1), the step of constructing a Markov chain model further includes: converting the sequencing fragments in the sequencing results of the healthy samples into haplotypes in the form of vectors of the methylation status of the CpG sites contained in the corresponding sequencing fragments, one sequencing fragment corresponding to one haplotype.

[0029] In one embodiment, in step 1), the sequencing fragments in the sequencing results of the healthy sample may include sequencing fragments of autosomes and / or sex chromosomes. Preferably, the sequencing fragments in the sequencing results of the healthy sample only include sequencing fragments of autosomes.

[0030] In one embodiment, in step 1), the step of constructing a Markov chain model further comprises: filtering out haplotypes that contain less than or equal to n-1 CpG sites.

[0031] In one embodiment, in step 1), the step of constructing the Markov chain model further includes: splitting each haplotype after filtering into sub-haplotypes, so that each sub-haplotype contains only n CpG sites.

[0032] In one embodiment, in step 1), each haplotype after filtering is split into sub-haplotypes by the following steps: a moving sliding window is performed on the haplotype with n CpG sites as a window, and the step length is 1 to n CpG sites (the step length is preferably 1 CpG site), and a group of n CpG sites in each window is a sub-haplotype.

[0033] In one embodiment, in step 1), the step of constructing a Markov chain model also includes: constructing a window corresponding to each CpG site in the rDNA portion of the human genome, counting the frequency of occurrence of sub-haplotypes in healthy samples within each window, where the window length is n CpG sites; and calculating the probability of occurrence of each sub-haplotype based on the frequency according to the chain rule of probability.

[0034] In one embodiment, in step 1), 3≤n≤30, for example, n can be 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29 or 30, or a sub-range consisting of any values ​​in these ranges, preferably 3≤n≤5, and most preferably n=3.

[0035] In one embodiment, in step 1), the step of constructing the Markov chain model further includes: calculating the abnormal score AbnormalScore of each sub-haplotype based on the calculated probability of occurrence of each sub-haplotype.

[0036] In one embodiment, in step 1), the abnormal score AbnormalScore of each sub-haplotype can be calculated as follows: in the window corresponding to each CpG site constructed in the rDNA portion of the human genome, for each specific sub-haplotype, the sum of the probabilities of all other sub-haplotypes in the window whose occurrence probability is less than or equal to that of the sub-haplotype is calculated, and the abnormal score AbnormalScore = 1-the aforementioned sum of probabilities.

[0037] In one embodiment, in step 1), a threshold value of the abnormal score AbnormalScore is determined according to the calculated distribution of the abnormal score AbnormalScore of each sub-haplotype.

[0038] In one embodiment, in step 1), at least 85% (e.g., at least 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99 or 100%, or a sub-range consisting of any value in these ranges) of the sub-haplotypes have an abnormal score AbnormalScore less than or equal to the threshold. Preferably, at least 95% of the sub-haplotypes have an abnormal score AbnormalScore less than or equal to the threshold.

[0039] In one embodiment, in step 1), the threshold is at least 0.85, for example, 0.85, 0.86, 0.87, 0.88, 0.89, 0.90, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98 or 0.99, or a sub-range consisting of any values ​​in these ranges). Preferably, the threshold is 0.95.

[0040] In one embodiment, in step 2), the sample to be tested is derived from a mammal, for example, a human, a cow, a horse, a dog, a sheep, a dog, a mouse, etc. Preferably, the sample to be tested is derived from a human.

[0041] In one embodiment, in step 2), the sample to be tested is derived from various tissues or organs, for example, skin, muscle, blood, brain, heart, liver, spleen, lung, etc. Preferably, the sample to be tested is derived from blood tissue.

[0042] In one embodiment, in step 2), the sample to be tested is derived from methylation sequencing.

[0043] In one embodiment, in step 2), the sample to be tested contains methylation sequencing results.

[0044] In one embodiment, in step 2), the methylation sequencing results of the sample to be tested are subjected to sequence alignment, for example, alignment using the rDNA sequence alignment method described in the present invention.

[0045] In one embodiment, in step 2), after sequence alignment, the methylation sequencing results of the sample to be tested are further processed to obtain the fragment and CpG site distribution, for example, by processing with Methyldackel and / or PANDaseq software.

[0046] In one embodiment, the step 2) comprises: performing sequence alignment on the methylation sequencing results of the sample to be tested, and obtaining the fragment and CpG site distribution through data processing based on the sequence alignment results.

[0047] In one embodiment, in step 2), the sequencing fragments of the sample to be tested are derived from the ribosomal DNA (rDNA) region.

[0048] In one embodiment, in step 2), the sequencing fragments in the sequencing results of the sample to be tested are converted into haplotypes in the form of vectors of the methylation status of the CpG sites contained in the corresponding sequencing fragments, and one sequencing fragment corresponds to one haplotype.

[0049] In one embodiment, in step 2), the sequencing fragments in the sequencing results of the sample to be tested may include sequencing fragments of autosomes and / or sex chromosomes. Preferably, the sequencing fragments in the sequencing results of the sample to be tested only include sequencing fragments of autosomes.

[0050] In one embodiment, in step 2), haplotypes containing only less than or equal to n-1 CpG sites are filtered out.

[0051] In one embodiment, in step 2), each haplotype after filtering is split into sub-haplotypes, so that each sub-haplotype contains only n CpG sites.

[0052] In one embodiment, in step 2), each haplotype after filtering is split into sub-haplotypes by the following steps: a moving sliding window is performed on the haplotype with n CpG sites as a window, and the step length is 1 to n CpG sites (the step length is preferably 1 CpG site), and a group of n CpG sites in each window is a sub-haplotype.

[0053] In one embodiment, in step 2), a window corresponding to each CpG site is constructed in the rDNA portion of the human genome, and the frequency of occurrence of sub-haplotypes in the sample to be tested within each window is counted, where the window length is n CpG sites; and according to the chain rule of probability, the probability of occurrence of each sub-haplotype is calculated based on the frequency.

[0054] In one embodiment, in step 2), 3≤n≤30, for example, n can be 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29 or 30, or a sub-range consisting of any values ​​in these ranges, preferably 3≤n≤5, and most preferably n=3.

[0055] In one embodiment, in step 2), the abnormal score AbnormalScore of each sub-haplotype is calculated based on the calculated probability of occurrence of each sub-haplotype.

[0056] In one embodiment, in step 2), the abnormal score AbnormalScore of each sub-haplotype is calculated as follows: in the window corresponding to each CpG site constructed in the rDNA portion of the human genome, for each specific sub-haplotype, the sum of the probabilities of all other sub-haplotypes in the window whose probability of occurrence is less than or equal to that of the haplotype is calculated, and the abnormal score AbnormalScore = 1-the above-mentioned sum of probabilities.

[0057] In one embodiment, in step 2), if the abnormal score AbnormalScore of a sub-haplotype of a sequencing fragment derived from the sample to be tested is greater than the threshold determined in step 2), the sequencing fragment corresponding to the sub-haplotype is retained; correspondingly, if the abnormal scores AbnormalScore of all sub-haplotypes of a sequencing fragment derived from the sample to be tested are less than or equal to the threshold determined in step 2), the sequencing fragment corresponding to the sub-haplotype is filtered.

[0058] In one embodiment, the sample to be tested is from a patient, for example, a patient suffering from cancer.

[0059] In one embodiment, the cancer includes, for example, melanoma, non-small cell lung cancer, prostate cancer, bile duct cancer, bladder cancer, pancreatic cancer, thyroid cancer, ovarian cancer, colorectal tumors, hairy cell leukemia, acute myeloid leukemia, multiple myeloma, liver cancer, breast cancer, esophageal cancer, head and neck cancer, or glioma.

[0060] In one embodiment, the healthy sample and the sample to be tested are derived from different mammalian populations or from the same mammalian population. Preferably, the two are derived from the same mammalian population, for example, the two are derived from the same population.

[0061] In a preferred embodiment, the method comprises:

[0062] 1) Using healthy samples to establish a signal baseline, where the healthy samples contain methylation sequencing results and the sequencing fragments of the healthy samples are derived from the rDNA region; and

[0063] 2) Using the established signal baseline to filter the samples to be tested, wherein the samples to be tested contain methylation sequencing results, and the sequencing fragments of the samples to be tested are derived from the rDNA region.

[0064] In a preferred embodiment, the method comprises:

[0065] 1) Using healthy samples, a signal baseline is established by constructing an n-order Markov chain model; and

[0066] 2) Filter the samples to be tested using the established signal baseline.

[0067] In a more preferred embodiment, the method comprises:

[0068] 1) Using healthy samples, a signal baseline is established by constructing an n-order Markov chain model, wherein the healthy samples contain methylation sequencing results and the sequencing fragments of the healthy samples are derived from the rDNA region; and

[0069] 2) Using the established signal baseline to filter the samples to be tested, wherein the samples to be tested contain methylation sequencing results, and the sequencing fragments of the samples to be tested are derived from the rDNA region.

[0070] In one aspect, the present invention provides a computer program product comprising a hardware processor and a non-transitory computer-readable storage medium storing executable instructions which, when executed by the hardware processor, cause the processor to perform operations including the above-described method.

[0071] In the present invention, a CpG site refers to a cytosine (C) and a guanine (G) connected by a phosphodiester bond, wherein C is located at the 5' end and G is located at the 3' end. For ease of description, a methylated CpG site can be referred to as "M x ", the unmethylated CpG site can be recorded as "U x ”, a CpG site only has “M x ” or “U x " There are two possibilities, where the subscript indicates the ordinal number of the site, which is a natural number starting from 1 and counted from the 5' end to the 3' end. Unless otherwise specified, the subscript count is based on the counting of all CpG sites in the rDNA part of the human genome from the 5' end to the 3' end.

[0072] In the present invention, the sequencing fragments in the sequencing results are converted into a vector of the methylation status of the CpG sites contained in the corresponding sequencing fragments, which can be expressed as, for example, <M x ,M x+1 ,...,U x+n >, where M indicates that the CpG site is a methylated site; U indicates that the CpG site is a non-methylated site; the subscript is the ordinal number of the CpG site, which is a natural number starting from 1 and is counted from the 5' end to the 3' end. Unless otherwise specified, the subscript count is based on the count of all CpG sites in the rDNA portion of the human genome from the 5' end to the 3' end. For example, a sequencing fragment contains 4 CpG sites, from the 5' end to the 3' end, the first three sites are methylated sites (the first site corresponds to the first CpG site counted from the 5' end to the 3' end of the rDNA portion of the human genome), and the last one is a non-methylated site. Converting it into a vector of the methylation status of the CpG sites contained in the fragment can be expressed as,<M1,M2,M3,U4> The haplotypes and sub-haplotypes in the present invention can be expressed in this form.

[0073] In the present invention, the splitting of haplotypes can be performed as follows, for example, haplotype<M1,M2,M3,U4> , if 3 CpG sites are used as a window and the step length is 1 CpG site for sliding window, it can be split into<M1,M2,M3> and<M2,M3,U4> Two sub-haplotypes.

[0074] In the present invention, unless otherwise specified, the human genome may include any currently available human genome, for example, T2T-CHM13v2.0.

[0075] In the present invention, the methylation state of a window corresponding to a CpG site refers to the combination of the methylation states of each CpG site within a window length from the 5' end to the 3' end of the CpG site and its downstream. For example, the window length is 3 CpG sites, and the first CpG site has the following 8 possible methylation states:<M1,M2,M3> ,<U1,M2,M3> ,<U1,U2,M3> ,<U1,U2,U3> ,<M1,U2,M3> ,<M1,U2,U3> ,<M1,M2,U3> or<U1,M2,U3> Understandably,<M1,M2> Depend on<M1,M2,M3> and<M1,M2,U3> constitute;<M1,U2> Depend on<M1,U2,M3> and<M1,U2,U3> constitute;<U1,U2> Depend on<U1,U2,M3> and<U1,U2,U3> constitute;<U1,M2> Depend on<U1,M2,M3> and<U1,M2,U3> constitute; <m1>Depend on<M1,M2> (Depend on<M1,M2,M3> and<M1,M2,U3> composition) and<M1,U2> (Depend on<M1,U2,M3> and<M1,U2,U3> constitute) constitute; <u1>Depend on<U1,U2> (Depend on<U1,U2,M3> and<U1,U2,U3> composition) and<U1,M2> (Depend on<U1,M2,M3> and<U1,M2,U3> constitute) constitute.

[0076] In the present invention, "haplotype" and "sub-haplotype" can be used interchangeably depending on the specific context, and can be used to represent a combination of methylation states of CpG sites (e.g., in the same sequencing fragment or the rDNA portion of the same genome (e.g., within a certain window length)). For example, a specific methylation state in the aforementioned "methylation state of a window corresponding to a certain CpG site" can also be referred to as a "haplotype" or "sub-haplotype".

[0077] In the present invention, a "window" refers to an interval containing a fixed number of n CpG sites, which is oriented from the 5' end to the 3' end, wherein the number of CpG sites n in the window is the length of the "window".

[0078] In the present invention, when splitting haplotypes, sliding the window refers to moving the "window" from the 5' end to the 3' end by a step length (1 to n CpG sites).

[0079] In the present invention, constructing a window corresponding to each CpG site in the rDNA portion of the human genome means that all CpG sites in the rDNA portion are numbered sequentially from 1 from the 5' end to the 3' end on the human genome, and each CpG site is used as a starting point to construct a window containing the CpG site and a fixed number of CpG sites downstream from the 5' end to the 3' end. For example, when the window length is 3, the window contains the CpG site and two CpG sites downstream of it, for a total of 3 CpG sites. This process can also be viewed as: all CpG sites in the rDNA portion of the human genome are numbered sequentially from 1 from the 5' end to the 3' end, and then starting from the first CpG site, with 1 CpG site as the step length, and n (for example, 3) CpG sites as the window length are used to slide the window.

[0080] According to the chain rule of probability, the probability of the methylation state of the window corresponding to a CpG site (also called the haplotype of the window corresponding to a CpG site) can be calculated as follows (in terms of haplotype<M1,M2,U3> For example): P(M1,M2,U3)=P(U3|M1,M2)*P(M2|M1)*P(M1), where P(M1,M2,U3) is the haplotype in the window<M1,M2,U3> The probability of occurrence, P(M1) is the methylation of the first site in the window (corresponding to the first CpG site counted from the 5' end to the 3' end of the rDNA part of the human genome) <m1>), P(M2|M1) is the probability that the second site is methylated when the first site is methylated in the window (corresponding to<M1,M2> ), P(U3|M1,M2) is the probability that the third site is unmethylated when the first and second sites are methylated in the window (corresponding to<M1,M2,U3> ) probability. If not otherwise specified, in the present invention, the probability of occurrence of various events is calculated based on the frequency of the corresponding events. For example, the aforementioned probabilities can be calculated as follows: P(M1)=counts(M1) / (counts(M1)+counts(U1)); P(M2|M1)=counts(M1,M2) / (counts(M1,M2)+counts(M1,U2)); P(U3|M1,M2)=counts(M1,M2,U3) / (counts(M1,M2,U3)+counts(M1,M2,M3)). Where counts(M1) or counts(U1) is the methylation of the first site in the window (corresponding to the first CpG site counted from the 5' end to the 3' end of the rDNA part of the human genome) (corresponding to <m1>) or unmethylated (corresponding to <u1>), counts(M1,M2) or counts(M1,U2) is the frequency of methylation of the second site when the first site is methylated in the window (corresponding to<M1,M2> ) or unmethylated (corresponding to<M1,U2> ), counts(M1,M2,U3) or counts(M1,M2,M3) is the frequency of the 3rd site being unmethylated (corresponding to<M1,M2,U3> ) or methylation (corresponding to<M1,M2,M3> ) frequency.

[0081] In the present invention, the term "abnormality score" may also be referred to as "abnormality index".

[0082] Compared with the existing technology, the optimized rDNA sequence alignment method of the present invention can:

[0083] -Effectively avoid multiple alignment problems and improve the accuracy of rDNA sequence alignment.

[0084] -Reduce false positive alignments and improve alignment efficiency by masking high similarity sequences.

[0085] - Uniform coverage depth and robust methylation site identification across whole-genome and all-bisulfite sequencing samples.

[0086] - Helps to discover disease-related rDNA methylation markers.

[0087] The method of the present invention effectively filters the background noise from normal tissues and effectively screens out the methylation signals (especially rDNA methylation signals) associated with the disease. BRIEF DESCRIPTION OF THE DRAWINGS

[0088] Figure 1 : fasta file for sequence similarity analysis.

[0089] Figure 2 : Comparison of the consistency between the masked regions of rDNA on five chromosomes and the existing rDNA annotations.

[0090] Figure 3 : The comparison results on the simulated samples showed that most positions in the entire rDNA genome region could be covered uniformly.

[0091] Figure 4 : The comparison results on whole genome sequencing WGS samples show that most of the rDNA regions in real samples can have good read length alignment and uniform coverage depth.

[0092] Figure 5 :The alignment results of whole genome bisulfate sequencing (WGBS) samples showed uniform sample coverage and stable methylation site recognition in the rDNA region.

[0093] Figure 6 : Differential analysis of methylation levels and signal screening in esophageal cancer.

[0094] Figure 7 It is the distribution of abnormal scores AbnormalScore of each sub-haplotype in normal sample blood cfDNA.

[0095] Figure 8 Figure 2 is the beta distribution of WGBS sequencing data (rDNA portion) from 23 blood samples from esophageal cancer patients before and after filtering. DETAILED DESCRIPTION

[0096] It should be noted that, in the absence of conflict, the features in the embodiments and examples of the present application can be combined with each other. The present invention will be described in detail below with reference to the accompanying drawings and in combination with examples.

[0097] definition

[0098] Unless defined otherwise, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which the invention belongs. For the purpose of the present invention, the terms used herein are defined below.

[0099] Unless expressly indicated to the contrary, the indefinite articles "a" and "an" as used herein in the specification and in the claims should be understood to mean "at least one". It should also be understood that, unless expressly indicated to the contrary, in any method claimed herein that includes more than one step or action, the order of the steps or actions of the method is not necessarily limited to the order in which the steps or actions of the method are listed.

[0100] As used herein, the term "about" or "approximately" refers to a quantity, level, value, number, frequency, percentage, dimension, size, amount, weight, or length that varies by up to 15%, 10%, 9%, 8%, 7%, 6%, 5%, 4%, 3%, 2%, or 1% compared to a reference quantity, level, value, number, frequency, percentage, dimension, size, amount, weight, or length. Where a range of values ​​is provided, each value between and including the upper and lower ends of the range is specifically contemplated and set forth herein.

[0101] Throughout this specification, unless the context requires otherwise, the terms "comprising", "including", "containing" and "having" should be understood to imply the inclusion of the steps or elements or groups of steps or elements, but not the exclusion of any other steps or elements or groups of steps or elements. In specific embodiments, the terms "comprising", "including", "containing" and "having" are used synonymously. "Consisting of" means including but limited to anything following the phrase "consisting of". Therefore, the phrase "consisting of" indicates that the listed elements are required or mandatory, and no other elements may be present. As used herein, the terms "optional" and "optionally" include both selection and non-selection.

[0102] Throughout this specification, references to "one embodiment," "some embodiments," "a specific embodiment," and the like mean that a particular feature, structure, or characteristic described in connection with the embodiment is included in at least one embodiment of the present invention. Thus, the appearances of the aforementioned phrases in various places throughout this specification are not necessarily all referring to the same embodiment. Furthermore, the particular features, structures, or characteristics may be combined in any suitable manner in one or more embodiments.

[0103] As used herein, the terms "optionally" and "optionally" include both the cases of being selected and not being selected. For example, "optionally modified" includes both the cases of being modified and not being modified.

[0104] As used herein, the term "individual" refers to a human individual. The term "healthy individual" refers to an individual who is assumed not to have cancer or disease. The term "subject" refers to an individual known to have or to be likely to have cancer or disease.

[0105] As used herein, the length units nt (nucleotide) and bp (base pair) can be used interchangeably and are only used as units to represent the length of a DNA fragment or the spacing between nucleotide (or base) sites, and are not used to limit the single-stranded or double-stranded state of a DNA chain.

[0106] As used herein, the term "read" or "read length" refers to any nucleotide sequence, including a nucleotide sequence obtained from a sequence read length obtained from an individual and / or a nucleotide sequence obtained from an initial sequence read length in a sample obtained from an individual.

[0107] As used herein, the term ribosomal deoxyribonucleic acid, also known as ribosomal DNA (rDNA), is a DNA sequence located on chromosomes in the cell nucleus that encodes for rRNA.

[0108] As used herein, the term "reference genome" refers to any specific known genome sequence, whether partial or complete, of any organism or virus that can be used to reference an identified sequence from a subject. For example, reference genomes for human subjects, as well as many other organisms, are found at the National Center for Biotechnology Information at ncbi.nlm.nih.gov. "Genome" refers to the complete genetic information of an organism or virus expressed as a nucleic acid sequence. In one embodiment, the reference sequence is a reference sequence of a full-length human genome.

[0109] Samples and free DNA

[0110] As used herein, the term "circulating free DNA" or "cfDNA" refers to a plurality of nucleic acid fragments that circulate in an individual's body (e.g., bloodstream) and originate from one or more cells. In addition, cfDNA may come from other sources, such as viruses, fetuses, etc.

[0111] In some embodiments, the isolated nucleic acid, particularly cfDNA, is extracted from a sample. In some embodiments, the sample is blood. In some embodiments, the sample is plasma, such as peripheral blood plasma.

[0112] The samples described in the present invention can all be obtained by standard techniques in the art. For example, a blood sample can be obtained by standard techniques, such as using a needle and a syringe. In some embodiments, the sample is a body fluid, for example, the sample is blood. In some embodiments, the blood sample is a peripheral blood sample. Alternatively, the blood sample can be a separated portion of peripheral blood, such as a plasma sample. In some embodiments, circulating free DNA (cfDNA) is extracted from a sample, such as blood or plasma. In some embodiments, the body fluid is selected from at least one of: blood, serum, gastric juice, intestinal fluid, saliva, bile, tumor fluid, cerebrospinal fluid, breast milk, semen, urine, vaginal fluid, interstitial fluid, and feces.

[0113] Library construction and sequencing

[0114] Those skilled in the art will appreciate that library construction and sequencing can be performed using any suitable method in the art. In some embodiments, the sequencing method is a whole genome sequencing (WGS) method. In some embodiments, the sequencing method is a whole genome bisulfate sequencing (WGBS) method.

[0115] Optimized ribosomal DNA sequence alignment method and its reference genome

[0116] In some embodiments of the present invention, an optimized ribosomal DNA (rDNA) sequence alignment method is provided, comprising the following steps:

[0117] a. Use the telomere-to-telomere (T2T) human genome (2022) and rDNA model sequences (2017) as base genomes;

[0118] b. Use sequence alignment tools to perform sliding window analysis on the rDNA pattern sequence (2017) and align it with the telomere-to-telomere (T2T) human genome (2022) to search for similar sequences;

[0119] c. Based on the comparison results, screen similarity sequences that meet preset conditions and the alignment positions corresponding to the similarity sequences that meet preset conditions (for example, construct a list of rDNA high similarity sequences), and mask the corresponding alignment positions in the telomere-to-telomere (T2T) human genome (2022) (using bedtools); and

[0120] d. The rDNA pattern sequence (2017) was connected to the masked telomere-to-telomere (T2T) human genome (2022) to construct a reference genome.

[0121] In the embodiments of the present invention, a high similarity sequence refers to a position in the telomere-to-telomere (T2T) human genome (2022) where the DNA sequence and the rDNA sequence are highly similar (each base position in the human genome has a coordinate, such as: chromosome 1, 2000 to 2500 bases, usually written as chr1:2000-2500. If the interval and a certain sequence of rDNA meet the similarity requirements, it can be called a high similarity sequence, and its position is the comparison position corresponding to the high similarity sequence. Further, in the telomere-to-telomere The list of positions of all high-similarity sequences in the (T2T) human genome (2022) is called the rDNA high-similarity sequence list. The alignment positions corresponding to the similarity sequences refer to the alignment positions corresponding to the sequences that meet the similarity conditions, i.e., the coordinates of the high-similarity sequences on the telomere-to-telomere (T2T) human genome (2022). Masking the corresponding alignment positions in the telomere-to-telomere (T2T) human genome (2022) prevents its own similar sequence regions (or homologous regions) from interfering with the rDNA alignment in subsequent algorithms.

[0122] In a preferred embodiment, in the alignment method of the present invention, the telomere to telomere (T2T) human genome (2022) is the human genome (T2T-CHM13v2.0), and the rDNA pattern sequence (2017) is the rDNA consensus (KY962518.1).

[0123] In a preferred embodiment, in step b of the alignment method, the sequence alignment tool is selected from BLAST and bwa-mem bowtie2, wherein BLAST is preferably used for alignment;

[0124] The sliding window analysis of the rDNA pattern sequence (2017) included sliding window splitting with a length of 100 nt and a step size of 1 nt for comparison with the telomere-to-telomere (T2T) human genome (2022).

[0125] In a further embodiment, the sequences with high similarity are determined by the Bit-score in the results of the BLAST comparison. Generally speaking, the larger the Bit-score, the more significant the similarity between the two sequences, and the higher the similarity. For example, 130 can be selected as a threshold, that is, Bit-score>=130 is a sequence with high similarity.

[0126] In a preferred embodiment, in step b of the alignment method, the parameters of the BLAST alignment are set to:

[0127] Mode: blastn; dust-mask False; softmasking: False,

[0128] And according to the BLAST results, the parameters of high similarity sequences were selected: Bit-score>=130.

[0129] In a further embodiment, the alignment method further comprises the step of: e. verifying the improvement effect of the optimized reference genome on rDNA alignment on simulated samples and real samples.

[0130] In some embodiments of the present invention, a reference genome for aligning rDNA sequences is provided, which is an rDNA masked genome and is constructed based on the telomere-to-telomere (T2T) human genome (2022) and the rDNA pattern sequence (2017) by masking out all real rDNA regions on the telomere-to-telomere (T2T) human genome (2022) as much as possible.

[0131] In a preferred embodiment, in the reference genome of the present invention, the telomere to telomere (T2T) human genome (2022) is the human genome (T2T-CHM13v2.0), and the rDNA pattern sequence (2017) is the rDNA consensus (KY962518.1).

[0132] In a preferred embodiment, the reference genome of the present invention can be constructed by the optimized ribosomal DNA (rDNA) sequence alignment method.

[0133] In a further preferred embodiment, in the reference genome, the masked regions on the telomere-to-telomere (T2T) human genome (2022) include: chr13:5964549-9123208; chr14:2109599-2707543; chr15:2518973-3206720; chr15:3206721-3941231; c hr15:3941232-4655338; chr21:3132131-3862134; chr21:3862135-5512426; chr22: 4796484-5509436;chr13:5770549-5964549;chr13:9123208-9347588;chr14:20995 39-2109599;chr14:2707543-2817149;chr15:2506444-2518973;chr15:4655338-4 707032;chr21:3108300-3132131;chr21:5512426-5612716;chr22:4793796-479648 4; chr22:5509436-5720604; chr1:91231861-91232222; chr5:72332175-72332371; chr8:70119668-70120039; chr12:20429886-20430034; and chr12:87082020-87082523.

[0134] Method and application of the present invention

[0135] In a preferred embodiment of the present invention, the optimized rDNA sequence alignment method can be performed as follows:

[0136] A combined strategy of the completed human genome (T2T-CHM13v2.0) + rDNA consensus (KY962518.1) published in 2022 was adopted to mask out all real rDNA regions on the human genome T2T as much as possible. Under this strategy, they are divided into three categories: 1) T2T-rDNA annotation model annotated regions, 2) rDNA annotation model upstream and downstream sequences, and 3) rDNA consensus highly consistent non-other gene coding regions. The rDNA model sequence (2017) sliding window was split (length 100nt, step length 1nt), and the genome that had been masked by the strategy was aligned, all possible alignments were collected, and all perfectly matched reads (i.e., regions consistent with the rDNA sequence other than the known rDNA annotation model were found) were selected to analyze the distribution and situation of the alignments, and finally 23 rDNA / rDNA-like intervals of 3 categories (rDNA-model X8 (8,986,220bp), rDNA-model upstream and downstream X10 (940,255bp), and other repeat regions X5 (1,584bp) were selected as the mask (HardMask) region. The obtained masked region is shown in Table 1 in Example 1.

[0137] Using the ribosomal DNA constructed by the present invention to compare the reference genome, the mining of rDNA methylation signals and the signal separation of blood samples:

[0138] a. The results of the comparison are used for methylation site detection (Calling) to identify the methylation status of rDNA. The paired tissue cancer and adjacent adjacent tissue samples are processed separately to identify the areas with changed methylation status as potential markers.

[0139] b. Construct an abnormal fragment screening model based on Markov chain: Use rDNA region fragments in healthy blood samples to establish a healthy signal baseline, filter the signal of the blood of cancer patients, extract abnormal fragments, and calculate their methylation levels for subsequent prediction of cancer. For example, we used 100 healthy human plasma samples for whole genome bisulfate sequencing (WGBS) and combined them together to form a healthy baseline group. According to the methylation status of the sequencing data fragments in the baseline group, a 3rd, 4th or 5th order Markov chain model is constructed. For the fragments of cancer samples, filtering is performed based on the Markov chain model. If the threshold of the fragment P value is less than M, it is considered an abnormal fragment (M can be 0.05, 0.001, 0.0005).

[0140] Therefore, in a further embodiment, the optimized rDNA sequence alignment method and the rDNA reference genome constructed therefrom described in the present invention can be used in a method for analyzing rDNA methylation levels.

[0141] In a preferred embodiment, a method for analyzing the rDNA methylation level of a sample comprises using the ribosomal DNA (rDNA) sequence alignment method optimized by the present invention, or using the reference genome described in the present invention, and

[0142] Further including converting the reference genome into methylation patterns, and

[0143] The methylation sequencing results of the samples were compared using a reference genome of methylation patterns.

[0144] In a preferred embodiment, the bwa-meth tool is used for comparison in the method of analyzing the rDNA methylation level of a sample.

[0145] In a further preferred embodiment, the method for analyzing the rDNA methylation level of a sample further comprises calculating the CpG site methylation level using the Methyldackel tool.

[0146] In a further embodiment, the optimized rDNA sequence alignment method and the rDNA reference genome constructed by the present invention can be used to identify disease-related rDNA methylation markers. Preferably, the rDNA methylation markers can be used as biomarkers for cancer or tumor detection.

[0147] In some embodiments of the present invention, a method for identifying rDNA methylation markers associated with a disease is provided, comprising using the ribosomal DNA (rDNA) sequence alignment method optimized by the present invention, or using the reference genome described in the present invention, and

[0148] Further including converting the reference genome into methylation patterns,

[0149] Comparison of methylation sequencing results of samples using a reference genome of methylation patterns, and

[0150] The results of the comparison were used to detect methylation sites and identify the methylation status of rDNA. The paired tissue lesions and adjacent lesion samples were processed separately to identify areas with changed methylation status as markers.

[0151] In a preferred embodiment, the results of the alignment using the Methyldackel tool are used to detect methylation sites.

[0152] In a further preferred embodiment, the disease is cancer. Further preferably, identification of disease-related rDNA methylation markers can be used for early diagnosis of cancer.

[0153] In some embodiments, the present invention provides a computer program product comprising a hardware processor and a non-transitory computer-readable storage medium storing executable instructions which, when executed by the hardware processor, cause the processor to perform operations including the methods described herein.

[0154] In some embodiments, the method disclosed herein can be implemented with a computer system and / or a computer program product, which includes a computer program mechanism embedded in a computer-readable storage medium. In addition, any method of the present invention can be implemented in one or more than one computer program product or computer system. Some embodiments of the present invention provide a computer system or a computer program product, which is encoded or has instructions for performing any or all methods disclosed herein. These methods / instructions can be stored in a non-temporary computer-readable storage medium, such as a CD-ROM, a DVD, a disk storage product, or any other computer-readable data or program storage product. These methods encoded in a computer program product can also be electronically distributed via the Internet or other means by digital or carrier transmission computer data signals.

[0155] The software tools involved in the present invention may be any existing version that can achieve the purpose of the present invention, for example, the version shown in the embodiments of the present invention. Unless otherwise specified, the software setting parameters adopt the parameters recommended by the software author.

[0156] Other applications of the present invention

[0157] In some embodiments, the methods, reference genomes and / or computer program products of the present invention can be used for diagnosis of a subject's disease, such as early diagnosis or screening, such as detecting the presence of cancer or a tumor. For example, the methods, reference genomes and / or computer program products described herein can be used to analyze free DNA, which can be used for diagnosis of a subject's disease or for determining a subject's disease state, enabling higher sensitivity and accuracy. In some embodiments, determining a subject's disease state is determining whether a subject has a certain disease or the risk or degree of a subject having a certain disease, such as for a cancer patient, the progression of the cancer or the stage of the cancer can be determined.

[0158] As described above, the methods, reference genomes and / or computer program products of the present invention can be used for early detection of cancer or tumors, for example, can be used to determine early samples of cancer or tumors from subjects. In some embodiments, cancer or tumor prediction is the possibility of whether the test sample is a cancer or tumor sample.

[0159] Cancers or tumors that can be detected using the methods, reference genomes and / or computer program products of the present invention can be malignant, benign, metastatic or precancers, including but not limited to acute and chronic leukemias, lymphomas, many solid tumors of mesenchymal or epithelial tissues. More specific examples of such tumors / cancers include but are not limited to squamous cell carcinoma, skin cancer, melanoma, retinoblastoma, lung cancer (including small cell lung cancer, non-small cell lung cancer, lung adenocarcinoma and lung squamous cell carcinoma), gastric cancer, pancreatic cancer, cervical cancer, ovarian cancer, liver cancer, hepatoma, bladder cancer, testicular cancer, breast cancer, brain cancer, glioma, colorectal cancer, uterine cancer, salivary gland cancer, kidney cancer, prostate cancer, vulvar cancer, thyroid cancer, anal cancer, penis cancer, head and neck cancer, esophageal cancer, nasopharyngeal cancer, laryngeal cancer, hematological malignancies, Kaposi's sarcoma, schwannoma, oligodendroglioma, neuroblastoma, rhabdomyosarcoma, osteosarcoma, leiomyosarcoma and urinary tract cancer.

[0160] Example

[0161] The present application will illustrate the beneficial effects of the present invention through the following examples. Those skilled in the art will recognize that these examples are illustrative and not restrictive. These examples do not limit the scope of the present invention in any way. The experimental techniques and experimental methods described in the following examples, unless otherwise specified, are all conventional technical methods. For example, the experimental methods in the following examples where specific conditions are not specified are usually carried out under conventional conditions or under conditions recommended by the manufacturer; unless otherwise specified, reagents and materials are all commercially available and can be obtained through regular commercial channels.

[0162] Example 1: Construction of rDNA-masked genome

[0163] 1. Download the basic genome

[0164] The download URL of Telomere-to-Telomere (T2T) Human Genome (2022) is as follows:

[0165] https: / / www.ncbi.nlm.nih.gov / datasets / genome / GCF_009914755.1 /

[0166] The download URL of rDNA model sequence (2017) is as follows:

[0167] https: / / www.ncbi.nlm.nih.gov / nuccore / KY962518

[0168] 2. Construction of rDNA-masked genome

[0169] a. Slide the rDNA pattern sequence KY962518 with a length of 100 and a step size of 1 to generate a fasta file for sequence similarity analysis, such as Figure 1 shown.

[0170] b. Next, the basic local alignment search tool BLAST (https: / / www.ncbi.nlm.nih.gov / blast / ) was used to perform a similar sequence search.

[0171] . First, use the T2T human genome to build a reference database:

[0172] Code: ncbi-blast-2.13.0+ / bin / makeblastdb-parse_seqids-in T2T.fasta-input_type'fasta′-dbtype′nucl′)

[0173] . Then run a BLAST search against the reference database using:

[0174] Code: ncbi-blast-2.13.0+ / bin / blastn -db T2T -query KY962518.fa -outKY962518blastn.out-outfmt 6

[0175] Collect the results of the comparison search, filter the similarity sequences and comparison positions that meet the preset conditions, and build a list of rDNA high similarity sequences

[0176] Parameters: BLAST two sequence similarity Bit-score>=130 (Bit-score is defined as the log2 generalization score of the probability of finding the current item by chance in the reference database https: / / www.metagenomics.wiki / tools / blast / evalue).

[0177] According to the r-DNA high similarity sequence coordinates and the UCSC genome browser https: / / genome.ucsc.edu / annotations, the sequence distribution coordinates of rDNA on the human genome are determined as shown in Table 1:

[0178] Table 1

[0179]

[0180] The T2T human genome was masked using the maskfasta tool of bedtools, and the bases at the coordinate positions in the above table were masked as N. The rDNA pattern sequence was then connected to the masked genome to construct the T2T_masked_rDNA.fa genome.

[0181] Comparison of the masked regions of rDNA on five chromosomes with the existing rDNA annotations Figure 2 shown.

[0182] Example 2: Verifying the effect of masking genome on rDNA alignment

[0183] a. Comparison effect on simulated samples:

[0184] After sliding the rDNA pattern sequence 100bp, a fastq file similar to single-end sequencing was constructed, and the bwa mem tool was used for alignment. The BAM file of the rDNA region on the T2T_masked_rDNA genome was selected and placed in IGV (Integrative Genomics Viewer) for visualization. It can be seen that most positions in the entire rDNA genome region can be uniformly covered, and there are a few alignment gaps in the first 22kb bp (above), such as Figure 3 shown.

[0185] b. Comparison results on whole genome sequencing WGS samples:

[0186] 30 healthy (12), colorectal cancer (10), and liver cancer (8) whole genome double-end sequencing samples were merged, and the T2T_masked_rDNA genome was aligned using the bwa mem tool. The reads aligned to the rDNA region were extracted for IGV visualization observation as follows. It can be seen that most of the rDNA regions in the real samples can have good read alignment and uniform coverage depth, such as Figure 4 shown.

[0187] c. Alignment results on whole genome bisulfate sequencing WGBS samples:

[0188] The bwa-meth tool was used to convert the T2T_masked_rDNA genome into the methylation pattern T2T_masked_rDNA-methyl, and the colorectal cancer samples were aligned and the methylation sites were identified. The reads aligned to the rDNA region were extracted for IGV visualization observation as follows. It can be seen that the sample coverage in the rDNA region is uniform and the methylation site identification is stable, as shown in the figure below. Figure 5 shown.

[0189] Example 3: Differential analysis and signal screening of methylation levels in whole genome bisulfate sequencing (WGBS) cancer select

[0190] Taking esophageal cancer as an example, the paired cancer samples and adjacent paracancerous samples of esophageal cancer (Xinjiang cohort, top, 12 pairs; Xiamen cohort, bottom, 29 pairs) were aligned using the bwa-meth tool on the T2T_masked_rDNA-methyl genome, and then the Methyldackel tool was used to calculate the methylation level of the CpG sites. After that, the CpG sites in the rDNA region of the cancer and adjacent paracancerous samples were extracted and compared with adjacent paracancerous samples (first row) and cancer samples (second row) as shown below. It can be seen that after the implementation of the new alignment method, it can be successfully observed that the methylation level of the cancer samples in the middle region is lower than that of the adjacent paracancerous samples, such as Figure 6 shown.

[0191] Example 4: Screening of methylation sequencing fragments

[0192] 1) Processing of healthy samples

[0193] Using WGBS sequencing data from 200 blood samples from healthy people, the sequencing fragments located in the rDNA region were extracted and mixed. Sequence alignment was performed as described in Example 1 and Example 2 to obtain a BAM file containing the alignment results, which was processed using PANDaseq software to obtain the corresponding fragments and CpG site distribution. The specific implementation steps and codes are as follows:

[0194] (1) Starting from the BAM file, first use bedtools to split the aligned BAM file into a paired-end read fastq format.

[0195] bedtools bamtofastq -i Sample.human.bam -fq sample.human.r1.fastq -fq2 sample.human.r2.fastq

[0196] (2) Use conda, the Python environment manager, to install PANDaseq.

[0197] conda install bioconda::pandaseq

[0198] (3) Use PANDaseq to merge the split double-end fastq files and save the results.

[0199] parameter

[0200] (i) Minimum overlap length (-o) 20, which means that at least 20 bp of base regions must overlap each other when the reads of the double-end sequencing are merged.

[0201] (ii) Minimum read length (-1) 50, reads shorter than 50 bp will be filtered.

[0202] (iii) Parallel number (-T) 12, calling 12 threads for parallel speed-up.

[0203] Code

[0204] pandaseq -f sample.human.r1.fastq -r sample.human.r2.fastq -uunaligned.txt -o 20 -1 50 -T 12>sample.pandaseq.stitch.res

[0205] The sequencing fragments in the sequencing results are converted into haplotypes of the corresponding sequencing fragments, and one sequencing fragment corresponds to one haplotype. The sequencing fragments in the sequencing results only contain sequencing fragments of autosomes.

[0206] Haplotypes containing less than or equal to 2 CpG sites were filtered out.

[0207] Each haplotype after filtering is split into sub-haplotypes, so that each sub-haplotype contains only 3 CpG sites. Specifically, it can be performed as follows: a moving sliding window is performed on the haplotype with 3 CpG sites as a window, the step length is 1 CpG site, and a group of 3 CpG sites in each window is a sub-haplotype.

[0208] A window corresponding to each CpG site was constructed in the rDNA part of the human genome T2T-CHM13v2.0, and the frequency of sub-haplotypes in healthy samples within each window was counted. The window length was 3 CpG sites. The probability of occurrence of each sub-haplotype was calculated based on the frequency according to the chain rule of probability.

[0209] Then, the abnormal score of each sub-haplotype is calculated based on the calculated probability of occurrence of each sub-haplotype. The abnormal score AbnormalScore of each sub-haplotype can be calculated as follows: in the window corresponding to each CpG site constructed from the rDNA part of the human genome, for each specific sub-haplotype, the probability sum of all other sub-haplotypes in the window whose occurrence probability is less than or equal to that sub-haplotype is calculated, and the abnormal score AbnormalScore = 1-the aforementioned probability sum.

[0210] According to the distribution of abnormal scores AbnormalScore of each sub-haplotype (such as Figure 7 As shown), the threshold is determined to be 0.95, and the abnormal score AbnormalScore of at least 95% of the sub-haplotypes is less than or equal to the threshold.

[0211] 2) Processing of patient samples

[0212] Using WGBS sequencing data from 23 blood samples from esophageal cancer patients (all from East Asian populations like the aforementioned healthy individuals), the sequencing fragments located in the rDNA region were extracted and mixed. Sequence alignment was performed as described in Examples 1 and 2 to obtain a BAM file containing the alignment results, and then the Methyldackel tool (using the default parameters recommended by the software author) was used to calculate the methylation level of the CpG sites. Afterwards, the CpG sites in the rDNA region of the blood samples from healthy individuals and the blood samples from esophageal cancer patients were extracted. Based on this, a plot was made. Figure 8 The upper half of the figure: methylation levels of blood samples from healthy people (first row), methylation levels of blood samples from esophageal cancer patients (second row), superposition of methylation levels of blood samples from healthy people and esophageal cancer patients (third row), and comparison of the difference in methylation levels of blood samples from healthy people and esophageal cancer patients (fourth row).

[0213] Return to the alignment result BAM file, and use PANDaseq software (using the default parameters recommended by the software author) to process the alignment result to obtain the corresponding fragments and CpG site distribution (see step 1 for specific implementation steps and code). Convert the sequencing fragments into haplotypes of the corresponding sequencing fragments, one sequencing fragment corresponds to one haplotype. The sequencing fragments in the sequencing results only contain sequencing fragments of autosomes.

[0214] Haplotypes containing less than or equal to 2 CpG sites were filtered out.

[0215] Each haplotype after filtering is split into sub-haplotypes, so that each sub-haplotype contains only 3 CpG sites. Specifically, it can be performed as follows: a moving sliding window is performed on the haplotype with 3 CpG sites as a window, the step length is 1 CpG site, and a group of 3 CpG sites in each window is a sub-haplotype.

[0216] A window corresponding to each CpG site was constructed in the rDNA part of the human genome T2T-CHM13v2.0, and the frequency of sub-haplotypes in patient samples within each window was counted. The window length was 3 CpG sites. According to the chain rule of probability, the probability of occurrence of each sub-haplotype was calculated based on the frequency.

[0217] Then, the abnormal score AbnormalScore of each sub-haplotype is calculated based on the calculated probability of occurrence of each sub-haplotype. The abnormal score AbnormalScore of each sub-haplotype can be calculated as follows: in the window corresponding to each CpG site constructed from the rDNA part of the human genome, for each specific haplotype, the probability sum of all other haplotypes in the window whose probability of occurrence is less than or equal to that haplotype is calculated, and the abnormal score AbnormalScore = 1-the aforementioned probability sum.

[0218] If the abnormal score AbnormalScore of a sub-haplotype of a sequencing fragment derived from the sample to be tested is greater than 0.95, the sequencing fragment corresponding to the sub-haplotype is retained; accordingly, if the abnormal scores AbnormalScore of all sub-haplotypes of a sequencing fragment derived from the sample to be tested are less than or equal to 0.95, the sequencing fragment corresponding to the sub-haplotype is filtered.

[0219] The methylation levels of the CpG sites located in the rDNA of the filtered samples were calculated one by one, and then the CpG sites in the rDNA region of the healthy blood samples and the esophageal cancer patient blood samples were extracted. The horizontal axis is the CpG sites distributed on the rDNA, which are plotted from left to right in the order on the genome: the methylation level of the healthy blood samples (first row), the methylation level of the esophageal cancer patient blood samples (second row), the superposition of the methylation levels of the healthy and esophageal cancer patient blood samples (third row), and the difference in methylation levels of the healthy and esophageal cancer patient blood samples (fourth row), so as to obtain Figure 8 The lower half of the picture.

[0220] Compare Figure 8 The schematic diagram before and after filtering shows that before filtering, the mean CpG methylation values ​​of healthy people and cancer patients are similar at each site, and the difference graph (i.e. Figure 8 There is no obvious difference in the fourth row of the upper half of the figure; after filtering, the methylation levels of the blood of cancer patients are separated, and the difference map (i.e. Figure 8 There are multiple obvious difference regions on the fourth row of the lower half of the figure. This shows that after screening by the method of the present invention, the rDNA signals from normal tissues are effectively screened out, and the cancer-related rDNA signals are well enriched.

[0221] The software versions and software open source addresses used in the examples are shown in the following table (unless otherwise specified, the default parameters recommended by the software author are used).

[0222]

[0223] The foregoing detailed description refers to the examples and drawings, which illustrate specific embodiments of the present disclosure. Other embodiments with different structures and operations do not depart from the scope of the present disclosure, and reference may be made to certain specific examples of many alternative aspects or embodiments of the applicant's invention set forth in this specification, and none of them is intended to limit the scope of the applicant's invention or claims. It will be understood that various details of the present disclosure may be changed without departing from the scope of the present disclosure. In addition, the foregoing description is for illustrative purposes only, not for limiting purposes.

Claims

1. A method for screening sequencing fragments, the method comprising: 1) Using healthy samples, a signal baseline was established by constructing an n-order Markov chain model, where the healthy samples contained methylation sequencing results and the sequencing fragments of the healthy samples were derived from the rDNA region; as well as 2) Using the established signal baseline to filter the sample to be tested, wherein the sample to be tested contains methylation sequencing results, and the sequencing fragments of the sample to be tested are derived from the rDNA region.

2. The method according to claim 1, wherein step 1) further comprises: The methylation sequencing results of healthy samples were sequenced and the fragment and CpG site distribution was obtained through data processing based on the sequence alignment results; and / or The step 2) further comprises: performing sequence comparison on the methylation sequencing results of the sample to be tested, and obtaining the fragment and CpG site distribution through data processing based on the sequence comparison results.

3. The method according to any preceding claim, wherein in step 1), the step of constructing a Markov chain model comprises: Methylation sequencing results data from multiple healthy samples were mixed.

4. The method according to any of the preceding claims, wherein in step 1), the sequencing fragments in the sequencing results of the healthy sample may include sequencing fragments of autosomes and / or sex chromosomes; preferably, the sequencing fragments in the sequencing results of the healthy sample only include sequencing fragments of autosomes; and / or In step 2), the sequencing fragments in the sequencing results of the sample to be tested may include sequencing fragments of autosomes and / or sex chromosomes; preferably, the sequencing fragments in the sequencing results of the sample to be tested only include sequencing fragments of autosomes.

5. The method according to any preceding claim, wherein in step 1), the step of constructing a Markov chain model further comprises: Converting the sequencing fragments in the sequencing results of the healthy samples into haplotypes in the form of vectors of the methylation status of the CpG sites contained in the corresponding sequencing fragments, where one sequencing fragment corresponds to one haplotype; and / or In step 2), the sequencing fragments in the sequencing results of the sample to be tested are converted into haplotypes in the form of vectors of the methylation status of the CpG sites contained in the corresponding sequencing fragments, and one sequencing fragment corresponds to one haplotype.

6. The method according to claim 5, in step 1), the step of constructing a Markov chain model further comprises: Filter out haplotypes that contain less than or equal to n-1 CpG sites; and / or In step 2), the haplotypes containing only less than or equal to n-1 CpG sites are filtered out.

7. The method according to claim 6, in step 1), the step of constructing a Markov chain model further comprises: Each haplotype after filtering is split into sub-haplotypes, so that each sub-haplotype contains only n CpG sites; and / or In step 2), each haplotype after filtering is split into sub-haplotypes, so that each sub-haplotype contains only n CpG sites.

8. The method of claim 7, in step 1), each haplotype after filtering is split into sub-haplotypes by the following steps: a sliding window is performed on the haplotype with n CpG sites as a window, the step length is 1 to n CpG sites (the step length is preferably 1 CpG site), and a group of n CpG sites in each window is a sub-haplotype; and / or In step 2), each haplotype after filtering can be split into sub-haplotypes by the following steps: a moving sliding window is performed on the haplotype with n CpG sites as a window, and the step length is 1 to n CpG sites (the step length is preferably 1 CpG site), and a group of n CpG sites in each window is a sub-haplotype.

9. The method according to any preceding claim, wherein in step 1), the step of constructing a Markov chain model further comprises: A window corresponding to each CpG site was constructed in the rDNA part of the human genome, and the frequency of occurrence of sub-haplotypes in healthy samples within each window was counted. The window length was n CpG sites. The probability of occurrence of each sub-haplotype was calculated based on the frequency according to the chain rule of probability. and / or In step 2), a window corresponding to each CpG site is constructed in the rDNA part of the human genome, and the frequency of occurrence of sub-haplotypes in the sample to be tested in each window is counted, and the window length is n CpG sites; and according to the chain rule of probability, the probability of occurrence of each sub-haplotype is calculated based on the frequency.

10. The method of any preceding claim, wherein in step 1), 3≤n≤30, for example, n can be 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29 or 30, or a sub-range consisting of any value in these ranges, preferably 3≤n≤5, most preferably n=3; and / or In step 2), 3≤n≤30, for example, n can be 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29 or 30, or a sub-range consisting of any values ​​in these ranges, preferably 3≤n≤5, and most preferably n=3.

11. The method according to any preceding claim, wherein in step 1), the step of constructing a Markov chain model further comprises: Calculate the abnormal score AbnormalScore of each sub-haplotype based on the calculated probability of occurrence of each sub-haplotype; and / or In step 2), the abnormal score AbnormalScore of each sub-haplotype is calculated based on the calculated probability of occurrence of each sub-haplotype.

12. The method according to claim 11, in step 1), the abnormal score AbnormalScore of each sub-haplotype is calculated as follows: in the window corresponding to each CpG site constructed from the rDNA part of the human genome, for each specific sub-haplotype, the probability sum of all other sub-haplotypes in the window whose occurrence probability is less than or equal to that of the sub-haplotype is calculated, and the abnormal score AbnormalScore = 1-the aforementioned probability sum; and / or In step 2), the abnormal score AbnormalScore of each sub-haplotype is calculated as follows: in the window corresponding to each CpG site constructed in the rDNA part of the human genome, for each specific sub-haplotype, the sum of the probabilities of all other sub-haplotypes in the window whose probability of appearance is less than or equal to that of the sub-haplotype is calculated, and the abnormal score AbnormalScore = 1-the above-mentioned sum of probabilities.

13. The method according to claim 12, wherein in step 1), a threshold value of the abnormal score AbnormalScore is determined according to the calculated distribution of the abnormal scores AbnormalScore of each sub-haplotype.

14. The method of claim 13, wherein in step 1), at least 85% of the sub-haplotypes have an abnormal score AbnormalScore less than or equal to the threshold; preferably, at least 95% of the sub-haplotypes have an abnormal score AbnormalScore less than or equal to the threshold.

15. The method according to claim 13, in step 1), the threshold value is at least 0.85; preferably, the threshold value is 0.

95.

16. The method of claim 13, wherein in step 2), if the abnormal score AbnormalScore of a sub-haplotype of a sequencing fragment derived from the sample to be tested is greater than the threshold value determined in step 2), the sequencing fragment corresponding to the sub-haplotype is retained; correspondingly, if the abnormal scores AbnormalScore of all sub-haplotypes of a sequencing fragment derived from the sample to be tested are less than or equal to the threshold value determined in step 1), the sequencing fragment corresponding to the sub-haplotype is filtered.

17. A method according to any preceding claim, wherein the sample to be tested is from a patient, for example a patient suffering from cancer.

18. The method of claim 17, wherein the cancer comprises, for example, melanoma, non-small cell lung cancer, prostate cancer, bile duct cancer, bladder cancer, pancreatic cancer, thyroid cancer, ovarian cancer, colorectal tumors, hairy cell leukemia, acute myeloid leukemia, multiple myeloma, liver cancer, breast cancer, esophageal cancer, head and neck cancer, or glioma.

19. A method as claimed in any preceding claim, wherein the healthy sample and the sample to be tested are derived from different mammalian populations or from the same mammalian population; preferably, the healthy sample and the sample are derived from the same mammalian population, for example, the healthy sample and the sample to be tested are derived from the same human population.

20. A computer program product, comprising a hardware processor and a non-transitory computer-readable storage medium, wherein the non-transitory computer-readable storage medium stores executable instructions, which, when executed by the hardware processor, cause the processor to perform operations including the method of any one of claims 1-19.

Images