Method and kit for analyzing susceptible site of smoking-related chronic obstructive pulmonary disease

Abstract

According to the analysis method and the kit for the susceptible site of the smoking-related chronic obstructive pulmonary disease, an initial comparison result is obtained through sample collection, library building and sequencing, data quality control and data comparison, and after the physical position of reads in a genome is obtained, SNP and genotyping are identified by using GATK4, and the susceptible site of the smoking-related chronic obstructive pulmonary disease is identified by using GATK4. The method comprises the following steps: carrying out GWAS analysis and PRS construction in combination with smoking factors by using Plink software, combining SNP data and using a logistic regression model to analyze a BED file to obtain GWAS overall data, and carrying out PRS construction based on a P + T method on the basis of taking the obtained GWAS overall data as baseline data and taking another group of independent samples as a data target, so that the method can be used for constructing the PRS for Chinese population through a bioinformatics method. The smoking-related chronic obstructive pulmonary disease genetic susceptible persons are accurately distinguished, so that potential high-risk people can prevent the chronic obstructive pulmonary disease in advance.

Description

Technical Field

[0001] The present invention relates to the technical field of bioinformatics, and in particular to a method for analyzing smoking-related chronic obstructive pulmonary disease susceptibility sites, and a kit utilizing the method for analyzing smoking-related chronic obstructive pulmonary disease susceptibility sites. Background Art

[0002] Chronic obstructive pulmonary disease (COPD) is a common respiratory disease characterized by persistent respiratory symptoms and airflow limitation, which severely impacts workability and quality of life. The China Adult Lung Health Study shows that there are approximately 99.9 million COPD patients in my country, with a prevalence of 13.7% among people aged 40 and above. COPD has become a chronic disease on par with hypertension and diabetes.

[0003] The pathogenesis of COPD is relatively complex and is currently classified as a polygenic disease. Its occurrence is affected by genetic factors, environmental factors (such as smoking) and the interaction between the two. Smoking is internationally recognized as the most important risk factor for COPD. The "China Smoking Health Report 2020" clearly pointed out that there is sufficient evidence that smoking can cause COPD, and even secondhand smoke exposure can cause a significant increase in respiratory symptoms and the incidence of COPD. Tobacco smoke can damage the lungs in various ways, such as affecting the defense function of the respiratory system and oxidative stress, resulting in a continuous decline in lung function. In addition, genetic causes play an important role in the occurrence and development of COPD. Current studies have confirmed that there is significant overlap between the risk loci for COPD and the population-based lung function-related loci, and that COPD-related gene loci are more abundant in areas involved in lung development.

[0004] In recent years, international research has been steadily developing using genome-wide association studies (GWAS) to identify genetic loci associated with COPD. To date, research using genetic risk scores to assess COPD risk is in its infancy. Summary of the Invention

[0005] In order to overcome the defects of the existing technology, the technical problem to be solved by the present invention is to provide an analysis method for smoking-related COPD susceptibility loci, which can accurately distinguish those who are genetically susceptible to smoking-related COPD in the Chinese population through bioinformatics methods, so that potential high-risk groups can prevent COPD early.

[0006] The technical solution of the present invention is: a method for analyzing smoking-related chronic obstructive pulmonary disease susceptibility sites, which comprises the following steps:

[0007] (1) Sample collection: Samples were collected from patients with smoking-related COPD and healthy individuals;

[0008] (2) Library construction and sequencing: DNA extraction, library construction, and sequencing were performed on all samples. The library construction used the Illumina whole genome amplification kit, the sequencing platform was Illumina's NovaSeq, and the target depth was 10x;

[0009] (3) Data quality control: Use the bbduk.sh script in the Trimmomatic and BBTools components to obtain high-quality valid data;

[0010] (4) Mutation identification: Perform data alignment to obtain the initial alignment results; after obtaining the physical location of the reads in the genome, use GATK4 to identify SNPs and genotype:

[0011] (4.1) Use the HaplotypeCaller module of GATK4 with default parameters to obtain the original mutation information;

[0012] (4.2) Use the SelectVariants module of GATK4 to separate SNVS and INDELs;

[0013] (4.3) Use the VariantFiltration module of GATK4 to filter low-quality mutation sites and obtain a vcf file containing accurate SNP information;

[0014] (4.4) Use the GenotypeGVCFs module of GATK4 to perform SNP annotation based on the dbSNP138 database, and label each SNP with an RS number;

[0015] (4.5) Use PLINK v1.90b6.21 64-bit to perform site and sample quality control on the obtained SNP vcf data;

[0016] (5) Plink was used to perform GWAS analysis and construct PRS in combination with smoking factors: Plink software was used in combination with SNP data, and the BED file was analyzed using a logistic regression model to obtain the overall GWAS data. The obtained overall GWAS data was used as baseline data, and another group of independent samples was used as data targets. PRS was constructed based on the P+T method.

[0017] The present invention collects samples, constructs libraries for sequencing, and performs data quality control to perform data comparison to obtain the initial comparison results. After obtaining the physical location of reads in the genome, GATK4 is used to identify SNPs and genotypes, and Plink is used to perform GWAS analysis and construct PRS in combination with smoking factors: Plink software is used in combination with SNP data, and the BED file is analyzed using a logistic regression model to obtain the overall GWAS data. Based on the obtained GWAS overall data as baseline data, another set of independent samples is used as data targets, and PRS is constructed based on the P+T method. Therefore, it is possible to accurately distinguish those who are genetically susceptible to smoking-related COPD in the Chinese population through bioinformatics methods, so that potentially high-risk groups can prevent COPD early.

