Patent Description:
High throughput next-generation sequencing (NGS) or massively parallel sequencing (MPS) technologies have significantly decreased the cost of DNA sequencing in the past decade. NGS has broad application in biology and dramatically changed the way of research or diagnosis methodologies. For example, RNA expression profiling or DNA sequencing can only be conducted with a few numbers of genes with traditional methods, such as quantitative PCR or Sanger sequencing. Even with microarrays, profiling the gene expression or identifying the mutation at the whole genome level can only be implemented for organisms whose genome size is relatively small. With NGS technology, RNA profiling or whole genome sequencing has become a routine practice now in biological research. On the other hand, due to the high throughput of NGS, multiplexed methods have been developed not just to sequence more regions but also to sequence more samples. Compared to the traditional Sanger sequencing technology, NGS enables the detection of mutation for much more samples in different genes in parallel. Due to its superiorities over traditional sequencing method, NGS sequencers are now replacing Sanger in routine diagnosis. In particular, genomic variations of individuals (germline) or of cancerous tissues (somatic) can now be routinely analyzed for a number of medical applications ranging from genetic disease diagnostic to pharmacogenomics fine-tuning of medication in precision medicine practice. NGS consists in processing multiple fragmented DNA sequence reads, typically short ones (less than <NUM> nucleotide base pairs). The resulting reads can then be compared to a reference genome by means of a number of bioinformatics methods, to identify small variants such as Single Nucleotide Polymorphisms (SNP) corresponding to a single nucleotide substitution, as well as short insertions and deletions (INDEL) of nucleotides in the DNA sequence compared to its reference.

In some pathologies, a specific gene variant has been associated with the illness, such as for instance the BRCA1, BRCA2 and TP53 tumor suppressor gene variants in certain forms of hereditary breast and ovarian cancers. Beyond inherent germline mutations, cancer cells often carry somatic gross chromosomal aberrations such as copy number variations, causing duplications or hemizygous deletions of certain alleles or genomic regions. Some of these variants may thus cause a loss of function of the tumor suppressor mechanisms, such as the homologous repair (HR) function, thus making the cancer more aggressive. Identifying those variants is of key importance in recent advances of personalized cancer medication, as the resulting tumors have been shown particularly responsive to certain therapeutic compounds such as the PARP (poly ADP-ribose polymerase) inhibitors. Rottenberg et al. reported for instance the high-sensitivity of BRCA1-deficient mammary tumors to the AZD2281 PARP inhibitor in "<NPL>, while Evers et al. identified it for BRCA2-deficient mammary tumors in "<NPL>. Another example of a PARP inhibitor that exploits tumor DNA damage response (DDR) pathway deficiencies to selectively destroy cancer cells is olaparib, currently approved in the US and in Europe for certain categories of ovarian cancers associated with BRCAgene mutations. Example of prior art methods to screen a genomic loss-of-function in a subject suffering from a cancer and to administer a PARP inhibitor according to the genomic loss-of-function status include genomic expression levels in cancer tissue samples, for instance <CIT>. There is therefore a clinical oncology interest to better characterize those genomic variants from somatic tumor samples.

Recent research has shown that beyond well-studied germline variants, somatic variants such as substitutions, insertions, deletions and rearrangement patterns can cause BRCAness (BRCA1 or BRCA2 loss of function) with different mutation profiles. In "<NPL>, Whole Genome Sequencing (WGS) is shown as a possible method to detect those loss-of-function mutation profiles, with a <NUM>% sensitivity, from mutational signatures in the whole genome sequencing data. In contrast to WGS, the authors report that when only Whole Exome Sequencing is applied, the sensitivity significantly drops, between <NUM>% and <NUM>% depending on specific adaptations of the bioinformatics algorithms.

In "<NPL>et al. describe a probe-based target enrichment assay coupled with Next Generation Sequencing (NGS) for deletion of CNVs and SNPs, including hemizygous status. This method assumes independent germline status testing from blood samples in addition to the somatic sample analysis. To provide enough sensitivity and specificity, it assumes a coverage of at least 500x by the dedicated probe-based assay, as well as sufficient DNA input from the FFPE samples. These requirements restrict its cost-effective routine use.

In "<NPL>, Dougherty et al. present preliminary clinical evidence that olaparib may benefit to ovarian cancer patients with somatic BRCA1/<NUM> mutations like it benefits to those with germline mutations. Moreover, they observe that somatic mutations had more than <NUM>% biallelic inactivation frequency and were predominantly clonal, suggesting that BRCA1/<NUM> loss occurs early in the development of these cancers. As explained by the authors, the conventional germline BRCA1 and BRCA2 mutational status by PCR and Sanger sequencing does not perform cost effectively in somatic applications. In particular it cannot accurately detect the low allele frequencies encountered in small amounts of low-quality DNA sequenced from FFPE tissue with variable tumour cellularity.

"<NPL>), describes methods to identify somatic copy number aberrations from cancer sequencing by directly accounting for tumor purity and clonality. These methods use a tumor-normal pair, i.e., data obtained from tumor (somatic) and normal (germline) samples from the same individual, and whole genome sequencing data preferably at a decent coverage of 30x or higher, which requires a significant genomic data processing, communication and storage cost when implemented in an automated genomic analysis workflow.

"<NPL> et al. reviews the challenges and possible solutions of NGS sequencing of FFPE samples for clinical oncology, using the TrueSeq Amplicon Cancer Panel, a targeted sequencing technology from Illumina, with various bioinformatics pipelines. They report the detection of mutations at allele frequencies down to <NUM>%, yet assuming an overall mean coverage of 2084x and median coverage of 2016x. The authors highlight the need for solutions to analyse the co-occurrence and relevance of somatic mutations, CNVs and germline variants and suggest the use of supervised machine learning models using paired tumour-normal samples obtained from the same patient to improve the prior art empirical methods of variant allele frequency cut-offs or nucleotide substitutions types. However, in routine clinical practice, these approaches are hindered by the extra cost of gathering tumor-normal samples from the patients.

There is therefore a need for dedicated NGS assays and methods development, which can be easily and cost-effectively deployed in automated NGS setups for routine diagnosis solely from tumour samples, without requiring manual data analysis and/or complementary biological assays such as germline (blood sample) genomic analysis.

