Methods for identifying chromosomal spatial instability in low-coverage next-generation sequencing data as homologous recombination deficiency

Abstract

A genomic data analyzer can be configured to detect and characterize the presence of genomic instability in a tumor sample with a machine learning model, such as a trained convolutional neural network. The genomic data analyzer can use whole-genome sequencing reads as input data, even at low sequencing coverage in high-throughput sequencing workflows often employed in various clinical oncology settings. The genomic data analyzer can arrange aligned read data coverage from a chromosome arm or a whole chromosome to form a coverage data signal array, which can be in the form of an image. The trained machine learning model can process the coverage data signal array to determine whether chromosomal spatial instability (CSI), such as genomic instability caused by homologous repair or recombination deficiency (HRD), is present in the tumor sample. An indication of the latter can guide the selection of a preferred anti-cancer treatment for the tumor.

Description

Technical Field

[0001] The methods described in this article generally involve genomic analysis, and more specifically, the application of genomic information in the detection and treatment of cancer. Background Technology

[0002] Tumor genomics as a predictor of response to anticancer therapy

[0003] In addition to inherent germline mutations, cancer cells often carry large-scale somatic chromosomal aberrations, such as copy number variations, duplications, or deletions in certain alleles or genomic regions. Some of these variants can lead to the loss of certain genomic functions, such as tumor suppressor mechanisms, particularly homologous recombination repair (HRR or HR), making the cancer more aggressive. Identifying these genomic mutations is crucial to recent advances in personalized cancer drug therapy, as they have been shown to predict responses to certain anticancer therapies in subjects with cell hyperproliferation disorders. Subjects can be humans or animals. Examples of anticancer treatments include alkylating agents, such as, but not limited to, platinum-based chemotherapeutic agents, carboplatin, cisplatin, isopropylplatin, nedaplatin, oxaliplatin, pyridine, nitrogen mustard, chlorambucil, melphalan, cyclophosphamide, ifosfamide, estradiol, carmustine, lomustine, formustine, streptozotocin, busulfan, piperobromane, procarbazine, dacarbazine, thiotepa, temozolomide, and / or other antitumor platinum complex compounds; DNA damaging agents or radiation; anthracyclines, such as, but not limited to, epirubicin or doxorubicin; topoisomerase I inhibitors, such as, but not limited to, camptothecin, topotecan, irinotecan; and / or PARP (poly-ADP-ribose polymerase) inhibitors. Examples of PARP inhibitors that selectively destroy cancer cells by exploiting tumor DNA homology repair defects (HRD) include olaparib, rucaparib, niraparib (MK 4827), and taprazole (BMN-673), which have been approved for certain cancer types in the United States and Europe; other examples of PARP inhibitors include iniparib and CEP. 9722(-14-methoxy-9-[(4-methylpiperazin-1-yl)methyl]-9,19-diazapentacyclo[10.7.0.0^{2,6}.0^{7,11}.0^{13,18}]nonadecan-1(12),2(6),7(11),13(18),14,16-hexen-8,10-dione), 3-aminobenzamide, velapirib, pamiparib or E7016(10-((4-hydroxypiperidin-1-yl)methyl)chromeno[4,3,2-de]phthalazin-3(2H)-one). PARP inhibitors indirectly induce a large number of DNA double-strand breaks (DSBs), which prevents the development of HRD tumor cells, while normal cells can generally repair these breaks through HR (Keung et al., 2019, Journal of Clinical Medicine 8(4), p.435).

[0004] HRD genomic analysis test

[0005] Therefore, a method for classifying DNA samples (such as those from tumors) could aid in the diagnosis of possible cancer types, or, due to genome sequencing and analysis, allow for the selection of the most suitable cancer treatment based on the characteristics of severe DNA copy number aberrations in a patient's tumor sample. In particular, determining whether a cancer is homology repair deficient (HRD) could be considerably helpful in treatment planning. Different genomic mutation signatures have been identified for this purpose over the past decade. European Patent EP2609216B1 discloses the use of Global Chromosomal Aberration Score (GCAS) to predict the outcome of anticancer treatments such as PARP inhibitors, radiotherapy, or chemotherapy including platinum-based drugs. European Patent EP2817630B1 proposes detecting the number of telomere allele imbalance (TAI) events, and selecting platinum-based therapies if this number exceeds a reference value for similar cancers known to be platinum-resistant. European patent EP2859118B1 by Stern, Mani, and Popova of the Curie Institute and INSERM discloses a method for predicting HRD by counting the number of segments spanning at least 3 megabases in each genome, the number of which corresponds to the large-scale transition of chromosome copy number (LST). European patent EP2794907B1 by Abkevich et al. from Myriad Genetics discloses counting the total number of loss of heterozygosity (LOH) regions longer than 11 megabases but shorter than the entire chromosome in at least one pair of human chromosomes and comparing this number with a reference value to predict patient response to various possible cancer treatments; EP2981624B1 discloses the use of LOH, TAI, and LST indicators, while EP3180447B1 discloses combining their values ​​into an HRD score and comparing it with a reference value to predict patient response to various possible cancer treatments. The latter detection method is currently used by Myriad Genetics' myChoice CDx assay to test the HRD in ovarian cancer patients, employing customized hybridization-based targeted enrichment on a number of genes with a median coverage of at least 500X. This HRD score has also recently been shown to be a possible prognostic predictor of progression to high-grade serous ovarian carcinoma (HGSOC) (Takaya et al., “Homologous recombination deficiency status-based classification of high-grade serous ovarian carcinoma”, Nature Research Scientific Reports (2020) 10:2757).

[0006] Similar to the Myriad myChoice assay, Foundation Medicine's commercial assay, FoundationFocus CDxBRCALOH, specifically analyzes the BRCA1 and BRCA2 genes to determine whether ovarian cancer patients will respond to rucapanib (a PARP inhibitor) therapy. The latter assay also uses a hybridization-based custom capture of all coding exons of many genes with a median coverage of 500X, but only uses LOH scores assessed from a combination of genome-wide copy number characteristics and small allele ratios of SNPs.

[0007] HRD Analysis Methods for WGS Data

[0008] More generally, WO2017191074 proposed characterizing the HRD status of tumor DNA samples based on probability scores calculated from the analysis of different genomic alterations, including base substitutions, rearrangements, and imprinted signatures (HRDetect score). As described in "HRDetect is a predictor of BRCA1 and BRCA2 deficiency based on mutational signatures," published online on March 13, 2017, by H. Davies et al., *Nature Medicine*, whole-genome sequencing (WGS) has proven to be a viable method for detecting HRD using the latter's HRDetect score, derived from mutational signatures observed across the entire genome, with a sensitivity of 98.7%. Compared to WGS, the sensitivity of HRDetect decreases significantly when only whole-exome sequencing (WES) is applied, ranging between 46.8% and 73%, depending on the specific tuning of the bioinformatics algorithm.

[0009] In addition to breast and ovarian cancer, as in "Genomic aberration based molecular signatures efficiently characterize homologous recombination deficiency in prostate cancer", Sztupinszki et al. have recently begun to study the use of LOH, TAI, and LST as indicators of HRD signatures for n prostate tumors in WGS and WES data (scarHRD).

[0010] Because the methods described above all rely on calls to SNVs and / or INDELs, they require high coverage depths (typically at least 30X) in next-generation sequencing (NGS) workflows. The high coverage requirements in WGS and large sample sizes significantly increase the analytical costs in clinical practice, including the overhead of wet-lab experiments and dry-lab data processing. Therefore, alternative methods using low-throughput WGS (LP-WGS - 1x to 5x) or even ultra-low-throughput WGS (ULP-WGS - as low as 0.1x) may be more advantageous for clinical oncology. In “ShallowHRD: Detection of Homologous Recombination Deficiency from shallow Whole Genome Sequencing”, Bioinformatics, Apr 21 2020, Eeckhoutte et al., ShallowHRD, a software-based method for characterizing LST status, is described, similar to the methods of Manie, Stern, and Pova, but specifically for shallow WGS data with coverage as low as approximately 1X. Their method simply uses the count of copy number changes within chromosome arms. This count can be estimated as the number of large-scale transitions between adjacent segments of at least 10 Mb after removing segments smaller than 3 Mb from the coverage data signal. According to the authors, fewer than 15 transitions can identify HRD-negative tumors, while more than 19 transitions can identify HRD-positive tumors, with sensitivity and specificity similar to the existing scarHRD method. While this method offers promise for a more cost-effective NGS workflow in clinical oncology, we note that its sensitivity remains lower than that achieved by measuring the LST index on HRDetect and SNP arrays. This is likely due to its inability to integrate LOH and TAI features from existing methods. Furthermore, it can only classify tumors with large and frequent chromosomal alterations sufficient to cause a dozen or so large-scale transitions, without considering how these transitions are distributed both within each chromosome arm and across multiple chromosomes in the human genome.

[0011] Recent research in oncology suggests that highly mutated dysfunctional DNA can cause heterogeneous alterations in different regions of chromosomes. In "Regulation of mitotic recombination between DNA repeats in centromeres," Nucleic Acids Research, 2017, Vol. 45, No. 19, Zafaret et al., and in "The dark side of centromeres, types, causes and consequences of structural abnormalities implicating centromeric DNA," Barra et al., Nature Communications (2018) 9:4340, Barra et al. reported a significant proportion of chromosomal rearrangements and breaks observed in cancer cell lines derived from colorectal and adenocarcinomas. This may be due to the inherent fragility of centromere regions in tumors, which is also common in other cancers. Barra et al. emphasized the lack of models and techniques to definitively analyze highly repetitive sequences in these regions.

[0012] Therefore, there is a need for improved genomic analysis methods that can be easily and cost-effectively deployed in automated NGS workflows within routine oncology clinical practice, while adapting to recent cancer genomics findings, particularly promising discoveries regarding the spatial characteristics of certain chromosomal aberrations in the tumor genome. Specifically, there is a need for improved genomic analysis methods that integrate recent cancer genomics findings and allow for the identification of features, including HRD, that are relevant to diagnosis, prognosis, treatment selection, and patient management. Even better, these methods and techniques should also be applicable at lower coverage levels to reduce the overall analytical costs in routine clinical settings, while maintaining reasonable accuracy in detecting HR events compared to existing technologies. Furthermore, improved genomic analysis methods can facilitate the characterization of HR-deficient tumors by employing machine learning techniques, rather than relying on explicit measurements and thresholds of scalar HRD indicators such as LOH, LST, or TAI indices in existing technologies, thereby enabling easier classification of tumors from a wide variety of cancers in humans and animals using a growing body of clinical data. Summary of the Invention

[0013] A method for determining the homologous recombination defect (HRD) status of a subject's DNA sample is proposed. The method includes: obtaining whole-genome sequencing reads of a set of subject DNA samples to be analyzed; aligning these sequencing reads with a reference genome, wherein the reference genome is divided into multiple bins, each bin belonging to the same genomic region along a chromosome arm of the chromosome to be analyzed; calculating and normalizing the alignment reads along each chromosome arm in each bin to obtain a chromosome arm coverage signal; and arranging the chromosome arm coverage signals into a sequence corresponding to the subject's DNA sample. A coverage data signal array of sample A is input into a trained machine learning model, wherein the model has been trained using a set of samples with known homologous recombination defect states to distinguish coverage data signal arrays from samples with positive homologous recombination defect states and from samples with negative homologous recombination defect states, thereby determining the homologous recombination defect score (HRD score) of the subject's DNA sample, and determining the negative, positive, or indeterminate homologous recombination defect (HRD) state of the subject's DNA sample based on the HRD score of the trained machine learning model. In a possible implementation, this set of sequencing reads may be obtained from whole-genome sequencing, wherein the read depth coverage is at most 30x; or this set of sequencing reads may be obtained from low-throughput whole-genome sequencing, wherein the read depth coverage is at least 0.1X and at most 5X. In a possible implementation, calculating and normalizing the numerical values ​​of the alignment readings in each block along each chromosome arm to obtain the chromosome arm coverage signal may include normalizing the coverage signal for each sample and / or normalizing by GC content to apply GC-offset correction. In a possible implementation, the chromosome arm coverage signal may be arranged into a one-dimensional coverage data signal vector or a two-dimensional coverage data signal image. In a possible implementation, the chromosome arm coverage signal may be arranged into a two-dimensional coverage data signal image by aligning the coverage data signal of each chromosome and the centromere block of each chromosome arm, i.e., the nearest block adjacent to the centromere region of the chromosome arm. In a possible implementation, the machine learning model may be pre-trained using a set of tumor data samples with known homologous recombination defect states as training labels. The training dataset may be augmented with artificial sample data generated by combining chromosome data from the data samples with labels of known homologous recombination defect states. The generated data-augmented samples may be intended to represent the purity-ploidy ratio distribution observed in the real sample dataset.In a possible implementation, the reference genome can be divided into a first group of blocks of up to 100 kbp, and further includes the step of collapsing the 100 kbp blocks into a second group of larger blocks of at least 500 kbp before arranging coverage signals on each chromosome arm. The blocks in the first group can have a uniform size of up to 100 kbp, while the blocks in the second group can have a size of 2.5 to 3.5 Mbp and are obtained by pooling 25 to 35 of the 100 kbp blocks in the first group.

[0014] An in vitro method for determining the homologous recombination defect (HRD) status of a patient's DNA sample is proposed, the method comprising:

[0015] - Provide DNA fragments from patient samples;

[0016] - Construct a library containing the fragments that overlap with a set of chromosomes;

