A personalized reference genome construction method and system for cancer genomics somatic mutation detection

Abstract

The application provides a personalized reference genome construction method and system for cancer genomics somatic mutation detection, and belongs to the field of cancer genomics. The construction method comprises five steps: 1. Using GATK and pindel to detect germ line mutations, and generating a high-confidence variation VCF file through merging, screening, deduplication and sorting; 2. Replacing the alleles of heterozygous sites to construct a "decoy" genome; 3. Constructing a homozygous reference genome based on homozygous mutations, and retaining the original bases of the heterozygous sites; 4. Aligning sequencing reads to the decoy genome, and re-aligning to the homozygous genome after correcting the alleles of the heterozygous sites; and 5. Outputting a standardized BAM file for use by downstream tools. The method can eliminate germ line variation mismatches, reduce false positives, improve somatic mutation detection accuracy, and has a stable process and does not require adjustment of downstream parameters.

Description

Technical Field

[0001] This invention relates to the fields of bioinformatics and cancer genomics, and in particular to a method for constructing a personalized reference genome. Background Technology

[0002] Somatic mutations play a crucial role in the occurrence and development of cancer, but the accuracy of their detection faces challenges. This is partly because the frequency of somatic mutations is extremely low, ranging from 0.1 to 100 mutations per million base pairs, and the allele frequency is not fixed.

[0003] However, the main reason is the high false positive rate. Existing tools (such as MuTect, VarScan2, etc.) rely on a universal reference genome, and the mismatch between the reads obtained from tumor sample sequencing and the reference genome leads to a high false positive rate.

[0004] To address the high false positive rate in detection, the shortcomings of a universal reference genome must be addressed. Current solutions include:

[0005] (1) Dual reference genome method: Construct two reference genomes, one containing reference alleles and the other containing variant alleles; or,

[0006] (2) Fragment insertion method: The variant allele and its flanking sequence are inserted as additional fragments at the end of the reference genome.

[0007] Existing methods have significant limitations. First, reduced coverage causes reads to be scattered across two locations, making it difficult for somatic mutation detection tools to accurately calculate mutation frequencies. Second, inconsistent coordinates can lead to alignment coordinate shifts caused by inserted fragments, affecting subsequent analyses such as mutation annotation. Furthermore, many personalized reference tools offer insufficient support for structural variations, only supporting SNPs and small indels while ignoring large insertions / deletions (>50bp), inversions, and other structural variations. Summary of the Invention

[0008] This invention provides a method and system for constructing a personalized reference genome to improve the accuracy of somatic mutation detection in cancer genomics.

[0009] This invention first provides a method for constructing a personalized reference genome for somatic mutation detection in cancer genomics, the method comprising the following steps:

[0010] Step 1: Germline mutation detection: Using patient-matched normal tissue samples, analyze their whole-genome sequencing data using standard variant detection tools to identify all germline mutations and generate a high-confidence germline mutation list; germline mutations include single nucleotide polymorphisms (SNPs) and insertion / deletion indels (insertion / deletion mutations);

[0011] Step 2: Constructing the decoy genome: In the universal reference genome, replace the reference alleles of all patient-specific heterozygous loci with alternative alleles;

[0012] Step 3: Construct a homozygous genome: Based on the germline mutation list from Step 1, construct a reference genome containing only homozygous mutations from the patient, while retaining the original reference bases for heterozygous sites;

[0013] Step 4: Sequencing reads correction and re-alignment: Align the original sequencing reads from tumor and normal samples to the bait genome to ensure that all reads carrying germline mutations are captured; restore the alternative alleles of heterozygous sites in the reads aligned to the bait genome to the reference alleles; re-align the corrected reads to the homozygous genome.

[0014] Step 5: Output standardized BAM files: Generate a pair of corrected BAM files for direct use by subsequent arbitrary somatic mutation detection tools.

[0015] Preferably, the germline mutation detection step includes: matching normal tissue, which is a healthy tissue sample from the same patient, to distinguish between germline mutations and somatic mutations.