Rather than sequencing the whole genome (WGS) from an individual sample, the genomic analysis can focus on the genome region associated with the illness, by targeting, with a set of region-specific DNA primers or probes, and enriching or amplifying, for instance with PCR (Polymerase Chain Reaction), the biological DNA sample specifically for sub-regions corresponding to the gene along the DNA strand. A number of next generation sequencing assays have now been developed along those principles as ready-to-use biological kits, such as for instance the Multiplicom MASTR™ (Dx (http://www. multiplicom. com/product/brca-mastr-dx) assay kits to facilitate DNA based diagnostics with next generation sequencers, such as for instance the Illumina MiSeq® sequencer, in medical research and clinical practice.

Target enrichment may be achieved from a small sample of DNA by means of probe-based hybridization (on arrays or in-solution) or highly multiplexed PCR-based targeted exon enrichment, so that both the gene coverage/read depth and the amplification specificity (amplifying the right region, as measured by further alignment to the desired target regions) are maximized. Examples of commercially available target enrichment systems include Agilent SureSelect™ Target Enrichment System, Roche NimbleGen SeqCap EZ, Illumina Nextera Rapid Capture, Agilent Haloplex™, and Multiplicom MASTR™.

In order to maximize the use of the massively-parallel processing NGS sequencer, a number of samples are multiplexed in the targeted NGS experiment - a pool of <NUM> or more target enrichment samples can thus be simultaneously input to the Illumina MiSeq sequencer for instance. Raw sequencing data out of the NGS sequencer may then be analyzed to identify specific subsequences, for instance by alignment to a reference genome. As a result, the amplification may produce more than a thousand reads for a given amplicon in a patient sample.

Next Generation Sequencing (NGS) enables in particular to detect and report small changes in the DNA sequence, such as single nucleotide polymorphisms (SNPs), insertions or deletions (INDELs), as compared to the reference genome, through bioinformatics methods such as sequencing read alignment, variant calling, and variant annotation. NGS workflows refer to the configuration and combination of such methods into an end-to-end genomic analysis application. In genomic research practice, NGS workflows are often manually setup and optimized using for instance dedicated scripts on a UNIX operating system, dedicated platforms including a graphical pipeline representation such as the Galaxy project, and/or a combination thereof. As clinical practice develops, NGS workflows may no longer be experimentally setup on a case-per-case basis, but rather integrated in SaaS (Software as a Service), PaaS (Platform as a Service) or IaaS (Infrastructure as a Service) offerings by third party providers. In that context, further automation of the NGS workflows is key to facilitate the routine integration of those services into the clinical practice.

While next generation sequencing methods have been shown more efficient than traditional Sanger sequencing in the detection of SNPs and INDELs, their specificity (rate of true positive detection for a given genomic variant) and sensitivity (rate of true negative exclusion for a given genomic variant) may still be further improved in clinical practice. The specificity and sensitivity of NGS genomic analysis may be affected by a number of factors:.

Due the above biases, the NGS data is particularly noisy and often requires complementary experiments on a case per case basis to detect variants of clinical relevance at sufficient sensitivity and specificity. As example in breast and ovarian cancer genomic analysis practice, most BRCA test service providers for germline mutations require complementary analysis such as Sanger sequencing or MLPA analysis of CNV mutations in parallel with the NGS genomic analysis of SNP and INDEL variants. This approach is costly and does not scale well as more and more tests have to be conducted in daily operation by a single multi-purpose genomic analysis platform. Patent application <CIT> addresses this limitation by describing a method to detect copy-number values (CNV) from a pool of DNA samples enriched with a target enrichment technology, based on a between-samples normalization of the genomic coverage signal, which can be applied by an NGS genomic data analyzer module in parallel with SNP/INDEL detection from the same target enrichment assay and next-generation sequencing experiment instead of requiring separate NGS and MLPA experiments and analyzers. This method however assumes sufficient homogeneity in the input DNA and resulting genomic alignment coverage data for its statistical model hypotheses to be valid. That is usually the case in germline cell analysis with good laboratory practice, but it cannot always reliably infer the somatic whole gene CNVs that may be present in tumor cell DNA from the coverage signal alone, especially when only low resolution coverage depth is available.

Still, recent works in somatic NGS analysis have highlighted the challenge of detecting large rearrangements as genomic alterations of pathologic relevance out from the tumor sample itself. For instance, <NPL> indicates that while MLPA may in theory be used to detect large rearrangements in FFPE tissue, in practice it raises data analysis challenges to cope with genomic instability in the tumor genome that may affect the control probes as well as to detect low level signals associated with somatic mutations.

The Tumor BRACAnalysis CDx commercial test from Myriad GmbH (Technical Specifications - Dec <NUM>, <NUM>) analyses copy number abnormalities indicative of deletion or duplication rearrangements based on the coverage signal (number of next generation sequencing reads that map to each nucleotide, normalized to the run median depth of coverage of the same nucleotide) and is therefore subject to certain limitations, in particular it cannot reliably detect certain specific insertions, inversions and regulatory mutations, the detection of large rearrangement deletions and duplications is subject to the quality of the FFPE sample, and not all of them can be assessed. Some duplications are reported as variants of uncertain significance, which have to be further analyzed manually.

Moreover, the current state of the art in automated CNV variant calling assumes a high-resolution depth of coverage analysis, as quoted for instance in <NPL>.

There is currently no known efficient genomic analysis method to automatically and reliably detect and categorize the biallelic loss of function status of tumor suppressor genes by analyzing FFPE tumor samples solely from next generation sequencing variant information without complementary biological or manual analyses.

There is therefore a need for improved genomic analysis methods to efficiently, automatically and reliably detect and categorize genomic alterations from FFPE tumor DNA samples to better determine the biallelic loss of function status of target genes in each individual sample, regardless of whether this loss of function is caused by germline or somatic mutations or chromosomal aberrations or a combination thereof.

It is proposed a method for identifying the presence of a chromosomal aberration in the tumor sample cells in at least one of several genes from the next generation sequencing data analysis of a pool of patient tumor samples comprising at least two analyzed genes, with a variant classification data processor module which operates solely on variant calling information from samples comprising a mix of germline and somatic cells and for each possible sample purity value p, in which a combination of at least two genomic alterations have caused its biallelic loss of function, the said data processor module performs a method comprising:.

The variant fraction of each putative germline heterozygous SNP may be compensated for capture-biases, by inferring and compensating the average loss a in capture efficiency due to the presence of a SNP in the tumor samples as the value a for which the compensated variant fraction distribution of all putative germline heterozygous SNP across all samples is maximally symmetric around <NUM>%. Possible chromosomal aberration in the tumor sample cells may comprise copy-neutral loss-of-heterozygosity, the deletion or the duplication of one allele.

The mixture model may be fitted to the data by minimizing a cost function based on the root mean squared error (RMSE) between the observed and the theoretical variant fraction expected for the most likely chromosomal aberration. Inferring chromosomal aberrations may then comprise minimizing the root mean squared error (RMSE) between the observed variant fraction and the theoretical variant fractions expected for a plurality of chromosomal aberrations.

A mixture model with mixture components corresponding to the different possible chromosomal aberrations using the Expectation-Maximization algorithm or using a Markov-Chain Monte-Carlo sampling algorithm may also be fitted to the data to estimate the sample purity. Inferring the gene chromosomal aberrations may then comprise applying Mixture Model Classifiers.

The method may further comprise detecting germline copy number variants (CNVs) in accordance with a Hidden Markov Model.

The method may further comprise reporting to the end user whether the detected alterations are biallelic, for each gene in each sequenced patient tumor sample.

A least one of the genes to be analyzed may be a tumor suppressor gene, such as for instance one of ATM, CHEK1, CHEK2, BARD1, BRCA1, BRCA2, BRIP1, RAD51C, RAD51D, FAM175A, MRE11A, NBN, PALB2, TP53, APC, DCC, DF1, NF2, PTEN, Rb, VHL, WT1, PP2A, LKB1, or INK4a/ARF.

The proposed methods to detect genomic alterations may also be used for determining if cancer or precancerous lesions or benign tumours in a patient will be responsive to treatment with a PARP inhibitor, comprising:.

<FIG> shows an exemplary genomic analysis system comprising a DNA enrichment assay <NUM>, a next generation sequencer <NUM> and a genomic data analyzer <NUM>.

In a NGS laboratory, a pool of DNA samples is processed by the DNA enrichment assay <NUM> to generate a library of pooled amplicons (for amplicon-based enrichment) or fragments (for probe-based enrichment) as DNA fragments input to the next generation sequencer <NUM>, each set of amplicons/fragments corresponding to a different sample. The number of amplicons/fragments is application dependent. In some genomic analysis experiments, target enrichment may require <NUM> primers to enrich <NUM> different regions to be targeted out of the sample genome, resulting in a set of <NUM> amplicons for each sample. The number of samples may also be adapted to the next-generation sequencing sequencer <NUM> parallel processing capability, for instance <NUM> samples in the form of a library of pooled amplicons may be sequenced in parallel by an Illumina MiSeq sequencer. Other NGS sequencer technologies may be used, such as for instance the Roche <NUM>™ GS Junior or GS FLX, Illumina MiSeq®, or Life Technologies Ion PGM™ sequencers.

The next-generation sequencer <NUM> analyses the input samples and generates sequence reads in a computer-readable file format representing raw NGS sequencing data. Depending on the NGS technology, one or more files may be output by the NGS sequencer <NUM>. For instance with Illumina sequencers, the FASTQ file format may be used with two different files for forward and reverse reads or as a single joined file. This text file typically starts with a sequence header marked by a '@' start character, followed by one line of sequence information represented as a string of 'A', 'T', 'C', 'G' nucleotide characters, then by a quality header marked by a '+' start character, followed by one line of quality metrics, one quality score matching each nucleotide read. The format for the quality metrics for each nucleotide in the sequence information string may depend on the sequencer. Some legacy sequencers output the raw sequencing data in the SFF (Standard Flowgram Format) binary file format, which comprises an informative header and the read data. Other outputs are also possible, for instance some legacy Roche sequencers output multiple FASTQ files for a single patient analysis, while other sequencers, for instance the Ion Torrent PGM sequencers, have migrated to the compressed unmapped BAM file format, as may be recognized from the. basecaller. bam file extension. As known to those skilled in the art of communication systems, the laboratory operates a computing infrastructure to store the resulting raw NGS sequencing data file in a laboratory biobank. The laboratory computing infrastructure connects, with authentication credentials, through a communication network, to the genomic data analyzer <NUM> and transmits a genomic analysis request comprising the raw NGS sequencing file to the genomic data analyzer <NUM>.

The genomic data analyzer <NUM> computer system (also "system" herein) <NUM> is programmed or otherwise configured to implement different genomic data analysis methods, such as receiving and/or combining sequencing data and/or annotating sequencing data.

The genomic data analyzer <NUM> 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, and the like. 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. The genomic data analyzer <NUM> may be integrated into a massively parallel system. IThe genomic data analyzer <NUM> may be directly integrated into a next generation sequencing system.

The genomic data analyzer <NUM> 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.

As illustrated on <FIG>, the genomic data analyzer <NUM> may comprise a sequence alignment module <NUM>, which compares the raw NGS sequencing data to a reference genome. The sequence alignment module <NUM> 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 algorithms 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 <NUM> requirements for storage optimization and/or genomic data privacy enforcement.

The resulting alignment data may be further filtered and analyzed by a variant calling module <NUM> to retrieve variant information such as SNP and INDEL polymorphisms. The variant calling module <NUM> may be configured to execute different variant calling algorithms. The resulting detected variant information may then be output by the genomic data analyzer module <NUM> as a genomic variant report for further processing by the end user, for instance with a visualization tool, and/or by a further variant annotation processing module (not represented).

The genomic data analyzer <NUM> may be adapted to automatically detect, with a processor, a set of characteristics that uniquely determine the input sequencing data and corresponding genetic context, the DNA enrichment context such as the sample type or laboratory process characteristics, the DNA enrichment technology such as the targeted enrichment target kit or capture probe assay characteristics, and/or the NGS sequencing technology. As will be apparent to those skilled in the art of next generation sequencing, these experimental characteristics may cause specific biases in the sequence alignment and/or the variant calling results.

The proposed genomic data analyzer system <NUM> may thus serve next generation sequencing genomic analysis requests from different labs that are independently operating different sequencer technologies and different DNA enrichment technologies on different samples for different genes. The proposed genomic data analyzer system <NUM> may automatically detect a set of characteristics from the input data and requests received from the laboratory and may adapt the configuration of the sequence alignment module <NUM> and the variant calling module <NUM> accordingly, without requiring a time consuming and costly manual setup to minimize the data biases possibly induced by each different biological workflow. As will be apparent to those skilled in the art, there may be dozens or even hundreds of different clinical laboratory setups for multiple sourcing laboratories operating with the same genomic analyzer <NUM>, and this number and diversity is likely to further increase with the deployment of additional technologies and assays as the NGS-based personalized medicine clinical practice develops.

Depending on the detected genomic experiment characteristics, the genomic data analyzer <NUM> may configure the sequence alignment module <NUM> to operate additional data processing steps and/or use different sets of configuration parameters such that the data biases caused by the genomic experiment characteristics are minimized.

Depending on the detected input characteristics, the genomic data analyzer may further configure the variant calling module <NUM> to operate additional data processing steps and/or use different sets of configuration parameters such that the data biases caused by the genomic experiment characteristics are minimized.

Depending on the results of the initial sequence alignment by the sequence alignment module <NUM>, the genomic data analyzer <NUM> may be further adapted to identify next generation sequencing data alignment biases that become apparent when analyzing the alignment data. The genomic data analyzer may accordingly configure the sequence alignment module <NUM> to execute an additional step of re-alignment of the raw NGS sequencing data. This re-alignment may be constrained by additional parameters as may be determined from the initial alignment results. In a possible embodiment, the re-alignment is applied specifically on a sub-region of the genomic sequence. The resulting re-alignment data may be further filtered and analyzed by the variant calling module <NUM> to output a more relevant genomic variant report with increased sensitivity and specificity for variant detection.

Depending on the results of the variant calling by the variant calling module <NUM>, the genomic data analyzer <NUM> may be further adapted to identify variant calling biases that become apparent when calling variants on the alignment data. The genomic data analyzer may accordingly configure the variant calling module <NUM> to execute an additional step of re-calling variants on all or part of the alignment data. This refined variant calling step may be constrained by additional parameters as may be determined from the former alignment and/or re-alignment and/or variant calling results.

In a possible embodiment, variants are specifically called on a subset of the aligned genomic data. The resulting refined variant calling data may be further combined with the standard variant calling results by the variant calling module <NUM> to output a more relevant genomic variant calling information report with increased sensitivity and specificity for variant detection. In a possible embodiment, some variant calling results may be excluded from the genomic variant report as identified possibly biased by the variant calling module <NUM>, so that a more relevant genomic variant report is generated by the genomic data analyzer <NUM> with increased sensitivity and specificity for variant detection.

As will be apparent to those skilled in the art of bioinformatics, the variant calling information reported by the variant calling module <NUM> may be a data file in the VCF format, comprising a list of genomic variants and the associated read depth coverage information. The variant fraction of each genomic variant in the next generation sequencing data can be derived with standard bioinformatics tools such as vcftools. For instance the list of called variants, the variant fraction information and/or the coverage information may be acquired by the genomic data analyzer <NUM> from different file formats.

As will be apparent to those skilled in the art of personalized medicine, the variant calling information is quite complex to interpret and currently requires a good understanding of the underlying methods and format specifics. In order to facilitate the interpretation of the data by non-bioinformatician experts such as biologists and doctors, a variant classification module <NUM> may further analyze, annotate, categorize and/or prioritize the variant calling information to report to the end user for each tumor suppressor gene a combination of SNPs, INDELs and/or CNV variants which potentially results in a pathogenic loss of function for the tumor suppressor gene. This report may then facilitate the oncogenomics diagnosis of the cancer origin. This report may also facilitate the selection of a medical treatment specifically targeting the analyzed tumor, such as for instance the use of PARP inhibitors in certain BRCA1 and/or BRCA2 loss of function mutations.

In the case of a somatic FFPE sample analysis, the variant classification module <NUM> may report germline and somatic alterations, for instance SNPs, INDELs, as well as chromosomal aberrations such as allelic imbalance (AI), loss of heterozygosity (LOH), as copy number variations such as aberrations (CNA), gain (CVG) or decrease (CND) as may be caused by large deletions or duplications of an allele genomic region.

In the case of a somatic FFPE sample, the variant classification module <NUM> may identify that the genome of germline cells is heterozygous for certain variants, and may report separately the variant fraction of different alteration events in cancer cells such as a copy-neutral loss of function, a large deletion resulting in a lower copy number for this variant, or a large duplication resulting in a larger copy number for this variant. The variant classification module <NUM> may thus identify, categorize and prioritize, in terms of loss-of-function, the combination of different genomic alterations such as germline and somatic alteration variants and chromosomal aberrations.

As will be apparent to those skilled in the art, this bioinformatics method then significantly facilitates detailed understanding of the cancer cells genomic alterations and the adaptation of a personalized medicine treatment to the specifics of the inferred cancer cell biology for the patient.

<FIG> shows accordingly a possible genomic analysis workflow for the genomic data analyzer <NUM>, comprising:.

An exemplary refined variant calling analysis method <NUM> will now be described in more detail with the example of genomic analysis of ovarian cancer samples, but as will be apparent to those skilled in the art of bioinformatics, the proposed embodiments are not specific to ovarian cancer samples and may similarly apply to other types of cancer tumors genomic analysis.

The fully automated genomic data analysis workflow of <FIG> operates on genomic data sourced from at least one next generation sequencing laboratory. The genomic data analyzer <NUM> receives <NUM> an NGS sequencing analysis request from the laboratory, comprising raw sequencing data. The genomic data analyzer <NUM> configures <NUM> the data alignment module <NUM> to execute <NUM> a first raw data alignment. As will be apparent to those skilled in the art of bioinformatics, the data alignment module <NUM> may first execute <NUM> pre-processing steps such as removing the assay specific adapters from the reads and/or merging pair-ended raw sequencing data files into a single merged file. The data alignment module <NUM> aligns <NUM> the raw sequencing data short reads to a reference genomic sequence such as the reference human genome (hg19), with a raw data alignment algorithm as known to those skilled in the art of bioinformatics, to produce a data alignment file. Standard 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 resulting data alignment file may be represented as one or several files in BAM or SAM format, but other steps are also possible, in particular the data alignment module <NUM> may also execute <NUM> post-processing steps such as compressing and/or encrypting the alignment data, for instance with an order-preserving encryption scheme, depending on the genomic data analyzer <NUM> requirements for storage optimization and/or genomic data privacy enforcement along the genomic analysis workflow processing.

The genomic data analyzer <NUM> may also further detect <NUM> the presence of alignment mismatches especially at the beginning and/or the end of the reads ("soft clipping"), as may be due to primer mispriming. As known to those skilled in the art of bioinformatics, soft clipping information may then be re-aligned in the genomic analysis workflow with specific algorithms in order to further detect structural variants of potential clinical relevance. The genomic data analyzer <NUM> may thus automatically identify <NUM> the short reads with soft clipping regions, from the results of the data alignment <NUM>, and configure <NUM> the data alignment module <NUM> to operate a further data re-alignment <NUM> on those reads specifically by taking into account the primer anchors information corresponding to the specific DNA enrichment technology in the alignment algorithm. As will be apparent to those skilled in the art of bioinformatics, a more robust algorithm than Bowtie2 or BWA may be used specifically on those regions, even if less computationally efficient, as only a subset of the whole sequencing data needs to be re-aligned this way.

Depending on the genomic context identifier, the genomic data analyzer <NUM> may also identify from the alignment data the presence of some regions that are particularly challenging to align, such as homopolymer regions, heteropolymer regions, or other regions with specific repeat patterns. Proper alignment of corresponding next generation sequencing reads is particularly challenging as those multiple repeats cause alignment ambiguities. The genomic data analyzer <NUM> may thus automatically identify <NUM> from the results of the raw data alignment <NUM> a specific genomic context requiring refinement on the reads overlapping those ambiguous regions. The genomic data analyzer <NUM> may accordingly configure <NUM> the data alignment module <NUM> to operate a further data re-alignment <NUM> on those reads to identify other possible alignment solutions, such as for instance by taking into account the PCR error rate and comparing reads to each other.

The genomic data analyzer <NUM> may then use the target enrichment technology identifier to configure <NUM> the variant calling module <NUM> to execute different variant calling algorithms in accordance with the genomic context (e.g. BRCA1, BRCA2 and TP53 genomic analysis). The variant calling module <NUM> calls <NUM> variants on the refined alignment data to produce a first VCF file. In particular, SNPs and short INDELs may be called and annotated by the variant calling module <NUM> configured to execute a variant calling method to produce a first variant calling information data, comprising a list of variants and their coverage information.

Moreover, in addition to SNPs and short INDELs variant calling or in parallel with SNPs and short INDELs variant calling, intragenic germline CNVs variants (germline CNVs that only affect part of a gene) may be called by the variant calling module <NUM> to produce a second variant calling data, comprising a list of intragenic germline CNVs of potential clinical relevance. SNPs/INDELs and CNVs variant calling may be executed in parallel or serially by the variant calling module <NUM>. In a possible embodiment, the method of patent application <CIT> may be used, but other CNV detection algorithms may also be used.

In somatic genomic analysis, the resulting variants may not be accurate enough to be reported by the genomic data analyzer <NUM>. The genomic data analyzer <NUM> may thus identify <NUM> a fourth set of characteristics such as the need for variant phasing from the first and or second variant calling results, and configure <NUM> the variant calling module <NUM> to refine <NUM> the variant calling data to be analyzed <NUM> and reported <NUM> by the genomic data analyzer <NUM>.

The variant classification module <NUM> may annotate and prioritize <NUM> the variant calling data by using information from cancer variant databases such as COSMIC or Clinvar, and/or pathogenicity prediction tools such as Align GVGD, Sift, Polyphen2, Alamut and others, as will be apparent to those skilled in the art of bioinformatics applied to personalized medicine research. The variant classification module <NUM> may categorize <NUM> the annotated variant calling data by detecting certain genomic alterations of particular interest, such as for instance biallelic loss of function due to a combination of genomic alterations and chromosomal aberrations in certain tumour suppressor genes.

Finally, the genomic data analyzer <NUM> may report <NUM> a set of refined classified genomic variants to an end user, as a result file and/or through a graphical user interface. The genomic data analyzer may report <NUM> the classified genomic variant information to a further data processor, through a communication network, for instance to an automated tool as a support to medical diagnosis or personalized medicine treatment recommendations. IThe genomic data analyzer may record <NUM> the classified genomic variant information into a medical knowledge database, for instance for reference in genomic research or in future clinical trials scientific data.

Other genomic analysis workflows for the genomic data analyzer <NUM> are also possible. In particular, as will be apparent to those skilled in the art, the proposed variant calling refinement <NUM> and variant analysis <NUM> methods may also be configured as a data post-processing stage to refine, analyze, annotate, categorize and/or prioritize genomic variants initially identified from targeted next generation sequencing variant calling <NUM> performed by a third party genomic analysis workflow. The genomic data analyzer <NUM> may then acquire the variant calling information from a third party genomic analyzer workflow. In a possible workflow (not represented), the genomic data analyzer <NUM> acquires the variant calling information from a third party genomic analyzer workflow as VCF file through a communication network, but other embodiments and file formats are also possible. In a further possible workflow the genomic data analyzer <NUM> may also acquire annotated variant information from a third party genomic analyzer workflow. Other workflows are also possible.

In a preferred embodiment which will now be described in further detail, the variant classification module <NUM> may be configured to analyze the variant fraction information as may be extracted from the variant calling data, in order to infer unknown variables such as the sample purity (ratio of somatic DNA material in the overall FFPE sample), and to detect the presence of germline heterozygous allelic variants and the presence of different somatic allelic variant fractions due to different chromosomal aberration events, such as those inducing genomic loss-of-heterozygosity (LOH) in the various tumor cells. Indeed, we observed that LOH does not necessarily require a change in copy number. In particular, copy-neutral LOH can occur in situations where the allele harboring a heterozygous mutation is duplicated and the non-mutated allele is deleted. This is a major limitation to state of the art bioinformatics approaches which primarily infer LOH from the coverage signal.

In a possible embodiment, the variant classification module <NUM> first extracts from the variant calling information data SNPs and their variant fraction, that is the ratio of the number of NGS reads harboring each SNP variant over the total number of reads. The variant classification module <NUM> further selects as putative germline heterozygous SNPs, for which one may expect a statistical distribution of the variant fraction centered at <NUM>%, across all observed samples as the SNPs which:.

In a possible embodiment, the reference database may be dbSNP to identify variants that are commonly found in human genetics, but other embodiments are also possible, for instance Exac or <NUM>-genomes may be used as alternative reference databases for identifying SNPs commonly found in human genetic variations. In another possible embodiment, clinical variant databases such as Clinvar may be used to exclude SNP variants that are not commonly found in human genetics.

In a possible embodiment, the heterozygous SNP selection range for the variant fraction values may be defined as [<NUM>%+Δh<NUM>, <NUM>%-Δh<NUM>], where the values of Δh<NUM> and Δh<NUM> may be chosen such that Δh<NUM> ≤ <NUM>% and Δh<NUM> ≤ <NUM>%. In a possible embodiment, the heterozygous SNP selection range may be between <NUM>% and <NUM>%, but other embodiments are also possible, for instance between <NUM>% and <NUM>% or between <NUM>% and <NUM>%.

In practice, we observed that the statistical distribution of the putative germline heterozygous SNPs over all samples may be biased towards a different value than the expected <NUM>% value, such as <NUM>% in some of our experiments. This may be caused by a bias due to the impact of SNPs on capture efficiency when using standard variant fraction estimator algorithms; indeed, we observed that compared to non-mutated fragments, DNA fragments harboring a SNP may have been captured with a lower efficiency, possibly down to <NUM>%.

As will be apparent to those skilled in the art of bioinformatics, the variant classification module <NUM> may accordingly compensate for capture-bias by estimating the average drop in capture efficiency due to the presence of a SNP in a DNA fragment. In a possible embodiment, the variant classification module <NUM> may apply an unbiased estimator of the variant fraction, such that: <MAT> where VF is the variant fraction value, NALT the number of fragments with SNP in the original DNA sample material, NREF the number of fragments without SNP in the original DNA sample material; the number of fragments without SNP in the captured material may be estimated as: <MAT> and the number of fragments with SNP in the captured material may be estimated as: <MAT> where the loss in capture efficiency due to the presence of a SNP may be defined as: <MAT>.

If the variant fraction estimation was unbiased, then we would expect the variant fraction distribution to have a symmetric mode centered at <NUM>%. The loss a in capture efficiency due to the presence of a SNP may then be experimentally inferred as the value for which the variant fraction distribution across all samples in the pool is maximally symmetric around <NUM>%.

In a possible embodiment, the variant classification module <NUM> may estimate the loss a in capture efficiency to unbias the variant fraction for each putative germline heterozygous SNP i from:.

by first calculating the unbiased variant fraction xi of a SNP i, according to the following formula: <MAT> where a is the ratio between the capture efficiency of reads including a SNP and the capture efficiency of reads matching the reference. The parameter a is unknown and has to be inferred from data. As will be apparent to those skilled in the art of statistics, various methods may be employed to this end. In a possible embodiment, for any possible value of a (e.g., <NUM> < a < <NUM>), the variant classification module <NUM> may compute the unbiased variant fractions {x<NUM>, x<NUM>,. , xN} of putative heterozygous SNPs and may then estimate distribution of variant fractions p(x) by using, e.g., a kernel density estimation technique: <MAT> where K(x) is an even function. In a possible embodiment, a Gaussian function with standard deviation <NUM>% may be used, but other embodiments are also possible. The degree of symmetry around <NUM>% of p(x) may then be calculated by various algorithms. For instance, in a possible embodiment the following RMSE-based score may be calculated by the variant classification module <NUM>: <MAT> where C is a parameter of choice that can be fixed e.g. to <NUM>% in an automated variant classification module <NUM>, as will be apparent to those skilled in the art of bioinformatics. Note that since xi (used to compute p(x)) depends on a, the RMSE-based score also depends on a. The optimal parameter a may then be obtained by finding the value under which E(a) is minimized. The variant classification module <NUM> may solve this problem by using standard optimization techniques, such as a line-search algorithm. The above embodiment has been presented by way of example and not limitation. Other embodiments are also possible; as will be apparent to those skilled in the art of statistics, the variant classification module <NUM> may be adapted to implement various methods and/or parametrizations to unbias the variant fraction distribution.

<FIG> illustrates the uncorrected and the corrected variant fraction distributions in an experiment where the value of a is measured at <NUM>, corresponding to a loss in capture efficiency of <NUM>%. In our experiments, we also observed that the resulting unbiased variant fraction distribution of putative germline heterozygous SNPs across all samples is characterized by both a high variability and a high degree of symmetry around the <NUM>% axis.

It is therefore of interest to map the experimental data against a statistical model jointly taking into account the above two hypotheses. With regards to the second hypothesis, <FIG> illustrates that different LOH-inducing chromosomal aberrations may result into specific variant fraction distributions peaks:.

With regards to the first hypothesis, <FIG> illustrates how the variant fraction distributions of a germline heterozygous SNP respectively affected by a copy-neutral LOH (black - CN-LOH), a deletion-induced LOH (dark grey - LOH-DEL) and a duplication-induced LOH (clear grey - LOH-DUP) may vary with sample purity. In particular, as pointed out for example purities ofp=<NUM>% or p=<NUM>% respectively, the above listed different LOH-inducing mechanisms result in clearly different variant fraction distribution symmetrical values relative to the <NUM>% variant fraction axis. In this situation, we observed that the observed variant fraction distributions can be used to reliably jointly infer sample purity and detect separate as well as combined genomic alteration events, as will now be described in further detail.

Thus, for each possible purity ratio of a patient sample, the variant classification module <NUM> may estimate the value of the theoretically symmetrical variant fraction distribution of somatic chromosomal aberration events. For instance, "NO-LOH", "CN-LOH", "LOH-DUP" and "LOH-DEL" events as described above all result in a different theoretical variant fraction value, but other events could be assessed similarly in more complex applications (e.g. frequently found subclonal alteration events in other cancers).

The variant classification module <NUM> may thus compare the observed variant fractions of the putative germline heterozygous SNPs corresponding to this with the theoretically symmetrical variant fraction values for different purity ratios. <FIG>) respectively illustrate how two exemplary patient sets of putative germline heterozygous SNPs map to those values, as represented by dashed lines for a purity ratio ofp=<NUM>%:.

Thus, for each chromosomal aberration inducing a theoretically symmetrical variant fraction distribution and for each possible purity ratio, the variant classification module <NUM> may estimate the theoretical variant fraction value for the putative germline heterozygous SNPs corresponding to this scenario for this purity ratio, and evaluate how closely the unbiased observed variant fraction distribution fits to the theoretical variant fraction distribution value corresponding to this scenario. As will be apparent to those skilled in the art of bioinformatics, this fit may be evaluated by different statistical methods. In a simple embodiment, the variant classification module <NUM> may calculate a distance from the unbiased observed variant fraction distribution to the theoretical variant fraction distribution value corresponding to this scenario. The variant classification module <NUM> may then calculate an overall distance and select the purity ratio for which the computed statistical distance is minimal over all scenarios.

As will be apparent to those skilled in the art of statistical analysis, various statistical algorithms may be used to this end. In a preferred embodiment, the variant classification module <NUM> may fit a mixture model to the observed variant fraction distribution of the overall putative germline heterozygous SNPs in the variant calling data for each possible value p of the sample purity, the modeled variant fraction distribution values being calculated as a function of the sample purity p. The mixture model may then be fitted to the data by minimizing a cost function measured between the observed and the theoretical variant fraction expected for the most likely chromosomal aberration scenario.

For instance, a simple embodiment will now be described in further detail:.

For each possible purity p between <NUM>% and <NUM>% the variant classification module <NUM> may compute a cost function model-prediction error E(p) as the root mean squared error (RMSE) of the SNP variant fraction distribution computed with respect to its closest theoretical variant fraction values according to the following equation: <MAT> where xi is the unbiased variant fraction of the putative heterozygous SNPs i in the set of N putative heterozygous SNPs identified by the variant classification module <NUM>, and D(p) is the set of theoretically symmetric variant fractions predicted by the variant classification module <NUM> for each possible scenario under the assumption of purity = p. Indeed, the position of these theoretical variant fractions depends on the purity p, as illustrated for instance in the <FIG> example, so the model-prediction error also depends on the purity variablep.

As will be apparent to those skilled in the art of statistics, the variant classification module <NUM> may then apply an optimization algorithm to infer from the experimental data the optimal purity poptimal for which E(p) is minimized. In a possible embodiment the variant classification module <NUM> may be adapted to apply a line-search algorithm to determine the optimal purity poptimal: <MAT>.

We note that other embodiments are also possible. For instance, the variant classification module <NUM> may fit a mixture model with theoretical value components corresponding to the different scenarios, using either the Expectation-Maximization (EM) algorithm or a Markov-Chain Monte-Carlo (sampling) approach, as will be apparent to those skilled in the art of bioinformatics. In a preferred embodiment, the mixture model may be a Gaussian Mixture Model (GMM) but other embodiments are also possible, depending on the specifics of the genomic variant distributions.

Once the purity has been inferred over the putative heterozygous SNPs throughout all genes in the sample, the variant classification module <NUM> may refine the analysis by detecting the most likely chromosomal aberrations for each gene separately, in order to facilitate the further identification and classification of loss-of-function genomic alterations of interest to report to the end user. In the exemplary application of biallelic loss-of-function genomic alterations detection, the possible chromosomal aberrations of interest at this stage may match the following hypotheses (<FIG>):.

In a possible embodiment, for each gene to analyze, for instance respectively BRCA1, BRCA2 and TP53 in the case of ovarian cancer, the variant classification module <NUM> may apply the following method:
For each possible event H<NUM> = { H<NUM>, H<NUM>, H<NUM>, H<NUM> }, a model-prediction error may be computed as the root mean squared error (RMSE) of the gene SNP variant fraction distribution computed with respect to its closest theoretical variant fraction values according to the following equation: <MAT> where xi is the unbiased variant fraction of the putative heterozygous SNPs I in the set of M putative heterozygous SNPs identified by the variant classification module <NUM> in the gene to analyze and D(poptimal, Hi) is the subset of theoretical variant fraction values predicted for the alteration event hypothesis Hi in the gene to analyze under the assumption of purity =poptimal.

In the example of <FIG>, there is one value for H<NUM> and two symmetrical values for the other hypotheses.

As will be apparent to those skilled in the art of bioinformatics, the variant classification module <NUM> may also apply other methods than the RMSE to identify the most likely alteration event hypothesis. For instance, when a mixture model has been used in step <NUM>, Mixture Model Classifiers as known in the art may be applied by the variant classification module <NUM> to identify the most likely alteration event hypothesis.

In a possible further disclosure, the information about the gene LOH status and sample purity may also further facilitate the identification of intragenic germline CNVs. The coverage signal of each patient sample in the multiplexed sequencing pool of samples may be normalized and a list of intragenic CNVs may be detected using for instance the method of patent application <CIT>. Due to the intrinsic heterogeneity of the cancer cell sub-clonal variants, the presence of a mix of germline and somatic cells in the tumor samples and the FFPE targeted sequencing coverage of generally low quality, the coverage signal is particularly noisy and it is difficult to directly detect the CNV variants without further knowledge on the underlying sample. In order to establish the relevance of each intragenic germline CNVs, the coverage signal may be further analyzed taking into account the information about the gene LOH status and sample purity p as formerly derived from the putative germline heterozygous SNP variant faction analysis.

While the above examples have been described for simple somatic event models of copy number values between <NUM> and <NUM>, it will be apparent to those skilled in the art of bioinformatics that more complex scenarios may also be modelled in a similar way, for instance with copy number values larger than <NUM>, in particular in cases where more genes are covered by the experimental panel.

The variant classification module <NUM> may then combine detected information on likely pathogenic SNPs, INDELs, intragenic germline CNVs and LOH status to classify and report them to facilitate their interpretation. As will be apparent to those skilled in the art of oncogenomics, in the case of tumor suppressor genes, the primary information of interest to the end user is whether the gene has lost its function. This loss of function may be due to a combination of at least two genomic alterations. Somatic variant databases such as Clinvar record known SNP and INDEL variants of likely pathogenicity. Some of those variants, regardless of their germline or somatic origin, may be combined with chromosomal aberrations in tumor cells, causing the biallelic loss of function of the gene. In the exemplary application of ovarian cancers targeted treatment with oliparib, as recently published by <NPL>, the secondary information of interest is whether the detected gene loss of function is due to:.

In a possible embodiment, the variant classification module <NUM> may annotate and prioritize <NUM> the variant calling data of each patient tumor sample by using information from cancer variant databases such as COSMIC or Clinvar, and/or pathogenicity prediction tools such as Align GVGD, Sift, Polyphen2, Alamut and others to identify SNPs and INDELs of likely pathogenicity, as will be apparent to those skilled in the art of bioinformatics applied to personalized medicine. The variant classification module <NUM> may further determine how these mutations have been combined in the cancer cells with the chromosomal arrangements as previously inferred to determine, for each gene to be analyzed, a combination of genomic alterations possibly causing a biallelic loss of function of the gene.

In the exemplary application of ovarian cancer analysis, for each gene in each sample, the analyzed genes may be classified in six different categories, based on their resulting genomic aberrations, as represented in <FIG>):.