[0017] - Sequencing the library to achieve a whole-genome sequencing coverage of up to 30X, preferably at least 0.1X and up to 5X;

[0018] Furthermore, the HRD status of patient samples is determined based on analysis using a machine learning model trained according to the methods disclosed herein. The patient DNA sample may be cell-free extracellular DNA (cfDNA), fresh frozen tissue (FFT), or formalin-fixed paraffin-embedded (FFPE) sample. The HRD score or HRD status of the patient sample can be a predictor of tumor response to cancer treatment regimens. The cancer treatment regimen can be selected from alkylating agents, platinum-based chemotherapy agents, carboplatin, cisplatin, isopropylplatin, nedaplatin, oxaliplatin, pyridine, nitrogen mustard, chlorambucil, melphalan, cyclophosphamide, ifosfamide, estradiol, carmustine, lomustine, fomomustine, streptozotocin, busulfan, piperobromane, procarbazine, dacarbazine, thiotepa, temozolomide, and / or other antitumor platinum complex compounds, DNA damaging agents, radiotherapy agents, anthracyclines, epirubicin, doxorubicin; topoisomerase I inhibitors, camptothecin, topotecan, irinotecan, PARP (poly-ADP-ribose polymerase) inhibitors, olaparib, rucapranib, nirapanib, taprazole, iniparib, CEP The group consisting of 9722, MK4827, BMN-673, 3-aminobenzamide, verapamil, pamiparib, or E7016.

[0019] A method for selecting cancer patients to be treated with platinum-based chemotherapy agents, DNA damaging agents, anthracyclines, topoisomerase I inhibitors, and PARP inhibitors is proposed, including the step of detecting HRD positivity in tumor patient samples according to the method disclosed herein. In possible embodiments, the patient may have cancer selected from high-grade serous ovarian cancer, prostate cancer, breast cancer, or pancreatic cancer.

[0020] A method for training a machine learning algorithm to determine the homologous recombination defect (HRD) status of a subject's DNA sample is proposed. The method includes inputting a coverage data signal array from samples known to be positive for homologous recombination defect and a coverage data signal array from samples known to be negative for homologous recombination defect into a supervised training algorithm.

[0021] The trained machine learning model can be a random forest model, a neural network model, a deep learning classifier, or a convolutional neural network model. The neural network model can be a convolutional neural network trained to produce single-label binary classifications of positive or negative HRD states, single-label multi-class classifications of positive, negative, or indeterminate HRD states, or scalar HRD scores representing HRD states. The machine learning model can be trained in a semi-supervised mode using a data augmentation set generated by sampling data from a set of real samples sharing the same HRD states, normalized purity, and ploidy ratio on chromosomes.

[0022] A method for characterizing DNA samples is proposed, comprising: obtaining a set of sequencing reads aligned to a reference genome from a patient sample; dividing the base pair positions (bp) of at least two chromosomes of the reference genome into a set of blocks, each block corresponding to a genomic region of up to 20 megabyte pairs (20 Mbp), each block containing coverage data from only a single chromosome arm; estimating the coverage data of each block of the sequenced patient sample from the aligned reads; arranging the coverage data into a multidimensional array comprising chromosomes or chromosome arms along one dimension and a set of blocks along another dimension of the chromosomes or chromosome arms, characterized in that the centromere blocks or telomere blocks of each chromosome or chromosome arm are arranged in a spatial arrangement of the multidimensional array; inputting the multidimensional array into a trained machine learning model; and generating a chromosome spatial instability (CSI) index at the output of the trained machine learning model. In a possible implementation, the CSI index may be an indicator of the high (HRD+) or low (HRD-) probability that the patient sample has a homologous recombination (HR) defect. Patient samples can be tumor samples, and the CSI index of patient samples can be a predictor of tumor response to cancer treatment regimens, including platinum-based chemotherapy agents, DNA damaging agents, anthracyclines, topoisomerase I inhibitors, or PARP inhibitors.

[0023] In one implementation, the machine learning model is trained in a supervised or semi-supervised mode using coverage data of labeled samples with known genomic instability states, such as HRD status (positive or negative).

[0024] In one implementation, the chromosome spatial instability (CSI) metric for genomic instability is based on the difference between a DNA sample in which at least one region of at least one chromosome arm has genomic instability and another DNA sample in which at least one region of at least one chromosome arm does not have genomic instability, which can be learned by a machine learning model.

[0025] In a possible implementation, the first set of blocks can be further collapsed into larger blocks by adjusting the size of the larger blocks relative to the length of each chromosome arm before arranging the multidimensional array, such that the block size remains constant along the chromosome arm. The first set of blocks can have a uniform size of up to 100 kbp, while the collapsed blocks can have a size of at least 500 kbp. Coverage data can be normalized and / or segmented to infer discrete copy values ​​before arranging the multidimensional array. In a possible implementation, empty elements of the multidimensional array can be filled with the value of the nearest block or a predetermined value. In a possible implementation, the trained machine learning model can be a random forest model or a neural network model, such as a convolutional neural network, trained to produce single-label binary classifications of positive or negative CSI states, or single-label multi-class classifications of positive, negative, or indeterminate CSI states, or scalar HRD scores representing CSI states at its output. Patient samples can be tumor samples, and the machine learning model can be trained in supervised or semi-supervised mode using a set of real samples labeled with CSI states according to the target application. In a possible implementation, the CSI state label can be the HRD state of artificial samples generated according to the HRDetect method and / or by sampling chromosomes from a set of real samples sharing the same HRD state and having similar tumor content. The artificial samples can be chosen to recreate the same purity-ploidy distribution as the real sample dataset.

[0026] In one specific implementation, the patient sample can be a tumor sample, and the machine learning model can be trained in a supervised or semi-supervised mode using labeled data with HRD status. The training data can consist of coverage data obtained from clinical samples with HRD status and arranged in a multidimensional array. The training data labels can be obtained from HRD detection methods, such as the HRDetect method.

[0027] In one particular implementation, data augmentation strategies that take into account biases introduced by sample-specific attributes (such as purity and ploidy) can be used to increase the diversity and number of multidimensional arrays used to train machine learning algorithms. Attached Figure Description

[0028] Figure 1 This refers to a next-generation sequencing system according to certain embodiments of the present disclosure.

[0029] Figure 2 This describes a genome analysis workflow according to certain embodiments of the present disclosure.

[0030] Figure 3The illustration schematically shows a possible implementation of a low-coverage data preparation preprocessing workflow, including collapsing and smoothing small coverage blocks into larger blocks to better demonstrate more specific spatial events observable along different chromosome arms.

[0031] Figure 4 An example of normalized readout coverage data signal along 22 non-sex chromosomes in a human DNA sample is depicted, which can be initially measured from the first set of small blocks in the coverage data signal of a low-pass whole genome sequencing experiment.

[0032] Figure 5 An example of a normalized reading coverage data signal is plotted, which can be derived from the coverage data signal using a second set of larger blocks. Figure 4 It was calculated in the middle.

[0033] Figure 6 The complete set of human genome chromosomes is shown (see x-axis) arranged to align their centromere regions, along with the lengths of their respective p-arms and q-arms (y-axis).

[0034] Figure 7 The illustrations illustrate some embodiments of the proposed coverage data preparation method, showing normalized coverage data that is spatially rearranged into a two-dimensional matrix of chromosomes arranged in rows and coverage blocks arranged in columns, such that centromere blocks are vertically arranged in a single column prior to chromosome spatial instability analysis, wherein a) is an HRD-negative tumor DNA sample and b) is an HRD-positive tumor DNA sample.

[0035] Figure 8 The internal data workflow architecture of a possible chromosome spatial instability analyzer (e.g., for HRD state analysis) is illustrated schematically.

[0036] Figure 9 Possible architectures of convolutional neural networks based on certain implementations of the proposed CSI analyzer data processing module are shown.

[0037] Figure 10a )and Figure 10b The image shows normalized coverage data for two samples from the 560 Breast Cancer Genome Database.

[0038] Figure 11 Normalized coverage data for 22 non-sex chromosomes are shown for 202 BAM file samples from the 560 Breast Cancer Genome Database.

[0039] Figure 12 The existing HRD score method is compared with the HRDetect score of the proposed CSI index classifier.

[0040] Figure 13 The results of HRD scores, proposed CSI index scores, and HRDetect scores for a series of 69 test samples with different BRCA missing (possibly related to HR missing) were plotted.

[0041] Figure 14 Exemplary alignments of coverage signal images for each chromosome arm (p-arm and q-arm) from at least two exemplary chromosomes (chromosome A and chromosome B) relative to at least one block adjacent to the centromere region or one block adjacent to the telomere region are shown. Image A is a coverage signal image for each chromosome arm (p-arm and q-arm) from at least two exemplary analyzed chromosomes (chromosome A and chromosome B), with dotted rectangles representing blocks adjacent to the telomere region and vertical striped rectangles representing blocks adjacent to the centromere region; Image B is an alignment of the coverage signal images with respect to blocks adjacent to the centromere region; Image C is an alignment of the coverage signal images with respect to blocks adjacent to the telomere region; Image D is an alignment of the coverage signal images with respect to blocks adjacent to the centromere region; Image E is an alignment of the coverage signal images with respect to blocks adjacent to the telomere region.