[0018] A kit for analyzing smoking-related COPD susceptibility loci is also provided. This kit is used for mining smoking-related COPD susceptibility loci. Genomic DNA is used as the material, DNA fragmentation and end-repair and A addition are performed, and adapters with single-molecule barcode sequences are ligated to both ends using DNA ligase. A library is obtained after amplification. The library is hybridized with biotin-labeled probes in liquid phase and captured and enriched using streptavidin-coated magnetic beads. Finally, primers with tag sequences are used for amplification to obtain a captured library. Sequencing data is obtained through high-throughput sequencing, and the final gene variation information is obtained after analysis using the corresponding data analysis software.

[0019] The components of Part A of the kit are as follows:

[0020] Reagent 1: The component name is DNA fragmentation and repair enzyme, and its ingredients are DNA fragmentation enzyme and repair enzyme. The specification and quantity are 96 μL × 1 tube;

[0021] Reagent 2: The component name is interruption repair solution, which includes Tris, and the specification and quantity are 96 μL × 1 tube;

[0022] Reagent 3: The component name is ligation buffer, which includes Tris and ATP. The specification and quantity are 192 μL × 1 tube;

[0023] Reagent 4: The component name is ligase, the ingredients include ligase, the specification and quantity are 96 μL × 1 tube;

[0024] Reagent 5: The component name is adapter, the ingredients include oligonucleotides, the specifications and quantity are 3 μL × 24 tubes;

[0025] Reagent 6: The component name is primer, the ingredients include oligonucleotides, the specification and quantity are 70 μL × 1 tube;

[0026] Reagent 7: The component name is amplification solution, which includes Tris, dNTPs, and DNA polymerase. The specification and quantity are 1390 μL × 1 tube;

[0027] Reagent 8: The component name is blocking solution, which includes oligonucleotides, and the specification and quantity are 88 μL × 1 tube;

[0028] Reagent 9: The component name is hybridization solution, which includes Tween and dextran sulfate. The specification and quantity are 336 μL × 1 tube;

[0029] Reagent 10: Hybridization Enhancer, containing formamide, 150 μL per tube.

[0030] Reagent 11: The component name is capture probe, including oligonucleotides, and the specification and quantity are 150 μL × 1 tube;

[0031] Reagent 12: The component name is washing solution S, which includes ethylenediaminetetraacetic acid disodium salt, and the specification and quantity are 576 μL × 1 tube;

[0032] Reagent 13: The component name is washing solution 1, which includes sodium dodecyl sulfate, and the specification and quantity are 480 μL × 1 tube;

[0033] Reagent 14: The component name is washing solution 2, which includes sodium dodecyl sulfate, and the specification and quantity are 576 μL × 1 tube;

[0034] Reagent 15: The component name is washing solution 3, which includes sodium dodecyl sulfate, and the specification and quantity are 576 μL × 1 tube;

[0035] Reagent 16: Magnetic bead cleaning solution, containing Tris, EDTA, and sodium chloride. Quantity: 720 μL x 4 tubes.

[0036] Reagent 17: The component name is label, the ingredients include oligonucleotides, the specifications and quantity are 4 μL × 24 tubes;

[0037] Reagent 18: The component name is positive control, the ingredients include deoxyribonucleic acid, the specification and quantity are 24μL×1 tube;

[0038] Reagent 19: The component name is negative control, the ingredients include deoxyribonucleic acid, the specification and quantity are 24μL×1 tube;

[0039] The components of Part B of the kit are as follows:

[0040] Reagent 20: The component name is capture magnetic beads, and its ingredients include streptavidin magnetic beads. The specification and quantity are 880 μL × 1 tube. BRIEF DESCRIPTION OF THE DRAWINGS

[0041] Figure 1 The figure is a flow chart of the method for analyzing smoking-related COPD susceptibility sites according to the present invention.

[0042] Figure 2 The number of significant SNPs obtained by screening with different P values ​​is shown.

[0043] Figure 3 Shows a SNP with a P value less than 5×10 -7 Then the Metascape enrichment analysis results of non-smoking SNPs were removed.

[0044] Figure 4 The ROC curves of the prior art and the present invention are shown.

[0045] Figure 5 Shows the correlation between smoking and COPD. DETAILED DESCRIPTION

[0046] like Figure 1 As shown, a method for analyzing smoking-related COPD susceptibility loci comprises the following steps:

[0047] (1) Sample collection: Samples were collected from patients with smoking-related COPD and healthy individuals;

[0048] (2) Library construction and sequencing: DNA extraction, library construction, and sequencing were performed on all samples. The library construction used the Illumina whole genome amplification kit, the sequencing platform was Illumina's NovaSeq, and the target depth was 10x;

[0049] (3) Data quality control: Use the bbduk.sh script in the Trimmomatic and BBTools components to obtain high-quality valid data;

[0050] (4) Mutation identification: Perform data alignment to obtain the initial alignment results; after obtaining the physical location of the reads in the genome, use GATK4 to identify SNPs and genotype:

[0051] (4.1) Use the HaplotypeCaller module of GATK4 with default parameters to obtain the original mutation information;

[0052] (4.2) Use the SelectVariants module of GATK4 to separate SNVS and INDELs;

[0053] (4.3) Use the VariantFiltration module of GATK4 to filter low-quality mutation sites and obtain a vcf file containing accurate SNP information;