In some samples, the detected information on likely pathogenic SNPs, INDELs, intragenic germline CNVs and LOH status may not be sufficient for the variant classification module <NUM> to classify and report the gene into a specific category. For instance, the variant classification module <NUM> may identify positive LOH status without confirming whether the likely pathogenic variant is of germline (class C) or somatic (class D) origin. The variant classification module <NUM> may accordingly categorize those variants as C/D (class C or class D).

In a possible embodiment, the above classes may be used by the variant classification module <NUM> to also determine the loss-of-function allelic type:.

In a possible embodiment, no distinction in made between germline-derived, somatically acquired and subclonal variants, but other embodiments are also possible, for instance by discriminating them into sub-classes A1, A2, A3. and B1, B2, B3. , subject to the theoretical variant fraction distribution scenarios used in the initial variant classification analysis steps.

The variant classification module <NUM> may then report <NUM> the detected genomic variant classes to the end user in various formats. In a possible embodiment, a list may be used to report for each gene the actual detection status (<NUM>) biallelic (<NUM>) potentially biallelic (<NUM>) monoallelic or (<NUM>) undetermined, and the class A, B, C, D, E, F or undetermined. In a possible embodiment illustrated on <FIG>), a table may be used to report for each gene (here BRCA1, BRCA2 and TP53) and for each sample in the pool with a color code the actual detection status as well as the detected class, marked into the table cell - for instance in <FIG>) dark cells represent monoallelic status, clear cells represent potentially bi-allelic status, marked medium grey cells marked represent bi-allelic status, and unmarked medium grey cells represent undetermined status. Grey crosses further mark samples in which homologous recombination repair is likely defective and sensitive to oliparib as a result of BRCA1/BRCA2 biallelic variants, as detected by the proposed methods.

