Patent Application
20230126920
The present invention pertains to the fields of genomics and bioinformatics, and relates to a classification method, device and use of urine sediment genomic DNA.
Urogenital tumors refer to tumors that occur in the urinary system. Common urogenital tumors include renal cancer (RC), bladder tumor (BT) and prostate cancer (PCA). The Cancer Statistics Report in 2018 shows that, among the top 20 common tumors in terms of new cases and death cases, there are three urogenital tumors and PCA is in top three.
Most of the patients with early-stage tumors can be radically cured by surgeries, but the prognosis and survival of patients are significantly reduced once metastases occur. Currently, the diagnosis of urogenital tumors mainly relies on tissue biopsies, while non-invasive diagnosis is immature, and the sensitivity and specificity in tumor detection are not high.
Renal cell carcinoma is also known as renal cancer, and a common subtype is kidney renal clear cell carcinoma, accounting for about 80-85% of renal cancer. The main types of renal cancer include kidney renal clear cell carcinoma, papillary renal cell carcinoma, and chromophobe renal cell carcinoma, which together account for about 95% of renal cancer. Due to lack of good markers for early diagnosis, renal cell carcinoma has progressed to advanced stages at the time of diagnosis in many patients.
Currently, the clinically recognized “gold standard” for the diagnosis and follow-up of BT relies on the combination of cystoscopy with pathological examination on shed cells in urine. The entire bladder can be examined by cystoscopy, but cystoscopy has a low diagnostic sensitivity (52%-68%) for high-grade bladder carcinoma in situ. In addition, the friction of the instrument against the urethra during the examination can easily lead to urothelial injury to a patient, resulting in a strong sense of pain to the patient. The diagnostic sensitivity of pathological examination on shed cells in urine is low, especially for BT with low pathological grade (4%-31%).
Prostate specific antibody (PSA) tests are widely used in the process of early diagnosis of prostate cancer. However, the PSA variation is susceptible to many factors, making its accuracy not high. Furthermore, prior to paracentesis, the selective use of multi-parameter parametric magnetic imaging (mpMRI) may improve the detection rate of prostate cancer (Gleason score >7). However, the use of mpMRI is controversial, and further diagnosis must rely on pathological diagnosis.
Liquid biopsy refers to a technique for detecting dynamic changes in tumors by using circulating tumor cells (CTCs), cell-free tumor DNAs, and exosomes released by tumor tissue into body fluids such as blood and urine. Due to its non-invasive or minimally invasive, real-time and dynamic characteristics, liquid biopsy has been widely used in the research of early diagnosis, metastasis, prognosis judgment, mechanisms of forming drug resistance and personalized treatment guidance of tumors. Currently, most of the studies on liquid biopsy mainly use blood as a carrier. In fact, the advantage of urine over blood is pronounced, i.e. truly non-invasive.
However, similar to liquid biopsy which uses blood as a carrier, urine-based liquid biopsy technology faces the problem of how to make use of a limited signal to trace the origin of a tumor tissue due to the low level of signal released by urogenital tumors. Currently, genomic variation tracing based on NGS technology has been reported, including driver gene mutations, and insertions and deletions. However, tumors are highly heterogeneous, and the driver gene variation may not be detected in shed cells. Furthermore, the identification of a mutation in a small number of tumor cfDNAs relies on targeted deep sequencing (>5000*) which may have sequencing errors.
At present, there is still a need to develop new means having good specificity and sensitivity for the detection of urogenital tumors. Such means is more convenient for multiple, long-term and prognostic monitoring, and reduces the suffering of patients.
With comprehensive research and efforts, the inventors of the present application developed, for the first time, a method of screening classification markers by detecting copy number variations (CNVs) and methylation haplotype load (MHL) of DNA methylation haplotype blocks (MHBs) in urine sediment genomic DNAs, and further developed a method of diagnosing urogenital tumors with high sensitivity and specificity, which can not only well distinguish tumor patients from healthy people, but also localize urogenital tumors. In addition, a prognostic survival model and corresponding 9 bladder cancer prognostic markers and 16 renal cancer prognostic markers were constructed by integrating clinical prognostic data from bladder cancer and renal cancer. Therefore, the following inventions are provided.
One aspect of the present application relates to a DNA classification method, comprising
calculating the MHL value or β mean of a DNA methylation haplotype block of a sample of interest and/or calculating the DNA copy number variation data of the sample of interest; and
calculating the similarity between the MHL value or β mean of the DNA methylation haplotype block of the sample of interest DNA and the MHL value or β mean of a DNA methylation haplotype block of a respective classification label, and/or calculating the similarity between the DNA copy number variation data of the sample of interest and the DNA copy number variation data of the respective classification label; and
determining the classification for the DNA in the sample of interest by using a classifier model and based on the similarity.
Preferably, the β mean is obtained by 450K chip data or 850K chip data.
In one or more embodiments of the present application, in the DNA classification method, the MHL value of the DNA methylation haplotype block and the DNA copy number variation data of a sample of interest are calculated; and the similarity between the MHL value of the DNA methylation haplotype block of the sample of interest and the MHL value of the DNA methylation haplotype block of a respective classification label, and the similarity between the DNA copy number variation data of the sample of interest and the DNA copy number variation data of a respective classification label are calculated.
In one or more embodiments of the present application, in the DNA classification method, the MHL value of the DNA methylation haplotype block of a sample of interest is calculated; and the similarity between the MHL value of the DNA methylation haplotype block of the sample of interest and the MHL value of the DNA methylation haplotype block of a respective classification label is calculated.
In one or more embodiments of the present application, in the DNA classification method, a β mean of a DNA methylation haplotype block of a sample of interest is calculated; and the similarity between the β mean of the DNA methylation haplotype block of the sample of interest and the β mean of the DNA methylation haplotype block of a respective classification label is calculated.
In one or more embodiments of the present application, in the DNA classification method, determining the classification for the DNA in the sample of interest comprises
determining a correlation between the MHL value of the DNA methylation haplotype block of a respective classification label and a human urogenital tumor, and/or a correlation between the DNA copy number variation data of a respective classification label and a human urogenital tumor by using a random forest model and based on the similarity; and
determining the classification for the DNA in the sample of interest by using the classifier model and based on the correlation.
In one or more embodiments of the present application, in the DNA classification method, determining the correlation between the MHL value of the DNA methylation haplotype block of a respective classification label and a human urogenital tumor comprises, based on the correlation, ranking the MHL value of the DNA methylation haplotype block to form a vector sequence, and inputting the vector sequence into the random forest model to determine a correlation between the MHL value of the DNA methylation haplotype block and a human urogenital tumor;
and/or
determining the correlation between the DNA copy number variation data of a respective classification label and a human urogenital tumor comprises, based on the correlation, ranking the DNA copy number variation data to form a vector sequence, and inputting the vector sequence into the random forest model to determine a correlation between the DNA copy number variation data of the classification label and a human urogenital tumor.
In one or more embodiments of the present application, in the DNA classification method, the human urogenital tumor is any one, any two (prostate cancer and urothelial cancer, urothelial cancer and renal cancer, or prostate cancer and renal cancer), or all three selected from the group consisting of prostate cancer, urothelial cancer, and renal cancer;
preferably, the renal cancer is a kidney renal clear cell carcinoma,
preferably, the urothelial cancer is upper tract urothelial cancer and/or bladder cancer,
preferably, the prostate cancer is prostate adenocarcinoma; and
preferably, the human urogenital tumor is diagnosed by biopsy from a surgery.
In one or more embodiments of the present application, in the DNA classification method, the random forest model includes at least three random forest binary classifiers and is selected from any one, any two, any three or all four of the following groups I-VI:
I. normal-vs-renal cancer, normal-vs-urothelial cancer, and normal-vs-prostate cancer;
II. renal cancer-vs-normal, renal cancer-vs-urothelial cancer, and renal cancer-vs-prostate cancer;
III. urothelial cancer-vs-normal, urothelial cancer-vs-renal cancer, and urothelial cancer-vs-prostate cancer; and
IV. prostate cancer-vs-normal, prostate cancer-vs-renal cancer, and prostate cancer-vs-urothelial cancer.
In one or more embodiments of the present application, the DNA classification method comprises voting for each group, and determining the group with the highest number of votes as the final classification, wherein if equal numbers of votes occur, the category with the highest prediction probability among the groups with the equal number of votes is determined as the final classification.
Since it is theoretically impossible for a female to be predicted to have prostate cancer, if a female sample is predicted to be prostate cancer, a sub-optimal prediction result is taken. For example, if the vote predicted to be renal cancer is second only to prostate cancer, the predictive label of the female sample is defined as renal cancer. If equal numbers of votes occur in groups, the probabilities in the groups are compared. The category with higher probability is determined as the final prediction result of the female sample.
In one or more embodiments of the present application, in the DNA classification method, the sample is a urine sample, preferably urina sanguinis, and more preferably, urine sediment of the urina sanguinis. Urine sediment can be obtained via technical means known to a person skilled in the art, for example, by centrifuging a urine sample and removing the supernatant; and preferably, the centrifugation is performed at a temperature less than or equal to 4° C.
In one or more embodiments of the present application, in the DNA classification method, the MHL value of the DNA methylation haplotype block of the sample of interest, the MHL value of the DNA methylation haplotype block of a respective classification label, the DNA copy number variation data of the sample of interest, and the DNA copy number variation data in a respective classification label are all calculated from the sequencing data of the DNAs in the urine sample;
preferably, the DNAs in the urine sample are urine sediment DNAs; and
preferably, the sequencing data is whole genome methylation sequencing data, such as whole genome bisulfite sequencing (WGBS) data; and preferably, the sequencing depth is 1×-5×.
In one or more embodiments of the present application, in the DNA classification method, the DNA methylation haplotype block of the sample of interest is the same as the DNA methylation haplotype block of a respective classification label; and/or
the DNA copy number variation regions of the sample of interest are the same as the DNA copy number variation regions of a respective classification label;
preferably, the methylation haplotype blocks and the copy number variation regions are those as shown in any one, any two, any three, any four, any five or all six of Tables 1-6, or as shown in Table 11 and/or Table 12.
In one or more embodiments of the present application, in the DNA classification method, the MHL value of the DNA methylation haplotype block of the sample of interest and the MHL value of DNA methylation haplotype block of a respective classification label are calculated by using MONOD2 software, and/or the DNA copy number variation data of the sample of interest and the DNA copy number variation data of a respective classification label are calculated by using Varbin;
preferably, the MHL value corresponding to the respective methylation haplotype block in the WGBS data is calculated by using MONOD2 software, and/or the copy number variation data corresponding to the respective copy number variation region in the WGBS data is calculated by using Varbin, wherein the methylation haplotype block and the copy number variation region are those as shown in any one, any two, any three, any four, any five, or all six of Table 1-6, or as shown in Table 11 and/or Table 12.
In one or more embodiments of the present application, in the DNA classification method, the DNA copy number variation data of the sample of interest and/or the DNA copy number variation data of a respective classification label are calculated in the following way.
In one or more embodiments of the present application, in the DNA classification method, the biomarker is a DNA segment from a start position S±m to a termination position T±n on a chromosome;
wherein S is a start site, T is a termination site, and the start and termination sites are those as shown in any one, any two, any three, any four, any five, or all six of Tables 1-6, or the start and termination sites are those as shown in Table 11 and/or Table 12; and
wherein m and n are independently non-negative integers less than or equal to 6000.
In one or more embodiments of the present application, in the DNA classification method, m and n are independently 5000, 4000, 3000, 2000, 1500, 1000, 500, 300, 200, 150, 100, 90, 80, 70, 60, 50, 40, 30, 20, 10, 5, or 0.
Another aspect of the present application relates to a method for the detection, diagnosis, classification, risk assessment or prognostic assessment of a human urogenital tumor, comprising
(1) obtaining a urine sample and extracting urine sediment DNAs;
(2) fragmenting the DNAs into fragments of 300-500 bp;
(3) constructing a whole genome library, preferably a whole genome methylation sequencing library, such as a whole genome bisulfite sequencing library, using the obtained DNA fragments; and
(4) classifying the DNA fragments in the library using any DNA classification method described in the present application, wherein the DNA fragments serve as the DNA in the sample of interest.
In one or more embodiments of the present application, in the method for the detection, diagnosis, classification, risk assessment, or prognostic assessment of a human urogenital tumor, the urogenital tumor is one or more selected from the group consisting of prostate cancer, urothelial cancer, and renal cancer; and preferably, the renal cancer is kidney renal clear cell carcinoma, the urothelial cancer includes upper tract urothelial cancer and bladder cancer, and the prostate cancer is prostate adenocarcinoma.
In one or more embodiments of the present application, in the method for the detection, diagnosis, classification, risk assessment or prognostic assessment of a human urogenital tumor, in step (1), the urine sample is urina sanguinis; and preferably, the urine sample is urine sediment of the urina sanguinis.
In one or more embodiments of the present application, in the method for the detection, diagnosis, classification, risk assessment or prognostic assessment of a human urogenital tumor, in step (2), the DNAs are fragmented into fragments of 350-450 bp.
A further aspect of the present application relates to a device for the detection, diagnosis, classification, risk assessment or prognostic assessment of a human urogenital tumor, comprising:
I. ‘normal decision unit’:
normal-vs-renal cancer, normal-vs-urothelial cancer, and normal-vs-prostate cancer;
II. ‘renal cancer decision unit’:
renal cancer-vs-normal, renal cancer-vs-urothelial cancer, and renal cancer-vs-prostate cancer;
III. ‘urothelial cancer decision unit’:
urothelial cancer-vs-normal, urothelial cancer-vs-renal cancer, and urothelial cancer-vs-prostate cancer;
IV. ‘prostate cancer decision unit’:
prostate cancer-vs-normal, prostate cancer-vs-renal cancer, and prostate cancer-vs-urothelial cancer,
preferably, the decision units can perform any DNA classification method described in the present application.
A further aspect of the present application relates to a device for the detection, diagnosis, classification, risk assessment or prognostic assessment of a human urogenital tumor, comprising
a memory; and
a processor coupled to the memory;
wherein program instructions which can be executed by the processor are stored in the memory, and the program instructions include any one, any two, any three, or all four decision units selected from the group consisting of
I. ‘normal decision unit’:
normal-vs-renal cancer, normal-vs-urothelial cancer, and normal-vs-prostate cancer;
II. ‘renal cancer decision unit’:
renal cancer-vs-normal, renal cancer-vs-urothelial cancer, and renal cancer-vs-prostate cancer;
III. ‘urothelial cancer decision unit’:
urothelial cancer-vs-normal, urothelial cancer-vs-renal cancer, and urothelial cancer-vs-prostate cancer;
IV. ‘prostate cancer decision unit’:
prostate cancer-vs-normal, prostate cancer-vs-renal cancer, and prostate cancer-vs-urothelial cancer;
wherein each decision unit comprises three random forest binary classifiers.
In one or more embodiments of the present application, for the device, the processor is configured to perform any classification method described in the present application based on the instructions stored in the memory.
In one or more embodiments of the present application, for the device, the urogenital tumor is one or more selected from the group consisting of prostate cancer, urothelial cancer, and renal cancer;
preferably, the renal cancer is a kidney renal clear cell carcinoma,
preferably, the urothelial cancer is upper tract urothelial cancer and/or bladder cancer, and
preferably, the prostate cancer is prostate adenocarcinoma.
A further aspect of the present application relates to the use of any one of the following items 1) to 3) in the preparation of a medicament for the detection, diagnosis, risk assessment or prognosis assessment of a human urogenital tumor:
1) the biomarkers described in the present application (i.e., the methylation haplotype blocks and/or the copy number variation regions);
2) DNAs in human urine, in particular in the urine sediment of human urine;
preferably, the urine is urina sanguinis, and
preferably, the DNAs are 300-500 bp, such as 350-450 bp, in length;
3) A DNA library prepared from item 2); preferably, the DNA library is a whole genome library, preferably a whole genome methylated sequencing library such as a whole genome bisulfite sequencing library;
preferably, the urogenital tumor is one or more selected from the group consisting of prostate cancer, urothelial cancer, and renal cancer;
preferably, the renal cancer is a kidney renal clear cell carcinoma,
preferably, the urothelial cancer is upper tract urothelial cancer and/or bladder cancer, and
preferably, the prostate cancer is prostate adenocarcinoma.
The present application also relates to a set of biomarkers (i.e., the methylation haplotype blocks and/or the copy number variation regions), wherein a biomarker is a DNA segment from a start position S±m to a termination position T±n on a chromosome;
wherein S is a start site, T is a termination site, and the start and termination sites are those as shown in any one, any two, any three, any four, any five, or all six of Tables 1-6, or the start and termination sites are those as shown in Table 11 and/or Table 12; and
wherein m and n are independently non-negative integers less than or equal to 6000.
In one or more embodiments of the present application, for the biomarkers, m and n are independently 5000, 4000, 3000, 2000, 1500, 1000, 500, 300, 200, 150, 100, 90, 80, 70, 60, 50, 40, 30, 20, 10, 5, or 0.
Some terms involved in the present application are explained below.
The term “bin” (section/region) is a generic description about artificially defining or dividing a genome by a certain length in the field of genomics. For example, if the human genome of about 3 billion base pairs is divided into 3000 bins on average, the size of each bin is about one million base pairs.
The term “coverage” refers to the proportion of a region of the genome that has been detected at least once accounting for the entire genome. Coverage is a term used to measure the extent to which the genome is covered by data. Due to the presence of complex structures (such as high GC and repeat sequences) in the genome, the final sequence obtained by sequencing, splicing and assembling often cannot cover all regions, and the regions which cannot be obtained are referred to as Gap. For example, when a bacterial genome is sequenced, and the coverage is 98%, 2% of the sequence region is not obtained by sequencing.
The term “sequencing depth” refers to the ratio of the total number of bases (bp) obtained by sequencing to the size of the genome, or it is understood as the average number of times that each base in the genome is sequenced. For example, assuming that the size of a gene is 2M and the obtained total amount of data is 20M, the sequencing depth is 20M/2M=10×.
The term “reads” or “read” refers to a read fragment, i.e., a read sequence.
The term “pair-end reads” refers to paired reads.
The term “copy number variations (CNVs)” refers to a deletion or duplication of a relatively large DNA fragment, typically an increase or a decrease in the copy number of DNA fragments of hundreds of bp to millions of bp. CNVs are caused by genomic rearrangements and are one of the important pathogenic factors of tumors. In one embodiment of the present application, the copy number variation is calculated in the following way.
The genome of a test sample is divided into 5,000-500,000 bins (e.g., 50,000 bins) of equal length or the same theoretical simulated copy number. The ratio A/B of the read number corresponding to each bin is calculated by software or algorithms such as Varbin, CNVnator, ReadDepth or SegSeq (A is the number of actual reads corrected for the GC content in a bin; B is the number of theoretical reads in the bin, which is obtained by dividing the total number of reads read in the sample by the total number of bins). The ratio A/B is the copy number variation.
The term “theoretical simulated copy number” involves dividing a genome into several regions of equal or unequal length by a software and/or method of calculating copy number, where theoretical copy number contained in each region is same by data simulation.
The term “MHB” refers to DNA methylation haplotype blocks, also referred to herein as DNA methylation haplotype region or DNA methylation haplotype modules, meaning a linkage region in which DNA co-methylation frequently occurs in the genome. The basic principle is based on the co-methylation linkage of adjacent CpG sites. The algorithm extends the concept of linkage disequilibrium (LD) in traditional genetics, which indicates the degree of co-methylation of adjacent CpG sites in DNA methylation, that is, the linkage condition of DNA methylation. The linkage condition of adjacent CpG sites is first calculated by DNA methylation haplotype, and the region with r2 not less than 0.5 in adjacent CpG sites is further defined as potential MHBs. The potential MHBs are then expanded according to the overlapping CpG sites in the MHB region, and final MHBs are obtained. They can be identified by using technical means known to a person skilled in the art, for example, by using MONOD2 software (http://genome-tech.ucsd.edu/public/MONOD_NG_TR44413/scripts_and_codes/) developed by Kun Zhang's Research Team.
The term “MHL” refers to DNA methylation haplotype load, which represents the heterogeneous distribution of different DNA methylation haplotypes in a given region, i.e., the proportion of CpG site methylation modifications.
The term “TNM” represents a tumor staging system in which:
“T” is the initial letter of the wording “tumor”, and refers to the size or direct extent of a primary tumor. With an increase in tumor volume and an increase in the extent of adjacent tissue involvement, it is represented by T1˜T4 in turn.
“N” is the initial letter of the wording “Node”, and refers to the involvement of regional lymph node. When the lymph node is not involved, it is represented by N0. With an increase of the degree and extent of lymph node involvement, it is represented by N1˜N3 in turn.
“M” is the initial letter of the wording “metastasis” and refers to distant metastasis (usually hematogenous metastasis). No distant metastasis is represented by M0 and the presence of distant metastasis is represented by M1. On this basis, a specific stage is delineated by the grouping of the three indicators of TNM.
One or more of the following technical effects are achieved in the present application.
(1) Non-invasive diagnosis in the true sense. Sampling is simple, which only requires obtaining a certain volume of urina sanguinis, and there is no trauma to the subjects. This is advantageous for sample collection, diagnosis, long-term monitoring and regular monitoring of prognosis.
(2) High success rate of library construction. The amount of urine sediment DNAs is much more than that of urine cell-free DNAs, so that the amount of starting DNAs for library construction is much more than that of cfDNAs for library construction. In addition, there are kits available for library construction and sequencing, which makes the operation easier and more stable and reliable.
(3) Low-depth high-throughput sequencing. In the present application, the integration of the information of DNA methylation and DNA copy number variation and the extraction of a tumor signal in a unit of a region by optimizing a modeling algorithm can not only maximumly retain the tumor signal, but also maximumly reduce sequencing cost. Theoretically, it is possible to obtain a result with high sensitivity and specificity at a sequencing depth of about 1× to 5×.
(4) High-accuracy diagnosis of a single tumor. The diagnosis and recurrence monitoring of common tumors of the urinary system (such as renal cancer, bladder cancer and prostate cancer) can be achieved using the constructed binary classifier model.
(5) Tumor localization. The use of the multi-stage classification system of the present application can not only determine whether a tumor is present or not, but also locate the potential tumor type of a tumor patient.
(6) Potential application in prognostic risk assessment. The prognostic markers screened by the present application can be potentially applied to the survival prognostic assay in a tumor patient.
The embodiments of the present application will be described in detail below in reference to Examples. It should be understood by a person skilled in the art that the following Examples are merely illustrative of the present application and are not intended to limit the scope of the present application. The experimental methods without specifying their protocols in the Examples are generally carried out according to conventional protocols, or according to protocols recommended by manufacturers. The used reagents or the instruments without specifying the manufacturer are commercially available conventional products.
In the present application,
The 450K chip data refers to the Illumina Infiium Human Methylation 450 BeadChip chip technology developed by Illumina, where 450K refers to the number of probes on the chip, which can detect the corresponding number of methylation sites.
The 850K chip data refers to the Illumina Infiium Human Methylation 850 BeadChip chip technology developed by Illumina, where 850K refers to the number of probes on the chip, which can detect the corresponding number of methylation sites.
The TCGA snp6.0 chip data is provided by a public database, which can be downloaded, for example, from http://firebrowse.org/?cohort=PRA or https://portal.gdc.cancer.gov/. The number of copy number variations in the area covered by the SNP6.0 chip can be detected.
The available clinical data of the TCGA is provided by a platform for tumor research, which is provided by the TCGA official website (https://www.cancer.gov/). A person skilled in the art can also obtain the available clinical data of the TCGA by other integration software and online platforms, such as http://firebrowse.org/and software such as TCGA download widgets.
1. Subject Population
Urine samples from a total of 313 subjects were collected, as shown in
2. Experimental Methods
(1) Fresh urine (urina sanguinis) from preoperative tumor patients and fresh urine (urina sanguinis) from healthy people were collected. The urines were collected in 50 ml centrifuge tubes with a volume of about 45-50 ml per urine sample.
(2) The collected urina sanguinis samples were centrifuged at 3500 rpm and 4° C. for 10 min, respectively. The supernatants were removed to obtain urine sediments.
(3) The urine sediments were washed twice with PBS buffer (500 ml of PBS buffer was added each time, and after centrifugation at 13000 g for 1 min, the supernatants were removed), and then the urine sediments were transferred to 1.5 ml EP tubes.
(3) Urine sediment genomic DNAs (urine sediment gDNAs) were extracted by using QIAamp DNA Mini Kit. After extraction, the concentration of the DNAs was measured with Qubit and the DNAs were stored at −80° C. for later use.
313 DNA samples were prepared.
50-200 ng of the DNA samples obtained in Example 1 were taken, respectively, as the start DNAs for library construction and lambda DNAs (all CpG sites included unmethylated C) and 5 mC DNAs (all CpG sites included methylated C) were added in a ratio of 3:1000. The DNAs were then fragmented with a Covaris sonicator such that the major length peaks of the fragments were in a range of 400 bp. The fragmented DNAs were then end repaired with NEBNext Ultra II End Repair/dA-Tailing Module 96 rxns (Cat. No. E7546) and were polyadenylated (polyA). Then, methylation PE linkers were added by using NEBNext Ultra II Ligation Module, 96 rxns unit (Cat. No. E7595L).
The resulting water-soluble DNAs with linkers ligated (i.e., the library) were subjected to a bisulfite treatment by using a EZ DNA methyhlation Gold kit (Zymo Research). The specific procedures were performed in accordance with the instructions for use of the kit. Afterwards, the DNAs were purified, amplified by PCR, and the concentration of the DNAs was determined by using the nucleic acid and protein quantitative analyzer Qubit2.0 of Life Tech, obtaining a DNA library.
The resulting DNA library was sent to Novogene for quality control of library fragmentation and concentration using Agilent 2100 and AB17500 Fluorescent quantitative PCR instruments, respectively. There was no problem in library examination, thereby obtaining a BS-seq library of 313 urine sediment gDNA samples for subsequent library sequencing.
1. Test Samples:
The BS-seq library of 313 urine sediment gDNAs prepared in the above Example 2.
2. Experimental Methods
Novogene sequencing company was entrusted to perform whole-genome sequencing on the BS-seq library of 313 urine sediment gDNAs.
3. Experimental Results
The data (i.e., a fastq raw file) on 150 bp pair-end reads of the BS-seq library of 313 urine sediment gDNAs was obtained for subsequent data preprocessing and tumor marker analysis.
The reads of the BS-seq library of 313 urine sediment gDNAs obtained by sequencing in Example 3 was first subjected to quality control by Trimmomatic (version: Trimmomatic-0.32), including removal of low-quality reads and linkers. Next, genomic alignment was performed using Bismark (version: bismark v0.14.5) alignment software and PCR repeat amplification reads (deduplication) were removed. Then, the overlap regions between reads were then removed using bamUtil (version: bamUtil_1.0.12) software. The resulting bam file was then used as a starting file for an analysis of DNA copy number and methylation. Finally, the output data coverage of each sample in the BS-seq library of 313 urine sediment gDNAs was approximately 1×-5×.
For the DNA methylation feature selection (shown in
(1) Calculation of Methylation Frequency (average methylation level): for a given region, if the number of reads covering the base C was defined as Nc and the number of reads covering the base T was defined as Nt, the methylation level of the region was Nc/(Nc+Nt).
Reference: Chen, K. et al. Loss of 5-hydroxymethylcytosine is linked to gene body hypermethylation in renal cancer. Cell Research. 26(1):103-118 (2016).
(2) Calculation of Methylation Entropy (ME):
wherein b denotes the number of corresponding CpG in a given region, n denotes the number of methylation haplotypes in a given region, and P (Hi) denotes the probability of observing a methylation haplotype in a given region.
Reference: Xie, H. et al. Genome-wide quantitative assessment of variation in DNA methylation patterns. Nucleic Acids Res. 39, 4099-4108 (2011).
(3) Calculation of Epi-polymorphism:
The probability of occurrence of methylation haplotype i for a given region was Pi, and the number of methylation haplotypes was n.
Reference: Landan, G. et al. Epigenetic polymorphism and the stochastic formation of differentially methylated regions in normal and cancerous tissues. Nat. Genet. 44, 1207-1214 (2012).
(4) Calculation of Methylation Haplotypes
For a given region, the methylation status of the corresponding CpG covering reads was the methylation haplotype.
Reference: Shoemaker, R., Deng, J., Wang, W. & Zhang, K. Allele-specific methylation is prevalent and is contributed by CpG-SNPs in the human genome. Genome Res. 20, 883-889 (2010).
Where an MHL value cannot be calculated for an MHB because the sequenced reads did not cover the MHB, the MHL value of the MHB was filled with the average MHL value of the sample itself. The average MHL value was calculated as follows.
For each sample, there were 147888 MHBs to calculate MHLs. The MHBs where MHLs cannot be calculated were NA, and the corresponding number was n(NA). The MHL values were calculated if the MHBs of the MHLs can be calculated. The corresponding number was 147888-n(NA). The sum of all MHLs of the corresponding MHBs for which MHL values can be calculated the was Sum, and the average MHL value for each sample was Sum/(147888-n(NA)).
Finally, almost 150,000 MHBs containing MHL values can be obtained for each sample. These MHBs were used as initial candidate features for DNA methylation analysis. In order to narrow the range of screening features, the inventors divided the features into two groups.
One group was candidate raw F1, representing that the MHL values of some MHBs were different for the urine sediment gDNAs not only between the tumor patients and healthy people (student t-test, p value<0.05) (the difference analysis can use statistical analysis languages such as limma R package, student t-test test, and filter features by limiting the p-value threshold; or statistical analysis software such as SPASS, SAS, Metalab or Origin; similarly hereinafter), but also between the solid tumor tissues and the corresponding pericarcinomatous tissues in the TCGA methylation 450 K data (student t-test, p value<0.05).
The other group was candidate raw F2, representing that the MHL values of some MHBs were different for the urine sediment gDNAs not only between the tumor patients and healthy people (student t-test, p value<0.05), but also between the solid tumor tissue and the corresponding pericarcinomatous tissue in the constructed Whole Genome Bisulfite Sequencing (WGBS) data (student t-test, p value<0.05).
Next, MHBs were gradually kicked out for raw F1 and raw F2, respectively, until the accuracy (obtained by 10-fold cross-validation) and the kappa coefficient (the kappa coefficient was used for consistency test, and can also be used to measure classification accuracy, which was calculated based on a hybrid matrix) of the corresponding random forest model no longer increased. At this time, the obtained MHBs corresponded to F1 and F2 (as shown in
In order to verify the reliability of the feature selection, the verification was performed by the inventors in combination with the TCGA methylation 450 K data. The verification method was as follows.
Firstly, using the screened F1 features, a β mean value of the F1 feature region corresponding to each sample was preliminarily calculated based on the TCGA 450K data (for a given region, if the number of 450K probes was n, and the sum of β values of all probes in the corresponding region was Sum β, then the average β value of the corresponding region was Sum_β/n), and then a hybrid matrix was constructed. Next, the samples were divided into a training set and a test set according to a ratio of 2:1. Then, the training set was modeled by a random forest algorithm, and the test set was used to test the predictive sensitivity and specificity of the model. Finally, the predictive performance of the model was displayed by combining the ROC curve.
The results showed that the selected feature could well distinguish a cancerous tissue from the corresponding pericarcinomatous tissue (as shown in
For the screening of subsequent feature of CNVs (F4) (as shown in
Similar to the F1 feature validation in Example 5, the inventors verified the F4 features using TCGA snp6.0 chip data. The results showed that the F4 features could well distinguish cancerous tissues from corresponding pericarcinomatous tissues (as shown in
In order to further improve the model performance, the F3 features and the F4 features were integrated with reference to the method in Example 6. The candidate features were gradually kicked out until the accuracy and the kappa value of the model prediction no longer increased, at which time the remaining features were used as F5, as shown in Tables 1 to 6 below, where the importance was a result of output with importance parameters after the model was built using randomForest R package.
F5 represented the features required for a hybrid model for integrating DNA methylation and copy number information, and the classification model constructed with F5 performs the best. In this way, the binary classification model was established.
This model can be used to distinguish tumor patients from healthy people.
As previously described, the inventors collected 100 samples of urothelial cancer (UC) (including bladder cancer and upper tract urothelial cancer), 65 samples of kidney renal clear cell carcinoma (KIRC) and 60 samples of prostate cancer (PRAD), and 88 samples of healthy people. Each sample included the feature information of F1 to F5. Taking the UC-vs-Healthy binary classifier as an example, the samples were first randomly rearranged so that the composite matrix of the samples had no preference, and then was split into a training set and a test set according to a ration of 5:1. Next, modeling was performed using the above-screened features (e.g., F5) combined with a support vector machine algorithm. Then, the test set was used to test the model performance, including accuracy, sensitivity, specificity, AUC and Kappa value. The above process was repeated 10 times, and the average accuracy, sensitivity, specificity, area under the curve (AUC) and Kappa coefficient of the ten results represented the stable classification performance of a binary classifier of urothelial cancer-vs-healthy. Other binary classifiers (Renal Cancer-vs-Healthy, Prostate Cancer-vs-Healthy) were constructed in a similar way.
The results were shown in Table 7 below.
The results showed that the accuracy of the 10-time repeated modeling and prediction of the corresponding classifier model was more than 90%. By feature selection and construction of the corresponding binary classifiers, the classifier model constructed by the inventors using the F5 features had the best performance, not only higher than the performance of the classifiers constructed only with DNA methylation information (F1, F2 and F3), but also higher than the performance of the classifier constructed with only DNA copy number information (F4).
For the tumor tissue typing model, the inventors constructed a multi-stage classification model (named as genitourinary cancers seek, abbreviated as GUseek) based on binary classifier models (shown in
The main aim of GUseek was to differentiate urothelial cancer (UC) (including bladder cancer and upper tract urothelial cancer), kidney renal clear cell carcinoma (KIRC), and prostate cancer (PRAD).
Based on the binary classification concept, there were six sets of binary classifiers, i.e., urothelial cancer-vs-healthy, urothelial cancer-vs-renal cancer, urothelial cancer-vs-prostate cancer, renal cancer-vs-healthy, renal cancer-vs-prostate cancer, and prostate cancer-vs-healthy, which can be combined into four sets of classification decision systems, i.e.:
a urothelial cancer decision system (including urothelial cancer-vs-healthy, urothelial cancer-vs-renal cancer and urothelial cancer-vs-prostate cancer),
a renal cancer decision system (including urothelial cancer-vs-renal cancer, renal cancer-vs-healthy and renal cancer-vs-prostate cancer),
a prostate cancer decision system (including urothelial cancer-vs-prostate cancer, renal cancer-vs-prostate cancer and prostate cancer-vs-healthy), and
a healthiness decision system (including urothelial cancer-vs-healthy, renal cancer-vs-healthy and prostate cancer-vs-healthy).
An unknown sample was first mapped to each decision system for predictive analysis, and the proportion of the prediction category of each decision system was provided accordingly. By integrating the scores of various types in the four decision systems, the category with the highest score was defined as the prediction category of the unknown sample. If there was more than one category with the highest score, the category with the highest score probability was selected as the final prediction category for the unknown sample. Considering that it was theoretically impossible for a female to be predicted to have prostate cancer, if a female sample was predicted to be prostate cancer, a sub-optimal prediction result was taken. For example, if the vote predicted to be renal cancer was second only to prostate cancer, the predictive label of the female sample was defined as renal cancer. If the numbers of votes were the same, then the probabilities were compared. The category with higher probability was taken as the final prediction result of the female sample.
The GUseek model can use the advantages of binary classification to the maximum, while a more powerful multi-stage classifier can be constructed by integrating multiple machine learning algorithms. By integrating the SVM algorithm, the GUseek constructed by the inventors can achieve 10-time repeated modeling and prediction accuracy up to nearly 90% (89.43%). The specific method was as follows.
The present inventors first randomly rearranged the collected 100 samples of urothelial cancer (UC) (including bladder cancer and upper tract urothelial cancer), 65 samples of kidney renal clear cell carcinoma (KIRC) and 60 samples of prostate cancer (PRAD), and 88 samples of healthy people and split the samples into a training set and a test set according to a ratio of 5:1 (see Table 8).
Six sets of binary classifiers were then constructed according to the above method of constructing binary classifiers, and were further combined to form four decision systems. For each sample in the test set, prediction was first performed in the binary classifiers and corresponding prediction categories and probabilities were obtained according to the input requirements of the binary classifiers of individual decision systems. The category of the predicted sample was determined by comparing the predicted times (the number of votes) of the sample by individual decision systems. If the numbers of votes for determining the decision category were comparable, the corresponding probabilities were further compared, and the category with the highest probability was taken as the final prediction category of the sample. In this way, the inventors can finally obtain the prediction classification of each test set sample, and can further obtain the prediction overall accuracy and Kappa coefficient of the GUseek model by constructing a hybrid matrix. The above process was repeated 10 times, and the obtained average accuracy was the stability performance of the GUseek. See
Using the integration algorithm GUseek proposed by the inventors, GUseek showed very high accuracies in 10-time remodeling and predictions (10-time average reached 89.43%, see
First, the training set that had been split according to a ratio of 5:1 by the GUseek analysis process was modeled according to the above algorithm in sequence, and then model evaluation was performed by using the test set. The assessment result was demonstrated by a hybrid matrix. The comparison results of one random time were shown in Tables 9-10, and the ten-time average accuracy was shown in
The algorithm developed by the present inventors can integrate the optimal conventional algorithm to achieve the optimal combination, i.e., each decision classification system, and can be constructed by selecting an algorithm with the best classification effect, which then can be combined into an overall optimal classification system.
Prognostic markers of bladder cancer and renal cancer were screened respectively by using available clinical data of TCGA. The specific steps were as follows.
Firstly, a statistical test was used to find the MHBs that can not only distinguish the tumor tissue from the corresponding pericarcinomatous tissue in the available clinical data of TCGA, but also distinguish the aforementioned 313 tumor patients from the healthy people in the urine sediment gDNAs. The specific procedure was shown in
These regions were then subjected to univariate and multivariate cox regression analysis. A statistically significant MHBs were selected for LASSO cox prognostic risk assessment to determine high-risk and low-risk groups and a combination of optimal prognostic risk features (resulting in a prognostic risk assessment model). The random forest algorithm was further used for these features, and the features were gradually kicked out until the accuracy of the prognostic model no longer increased. The MHBs (9 MHBs for the prognosis of bladder cancer and 16 MHBs for the prognosis of renal cancer) closely related to the prognosis of bladder cancer and renal cancer were finally found, which can potentially be applied to prognostic survival analysis of tumor patients.
The R packages used in the selection of model features include survival, survminer, glmnet and glmSparseNet. After the features for constructing a model were selected, there were many relevant R packages in R that can be used to analyze ROC curve and K-mean survival. For example, in the Example, the R package used in constructing the ROC curve was ROCR and the R package used in analyzing the K-mean survival was glmSparseNet.
The markers for bladder cancer and renal cancer prognosis were shown in Tables 11 and 12 below.
The AUC value of the ROC curve of the prognostic survival model constructed by the present inventors was very high (
The above experimental results showed that the present inventors have developed, for the first time, a model for the diagnosis, localization and prognosis of urogenital tumors that integrates the methylation haplotype and copy number information of urine sediment genomic DNAs. The model can be used to not only predict with high accuracy whether an unknown sample is a tumor or healthy, but also determine the tissue origin of the tumor if the sample is a tumor. By comparing the multivariate classifier algorithms, the GUseek system constructed by the inventors is significantly superior to other commonly used machine algorithm models, including SVM, LASSO, LDA, knn, RandomForest, and Bayes algorithms (
On the first day, the test subjects were enrolled, and a 50 ml of urina sanguinis collection tube was distributed to each subject. The test subjects were then required to collect 50 ml of urina sanguinis in the following morning and send it to the urine collection site of the clinic. The urine was then centrifuged to obtain the corresponding urine sediment. Next, the urine sediment DNAs were extracted and a WGBS library was constructed and sequenced to obtain data information of the F5 features in WGBS. For example, MHL values corresponding to the F5 features in WGBS were calculated using MONOD2 software, and copy number variation data corresponding to the F5 features in WGBS were calculated by using Varbin. The basic protocols can follow those in the above Examples 1-4 and Example 7.
The acquired data information of the F5 features in WGBS was then imported into the classifier model constructed according to Example 7 or 8 of the present application. The model can output a possible category of an unknown subject, such as healthy or unhealthy, in particular which type of tumor it is where the subject is unhealthy. If a patient has developed a tumor and undergone surgery, testing at this time was similar to regular follow-up of the patient after surgery.
The prognosis model is only for tumor patients. The tumor patients with good prognosis and survival are expressed as a low-risk group, and the tumor patients with poor prognosis and survival are expressed as a high-risk group. The purpose of the prognostic model of the present application is to divide the high-risk and low-risk groups of patients.
On the first day, the test patients with renal or bladder cancer were enrolled, and a 50 ml of urina sanguinis collection tube was distributed to each patient. The test subjects were then required to collect 50 ml of urina sanguinis in the following morning and send it to the urine collection site of the clinic. The urine was then centrifuged to obtain the corresponding urine sediment. Next, the urine sediment DNAs were extracted and sent to a company to measure the 450 K or 850 K chip data of the sample. The data information of the prognostic marker characteristics in Table 11 and/or Table 12 in the 450 K or 850 K chip data was then obtained, such as the corresponding β mean (the mean of probe signals, which is positively correlated with the methylation level) of the prognostic markers in Table 11 and/or Table 12 in the 450 K or 850 K chip data. The acquired data information of the feature candidate prognostic markers in the 450 K or 850 K chip was then imported into the prognostic risk assessment model constructed in Example 9 of the present application. The model can output a possible category of a patient with unknown risk category, such as a high-risk group or a low-risk group. If a patient has developed a tumor and undergone surgery, testing at this time was similar to regular follow-up of the patient after surgery.
Although specific embodiments of the present application have been described in detail, a person skilled in the art will appreciate that various modifications and substitutions can be made to those details from the teachings of the disclosure, all of which are within the scope of the present application. The full scope of the present application is covered by the appended claims and any equivalents thereof.
| Number | Date | Country | Kind |
|---|---|---|---|
| 201911088433.6 | Nov 2019 | CN | national |
| Filing Document | Filing Date | Country | Kind |
|---|---|---|---|
| PCT/CN2020/122821 | 10/22/2020 | WO |