[0054] (4.4) Use the GenotypeGVCFs module of GATK4 to perform SNP annotation based on the dbSNP138 database, and label each SNP with an rs number (rs number, also known as rs number, is a systematic identifier applied in human genomics to locate a single position in genotype analysis);

[0055] (4.5) Use PLINK v1.90b6.21 64-bit to perform site and sample quality control on the obtained SNP vcf data;

[0056] (5) Using Plink for GWAS analysis and PRS construction combined with smoking factors: Using Plink software and SNP data, the BED file was analyzed using a logistic regression model to obtain the GWAS overall data. The obtained GWAS overall data was used as the baseline data, and another set of independent samples was used as the data target. Based on the P+T method, PRS was constructed. The P+T method is a p-value clumping + thresholding method, which only includes a part of the SNPs in the calculation of the PRS. That is, clustering (pruning based on the p-value) is first performed to screen out the SNPs with the lowest p-value in each module, and then the SNPs included are selected based on a certain threshold of the p-value.

[0057] The present invention collects samples, constructs libraries for sequencing, and performs data quality control to perform data comparison to obtain the initial comparison results. After obtaining the physical location of reads in the genome, GATK4 is used to identify SNPs and genotypes, and Plink is used to perform GWAS analysis and construct PRS in combination with smoking factors: Plink software is used in combination with SNP data, and the BED file is analyzed using a logistic regression model to obtain the overall GWAS data. Based on the obtained GWAS overall data as baseline data, another set of independent samples is used as data targets, and PRS is constructed based on the P+T method. Therefore, it is possible to accurately distinguish those who are genetically susceptible to smoking-related COPD in the Chinese population through bioinformatics methods, so that potentially high-risk groups can prevent COPD early.

[0058] Preferably, in step (3), the filtering conditions are as follows: remove adapters; remove bases with a 3' and 5' quality score less than 3; set a 4-base sliding window and remove bases with an average quality score less than 15; remove reads less than 36 bp in length; and then check the quality of the sequencing data after quality control using FastQC software. Reads are short fragment sequences obtained during the DNA sequencing process, typically 50 to 500 base pairs in length. These short fragments can be spliced ​​into a complete DNA sequence to determine the sequence of the original DNA molecule.

[0059] Preferably, the data comparison in step (4) includes: using the MEM module in BWA software to compare, using default parameters, and comparing the high-quality data obtained by quality control to the reference genome, the reference genome is hg38 of UCSC; sorting the original comparison results by the sort command of SAMtools to obtain a file in BAM format; using the default parameters of Picard to remove optical and PCR duplicates to obtain the initial comparison results.

[0060] Preferably, in step (4.3), the parameters are as follows:

[0061] QD<2.0, FS>60.0, MQ<40.0, MQRankSum<-12.5, ReadPosRankSum<-8.0, indels: QD<2.0, FS>200.0, ReadPosRankSum<-20.0.