<FIG> illustrates a possible automated workflow of a genomic data analyzer <NUM> adapted to analyze <NUM> with a variant classification module <NUM> the next generation sequencing variant calling information from somatic patient samples to report likely pathogenic biallelic genomic alterations of analyzed genes in a genomic alteration report, comprising:.

In a further disclosure, the genomic data analyzer may also be adapted to identify <NUM> intragenic germline CNVs for each sample and also take them into account in identifying <NUM> biallelic genomic alterations for each gene of each sample.

The genomic data analyzer <NUM> does not need germline reference samples to operate the variant classification analysis and instead can operate solely on variant calling information from samples comprising a mix of germline and somatic cells of unknown purity, such as FFPE tumour samples. This is a significant advantage over prior art published methods which separately analyze germline and somatic variants, both in terms of routine diagnosis costs and in terms of genetic counselling as may be required for germline analysis due to the possible relationship of germline mutations with other pathologies in the patient and/or its siblings.

As will be apparent to those skilled in the art, another limitation of prior art methods is that the genomic data signals out from sequencing such samples are of low quality and particularly challenging to analyze without proper data pre-processing. In particular, as we experimented with prior art bioinformatics analysis methods applied to somatic FFPE tumour samples sequencing data without germline sample references, the coverage signal alone was not sufficient for the genomic data analyser <NUM> to reliably infer the presence of somatic CNVs that affected entire genes. The reasons for this are: i) the frequent presence of somatic CNVs in ovarian cancer patients; ii) the fact that, in some samples, the coverage signal was affected by an important presence of noise, which was likely induced by formalin-fixed paraffin-embedding; iii) the fact that somatic samples do not simply comprise a homogeneous population of cancer cells, but also include an unknown fraction of germline cells. Together with the absence of reference samples, these factors prevented us from performing an appropriate coverage normalization. In contrast, the proposed methods grounded on statistical modelling and proper processing, with a variant classification data processor module, of the variant fraction distribution data enables to overcome these limitations and to detect chromosomal aberration events such as LOH induced by somatic deletions or duplications. Importantly, since our approach fully relies on processing the variant fraction information rather than the coverage signals, copy-neutral LOH events can also be detected.