[0042] Figure 15 shows exemplary output data obtained from the combined lpWGS and directional capture as described in Example 2. Figure 15A The raw coverage obtained from NGS data generated from a workflow combining low-pass WGS (lpWGS) and targeted sequencing is shown, calculated at 100kB intervals. For illustrative purposes, raw coverage of genomic regions enriched in targeted experiments is represented by empty markers, while lpWGS regions are represented by solid markers. Figure 15B The image shows the area outside the target (lpWGS only) (left, right). Figure 15B B1) and target (capture enrichment) region (right, Figure 15B The arrangement of the mapping readings of B2). Figure 15C The normalized coverage is shown only for the lpWGS region (left figure). Figure 15C C1), and allele portions of variants identified in the capture enrichment region (right figure, Figure 15C (C2).

[0043] Figure 16The preparation of the artificial sample as described in Example 3 is illustrated. Image A shows the normalized coverage (Y-axis) of the genome (X-axis, normalized coverage map) across three randomly selected original tumor samples used to generate the artificial sample. The tumor content (or purity, T), ploidy (P), and purity-ploidy ratio (PPR) for each sample are indicated above the corresponding normalized coverage map of the sample. Image B shows the normalized coverage map of the sample after in-silico dilution to match the lowest PPR observed in the original samples of this group. Image C shows the normalized coverage map of the artificial sample assembled by combining the data from the coverage maps of the diluted samples. The colors of the normalized coverage points for the original samples, diluted samples, and artificial samples are preserved. Detailed Implementation

[0044] This disclosure is based in at least part on the following findings: the machine learning-trained analyzer of this disclosure, designed to process sequencing data coverage of tumor samples related to basic chromosome arrangement, can extract indicators of general chromosomal spatial instability (CSI), particularly indicators of homology repair deficiency (HRD) status in tumor samples.

[0045] The proposed methods and systems will now be described with reference to more detailed embodiments. However, the proposed methods and systems may be embodied in different forms and should not be construed as limited to the embodiments described herein. Rather, these embodiments are provided to make this disclosure thorough and complete, and to fully express its scope to those skilled in the art.

[0046] definition

[0047] A "DNA sample" refers to a nucleic acid sample derived from an organism, such as a solid tumor or a liquid. This organism can be human, animal, plant, fungus, or microorganism. Nucleic acids can be found in solid samples, such as formalin-fixed paraffin-embedded (FFPE) samples. Alternatively, nucleic acids can be found in limited quantities or low concentrations, such as circulating tumor DNA in blood or plasma.

[0048] A "DNA fragment" refers to a short DNA segment resulting from the fragmentation of high-molecular-weight DNA. Fragmentation can occur naturally within the sample organism or be artificially generated through DNA fragmentation methods applied to the DNA sample, such as mechanical shearing, sonication, enzymatic fragmentation, and other methods. After fragmentation, DNA fragments can undergo end repair to ensure that each molecule has a blunt end. To improve ligation efficiency, adenine can be added to each 3' blunt end of the fragmented DNA, enabling the DNA fragment to ligate to adapters with complementary dT-protrusions.

[0049] "DNA products" refers to engineered DNA fragments produced by manipulating, extending, ligating, doubling, amplifying, replicating, editing, and / or cutting DNA fragments to adapt them to next-generation sequencing workflows.

[0050] "DNA adapter products" refer to DNA products created by connecting DNA fragments to DNA adapters to make them compatible with next-generation sequencing workflows.

[0051] A “DNA library” refers to a collection of DNA products or DNA adapter products that make DNA fragments compatible with next-generation sequencing workflows.

[0052] A "pool" refers to multiple DNA samples (e.g., 48, 96, or more) from the same or different organisms that can be multiplexed into a single high-throughput sequencing analysis. Each sample can be identified in the pool by a unique sample barcode.

[0053] A "nucleotide sequence" or "polynucleotide sequence" refers to any polymer or oligomer of nucleotides, such as cytosine (represented by the letter C in the sequence string), thymine (represented by the letter T in the sequence string), adenine (represented by the letter A in the sequence string), guanine (represented by the letter G in the sequence string), and uracil (represented by the letter U in the sequence string). It can be DNA or RNA, or a combination thereof. It can exist permanently or temporarily in single-stranded or double-stranded form. Unless otherwise specified, nucleic acid sequences are written from left to right in a 5' to 3' orientation.

[0054] "Amplification" refers to a polynucleotide amplification reaction to produce multiple polynucleotide sequences that are replicated from one or more parent sequences. Amplification can be achieved through various methods, such as polymerase chain reaction (PCR), linear polymerase chain reaction, nucleic acid sequence-based amplification, rolling cycle amplification, and other methods.

[0055] "Sequencing" refers to reading the sequence of nucleotides from a DNA library to produce a set of sequencing reads that can be processed by a bioinformatics computer in a bioinformatics workflow. High-throughput sequencing (HTS) or next-generation sequencing (NGS) refers to the real-time parallel sequencing of multiple sequences, typically between 50 and several thousand base pairs per sequence. Typical NGS technologies include those from Illumina, Ion Torrent Systems, Oxford Nanopore Technologies, Complete Genomics, Pacific Biosciences, BGI, and others. Depending on the specific technology, NGS sequencing may require sample preparation with sequencing adapters or primers to facilitate further sequencing steps, as well as amplification steps to sequence multiple instances of a single parent molecule, such as PCR amplification before delivery to a flow cell in the case of synthetic sequencing.

[0056] “Sequencing depth,” “sequencing coverage,” or “depth of sequencing” refers to the number of times a genome is sequenced. In target enrichment workflows, only a small fraction of the entire genome is sequenced, so increasing sequencing depth may be justifiable without incurring significant data storage and processing overhead. In some genomic analysis applications that do not require high resolution along the genome, such as in detecting copy number changes, low-pass (LP) coverage (1x–10x) or even ultra-low-pass (ULP) coverage (<1x – not all locations are sequenced) may be more efficient in terms of IT infrastructure costs, but these workflows require more sophisticated bioinformatics methods and techniques to handle less reliable data input from sequencers and alignment machines. Furthermore, in addition to the higher costs associated with data storage and processing, the operational costs of experimental NGS runs need to be optimized by balancing coverage depth with the number of samples that might be tested in parallel in a routine clinical workflow – i.e., loading the sequencer with sequencing samples. Indeed, the total number of reads a next-generation sequencer can produce in a single experiment (i.e., in a given run) remains limited. Lower coverage means fewer reads per sample for genomic analysis, but more samples can be reused in a next-generation sequencing run.

[0057] "Alignment," "alignment," or "aligner" refers to the base-by-base mapping and alignment of sequencing reads with a reference genome sequence within a bioinformatics workflow, depending on the application. For example, in targeted enrichment applications, sequencing reads are expected to map to specific target genomic regions based on the hybridization capture probes used during experimental amplification. Alignment can be specifically searched relative to the corresponding sequence, defined by genomic coordinates such as chromosome number, start position, and end position in the reference genome. As is known in bioinformatics practice, in some implementations, the "alignment" method described herein may also include certain preprocessing steps to facilitate the mapping of sequencing reads and / or remove irrelevant data from the reads, such as by removing unpaired reads, and / or by pruning adapter sequences as ends of reads, and / or other read preprocessing filtering methods. Exemplary bioinformatics data representation methods with different coordinate systems (absolute or relative position indexes, 0-based or 1-based, etc.) include BED format, GTF format, GFF format, SAM format, BAM format, VCF format, BCF format, Wiggle format, GenomicRanges format, BLAST format, GenBank / EMBL feature table format, and others.

[0058] "Coverage," or "sequence read coverage," refers to the number of sequencing reads aligned with a genomic location or a set of genomic locations. Generally, genomic regions with high coverage are associated with higher reliability in downstream genome characterization, especially when calling variants.

[0059] "Bin," "genome fragment," "region," "genome section," or "chromosomal section" refers to a contiguous region of interest within the genome. Because such regions may include variants, a bin generally refers to a location or region of the genome, rather than a fixed DNA sequence. In bioinformatics methods and processes, a bin can be identified by its genomic coordinates along the start and end points of a reference genome, and the length of the bin can be measured in bases (b, kb, Mb) or base pairs (bp, kbp, Mbp) from the start to the end coordinates. Generally, a bin can correspond to an entire chromosome, a segment of a chromosome, a segment of the reference genome, a portion of multiple chromosomes, multiple chromosomes, portions of multiple chromosomes, and / or combinations thereof. Preferably, a bin is a portion of a chromosome obtained by dividing the reference genome into genomic sections (e.g., by size, fragment, contiguous region, contiguous region of arbitrarily defined size, etc.). Genome sections can be selected, classified, filtered, and / or removed from consideration using any suitable criteria known in bioinformatics. Along the genome, bins can have the same, uniform length or different, variable lengths.

[0060] "Coverage block count" or "coverage count" or "count" refers to the number of sequencing reads or end-of-line reads (representing DNA fragments) mapped to or partially overlapping with blocks. Counts can be derived from mapping or aligning some or all of the raw sequence reads and / or preprocessed sequence reads to (i.e., the relevant) genomic portion. According to various bioinformatics methods known in the art, some sets of reads can be weighted, removed, filtered, normalized, adjusted, collapsed, merged, added, and / or subtracted, or by a combination thereof, on a read-by-read or end-of-line basis before counting. In some embodiments, the count may be associated with uncertainty or quality values. In some embodiments, reads or end-of-line reads may not be entirely contained within a single block but overlap within two adjacent blocks; in this case, the reads may be counted within the block with the greatest overlap.

[0061] "Coverage data signal" or "coverage signal" or "coverage data signal array" refers to a collection of blocks associated with their respective counts, which can be arranged as a one-dimensional vector (thus obtaining a coverage signal vector or coverage signal array or one-dimensional heatmap), a two-dimensional matrix (thus obtaining a coverage signal matrix or coverage signal image or two-dimensional heatmap), or any suitable topology. Coverage data signals can be selectively normalized to eliminate technical bias. Coverage data signal arrays can be arranged in a two-dimensional representation, referred to herein as a coverage data signal image or a coverage data signal multidimensional array.

[0062] A "chromosomal arm" refers to either of the two arms of a chromosome, joined together by the centromere. The p-arm is the shortest arm, and the q-arm is the longest. Each arm terminates at a telomere. In mesocentric chromosomes, the p-arm and q-arm are similar in size. In submesocentric chromosomes, the p-arm is shorter than the q-arm. In acrocentric chromosomes, the p-arm is very short. In telomere chromosomes, the p-arm is absent or too short to be visible when examining the chromosome. Normal human cells do not carry telomere chromosomes; however, they can be found in some tumor cells.

[0063] The "centromere region" or "centromere area" of a chromosome refers to the region immediately adjacent to the centromere of the chromosome.

[0064] The “centromere block” in chromosome coverage data refers to the coverage block that is closest to the centromere region in the set of coverage blocks associated with chromosome arms. At most two centromere blocks can be associated with any given human chromosome's genome coverage data: one to the left of the centromere and one to the right along the genome coordinate system of the sequencing data.

[0065] The "telomere region" of a chromosome refers to the area next to the telomeres of a chromosome.

[0066] A “telomere block” in chromosome coverage data refers to the coverage block in the set of coverage blocks associated with that chromosome that is closest to either of the two telomere regions at the beginning or end of that chromosome. For any given human chromosome, there can be at most two independent, distant telomere blocks associated with genome coverage data, one on the p-arm and the other on the q-arm.

[0067] Chromosomal instability, or CIN, refers to a type of genomic instability in which chromosomes are unstable, resulting in the duplication or deletion of entire chromosomes or parts of chromosomes. More specifically, CIN refers to an increased rate of increase or loss of entire chromosomes or parts of chromosomes. This type of genomic alteration can occur particularly in tumor cells, involving the increase or loss of entire chromosomes or structural aberrations, such as severe chromosomal rearrangements.

[0068] Chromosomal spatial instability, or "CSI," refers to chromosomal instability where the rate of increase or decrease in chromosomes can be characterized by the spatial distribution of events across the genome. Understandably, CSI information can be obtained from the output of a trained machine learning model, either directly as a CSI score calculated by the machine learning model for a given phenotype (e.g., a scalar number), or as a CSI state, which can be labeled as a positive "CSI+", a negative "CSI-", or a potentially uncertain or indeterminate value "Undetermined CSI". The positive or negative status of CSI can also be inferred from the CSI score.

[0069] "Homologous recombination defective state" or "HRD state" refers to the classification of the homologous recombination pathway and is associated with any cellular state / event that leads to a defect in the homologous recombination pathway. HRD state can be classified as positive (HRD+), in which the homologous recombination pathway is defective; or negative (HRD-), in which the homologous recombination pathway is not defective; or other undetermined states (HRD indeterminate, HRD unknown).

[0070] A "machine learning model" refers to a data model or data classifier trained using supervised, semi-supervised, or unsupervised learning techniques known in data science, rather than an explicit statistical model. Data input can be represented as a one-dimensional signal (vector), a two-dimensional signal (matrix), or more generally, a multidimensional array of signals (e.g., a tensor, or an RGB color image represented as a 3*2D signal with red, green, and blue color decomposition planes—three matrices), and / or combinations thereof. A multidimensional array is mathematically defined as a data structure arranged along at least two dimensions, each recording a value greater than 1.

[0071] In the case of deep learning classifiers, data input is further processed through a series of data processing layers to implicitly capture hidden data structures, data signatures, and basic patterns. Due to the use of multiple data processing layers, deep learning facilitates the generalization of automated data processing, enabling a variety of complex pattern detection and data analysis tasks. Machine learning models can be trained within supervised, semi-supervised, or unsupervised learning frameworks. Within a supervised learning framework, the model learns a function based on instances of input and matching output to map the output from the input dataset. Examples of machine learning models used for supervised learning include: Support Vector Machines (SVM), regression analysis, linear regression, logistic regression, and Naive Bayes (Naive Bayes). Unsupervised learning frameworks include Bayes, linear discriminant analysis, decision trees, k-nearest neighbors, random forests, artificial neural networks (ANNs) such as convolutional neural networks (CNNs), recurrent neural networks (RNNs), fully connected neural networks, long short-term memory (LSTM) models, and others; and / or combinations thereof. Models trained within unsupervised learning frameworks infer functions that determine the hidden structure of a dataset without requiring prior knowledge of the data. Examples of unsupervised machine learning models known in the art include: clustering, such as k-means clustering, mixture model clustering, hierarchical clustering; anomaly detection methods; principal component analysis (PCA), independent component analysis (ICA), stochastic neighbor embeddings of t-shaped distributions (t-SNE); generative models; and / or unsupervised neural networks; autoencoders; and / or combinations thereof. Semi-supervised learning (SSL) is a machine learning framework in which models can be trained using both labeled and unlabeled data. Data augmentation methods can be used to generate artificial data samples from sparse groups of real data samples and increase the quantity and diversity of data used for model training. Unlabeled data, when used with a small amount of labeled data, can produce a considerable improvement in learning accuracy compared to other frameworks. This approach is particularly interesting when only a portion of the available data is labeled.

[0072] A "convolutional neural network" or "CNN" refers to a machine learning model that uses multiple data processing layers, called convolutional layers, to represent input data in a way best suited for solving classification or regression tasks. During training, optimization algorithms known in the art, such as backpropagation, are used to perform stochastic gradient descent, optimizing the weight parameters for each CNN layer. At runtime, the resulting trained CNN can process input data very efficiently, for example, classifying it as the correct data output label in the case of learning a classification task, and minimizing false positives and false negatives. Convolutional neural networks can also be combined with recurrent neural networks to produce deep learning classifiers.

[0073] Genome analysis system

[0074] Now refer to Figure 1 Exemplary genome analysis systems and workflows are described in further detail. It will be apparent to those skilled in the art of DNA analysis that a genome analysis workflow comprises: preliminary experimental steps performed in a laboratory (also known as a “wet laboratory”) to generate DNA analysis data, such as raw sequencing reads in a next-generation sequencing workflow; and subsequent data processing steps using a bioinformatics system (also known as a “dry laboratory”) to further determine information of interest to the end user, such as detailed identification of DNA variants and associated annotations. Various implementations of the DNA analysis workflow are possible depending on the specific application, laboratory setup, and bioinformatics platform. Figure 1 Examples of NGS systems are described, including a wet laboratory system in which DNA samples are first prepared experimentally using a DNA library preparation protocol 100, which generates, adapts to sequencing, and amplifies DNA fragments for processing by an NGS sequencer 110. In the next-generation sequencing workflow, the resulting DNA analysis data can be generated as a data file of raw sequencing reads in FASTQ format. This workflow can then further include a dry laboratory genomic data analysis system 120, which takes raw sequencing reads from a DNA sample pool prepared according to a proposed method as input and applies a series of data processing steps to characterize certain genomic features of the input sample. An exemplary genomic data analysis system 120 is the Sophia Data-Driven Medicine Platform (Sophia DDM), which was used by over 1000 hospitals worldwide in 2020 for the automated identification and characterization of genomic variants and reporting to end users, but other systems may also be used. For example, different detailed possible implementations of the data processing steps that can be applied by the genomic data analysis system 120 for genomic variant analysis are described in International PCT patent application WO2017 / 220508, but other implementations are also possible.

[0075] like Figure 1As shown, the genome data analyzer 120 may include a sequence alignment module 121 that compares raw NGS sequencing data with a reference genome, such as a human genome in medical applications or an animal genome in veterinary applications. In conventional genome data analysis systems, the resulting alignment data may be further filtered and analyzed by a variant invocation module (not shown) to retrieve variant information, such as SNPs and INDEL polymorphisms. The variant invocation module may be configured to execute different variant invocation algorithms. The detected variant information may then be output by the genome data analysis module 120 as a genome variant report for further processing by the end user, such as with visualization tools, and / or by a further variant annotation processing module (not shown). In a possible implementation, the genome data analysis system 120 may include automated data processing modules, such as a coverage data preparation module 122 to prepare spatially arranged coverage data signals, and a chromosome spatial instability (CSI) analysis module 123 to analyze the spatially arranged coverage data signals and obtain CSI information, which may then be reported to the end user, for example, with visualization tools, or to another downstream process (not shown). In a possible implementation, the chromosome spatial instability (CSI) analysis module 123 may be adapted to analyze the coverage data signal of spatial arrangement and derive a CSI score, such as an HRD score, or a CSI state, such as an HRD state, from the CSI information. The HRD score, HRD state, and / or CSI information may then be reported to the end user, for example, using a visualization tool, or to another downstream process (not shown).

[0076] Data processing workflow Procedure

[0077] The Genome Data Analyzer 120 can process sequencing data by employing and combining different data processing methods to generate genome data analysis reports.

[0078] The sequence alignment module 121 can be configured to execute different alignment algorithms. Standard raw data alignment algorithms, such as Bowtie2 or BWA optimized for rapid processing of large volumes of genomic sequencing reads, can be used, but other implementations are also possible. The alignment results can be represented as one or more BAM or SAM format files, as is known to those skilled in bioinformatics, but other formats, such as compressed formats or formats optimized for sequence-preserving encryption, can also be used, depending on the requirements of the genomic data analyzer 120 for storage optimization and / or genomic data privacy enforcement.

[0079] Figure 2A possible workflow of an automated data processing module is illustrated to generate a CSI report from the analysis of alignment files (such as BAM files) from bioinformatics sources. In a possible implementation, the coverage data preparation module 122 can be configured to generate 200 coverage data signals from the aligned readings as the raw number of readings aligned to each genomic block or the opposite end readings, and arrange the resulting data 210 into a spatially arranged coverage data signal. The spatially arranged coverage data can be represented as a matrix, i.e., a set of matrices, or as coverage multidimensional array data (tensors) to facilitate the processing of two-dimensional images, i.e., multiplanar two-dimensional images, or multidimensional array (tensor) input signals by the CSI analyzer 123. The chromosome spatial instability (CSI) analyzer module 123 can then further analyze, classify, and / or categorize 220 the spatially arranged coverage data information to report 230 to the end user whether the data exhibits certain characteristics that are features of possible chromosome spatial instability (CSI) in the tumor sample. In a possible implementation, the Chromosomal Spatial Instability (CSI) analysis module 123 can generate a CSI score. In another possible implementation, the Chromosomal Spatial Instability (CSI) analysis module 123 can generate a CSI status, such as a negative status, a positive status, or an indeterminate status. In a possible implementation, the CSI status can be inferred from the CSI score, but other implementations are also possible. In a preferred implementation, the CSI can characterize homology repair deficiency (HRD) of a tumor sample as an HRD score or HRD status, but it can also characterize other genomic features associated with chromosomal spatial instability, such as HRD caused by BRCA1 or BRCA2 deficiency. This report can then aid in tumor genomic diagnosis of cancer subtypes and / or their prognosis. The report also helps in selecting medical approaches specifically for the analyzed tumor, such as using certain cancer treatments that exploit HRR pathway defects to reduce or stop tumor cell proliferation.

[0080] The genome data analyzer 120 may be a computer system or part of a computer system, including a central processing unit (CPU, "processor" or "computer processor" herein), memory such as RAM and storage units such as hard disks, and a communication interface for communicating with other computer systems via a communication network, such as the Internet or a local area network. Examples of genome data analyzer computing systems, environments, and / or configurations include, but are not limited to: personal computer systems, server computer systems, thin clients, thick clients, handheld or laptop devices, multiprocessor systems, microprocessor-based systems, set-top boxes, programmable consumer electronics, network PCs, microcomputer systems, mainframe computer systems, graphics processing units (GPUs), etc. In some embodiments, the computer system may include one or more computer servers that operate alongside numerous other general-purpose or special-purpose computing systems and enable distributed computing, such as cloud computing in a genome data farm. In some embodiments, the genome data analyzer 120 may be integrated into a large-scale parallel system. In some embodiments, the genome data analyzer 120 may be directly integrated into a next-generation sequencing system.

[0081] The 120 genomic data analyzer computer system can be modified in general, based on executable instructions (such as program modules executed by the computer system). Generally, program modules can include routines, programs, objects, components, logic, data structures, etc., that perform specific tasks or implement specific abstract data types. As those skilled in computer programming know, program modules can utilize local operating system and / or file system functionalities; standalone applications; browser or application plugins, applets, etc.; commercial or open-source libraries and / or library tools, as they can be programmed in Python, Biopython, C / C++, or other programming languages; and custom scripts, such as Perl or Bioperl scripts.

[0082] Instructions can be executed in a distributed cloud computing environment, where tasks are performed by remote processing devices connected via a communication network. In a distributed cloud computing environment, program modules can reside on storage media, including memory storage devices, in both local and remote computer systems.

[0083] Therefore, it can be understood that the method described here is a computer-implemented method.

[0084] Coverage data preparation - blockization

[0085] In a possible implementation, the reference genome coordinates are first divided into a set of P blocks. Alignment reads from a BAM or SAM file (after sequencing read alignment) are used as input. Bioinformatics methods known in the art are used to count the number of DNA fragments mapped to each block, generating 200 coverage data signals along the P blocks. It will be apparent to those skilled in bioinformatics that in some implementations, some blocks may be ignored, filtered, or merged with other blocks to improve the quality of the input data.

[0086] In a possible simple implementation, blocks along the reference genome can use a uniform size, but heterogeneous sizes can also be chosen. The block size can be selected to preserve sufficient spatial detail for the analysis of coverage data signals for each chromosome, while facilitating coverage data signal analysis using the machine learning analysis methods described herein. Preferably, the block size can be from 300 bp to 20 Mpb, but other implementations are also possible.

[0087] In the case of the reference human genome, chromosomes are ordered from the largest chromosome (chr1) to the smallest chromosome (chr22), possibly excluding sex chromosomes. Human chromosomes are numbered in descending order of length, from chr1 (249 Mbp) to chr22 (51 Mbp). For example, by dividing the reference genome into 3 Mbp blocks, a one-dimensional array (vector) can be obtained, where the longest chromosome 1 has 84 blocks and the shortest chromosome 22 has 4 blocks. Table 1 shows the theoretical starting and ending points of the centromeres, the total chromosome length and its type, and the length (Mbp) of the p-arm (short arm) and q-arm (long arm) of each chromosome in the human genome, referring to the GRCh37 reference genome (hs37d5 version) in the coordinate system of the reference genome used for read alignment. Please note that with current high-throughput sequencing technologies, not all human genome blocks can contain alignment reads from BAM files. This is because it is difficult to accurately map short NGS reads to certain repetitive sequences or problematic regions of the human genome, such as the p arms of acrocentromeres 13, 14, 15, 21, 22, or Y.

[0088]

[0089] Table 1. Centromere positions, chromosome length (bp), and start and end positions of chromosome arm length (Mbp) of human chromosomes as represented in the GRCh37 reference genome coordinate system (hs37d5 version).

[0090] It is obvious to those skilled in bioinformatics that, in possible implementations, the block size and location along the whole-genome sequencing data can also be adjusted based on the actual coverage of the BAM file. Generally, we can consider each block in a set of P blocks to be defined by two values ​​corresponding to the start and end positions of the block {block_start, block_end} in the reference coordinate system, where 1 <= block_index <= P. Therefore, the block size can be variable, i.e., block_size = block_end - block_start. In possible implementations, the block size can be 0.5 Mbp to 1.5 Mbp, or from 1.5 Mbp to 2.5 Mbp, or from 2.5 Mbp to 3.5 Mbp, or from 3.5 Mbp to 4.5 Mbp, or from 4.5 Mbp to 5.5 Mbp, or from 5.5 Mbp to 6.5 Mbp, or from 6.5 Mbp to 7.5 Mbp, or from 7.5 Mbp to 8.5 Mbp, or from 8.5 Mbp to 9.5 Mbp. In another possible implementation, the block size can be selected to obtain a fixed number of blocks per chromosome or per chromosome arm. Other implementations are also possible.

[0091] Coverage signals for each chromosome or chromosome arm can then be obtained by counting the number of alignment reads in each block of a fixed set of blocks. In a possible implementation, the resulting array of high-resolution coverage data signals for each chromosome or chromosome arm can serve as high-dimensional data input for downstream CSI analysis. In a possible alternative implementation, a first size of 100 kbp can be used to generate a first set of high-resolution coverage blocks. The resulting high-resolution coverage data signals can be further collapsed into a second set of larger blocks of 1 to 20 Mbp, for example, by combining approximately 30 large initial blocks of 100 kbp to generate larger blocks (2.5 Mbp to 3.5 Mbp), thereby reducing the dimensionality of the data for downstream CSI analysis. Various methods known to those skilled in bioinformatics can be employed for this purpose, such as smoothing, using the median or mean to reduce dimensionality, pooling, sampling, and other methods. In a preferred embodiment, the collapse uniquely assigns each block to a specific chromosome arm without mixing coverage data from different chromosome arms into the blocks. Then, the coverage data is clearly separated from one chromosome arm to the next, which may help in analyzing chromosome spatial instability events at the chromosome arm level.

[0092] In a possible implementation, collapse may include adjusting the size of the collapsed blocks according to the length of each chromosome arm to reduce boundary effects at chromosome arm boundaries. For example, if the length of a chromosome arm cannot be divided by the initial block size under consideration, the resulting collapsed blocks will typically be shorter than others, causing boundary effects (e.g., centromere or telomere blocks at either end of a chromosome arm may not include the entire 3 Mbp region, potentially leading to distorted coverage of these regions, which could be critical for analyzing CSI events). A strategy to reduce this effect is to adjust the block size relative to the chromosome arm length so that the considered chromosome arm length can be divided by the block size with negligible remainder. This ensures that all collapsed blocks on the considered chromosome arm have the same actual size.

[0093] In a possible further implementation, collapse may include individually adjusting the start or end position of the collapsed blocks along the chromosome arm to adjust the actual size so that the sum of the variable sizes of the collapsed blocks is still equal to the length of the chromosome arm, while being more or less uniformly distributed along the chromosome arm.

[0094] Figure 3 Chromosome 300 is schematically illustrated, including a central position 320 corresponding to the centromere, a starting position 310 corresponding to the p-arm telomere, and a terminal position 330 corresponding to the q-arm telomere. This chromosome includes a first telomere region 315 starting from the starting position 310, a centromere region 325 surrounding the centromere, and a second telomere region 335 ending at the terminal position 330. It will be apparent to those skilled in bioinformatics that a portion 351 of the centromere region 325 has no corresponding reference sequence in the existing reference genome covering this region 325 and therefore cannot be mapped. In possible implementations, high-resolution blocks corresponding to regions with limited coverage or low quality can also be removed based on genomic knowledge of the corresponding region (e.g., Figure 3 (Region 370 in the BAM file). In a preferred embodiment, a different set of first high-resolution coverage blocks 350, 352 can be constructed separately for each chromosome arm from the BAM file. Each set of high-resolution blocks can then be further collapsed into a second set of larger-range blocks 360, 362 for each arm. The blocks can also span shorter or larger genomic regions based on genomic knowledge of the corresponding region (e.g., a 2.5 Mbp block 371, then a 3.5 Mbp block 372, averaging to a 3 Mbp block), and / or optimize the blocks at the ends of the two chromosome arms (telomere blocks 372, 375, which are the closest blocks to telomeres 310, 330; or centromere blocks 373, 374, which are the closest blocks to centromere 320).

[0095] For example, initial 100kbp blocks from each chromosome arm can be smoothed independently to obtain 3Mbp blocks (i.e., a cluster of 30 100kbp blocks). Since the total number of blocks per chromosome arm is unlikely to be divisible by 30, the final smoothed block of a chromosome arm typically contains fewer than 30 of the original blocks. Retaining these regions could give them too much importance, while deleting them would result in a loss of too much information. To reduce boundary effects, larger, variable-sized blocks can be collapsed independently for each chromosome arm. In possible implementations, the block size that results in the minimum amount of incomplete blocks and is closest to the target size can be selected. In possible implementations, if the last block is incomplete, it can be rejected (because it can now be ignored), but other implementations are also possible. Column 3 of Table 2 shows an exemplary distribution of proposed collapsed block sizes for each chromosome arm, with the initial set of 100kbp blocks collapsing along the 2*22 chromosome arms into larger blocks of the target 3Mbp size. Column 4 shows the size of the resulting rejected boundary block, compared to the number of remaining high-resolution blocks by default in Column 2. When all arms have a fixed size of 3 Mpb, the last collapsed block of the arm boundary will be formed.

[0096]

[0097]

[0098] Table 2. Exemplary block collapse distributions, from 100kbp blocks to target collapse block sizes of 3MBp.

[0099] Coverage data preparation - normalization

[0100] It will be apparent to those skilled in the art of bioinformatics that, in some implementations, the raw coverage data signal can be further normalized, filtered, or smoothed using methods known in the art. In possible implementations, normalization may include sample-by-sample normalization by dividing the coverage data signal by the average coverage signal for the entire sample to eliminate biases that may cause differences in coverage between samples, as associated with the sample sequencing experiment. Normalization may also include applying GC-bias correction by normalizing the GC content. Other implementations are also possible, such as block-by-block normalization, linear and nonlinear least squares regression, GC LOESS (GC normalization), LOWESS, PERUN (sample-by-sample normalization), RM, GCRM, cQn, and / or combinations thereof.

[0101] In a possible further implementation, the normalized data can be further segmented into sub-regions within blocks to identify homogeneous regions in the coverage data signal, which may correspond to different copy numbers (typically two per chromosome). Coverage segmentation can infer the copy number of discrete segments relative to other genomic regions from the coverage curve of a given sample. For this purpose, the coverage curve is decomposed into segments (i.e., parts of the genome) where coverage is considered constant. Each segment defined in the previous step can then be associated with discrete levels defined by the amplitude of the coverage signal. At this stage, copy number variation (CNV) events can be detected, but the absolute copy number associated with each coverage level remains unknown.

[0102] When the tumor content and ploidy of a sample are known or can be inferred from the data, the discrete levels determined by segmentation can be mapped to integers reflecting the copy number present in the tumor (e.g., CN = 1, CN = 2, ...). It will be apparent to those skilled in bioinformatics that this segmentation method reduces noise and facilitates identification of the copy number profile of a sample, regardless of the availability of an absolute copy number reference. However, in the presence of noise, segmentation can produce erroneous results because it suppresses coverage fluctuations caused by copy number variations, particularly in low-coverage-depth sequencing. Therefore, in a preferred embodiment, the genomic analysis method described herein can be run without a segmentation step. In another embodiment, a segmentation step can be applied to preprocess the data before applying the machine learning methods disclosed herein.

[0103] Following possible normalization and / or segmentation steps, the coverage data typically consists of a set of one-dimensional vectors, where each vector corresponds to a chromosome. In a possible implementation, the coverage data may also consist of a set of one-dimensional vectors, where each vector corresponds to a chromosome arm. In a possible implementation, the size of the vectors can be variable, reflecting the different sizes of different chromosomes. In another possible implementation, the size of the blocks is variable so that each chromosome vector (or each chromosome arm vector) has the same size.

[0104] Figure 4 Exemplary coverage data 1D signal of a whole-genome NGS sequencing experiment was depicted, which was performed with 1X coverage on the first 22 chromosomes of the human genome (excluding sex chromosomes), 100 kbp blocks, normalized by sample and GC content. Figure 5 The same coverage data signal is shown in the arrangement on each chromosome, with the collapsed block being 3 Mbp.

[0105] Spatial arrangement of coverage data

[0106] In a preferred embodiment, the coverage data can therefore be further spatially arranged 210 to form a multidimensional data array, wherein the coverage data is organized to align each chromosome or chromosome arm with blocks of their telomeres or centromeres along any of these blocks in the array structure. This multidimensional array structure may be particularly well-suited to benefiting from efficient genomic analysis machine learning training models recently developed for non-genomic applications such as image classifiers, rather than relying on simple human-engineered features as in existing HRD detection methods.

[0107] In other words, and consistent with the provided definition of telomere or centromere blocks, the coverage data signal vector can be further arranged in space 210 to form a coverage data signal image, wherein the coverage data is organized into blocks aligned with each chromosome or chromosome arm along at least one centromere block adjacent to a centromere region or at least one telomere block adjacent to a telomere region in a two-dimensional array structure.

[0108] In a possible implementation, the coverage data signal vector can be further arranged in space 210 to form a one-dimensional coverage data array, wherein the coverage data is organized as tail-to-head connections of blocks along each chromosome or chromosome arm of the genome.

[0109] In another preferred embodiment, the coverage data signal vector may be further arranged in space 210 to form a multidimensional coverage data signal array, wherein the coverage data is organized into blocks aligned with each chromosome or chromosome arm along at least two blocks adjacent to the centromere region in the array structure, one block being to the left of the centromere and the other block being to the right of the centromere, and the multidimensional data array is a two-dimensional array, specifically a two-dimensional image.

[0110] In another preferred embodiment, the coverage data can be further arranged spatially 210 to form a multidimensional data array, wherein the coverage data is organized to align blocks of each chromosome or chromosome arm along at least two blocks adjacent to telomere regions in the array structure, wherein one block is located at the chromosome origin within the set of coverage blocks associated with the chromosome, and the other block is located at the chromosome end within the set of coverage blocks associated with the chromosome, and the multidimensional data array is a two-dimensional array, particularly a two-dimensional image.

[0111] In another preferred embodiment, the coverage data may be further arranged in space 210 to form a multidimensional data array, wherein the coverage data is organized to compare each chromosome or chromosome arm along at least one or at least two blocks adjacent to the centromere region or at least one or at least two blocks adjacent to the telomere region in the array structure, wherein this multidimensional data array is a two-dimensional array, specifically a two-dimensional image.

[0112] In this case, the block adjacent to the centromere or telomere region means that it is the block closest to these regions.

[0113] In a further preferred embodiment, the blocks of each chromosome arm are calculated independently.

[0114] In one implementation, the method may include, but is not limited to, the step of further normalizing the coverage data for each chromosome arm by taking the median normalized coverage over a window of 2.5 Mbp to 3.5 Mbp. For each chromosome arm, the exact block resolution can be independently chosen between these two values, resulting in the minimum possible incomplete blocks at the ends of the chromosome arms. The final incomplete blocks are discarded as input to a machine learning model, such as a CNN. Therefore, the coverage signal contained in the last block of each p-arm is geographically closest to its corresponding centromere. Similarly, the coverage signal contained in the first block of each q-arm is geographically closest to its corresponding centromere. Simultaneously, the coverage signal contained in the first block of each p-arm is geographically closest to its corresponding p-arm telomere region. Likewise, the coverage signal contained in the last block of each q-arm is geographically closest to its corresponding q-arm telomere region.

[0115] Therefore, in Figure 3 It can be understood and described that the centromere region may not be included in these blocks. Therefore, the last block of the p-arm and the first block of the q-arm actually mark the centromere region of a chromosome.

[0116] Tissue coverage data, when compared to the blocks described in this paper, allows for the detection of spatial relationships between detected events, such as genomic instabilities within chromosome arms. Machine learning models, such as CNNs, can utilize the spatial distribution information of detected events, thus potentially improving their output predictions when accurate state detection benefits from this type of spatial arrangement.

[0117] Figure 6 A graphical representation of chromosomes arranged according to centromere position is shown. (From Table 1 and...) Figure 6It can be seen that the longest chromosome chr1 is 249 Mbp in length, the longest p-arm is 124.7 Mbp in length for chr1, while the longest q-arm is actually even longer in chr2, at 147.9 Mbp. In a possible implementation, within a given constant block size, the coverage data can be adapted to a two-dimensional array of (upper bound (chr_p-arm_max_len) / block_size + upper bound (chr_q-arm_max_len) / block_size)*N, where N represents the number of chromosomes, chr_p-arm_max_len is the longest p-arm length, and chr_q-arm_max_len is the longest q-arm length among the N chromosome sets to be jointly analyzed. As an illustrative example, whole-genome sequencing coverage data for the arms of 22 non-sex chromosomes can be arranged in centromere blocks with a constant block size of 3 Mbp. An array of 92*22 (or 22*92) blocks can be spatially arranged 210 by the genome data analyzer 120 as data input to the CSI analyzer 123. The coverage data for the longest p-arm, 124.7 Mbp (chr1), might indeed fit into 42 3 Mbp blocks, while the longest q-arm, 147.9 Mbp (chr2), might fit into 50 3 Mbp blocks. In a possible implementation, the multidimensional array can be arranged into blocks representing the longest chromosome arm among the 2*N chromosome arms in the N chromosome sets to be analyzed. For example, the coverage data can be spatially arranged into a 44*50 array for analyzing 22 chromosomes with a constant block size of 3 Mbp. The coverage data for each chromosome arm can be arranged in rows (or columns) in the array so that the centromere blocks of the chromosome arms remain consistent throughout the array. In a possible implementation, the following can be applied: Figure 3 The above schematically illustrates adaptive blocks to remove irrelevant high-resolution blocks and adjust the size and position of variable blocks along the chromosome genome. For example, an array of size 84*22 (or 22*84) can be spatially arranged 210 by the genome data analyzer 120 as data input to the CSI analyzer 123. Other implementations are also possible.

[0118] Because the number of blocks varies with the length of the chromosome and / or chromosome arms, some elements of the array need to be filled with dummy data that is not present in the coverage data. For those skilled in data science, there are different options for filling empty entries (or elements) in a multidimensional array. In one possible implementation, empty elements can be filled with constant values. In another possible implementation, empty elements can be filled from a mask array of predefined constant values. In yet another possible implementation, empty elements can be filled based on the coverage data of the portions of rows or columns that are filled with real data. In a possible implementation, empty elements can be filled by repeating coverage values, moving left or right (above or below, respectively) from the last available block in the column (or row) of a block along a chromosome or chromosome arm row (or column).

[0119] In another possible implementation, the individual size of the blocks can be selected to facilitate dense array filling. Smaller block sizes (corresponding to better resolution) can be selected in regions with specific genomic significance to increase the total number of blocks for certain chromosome arms.

[0120] Figure 7 Examples of normalized coverage data spatially rearranged into two-dimensional array images suitable for CSI analysis using machine learning models are shown, from a) HRD-negative tumor DNA samples and b) HRD-positive tumor DNA samples, respectively. Chromosomes are arranged in rows, and coverage blocks are arranged along columns, such that, according to some embodiments of the proposed coverage data preparation method, centromere blocks are vertically arranged in a single column prior to chromosome spatial instability analysis. The brighter the pixel, the higher the normalized coverage count value C of the plotted block. 区块 The higher the value, the darker the pixel; conversely, the lower the value, the more likely it is to be.