[0016] Preferably, the decoy genome construction step includes: retrieving the patient germline mutation information from step 1, replacing the reference allele with the alternative allele at the patient-specific heterozygous sites in the universal reference genome, integrating all the patient's germline variation information, and then calibrating the sequence to obtain the decoy genome.

[0017] Preferably, the steps of constructing a homozygous genome include: based on the patient germline mutation information in step 1, replacing the patient homozygous bases in the universal reference genome, retaining the original bases of the heterozygous sites, and obtaining a homozygous genome containing only patient homozygous variants through sequence calibration.

[0018] Preferably, the sequencing reads correction and realignment steps include: mapping the original sequencing reads of tumor and normal samples to a decoy genome to capture reads carrying heterozygous variants, and replacing the bases at heterozygous sites with bases from a universal reference genome.

[0019] Preferably, the steps for outputting a standardized BAM file include: remapping the corrected reads to a homozygous genome containing only homozygous variants, and retaining the original reference genome coordinates in the final BAM file.

[0020] Preferably, when processing indels, bases are deleted or inserted to maintain a coordinate system consistent with the reference genome.

[0021] This invention also discloses a personalized reference genome construction system for somatic mutation detection in cancer genomics, the system comprising at least a processor and a memory connected to the processor, wherein:

[0022] The memory is used to store computer programs or instructions;

[0023] The processor is used to acquire and execute computer programs or instructions in memory to enable the system to execute a personalized reference genome construction method for somatic mutation detection in cancer genomics proposed in this invention.

[0024] Beneficial effects:

[0025] The method of this invention generates a final BAM file (tumor / normal sample) for direct use by any subsequent somatic mutation detection tool. The file coordinate system is consistent with the universal reference genome, seamlessly adapting to existing analysis workflows. The original reference genome coordinates are retained in the final BAM file. This allows the final BAM file to more accurately locate the correct positions, reducing mismatches caused by germline variations and thus improving the accuracy of subsequent somatic mutation detection tools. The MapVariants command of this invention can be used for coordinate transformation, facilitating the detection of somatic mutations on the final BAM using tools such as Mutect2 and VarScan, generating a somatic mutation VCF file. The ViewFasta command can also be used to view the reference genome region and the SomaticMutationsOnGermlineInsertion command, marking the relative positions of somatic mutations located in the germline insertion region. Attached Figure Description

[0026] Figure 1 This is a flowchart of the present invention. Detailed Implementation

[0027] The following embodiments are used to illustrate the present invention, but are not intended to limit the scope of the invention. The present invention will be further described in detail below with reference to the accompanying drawings and embodiments.

[0028] Preprocessing: GATK is mainly used for SNP and indel detection. The workflow includes alignment, sorting, deduplication, local repetition alignment, BQSR, variant detection, and filtering. Pindel is specifically used for indel detection, with parameter settings focusing on supporting read count, indel length, and paired-end read information. This invention provides details on parameter settings for GATK and Pindel when detecting germline mutations. The details are as follows:

[0029] The GATK workflow is used to detect germline mutations and mainly includes the following steps:

[0030] 1. Comparison

[0031] The `-bwa mem` option is used for comparison, and its parameters include:

[0032] bwa mem -M -R

[0033] @RG\tID:sample_1\tLB:sample_1\tPL:ILLUMINA\tPM:HISEQ\tSM:sample_1' $ref 1.fq 2.fq > aligned_reads.sam

[0034] -M: Marks shorter split hits as secondary.

[0035] -R: Add read group information (ID, LB, PL, PM, SM)

[0036] 2. Sorting and deduplication

[0037] - Use the Picard tool to sort and mark duplicate reads:

[0038] java -jar $picard SortSam INPUT=aligned_reads.sam OUTPUT=sorted_reads.bam SORT_ORDER=coordinate

[0039] java -jar $picard MarkDuplicates INPUT=sorted_reads.bam OUTPUT=dedup_reads.bam METRICS_FILE=metrics.txt

[0040] java -jar $picard BuildBamIndex INPUT=dedup_reads.bam