The genomic variant classification report may further facilitate personalized diagnosis and treatment of pathologies such as cancers by identifying whether the sample alterations are biallelic solely from next generation sequencing variant calling information without requiring complementary assays, even with low coverage signals and without the need for dedicated variant phasing bioinformatics methods in the genomic data analyzer <NUM> workflow, thus enabling lower cost genomic analysis than prior art methods and assays.

While the proposed methods have been described as a standalone workflow, it will be apparent to those skilled in the art of bioinformatics that for certain applications they may also be combined or complemented with other methods, in parallel with the proposed workflow or as additional steps, in order to further improve the specificity and sensitivity of the overall genomic analysis results. For instance, the proposed method may complement the method of patent application <CIT> to better detect germline CNVs from somatic samples, but other embodiments are also possible.

<FIG> shows the results obtained with the proposed variant detection method out of variant calling information resulting from the genomic analysis of an FFPE ovarian cancer patient sample with targeted next generation sequencing on the BRCA1, BRCA2 and TP53 genes. Putative germline heterozygous SNPs were first identified from the variant calling information and their variant fraction distribution was unbiased, as shown in <FIG>) where the black dots represent the individual putative germline heterozygous SNPs and the grey area their variant fraction distribution. Vertical dashed lines show the theoretical variant fraction values predicted for germline-derived SNPs affected by four different alteration event hypotheses, namely: no loss of heterozygoty ("NO-LOH"), copy-neutral LOH ("CN-LOH"), deletion-induced LOH ("LOH-DEL") and duplication-induced LOH ("LOH-DUP").