[0062] Among them, QD is QaulByDepth (the credibility of the variant site divided by the number of unfiltered non-reference reads), FS is FisherStrand (Fisher's exact test assesses the possibility that the current variant is a strand bias), MQ is MappingQuality (the square root of the alignment quality in all samples), MQRankSum is MappingQualityRankSumTest (the credibility is assessed based on the alignment quality of the REF and ALT reads), and ReadPosRankSum is ReadPosRankSumTest (the credibility of the variant is assessed by the position of the variant in the read, and the error rate is usually higher at the two ends of the read).

[0063] Preferably, the step (4.5) comprises the following sub-steps:

[0064] (4.5.1) Delete SNP sites with a missing rate greater than 5% in the sample;

[0065] (4.5.2) Delete samples with missing rate greater than 2%, a total of 2 samples were removed;

[0066] (4.5.3) Filtering was performed using Hadwinberg equilibrium, with a threshold of P value < 10 -6 ;

[0067] (4.5.4) Filter by minimum allele frequency, with a threshold of MAF < 0.01;

[0068] (4.5.5) Use snpEff software to perform gene and function annotation for all sites.

[0069] Preferably, in step (5), the PRS construction includes the following sub-steps:

[0070] (5.1) Use the coord module of LDpred software to keep the genomic positions of the two loci and genotype data consistent;

[0071] (5.2) Using the P+T module and score module of the LDpred software, based on the LD relationship between genetic variants, several P+T standards were set to generate corresponding PRSs. The PRS with the greatest association with COPD was considered optimal. The logistic regression model based on the RMS package was used to analyze the PRS.

[0072] Preferably, in step (5.2), there are 13 P+T standards, namely: P-value<3×10 -1 ,1×10 -1 ,3×10 -2 ,1×10 -2 ,3×10 -3 ,1×10 -3 ,3×10 -4 ,1×10 -4 ,3×10 -5 ,1×10 -5 ,1×10 -6 ,1×10 -7 ,1×10 -8 .

[0073] Preferably, in the PRS of step (5.2), the smoking factors include: age at first smoking, years of smoking, and number of cigarettes smoked per day.

[0074] A kit for analyzing smoking-related COPD susceptibility loci is also provided. This kit is used for mining smoking-related COPD susceptibility loci. Genomic DNA is used as the material, DNA fragmentation and end-repair and A addition are performed, and adapters with single-molecule barcode sequences are ligated to both ends using DNA ligase. A library is obtained after amplification. The library is hybridized with biotin-labeled probes in liquid phase and captured and enriched using streptavidin-coated magnetic beads. Finally, primers with tag sequences are used for amplification to obtain a captured library. Sequencing data is obtained through high-throughput sequencing, and the final gene variation information is obtained after analysis using the corresponding data analysis software.

[0075] The components of Part A of the kit are as follows:

[0076] Reagent 1: The component name is DNA fragmentation and repair enzyme, and its ingredients are DNA fragmentation enzyme and repair enzyme. The specification and quantity are 96 μL × 1 tube;

[0077] Reagent 2: The component name is interruption repair solution, which includes Tris, and the specification and quantity are 96 μL × 1 tube;

[0078] Reagent 3: The component name is ligation buffer, which includes Tris and ATP. The specification and quantity are 192 μL × 1 tube;

[0079] Reagent 4: The component name is ligase, the ingredients include ligase, the specification and quantity are 96 μL × 1 tube;

[0080] Reagent 5: The component name is adapter, the ingredients include oligonucleotides, the specifications and quantity are 3μL × 24 tubes;

[0081] Reagent 6: The component name is primer, the ingredients include oligonucleotides, the specification and quantity are 70 μL × 1 tube;

[0082] Reagent 7: The component name is amplification solution, which includes Tris, dNTPs, and DNA polymerase. The specification and quantity are 1390 μL × 1 tube;

[0083] Reagent 8: The component name is blocking solution, which includes oligonucleotides, and the specification and quantity are 88 μL × 1 tube;

[0084] Reagent 9: The component name is hybridization solution, which includes Tween and dextran sulfate. The specification and quantity are 336 μL × 1 tube;

[0085] Reagent 10: Hybridization Enhancer, containing formamide, 150 μL per tube.

[0086] Reagent 11: The component name is capture probe, including oligonucleotides, and the specification and quantity are 150 μL × 1 tube;

[0087] Reagent 12: The component name is washing solution S, which includes ethylenediaminetetraacetic acid disodium salt, and the specification and quantity are 576 μL × 1 tube;

[0088] Reagent 13: The component name is washing solution 1, which includes sodium dodecyl sulfate, and the specification and quantity are 480 μL × 1 tube;

[0089] Reagent 14: The component name is washing solution 2, which includes sodium dodecyl sulfate, and the specification and quantity are 576 μL × 1 tube;

[0090] Reagent 15: The component name is washing solution 3, which includes sodium dodecyl sulfate, and the specification and quantity are 576 μL × 1 tube;

[0091] Reagent 16: Magnetic bead cleaning solution, containing Tris, EDTA, and sodium chloride. Quantity: 720 μL x 4 tubes.

[0092] Reagent 17: The component name is label, the ingredients include oligonucleotides, the specifications and quantity are 4 μL × 24 tubes;

[0093] Reagent 18: The component name is positive control, the ingredients include deoxyribonucleic acid, the specification and quantity are 24μL×1 tube;

[0094] Reagent 19: The component name is negative control, the ingredients include deoxyribonucleic acid, the specification and quantity are 24μL×1 tube;

[0095] The components of Part B of the kit are as follows:

[0096] Reagent 20: The component name is capture magnetic beads, and its ingredients include streptavidin magnetic beads. The specification and quantity are 880 μL × 1 tube.

[0097] Preferably, 200 bp regions upstream and downstream of the SNP were selected to design probes. The web-based probe design software was used to design and evaluate the probes. The probe information is as follows:

[0098] Site 1: rs763697172, located at chr12:104546413-104546553, probe coverage area: 141 bp, CHST11;

[0099] Site 2: snp5_682363, located at chr5:682303-682423, probe coverage area: 121 bp, TPPP;

[0100] Site 3: rs56408533, located at chr2:63534400-63534580, probe coverage area: 181 bp, WDPCP;

[0101] Site 4: rs1232553, located at chr11:60866012-60866211, probe coverage area: 200 bp, ZP1;

[0102] Site 5: snp11_26449753, located at chr11:26449682-26449802, probe coverage area: 121 bp, ANO3;

[0103] Site 6: snp17_44011024, located at chr17:44011004-44011124, probe coverage area: 121 bp, NAGS;

[0104] Site 7: snp14_10246376, located at chr14:102463246-102463386, probe coverage area: 141 bp, TECPR2;

[0105] Site 8: rs573468786, located at chr7:14415389-14415609, probe coverage area: 221bp, DGKB

[0106] Site 9: rs386419685, located at chr2:108591764-108591884, probe coverage area: 121 bp, LIMS1;

[0107] Site 10: rs7936710, located at chr11:94375411-94375654, probe coverage region: 244 bp, GPR83;

[0108] Site 11: rs7671167, located at chr4:88962808-88963108, probe coverage region: 301 bp, FAM13A;

[0109] Site 12: rs11347214, located at chr4:88962812-88963064, probe coverage area: 253bp, FAM13A.

[0110] The embodiments of the present invention are described in detail below.

[0111] Example 1

[0112] Data collection and library construction and sequencing: samples from patients with smoking-related COPD and normal individuals were collected. DNA extraction, library construction and sequencing were performed on all samples. The library was constructed using the Illumina whole genome amplification kit, the sequencing platform was Illumina's NovaSeq, and the target depth was 10x. Quality control was performed at every step from sample collection to library construction and sequencing, and unqualified samples were supplemented or removed in a timely manner. The total sample library contains 2795 samples, including 1415 samples from patients with smoking-related COPD and 1380 normal samples (see Table 1).

[0113] Table 1

[0114]

[0115] Data Preprocessing and Mutation Calling: All off-machine FastQ data were first subjected to preliminary quality control (QC) for each sample, including data quality assessments such as average base quality, GC content, N content, and Ts / Tv ratio. Results showed that all 2,795 samples had data volumes exceeding 30 GB, Q30 values ​​exceeding 85%, and an average sequencing depth exceeding 10×, meeting relevant standards. BAM files were then generated using BWA (Brain Mapping and Analysis) to align against the hg38 human reference genome downloaded from UCSC. Finally, mutation calling was performed across multiple samples using GATK HaplotypeCaller. Given the low sequencing depth of whole genomes, simultaneous mutation calling across multiple samples can reduce false positives. To ensure the accuracy of SNP and indel analysis, GATK VariantFiltration was used for further mutation QC, ultimately identifying over 8,085,738 SNPs and indels.

[0116] GWAS analysis was performed on the VCF files of all samples using Plink software. Plink integrates common GWAS functions, including chi-square tests for individual SNPs / indels, Wald tests, and Bonferroni and Benjamin Hochber P-value corrections for multiple testing. Manhattan plots were generated using the qqman module (0.1.8) in the R language (version 4.1.1). The concordance of loci obtained from the entire sample and smoker-only samples at different P values ​​was then compared.

[0117] In order to further verify the accuracy of the SNPs reported in this paper, the sites related to smoking-related COPD found in this paper were compared with SNPs reported in previous literature. The results showed that there were a large number of overlapping sites between the smoking-related COPD genes obtained in this paper and those reported previously (P-value < 5e -4 ), which further proves the accuracy of the SNP reported in the present invention (see Figure 2 ).

[0118] Enrichment analysis: For the significant SNPs obtained, we used the enrichment analysis software (Metascape, https: / / metascape.org / gp / index.html# / main / ) to perform enrichment analysis on the genes where the SNPs were located. It was found that genes with more significant P values ​​were more enriched in pathways related to lung function and human lifespan (see Figure 3 After removing the sites of non-smokers, the enrichment analysis included the results of pathways related to lung disease and smoking.

[0119] PRS construction and example testing, using the P+T (pruning+thresholding) method of LDpred (version 1.0.10) software to construct a polygenic risk score:

[0120] 1) Using the coord module of the LDpred software, we extracted the common sites between the baseline and target genotype data. The target data included 267 disease samples and 483 normal samples, and kept the genomic locations of the genotype data consistent.

[0121] 2) Using the P+T module and score module of LDpred software, 13 P+T standards were set based on the LD relationship between genetic variants (--ldr 200), namely: P-value < 3×10 -1 ,1×10 -1 ,3×10 -2 ,1×10 -2 ,3×10 -3 ,1×10 -3 ,3×10 -4 ,1×10 -4 ,3×10 -5 ,1×10 -5 ,1×10 -6 ,1×10 -7 ,1×10 -8 The corresponding PRS was generated, and the PRS with the greatest correlation with COPD was the best;