[0041] 3. Local re-comparison

[0042] - Use GATK's RealignerTargetCreator and IndelRealigner for indel realignment:

[0043] java -jar $gatk -T RealignerTargetCreator -R $ref -I dedup_reads.bam-o realignment_targets.list

[0044] java -jar $gatk -T IndelRealigner -R $ref -I dedup_reads.bam -targetIntervals realignment_targets.list -o realigned_reads.bam

[0045] 4. Base quality correction

[0046] - Use BaseRecalibrator and PrintReads for BQSR:

[0047] java -jar $gatk -T BaseRecalibrator -R $ref -I realigned_reads.bam -knownSites $dbsnp -o recal_data.table

[0048] java -jar $gatk -T PrintReads -R $ref -I realigned_reads.bam -BQSRrecal_data.table -o recal_reads.bam

[0049] 5. Mutation detection

[0050] - Use HaplotypeCaller to detect variants:

[0051] java -jar $gatk -T HaplotypeCaller -R $ref -I recal_reads.bam -o raw_variants_recal.vcf

[0052] 6. Filtration

[0053] - Filter SNPs and Indels separately:

[0054] SNP filtering:

[0055] java -jar $gatk -T VariantFiltration -R $ref -V raw_snps_recal.vcf--filterExpression 'QD < 2.0 || FS > 60.0 || MQ < 40.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0 || SOR > 4.0' --filterName "basic_snp_filter" -ofiltered_snps_final.vcf

[0056] Indel filtering:

[0057] java -jar $gatk -T VariantFiltration -R $ref -V raw_indels_recal.vcf --filterExpression 'QD < 2.0 || FS > 200.0 || ReadPosRankSum < -20.0 || SOR >10.0' --filterName "basic_indel_filter" -o filtered_indels_final.vcf

[0058] Pindel is used to detect indel variants, with the following parameters:

[0059] pindel2vcf -e 4 -ir 2 -il 10 -pr 1 -pl 10 -as 250

[0060] -e 4: Minimum number of reads supported (4)

[0061] -ir 2: Minimum number of reads for the indel region (2)

[0062] -il 10: Maximum length of indel (10bp)

[0063] -pr 1: Minimum number of supported pairs of reads (1)

[0064] -pl 10: Maximum insertion length (10bp) for paired-end reads

[0065] -as 250: Minimum length of the anchor sequence (250bp)

[0066] Combination Figure 1 A personalized reference genome construction method for somatic mutation detection in cancer genomics, specifically including the following steps:

[0067] S1: Detecting germline mutations

[0068] To clearly distinguish between a patient's hereditary variants and potential somatic mutations, whole-genome sequencing data was analyzed using patient-matched normal tissue samples (such as blood or healthy tissue) and standard variant detection tools (GATK and Pindel) to identify all germline mutations (including SNPs and indels), generating a high-confidence list of germline mutations. Specifically, GATK and Pindel were used to identify germline mutations, followed by processing the variant VCF file using the commands CombineVariants, SelectGenotype, RemoveOverlaps, and SortVariants. These commands were used to merge, filter, remove duplicates, and sort variants, generating a clean variant VCF file for subsequent generation of a personalized reference genome.

[0069] The RemoveOverlaps command described in this invention is based on the Bedtools intersect tool to remove overlapping mutation sites, and the SelectGenotype command is based on the GATK SelectVariants command (parameter: --genotype-filter) to screen for homozygous and heterozygous mutations; SortVariants and CombineVariants are native GATK commands, and their parameter settings follow the official GATK specifications.

[0070] S2: Constructing the "decoy" genome

[0071] To capture sequencing reads carrying germline mutations and prevent alignment failures or shifts due to reference genome mismatches, all patient-specific heterozygous sites (e.g., reference base A, patient is A / C heterozygous) in the universal reference genome are replaced with alternative alleles (base C). This ensures that reads carrying alternative alleles can correctly align to the decoy genome, avoiding filtering or alignment to incorrect positions due to mismatches with the universal reference genome (base A). Specifically, the `MakePersonalizedReference` command is used to generate a decoy genome file (decoy_ref.fa), which can generate two haplotype reference genomes. The "decoy" reference genome contains all variants, while the homozygous reference genome contains only homozygous variants; therefore, different `-genotype` parameters are needed when generating two reference genomes during `MakePersonalizedReference`. This step uses the `-genotype -heter` parameter.