The unbiased variant fraction distribution of <FIG>) exhibits variant fraction peaks which can be mapped to a scenario of a sample purity of <NUM>% as the value minimizing the prediction-error of a statistical model of the variant fraction distribution, as illustrated in <FIG>). This corresponds to a somatic alteration scenario characterized by the presence of LOH in all genes BRCA1, BRCA2 and TP53 as illustrated in <FIG>). While in BRCA1 and BRCA2 LOH was induced by a somatic deletion (multiple SNPs mapped around the symmetrical <NUM>% and <NUM>% variant fraction peaks corresponding to the "LOH-DEL" theoretical variant fraction model for a purity of <NUM>%), we found that TP53 underwent copy-neutral LOH (multiple SNPs mapped around the symmetrical <NUM>% and <NUM>% variant fraction peaks corresponding to the "CN-LOH" theoretical variant fraction model for a purity of <NUM>%).

This information may be further used for detecting the most relevant variants of likely pathogenicity in the cancer cells as follows. First, most of the variant fractions observed in this particular sample correspond to theoretical scenario predictions made for germline-derived variants (<FIG>) solid gray lines). These variants may thus be classified as germline-derived. Furthermore, a positive LOH status may be attributed to germline-derived variants whose variant fraction is larger than <NUM>% (<FIG>) variants v1, v3 and v5). In contrast, germline-derived variants with variant fraction smaller than <NUM>% (<FIG>) variants v2, v4 and v6) correspond to lost variants, carried by the allele that was somatically deleted. Since these variants are no longer present in the somatic cells, their LOH status may be classified as negative. Interestingly, the variant fraction of a likely pathogenic TP53 SNP2 (<FIG>), v7) does not match any of the theoretical predictions made for germline-derived variants, but rather corresponds with the theoretical value expected in cases where a somatically acquired SNP was subsequently affected by copy-neutral LOH (<FIG>) black line). Overall, this analysis allows us to infer that cancer cells from this sample are likely affected by a biallelic TP53 variant.