[0122] 3) The R packages used for the ROC graph are as follows: R version 4.1.1, R package ggplot2 (3.4.0), R package rms (1.18.0), R package pROC (1.18.0), and R package forestmodel (0.6.2). Finally, a logistic regression model based on the R package rms (1.18.0) was used to analyze PRS+smoking factors, where smoking factors included smoking age and years of smoking. The results showed that the use of susceptibility loci and PRS models of European populations could not accurately predict whether the Chinese population had smoking-related COPD, but the use of susceptibility loci and the constructed PRS model of the Chinese population could provide a relatively accurate prediction for the Chinese population (see Figure 4 ).

[0123] This study, the first large-scale GWAS study of smoking-related COPD in Chinese samples, used thousands of GWAS data sets as the sample set and identified SNPs significantly associated with the onset of smoking-related COPD. Furthermore, the PRS model constructed by this study (incorporating smoking factors) can also relatively accurately predict individuals at high risk of smoking-related COPD. The list of selected SNPs is shown in Table 2.

[0124] Table 2

[0125]

[0126] Example 2

[0127] 1. Genomic DNA Extraction

[0128] Genomic DNA was extracted from fresh or frozen whole blood (blood treated with anticoagulants such as citrate, EDTA, or heparin) using a commercial kit.

[0129] Total sample volume: Applicable to starting amounts of 1 ng to 1000 ng DNA. To improve data quality, a minimum of 50 ng DNA is recommended. DNA concentration is determined using Qubit fluorescence quantification.

[0130] Extracted sample quality and purity: DNA integrity and protein residue can be assessed by agarose gel electrophoresis. Nanodrop A260 / A280 = 1.8-2.0; A260 / A230 > 2.0.

[0131] 2. DNA fragmentation, end repair and A addition

[0132] After the extraction of genomic DNA, the fragmentation repair solution and fragmentation repair enzyme were added, mixed and placed in a PCR instrument for incubation; the PCR program settings were: 37°C, 30 min; 65°C, 30 min; 4°C, ∞.

[0133] Note: Take 5uL of the positive control and negative control of the kit and test them simultaneously with the samples to be tested.

[0134] 3. Add connector

[0135] Add adapters, ligase, and ligation buffer to the end-repair plus A product, mix well, and place in a PCR instrument for incubation at 20°C for 15 minutes.

[0136] 4. Purification of Ligation Products