[0121] Chromosome Spatial Instability (CSI) Analyzer

[0122] The genomic data analysis system 250 can further apply the chromosome spatial instability analyzer 123 to spatially arranged coverage data multidimensional arrays to automatically analyze and report one or more genomic characteristics representing certain spatial instability patterns of two or more chromosomes or chromosome arms in each DNA sample. In the case of tumor sample analysis, the CSI analyzer module 123 can report the burden of chromosomal aberrations, such as large deletions or duplications of alleles in specific genomic regions along each chromosome arm.

[0123] In a possible implementation, the CSI analyzer module 123 is a computer-implemented algorithm trained to distinguish between DNA samples with genomic instability (marked as HRD positive) and another DNA sample without genomic instability (marked as HRD negative).

[0124] In a possible implementation, in the case of a tumor sample, the CSI analyzer module 123 can identify a signature that can identify samples with similar characteristics, which may lead to a good response to a given treatment.

[0125] CSI analyzer module 123 can accordingly identify, classify, and report combinations of indicators, scalar scores, or characteristic indicators to describe one or more genomic characteristics of a tumor sample. In a preferred embodiment, CSI analyzer module 123 can accordingly identify, classify, and report the status of the tumor sample, such as HRD-negative or HRD-positive, or possibly undetermined (uncertain). In a possible embodiment, CSI analyzer module 123 can report scalar scores as an indicator of the likelihood of HRD in the tumor sample. It will be apparent to those skilled in the art that this bioinformatics approach subsequently greatly facilitates a detailed understanding of genomic alterations in cancer cells within tumor samples and enables personalized medicine to be tailored to the specific circumstances of a patient's inferred cancer cell biology.