[0072] The MakePersonalizedReference command described in this invention implements the core base substitution logic based on the consensus module of bcftools (version 1.15+), and completes the genome index construction by combining the faidx command of samtools (version 1.15+). By encapsulating the operation process of the above-mentioned existing tools, it realizes one-click generation of personalized reference genomes (bait genomes and homozygous genomes); wherein, the -genotype parameter -heter / homo corresponds to calling the -H2 parameter (heterozygous site replacement) and the default parameter (homozygous site replacement only) of bcftoolsconsensus, respectively.

[0073] S3: Constructing a homozygous genome

[0074] Based on the S1 germline mutation list, a reference genome is constructed containing only homozygous mutations from the patient (e.g., reference base A, where the patient is T / T homozygous), while heterozygous sites retain their original reference bases (e.g., base A). This avoids read alignment dispersion caused by heterozygous sites (e.g., reads at the same A / C heterozygous position may be assigned to different coordinates) and provides stable coverage, offering high-quality input data for subsequent somatic mutation detection tools. Specifically, the MakePersonalizedReference command is used to generate a homozygous genome file (homo_ref.fa), which can generate two haplotype reference genomes. The homozygous reference genome contains only homozygous variants, therefore the -genotype parameter homo is used in MakePersonalizedReference.

[0075] Meanwhile, after the decoy genome and homozygous genome are generated, a personalized database file is generated using the MakePersonalizedVariantsDB command. This file provides the necessary database support for downstream tools such as read correction, re-alignment, and somatic mutation detection, thus completing the construction of the reference genome and supporting database.

[0076] S4: Sequencing read correction and repetition alignment

[0077] First, the data is mapped to the decoy genome: raw sequencing reads from tumor and normal samples are aligned to the decoy genome, ensuring that all reads carrying germline variations are captured. Then, in the reads aligned to the decoy genome, surrogate alleles at heterozygous sites are restored to the reference alleles.

[0078] Secondary alignment to homozygous genome: The corrected reads are re-aligned to the homozygous genome to generate the final BAM file. By mapping to the "decoy" and homozygous genomes, variant reads are captured while heterozygous site interference is eliminated, ensuring high consistency in alignment location and coverage. The corrected reads carry only somatic mutation signals; germline mutations have been eliminated, reducing false positives.

[0079] The specific steps involve using the BWA tool to map the raw reads to the decoy genome, generating a raw BAM file. The ReplaceGenotype command then modifies the reads based on heterozygous mutations, replacing substituted alleles with reference alleles. The modified reads are then mapped again using BWA to a reference genome containing only homozygous variants, completing the core workflow of "raw data mapping - sequence correction - remapping" to generate a BAM file for somatic mutation detection. Additionally, the MakePersonalizedVariantsDB command generates a personalized database file, preparing the necessary database files for downstream tools.

[0080] The ReplaceGenotype command described in this invention is a customized naming for the function of "correcting sequencing read bases based on heterozygous mutations". Its core logic is based on the built-in functions of GATK's FastaAlternateReferenceMaker and ReadsPipelineSpark. By identifying the substitution alleles of heterozygous mutation sites in reads, it restores them to the reference alleles of the universal reference genome, thus completing the accurate correction of read sequences.

[0081] The MakePersonalizedVariantsDB command described in this invention is not a native command of existing tools, but a customized name for the "build a personalized variant database" function. Its core logic is based on the built-in functions of GATK's VariantDB and IndexFeatureFile. By integrating homozygous and heterozygous mutation VCF files, it builds a structured variant database and generates an index, providing necessary database support for downstream read correction and somatic mutation detection tools.

[0082] S5: Output standardized BAM file