[0137] Add 1.0 times the purified magnetic beads to the ligation product, mix well and let it stand; place the PCR tube in a magnetic rack, let it stand, and discard the liquid; wash twice with freshly prepared 80% ethanol; finally, add nuclease-free water for eluation.

[0138] 5. Amplification of Ligation Products

[0139] Add primers and amplification solution to the purified ligation product, mix well and place in a PCR instrument for amplification. The PCR program is set as follows: 98°C, 45s; (98°C, 15s; 65°C, 30s; 72°C, 30s) 11-13 cycles; 72°C, 1min; 4°C, ∞.

[0140] 6. PCR Product Purification

[0141] Add 1.5 times the purified magnetic beads to the amplified product, wash twice with freshly prepared 80% ethanol, and elute with nuclease-free water.

[0142] The library concentration was detected using Qubit, and the library fragment length was detected using a fragment analyzer.

[0143] 7. Library and Probe Hybridization

[0144] The library was concentrated in a vacuum concentrator. After concentration, the hybridization reaction solution (probe, blocking solution, and hybridization buffer) was added, vortexed to ensure that the DNA dried at the bottom of the tube was dissolved, and briefly centrifuged; hybridization was performed overnight at 65°C on a PCR instrument.

[0145] 8. Preparation of capture magnetic beads: Mix the capture magnetic beads equilibrated at room temperature thoroughly, take the mixed capture magnetic beads and add 1× magnetic bead cleaning solution, mix well, centrifuge briefly, place the PCR tube on a magnetic rack until the solution is clear, and remove the supernatant; repeat the wash three times, discard the supernatant and centrifuge for the last time, and place it on a magnetic rack until the liquid is clear, and remove the residual liquid; remove the PCR tube and immediately add magnetic bead resuspension buffer (prepared with hybridization solution and hybridization enhancer), centrifuge briefly, vortex mix at low speed, if there are bubbles, flick to remove the bubbles, centrifuge again, place the PCR tube in a 65℃ metal bath or PCR instrument to preheat and set aside.

[0146] 9. Target region DNA capture

[0147] Mix the capture beads by slowly pipetting them up and down. Quickly transfer them to the hybridization product (performed on a PCR instrument). Pipet them up and down 10 times to mix them again. Centrifuge briefly, quickly return the PCR tube to the instrument, and run the enrichment program. Set the instrument to 65°C, the heated lid to 70°C, and the reaction volume to 100 μL.

[0148] Take out the sample every 10-15 minutes, vortex at low speed to avoid bubbles, centrifuge briefly (keep the magnetic beads suspended), and quickly return it to the PCR instrument.

[0149] 10. Preparation of cleaning solution

[0150] The washing solution was placed in a water bath and completely dissolved to prepare a 1× working solution;

[0151] Preheat Wash Buffer 1 and 1× Wash Buffer S in a 65°C metal bath or PCR instrument.

[0152] 11.65℃ washing

[0153] After the enrichment reaction is completed, open the PCR tube cap of the captured product directly on the PCR instrument (the 65℃ Hold program continues to run), quickly add 1× Wash Solution 1 preheated at 65℃ to the captured sample, pipette to mix, centrifuge briefly, quickly place the PCR tube on a magnetic rack until the solution is clear, and discard the supernatant.

[0154] Place the PCR tube back on the PCR instrument, quickly add 1× Wash Buffer S preheated at 65°C, mix well with a pipette, centrifuge briefly, place the PCR tube back on the PCR instrument, and incubate for 5 minutes; repeat this step three times.

[0155] 12. Washing at room temperature

[0156] Remove the PCR tube from the magnetic rack, add 1× Wash Buffer 1 at room temperature, vortex to mix, centrifuge briefly, place the PCR tube on the magnetic rack until the solution is clear, and discard the supernatant.

[0157] Remove the PCR tube from the magnetic rack, add 1× Wash Buffer 2 at room temperature, vortex to mix, centrifuge briefly, place the PCR tube on the magnetic rack until the solution is clear, and discard the supernatant;

[0158] Remove the PCR tube from the magnetic rack, add 1× Wash Buffer 3 at room temperature, vortex to mix, centrifuge briefly, place the PCR tube on the magnetic rack until the solution is clear, and discard the supernatant; centrifuge again, place the PCR tube on the magnetic rack, and use a pipette to completely remove the remaining supernatant.

[0159] Remove the PCR tube from the magnetic stand, add deionized water, vortex to mix, centrifuge, retain the magnetic beads, and proceed to the next reaction.

[0160] 13. Post-capture PCR Amplification

[0161] Sequencing primers and amplification solution were added to the resuspended capture magnetic beads, mixed well, and placed in a PCR instrument for amplification. The PCR program was set as follows: 98°C, 45s; (98°C, 15s; 60°C, 30s; 72°C, 30s) 11-13 cycles; 72°C, 1min; 4°C, ∞.

[0162] 14. Post-amplification purification

[0163] Remove the purified magnetic beads and equilibrate to room temperature. Add 1.1 times the volume of purified magnetic beads to the amplified product, mix well, and let it stand. Place the PCR tube in a magnetic rack, let it stand, and discard the liquid. Wash twice with freshly prepared 80% ethanol. Finally, add nuclease-free water for eluent.

[0164] 15. Library Quality Control

[0165] Take 1 μL of library and use Qubit dsDNA HS Assay Kit reagent to measure the library concentration on Qubit 4.0 Fluorometer and record the library concentration; take 1 μL of library and use fragment analyzer to perform fragment quality inspection. The fragment size should be basically consistent with the pre-library size.