[0126] In one implementation, the genome characterization report can be displayed to the end user on a graphical user interface. In another possible implementation, the genome characterization report can be prepared as a text file for further automated processing. Other implementations are also possible.

[0127] In a preferred embodiment, the CSI analyzer module 123 can be configured to apply a machine learning model to a spatially arranged coverage data multidimensional array. In a possible embodiment, the CSI analyzer module 123 may include a trained neural network or a combination of trained neural networks as a data processing module to analyze the spatially arranged coverage data multidimensional array 220 as input to the trained neural network and produce CSI results as outputs of the trained neural network.

[0128] Now refer to Figure 8 The exemplary CSI analyzer 123 and its data processing workflow are described in further detail. It will be apparent to those skilled in the art that various methods and architectures known in the field can be used to process the prepared coverage data multidimensional array input and derive the CSI status or CSI score of the tumor sample.

[0129] Exemplary CNN machine learning model

[0130] To classify multidimensional arrays (such as two-dimensional images) where chromosome arms are aligned along their respective centromere positions to highlight visual patterns characterizing the CSI state of the chromosome set, convolutional neural networks (CNNs) can be built and trained using methods and techniques known in the field of image pattern recognition, such as TensorFlow, Caffe, Caffe2, PyTorch, Theano, or others.

[0131] Figure 8 An example of a convolutional neural network (CNN) is illustrated, which can process a matrix image of the spatial arrangement coverage data of a patient's tumor DNA sample to classify it as HRD positive (HRD+) or HRD negative (HRD-). Figure 8 In one example, the trained CNN takes as input a spatially arranged coverage data matrix image 800, which is arranged in rows of 22 non-sex chromosomes, from chr1 at the top to chr22 at the bottom, and in columns of blocks (in... Figure 8 In the exemplary illustration, the image is approximately 3 Mbp in size, arranged along the centromere region of the chromosome, with essentially zero coverage (although other implementations are possible). The input image 800 is fed into a first convolutional neural network 810 with several processing layers to quantify a set of intermediate features, which are themselves fed into a second neural network classifier 820, responsible for extracting two scores (one for HRD- and one for HRD+) as its two outputs. The CSI analyzer can then determine the HRD state (HRD-, HRD+, or indeterminate) by comparing and thresholding the output values ​​of HRD- and HRD+.