Variant fraction analysis for the same scenario models as illustrated in <FIG> was performed on a pool of sixty individual samples, resulting in the identification and classification of <NUM> variants most likely pathogenic. The number of likely pathogenic variants identified in each individual sample varied between <NUM> and <NUM>. The majority of these variants were found to be missense SNPs and frameshift INDELs. Variant origin and LOH status could be determined respectively in <NUM>% and <NUM>% of likely pathogenic variants. Our results revealed the presence of <NUM> germline-derived variants affected by LOH, <NUM> somatic variants affected by LOH as well as <NUM> somatic subclonal variants of likely pathogenicity. Amongst the <NUM> variants whose origin could not be identified, <NUM> had positive LOH status. We found that TP53 carried at least one likely pathogenic in about <NUM>% of the samples, while <NUM>% of the samples carried a likely pathogenic BRCA2 mutation of germline origin. Finally, among the sixty samples that were analyzed, <NUM> carried a germline- derived, likely pathogenic BRCA1 variant.

A further analysis also enabled to identify germline CNVs on the same patient as illustrated on <FIG>). <FIG>) shows an example of a normalized coverage signal for BRCA2, TP53 and BRCA1 for the variant calling information for this patient sample. Each dot shows the average normalized coverage observed in a single exon. Gray and white areas correspond to different exons. The normalized coverage has arbitrary units (a. ) and cannot be directly interpreted as the total number of alleles (i.e., copy number), but it still enables to infer relevant information: the exon17 genomic region of gene BRCA1 shows a significant coverage drop. The normalized coverage ratio between exon <NUM> and the other exons in BRCA1 is measured as Crel = <NUM>%. Under the assumption that this sample comprises a mix of a homogenous population of cancer cells and a homogeneous population of germline cells, the relative coverage observed in exon <NUM> is given by: <MAT> where Sex17 and Gex17 are the number of copies of exon <NUM> present in somatic and germline cells respectively, SBRCA1 and GBRCA1 are the number of copies of the other BRCA1 exons present in somatic and germline cells, respectively, and p denotes the sample purity. The sample purity may be inferred as the value minimizing the prediction error of the mathematical model accounting for the variant fraction of putative heterozygous germline SNPs found in BRCA1, BRCA2 and TP53, in this experiment p=<NUM>% as shown in <FIG>).