[0083] The modified reads are remapped to homozygous genomes containing only homozygous variants, generating a pair of corrected BAM files (tumor / normal sample) for direct use by subsequent arbitrary somatic mutation detection tools.

[0084] After the standardized BAM file generated in this step enters the somatic mutation detection stage, the coordinate transformation can be completed using the MapVariants command to convert the alignment coordinates of the personalized reference genome into the coordinates of the universal reference genome, ensuring compatibility with the coordinate systems of downstream somatic mutation detection tools (MuTect2, VarScan, etc.). After the coordinate transformation, somatic mutation detection can be carried out, and a somatic mutation VCF file will be generated in the end.

[0085] The MapVariants command described in this invention is a customized name for the "Personalized Reference Genome BAM File Coordinate Transformation" function. Its core logic is to convert the alignment results of the personalized reference genome into the standard coordinates of the universal reference genome (hg38) by resetting the reference genome header information of the BAM file and correcting the alignment coordinates, so as to adapt to downstream somatic mutation detection tools.

[0086] The effectiveness of this invention has been verified on the DREAM Challenge, COLO-829, and simulated datasets. The specific data and results are as follows:

[0087] The DREAM Challenge dataset is a publicly available dataset from ICGC-TCGA, simulating tumor-normal paired sequencing data, containing known somatic mutations, with sequencing depths at different stages (20X-60X), and read lengths of 101 bp. The COLO-829 dataset was obtained from sequencing conducted by the Wellcome Trust Sanger Institute and is a validated set of somatic mutations in melanoma cell lines (containing 35,543 SNVs), with whole-genome sequencing coverage of 40X and 32X for tumor and matched normal tissues, respectively, and read lengths of 75 bp. The simulated data (VarSim) is derived from the NA12878 genome and DREAM Stage 3 mutations. It is artificially synthesized data, containing real germline variations and simulated somatic mutations, with sequencing depths at different stages (20X-60X), and read lengths of 50 / 75 / 100 / 150 bp.

[0088] DREAM Challenge Data: Using a personalized reference genome constructed with this method as a reference, four testing tools (MuTect, MuSE, Somatic-Sniper, and VarScan2) were applied to the DREAM Challenge data. The results showed that: FP reduction ranged from 8.5% to 55.5% (MuSE performed best, with a 55.5% reduction in FP); F1 scores improved, such as MuSE in Stage 3 from 0.847 to 0.894; the best performance was MuSE in Stage 3, where FP decreased from 1611 to 740, a reduction of 54.1%.

[0089] The TP (true positives), FP (false positives), and F1 (F1 score) of this invention on the DREAM Challenge dataset are shown in Table 1:

[0090] Table 1: Performance metrics of the present invention on the DREAM Challenge dataset

[0091]

[0092] COLO-829 Data: Using a personalized reference genome constructed using this method, four testing tools (MuTect, MuSE, Somatic-Sniper, and VarScan2) were applied to the COLO-829 data. The results showed: FP reduction from 7.2% to 23.5% (MuSE performed best, with a 23.5% decrease in FP); F1 score improvement, such as MuSE increasing from 0.894 to 0.913; and a new finding: a somatic point mutation (located in the germline insertion region) in the PDE1A gene (melanoma-related), previously undetectable by other tools, was detected.

[0093] The TP (true positives), FP (false positives), and F1 (F1 score) of this invention on the COLO-829 dataset are shown in Table 2:

[0094] Table 2: Performance metrics of the present invention on the COLO-829 dataset

[0095]

[0096] Simulated data (VarSim): Using a personalized reference genome constructed using this method as a reference, four testing tools (MuTect, MuSE, Somatic-Sniper, VarScan2) were applied to the simulated data. The results showed that: at different sequencing depths (20X-60X), the processing effect was more significant for short reads (50 bp) at a sequencing depth of 20X; the FP reduction rate of the four testing tools ranged from 59.3% to 96.8%, with MuSE showing the best effect, reducing FP by 96.8%; and the mapping quality was improved: compared to the universal reference genome, the personalized reference genome constructed using this method significantly reduced the number of mismapped reads in highly repetitive regions of the genome. This method also showed good improvement effects on low mappability regions.