[0132] More generally, it will be apparent to those skilled in the art that the CSI analyzer 123 can be adapted to employ various CNN architectures. The convolutional network may include a first series of sequentially connected convolutional data processing layers 810, which may be further fed into a second series of fully connected layers 720 to output various predictions. For example, the CNN architecture may include one or more convolutional layers arranged to process one-dimensional, two-dimensional, or multi-dimensional data; optionally, one or more ensemble layers of maximum, median, or average values ​​arranged to process one-dimensional, two-dimensional, or multi-dimensional data; one or more multiple filtering layers for regularization during training; optionally, one or more batch regularization layers; flattening layers followed by one or more fully connected layers; or combinations thereof. Various activation methods known in the art, such as ReLU, SELU, or more generally, any activation function that prevents gradient vanishing, may be used. Different output activations can also be generated using methods known in the art, such as the sigmoid function method to produce probabilities associated with single-label binary classification (e.g., the probability of a positive state, the probability of a negative state = 1 - the probability of a positive state); or using the softmax method to generate probabilities associated with single-label multi-class classification outputs (e.g., mutually exclusive probabilities of positive, negative, or uncertain states); or using regression methods, such as linear methods or ReLU classification methods, to generate scalar values ​​that are not intended to be interpreted as probabilities, such as HRD scores.

[0133] In a preferred embodiment, the CSI analyzer 123 generates a final output decision, such as HRD+ or HRD- status analysis for each sample, by thresholding the output of the CNN. The threshold can be adjusted according to the application to optimize the sensitivity and / or specificity of the CSI analysis based on the end-user's needs, such as for diagnosis, prognosis, or guiding the selection of the most effective cancer treatment. In a possible embodiment, the CSI analyzer 123 may also further report the output values ​​of the 230CNN as needed by the end-user.

[0134] Figure 9An example of a possible detailed CNN architecture is shown, which can be used by a CSI analyzer 123 to analyze coverage data with a 220 spatial arrangement. For each convolutional layer within the network, the number of filters and the filter size are represented as [N x H x W], where N is the number of filters, H is the height of the filters, and W is the width of the filters. Each convolutional layer is activated by an activation function (Regularized Linear Unit, ReLU) and then batch normalized (BN). Intermediate convolutional layers can be interleaved with max-pooling layers, and then, following the typical CNN structure, are interleaved with average pooling layers and flattening layers. In this example, the CNN's output layer is flattened, passed through batch normalization layers, dropout layers, and fed into a single-node dense layer (i.e., a fully connected layer) using a sigmoid activation function.

[0135] Chromosome Spatial Instability (CSI) Analyzer Training

[0136] For supervised and semi-supervised CSI analyzer training, a set of labeled input-output pairs can be used as the training set. For example, to train a CSI analyzer to classify HRD states, a subset of data samples labeled with public-domain HRDs can be used, allowing the network to adjust its parameters (such as weights) along its different layers through an optimization process (such as backpropagation) to minimize its output error; this optimization process is specifically applied during the training phase. For those skilled in machine learning, it is clear that different loss functions can be used to measure the error; for example, binary cross-entropy can be used for models with single-label binary classification outputs, classification cross-entropy can be used for models with single-label multi-class classification outputs, or regression methods based on root mean square error, absolute error, or other error measures. Optimization can employ existing stochastic gradient descent-based optimizers such as ADAM or RMPSprop. Regularization can employ dropout methods, early stopping methods, and / or L-1 or L-2 parameter regularization. Other implementations are also possible.

[0137] As more clinical conditions related to chromosomal instability are categorized in current and future oncology research, additional datasets can be used to train machine learning models designed to detect genomic scarring caused by HRD or other phenomena, such as response to specific treatments. In a possible implementation, binary states of HRD+ and HRD-, as well as an optional “uncertain” state, can be used as multi-class underlying ground truth labels. In another possible implementation, HRD scores from existing methods (such as HRDetect) can be used as scalar output labels. According to any possible implementation of the methods disclosed herein, the input coverage signal data array can also be spatially arranged as multidimensional data, provided that the choice of the multidimensional array format implementation remains unchanged during training and runtime.

[0138] In a possible implementation, the additional dataset used to train the machine learning model may include an array of coverage data signals from patient DNA samples, wherein these patients are known to have a homologous recombination defect, and / or wherein these patients have received a cancer treatment regimen and the treatment outcome is known (response vs. non-response).

[0139] In a possible implementation, additional datasets for training machine learning models may include an array of coverage data signals from DNA samples generated from primary or immortalized cell line experiments or tumor-derived organoid experiments, particularly experiments that determine whether organoids predicted to have HRD are eliminated or have reduced growth by cancer treatment regimens.

[0140] In a possible implementation, a semi-supervised learning framework can be employed to train the machine learning model. Indeed, acquiring labeled data can be challenging and costly, especially in the context of HRD state prediction, making semi-supervised learning a promising framework for this application. Since the degree of overfitting of a model is determined by its complexity and the amount of training it receives, providing the model with more training instances can help reduce overfitting. Data augmentation involves creating new artificial training samples from existing data. In a possible implementation, artificial samples can be generated by extracting chromosomes from existing real samples that share the same HRD state and have similar tumor content (or biomimetic preprocessing to mimic similar tumor content) and ploidy.

[0141] In a possible implementation, artificial samples can be generated by sampling chromosomes from available real samples that share the same HRD state and the same normalized purity and ploidy ratio, where purity is the percentage of tumor cells in the sample and ploidy is the average number of chromosome copies in the entire genome, thus characterizing how many complete sets of chromosomes are present on average in tumor cells.

[0142] In one implementation, the machine learning model is trained on a coverage data signal array data from labeled real samples.

[0143] In another implementation, the machine learning model is trained on a coverage data signal array, which consists of labeled real samples and labeled data-enhanced samples generated by sampling chromosomes from available real samples with the same HRD state and ensuring that the sampled chromosomes have the same purity / ploidy ratio.

[0144] Therefore, any machine learning model can learn to quantify the difference between a DNA sample with an HRD positive status and another DNA sample with an HRD negative status as a CSI score, and then classify the tested DNA samples based on the CSI score results.

[0145] In possible implementations, the machine learning model is trained in a supervised or semi-supervised mode, using at least one set of real samples labeled with HRD states according to the HRDetect method. Preferably, the learning model is trained in a supervised mode using at least one set of real samples labeled with HRD states according to the HRDetect method.

[0146] In one possible implementation, the machine learning model is trained in a supervised or semi-supervised mode using artificial samples generated by sampling chromosomes from a set of real samples sharing the same HRD state, taking into account differences in purity ploidy ratio. In another possible implementation, the machine learning model is trained in a supervised mode using artificial samples generated by sampling chromosomes from a set of real samples sharing the same HRD state, taking into account differences in purity ploidy ratio.

[0147] In one implementation, a trained machine learning model calculates the CSI score of a subject's DNA sample and uses the CSI score to classify the probability that the subject's DNA sample belongs to a sample group with a certain CSI state.

[0148] In one implementation, a trained machine learning model calculates the CSI score of a subject's DNA sample and uses the CSI score to determine the HRD state of the subject's DNA sample belonging to a set of samples.

[0149] In one implementation, a computer-based method is provided for determining the CSI status of a DNA sample, wherein the CSI status is the HRD status.

[0150] In one implementation, a computer-based method is provided for determining the CSI status of a DNA sample, wherein the DNA sample is DNA from a cancer sample.

[0151] In one embodiment, a computer-based method is provided for determining the CSI status of a DNA sample, wherein the low-pass whole-genome sequencing coverage is at least 0.1X and at most 30X, such as 1X to 10X, such as 0.1X to 5X, or such as 0.1X to 1X.

[0152] In one embodiment, a computer-based method is provided for determining the CSI status of a DNA sample, wherein low-throughput genome sequencing is performed at a rate of at least 0.1X and at most 30X, such as 1X to 10X, such as 0.1X to 5X, or such as 0.1X to 1X, to obtain a set of sequencing reads of a subject's DNA sample containing reads from at least 2 and at most 22 chromosomes.

[0153] In one implementation, a computer-based method is provided for determining the CSI status of a DNA sample, wherein the process involves comparing sequencing reads of the DNA sample with a human reference genome.

[0154] In one implementation, a computer-based method is provided for determining the CSI status of a DNA sample, wherein normalization by GC content is performed to apply GC-offset correction.

[0155] In one implementation, a computer-based method is provided for determining the CSI status of a DNA sample, wherein the CSI status is an HRD status.

[0156] In one implementation, a computer-based method is provided for determining the CSI state of a DNA sample, wherein a machine learning model is pre-trained using a training dataset of samples with known genomic instability states such as HRD+ / HRD- states, and thus the machine learning model is trained to distinguish samples with a certain CSI state such as HRD state based on coverage spectral characteristics.

[0157] Diagnostic / Prognostic Methods

[0158] In one embodiment, an in vitro method is provided for characterizing a DNA sample obtained from a patient, the method comprising:

[0159] - DNA fragments were isolated from patient samples;

[0160] - Construct a sequencing library containing the fragments that overlap with a set of chromosomes;

[0161] - Sequencing the library to achieve up to 30X sequencing coverage;

[0162] - Compare the obtained sequencing reads with the human reference genome;

[0163] - and generate the CSI score of the subject sample using the method of the present invention.

[0164] In one embodiment, an in vitro method is provided for characterizing a DNA sample obtained from a patient, the method comprising:

[0165] - DNA fragments were isolated from patient samples;

[0166] - Construct a sequencing library containing the fragments that overlap with a set of chromosomes;

[0167] - Sequencing the library to achieve up to 30X sequencing coverage;

[0168] - Compare the obtained sequencing reads with the human reference genome;

[0169] - And based on the CSI score obtained using the method of the present invention, the HRD status of the patient sample is determined.

[0170] A method for selecting cancer patients to be treated with PARP inhibitors includes the step of detecting HRD positivity in tumor patient samples using any of the methods described in this invention.

[0171] In one implementation, a patient with cancer can be classified as HRD positive (HRD+) or HRD negative (HRD-) using any of the methods described herein.

[0172] In one embodiment, an in vitro method is provided for characterizing a DNA sample obtained from a patient, wherein the patient sample is a tumor sample, and an indicator of chromosomal spatial instability of the patient sample is a predictor of the tumor’s response to a cancer treatment regimen, including platinum-based chemotherapy agents, DNA damaging agents, anthracyclines, topoisomerase I inhibitors, or PARP inhibitors.

[0173] In one embodiment, an in vitro method for characterizing DNA samples obtained from a patient is provided, including a step of identifying whether a tumor patient sample has a homologous recombination (HR) defect based on a chromosomal spatial instability index, wherein a high probability of a homologous recombination (HR) defect indicates that a PARP inhibitor may be used to treat cancer. PARP inhibitors may be used alone or in combination with other treatment methods.

[0174] In one embodiment, a method is provided for selecting cancer patients to receive treatment with platinum-based chemotherapy agents, DNA damaging agents, anthracyclines, topoisomerase I inhibitors, or PARP inhibitors, including the step of detecting that the patient's DNA cancer sample is HRD-positive according to a computer-based method for determining the CSI status of a DNA cancer sample as described herein. In one embodiment, the method for selecting cancer subjects for treatment is an in vitro method.

[0175] In one embodiment, a computer-based method is provided for determining the CSI status of a DNA cancer sample described herein, wherein the method further includes:

[0176] When the DNA cancer sample is determined to be defective due to homologous recombination, the cancer is treated by administering a poly-ADP-ribose polymerase (PARP) inhibitor to the test subject.

[0177] According to one aspect, the sample is a patient's sample, which may be in the form of tissue, fresh frozen tissue (FFT), blood or any bodily fluid, or cytological specimens / preparations (FFPE, smears), etc. According to another aspect, the sample is a patient's tumor sample, including FFPE samples.

[0178] In another specific implementation, the patient has cancer, particularly high-grade serous ovarian cancer, prostate cancer, breast cancer, pancreatic cancer, etc.

[0179] According to one embodiment, a method is provided for diagnosing cancer type HDR+ or HDR- according to the method described herein, wherein a sample is obtained from a subject who has or is suspected of having cancer.

[0180] Wet laboratory workflow Procedure Integration

[0181] - CSI Analyzer Data Generation - WGS Workflow

[0182] In a possible implementation, a method for determining the CSI status of a subject's DNA sample is provided, the method comprising:

[0183] a) Provide sample materials including nucleic acid;

[0184] b) Prepare the first nucleic acid sequencing library for whole genome sequencing;

[0185] c) Sequencing the first nucleic acid sequencing library obtained in step b);

[0186] d) Analyze the nucleic acid sequences obtained in step c), where,

[0187] The sequence of nucleic acid (obtained in step c) is analyzed according to the computer-based method for determining the CSI status of a subject's DNA sample as described herein.

[0188] - CSI Analyzer Data Generation - WGS Workflow and Variant Call Program Target Enrichment Workflow

[0189] In a possible implementation, a method for determining the CSI status of a subject's DNA sample is provided, the method comprising:

[0190] a) Provide sample materials including nucleic acid;

[0191] b) Prepare the first nucleic acid sequencing library for whole genome sequencing;

[0192] c) Sequencing the nucleic acid sequencing library obtained in step b);

[0193] d) Analyze the nucleic acid sequences obtained in step c), where,