The BRCA1 status of this sample <NUM> is represented schematically in <FIG>). In germline cells (top), BRCA1 was present in two copies (GBRCA1=<NUM>), except for exon <NUM>, which had only one copy (Gex17=<NUM>, dark grey) as a result of a heterozygous germline deletion (soft grey). In somatic cells (bottom), BRCA1 was only present in one allele (SBRCA1=<NUM>), as a result of a somatic full-gene deletion of the non-mutated allele (black dotted line), as formerly inferred from the putative germline heterozygous SNP analysis. Moreover, the only allele left in cancer cells carried an exon <NUM> deletion (Sex17=<NUM>). Biologically, this means that BRCA1 was affected by a germline heterozygous deletion of exon <NUM> as well as by a somatic deletion of the non-mutated allele. Overall, these results indicate that in this patient cancer cells, both alleles list the BRCA1 tumor suppressor function, as reported in the Sample <NUM> column of the variant classification report table for the pool of <NUM> samples (<FIG>).

Overall, in the experiment the proposed method allowed us to identify <NUM> SNPs, <NUM> INDELs and <NUM> intragenic germline CNVs, which could explain the underlying genetic basis of the disease. Inferring sample purity, further allowed us to identify variant origin (i.e., germline vs. somatic) and perform LOH detection. Together with variant annotation information obtained from public databases to identify likely pathogenic variants, these results allowed us the establish the gene mutation status of BRCA1, BRCA2 and TP53 in each individual sample. Amongst the <NUM> patients analysed in our experiment, <NUM> biallelic pathogenic variants were identified in the BRCA1 gene (<NUM>%), <NUM> in the BRCA2 gene (<NUM>%) and <NUM> in the TP53 gene (<NUM>%). Moreover, in <NUM> samples, none of the genes that were analysed carried biallelic pathogenic variants. We also observed that in the BRCA1/BRCA2 genes, biallelic pathogenic variants were in most cases germline-derived variants which underwent LOH while in the TP53 gene, they were in most cases somatic variants which underwent LOH. Overall, this analysis allowed us to identify <NUM> cases in which patients would likely benefit from a PARP inhibition therapy such as olaparib to treat their ovarian cancer tumour.

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 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 of BRCA1, BRCA2 and TP53 genes, it will be apparent to those skilled in the art of bioinformatics that they also apply to the detection and classification of variants of other tumor suppressor or homologous repair genes such as for instance ATM, CHEK1, CHEK2, BARD1, BRIP1, RAD51C, RAD51D, FAM175A, MRE11A, NBN, PALB2, APC, DCC, DF1, NF2, PTEN, Rb, VHL, WT1, PP2A, LKB1, or INK4a/ARF, and others which may be identified in the future as the scientific oncogenomics knowledge develops.

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 structures", "data fields", "file", or "stream" may thus be used indifferently throughout this specification.

As will be apparent to those skilled in the art statistics, the methods described herein may be indifferently applied to various statistical methods such as probability representations and statistical measurements. The terms "distribution", "likelihood", "probability" 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.

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. In fact, after reading the above description, it will be apparent to one skilled in the relevant art(s) how to implement alternative 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.

Although the term "at least one" may often be used in the specification, claims and drawings, the terms "a", "an", "the", "said", etc. also signify "at least one" or "the at least one" in the specification, claims and drawings.

Certain embodiments are described herein as including logic or a number of components, modules, units, or mechanisms. Modules or units may constitute either software modules (e.g., code embodied on a machine-readable medium or in a transmission signal) or hardware modules. A "hardware module" is a tangible unit capable of performing certain operations and may be configured or arranged in a certain physical manner. In various example embodiments, one or more computer systems (e.g., a standalone computer system, a client computer system, or a server computer system) or one or more hardware modules of a computer system (e.g., a processor or a group of processors) may be configured by software (e.g., an application or application portion) as a hardware module that operates to perform certain operations as described herein.

In some embodiments, a hardware module may be implemented mechanically, electronically, or any suitable combination thereof. For example, a hardware module may include dedicated circuitry or logic that is permanently configured to perform certain operations. For example, a hardware module may be a special-purpose processor, such as a field-programmable gate array (FPGA) or an ASIC. A hardware module may also include programmable logic or circuitry that is temporarily configured by software to perform certain operations. For example, a hardware module may include software encompassed within a general-purpose processor or other programmable processor.

Claim 1:
A method for identifying the presence of a chromosomal aberration in the tumor sample cells in at least one of several genes from the next generation sequencing data analysis of a pool of patient tumor samples comprising at least two analyzed genes, with a variant classification data processor module which operates solely on variant calling information from samples comprising a mix of germline and somatic cells and for each possible sample purity value p, in which a combination of at least two genomic alterations have caused its biallelic loss of function, the said data processor module performs a method comprising:
- acquiring a variant calling information of the next generation sequencing of the pool of patient tumor samples comprising germline and somatic cells, the variant calling information comprising for each patient a list of SNPs genomic variants in at least two genes to be analyzed;
- selecting, in the list of genomic variants, a set of putative germline heterozygous SNPs by comparing the list of genomic variants to a reference human genome variant database;
- for each possible value of sample purity p, a sample purity p being a ratio of somatic DNA material to germline DNA material in the overall FFPE sample, fitting a mixture model to the observed variant fraction distribution of said set in the variant calling data, the modeled variant fraction distribution values being calculated as a function of the sample purityp;
- inferring an optimal sample purity poptimal as the possible sample purity p for which the mixture model best fits the observed variant fraction distribution of said set;
- for each gene to be analyzed:
o selecting, in the plurality of putative germline heterozygous SNPs, a subset of putative germline heterozygous SNPs located in the gene to be analyzed,
o inferring the presence of a chromosomal aberration in said tumor sample cells by comparing, on at least said subset, the observed variant fraction distribution with the predicted variant fraction distribution for each possible chromosomal aberration, the predicted variant fraction distribution values being calculated as a function of the optimal sample purity poptimal.