[0166] 16. Sequencing

[0167] The library was sequenced using the NGS sequencing platform, and sequencing was performed in PE150 mode according to the instrument instructions.

[0168] 17. Data preprocessing and mutation detection

[0169] For all fastq data after being processed, we first performed preliminary quality control on each sample using fastqc to check the data quality, including average base quality, GC content, N content, and Ts / Tv ratio. The results showed that the sample data volume after quality control was greater than 30GB, the Q30 was greater than 85%, and the average sample sequencing depth was greater than 10×, meeting the corresponding standards. We then used BWA to align the data to generate BAM files. The genome was aligned to the hg38 human reference genome downloaded from UCSC. We used GATK VariantFiltration for mutation quality control and GATK haplotypeCaller for multi-sample mutation search. Finally, we determined the sample genotypes.

[0170] 18. Sample genotype interpretation

[0171] Table 3

[0172]

[0173] By inputting each sample's typing into the PRS model, a PRS score was obtained for each sample. The scores for these five samples were 0.9, 0.12, 0.27, 0.55, and 0.13, respectively. A PRS score close to 1 indicates a high risk of smoking-related COPD, while a PRS score close to 0 indicates a low risk (see Table 3).

[0174] The above description is merely a preferred embodiment of the present invention and does not constitute any form of limitation to the present invention. Any simple modifications, equivalent changes and modifications made to the above embodiments based on the technical essence of the present invention are still within the scope of protection of the technical solution of the present invention.

Claims

1. A method for analyzing smoking-related COPD susceptibility loci, characterized by: It includes the following steps: (1) Sample collection: Samples were collected from patients with smoking-related COPD and healthy individuals; (2) Library construction and sequencing: DNA extraction, library construction, and sequencing were performed on all samples. The library construction used the Illumina whole genome amplification kit, the sequencing platform was Illumina's NovaSeq, and the target depth was 10x; (3) Data quality control: Use the bbduk.sh script in the Trimmomatic and BBTools components to obtain high-quality valid data; (4) Mutation identification: Perform data alignment to obtain the initial alignment results; after obtaining the physical location of the reads in the genome, use GATK4 to identify SNPs and genotype: (4.1) Use the HaplotypeCaller module of GATK4 with default parameters to obtain the original mutation information; (4.2) Use the SelectVariants module of GATK4 to separate SNVS and INDELs; (4.3) Use the VariantFiltration module of GATK4 to filter low-quality mutation sites and obtain a vcf file containing accurate SNP information; (4.4) Use the GenotypeGVCFs module of GATK4 to perform SNP annotation based on the dbSNP138 database, and label each SNP with an RS number; (4.5) Use PLINK v1.90b6.21 64-bit to perform site and sample quality control on the obtained SNP vcf data; (5) Plink was used to perform GWAS analysis and construct PRS in combination with smoking factors: Plink software was used in combination with SNP data, and the BED file was analyzed using a logistic regression model to obtain the overall GWAS data. The obtained overall GWAS data was used as baseline data, and another group of independent samples was used as data targets. PRS was constructed based on the P+T method.

2. The method for analyzing smoking-related COPD susceptibility loci according to claim 1, characterized in that: In step (3), the filtering conditions are as follows: remove the linker; remove bases with a quality value lower than 3 at the 3' and 5' ends; set a sliding window of 4 bases and remove bases with an average quality value lower than 15; Reads shorter than 36 bp were removed, and the quality of the sequencing data after quality control was checked using FastQC software.

3. The method for analyzing smoking-related COPD susceptibility loci according to claim 2, characterized in that: The data comparison in step (4) includes: using the MEM module in the BWA software to compare, using the default parameters, and comparing the high-quality data obtained by quality control to the reference genome, which is the hg38 of UCSC; sorting the original comparison results by the sort command of SAMtools to obtain a file in BAM format; using the default parameters of Picard to remove optical and PCR duplicates to obtain the initial comparison results.

4. The method for analyzing smoking-related COPD susceptibility loci according to claim 3, characterized in that: In the step (4.3), the parameters are as follows: QD<2.0, FS>60.0, MQ<40.0, MQRankSum<-12.5, ReadPosRankSum<-8.0, indels: QD<2.0, FS>200.0, ReadPosRankSum <-20.

0.

5. The method for analyzing smoking-related COPD susceptibility loci according to claim 4, characterized in that: The step (4.5) comprises the following sub-steps: (4.5.1) Delete SNP sites with a missing rate greater than 5% in the sample; (4.5.2) Delete samples with missing rate greater than 2%, a total of 2 samples were removed; (4.5.3) Filtering was performed using Hadwinberg equilibrium, with a threshold of P value < 10 -6 ; (4.5.4) Filter by minimum allele frequency, with a threshold of MAF < 0.01; (4.5.5) Use snpEff software to perform gene and function annotation for all sites.

6. The method for analyzing smoking-related COPD susceptibility loci according to claim 5, characterized in that: In step (5), the PRS construction includes the following sub-steps: (5.1) Use the coord module of LDpred software to keep the genomic positions of the two loci and genotype data consistent; (5.2) Using the P+T module and score module of the LDpred software, based on the LD relationship between genetic variants, several P+T standards were set to generate corresponding PRSs. The PRS with the greatest association with COPD was considered optimal. The logistic regression model based on the RMS package was used to analyze the PRS.