[0194] The sequence of nucleic acid (obtained in step c) is analyzed according to the computer-based method for determining the CSI status of a subject's DNA sample as described herein.

[0195] The method also includes:

[0196] e) Prepare a second nucleic acid sequencing library from the same sample material;

[0197] f) Targeted enrichment of the second nucleic acid sequencing library;

[0198] g) Sequencing the target-enriched library obtained in step f);

[0199] h) Analyze the nucleic acid sequences obtained in step g), where,

[0200] Analyze the sequences from the enriched nucleic acid library (obtained in step f) according to any known variant invocation method.

[0201] It is obvious to those skilled in the field of genomic analysis that the first and second libraries can be prepared independently from the same sample at the same time or at different times, in any order. In a possible implementation, the first and second libraries can be obtained by splitting a library into two subsets, one subset (the first library - the WGS library for CSI analysis) proceeding through steps c) and d), while the other subset (the second library - the targeted enrichment library for variant call analysis) proceeds through subsets e) through h).

[0202] In another embodiment, a method for determining the CSI status of a subject's DNA sample is provided, the method comprising:

[0203] a) Provide sample materials including nucleic acid;

[0204] b) Prepare the first nucleic acid sequencing library;

[0205] c) Targeted enrichment of the nucleic acid sequencing library to obtain a second enriched nucleic acid sequencing library;

[0206] d) Sequencing the nucleic acid sequencing libraries obtained in steps b) and c);

[0207] e) Analyze the nucleic acid sequences obtained in step d), where,

[0208] After excluding sequences aligned with the enriched regions in c), sequences from non-targeted genomic regions of the first non-enriched nucleic acid sequence library (obtained in step b) were analyzed according to the computer-based method for determining the CSI status of subject DNA samples described herein.

[0209] and

[0210] Analyze the sequences from the second enriched nucleic acid library (obtained in step c) according to any known variant invocation method.

[0211] It is obvious to genomics professionals that the above workflow can be achieved using different genomic analysis steps.

[0212] The steps for providing sample materials including nucleic acids can be performed using any known DNA extraction method. If multiple samples are sequenced together in one experiment, genome barcoding methods can be used to identify the sample materials for downstream analysis.

[0213] The steps for preparing a nucleic acid sequencing library can be performed according to any known method for preparing a sequencing library. Optionally, step b) of preparing a nucleic acid sequencing library may include a further nucleic acid sequence amplification step b.0. Nucleic acid amplification can be performed according to any known method, such as polymerase chain reaction (PCR).

[0214] The targeted enrichment steps can be performed using any known targeted enrichment DNA-sequencing or RNA-sequencing method, such as targeted hybridization capture (i.e., targeted sequencing based on hybridization capture) or amplicon-based methods (amplicon sequencing).

[0215] According to one embodiment, targeted enrichment is capture-based targeted enrichment. According to this embodiment, the steps of performing targeted enrichment include hybridizing at least one probe or probe set with a target nucleic acid (“target region” or “enriched region”) from a genomic region known to potentially carry variants of interest (e.g., but not limited to BRCA1 and BRCA2 regions), washing away non-target nucleic acids, and enriching the nucleic acid. It is understood that the term “probe” or the interchangeable terms “decoy” or “(probe) nucleic acid molecule” or “capture probe” or “(DNA / RNA) oligonucleotide (capture) probe” refer to a nucleic acid molecule that can hybridize with a target nucleic acid molecule. Any known probe design can be used. The term “target nucleic acid” refers to a nucleic acid region in a gene or transcript that can be captured by the probe used. Preferably, the target nucleic acid is a gene associated with the HRR pathway, such as selected from, but not limited to, BRCA1 and BRCA2.

[0216] It is understandable that the nucleic acid sequencing library obtained in step b) can be divided into at least two parts. One part will not be targeted enriched and will therefore remain a non-enriched nucleic acid library; while the other part will be targeted and form an enriched nucleic acid library.

[0217] Understandably, non-enriched and enriched nucleic acid libraries can be loaded onto the sequencer together or separately. Therefore, depending on the chosen workflow, nucleic acid sequences can be obtained individually or together.

[0218] When non-enriched and enriched nucleic acid libraries are loaded together into a sequencer, these libraries are loaded at predetermined concentrations. These relative concentrations are a function of the desired coverage of the targeted (enriched) and non-targeted regions, such that the nucleic acid coverage of the non-enriched library is 0.1X to 10X, 0.5X to 10X, at least 0.2X, or preferably at least 0.1X and at most 30X; while the nucleic acid coverage of the enriched library is at least 30X, or 30X to 100X, or 100X to 500X, or 500X to 1000X, or 1000X to 5000X, preferably at least 4000X.

[0219] The sequencing steps for nucleic acid libraries can be performed using any known method and any known sequencer.

[0220] In a possible implementation, sequence analysis from non-enriched nucleic acid libraries and sequence analysis from enriched nucleic acid samples are performed separately.

[0221] When two libraries are sequenced together, analysis of sequences from the non-enriched nucleic acid library to obtain CSI status using the computer-based method described herein may include further steps to filter out regions covered in the targeted enriched library. These regions can be masked or filtered by a CSI analyzer according to any known masking strategy, such as:

[0222] m1) This involves simply aligning sequencing reads to a reference genome at genomic regions not covered by the targeted enrichment library; or

[0223] m2) By selecting only the aligned genomic regions not covered by the target-enriched library in the alignment file (e.g., BAM or SAM file); or

[0224] m3) This is achieved by excluding high-resolution blocks from the CSI score determination, i.e., blocks that overlap with the target genomic region (e.g., by masking blocks that overlap with the genomic region covered by the target enriched library).

[0225] In the joint sequencing of non-enriched and enriched libraries, the step of masking the targeted enriched regions of the genome obtained from the enriched library ensures that the coverage differences introduced by targeted enrichment do not affect CSI analysis.

[0226] In a possible implementation, analysis of sequences from an enriched nucleic acid library includes variant calling according to known methods. Data on variant allele fractions (VAFs) obtained from this analysis are additional information that can be used to further characterize the patient DNA sample. In a possible implementation, the HRD status of a sample can be derived from the CSI status, variant calling results (e.g., if one or more variants in genes associated with homology repair deficiency are identified, such as, but not limited to, certain variants in the BRCA1 or BRCA2 genes), or any combination thereof (e.g., if the CSI status is identified as positive or if one or more variants in genes associated with homology repair deficiency are identified). In a possible implementation, variant calling analysis may be performed first, and CSI analysis of the sample may be performed according to the computer-based methods described herein only if the variant calling results do not identify the HRD status of the sample, in order to identify possible HRD statuses from genomic instabilities across the sample's genome.

[0227] Experiment - Example 1

[0228] In the first experiment, Figure 9The CNN was trained using low-pass WGS coverage data measured in a subset of 133 fresh frozen tumor samples from 560 breast cancer patients, where BRCA deficiency state and HRD state (defined by HRDetect score) were publicly available. See Nik-Zainal et al., “Landscape of somatic mutations in 560 breast cancer whole-genome sequences”, Nature 534, 47–54 (2016) – from datasets from the Wellcome Trust Sanger Institute and the International Cancer Genome Consortium ICGC, available at the European Genome-phenome Archive EGA (https: / / www.ebi.ac.uk / ega / studies / EGAS00001001178), and see “HRDetect is a predictor of BRCA1 and BRCA2 deficiency based on mutational signatures”, H. Davies et al., Nature. Published online in Medicine on March 13, 2017. The original training dataset, consisting of 133 publicly available real samples, was augmented to 3083 artificial training samples. Data augmentation was performed to maintain the same purity-ploidy distribution observed in the original dataset. The artificial samples obtained through data augmentation were used only for training and validation. A large portion of the original dataset (69 out of 202 samples from 560 breast cancer studies) was reserved for testing and did not participate in the data augmentation process.

[0229] The original BAM file for each tumor sample was first downsampled to 10 million paired-end reads to simulate a low-pass WGS. The coverage signal was first pre-computed using a first set of 100kbp high-resolution coverage blocks and then preprocessed by performing sample and GC normalization. The resulting normalized coverage signal (y-axis) is shown in the figure. Figure 10a HRD-samples in ) and Figure 10b The whole genome (excluding sex chromosomes) of the HRD+ sample in the image is shown on the x-axis. The x-axis represents the genomic coordinates of the vertical dashed lines at the boundaries of each chromosome.

[0230] use Figure 3 The adaptive blockization strategy further collapses the normalized coverage data into an average target block size of 3Mbp, ranging from 2.5Mpb to 3.5Mbp. Figure 11The diagram shows the coverage data obtained by sorting 133 original training samples by their CSI status as a one-dimensional signal before applying the spatial arrangement method according to some embodiments of this disclosure, with each sample on a separate line and HR-deficient samples at the top of the Y-axis. The X-axis shows the genomic coordinates of the vertical dashed lines at each chromosome boundary.

[0231] The normalized coverage data of the samples in the training dataset were further spatially arranged into a two-dimensional array of 84 blocks * 22 chromosomes, used to train the CNN model to predict HRD states. The coverage blocks of the 22 chromosomes were plotted in 22 rows, from chr1 to chr22, arranged relative to their centromere blocks. For shorter chromosome arms, empty blocks were filled by copying the value of the nearest telomere block in the same row. Figure 7 a) shows an array of the resulting HRD samples. Figure 7 b) shows the array of the resulting HRD+ samples.

[0232] Figure 12 The left side shows the results of applying the trained CNN to the full dataset of 202 samples, including 69 test samples not used for training (represented as “BRCAness” state indicators measured at the CNN output and then converted to probabilities via softmax), and the right side shows the results of the state-of-the-art HRD scores (obtained by averaging LOH, LST, and TAI scores) against the HRDetect score of the classifier. Figure 12 As shown in the histogram, the proposed method outperforms existing HRD score methods in predicting HRDetect positive samples, whether tested on the entire set of 202 real samples or only on the 69 samples not used to train the machine learning model.

[0233] Figure 13 The image shows the three HRDetect scores (top), HRD scores (bottom), and SOPHiA CSI (represented here as BRCAness) scores obtained at the output of the trained CNN (middle) for each of the 69 samples in the test set. In each individual image, samples are sorted from left to right in ascending order of score. Samples without BRCA defects are shown in light gray. Samples with BRCA defects are shown in dark gray (different gray levels are used to represent the mutation types that cause BRCA defects). Figure 13 The results confirm that all BRCA defect samples had very high HRDetect scores and SOPHiA BRCAness CSI scores. This result indicates that the CSI index is superior to the HRD score in identifying BRCA defect samples.

[0234] Experimental Example 2 - High-grade ovarian serous adenocarcinoma Ovkate in a single wet laboratory workflow Cells undergo lpWGS combined with capture-based target plates

[0235] Low-throughput whole-genome sequencing and capture-based target plates were performed on OVKATE (RRID: CVCL_3110) cells (https: / / web.expasy.org / cellosaurus / CVCL_3110). DNA was extracted using the DNeasy Blood and Tissue Kit (Qiagen) according to the manufacturer's instructions. Whole-genome sequencing libraries were then prepared using the SOPHiA Gene Library Preparation Kit. A portion of the whole-genome library was used for target enrichment using a probe group covering genes related to the HRR pathway and the SOPHiA gene capture protocol. The enriched library and the whole-genome library were then sequenced together in a Nextseq Mid flow cell (Illumina). Both types of libraries were loaded and balanced on flow cytometry to achieve approximately 1-2X coverage on the genome and over 1000X coverage in the target regions of the capture plate.

[0236] Figure 15 illustrates an exemplary workflow for analyzing data obtained after lpWGS and directional capture.

[0237] DNA extracted from cultured Ovkate cells was processed using a combined low-pass WGS and targeted enrichment NGS workflow. The NGS data underwent bioinformatics processing to generate BAM files containing reads aligned to a reference genome. Reads mapped to targeted regions were used for variant recall and enrichment of regions. Figure 15C (right, Figure 15C The proportion of alleles in C2). After masking regions targeted by the capture enrichment method, reads from the lp-WGS region can be used to calculate a genome-wide coverage profile. The latter is normalized to eliminate bias. Figure 15C (Left, Figure 15C The C1 group was processed using a CNN to calculate the HRD score. Ovkate cells were classified as HRD positive (score 15.8).

[0238] Therefore, combining lpWGS and capture-based target sequencing allows for the parallel detection of sample CSI status and mutation HRD status of genes of interest.

[0239] Experiment - Example 3 - Obtaining Enhanced Training Data

[0240] In one experiment, 169 real samples of fresh frozen tissue with known HRD states were first selected from publicly available data from Nik-Zainal et al., “Landscape of somatic mutations in 560 breast cancer whole-genome sequences”, Nature 534, 47–54 (2016). For each artificial sample, chromosomes from a random number of real samples with the same HRD state were randomly combined to create data augmentation (DA) training samples. The number of samples used to create the artificial samples can be a random number drawn from an exponential distribution, using: N = K * exp(-K * x), where x is a random number and K = 1 / 3. Using this method ensures that only a finite number of sample data are combined into artificial samples.

[0241] Next, the sample with the lowest purity / ploidy ratio in the pool is identified, and the purity of all other samples is reduced until the purity / ploidy ratio of all samples in the pool equals this ratio. This prevents the introduction of differences in coverage amplitude between chromosomes randomly selected from samples with different tumor purities and ploidies. Figure 16 A). To achieve this, purity is reduced in a computer (in-silico) by adding sequencing reads from normal samples. Figure 16 B). This strategy ensures that, for data-enhanced samples, the magnitude of the difference in coverage for a given ploidy is constant across all chromosomes in all DA samples. Figure 16 C), but this introduces a bias in the purity / ploidy ratio of the DA training data. In particular, the purity and ploidy distributions of these samples tend to be lower than observed in the original sample set. This potential confounding factor was accounted for by applying the Metropolis-Hastings and Gibbs sampling methods, and it was ensured that the purity / ploidy distributions of the resulting 4403 retained DA samples matched those of the original samples.