[0097] The TP (true positives), FP (false positives), and F1 (F1 score) of this invention on the simulated dataset (VarSim) are shown in Table 3:

[0098] Table 3: Performance metrics of the present invention on the simulated dataset (VarSim)

[0099]

[0100] This invention has demonstrated its effectiveness on DREAM Challenge, COLO-829, and simulated data: it significantly reduces the false positive rate, by up to 96.8%, with particularly good results on short reads; it improves the mapping accuracy of low-mappability regions and reduces "dark region" errors; it discovers new mutations, such as somatic mutations in the PDE1A gene; it is parameter-independent, and performance improvements do not require parameter tuning, making it suitable for different datasets and tools (such as MuTect and MuSE).

[0101] The specific embodiments described herein are merely illustrative of the spirit of the invention. Those skilled in the art to which this invention pertains may make various modifications or additions to the described specific embodiments or use similar methods to substitute them, without departing from the spirit of the invention or exceeding the scope defined by the appended claims.

Claims

1. A method for constructing a personalized reference genome for somatic mutation detection in cancer genomics, characterized in that, The construction method includes the following steps: Step 1: Germline mutation detection: Using patient-matched normal tissue samples, whole-genome sequencing data are analyzed using standard variant detection tools to identify all germline mutations and generate a high-confidence germline mutation list; germline mutations include single nucleotide polymorphisms (SNPs) and insertion / deletion indels. Step 2: Constructing the decoy genome: In the universal reference genome, replace the reference alleles of all patient-specific heterozygous loci with alternative alleles; Step 3: Construct a homozygous genome: Based on the germline mutation list from Step 1, construct a reference genome containing only homozygous mutations from the patient, while retaining the original reference bases for heterozygous sites; Step 4: Sequencing reads correction and re-alignment: Align the original sequencing reads from tumor and normal samples to the bait genome to ensure that all reads carrying germline mutations are captured; restore the alternative alleles of heterozygous sites in the reads aligned to the bait genome to the reference alleles; re-align the corrected reads to the homozygous genome. Step 5: Output standardized BAM files: Generate a pair of corrected BAM files for direct use by subsequent arbitrary somatic mutation detection tools.

2. The method according to claim 1, characterized in that, The steps for detecting germline mutations include: matching normal tissue, which is a healthy tissue sample from the same patient, to distinguish between germline mutations and somatic mutations.

3. The method according to claim 1, characterized in that, The steps for constructing the decoy genome include: retrieving the patient germline mutation information from step 1, replacing the reference allele with the alternative allele at the patient-specific heterozygous sites in the universal reference genome, integrating all the patient's germline variation information, and then calibrating the sequence to obtain the decoy genome.

4. The method according to claim 1, characterized in that, The steps for constructing a homozygous genome include: based on the patient germline mutation information in step 1, replacing the patient homozygous bases in the universal reference genome, retaining the original bases at heterozygous sites, and obtaining a homozygous genome containing only patient homozygous variants after sequence calibration.

5. The method according to claim 1, characterized in that, The sequencing reads correction and realignment steps include mapping the raw sequencing reads of tumor and normal samples to the decoy genome to capture reads carrying heterozygous variants and replacing the bases at heterozygous sites with bases from the universal reference genome.

6. The method according to claim 1, characterized in that, The steps for outputting a standardized BAM file include: remapping the corrected reads to a homozygous genome containing only homozygous variants, and preserving the original reference genome coordinates in the final BAM file.

7. The method according to claim 1, characterized in that, When processing indels for base deletion or insertion, maintain a coordinate system consistent with the reference genome.

8. A personalized reference genome construction system for somatic mutation detection in cancer genomics, characterized in that: The system includes at least one processor and a memory connected to the processor, wherein: The memory is used to store computer programs or instructions; The processor is configured to acquire and execute computer programs or instructions in memory to cause the system to perform the method as described in any one of claims 1-7.

Images