Patent Application
20220084626
Methods described herein relate to genomic analysis in general, and more specifically to the use of genomic information for detecting and treating cancer.
Beyond inherent germline mutations, cancer cells often carry somatic gross chromosomal aberrations such as copy number variations, duplications or deletions of certain alleles or genomic regions. Some of these variants may cause a loss of some genomic functions such as for instance the tumor suppressor mechanisms, in particular the homologous recombination repair (HRR or HR) function, thus making the cancer more aggressive. Identifying those genomic mutations is of key importance in recent advances of personalized cancer medication, as they have been shown to be predictive of the response of a subject with a cell hyperproliferative disorder to certain anti-cancer therapies. The subject may be a human or an animal. Examples of anti-cancer treatments comprise such alkylating agents such as for instance, but not limited to, platinum-based chemotherapeutic agents, carboplatin, cisplatin, iproplatin, nedaplatin, oxaliplatin, picoplatin, chlormethine, chlorambucil, melphalan, cyclophosphamide, ifosfamide, estramustine, carmustine, lomustine, fotemustine, streptozocin, busulfan, pipobroman, procarbazine, dacarbazine, thiotepa, temozolomide and/or other antineoplastic platinum coordination compounds; a DNA damaging agent or radiation; an anthracycline, such as for instance, but not limited to, epirubincin or doxorubicin; a topoisomerase I inhibitor, such as for instance, but not limited to, campothecin, topotecan, irinotecan; and/or a PARP (poly ADP-ribose polymerase) inhibitor. Examples of a PARP inhibitors that exploits the tumor DNA homologous repair deficiency (HRD) to selectively destroy cancer cells are olaparib, rucaparib, niraparib (MK 4827), and talazoparib (BMN-673) which have been approved in the US and in Europe for certain cancer types; other examples of PARP inhibitors include iniparib, CEP 9722 (-14-methoxy-9-[(4-methylpiperazin-1-yl)methyl]-9,19-diazapentacyclo[10.7.0.0{circumflex over ( )}{2,6}.0{circumflex over ( )}{7,11}.0{circumflex over ( )}{13,18}]nonadeca-1(12),2(6),7(11),13(18),14,16-hexaene-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 multitude of DNA double strand breaks (DSBs) that prevents the development of HRD tumor cells while normal cells are generally able to repair those breaks via HR (Keung et al., 2019, Journal of Clinical Medicine 8 (4), p. 435).
Accordingly, a method of classifying DNA samples, such as samples from a tumor, may therefore facilitate the diagnosis of possible cancer types or may allow to adapt the most suitable cancer treatment to patients in accordance with the characterization of gross DNA copy number aberrations in their tumor samples thanks to genomic sequencing and analysis. In particular, identification of whether a cancer is homologous repair deficiency (HRD) may be of considerable assistance in the planning of treatment. In the past decade, different genomic mutational signatures have been identified to this end. European patent EP2609216B1 discloses the use of a global chromosome aberration score (GCAS) to predict the outcome of an anti-cancer treatment such as a PARP inhibitor, radiation therapy or a chemotherapy comprising platinum-based chemotherapeutic agents. European patent EP2817630B1 proposes the detection of the number of telomeric allelic imbalance (TAI) events and the selection of a platinum-comprising therapy if this number is above a reference value from similar cancers known to be platinum-resistant. European patent EP2859118B1 by Stern, Manie and Popova from Institut Curie and INSERM discloses a method to predict HRD by counting the number, per genome, of segments spanning at least 3 megabases, which corresponds to large-scale transitions (LST) in chromosomal copy number. European patent EP2794907B1 by Abkevich et al. from Myriad Genetics discloses counting the total number of loss-of-heterozygosity (LOH) regions that are longer than 11 megabases but shorter than the whole chromosome in at least one pair of human chromosomes and comparing this number to a reference number to predict the patient response to a diverse possible cancer treatments; EP2981624B1 discloses using the LOH, TAI and LST indicators and EP3180447B1 discloses summing their values into an HRD score to be compared to a reference number in order to predict the patient response to a diverse possible cancer treatments. The latter test method is currently used as a test for HRD in ovarian cancer patients by the Myriad Genetics myChoice CDx assay, which uses a custom hybridization-based target-enrichment on a number of genes at a median coverage of at least 500×. This HRD score was also recently shown as a possible prognosis predictor regarding the progression of 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).
Similar to the Myriad myChoice assay, the Foundation Medicine commercial assay FoundationFocus CDxBRCALOH specializes on the analysis of the BRCA1 and BRCA2 genes to determine whether an ovarian cancer patient will be responsive to responsive to rucaparib, a PARP inhibitor therapy. The latter assay also uses a custom hybridization-based capture of all coding exons of a number of genes at a median coverage of 500×, yet solely uses an LOH score assessed from the combination of a genome-wide copy number profile and minor allele fractions of SNPs.
More generally, WO2017191074 proposes to characterize the HRD status of a tumor DNA sample according to a probabilistic score calculated from the analysis of different genomic alterations, including base substitutions, rearrangements and indel signatures (HRDetect score). As described in “HRDetect is a predictor of BRCA1 and BRCA2 deficiency based on mutational signatures”, H. Davies et al., Nature Medicine, published online Mar. 13, 2017, with the latter HRDetect score, Whole Genome Sequencing (WGS) is shown as a possible method to detect HRD, with a 98.7% sensitivity, from mutational signatures observed throughout the whole genome. In contrast to WGS, when only Whole Exome Sequencing (WES) is applied, the sensitivity of HRDetect significantly drops, between 46.8% and 73% depending on specific adaptations of the bioinformatics algorithms.
Beyond breast and ovarian cancer, as reported in “Genomic aberration based molecular signatures efficiently characterize homologous recombination deficiency in prostate cancer”, Sztupinszki et al. also recently started to investigate the use of the LOH, TAI and LST as indicators of HRD signatures in prostate tumors out of WGS and WES data (scarHRD).
Since all the above methods rely on SNV and/or INDEL calling, they require a high depth of coverage (typically at least 30×) in next generation sequencing (NGS) workflows. High coverage requirements in WGS and large panels significantly increases the cost of the analysis in clinical practice, both in terms of wet lab experiments and dry lab data processing overhead. Alternative methods using low-pass WGS (LP-WGS—1× to 5×) or even ultra-low-pass WGS (ULP-WGS—down to 0.1×) may thus 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. describe ShallowHRD, a software method to characterize the LST status similarly to the method of Manie, Stern and Pova, yet specifically on shallow WGS data, down to a coverage of about 1×. Their approach simply uses the count of intra-chromosome arm copy number alterations. This count may be estimated as the number of large-scale transitions between adjacent segments of at least 10 megabases after removing segments of less than 3 Mb in the coverage data signal. According to the authors, less than 15 transitions enables to identify HRD negative tumors while more than 19 enables to identify HRD positive tumors with similar sensitivity and specificity as the prior art scarHRD method. While this method offers the promise for the use of more cost-effective NGS workflows in clinical oncology, we observe that its sensitivity however still remains lower than the results achieved by HRDetect and the measurement of the LST indicator on SNP arrays. This may be due to the fact that it cannot integrate the LOH and TAI signatures of the prior art methods. Moreover, it only enables to classify tumors for which the chromosomal alterations are large and frequent enough to cause more than a dozen large scale transitions, regardless of how those transitions are distributed within each chromosome arm on the one hand and among the multiple chromosomes of the human genome on the other hand.
A number of recent research works in oncology have suggested that the hypermutability of dysfunctional DNA may cause heterogenous alterations in different areas of the chromosomes. In “Regulation of mitotic recombination between DNA repeats in centromeres”, Nucleic Acids Research, 2017, Vol. 45, No. 19, Zafar et al. 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. report that a significant ratio of chromosome rearrangements and breaks are observed in pericentromeric and centromeric regions in cancer cell lines derived from colorectal carcinomas and adenocarcinomas, probably due to an inherent fragility of the centromere regions in tumors and this is also common in other cancers. Barra et al. highlight the lack of models and technologies to explicitly analyze the highly repetitive sequences in those regions.
There is therefore a need for improved genomic analysis methods, which can be easily and cost-effectively deployed in automated NGS workflows for routine oncology clinical practice while being adaptive to recent cancer genomics discoveries and in particular to the promising discoveries of certain spatial features of chromosomal aberrations in the tumour genomes. There is a particular need for improved genomic analysis methods that integrate recent cancer genomics discoveries and allows to identify features, including HRD, relevant for diagnostic, prognostic, treatment selection and patient management. Preferably, those methods and technologies may also be applicable at a lower coverage in order to decrease the overall analysis cost in routine clinical setups, while preserving a reasonable accuracy in the detection of HR events compared to prior art techniques. Moreover, improved genomic analysis methods may also facilitate the characterization of HR-deficient tumours by employing machine learning technologies instead of relying upon explicit measurement and thresholding of scalar HRD indicators such as the LOH, LST or TAI indicators of the prior art, so that they can more easily exploit the increasing clinical data availability in the classification of tumours for a diversity of different cancers in humans and animals.
Disclosed herein is a method for determining a homologous recombination deficiency (HRD) status of a subject DNA sample, the method comprising obtaining a set of sequencing reads of the whole genome of the subject DNA sample to be analyzed, aligning the set of sequencing reads of the subject DNA sample to a reference genome, wherein the reference genome is divided into a plurality of bins, each bin belonging to the same genomic region from a chromosome arm in the whole genome chromosomes to be analyzed, counting and normalizing the number of aligned reads in each bin along each chromosome arm to obtain a coverage signal on the chromosome arm, arranging the coverage signals of the chromosome arms into a coverage data signal array for the subject DNA sample, inputting the coverage data signal array to a trained machine learning model, wherein the model has been trained using a set of samples of known homologous recombination deficiency status to distinguish between the coverage data signal array from samples with a positive homologous recombination deficiency status and the coverage data signal array from samples with a negative homologous recombination deficiency status, thereby determining a homologous recombination deficiency score (HRD score) of the subject DNA sample, and determining a negative, a positive or an uncertain homologous recombination deficiency (HRD) status of the subject DNA sample according to the HRD score out of the trained machine learning model. In a possible embodiment, the set of sequencing reads may be obtained from a whole genome sequencing wherein a read depth coverage is at most 30×, or from a low-pass whole genome sequencing wherein a read depth coverage is at least 0.1× and at most 5×. In a possible embodiment, counting and normalizing the number of aligned reads in each bin along each chromosome arm to obtain a coverage signal on the chromosome arm may comprises normalizing the coverage signal per sample, and/or normalizing by GC content to apply a GC-bias correction. In a possible embodiment, the coverage signals of the chromosome arms may be arranged into a 1D coverage data signal vector or a 2D coverage data signal image. In a possible embodiment, the coverage signals of the chromosome arms may be arranged into a 2D coverage data signal image by aligning in rows the coverage data signal for each chromosome with respect to the centromeric bin of each chromosome arm, that is the closest bin adjacent to the centromere region of the chromosome arm. In a possible embodiment, the machine learning model may have been previously trained using a set of tumor data samples with a known homologous recombination deficiency status as training label. The training dataset may be augmented with artificial sample data generated by combining data from chromosomes of data samples with a known homologous recombination deficiency status label. The data augmented samples may be generated in order to represent the purity-ploidy ratio distribution as observed in the real samples dataset. In a possible embodiment, the reference genome may be divided into a first set of at most 100 kbp bins and further comprising a step of collapsing the 100 kbp bins into a second set of larger bins of at least 500 kbp prior to arranging the coverage signals on each chromosome arm. The bins of the first set of bins may have a uniform size of at most 100 kbp and the bins of the second set of bins may have a size of between 2.5 to 3.5 Mbp and are obtained by pooling between 25 to 35 of 100 kbp bins from the first set of bins.
An in vitro method is provided for determining a homologous recombination deficiency (HRD) status of a patient DNA sample, the method comprising:
A method is provided for selecting a cancer patient for treatment with a platinum-based chemotherapeutic agent, a DNA damaging agent, an anthracycline, a topoisomerase I inhibitor, a PARP inhibitor comprising a step of detecting that a tumor patient sample is HRD positive according to the methods disclosed herein. In a possible embodiment, the patient may have a cancer selected from high grade serous ovarian cancer, prostate cancer, breast cancer or pancreatic cancer.
A method is provided for training a machine learning algorithm for determining the homologous recombination deficiency (HRD) status of a subject DNA sample, the method comprising inputting to a machine learning supervised training algorithm a coverage data signal array from samples with known positive homologous recombination deficiency status and a coverage data signal array from samples with known negative homologous recombination deficiency status.
The trained machine learning model may be a random forest model, a neural network model, a deep learning classifier or a convolutional neural network model. The neural network model trained machine learning model may be a convolutional neural network trained to produce at its output a single label binary classification of the positive or negative HRD status, or a single label multiclass classification of the positive, negative or uncertain HRD status, or a scalar HRD score representative of the HRD status. The machine learning model may be trained in semi-supervised mode using a data augmented sets generated by sampling data from chromosomes of a set of real samples sharing the same HRD status and the same normalized purity and ploidy ratios.
A method is proposed for characterizing a DNA sample, the method comprising obtaining a set of sequencing reads from a patient sample aligned to said reference genome; dividing the base pair positions (bp) for at least two chromosomes of a reference genome into a set of bins, each bin corresponding to a genomic region of at most 20 megabase pairs (20 Mbp), each bin solely containing coverage data from a single chromosome arm; estimating, from the aligned reads, the coverage data for each bin of the sequenced patient sample; arranging the coverage data into a multidimensional array comprising either the chromosomes or the chromosome arms along one dimension and the set of bins for said chromosome or chromosome arm along another dimension, characterized in that either the centromeric bins or the telomeric bins for each chromosome or chromosome arm are aligned into the multidimensional array spatial arrangement; inputting the multidimensional array to a trained machine learning model; and generating, at the output of the trained machine learning model, a chromosomal spatial instability (CSI) indicator. In a possible embodiment, the CSI indicator may be an indicator of whether said patient sample has a high (HRD+) or low (HRD−) likelihood of being homologous recombination (HR)-deficient. The patient sample may be a tumor sample and the CSI indicator of the patient sample may be a predictor of the tumor response to a cancer treatment regimen comprising a platinum-based chemotherapeutic agent, a DNA damaging agent, an anthracycline, a topoisomerase I inhibitor, or a PARP inhibitor.
In one embodiment, a machine learning model is trained in supervised or semi-supervised mode using labelled coverage data from samples with a known genomic instability status, such as HRD status (positive or negative).
In one embodiment, a chromosomal spatial instability (CSI) indicator of genomic instability is based on the difference between a DNA sample with a genomic instability in at least one region of at least one of the chromosome arms in said DNA sample and another DNA sample without a genomic instability in at least one region of at least one of the chromosome arms in said DNA sample, as may be learnt by a machine learning model.
In a possible embodiment, the first set of bins may be further collapsed into larger bins prior to arranging the multidimensional array, by adapting the larger bin size with respect to each chromosome arm length such so that the bin size remains constant along the chromosome arm. The first set of bins may have a uniform size of at most 100 kbp while the collapsed bins may have a size of at least 500 kbp. The coverage data may be normalized and/or may be segmented to infer discrete copy number values prior to arranging the multidimensional array. In a possible embodiment, the empty elements of the multidimensional array may be padded with either the value of the closest bin or with a predefined value. In a possible embodiment, the trained machine learning model may be a random forest model or a neural network model, such as for instance a convolutional neural network trained to produce at its output a single label binary classification of the positive or negative CSI status, or a single label multiclass classification of the positive, negative or uncertain CSI status, or a scalar CSI score representative of the CSI status. The patient sample may be a tumor sample and the machine learning model may be trained in supervised or semi-supervised mode using a set of real samples labelled with a CSI status according to the target application. In a possible embodiment, the CSI status label may be an HRD status according to the HRDetect method and/or using artificial samples generated by sampling chromosomes of a set of real samples sharing the same HRD status and having similar tumor content. The artificial samples may be chosen in order to recreate the same purity-ploidy distribution as with the real samples dataset.
In one particular embodiment, the patient sample may be a tumor sample and the machine learning model may be trained in supervised or semi-supervised mode using labelled data with an HRD status. Training data can be composed of coverage data obtained from clinical samples with an HRD status and arranged as multidimensional arrays. Training data labels can be obtained from HRD detection methods for example HRDetect method.
In one particular embodiment, data augmentation strategies that account for biases introduced by sample specific properties, such as purity and ploidy, can be used to increase the diversity and number of multidimensional arrays used to train the machine learning algorithm.
Color of normalized coverage points in original, diluted and artificial samples was preserved.
The present disclosure is based, at least in part, on the discovery that the presently disclosed machine learning trained analyzer, designed to process a tumor sample sequencing data coverage in relation with the underlying chromosomal arrangements, can extract an indicator of a chromosomal spatial instability (CSI) in general and in particular an indicator of a homologous repair deficiency (HRD) status for the tumor sample.
The proposed methods and systems will now be described by reference to more detailed embodiments. The proposed methods and systems may, however, be embodied in different forms and should not be construed as limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope to those skilled in the art.
The proposed methods and systems will now be described by rAn exemplary genomic analysis system and workflow will now be described with further detail with reference to
As illustrated on
The Genomic Data Analyzer 120 may process the sequencing data to produce a genomic data analysis report by employing and combining different data processing methods.
The sequence alignment module 121 may be configured to execute different alignment algorithms. Standard raw data alignment algorithms such as Bowtie2 or BWA that have been optimized for fast processing of numerous genomic data sequencing reads may be used, but other embodiments are also possible. The alignment results may be represented as one or several files in BAM or SAM format, as known to those skilled in the bioinformatics art, but other formats may also be used, for instance compressed formats or formats optimized for order-preserving encryption, depending on the genomic data analyzer 120 requirements for storage optimization and/or genomic data privacy enforcement.
The Genomic 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 a hard disk, and communication interfaces to communicate with other computer systems through a communication network, for instance the internet or a local network. Examples of genomic 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, minicomputer systems, mainframe computer systems, graphical processing units (GPU), and the like. In some embodiments, the computer system may comprise one or more computer servers, which are operational with numerous other general purpose or special purpose computing systems and may enable distributed computing, such as cloud computing, for instance in a genomic data farm. In some embodiments, the genomic data analyzer 120 may be integrated into a massively parallel system. In some embodiments, the genomic data analyzer 120 may be directly integrated into a next generation sequencing system.
The Genomic Data Analyzer 120 computer system may be adapted in the general context of computer system-executable instructions, such as program modules, being executed by a computer system. Generally, program modules may include routines, programs, objects, components, logic, data structures, and so on that perform particular tasks or implement particular abstract data types. As is well known to those skilled in the art of computer programming, program modules may use native operating system and/or file system functions, standalone applications; browser or application plugins, applets, etc.; commercial or open source libraries and/or library tools as may be programmed in Python, Biopython, C/C++, or other programming languages; custom scripts, such as Perl or Bioperl scripts.
Instructions may be executed in distributed cloud computing environments where tasks are performed by remote processing devices that are linked through a communications network. In a distributed cloud-computing environment, program modules may be located in both local and remote computer system storage media including memory storage devices.
It is thus understood that methods described herein are computer-implemented methods.
In a possible embodiment, the reference genome coordinates are first divided into a set of P bins and the number of DNA fragments mapping to each bin is counted taking as input the aligned reads from the BAM or SAM file (after sequencing read alignment), using bioinformatics methods as known in the art, to produce 200 the coverage data signal along the P bins. As will be apparent to those skilled in the art of bioinformatics, in some embodiments, some of the bins may be ignored, filtered, or merged with other bins to improve the input data quality.
In a possible simple embodiment, a uniform size may be used for the bins along the reference genome, but a choice of heterogeneous sizes is also possible. The size of the bins may be chosen such that it preserves enough spatial detail for the analysis of the coverage data signal for each chromosome while facilitating the coverage data signal analysis with the machine learning analysis methods described herein. Preferably, the size of the bins may be in the range from 300 bp to 20 Mpb, but other embodiments are also possible.
In the case of the reference human genome, chromosomes are sequentially ordered from the largest chromosome (chr1) to the smallest one (chr22), possibly with the exclusion of the sexual chromosomes. Human chromosomes are numbered by decreasing length, from 249 Mbp for chr1 down to 51 Mbp for chr22. By dividing the reference genome in bins of 3 Mbp for instance, it is therefore possible to obtain a 1D array (vector) of 84 bins for the longest chromosome 1 down to 4 bins for the shortest chromosome 22. Table 1 shows the theoretical respective start and end position of the centromeric region, the total chromosome length and its type, as well as the respective lengths in Mbp for the p-arm (short arm) and the q-arm (long arm) for each chromosome in the human genome, with reference to the GRCh37 reference genome (hs37d5 version), in the coordinate system of the reference genome used for the read alignment. Note that not all bins of the human genome may contain aligned reads in the BAM file with the current high throughput sequencing technologies, due to the difficulty to accurately mapping the short NGS reads to certain repetitive sequences or problematic regions of the human genome, for instance the p-arm of the acrocentric chromosomes 13, 14, 15, 21, 22 or Y.
As will be apparent to those skilled in the arts of bioinformatics, in a possible embodiment, the bin sizes and positions along the whole genome sequencing data may also be adapted to the actual coverage from the BAM file. In general, we may consider that each bin in the set of P bins may be defined by 2 values corresponding to the start and end positions of the bin {bin_start, bin_end}1<=bin_index<=P in the reference coordinate system. The size of the bin may thus be variable as bin_size=bin_end−bin_start. In a possible embodiment, the size of the bin may be in the range from 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 embodiment, the size of the bins may be selected such that a fixed number of bins may be obtained for each chromosome or for each chromosome arm. Other embodiments are also possible.
A coverage signal for each chromosome or for each chromosome arm may be then obtained by counting the number of aligned reads in each bin of the fixed set of bins. In a possible embodiment, the resulting high resolution coverage data signal array for each chromosome or for each chromosome arm may be used as a data input of a high dimensionality for the downstream CSI analysis. In a possible alternate embodiment, a first size of 100 kbp may be used to produce a first set of high-resolution coverage bins. The resulting high resolution coverage data signal may be further collapsed into a second set of larger bins of 1 to 20 Mbp, for instance by combining about 30 100 kbp-large initial bins to produce a bin of a larger size, between 2.5 Mbp and 3.5 Mbp, thus reducing the dimensionality of the data for the downstream CSI analysis. Various methods known to those skilled in the art of bioinformatics may be used to this end, such as smoothing, using the median or the mean to reduce the dimensionality, pooling, sampling, and other methods. In a preferred embodiment, collapsing uniquely assigns each bin to a specific chromosome arm without mixing coverage data from different chromosome arms into the bin. The coverage data is then clearly separated from one chromosome arm to the next one, which may facilitate the analysis of chromosomal spatial instability events at the chromosome arm level.
In a possible embodiment, collapsing may comprise adapting the size of the collapsed bins to each chromosome arm length in order to reduce the border effect at the chromosome arm boundaries. For instance, in a situation where the chromosome arm length is not dividable by the initial bin size considered, the last collapsed bin obtained will usually be shorter than the other ones, causing a border effect (for instance the centromeric bin or the telomeric bins on any end of the chromosome arm may not comprise an entire area of 3 Mbp, and may result in a distorted coverage for those areas while possibly key in analyzing CSI events). A strategy to reduce this effect is to adapt the binning size with respect to chromosome arm length such that the chromosome arm length considered is dividable by the binning size with a neglectable remainder. This ensures that all the collapsed bins over the chromosome arm considered will have the same actual size.
In a possible further embodiment, collapsing may comprise adapting the start or the end position, and thus the actual size, of the collapsed bins individually along a chromosome arm such that the sum of the variable sizes of the collapsed bins still equals the chromosome arm length while being more or less evenly distributed along the chromosome arm.
For example, the initial 100 kbp bins from each chromosome arm may be smoothed independently to obtain bins of 3 Mbp (i.e aggregation over 30 bins of 100 kbp). As the total number of bins in each chromosome arm may not be exactly dividable by 30, the last smoothed bin of chromosome arms usually contains less than 30 original bins. This may either give those regions too much importance if keeping them, or conversely if removing them, losing too much information. To reduce the border effect, larger bins of variable size may be collapsed for each chromosome arms independently. In a possible embodiment, the bin size that leads to the smallest non-complete bin remainder and that is closest to the target size may be selected. In a possible embodiment, the last bin if not complete may then be rejected as now neglectable, but other embodiments are also possible. Table 2 column 3 shows an exemplary distribution of the proposed collapsed bin sizes for each chromosome arm for an initial set of 100 kbp bins collapsed into larger bins at a target 3 Mbp size along the 2*22 chromosome arms. Column 4 shows the resulting sizes of the rejected border bins, compared to the default number, in column 2, of remainder higher resolution bins that would form the last collapsed bin at the arm border when dividing by a fixed size of 3 Mpb for all arms.
As will be apparent to those skilled in the art of bioinformatics, in some embodiments, the raw coverage data signal may be further normalized, filtered or smoothed by a method known in the art. In a possible embodiment, normalization may comprise normalizing per sample, by dividing the coverage data signal by the mean coverage signal for the whole sample, in order to get rid of biases associated with the sample sequencing experiment possibly causing differences in the coverage per sample. Normalization may also comprise normalization by GC content to apply a GC-bias correction. Other embodiments are also possible, such as bin-wise normalization, linear and nonlinear least squares regression, GC LOESS (GC normalization), LOWESS, PERUN (normalization per sample), RM, GCRM, cQn and/or combinations thereof.
In a possible further embodiment, the normalized data may be further segmented into sub-regions within the bins to identify homogeneous areas within the coverage data signal, possibly corresponding to different copy numbers (normally 2 per chromosome). Coverage segmentation allows to infer copy number at discrete segments relative to the other genomic regions from the coverage profile of a given sample. To do so, the coverage profile is decomposed into segments (i.e., portion of the genome) within which the coverage is believed to be constant. Each of the segments defined in the previous step may then be associated to discrete levels defined in terms of coverage signal amplitude. At this stage, Copy Number Variation (CNV) events can already be detected but the absolute Copy Number associated to each coverage level remains unknown.
In cases where the sample tumor content and ploidy are known or can be inferred from the data, the discrete levels identified by segmentation can be mapped to integer numbers reflecting the number of copies present in the tumor (e.g, CN=1, CN=2, . . . ). As will be apparent to those skilled in the art of bioinformatics, such segmentation methods with or without an absolute copy number reference may reduce the noise and facilitate the identification of the copy number profile for the sample. However, in the presence of noise, segmentation may produce erroneous results by suppressing coverage fluctuations which results from changes in copy number, in particular when sequencing at low coverage depths. In a preferred embodiment the genomic analysis method as described herein may therefore operate without the segmentation step. In an alternate embodiment, a segmentation step may be applied to pre-process the data prior to applying the machine learning methods as disclosed herein.
After possible steps of normalization and/or segmentation, the coverage data generally consists of a set of 1D vectors, each of them corresponding to a chromosome. In a possible embodiment, the coverage data may also consist of a set of 1D vectors, each of them corresponding to a chromosome arm. In a possible embodiment, the size of the vectors may be variable, reflecting the fact that different chromosomes have different sizes. In another possible embodiment, the bins have a variable size so that each chromosome vector (or each chromosome arm vector) has the same size.
An exemplary coverage data 1D signal is plotted on
In a preferred embodiment, the coverage data may therefore be further spatially arranged 210 to form a multidimensional data array in which the coverage data is organized to align the bins for each chromosome or chromosome arm along any of their telomeric or centromeric bins in the array structure. A multidimensional array structure may be particularly well suited to benefit from efficient genomic analysis machine learning trained models recently developed for non-genomic applications such as image classifiers, instead of relying upon simple human engineered features as in the prior art HRD detection methods.
In other words, and in line with provided definition of telomeric or centromeric bins, the coverage data signal vectors may be further spatially arranged 210 to form a coverage data signal image in which the coverage data is organized to align the bins for each chromosome or chromosome arm along at least one centromeric bin adjacent to a centromere region or at least one telomeric bin adjacent to a telomere region in a 2D array structure.
In a possible embodiment, the coverage data signal vectors may be further spatially arranged 210 to form a 1D coverage data array in which the coverage data is organized as the tail to head concatenation of the bins for each chromosome or chromosome arm along the genome.
In another preferred embodiment, the coverage data signal vectors may be further spatially arranged 210 to form a multidimensional coverage data signal array in which the coverage data is organized to align the bins for each chromosome or chromosome arm along at least two bins adjacent to a centromere region in the array structure, wherein one bin is on the left of a centromere, and another bin is on the right of the centromere and wherein this multidimensional data array is a 2D array, and specifically a 2D image.
In another preferred embodiment, the coverage data may be further spatially arranged 210 to form a multidimensional data array in which the coverage data is organized to align the bins for each chromosome or chromosome arm along at least two bins adjacent to telomeric regions in the array structure, wherein one bin is at the start of the chromosome within the set of coverage bins associated with this chromosome, and another bin is at the end of the chromosome within the set of coverage bins associated with this chromosome and wherein this multidimensional data array is a 2D array, and specifically a 2D image.
In a further preferred embodiment, the coverage data may be further spatially arranged 210 to form a multidimensional data array in which the coverage data is organized to align the bins for each chromosome or chromosome arm along at least one or at least two bins adjacent to a centromere region or at least one or at least two bins adjacent to a telomere region in the array structure, wherein this multidimensional data array is a 2D array, and specifically a 2D image.
In this context a bin adjacent to a centromere or telomere regions means that it is a bin that is the closest to these regions.
In a further preferred embodiment, the bins for each chromosome arm are computed independently.
In one embodiment, a method may comprise, but is not limited to, a step wherein for each of the chromosome arm the coverage data may be further normalized by taking the median normalized coverage over a window of 2.5 Mbp to 3.5 Mbp. The exact bin resolution may be chosen between those two values independently for each chromosome arm such that a smallest possible non-complete bin at the end of the chromosome arm is obtained. The last non-complete bin is discarded from the input of the machine learning model, such as CNN. Therefore, the last bin for each p-arm contains coverage signal that is the closest regarding its position to its corresponding chromosome centromere. Similarly, the first bin for each q-arm is the one that contains coverage signal that is the closest regarding its position to its corresponding chromosome centromere. In parallel, the first bin for each p-arm contains coverage signal that is the closest regarding its position to its corresponding p-arm chromosome telomere region. Similarly, the last bin for each q-arm is the one that contains coverage signal that is the closest regarding its position to its corresponding q-arm chromosome telomere region.
It is thus understood and described in
The effect of organizing of the coverage data to align the bins as described herein allows to detect a spatial relation between detected events, such as genomic instability within chromosome arms. The machine learning model, such as CNN, may use information on spatial distribution of detected events and thus potentially improve its output prediction when accurate status detection benefits from this type of spatial arrangement.
As the number of bins changes with the length of the chromosome and/or the chromosome arms, some elements of the array need to be filled in by padding virtual data not available from the coverage data. As will be apparent to those skilled in the art of data science, different options are possible to pad the empty entries (or elements) of the multidimensional array. In a possible embodiment, the empty elements may be filled by a constant value. In another possible embodiment, the empty elements may be filled from a mask array of predefined constant values. In another possible embodiment, the empty elements may be filled as a function of the coverage data in the part of the row or the column which has been filled by real data. In a possible embodiment, the empty elements may be padded by repeating the coverage value left or right (resp. above or below) from the last available bin in the column (resp. row) of bins along a chromosome or chromosome arm row (resp. column).
In another possible embodiment, the individual size of the bins may be selected so as to facilitate a denser filling of the array. Smaller bin sizes (corresponding to a better resolution) may be selected in areas of particular genomic interest to increase the total number of bins for some chromosome arms.
The Genomic Data Analyzer system 250 may further apply a chromosomal spatial instability analyzer 123 to the spatially arranged coverage data multidimensional array to automatically analyze and report one or more genomic properties characterizing certain spatial instability patterns for two of more chromosomes or chromosome arms in each DNA sample. In the case of a tumor sample analysis, the CSI analyzer module 123 may report the burden of chromosomal aberrations such as large deletions or duplications of an allele in a given genomic region along each chromosome arm.
In a possible embodiment, the CSI analyzer module 123 is a computer-implemented algorithm trained to discriminate between a DNA sample with a genomic instability, labeled with an HRD positive status, and another DNA sample without a genomic instability, labeled with an HRD negative status.
In a possible embodiment, in the case of a tumor sample, the CSI analyzer module 123 may identify a signature which may identify samples that share a similar features and may thus respond well to a given treatment.
The CSI analyzer module 123 may accordingly identify, categorize and report an indicator, a scalar score or a combination of feature indicators to characterize one or more genomic properties of the tumor sample. In a preferred embodiment, the CSI analyzer module 123 may accordingly identify, categorize and report a status of the tumor sample, such as HRD-negative or HRD-positive or, possibly, undetermined (uncertain). In a possible embodiment, the CSI analyzer module 123 may report a scalar score as an indicator of the HRD likelihood for the tumor sample. As will be apparent to those skilled in the art, this bioinformatics method then significantly facilitates detailed understanding of the cancer cells genomic alterations in the tumor sample and the adaptation of a personalized medicine treatment to the specifics of the inferred cancer cell biology for the patient.
In one embodiment, the genomic property report may be displayed to the end user on a graphical user interface. In another possible embodiment, the genomic property report may be produced as a text file for further automated processing. Other embodiments are also possible.
In a preferred embodiment, the CSI analyzer module 123 may be adapted to apply a machine learning model to the spatially arranged coverage data multidimensional array. In a possible embodiment, the CSI analyzer module 123 may comprise a trained neural network or a combination of trained neural networks as a data processing module to analyze 220 the spatially arranged coverage data multidimensional array as the trained neural network input and to produce a CSI result as the trained neural network output.
An exemplary CSI Analyzer 123 and its data processing workflow will now be described with further detail with reference to
In order to classify a multidimensional array such as a 2D image in which chromosome arms have been aligned along their respective centromeric positions to highlight visual patterns which characterize the CSI status of this chromosome set, a convolutional neural network (CNN) may be architected and trained by employing methods and technologies known in the art of image pattern recognition such as TensorFlow, Caffe, Caffe2, Pytorch, Theano, or others.
More generally, as will be apparent to those skilled in the art of deep learning, a CSI Analyzer 123 may be adapted to employ various CNN architectures. The convolutional network may comprise a first series of serially connected convolutional data processing layers 810, which may further feed a second series of fully connected layers 720 to output various predictions. The CNN architecture may for instance comprise one or more convolutional layers arranged to process 1D, 2D or multidimensional data; optionally, one or more maximum-, median- or average-pooling layers arranged to process 1D, 2D or multidimensional data; one or more multiple dropout layers, to facilitate regularization during training; optionally, one or more batch normalization layers; a flattening layer followed by one or more fully connected layers; or a combination thereof. Different activation methods may be used as known in the art, such as for instance ReLU, SELU, or more generally any activation function that prevents vanishing gradient. Different output activation may also be used as known in the art, for instance a sigmoid method, to produce the probability associated to a single label binary classification (for instance, probability of positive status, probability of negative status=1—probability of positive status), or a softmax method to produce the probabilities associated to a single label multiclass classification output (for instance, mutually exclusive probabilities of a positive status, negative status, or an uncertain status), or a regression method such a linear method or the ReLU classification method to produce a scalar value which is not meant to be interpreted as a probability, such as for instance the HRD score.
In a preferred embodiment, the CSI Analyzer 123 produces a final output decision such as an HRD+ or an HRD− status analysis for each sample, by thresholding the output from the CNN. The threshold may be adapted depending on the application, so as the optimize the sensitivity and/or the specificity of the CSI analysis according to the end user needs, for instance for diagnosis, prognosis, or to guide the choice of the most efficient cancer treatment. In a possible embodiment, the CSI Analyzer 123 may also further report 230 the CNN output value(s) in accordance with the end user needs.
For supervised and semi-supervised CSI Analyzer trainings, a set of labelled input-output pairs may be used as a training set. For instance, to train the CSI Analyzer to classify the HRD status, a subset of the public domain HRD labelled data samples may be used to enable the network to minimize its output errors by adjusting its parameters (e.g., weights) along its different layers through an optimization process such as back propagation which is specifically applied during the training phase. As will be apparent to those skilled in the art of machine learning, different loss functions may be used to measure the error, for instance a binary cross entropy measurement for a model with a single label binary classification output, a categorical cross entropy measurement for a model with a single label multiclass classification output, or a regression method based on the root mean square error, the absolute error, or other error measurements. Optimization may employ state of the art stochastic descent gradient based optimizers, such as ADAM or RMPSprop. Regularization may employ a dropout method, an early-stop method, and/or an L-1 or L-2 parameter regularization. Other embodiments are also possible.
As more clinical conditions linked to chromosomal instability get classified by current and future oncology research, additional datasets may be used for training machine learning models aimed at detecting genomic scars induced by HRD or other phenomena such a response to a given treatment. In a possible embodiment, the HRD+ and HRD− binary statuses, and optionally an “uncertain” status may used as multiclass ground truth labels. In another possible embodiment, the HRD score from a prior art method such as HRDetect may be used as a scalar output label. The input coverage signal data array may also be spatially arranged as a multidimensional data according to any of the possible embodiments of the methods herein disclosed, as long as the choice of the multidimensional array format embodiment remains the same at training time and at runtime.
In a possible embodiment, additional data sets for training machine learning model may include a coverage data signal array from patient DNA samples, wherein these patients are of known homologous recombination deficiency status, and/or wherein these patients received a cancer treatment regimen and the outcome of the treatment is known (response vs. no response).
In a possible embodiment, additional datasets for training machine learning model may include a coverage data signal array from DNA samples generated from primary or immortalized cell line experiments or tumor-derived organoid experiments, especially experiments that determine whether organoids predicted to have HRD are eliminated or have their growth reduced by a cancer treatment regimen.
In a possible embodiment, a semi-supervised learning framework may be employed for training the machine learning model. Indeed, access to labeled data can be challenging and costly, particularly in the context of HRD status prediction which makes semi-supervised learning a promising framework for this application. Since the degree of model overfitting is determined by both its complexity and the amount of training it receives, providing a model with more training examples can help reduce overfitting. Data augmentation consists in the creation of new artificial training samples from existing data. In a possible embodiment, artificial samples may be generated by sampling chromosomes from available real samples sharing the same HRD status and having similar tumor content (or bionomically preprocessed to mimic similar tumor content) and ploidy.
In a possible embodiment, artificial samples may be generated by sampling chromosomes from available real samples sharing the same HRD status and the same normalized purity and ploidy ratios, wherein purity is the percentage of tumor cells in the sample and ploidy is the average chromosome copy number across the genome, thus characterizing how many complete chromosome sets on average tumor cells have.
In one embodiment, a machine learning model is trained on the coverage data signal array data from labelled real samples.
In another embodiment, a machine learning model is trained on the coverage data signal array data from both labelled real samples and labelled data augmented samples generated by sampling chromosomes from available real samples sharing the same HRD status and ensuring that the purity/ploidy ratios of sampled chromosomes is the same.
Any machine learning model may thus learn to quantify as a CSI score a difference between a DNA sample with a positive HRD status and another DNA sample a negative HRD status and then to classify tested DNA samples based on the CSI score result.
In a possible embodiment, the machine learning model has been trained in supervised or semi-supervised mode using at least a set of real samples labelled with an HRD status according to the HRDetect method. Preferably, the learning model has been trained in supervised mode using at least a set of real samples labelled with an HRD status according to the HRDetect method.
In a possible embodiment, the machine learning model has been trained in supervised or semi-supervised mode using artificial samples generated by sampling chromosomes of a set of real samples sharing the same HRD status and accounting for differences in purity ploidy ratio. In a possible embodiment, the machine learning model has been trained in supervised mode using artificial samples generated by sampling chromosomes of a set of real samples sharing the same HRD status and accounting for differences in purity ploidy ratio.
In one embodiment, the trained machine learning model computes a CSI score of the subject DNA sample and uses this CSI score to classify the probability of the subject DNA sample belonging to a group of samples with a certain CSI status.
In one embodiment, the trained machine learning model computes a CSI score of the subject DNA sample and uses this CSI score to determine a HRD status of the subject DNA sample belonging to a group of samples.
In one embodiment provided is a computer-based method of determining a CSI status of a DNA sample, wherein a CSI status is a HRD status.
In one embodiment provided is computer-based method of determining a CSI status of a DNA sample, wherein a DNA sample is DNA from a cancer sample.
In one embodiment provided is computer-based method of determining a CSI status of a DNA sample, wherein a low-pass whole genome sequencing coverage is at least 0.1× and up to 30×, such as between 1× to 10×, such as between 0.1× to 5× or such as between 0.1× to 1×.
In one embodiment provided is computer-based method of determining a CSI status of a DNA sample, wherein a low-pass genome sequencing is at least 0.1× and up to 30×, such as between 1× to 10×, such as between 0.1× to 5× or such as between 0.1× to 1× and obtained is a set of sequencing reads of a subject DNA sample that contain reads from at least 2 and up to 22 chromosomes.
In one embodiment provided is computer-based method of determining a CSI status of a DNA sample, wherein performed is an alignment of the sequencing reads of the DNA sample to a human reference genome.
In one embodiment provided is computer-based method of determining a CSI status of a DNA sample, wherein performed is a normalization by GC content to apply a GC-bias correction.
In one embodiment provided is computer-based method of determining a CSI status of a DNA sample that is a HRD status.
In one embodiment provided is computer-based method of determining a CSI status of a DNA sample, wherein a machine learning model, was previously trained using a training data set of samples of known genomic instability status such as HRD+/HRD− status and thus this machine learning model is trained to distinguish samples with coverage profiles characteristic of a samples with a certain CSI status such as a HRD status.
In one embodiment provided is an in vitro method of characterizing a DNA sample obtained from a patient, the method comprising:
In one embodiment provided is computer-based method of determining a CSI status of a DNA sample, wherein performed is a normalization by GC content to apply a GC-bias correction.
A method for selecting a cancer patient for treatment with a PARP inhibitor comprising a step of detecting that a tumor patient sample is HRD positive with any of the methods described herein.
In an embodiment, patients are suffering from a cancer, that can be classified with the use any of the methods described herein as either being HRD positive (HRD+) or being HRD negative (HRD−).
In one embodiment provided is an in vitro method of characterizing a DNA sample obtained from a patient, wherein the patient sample is a tumor sample and a chromosomal spatial instability indicator of the patient sample is a predictor of the tumor response to a cancer treatment regimen comprising a platinum-based chemotherapeutic agent, a DNA damaging agent, an anthracycline, a topoisomerase I inhibitor, or a PARP inhibitor.
In one embodiment provided is an in vitro method of characterizing a DNA sample obtained from a patient, comprising a step of identifying based on a chromosomal spatial instability indicator, whether a tumor patient sample is homologous recombination (HR)-deficient and wherein a high likelihood of being homologous recombination (HR)-deficient indicates that a PARP inhibitor may be for use in a method for treating cancer. A PARP inhibitor may be used alone or in combination with other treatments.
In one embodiment provided is a method for selecting a cancer patient for a treatment with a platinum-based chemotherapeutic agent, a DNA damaging agent, an anthracycline, a topoisomerase I inhibitor, or a PARP inhibitor, comprising a step of detecting that the patient DNA cancer sample is HRD positive according to the computer-based method of determining a CSI status of a DNA cancer sample as described herein. In one embodiment methods for selecting a cancer patient for a treatment are an in vitro methods.
In one embodiment provided is computer-based method of determining a CSI status of a DNA cancer sample as described herein, wherein the method further comprises: when it is determined that the DNA cancer sample is homologous recombination deficient, treating the cancer by administering a poly ADP ribose polymerase (PARP) inhibitor to the test subject.
According to one aspect, a sample is a patient sample and is in a form of tissue, fresh frozen tissue (FFT), blood or any body fluid, or cytological specimens/preparation (FFPE, smears), and the like. According to one aspect, a sample is a patient tumour sample including FFPE samples.
In another particular embodiment, patients are suffering from a cancer, in particular high grade serous ovarian cancer, prostate cancer, breast cancer, pancreatic cancer and the like.
According to one embodiment, provided is a method of diagnosing of a cancer type being HDR+ or HDR− according to the methods described herein, wherein a sample is obtained from a subject that has or is suspected of having a cancer.
In a possible embodiment is provided a method of determining a CSI status of a subject DNA sample, the method comprising:
In a possible embodiment is provided a method of determining a CSI status of a subject DNA sample, the method comprising:
As will be apparent to those skilled in the art of genomic analysis, the first and the second library may be independently prepared from the same sample at the same time, or at different times, in any order. In a possible embodiment, the first and second library may be obtained by dividing a library into two subsets, one subset (the first library—WGS library for CSI analysis) going through steps c) and d) and the other subset (the second library—targeted enrichment library for variant calling analysis) going through subsets e) to h).
In another possible embodiment is provided a method of determining a CSI status of a subject DNA sample, the method comprising:
As will be apparent to those skilled in the art of genomics, the above workflows may use different steps of genomic analysis.
The step of providing a sample material comprising nucleic acids may be performed according to any known DNA extraction method. A genomic barcoding method may be employed to identify the sample material in the downstream analysis if multiple samples are to be sequenced together in a single experiment.
The step of preparing a nucleic acid sequencing library may 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 step b.0.) of amplifying nucleic acid sequences. Amplifying nucleic acid may be performed according to any known method such as polymerase chain reaction (PCR).
The step of performing targeted enrichment may be performed according to any known method of targeted enrichment DNA-sequencing or RNA-sequencing such as target hybridization capture (i.e., hybridization capture-based targeted sequencing) or amplicon-based approaches (amplicon sequencing).
According to one embodiment, targeted enrichment is a capture-based targeted enrichment. According to this embodiment, the step of performing targeted enrichment comprises steps of hybridizing at least one probe or a probe set to target nucleic acids from genomic regions known to possibly carry variants of interest (for instance but not limited to the BRCA1 and BRCA2 regions) (“targeted regions”, or “enriched regions”), washing away non-target nucleic acids and enriching nucleic acids. It is understood that the term “probe” or interchangeably a “bait” or “(probe) nucleic acid molecule” or “capture probe” or “(DNA/RNA) oligonucleotide (capture) probe” refers to a nucleic acid molecule that can hybridize to a target nucleic acid molecule. Any known probe design can be used. The term “target nucleic acid” refers to a nucleic acid region within a gene or a transcript that may be captured by the used probe. Preferably target nucleic acids are genes related to HRR pathway such as selected from, but not limited to, BRCA1 and BRCA2.
It is understood that a nucleic acid sequencing library obtained in step b) may be divided into at least two parts, wherein one part will not be subjected to a targeted enrichment and thus will remain a non-enriched nucleic acid library, and wherein another part will be subjected to a targeted and will form an enriched nucleic acid library.
It is understood that a non-enriched nucleic acid library and an enriched nucleic acid library may be loaded together or separately to a sequencer. Therefore, nucleic acids sequences may be obtained separately or together depending on the chosen workflow.
When a non-enriched nucleic acid library and an enriched nucleic acid library are to be loaded together to a sequencer, these libraries are loaded at pre-determined concentrations. These relative concentrations are a function of the desired coverage at targeted (enriched) and untargeted regions such that a coverage of nucleic acids of a non-enriched library is from between 0.1× to 10×, between 0.5× to 10×, at least 0.2×, or preferably at least 0.1×, and at most 30×; and a coverage of nucleic acids of an enriched library is at least 30×, or between 30× and 100×, or between 100× and 500×, or between 500× and 1000×, or between 1000× and 5000×, preferably at least 4000×.
The step of sequencing nucleic acid libraries may be performed according to any known method and using known sequencers.
In a possible embodiment, the analysis of sequences from a non-enriched nucleic acid library and the analysis of sequences from an enriched nucleic acid sample are performed separately.
When the two libraries are sequenced together, the analysis of sequences from a non-enriched nucleic acid library to obtain a CSI status according to the computer-based methods as described herein may comprise a further step of filtering out the regions covered in the targeted enrichment library. These regions may be masked or filtered out by the CSI analyser according to any known masking strategy, for instance:
The step of masking the targeted enriched regions of the genome obtained from the enriched library in the combined sequencing of the non-enriched and enriched libraries ensures that coverage differences introduced by targeted enrichment do not impact the CSI analysis.
In a possible embodiment, the analysis of sequences from an enriched nucleic acid library include variant calling according to known methods. Data obtained from this analysis regarding variant allelic fraction (VAF) are additional information that may be used to further characterize a patient DNA sample. In a possible embodiment, a HRD status for the sample may be derived from the CSI status, from the variant calling results (for instance if one or more variants in genes associated with the homologous repair deficiency have been identified, for instance, but not limited to, certain variants in the BRCA1 or the BRCA2 genes), or from any combination thereof (for instance if either a CSI status has been identified as positive or if one or more variants in genes associated with the homologous repair deficiency have been identified). In a possible embodiment, the variant calling analysis may be conducted first and the CSI analysis for the sample according to the computer-based methods as described herein may be analysed only if the variant calling results did not identify a HRD status of the sample, so as to identify a possible HRD status from a genomic instability along the whole genome of the sample.
In a first experiment, the CNN of
The original BAM file of each tumor sample was first downsampled to 10 million pair-end reads to mimic low-pass WGS, the coverage signal was first pre-computed using a first set of high-resolution coverage bins of 100 kbp, and was preprocessed by performing sample and GC normalization. The resulting normalized coverage signal (y-axis) is displayed respectively on
The normalized coverage data was further collapsed into an average 3 Mbp target bin size in the 2.5 Mpb-3.5 Mbp range using the adaptive binning strategy of
The normalized coverage data for samples in training dataset was further spatially arranged to form a 2D array of 84 bins*22 chromosomes that was used to train a CNN model to predict HRD status. The coverage bins for 22 chromosomes are plotted as 22 rows from chr1 down to chr 22, aligned with respect to their centromeric bins. For shorter chromosome arms, empty bins are filled by copying the value of nearest telomeric bin present in the same row.
Low pass Whole Genome Sequencing and capture-based target panel was performed on OVKATE (RRID:CVCL_3110) cells (https://web.expasy.org/cellosaurus/CVCL_3110). DNA was extracted using the DNeasy Blood & Tissue Kits (Qiagen) according to the manufacturer's instructions. Whole-genome sequencing libraries were then prepared using SOPHiA Genetics library preparation kit. A portion of the whole-genome libraries was used to perform target enrichment using a probe panel covering genes related to the HRR pathway with the SOPHiA Genetics capture protocol. The enriched libraries and the whole-genome libraries were then sequenced together on a Nextseq Mid flow-cell (Illumina). The two types of libraries were loaded and balanced on the flow-cell such as to attain a coverage of approximately 1-2× on the genome and of more than 1000× on the target regions of the capture panel.
The DNA extracted from cultured Ovkate cells was processed following the NGS workflow combining low-pass WGS and targeted enrichment. The NGS data was processed bioinformatically to produce a BAM file containing information of reads aligned to a reference genome. The reads mapping to targeted regions allowed to perform variant calling and measure allelic fractions in enriched regions
Therefore, combination of lpWGS and capture-based target sequencing allows parallel detection of a sample CSI status and the mutational HRD status at genes of interest.
In one experiment, first selected were 169 real samples of fresh frozen tissue with known HRD status from the publicly available data of Nik-Zainal et al., “Landscape of somatic mutations in 560 breast cancer whole-genome sequences”, Nature 534, 47-54 (2016). For each artificial sample, the chromosomes from a random number of the real samples with the same HRD status may be randomly combined to create the data augmented (DA) training sample. The number of samples to create an artificial sample may be a random number drawn from an exponential distribution using: N=K * exp(−K*x) with x a random number and K=⅓. Using such an approach ensures that data from only a limited number of samples is combined to assemble an artificial sample.
Next, identified were the sample(s) in the pool with the lowest purity/ploidy ratio and decreased the purity of all the other samples so that the purity/ploidy ratio of all samples in the pool was the equal to this ratio. This enables to prevent the introduction of differences in the amplitude of coverage between chromosomes randomly selected from samples with different tumor purity and ploidy (
Such set of DA training data may be used for training a machine learning model. As will be apparent to those skilled in the art of machine learning, such a data augmentation strategy may facilitate the training of a CSI classifier machine learning model when the large volumes of data that may be needed to support the performance of some classes of methods such as HRD are not available. Data augmentation may be prepared in-silico to increase the size and diversity of the training datasets. Successful data augmentation strategies that preserve the key properties of the augmented dataset result in leading to artificial and real data that cannot be distinguished.
The proposed machine learning methods further bring additional advantages as they are inherently adaptive to the increasing availability of personalized medicine data, in particular in oncology practice. They do not require explicit biology models (for instance of suspected events occurring preferably in certain chromosome arm regions, such as around the telomeres or the 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 less sensitive to noisy data or low tumor content in the DNA sample. As more data becomes available, the model may be retrained without the need to rearchitect the runtime genomic analyzer system and workflows, as only the model parameters change. Trained models initially developed for a specific application (e.g. breast cancer) may be transferred to other applications (e.g. ovarian cancer, prostate cancer or pancreatic cancer) and/or other sample types (FFPE, FFT, cfDNA, ctDNA). They may also be used to predict different statuses, including responses to different therapies and treatments.
The methods proposed herein are suitable for use on lpWGS sequencing data and are thus cheaper and easier to implement and deploy in standard clinical practice than prior art methods that use high coverage (>30×) WGS or SNP arrays.
While various embodiments have been described above, it should be understood that they have been presented by way of example and not limitation. It will be apparent to persons skilled in the relevant art(s) that various changes in form and detail can be made therein without departing from the spirit and scope. In fact, after reading the above description, it will be apparent to one skilled in the relevant art(s) how to implement alternative embodiments.
While exemplary embodiments and applications of the proposed methods have been described in relation with the targeted next generation sequencing genomic analysis, it will be apparent to those skilled in the art of bioinformatics that they may also be adapted to apply to the detection and classification of HR deficiency out of alternate tumor sample genomic analysis workflows, for instance using genomic data out of SNP arrays and array CGH wet lab workflows. An SNP array, an array CGH, or a next-generation-sequencing (NGS) technology may be used to produce a count that varies with the copy number. Moreover, low pass WGS (0.1× to 5×) or large targeted enriched sequencing libraries as may be obtained for instance from WES (whole exome sequencing) or CES (capillary electrophoresis) panels may be solely used or combined to generate the input data to the proposed machine learning methods. Alternatively, low-pass WGS can be combined with small targeted panels (amplicon-based or capture-based) in a single or separate wetlab workflows. Possibly, off-target reads issued from targeted sequencing methods may also be used as main or complementary input data into the multidimensional array, while remaining well suited for a trained neural network as the input data processing component in the CSI Analyzer 123 architecture.
While exemplary embodiments and applications of the proposed methods have been described in relation with a coverage data signal array arranged as an image of the chromosome coverage bins arranged in rows and vertically aligned by their centromeric bins, it will be apparent to those skilled in the art of bioinformatics that various other embodiments are also possible.
In a possible embodiment, the chromosome arms may be represented as rows in the array and the centromeric bins for all arms may be aligned along the first column of the array. In another possible embodiment, the chromosome arms may be represented as rows in the array and the telomeric bins for all arms may be aligned along the first column of the array. In another possible embodiment, the chromosome arms may be represented as rows in the array and the centromeric bins for all arms may be aligned along the last column of the array. In another possible embodiment, the chromosome arms may be represented as rows in the array and the telomeric bins for all arms may be aligned along the last column of the array.
In a possible embodiment, the chromosome arms may be represented as columns in the array and the centromeric bins for all arms may be aligned along the first row of the array. In another possible embodiment, the chromosome arms may be represented as columns in the array and the telomeric bins for all arms may be aligned along the first row of the array. In another possible embodiment, the chromosome arms may be represented as columns in the array and the centromeric bins for all arms may be aligned along the last row of the array. In another possible embodiment, the chromosome arms may be represented as columns in the array and the telomeric bins for all arms may be aligned along the last row of the array.
In a possible embodiment, the whole chromosomes may be represented as columns in the array and the centromeric bins for all chromosomes may be aligned along a row of the array. The alignment row for the centromeric bins may be at the center of the array or it may be shifted up or down from the center according to the respective bin lengths for the p-arms and the q-arms of the set of chromosomes. The bins associated with the p-arm and the bins associated with the q-arm for a chromosome may be on above and below respectively from the centromeric bin row, or conversely below and above from the centromeric bin row.
In a possible embodiment, the whole chromosomes may be represented as rows in the array and the centromeric bins for all chromosomes may be aligned along a column of the array. The alignment column for the centromeric bins may be at the middle of the array or it may be shifted right or left from the middle according to the respective bin lengths for the p-arms and the q-arms of the set of chromosomes. The bins associated with the p-arm and the bins associated with the q-arm for a chromosome may be right and left respectively from the centromeric bin column, or conversely left and right from the centromeric bin column.
In a possible embodiment, the whole chromosomes may be represented as rows in the array and one of the two telomeric bins for all chromosomes may be aligned along a column of the array. In a possible embodiment, the alignment column for the telomeric bins may be the first column of the array. In another possible embodiment, the alignment column for the telomeric bins may be the last column of the array. In a possible embodiment, the p-arm telomeric bins may be used for the alignment. In another possible embodiment, the q-arm telomeric bins may be used for the alignment. In another possible embodiment, the choice of the p-arm or the q-arm telomeric bin for the alignment may be selected individually for each chromosome.
In a possible embodiment, the whole chromosomes may be represented as columns in the array and one of the two telomeric bins for all chromosomes may be aligned along a row of the array. In a possible embodiment, the alignment row for the telomeric bins may be the first row of the array. In another possible embodiment, the alignment row for the telomeric bins may be the last row of the array. In a possible embodiment, the p-arm telomeric bins may be used for the alignment. In another possible embodiment, the q-arm telomeric bins may be used for the alignment. In another possible embodiment, the choice of the p-arm or the q-arm telomeric bin for the alignment may be selected individually for each chromosome.
The Genomic Data Analyzer 120 computer system (also “system” herein) 120 may be programmed or otherwise configured to implement different genomic data analysis methods in addition to the CSI and HRD analyzer systems and methods as described herein, such as receiving and/or combining sequencing data, calling copy number alterations and/or annotating variants to further characterize the tumor samples. The machine learning model may employ extended multidimensional array inputs further comprising information such as GC content, bin size, mappability, mapping quality, variant allele fraction (VAF) in addition to the normalized coverage data 2D array input. The machine learning model may also employ additional scalar inputs providing tumor content information (purity), sample ploidy information. The Genomic Data Analyzer system may also be used to assess the quality of a sample (i.e. measure the degree of FFPE degradation). The Genomic Data Analyzer system may also be used for classification of cancers based on their CSI status level. The Genomic Data Analyzer system may also be used for stratification of cancers based on their CSI status and association with immune escape events (Bakhoum et al., 2018, Cell 174 (6), p. 1347-1360).
While exemplary embodiments and applications of the proposed methods have been described in relation with a CNN as the trained machine learning model, in a possible embodiment (not illustrated), a random forest machine learning model may be alternately employed. Random forests or random decision forests are an ensemble learning method for classification, regression and other tasks that operate by constructing a multitude of decision trees at training time and outputting the class that is the mode of the classes (classification) or mean prediction (regression) of the individual trees. Random decision forests correct for decision trees' habit of overfitting to their training set which makes it a good candidate for high dimensional input data.
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 belongs. The terminology used in the description herein is for describing particular embodiments only and is not intended to be limiting. As used in the description and the appended claims, the singular forms “a,” “an,” and “the” are intended to include the plural forms as well, unless the context clearly indicates otherwise.
Unless indicated to the contrary, the numerical parameters set forth in the following specification and attached claims are approximations that may vary depending upon the desired properties sought to be obtained and thus may be modified by the term “about”. At the very least, and not as an attempt to limit the application of the doctrine of equivalents to the scope of the claims, each numerical parameter should be construed in light of the number of significant digits and ordinary rounding approaches.
Notwithstanding that the numerical ranges and parameters setting forth the broad scope are approximations, the numerical values set forth in the specific examples are reported as precisely as possible. Any numerical value, however, inherently contains certain errors necessarily resulting from the standard deviation found in their respective testing measurements. Every numerical range given throughout this specification will include every narrower numerical range that falls within such broader numerical range, as if such narrower numerical ranges were all expressly written herein.
As will be apparent to those skilled in the art of digital data communications, the methods described herein may be indifferently applied to various data structures such as data files or data streams. The terms “data”, “data set”, “data structures”, “data fields”, “file”, or “stream” may thus be used indifferently throughout this specification.
Although the detailed description above contains many specific details, these should not be construed as limiting the scope of the embodiments but as merely providing illustrations of some of several embodiments.
In addition, it should be understood that any figures which highlight the functionality and advantages are presented for example purposes only. The disclosed methods are sufficiently flexible and configurable such that they may be utilized in ways other than that shown.
| Number | Date | Country | Kind |
|---|---|---|---|
| 20187813.9 | Jul 2020 | EP | regional |
| Number | Date | Country | |
|---|---|---|---|
| Parent | 17386255 | Jul 2021 | US |
| Child | 17534368 | US |