[0242] Such a set of data augmentation (DA) training data can be used to train machine learning models. It is evident to those skilled in machine learning that such data augmentation strategies can facilitate the training of CSI classifier machine learning models when a large amount of data that might be needed to support the performance of certain classes of methods, such as HRD, is unavailable. Data augmentation can be prepared on a computer to increase the size and diversity of the training dataset. Successful data augmentation strategies can preserve the key attributes of the augmented dataset, making it indistinguishable from real data.

[0243] Further advantages and benefits of the proposed method

[0244] The proposed machine learning methods offer further advantages because they are inherently adaptable to the increasing availability of personalized medical data, particularly in oncology practice. They do not require explicit biological models (e.g., preferably suspected events occurring in certain chromosomal arm regions, such as telomeres or around centromeres) to establish predictive features suitable for guiding diagnosis, prognosis, and / or treatment. When employing a semi-supervised training framework, data augmentation facilitates the development of more robust data models that are less sensitive to noisy data or low tumor content in DNA samples. As more data becomes available, models can be retrained without requiring a complete overhaul of the existing genomic analysis system and workflow, as only the model parameters change. Training models initially developed for a specific application (e.g., breast cancer) can be transferred to other applications (e.g., ovarian cancer, prostate cancer, or pancreatic cancer) and / or other sample types (FFPE, FFT, cfDNA, ctDNA). They can also be used to predict different states, including responses to different therapies and treatments.

[0245] The method presented in this paper is applicable to lpWGS sequencing data, and is therefore cheaper and easier to implement and deploy in standard clinical practice compared to existing techniques using high-coverage (>30X) WGS or SNP arrays.

[0246] Other implementation methods and applications

[0247] While various implementation methods have been described above, it should be understood that these methods are presented by way of example and not as limiting. It will be apparent to those skilled in the art that various changes in form and detail can be made without departing from the spirit and scope. In fact, after reading the above description, it will be obvious to those skilled in the art how to implement alternative implementation methods.

[0248] While exemplary implementations and applications of the proposed methods have been associated with targeted next-generation sequencing (NGS) genomic analysis, it will be apparent to those skilled in the art of bioinformatics that they can also be applied to detect and classify HR defects from other tumor sample genomic analysis workflows, such as those using SNP arrays and array CGH wet lab workflows. SNP arrays, array CGH, or next-generation sequencing (NGS) technologies can be used to generate copy number-dependent counts. Furthermore, low-throughput WGS (0.1X to 5X) or large targeted enrichment sequencing libraries, such as those obtained from WES (whole exome sequencing) or CES (capillary electrophoresis) plates, can be used alone or in combination to generate input data for the proposed machine learning methods. Additionally, low-throughput WGS can be combined with small targeted plates (amplifier-based or capture-based) in single or separate wet lab workflows. Off-target reads from targeted sequencing methods can also be fed as primary or supplementary input data into multidimensional arrays, while still remaining well-suited for trained neural networks as input data processing components in the CSI Analyzer 123 architecture.

[0249] While exemplary instances and applications of the proposed method have been described in relation to coverage data signal arrays arranged as images of chromosome coverage blocks in rows and vertically arranged according to their centromere blocks, it will be apparent to those skilled in the art of bioinformatics that various other embodiments are also possible. Figure 14 An example alignment of coverage signal images of each chromosome arm (p arm and q arm) from at least two analyses (chromosome A and chromosome B) relative to at least one block adjacent to a centromere region or telomere block is shown.

[0250] In one possible implementation, chromosome arms can be represented as rows in an array, with the centromere regions of all arms aligned along the first column of the array. In another possible implementation, chromosome arms can be represented as rows in an array, with the telomere regions of all arms aligned along the first column of the array. In yet another possible implementation, chromosome arms can be represented as rows in an array, with the centromere regions of all arms aligned along the last column of the array. In yet another possible implementation, chromosome arms can be represented as rows in an array, with the telomere regions of all arms aligned along the last column of the array.

[0251] In one possible implementation, chromosome arms can be represented as columns in an array, with the centromere regions of all chromosome arms aligned along the first row of the array. In another possible implementation, chromosome arms can be represented as columns in an array, with the telomere regions of all arms aligned along the first row of the array. In yet another possible implementation, chromosome arms can be represented as columns in an array, with the centromere regions of all arms aligned along the last row of the array. In yet another possible implementation, chromosome arms can be represented as columns in an array, with the telomere regions of all arms aligned along the last row of the array.

[0252] In a possible implementation, the entire chromosome can be represented as columns in an array, with centromere blocks of all chromosomes aligned along rows of the array. The alignment rows of the centromere blocks can be at the center of the array, or they can be moved upwards or downwards from the center depending on the block lengths of the p-arms and q-arms of the chromosome set. Blocks associated with the p-arms and blocks associated with the q-arms of the chromosomes can be above and below the centromere block rows, respectively, or conversely, below and above the centromere block rows.

[0253] In a possible implementation, the entire chromosome can be represented as rows in an array, with centromere blocks of all chromosomes aligned along columns of the array. The alignment columns of the centromere blocks can be in the middle of the array, or they can be offset to the right or left from the middle according to the block lengths of the p-arms and q-arms of the chromosome set. Blocks associated with the p-arms and blocks associated with the q-arms of the chromosomes can be moved to the right and left from the centromere block columns, respectively, or vice versa.

[0254] In one possible implementation, the entire chromosome can be represented as rows in an array, and one of the two telomere blocks of all chromosomes can be aligned along a column of the array. In another possible implementation, the alignment column for the telomere blocks can be the first column of the array. In yet another possible implementation, the alignment column for the telomere blocks can be the last column of the array. In one possible implementation, p-arm telomere blocks can be used for alignment. In yet another possible implementation, q-arm telomere blocks can be used for alignment. In yet another possible implementation, the alignment of either the p-arm or q-arm telomere blocks can be selected individually for each chromosome.

[0255] In one possible implementation, the entire chromosome can be represented as columns in an array, and one of the two telomere blocks of all chromosomes can be aligned along a row of the array. In another possible implementation, the alignment row for the telomere blocks can be the first row of the array. In yet another possible implementation, the alignment row for the telomere blocks can be the last row of the array. In one possible implementation, p-arm telomere blocks can be used for alignment. In yet another possible implementation, q-arm telomere blocks can be used for alignment. In yet another possible implementation, the alignment of either the p-arm or q-arm telomere blocks can be selected individually for each chromosome.

[0256] The genomic data analyzer 120 computer system (also referred to herein as the "system") 120 can be programmed or otherwise configured to implement different genomic data analysis methods beyond those described herein for CSI and HRD analyzer systems and methods, such as receiving and / or combining sequencing data, invoking copy number alterations and / or annotation variants to further characterize tumor samples. Machine learning models can employ extended multidimensional array inputs, including information such as GC content, block size, mapping, mapping quality, and variant allele fractions (VAF), in addition to normalized coverage data in a two-dimensional array. Machine learning models can also employ additional scalar inputs, providing tumor content information (purity) and sample ploidy information. The genomic data analysis system can also be used to assess sample quality (i.e., measure the extent of FFPE degradation). The genomic data analysis system can also be used to classify cancers according to their CSI status level. The genomic data analyzer system can also be used to stratify cancers based on their CSI status and association with immune escape events (Bakhoum et al., 2018, Cell 174(6), p.1347–1360).

[0257] While exemplary implementations and applications of the proposed method have been described using CNNs as trained machine learning models, random forest machine learning models can be used interchangeably in possible implementations (not illustrated). Random forests, or random decision forests, are ensemble learning methods for classification, regression, and other tasks that operate by building a large number of decision trees during training and outputting patterns of classes (classification) or the average prediction of a single tree (regression). Random decision forests correct the tendency of decision trees to overfit their training set, making them a good candidate for high-dimensional input data.

[0258] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains. The terminology used in this description is for the purpose of describing particular embodiments only and is not intended to be limiting. In the specification and appended claims, unless the context clearly indicates otherwise, the singular forms “a,” “an,” and “the” also include the plural forms.

[0259] Unless otherwise stated, the numerical parameters listed in the following description and appended claims are approximate values ​​and may vary depending on the desired characteristics sought to be obtained, and are therefore subject to modification by the term "about". At least, and not to limit the application of equivalence theory within the scope of the claims, each numerical parameter should be interpreted in accordance with the number of significant figures and ordinary rounding methods.

[0260] While the ranges of numbers and parameters presented are approximate, the values ​​given in specific examples are reported as precisely as possible. However, any value will inevitably contain some error, caused by the standard deviation found in the respective test measurements. Each range of values ​​given in this specification will include every narrower range of values ​​falling within this wider range, as these narrower ranges are explicitly stated here.

[0261] It will be obvious to those skilled in the art of digital data communications that the methods described herein can be applied indiscriminately to various data structures, such as data files or data streams. Therefore, the terms “data,” “dataset,” “data structure,” “data field,” “file,” or “stream” can be used interchangeably in this specification.

[0262] Although the above detailed description contains many specific details, these should not be construed as limiting the scope of the invention, but merely as providing some illustrations of a few embodiments.

[0263] Furthermore, it should be understood that any diagrams highlighting features and advantages are for illustrative purposes only. The disclosed methods are flexible and configurable enough that they can be utilized in ways other than those shown.

Claims

1. A computer-based in vitro method for processing sequencing data of isolated subject DNA samples to generate a homologous recombination defect (HRD) status of said isolated subject DNA samples, said method comprising the following steps: (a) Extracting and separating DNA fragments from the isolated subject DNA sample; (b) Constructing a sequencing library containing DNA fragments from said isolated subject DNA samples that overlap with a set of chromosomes; (c) Sequencing the sequencing library to a read depth of at least 0.1X and at most 5X via whole-genome sequencing; and (d) Determining the HRD status of the isolated subject DNA sample, wherein determining the HPD status of the isolated subject DNA sample includes the following steps: (i) Obtain sequencing reads of the isolated subject DNA sample to be analyzed, wherein the sequencing reads are obtained by low-throughput whole-genome sequencing of the isolated subject DNA sample. (ii) The sequencing reads of the isolated subject DNA sample are compared with a reference genome, wherein the reference genome is divided into multiple blocks, each block belonging to the same genomic region of a chromosome arm in the reference genome to be analyzed. (iii) Count and normalize the numerical values ​​of the alignment readings in each block along each chromosome arm; (iv) Based on the steps of counting and normalizing the values ​​of the alignment readings in each block along each chromosome arm, the coverage signal of each chromosome arm is obtained. (v) Arranging the coverage signals of each chromosome arm into a coverage data signal array of the separated subject DNA sample, including arranging the coverage signals of each chromosome arm into a one-dimensional coverage data signal vector or a two-dimensional coverage data signal array; arranging the coverage signals of the chromosome arms into the two-dimensional coverage data signal array, including aligning the coverage data signals of each chromosome row by row to the centromere region of each chromosome arm, i.e., aligning to the region adjacent to the centromere region of the chromosome arm; (vi) The coverage data signal array is input into a trained machine learning model, wherein the model is trained using coverage signal arrays of samples with known homologous recombination defect states to distinguish coverage data signal arrays from samples with positive homologous recombination defect states and coverage data signal arrays from samples with negative homologous recombination defect states. (vii) The homologous recombination defect score (HRD score) of the isolated subject DNA sample is determined and output via the trained machine learning model. (viii) Determine the negative, positive, or indeterminate homologous recombination defect status of the isolated subject DNA sample based on the HRD score output by the trained machine learning model.

2. The method according to claim 1, wherein, Counting and normalizing the numerical values ​​of the alignment readings in each block along each chromosome arm to obtain the coverage signal of the chromosome arm includes normalizing the coverage signal of each sample and / or normalizing by GC content to apply GC-offset correction.

3. The method according to claim 1, wherein the plurality of samples with known homologous recombination defect states are tumor data samples with known homologous recombination defect states.

4. The method according to claim 3, wherein, The tumor data sample with a known homologous recombination defect state is augmented by artificial sample data generated by combining the chromosomal data of the tumor data sample with the known homologous recombination defect state label, forming a data augmentation sample.

5. The method according to claim 4, wherein, The data augmentation samples are generated to represent the purity-ploidy ratio distribution observed in the tumor data samples.

6. The method according to claim 1, wherein, The reference genome is divided into a first group of blocks of up to 100 kbp, and further includes the step of collapsing the 100 kbp blocks into a second group of larger blocks of at least 500 kbp before arranging coverage signals on each chromosome arm.

7. The method according to claim 6, wherein, The first group of blocks has a uniform size of up to 100kbp, while the second group of blocks has a size of 2.5 to 3.5Mbp and is obtained by pooling 25 to 35 100kbp blocks from the first group.

8. The method according to claim 1, wherein, The isolated subject DNA samples are cell-free extracellular DNA (cfDNA), fresh frozen tissue (FFT), or formalin-fixed paraffin-embedded samples (FFPE).

Images