7. The method for analyzing smoking-related COPD susceptibility loci according to claim 6, characterized in that: In the step (5.2), there are 13 P+T standards, namely: P-value < 3 × 10 -1 ,1×10 -1 ,3×10 -2 ,1×10 -2 ,3×10 -3 ,1×10 -3 ,3×10 -4 ,1×10 -4 ,3×10 -5 ,1×10 -5 ,1×10 -6 ,1×10 -7 ,1×10 -8 .

8. The method for analyzing smoking-related COPD susceptibility loci according to claim 7, characterized in that: In the PRS of step (5.2), smoking factors include: age at first smoking, years of smoking, and number of cigarettes smoked per day.

9. A kit for analyzing smoking-related COPD susceptibility loci according to claim 8, characterized in that: This kit is used to explore smoking-related COPD susceptibility loci. Genomic DNA is used as the material, DNA fragmentation is performed, end-repair and A addition is performed, and adapters with single-molecule barcode sequences are ligated to both ends using DNA ligase. The library is then amplified. The library is hybridized with a biotin-labeled probe in liquid phase and captured and enriched using streptavidin-coated magnetic beads. Finally, primers with tag sequences are used for amplification to obtain the captured library. The sequencing data is obtained through high-throughput sequencing, and the final gene variation information is obtained after analysis by the supporting data analysis software; The components of Part A of the kit are as follows: Reagent 1: The component name is DNA fragmentation and repair enzyme, and its ingredients are DNA fragmentation enzyme and repair enzyme. The specification and quantity are 96 μL × 1 tube; Reagent 2: The component name is interruption repair solution, which includes Tris, and the specification and quantity are 96 μL × 1 tube; Reagent 3: The component name is ligation buffer, which includes Tris and ATP. The specification and quantity are 192 μL × 1 tube; Reagent 4: The component name is ligase, the ingredients include ligase, the specification and quantity are 96 μL × 1 tube; Reagent 5: The component name is adapter, the ingredients include oligonucleotides, the specifications and quantity are 3 μL × 24 tubes; Reagent 6: The component name is primer, the ingredients include oligonucleotides, the specification and quantity are 70 μL × 1 tube; Reagent 7: The component name is amplification solution, which includes Tris, dNTPs, and DNA polymerase. The specification and quantity are 1390 μL × 1 tube; Reagent 8: The component name is blocking solution, which includes oligonucleotides, and the specification and quantity are 88 μL × 1 tube; Reagent 9: The component name is hybridization solution, which includes Tween and dextran sulfate. The specification and quantity are 336 μL × 1 tube; Reagent 10: Hybridization Enhancer, containing formamide, 150 μL x 1 tube. Reagent 11: The component name is capture probe, including oligonucleotides, and the specification and quantity are 150 μL × 1 tube; Reagent 12: The component name is washing solution S, which includes disodium ethylenediaminetetraacetic acid, and the specification and quantity are 576 μL × 1 tube; Reagent 13: The component name is washing solution 1, which includes sodium dodecyl sulfate, and the specification and quantity are 480 μL × 1 tube; Reagent 14: The component name is washing solution 2, which includes sodium dodecyl sulfate, and the specification and quantity are 576 μL × 1 tube; Reagent 15: The component name is washing solution 3, which includes sodium dodecyl sulfate, and the specification and quantity are 576 μL × 1 tube; Reagent 16: Magnetic bead cleaning solution, containing Tris, EDTA, and sodium chloride, in 4 tubes (720 μL). Reagent 17: The component name is label, the ingredients include oligonucleotides, the specifications and quantity are 4 μL × 24 tubes; Reagent 18: The component name is positive control, the ingredients include deoxyribonucleic acid, the specification and quantity are 24μL×1 tube; Reagent 19: The component name is negative control, the ingredients include deoxyribonucleic acid, the specification and quantity are 24μL×1 tube; The components of Part B of the kit are as follows: Reagent 20: The component name is capture magnetic beads, and its ingredients include streptavidin magnetic beads. The specification and quantity are 880 μL × 1 tube.

10. The kit according to claim 9, characterized in that: Probes were designed from 200 bp upstream and downstream of the SNP. The web-based probe design software was used to design and evaluate the probes. The probe information is as follows: Site 1: rs763697172, located at chr12:104546413-104546553, probe coverage area: 141 bp, CHST11; Site 2: snp5_682363, located at chr5:682303-682423, probe coverage area: 121 bp, TPPP; Site 3: rs56408533, located at chr2:63534400-63534580, probe coverage area: 181 bp, WDPCP; Site 4: rs1232553, located at chr11:60866012-60866211, probe coverage area: 200 bp, ZP1; Site 5: snp11_26449753, located at chr11:26449682-26449802, probe coverage area: 121 bp, ANO3; Site 6: snp17_44011024, located at chr17:44011004-44011124, probe coverage area: 121 bp, NAGS; Site 7: snp14_10246376, located at chr14:102463246-102463386, probe coverage area: 141 bp, TECPR2; Site 8: rs573468786, located at chr7:14415389-14415609, probe coverage area: 221bp, DGKB Site 9: rs386419685, located at chr2:108591764-108591884, probe coverage area: 121 bp, LIMS1; Site 10: rs7936710, located at chr11:94375411-94375654, probe coverage region: 244 bp, GPR83; Site 11: rs7671167, located at chr4:88962808-88963108, probe coverage region: 301 bp, FAM13A; Site 12: rs11347214, located at chr4:88962812-88963064, probe coverage area: 253bp, FAM13A.

Images