Richard Iggo
email: r.iggo@bordeaux.unicancer.fr

Justine Rudewicz

Elodie Monceau

Nicolas Sevenet

Jonas Bergh

Tobias Sjoblom

Hervé Bonnefoi

Validation of a yeast functional assay for p53 mutations using clonal sequencing

Keywords: p53 yeast functional assay, breast cancer, next generation sequencing

We have previously tested biopsies from 1469 breast tumours with a p53 functional assay in the context of a prospective clinical trial (EORTC 10994/BIG 1-00). The goal of the trial was to determine whether p53 status could be used to select patients who would benefit from inclusion of taxanes in anthracycline-based chemotherapy. The results of the trial were negative. To test whether this was because the functional assay misclassified the tumours, we have reanalysed two groups of biopsies by Sanger sequencing and Roche 454 next generation sequencing (NGS). Comparison of yeast data with pooled cDNA sequencing data in an initial cohort of 69 biopsies showed that conventional sequencing is insensitive when the mutant p53 content is low. A second cohort of 48 biopsies was used to compare directly the yeast assay with Sanger and NGS technology. The mutant sequence was difficult to detect in sequence chromatograms of pooled cDNA, whereas NGS unequivocally identified mutations in every case classified as mutant by the functional assay. The NGS data showed that small deletions, probably caused by PCR splicing, account for most of the unexplained background in the yeast assay. We conclude that mutation detection techniques that test multiple clones, such as the p53 functional assay and NGS, are more reliable than Sanger sequencing of pooled DNA; that the high p53 mutation rate (44%) seen with the yeast assay in the EORTC 10994/BIG 1-00 trial reflects this high sensitivity; and that NGS with Roche 454 technology could be used to identify the p53 mutations in the remaining tumours previously tested in yeast in the EORTC10994/BIG 1-00 trial.

Introduction

The p53 gene (TP53 ) is one of the most commonly mutated genes in human tumours [START_REF] Petitjean | Impact of mutant p53 functional properties on TP53 mutation patterns and tumor phenotype: lessons from recent developments in the IARC TP53 database[END_REF][START_REF] Goh | The role of mutant p53 in human cancer[END_REF]. A longstanding model posits that specific p53 mutations are linked to drug response [START_REF] Aas | Specific P53 mutations are associated with de novo resistance to doxorubicin in breast cancer patients[END_REF]. Given the established role that p53 plays in the response to chemotherapeutic drugs in experimental models, we organized a multi-centric phase III clinical trial in 1999 to test whether p53 status influences the response of human breast tumours to chemotherapy. Our hypothesis was that p53 mutant tumours would be resistant to anthracyclines but would remain sensitive to taxanes. The trial, EORTC 10994/BIG 1-00, compared a regimen containing taxanes with one without taxanes [START_REF] Bonnefoi | TP53 status for prediction of sensitivity to taxane versus non-taxane neoadjuvant chemotherapy in breast cancer (EORTC 10994/BIG 1-00): a randomised phase 3 trial[END_REF]. The results of the clinical trial did not confirm our hypothesis. We used a yeast functional assay to determine p53 status in the clinical trial. The yeast assay distinguishes inactivating mutations from biologically unimportant passenger mutations [START_REF] Flaman | A simple p53 functional assay for screening cell lines, blood, and tumors[END_REF] but provides no information on the identity of the mutations, so we were unable to ask whether specific p53 mutations play a role in chemotherapy resistance. To address this question, we have now performed a feasibility study to determine whether next generation sequencing (NGS) technology could be used to identify the p53 mutations present in the cDNA that was originally used to test p53 status in the clinical trial.

We show that NGS is more reliable than direct sequencing when the mutant p53 mRNA fraction in the tumour is low. The reason for the high sensitivity of the yeast assay and NGS is that they test clones, which means that wild-type sequences derived from normal tissue do not obscure the mutant signal.

Materials and methods

Consent

p53 sequencing was a secondary endpoint of the original EORTC 10994 study for which all patients gave signed informed consent as a condition for enrolment in the study [START_REF] Bonnefoi | TP53 status for prediction of sensitivity to taxane versus non-taxane neoadjuvant chemotherapy in breast cancer (EORTC 10994/BIG 1-00): a randomised phase 3 trial[END_REF].

DNA sequencing

The PCR products used to test p53 status in yeast in the original study [START_REF] Bonnefoi | TP53 status for prediction of sensitivity to taxane versus non-taxane neoadjuvant chemotherapy in breast cancer (EORTC 10994/BIG 1-00): a randomised phase 3 trial[END_REF] were split into two fragments for NGS. Tags and 454 sequencing primers (Supplementary Table 1) were added by PCR as shown in Supplementary Figure 1 (see the Supplementary materials and methods for details). A pool containing 96 tagged PCR products was sequenced in a single batch on a GS FLX Amplicon sequencer using Titanium series chemistry (Roche 454 Life Sciences, Branford, CT, USA). Plasmids and PCR products were analysed by conventional Sanger sequencing as described in refs 4 and 6 and in the Supplementary materials and methods.

The sequencing data (fastq files) are available in the NCBI SRA database under accession number SRP020456 (BioProject PRJNA193388).

NGS analysis pipeline

The pipeline for calling variants is shown in Figure 1. The fasta and qual files were merged to give fastq files that have been deposited in the NCBI SRA database under accession number SRP020456 (BioProject PRJNA193388). The alignment and variant calling pipeline is described in the Supplementary materials and methods.

Yeast functional assay

Several mutants identified in this study were retested by yeast functional assay as previously described [START_REF] Flaman | A simple p53 functional assay for screening cell lines, blood, and tumors[END_REF]. To test for temperature sensitivity, yeasts were incubated at 25 • C, 32 • C, and 36 • C for 3 days.

Results

Comparison of next generation sequencing to Sanger sequencing

We initially attempted to identify p53 mutations by Sanger sequencing of the RT-PCR products used for the yeast assay. Most mutations were missed by standard mutation detection software when the mutant fraction estimated from the yeast assay was below 40% (Supplementary Table 2). To explore the reasons for this failure, we performed a direct comparison of mutation detection techniques in a fresh cohort of 48 biopsies taken from 31 tumours. The tumours were not selected randomly or sequentially; instead, they were chosen to give a group of nine 'positive control' tumours with 59-92% red colonies in the yeast assay and a group of 12 'negative control' tumours scored as wild type in the yeast assay. To challenge the NGS technology, we selected a third group of ten 'difficult' tumours that gave 21-48% red colonies.

To check for reproducibility, we tested two biopsies per tumour for 17 of the 19 tumours in the positive control and difficult groups. The mutations identified in the plasmids in the 'positive control' group were all missense or in-frame mutations, whereas the mutations in the 'difficult' group were all nonsense, frameshift or splicing mutations, apart from one in-frame deletion Full details of the mutations are given in Supplementary Tables 3 and4. ¶ The plasmids contained a mixture of abnormally spliced products; the same splice variant was seen in one plasmid as in the pooled cDNA and NGS reads (the causal mutation was only visible in the plasmid sequences). * The causal mutation was not identified. † A probable PCR splicing artefact was also identified in biopsy 1. ‡ A splice donor mutation was present in one plasmid and intron retention in the other. § A P191A mutation encoding functional p53 was also present in this case (see text for details). nt = not tested; na = not applicable; wt = no mutations identified; fs = frameshift (the mutant codon is indicated); sd = splice donor.

(Table 1). This is expected because nonsense-mediated mRNA decay (NMD) reduces the level of mRNA containing premature stop codons.

NGS sequence analysis: common artefacts

The base quality, read depth, and coverage of the NGS sequences are summarized in Figure 2 and Supplementary Figure 2. As expected for 454 technology, low phred scores were mainly seen at homopolymers. The sites most commonly affected are indicated by black circles in Figure 2. We previously reported that PCR splicing artefacts can give rise to significant background mutations in the yeast assay when reverse transcription is primed with hexamers [START_REF] Waridel | Field cancerisation and polyclonal p53 mutation in the upper aero-digestive tract[END_REF]. Before starting the EORTC 10994 study, we showed that use of a p53 -specific reverse transcription primer substantially reduced this problem. The NGS data show that small deletions probably caused by PCR splicing still occurred to some extent in all samples in exons 4-6 (Figure 3A). To quantify the phenomenon, we generated a 'PCR splicing score' by counting the number of reads containing deletions present in at least two tumours (Supplementary Table 1). To visualize the deletions, we plotted deletion size against position in the p53 cDNA (Figures 3B and3C).

The size of the plotted points reflects the number of reads containing the deletion. Presented in this way, it is easy to distinguish high-quality samples (Figure 3B) from those with frequent PCR splicing artefacts (Figure 3C). These artefacts were easily removed by the filters in the analysis script and did not generate false positives in the final analysis.

NGS allele depth versus % red colonies in the yeast functional assay

Figure 4 shows the fraction of mutant reads compared with the fraction of red colonies in the yeast assays.

There is excellent agreement between the two measures of mutant fraction (r = 0.97). The y-intercept is 15.1, indicating that there are no mutant reads in the NGS data when the % red colonies is 15.1%. This is slightly higher than our estimate of the background based on the % red colonies in cases scored as wild type in the published study (11 ± 4%) [START_REF] Bonnefoi | TP53 status for prediction of sensitivity to taxane versus non-taxane neoadjuvant chemotherapy in breast cancer (EORTC 10994/BIG 1-00): a randomised phase 3 trial[END_REF], but nevertheless remarkably similar given the fact that the current estimate is derived from NGS of a small selected population of mutant cases, whereas the previous estimate was derived from yeast assays of the full population of wild-type cases.

NGS sequence analysis: 'positive control' group

In the group of samples with high mutant RNA content, the same mutations were found in the plasmid sequences as in the NGS sequences, and the same mutations were found in both biopsies of a pair (Table 1 and Supplementary Table 3). We conclude that NGS faithfully detects abundant mutations, as expected.

NGS sequence analysis: 'negative control' group

The analysis pipeline reported no mutations in 11 of the 12 samples in the negative control group. The final sample had a P222L variant in 81% of the reads (Table 1 and Supplementary Table 3). The average base phred score was 39, indicating that the variant was called with high confidence. This variant is labelled as an SNP (rs146340390) but the dbsnp minor allele frequency is so low (0.000) that it could be a mutation. The reason for using the yeast assay was to distinguish functionally important mutations from polymorphisms and passenger mutations. To test whether the P222L mutant was partially active, we tested it in yeast at different temperatures. The colonies were white at 25-36 • C, but after storage at 4 • C they gradually became slightly pink (Figure 5), indicating that the P222L mutant protein has a mild transcriptional defect.

We conclude that NGS identifies wild-type samples with high confidence but additional filtering is required to distinguish functionally important mutations from polymorphisms and passenger mutations.

NGS sequence analysis: 'difficult' group

The 'difficult' tumours included four with nonsense, frameshift or in-frame deletion mutations and two for which we could not confidently identify the mutation in the plasmid sequences. We identified nonsense or frameshift mutations by NGS in all ten biopsies from these tumours (Table 1 and Supplementary Table 3). Nine of these biopsies had yeast results in the range 21-29% red colonies, close to the 20% cut-off used to score tumours as mutant. NGS thus demonstrates that tumours just above the cut-off in the yeast study were correctly scored as mutant. The remaining four tumours in the 'difficult' group had splicing mutations. These mutations commonly lead to a mixture of intron retention and exon skipping; the actual mutation is typically only visible in the former. This is exactly what was seen in the NGS data. The mutant base was present in the plasmid sequences from all four tumours but only present in the NGS data in three of them. In the fourth tumour (183, Supplementary Table 3), the mutation was located five nucleotides after the splice site (cgˆgtcag>cgˆgtcac). It was visible in the three plasmids with intron retention. In contrast, it was not present in the NGS sequences because none contained the retained intron, which increases the size of the fragment by 757 bp. Instead, there were multiple reads using a cryptic donor site in exon 4; this variant was also seen in one plasmid. The cryptic donor had the sequence agˆgt. Since this corresponds to a normal splice donor, the alignment software labelled it as a normal intron. To solve this problem, we aligned the NGS sequences to specific p53 mRNA isoforms and forced the software to report unaligned bases as deletions. The mutation was efficiently called with this strategy (Table 1 and Supplementary Table 3). The alternate splice in case 83 is caused by a mutation in the donor site at the end of the alternative exon that improves the donor site (ttˆgt>tgˆgt), leading to more frequent inclusion of the alternative exon. All four plasmids contained the mutation but it was only called in one of the biopsies in the NGS data. The obstacle to calling the mutation in the other biopsy was the low allele depth (two reads with the beta and two with the gamma isoform). The splicing mutations in the two remaining cases were correctly identified (Table 1 and Supplementary Table 3). In both cases, the mutation affected the intron 6 splice donor site, resulting in the use of a cryptic donor (tgˆgt) five nucleotides after the normal donor. The pathologically spliced forms were not called as such after alignment to the genomic sequence, although the mutation itself was called, but they were correctly identified in the alignment to mRNA. In summary, all of the mutations in the 'difficult' group were correctly identified by NGS.

Direct comparison of NGS with Sanger sequencing

To allow direct comparison of NGS with Sanger sequencing, we performed Sanger sequencing of the pooled cDNA PCR products in the second cohort and asked an accredited diagnostic lab geneticist to identify the mutations by visual inspection of the chromatograms. Despite the low mutant content in some samples, all but two of the mutations were identified correctly (Table 1 and Supplementary Table 4). The two mutations missed were present in the NGS data in 5% of the reads. Figure 6 shows three examples of chromatograms with low mutant content in which the mutation was nevertheless correctly identified. The phred scores calculated with Tracetuner software are shown beneath the mutant peaks. They indicate that the probability of error in automated calling of the mutant bases was commonly in the range 10-30%. This should be compared with the phred = 30 cut-off for calling mutant peaks in the NGS data, which corresponds to a probability of error of 0.1%. We conclude that NGS allows much more confident identification of mutations than Sanger sequencing in samples heavily contaminated with normal tissue.

Discussion

The main conclusion from this study is that mutation detection techniques that test multiple clones, such as the yeast assay or NGS, are more sensitive than direct sequencing of bulk tumour DNA. The NGS data also indicate that the tumours in the EORTC 10994 clinical trial scored as mutant by the yeast assay probably do contain authentic p53 mutations. The mutation rate based on the yeast assay was 44%, a higher figure than was anticipated by the study design [START_REF] Bonnefoi | TP53 status for prediction of sensitivity to taxane versus non-taxane neoadjuvant chemotherapy in breast cancer (EORTC 10994/BIG 1-00): a randomised phase 3 trial[END_REF]. Analysis of the distribution of red colonies showed two mutant peaks, one at 76% and the other at 29% [START_REF] Bonnefoi | TP53 status for prediction of sensitivity to taxane versus non-taxane neoadjuvant chemotherapy in breast cancer (EORTC 10994/BIG 1-00): a randomised phase 3 trial[END_REF]. The latter was uncomfortably close to the cut-off for scoring tumours as mutant (20% red colonies). NGS shares with the yeast assay the fact that it tests multiple clones, but has the advantage that background mutant clones (reads) can be distinguished from one another by their sequence. This means that the cut-off can be pushed much lower in studies using NGS than in studies using the yeast assay. The NGS results indicate that samples near the threshold in the yeast assay contained genuine mutations. Taken together, the results indicate that NGS is more reliable than Sanger sequencing of pooled cDNA when the mutant RNA content is low, as is often the case in clinical samples. We suspect that many studies in the older literature may have underestimated the p53 mutation rate for this reason.

In the conditions used for the clinical trial, we expected PCR mutations to generate a background of about 3.5% red colonies [START_REF] Flaman | A rapid PCR fidelity assay[END_REF] but the observed background was 11% ± 4% [START_REF] Bonnefoi | TP53 status for prediction of sensitivity to taxane versus non-taxane neoadjuvant chemotherapy in breast cancer (EORTC 10994/BIG 1-00): a randomised phase 3 trial[END_REF]. The NGS data show us the likely source of much of the unexplained background. The PCR analogue of DNA repair by microhomology-mediated end joining is called PCR splicing. It occurs when template breaks terminate replication before the polymerase reaches the opposing primer site. The termini of the resulting fragments anneal to one another in ensuing cycles at sequences with a few bases of homology. We previously reported that this was a major source of background in the yeast assay when reverse transcription was primed with hexamers [START_REF] Waridel | Field cancerisation and polyclonal p53 mutation in the upper aero-digestive tract[END_REF]. It is substantially reduced by priming reverse transcription with a p53 -specific primer, but the NGS data show that it still occurred. The diversity of the variants generated by PCR splicing means that they were almost all removed automatically by the filters in Varscan, but to make the pipeline more robust we used a Perl script to explicitly remove the commonest PCR splice variants. The Varscan filters in the script are stringent for base quality (phred ≥ 30) but relatively permissive for mutant read depth (≥ 5 mutant reads, ≥ 4.5% mutant reads) because we were interested in identifying mutants that scored near the background in the yeast assay. At these thresholds, there was one case (276) in which the script returned two mutations. The second mutation (P191A) was present in 11% of the NGS reads, but is probably a passenger mutation, a PCR mutation, or a sequencing artefact. The mutant protein is wild type in the Ishioka functional assay database [START_REF] Kato | Understanding the function-structure and function-mutation relationships of p53 tumor suppressor protein by high-resolution missense mutation analysis[END_REF], and the mutation is adjacent to a deletion with a low phred score, indicating that the NGS data are potentially unreliable at this base.

One argument for using the yeast assay in the clinical trial was that it would distinguish functionally important mutations from passenger mutations and polymorphisms. Interestingly, the negative control group in the present study contained a mutation with exactly those properties. To exclude the possibility that the sample was wrongly scored, we retested it in yeast and it to have a very mild transcriptional defect. According to a p53 function predictor [START_REF] Soussi | MUT-TP53 2.0: a novel versatile matrix for statistical analysis of TP53 mutations in human cancer[END_REF], the P222L mutation should have 24% activity, a value that would place it in the mutant camp, but the Ishioka database shows it to be active on most binding sites [START_REF] Kato | Understanding the function-structure and function-mutation relationships of p53 tumor suppressor protein by high-resolution missense mutation analysis[END_REF], consistent with our results (Figure 5).

We chose to sequence the p53 cDNA PCR products used for the clinical trial because genomic DNA was not available. Alignment to genomic DNA with GMAP worked perfectly but the variant calling software was unable to identify pathological splice variants in the GMAP output because the cryptic donor sites were similar or identical to normal splice sites, so the corresponding intron bases were set by GMAP to N in the CIGAR string. The variants were definitely abnormal because we could see the causal mutations in the plasmids and they do not correspond to any of the rare variants reported in a high-depth RNA-seq study on p53 [START_REF] Mercer | Targeted RNA sequencing reveals the deep complexity of the human transcriptome[END_REF]. To identify them, we aligned the NGS reads to the p53 mRNA and forced GMAP to report unaligned bases as deletions. A more pernicious problem was the loss of fragments containing retained introns or additional exons. This reduced the mutant allele depth for one mutant (the p53 beta splice donor mutation in case 83) to the range that overlaps with trivial NGS artefacts. It remains to be seen whether this problem is tractable because allele depth is such a useful criterion for filtering NGS data.

Our ultimate goal is to sequence p53 from the 1469 tumours tested in yeast in the EORTC 10994 trial, to find out whether specific p53 mutations predict the response of tumours to chemotherapy. We chose to test 454 NGS technology because we could incorporate tags for 96-fold multiplexing in a simple PCR step starting from our existing PCR products. While we were able to identify all of the mutations that we sought, few scientists could view the profile in Figure 2 with total equanimity. In addition to its inability to sequence homopolymers, 454 technology has a rather low capacity by current standards. Were the read length of competing technologies to reach ∼550 bp, libraries based on our existing strategy could easily be further multiplexed to sequence all 1469 tumours in one run.

The three main conclusions from this study are that it would be possible to identify the p53 mutations in the tumours from the EORTC 10994 clinical trial with the protocol described here; that the tumours from the trial that were scored as mutant in the yeast assays [START_REF] Bonnefoi | TP53 status for prediction of sensitivity to taxane versus non-taxane neoadjuvant chemotherapy in breast cancer (EORTC 10994/BIG 1-00): a randomised phase 3 trial[END_REF] probably do contain authentic p53 mutations; and that studies from the pre-NGS era probably underestimate p53 mutation rates.

Figure 1 .

 1 Figure 1. Pipeline for sequence analysis. Plasmids rescued from yeast were sequenced to identify the mutations present in the biopsies. The plasmid data were compared with data from direct sequencing of the PCR products or NGS from the PCR products.

Figure 2 .

 2 Figure 2. Sequence quality showing the effect of homopolymers in p53. Mean base quality (phred scale) and read depth are plotted against position in p53. The continuous line shows base quality. The grey silhouette shows read depth. The domain structure of p53 is shown below the plot: TAD1, transactivation domain 1; Pro, proline-rich domain; DBD, DNA binding domain; N, nuclear localization signal; Tet, tetramerization domain; RD, regulatory domain. The region sequenced spanned codons 54-348. The black circles show 16 bases in homopolymers where the mean quality was below 25, affecting codons 89, 151, 152, 153, 176, 177, and 178.

Figure 3 .

 3 Figure 3. Artefactual deletions caused by PCR splicing. PCR splicing creates artefactual deletions in exons 4-6. (A) Reads from a sample classified as wild type (158) displayed in the Integrated Genome Viewer (IGV). The thick bars in the bottom panel show p53 coding exons. The light grey bars in the middle panel show sequencing reads. Gaps in reads that do not begin at exon boundaries are probably caused by PCR splicing. (B, C) Plots showing the position and length of deletions in the biopsies from tumour 340. Biopsy 2 has few deletions (B), whereas biopsy 1 has frequent deletions (C). The size of the dots is proportional to the number of mutant reads. Deletions present in at least two tumours were used to generate a PCR splicing score (see Supplementary Table1).

1

 1 Figure 3. Artefactual deletions caused by PCR splicing. PCR splicing creates artefactual deletions in exons 4-6. (A) Reads from a sample classified as wild type (158) displayed in the Integrated Genome Viewer (IGV). The thick bars in the bottom panel show p53 coding exons. The light grey bars in the middle panel show sequencing reads. Gaps in reads that do not begin at exon boundaries are probably caused by PCR splicing. (B, C) Plots showing the position and length of deletions in the biopsies from tumour 340. Biopsy 2 has few deletions (B), whereas biopsy 1 has frequent deletions (C). The size of the dots is proportional to the number of mutant reads. Deletions present in at least two tumours were used to generate a PCR splicing score (see Supplementary Table1).

  Figure 3. Artefactual deletions caused by PCR splicing. PCR splicing creates artefactual deletions in exons 4-6. (A) Reads from a sample classified as wild type (158) displayed in the Integrated Genome Viewer (IGV). The thick bars in the bottom panel show p53 coding exons. The light grey bars in the middle panel show sequencing reads. Gaps in reads that do not begin at exon boundaries are probably caused by PCR splicing. (B, C) Plots showing the position and length of deletions in the biopsies from tumour 340. Biopsy 2 has few deletions (B), whereas biopsy 1 has frequent deletions (C). The size of the dots is proportional to the number of mutant reads. Deletions present in at least two tumours were used to generate a PCR splicing score (see Supplementary Table1).

Figure 4 .

 4 Figure 4. Estimating the background in the yeast assay. The percentage of NGS reads containing the mutations in Table1is plotted against the percentage of red yeast colonies. Reads containing P191A in sample 276 and P222L in sample 256 were excluded these variants do not give red colonies. Cases 83 and 183 were excluded because mutant reads containing splice variants that increased the length of the template were lost from the NGS data. The NGS data are derived from the alignment to the prototypic p53α mRNA.

Figure 5 .

 5 Figure 5. Yeast functional assay for the P222L p53 mutant. The yeast reporter strain contains an ADE2 gene expressed from a p53-responsive promoter. (A) Yeast transfected with the p53 PCR product from case 256. The pink circle shows a pink colony containing the P222L mutant; the white circle shows a white colony containing wild-type p53 derived from normal tissue; the red circle shows a red colony containing mutant p53 inactivated by a PCR mutation. (B-D) Yeasts expressing some of the mutants identified in this study were grown at the temperatures indicated. The P222L mutant colonies are slightly cold-sensitive. The R110P mutation is labelled as an SNP (rs11540654) but is clearly inactive.

Figure 6 .

 6 Figure 6. Sequence chromatograms showing secondary peaks. Sanger sequencing was performed on pooled cDNA from biopsies 215.2 (A), 272.2 (B), and 326.2 (C). The mutant content of these biopsies estimated from the NGS data was 14%, 7%, and 13%, respectively. The base calls and quality (phred) scores are shown for the wild-type and mutant sequences above and below the chromatograms, respectively. An asterisk in the mutant track indicates that the base was identical in both sequences.

Validation of a p53 functional assay by NGS 443Table 1 .

 1 Mutations identified by NGS compared with mutations identified by Sanger sequencing

EORTC ID Red yeast colonies biopsy 1 (%) Red yeast colonies biopsy 2 (%) No of plasmids with the indicated mutation Sanger plasmid mutations No of biopsies tested by NGS/Sanger cDNA pool sequencing Sanger cDNA pool mutations NGS mutations Mutation type

	169	85	77	3/3	Y163C	2	Y163C	Y163C	Missense
	207	62	73	2/2	R110P	2	R110P	R110P	Missense
	221	84	67	4/4	Y220C	2	Y220C	Y220C	Missense
	279	79	72	4/4	C275R	2	C275R	C275R	Missense
	316	92	92	4/4	F134L	2	F134L	F134L	Missense
	318	67	33	nt	na	2	R282W	R282W	Missense
	320	83	73	4/4	R110 RLGF	2	R110 RLGF	R110 RLGF	In-frame deletion
	340	66	77	4/4	I195S	2	I195S	I195S	Missense
	341	59	86	3/3	R273H	2	R273H	R273H	Missense
	83	37	22	4/4	p53-beta	2	p53-beta*	p53-beta	Splicing
	183	48	45	4/4	intron 4 sd ¶	2	Intron 4 sd ¶	Intron 4 sd ¶	Splicing
	192	23	29	1/4	P64fs	2	P64fs	P64fs	Frameshift
	215	38	25	2/3	S241fs	2	S241fs †	S241fs	Frameshift
	269	40	43	3/4	Intron 6 sd	2	Intron 6 sd	Intron 6 sd	Splicing
	272	36	25	2/3 †	Intron6sd	2	Intron6sd	Intron6sd	Splicing
	276	27	na	2/4	G245 G	1	wt	G245 G §	In-frame deletion
	290	22	21	2/4	S166fs	2	S166fs	S166fs	Frameshift
	323	22	na	2/4	Q165Z	1	wt	Q165Z	Nonsense
	326	22	23	1/4	L145fs	2	L145fs	L145fs	Frameshift
	158	7	na	na	na	1	wt	wt	na
	193	7	na	na	na	1	wt	wt	na
	256	8	na	na	na	1	P222L	P222L	Missense (SNP)
	267	8	na	na	na	1	wt	wt	na
	284	5	na	na	na	1	wt	wt	na
	285	13	na	na	na	1	wt	wt	na
	288	10	na	na	na	1	wt	wt	na
	306	6	na	na	na	1	wt	wt	na
	311	10	na	na	na	1	wt	wt	na
	312	7	na	na	na	1	wt	wt	na
	319	10	na	na	na	1	wt	wt	na
	322	8	na	na	na	1	wt	wt	na

Acknowledgments

We thank Monica Wallin for technical assistance; the bioinformaticians who generously made available their software through public websites; the members of the Institut Bergonie breast cancer group for helpful advice; and Macha Nikolski for critical reading of the manuscript. We also thank the French Cancer League ('Equipe Labellisée Ligue Contre le Cancer 2011'), the French National Research Agency (ANR-08-CEXC-016-01), and the Fondation Widmer pour la Lutte Contre le Cancer for financial support.

Author contribution statement

The authors contributed in the following way: TS, NS, and EM: Sanger sequencing; EM: preparation of libraries; TS, JR, and RI: data analysis; HB, JB, and RI: study design; RI: yeast assays and writing of the manuscript. All of the authors approved the final version of the manuscript.

SUPPORTING INFORMATION ON THE INTERNET

The following supporting information may be found in the online version of this article.

Supplementary materials and methods.

Table S1. Primers used for multiplexing showing the PCR splicing score for each biopsy and the number of reads per tag.

Table S2. Mutations identified by Sanger sequencing in cohort 1.

Table S3. Mutations identified by NGS in cohort 2.

Table S4. Mutations identified by Sanger sequencing in cohort 2.


Philippe Owezarski

Unsupervised Classification and Characterization of Honeypot Attacks

Keywords: Honeypot attack detection, Anomaly classification, unsupervised machine learning, autonomous security systems

Monitoring communication networks and their traffic is of essential importance for estimating the risk in the Internet, and therefore designing suited protection systems for computer networks. Network and traffic analysis can be done thanks to measurement devices or honeypots. However, analyzing the huge amount of gathered data, and characterizing the anomalies and attacks contained in these traces remain complex and time consuming tasks, done by network and security experts using poorly automatized tools, and are consequently slow and costly. In this paper, we present an unsupervised method for classification and characterization of security related anomalies and attacks occurring in honeypots. This as automatized as possible method does not need any attack signature database, learning phase, or labeled traffic. This corresponds to a major step towards autonomous security systems. This paper also shows how it is possible from anomalies characterization results to infer filtering rules that could serve for automatically configuring network routers, switches or firewalls.

I. INTRODUCTION

Monitoring communication networks and their traffic is an important advance for protecting computers systems. Several methods exist. For example, using monitoring devices at the interconnection points of networks allows the analysis of all incoming and outgoing flows. Another method consists in installing honeypots in the network. Honeypots provide more or less emulated services (according to the expected interaction level with the attackers), and allows supervising the use of the provided services by attackers. Thus, gathered data allows security experts on one hand to detect and analyze computer systems weaknesses, and on the other hand to collect useful information on the attacker activities that allows the analysis of their attacking methods, objectives, and strategies. Estimating and analyzing the risk related to illegitimate activities on the Internet are essential for security experts to design and develop adapted and efficient defense and protection systems with regard to the actual risk.

Characterizing and classifying current anomalies and attacks of the Internet are very complex and time consuming tasks, done by experts. Such tasks are therefore slow and costly. The main difficulty related to the identification and analysis of the different classes of illegitimate traffic is their fast evolution, amplification, and renewing capabilities. Designing autonomous classification and characterization processes is crucial for easy to deploy and to use defense systems. Hence, modern classification systems must not rely on human expert knowledge, and must be able to autonomously adapt to the evolution of all traffic components, be they legitimate or not (this paper obviously focuses on illegitimate ones). For this purpose, we propose to take advantage of unsupervised machine learning techniques that do not require the help of any human expert. In this paper we then present an unsupervised method for classifying and characterizing anomalies and attacks contained in traffic gathered on honeypot networks. This method does not take advantage of any attack signature, learning stage, or labeled traffic, what constitutes a major advance towards autonomous security systems. This approach uses robust clustering techniques, as sub-space clustering, density-based clustering, and evidence accumulation for classifying flow ensembles in traffic classes, and builds easily understandable associated signatures. This method was first proposed in our previous work [START_REF] Mazel | Sub-space clustering, interclustering results association & anomaly correlation for unsupervised network anomaly detection[END_REF] that aims at detecting, classifying and characterizing anomalies in the full network traffic. This method has been adapted and extended in order to cope with the specific honeypot network traffic. Indeed, honeypot network traffic only contains illegitimate traffic (attacks, anomalies, intrusion attempts, ...) what significantly changes the way classification and characterization can be done. This paper mainly deals with taking into account the illegitimate nature of traffic to be analyzed, and designing adapted and optimized classification and characterization algorithms.

These algorithms are illustrated by the analysis of honeypot traffic gathered at the University of Maryland, at the Netflow format [START_REF] Berthier | Nfsight: Netflow-based network awareness tool[END_REF]. The University of Maryland put measurement devices on routers at the border of the honeypots networks. This allows the capture of the traffic exchanged between honeypots and the Internet by monitoring only very few network devices. We nevertheless do not have the communication traces between the honeypots of the same network. The traffic traces we have been working on have been captured on a duration of more than one year.

The rest of the paper is as follows: section II presents related work in the domain of attack and anomalies detection and characterization, covering the full set of approaches from signature to statistical profile based, and finishing with approaches based on data mining and machine learning, especially unsupervised ones. This section also presents supervised and semi-supervised classification methods. Section III presents our previous unsupervised anomaly detection algorithm. This algorithm has already been published in [START_REF] Mazel | Sub-space clustering, interclustering results association & anomaly correlation for unsupervised network anomaly detection[END_REF]. We nevertheless add a presentation of this algorithm in this paper to make it self-contained. Interested readers can find an extended version of the algorithms description in [START_REF] Mazel | Unsupervised network anomaly detection[END_REF]. However this algorithm originally works on full traffic, i.e. normal traffic on Internet links, and the objective is then to detect anomalies and attacks.

Given the nature of the traffic we are considering in this paper -i.e. traffic between honeypots and the Internet -we know that it is completely illegitimate. This is a strong hypothesis that can help going further for characterizing and classifying the different component of this illicit traffic, and estimating the risk in the Internet. This section III then presents the sub-space clustering approach that aims at increasing the robustness of clustering algorithms, limiting its sensitivity, etc. It also presents the recombination mechanisms based on evidence accumulation and inter-clustering associations. As a strong improvement of the previous algorithm, section IV describes how we can show correlation existing between clusters in different subspaces, these correlations exhibiting that such clusters could correspond to the same attack or anomaly. This is a strong help for more accurately and completely identifying and characterizing all attacks contained in the honeypots traffic. Section V then proposes a way to automatically characterize the identified anomalies and attacks. It relies on the clustering and correlation results. Section VI proposes a method for ranking the risk an attack represents. Such ranking is aimed at helping security experts or network operators to prioritize their work, focusing first on the most dangerous ones. Section VII illustrates the results of this characterization algorithm by showing a set of attacks and anomalies that have been identified in the real traffic traces of the University of Maryland, and the characterizing signatures that have been issued. Last, section VIII concludes this paper.

II. RELATED WORK

The problem of network anomaly detection has been extensively studied during the last decade. Most of the approaches analyze statistical variations of traffic volume (e.g. number of packets, bytes or new flows) and/or traffic features (e.g. IP addresses and ports), using either single-link measurements or network-wide data. A non-exhaustive list of standard methods includes the use of signal processing techniques (e.g. ARIMA -Autoregressive Integrated Moving Average -modeling, wavelets-based filtering) on single-link traffic measurements [START_REF] Barford | A signal analysis of network traffic anomalies[END_REF], [START_REF] Brutlag | Aberrant behavior detection in time series for network monitoring[END_REF], PCA (Principal Component Analysis) for networkwide anomaly detection [START_REF] Lakhina | Characterization of networkwide anomalies in traffic flows[END_REF]- [START_REF] Lakhina | Mining anomalies using traffic feature distributions[END_REF], and Sketches applied to IPflows [START_REF] Krishnamurthy | Sketch-based change detection: Methods, evaluation, and applications[END_REF], [START_REF] Li | Detection and identification of network anomalies using sketch subspaces[END_REF].

The simultaneous detection and characterization of traffic anomalies has also received quite a lot of attention in the past, but results are few and present important limitations, either because they rely on some kind of training data and/or anomaly signatures, or because they do not provide meaningful and tractable information to a human network operator, who has to take the final decision about the nature of the detected problem. Authors in [START_REF] Lakhina | Characterization of networkwide anomalies in traffic flows[END_REF] characterize network-wide anomalies in highly aggregated traffic (Origin-Destination flows or OD flows for short), using PCA and the sub-space approach [START_REF] Lakhina | Diagnosing network-wide traffic anomalies[END_REF]. An important limitation of this approach is that the information obtained from OD flow data is too coarse-grained to provide meaningful information to the network operator. Papers like Lakhina et al. [START_REF] Lakhina | Mining anomalies using traffic feature distributions[END_REF] and Biang et al. [START_REF] Li | Detection and identification of network anomalies using sketch subspaces[END_REF] detect and characterize anomalies using finer-grained traffic information, basically applying the same PCA approach to the sample entropy of the empirical distribution of specific traffic features. One clear limitation of these approaches is that the information they provide is not immediately usable and easy-to-understand by the network operator, who may not even be familiar with concepts distant from his tasks such as sample entropy. Besides, the PCA approach is highly sensitive to noise when used for anomaly detection [START_REF] Ringberg | Sensitivity of pca for traffic anomaly detection[END_REF], [START_REF] Casas | Optimal volume anomaly detection and isolation in large-scale ip networks using coarsegrained measurements[END_REF], requiring in practice a fine-tuning and data-dependent calibration step to work. UNADA (Unsupervised Network Anomaly detection Algorithm) [START_REF] Casas | Unada: Unsupervised network anomaly detection using sub-space outliers ranking[END_REF] falls within the unsupervised anomaly detection domain, a novel research area that has drawn quite a lot of interest in the research community, but that still represents a rather immature field. Most work on unsupervised network anomaly detection has been devoted to the IDS field, generally targeting the detection of network intrusions in the very well known KDD'99 dataset. The great majority of the detection schemes proposed in the literature are based on clustering techniques and outliers detection, being [START_REF] Portnoy | Intrusion detection with unlabeled data using clustering[END_REF]- [START_REF] Leung | Unsupervised anomaly detection in network intrusion detection using clustering[END_REF] some examples. The objective of clustering is to partition a set of unlabeled patterns into homogeneous groups of "similar" characteristics, based on some similarity measure. Outliers detection consists in identifying those patterns that do not belong to any of these clusters. In [START_REF] Portnoy | Intrusion detection with unlabeled data using clustering[END_REF], authors use a simple single-linkage hierarchical clustering method to cluster data from the KDD'99 dataset, based on the standard Euclidean distance for interpattern similarity. Eskin et al. [START_REF] Eskin | A geometric framework for unsupervised anomaly detection: Detecting intrusions in unlabeled data[END_REF] reports improved results in the same dataset, using three different clustering algorithms: the Fixed-Width clustering algorithm, an optimized version of the k-NN algorithm, and the one class SVM algorithm. Leung and Leckie [START_REF] Leung | Unsupervised anomaly detection in network intrusion detection using clustering[END_REF] present a combined density-based and grid-based clustering algorithm to improve computational complexity, obtaining similar detection results.

Previous work of our own permits to automatically characterize network traffic anomalies [START_REF] Fernandes | Automated classification of network traffic anomalies[END_REF], but using a-priori welldefined anomaly signatures. Closer to our current work, authors in [START_REF] Silveira | Rca: Pulling anomalies by their root causes[END_REF] present URCA (Unsupervised Root Cause Analysis), a two-steps algorithm to characterize network anomalies in an unsupervised fashion. URCA uses as input the traffic in the anomalous time slots detected by any generic time-slotbased detection algorithm [START_REF] Cormode | What's new: Finding significant differences in network data streams[END_REF]. In the first step, it identifies the anomaly by iteratively removing from the anomalous time slots those flows that seem normal. In the second step, the algorithm uses a hierarchical clustering method to characterize the particular flows identified as anomalous. We identify some serious drawbacks and omissions in URCA: authors claim that the approach is unsupervised, which is not true, simply because it uses previously labeled anomalous events for the characterization. As in previous works, the algorithm uses difficult-to-interpret traffic descriptors for the clustering step (e.g. sample entropy of the distribution of IP addresses, aggregated at different levels), obscuring the comprehension of the network operator. Finally, the algorithm removes those flows that seem normal before the characterization step, which drags possible errors to the clustering step.

Our Unsupervised Anomaly Detection and Characterization algorithm [START_REF] Mazel | Sub-space clustering, interclustering results association & anomaly correlation for unsupervised network anomaly detection[END_REF] presents several advantages w.r.t. current state of the art. First and most important, it works in a completely unsupervised fashion, which means that it can be directly plugged into any monitoring system and start to work from scratch. Secondly, we perform anomaly detection based not only on outliers detection, but also by identifying small-clusters. This is achieved by using different levels of traffic aggregation, both at the source and destination of the traffic; this additionally permits to discover low-intensity and distributed anomalies. Thirdly, we avoid the lack of robustness of general clustering approaches, by combining the notions of Sub-Space Clustering [START_REF] Parsons | Subspace clustering for high dimensional data: a review[END_REF] and multiple Evidence Accumulation [START_REF] Fred | Combining multiple clusterings using evidence accumulation[END_REF]. In particular, our algorithm is immune to general clustering problems like sensitivity to initialization, specification of number of clusters, or structure-masking by irrelevant features. Fourthly, the algorithm performs clustering in low-dimensional feature spaces, using simple traffic descriptors like number of source IP addresses or fraction of SYN packets. This simplifies the characterization of the anomaly, and avoids well-known clustering problems when working with high-dimensional data [START_REF] Jain | Data clustering: 50 years beyond k-means[END_REF]. Our algorithm ranks the multiple evidence of an anomaly detected in different subspaces, combining the most relevant traffic descriptors into a compact and easy-to-interpret signature that characterizes the problem. This permits to reduce the time spent by the network operator to understand the nature of the anomaly. Finally, this algorithm is designed to work in an on-line fashion, analyzing traffic from consecutive time slots in near real time. This is possible even when working with large number of traffic descriptors, because the sub-space clustering and the evidence accumulation algorithms are perfectly adapted for parallelization (see [START_REF] Mazel | Sub-space clustering, interclustering results association & anomaly correlation for unsupervised network anomaly detection[END_REF]).

To the best of our knowledge, there is no paper in the literature on the use of unsupervised classification and characterization algorithms on honeypot traffic, or for intrusion or attacks characterization on big data sets of attack traces. Some very recent work exists, especially by Symantec and the TRIAGE project (Data Analytics Framework for Intelligence Analysis) that aims to use autonomous data mining techniques for the analysis of all the gathered traces of attacks. However, TRIAGE aims at designing visualization techniques for the experts to make decision, while in this paper we propose to autonomously apply countermeasures for cheaper and faster defense.

III. UNSUPERVISED ANOMALY DETECTION

Our anomaly detection works on single-link packet-level traffic captured in consecutive time-slots of fixed length ∆ T . The first analysis stage consists in change detection. At each time-slot, traffic is aggregated in 9 different flow levels l i . These include (from finer to coarser-grained resolution): source IPs (l 1 : IPsrc), destination IPs (l 2 : IPdst), source Network Prefixes (l 3,4,5 : IPsrc/24, /16, /8), destination Network Prefixes (l 6,7,8 : IPdst/24, /16, /8), and traffic per Time Slot (l 9 : tpTS). Time series Z li t are built for basic traffic metrics such as number of bytes, packets, and IP flows per time slot, using the 9 flow resolutions l 1...9 . Analyzing honeypot traffic at multiple aggregation levels permits to detect both single source-destination and distributed attacks of very different intensities.

The unsupervised anomaly detection stage takes as input all the flows in the time slot flagged as anomalous, aggregated according to one of the different levels used in the first stage. An anomaly will generally be detected in different aggregation levels, and there are many ways to select a particular aggregation to use in the unsupervised stage; for the sake of simplicity, we shall skip this issue, and use any of the aggregation levels in which the anomaly was detected. Without loss of generality, let Y = {y 1 ,..,y F } be the set of F flows in the flagged time slot, referred to as patterns in more general terms. Each flow y f ∈ Y is described by a set of A traffic attributes or features. In this paper, we use a list of common traffic attributes. The list includes A =9traffic features: number of source/destination IP addresses and ports, ratio of number of sources to number of destinations, packet rate, ratio of packets to number of destinations, and fraction of ICMP and SYN packets. According to our previous work on signature-based anomaly characterization [START_REF] Fernandes | Automated classification of network traffic anomalies[END_REF], such simple traffic descriptors permit characterization of general traffic anomalies in easy-tointerpret terms. The list is therefore by no means exhaustive, and more features can be easily plugged-in to improve results. Let x f =(x f (1),..,x f (A)) ∈ R A be the corresponding vector of traffic features describing flow y f , and X =(x 1 ; ..; x F ) the complete matrix of features, referred to as the feature space.

The unsupervised detection algorithm is based on clustering techniques applied to X. The objective of clustering is to partition a set of unlabelled patterns into homogeneous groups of similar characteristics, based on some measure of similarity. Table I explains the characteristics of each anomaly in terms of type, distributed nature, aggregation type and netmask used, and impact on traffic features. On one hand, a SYN DDoS which targets one machine from a high number of hosts located in several /24 addresses will constitute a cluster if flows are aggregated in l 3 . In fact, each of these /24 addresses will have traffic attributes values different from the ones of normal traffic: a high number of packet, a single destination and many SYN packets. It is the whole set of these flows that will create a cluster. On the other hand, if flows are aggregated in l 6 , the only destination address will be an outlier characterized by many sources and a high proportion of SYN packets.

Our particular goal is to identify and to isolate the different flows that compose the anomaly flagged in the first stage, both in a robust way. Unfortunately, even if hundreds of clustering algorithms exist [START_REF] Jain | Data clustering: 50 years beyond k-means[END_REF], it is very difficult to find a single one that can handle all types of cluster shapes and sizes, or even decide which algorithm would be the best for our particular problem. Different clustering algorithms produce different partitions of data, and even the same clustering algorithm provides different results when using different initializations and/or different algorithm parameters. This is in fact one of the major drawbacks in current cluster analysis techniques: the lack of robustness.

To avoid such a limitation, we have developed a divide and conquer clustering approach, using the notions of clustering ensemble [START_REF] Strehl | Cluster ensembles -a knowledge reuse framework for combining multiple partitions[END_REF] and multiple clusterings combination. A clustering ensemble P consists of a set of N partitions P n produced for the same data with n =1 ,..,N. Each of these partitions provides a different and independent evidence of data structure, which can be combined to construct a global clustering result for the whole feature space. There are different ways to produce a clustering ensemble. We use Sub-Space Clustering (SSC) [START_REF] Parsons | Subspace clustering for high dimensional data: a review[END_REF] to produce multiple data partitions, applying the same clustering algorithm to N different subspaces U n ⊂ X of the original space.

A. Clustering Ensemble and Sub-Space Clustering

Each of the N sub-spaces U n ⊂ X is obtained by selecting R features from the complete set of A attributes. The number of sub-spaces N hence is equal to R-combinations-obtainedfrom-A. To set the sub-space dimension R, we take a very useful property of monotonicity in clustering sets, known as the downward closure property: "if a collection of points is a cluster in a d-dimensional space, then it is also part of a cluster in any (d -1) projections of this space" [START_REF] Agrawal | Automatic subspace clustering of high dimensional data for data mining applications[END_REF]. This directly implies that, if there exists any evidence of density in X, it will certainly be present in its lowest-dimensional sub-spaces.

Using small values for R provides several advantages: firstly, doing clustering in low-dimensional spaces is more efficient and faster than clustering in bigger dimensions. Secondly, density-based clustering algorithms provide better results in low-dimensional spaces [START_REF] Agrawal | Automatic subspace clustering of high dimensional data for data mining applications[END_REF], because high-dimensional spaces are usually sparse, making it difficult to distinguish between high and low density regions. We shall therefore use R =2 in our SSC algorithm, which gives

N = C A R = A(A -1)/2 partitions.

B. Combining Multiple Partitions

Having produced the N partitions, we now explore different methods to combine these partitions in order to build a single partition where anomalous flows are easily distinguishable from normal-operation traffic: the classical Evidence Accumulation (EA) and the new Inter-Clustering Result Association (ICRA) method.

1) Combining Multiple Partitions using Evidence Accumulation:

A possible answer is provided in [START_REF] Fred | Combining multiple clusterings using evidence accumulation[END_REF], where authors introduced the idea of multiple-clusterings Evidence Accumulation (EA). By simple definition of what it is, an anomaly may consist of either outliers or small-size clusters, depending on the aggregation level of flows in Y (cf table I). EA then uses the cluster ensemble P to build two inter-pattern similarity measures between the flows in Y. These similarity measures are stored in two elements: a similarity matrix S to detect small clusters and a vector D used to rank outliers. S(p, q) represents the similarity between flows p and q. This value increases when the flows p and q are in the same cluster many times and when the size of this cluster is small. These two parameters allows the algorithm to target small clusters. D(o) represents the abnormality of the outlier o. This value increases when the outlier has been classified as such several times and when the separation between the outlier and the normal traffic is important. As we are only interested in finding the smallestsize clusters and the most dissimilar outliers, the detection consists in finding the flows with the biggest similarity in S and the biggest dissimilarity in D. Any clustering algorithm can then be applied on the matrix S values to obtain a final partition of X that isolates small-size clusters of close similarity values. A variable detection threshold over the values in S is also able to detect small-size cluster. Concerning dissimilar outliers, they can be isolated though a threshold applied on the values in D.

2) Combining Multiple Partitions using Inter-Clustering result Association: However, by reasoning over the similarities between patterns (here flows), EA introduces several potential errors. Let us consider two pattern sets P i and P j , if the cardinality of these pattern sets is close and if they are present in a similar number of sub-spaces, then EA will produce a very close (potentially the same) similarity value for both flow sets. They will then likely be falsely considered as belonging to the same cluster. This possibility has to be considered very seriously as it can induce a huge error: different anomalies will be merged together and will then likely be wrongly identified and characterized. Another source of potential error when using a clustering algorithm over S values is the algorithm sensitivity to wrong parameters. Furthermore, the use of a threshold over S and/or D can decrease the system performance in case of a wrong value used.

In order to avoid the previously exposed sources of error, we introduce a new way of combining clustering results obtained from sub-spaces: Inter-Clustering Results Association. The idea is to address the problem in terms of cluster of flows and outlier of flow similarity instead of pattern (or flow) similarity. Hence, we shift the similarity measure from the patterns to the clustering results. The problem can then be split in two sub-problems: correlate clusters through Inter-CLuster Association (ICLA), and correlate outlier through Inter-Outlier Association (IOA).

In each case, a graph is used to express similarity between either clusters or outliers. Each vertex is a cluster/outlier from any sub-space U n and each edge represents the fact that two connected vertices are similar. The underlying idea is straightforward: identify clusters or outliers present in different sub-spaces that contain the same flows. To do so, we first define a cluster similarity measure called CS between two clusters C r and C s : CS(C r ,C s )= card(Cr∩Cs) max(card(Cr),card(Cs) , card being the function that associates a pattern set with its cardinality, and C r ∩ C s the intersection of C r and C s . Each edge in the cluster similarity graph between two C r and C s means CS(C r ,C s ) > 0.9, being this an empirically chosen value. The value 0.9 guarantees that the vast majority of patterns are located in both clusters with a small margin of error. IOA uses an outlier similarity graph built by linking every outlier to every other outlier that contains the same pattern. Once these graphs are built, we need to find cluster sets where every cluster contains the same flows. In terms of vertices, we need to find vertex sets where every vertex is linked to every other vertex. In graph theory, such vertex set is called a clique. The clique search problem is a NP-hard problem. Most existing solutions use exhaustive search inside the vertex set which is too slow for our application. We then make the hypothesis that a vertex can only be part of a single clique. A greedy algorithm is then used to build each clique. Anomalous flow set are finally identified as the intersection of all the flow sets present in the clusters or outliers within each clique.

IV. CORRELATING ANOMALOUS TRAFFIC CLASSES

A. Address related correlation

Thanks to previous algorithm, we can detect several classes of illegitimate traffic, but these classes can appear in different aggregation levels. We know that such classes present at different levels can be related to each other. Indeed, two traffic classes in two different aggregation levels are related for example when their flows come from the same sources and go towards the same destinations. For a better characterization of traffic classes, it is then important to link classes corresponding to the same single anomaly. Correlating illegitimate traffic classes is a solution for that purpose; it is able to determine the similarity between classes at each aggregation level, for possibly grouping them if they belong to the same anomaly.

For estimating the similarity between two illegitimate traffic classes on two different aggregation levels, we use a comparison function. It relies on IP addresses comparison [START_REF] Mazel | Sub-space clustering, interclustering results association & anomaly correlation for unsupervised network anomaly detection[END_REF]. The comparison method then uses the IP source and destination addresses of the different traffic classes a 1 and a 2 . It then compares source IP addresses with each other, and destination IP addresses with each other, using function [START_REF] Mazel | Sub-space clustering, interclustering results association & anomaly correlation for unsupervised network anomaly detection[END_REF].

Sim @ (@ 1 , @ 2 )= |@ 1 ∩ @ 2 | max(|@ 1 |, |@ 2 |) (1)
where @ n is a set of IP addresses, and |@ n | the number of addresses in this set.

It exists a similarity between two illegitimate traffic classes if equation ( 2) is true. tT SAddrSims is a threshold to be defined.

Sim src (@ 1 , @ 2 ) > tT SAddrSims ∧ Sim dest (@ 1 , @ 2 ) > tT SAddrSims (2)
This method as defined in our previous work [START_REF] Mazel | Sub-space clustering, interclustering results association & anomaly correlation for unsupervised network anomaly detection[END_REF] has some lacks, and especially because it does not consider time. It appears when testing this method on the Maryland data that some flows of the same class are separated in time by several months. It is then clear that it is not satisfactory to consider only IP addresses. For example, illegitimate traffic classes detailed in tables II and III are completely different from a behavior point of view. For the first class, the number of sent packets in each flow is around 54, whereas it is around 4 in the second class. In addition, the time difference between these sendings is around two months. The only similarity between these two classes is related to their source IP addresses (saddr), and destination IP addresses (daddr). In this case the similarity value between the two classes with the previous function is 66%, what is significantly high, and would indicate a strong link between them. Given the time at which they happened it is certainly not true, and malware infecting the machines are certainly not the same. As a consequence, we added temporal features to the correlation function.

Let us consider two illegitimate traffic classes a 1 and a 2 . T 1 and T 2 are the set of time intervals in which traces of a 1 and a 2 appear. @ n 1 and @ n 2 are the set of source and destination IP addresses of classes a 1 and a 2 for the interval t n ∈ T 1 . The new similarity function is then defined as:

SimT ime @ = X t∈T1 Sim @ (@ t 1 , @ t 2 ) max(|@ t 1 |, |@ t 2 |) (3)
However, if it exists time intervals between two classes, the similarity function equals a value very close from zero, whereas these two classes can be highly similar for other features. We then still need to keep in the computing of the similarity function the IP addresses of the anomaly classes a 1 and a 2 . But the addresses must be computed intependentely from any temporal feature.

SimGlobal @ = Sim @ (@ 1 , @ 2 ) max(|@ 1 |, |@ 2 |) (4)
We then obtain the simlarity function defined by equation [START_REF] Brutlag | Aberrant behavior detection in time series for network monitoring[END_REF].

SimAnomalies(a1,a2) = (SimT ime Src (a 1 ,a 2 ) > λ 1 ) ∧ (SimT ime Dest (a 1 ,a 2 ) > λ 1 ) ∧ (SimGlobal Src (a 1 ,a 2 ) > λ 2 ) ∧ (SimGlobal Dest (a 1 ,a 2 ) > λ 2 ) (5)
The result of this function is a boolean value which is true if a 1 and a 2 are similar.

B. Time related correlation

Let's take again the examples of classes detailed in tables II and III. It is clear that there is no link between these two classes because they happened at very different times. Let's apply equation ( 5). Let's consider IP addresses of the first class on February 1st at 2h40 am. Source IP addresses are 192.168.0.1 and 192.168.0.2. At that time, the second class is empty (it happened on April 10th to 13th). The similarity is then correctly estimated as false. By continuing with other times and the same method, we always obtain a false value for the similarity. The new similarity function then correctly does not find any similarity between these two illegitimate traffic classes.

Nevertheless, it is required with this new method to fix correct threshold values. If thresholds are not well selected returned values could be erroneous. With a small λ 1 value, each time a small temporal similarity will appear between two illegitimate classes, it will be the global similarity value that will determine the final similarity result between the two classes. If λ 1 is high, the two classes will be considered as similar if they sent packet almost at the same time, from the same source and to the same destinations. This is a strong constraint and forbids any time difference in the sendings. The global similarity will then be of less importance. Based on our experience with the traces of University or Maryland, we empirically recommend to select for λ 1 a low threshold value, between 10% and 30%, and for λ 2 a threshold value greater than 25%.

Illegitimate traffic classes that appear as different after the sub-space clustering phase can then be grouped, as it is shown that they correspond to the same anomaly.

V. AUTOMATIC CHARACTERIZATION OF ANOMALIES

At this stage, the global traffic from and towards honeypots has been decomposed into traffic classes that exhibit different behaviors. Classification techniques of the Internet traffic have now to identify the type of each of these traffic classes, the generating application, or the type of attack or anomaly, each of these clusters is related to. Most advanced classification techniques, i.e. semi-supervised techniques, take advantage of signatures that specifically identify one of the possible traffic families, applications, or attacks. These signatures come either from a previous knowledge or expertise in this domain, either on a training stage on a known traffic, whose application components have already been labeled. This is obviously a strong requirement for classification purposes, but it remains a severe limit for designing a fully autonomous method. However, in this work, at the opposite of what has been already done in the research area of autonomous traffic classification, this constraint does not exist because of the traffic nature that we have to analyze: illegitimate traffic. It is thus not needed to perfectly identify the anomaly kind (or attack), or to name it. Indeed, all traffic classes that have been isolated by the sub-space clustering algorithm, and evidence accumulation are anomalies or attacks. Therefore, the computing to be performed for each traffic anomaly or attack, in fine, consists in discarding them, after having identified any necessary feature required for instance for estimating the risk they represent in the Internet.

We then propose for this purpose to automatically generate the rules characterizing the anomaly classes. Based on these rules, it is easy to understand the anomalies characteristics and to infer the countermeasures to be performed (and of course, this can be done by a computer process).

At this stage, the sub-space clustering / evidence accumulation / correlation algorithm has identified several correlated anomalies containing a set of traffic flows in Y far out the rest of the traffic. The following task is to produce the appropriate filtering rules to correctly isolate and characterize each of these anomalies.

In order to produce filtering rules, the algorithm selects those sub-spaces U n where the separation between the considered anomalous flows and the rest of the traffic is the biggest. We define two different classes of filtering rule: absolute rules FR A (Y) and relative rules FR R (Y). Absolute rules do not depend on the separation between flows, and correspond to the presence of dominant features in the considered flows. An absolute rule for a certain feature j characterizing a certain flow set Y g has the form

FR A (Y g ,a)={∀y f ∈ Y g ⊂ Y : x f (a)==λ}.
For example, in the case of an ICMP flooding attack, the vast majority of the associated flows use only ICMP packets, hence the absolute filtering rule {nICMP/nPkts == 1} makes sense. On the contrary, relative filtering rules depend on the relative separation between anomalous and normaloperation flows. Basically, if the anomalous flows are well separated from the normal cluster in a certain partition P n , then the features of the corresponding sub-space U n are good candidates to define a relative filtering rule. A relative rule has the form

FR R (Y g ,a)={∀y f ∈ Y g ⊂ Y : x f (a) < λ ∨ x f (a) > λ}.
We shall also define a covering relation between filtering rules: we say that rule

f 1 covers rule f 2 ⇔ f 2 (Y) ⊂ f 1 (Y).
If two or more rules overlap (i.e., they are associated to the same feature), the algorithm keeps the one that covers the rest.

In order to construct a compact signature of the anomaly, we have to devise a procedure to select the most discriminant filtering rules. Absolute rules are important, because they define inherent characteristics of the anomaly. As regards relative rules, their relevance is directly tied to the degree of separation between anomalous and normal flows. In the case of outliers, we select the K features for which the Mahalanobis distance to the normal-operation traffic is among the top-K biggest distances. In the case of small-size clusters, we rank the relatives rules according to the degree of separation to the normal anomaly using the well-known Fisher Score (FS) which uses the variance in each cluster (normal and anomalous). To finally construct the signature, the absolute rules and the top-K relative rules are combined into a single inclusive predicate, using the covering relation in case of overlapping rules. On Figure 1, there exists a relative rule between cluster 1 and cluster 4. For generating it, we draw the median between clusters 1 and 4. Based on this median value, we can create the relative rule {nDif f DestAddr < 9}.

VI. RISK BASED ANOMALIES RANKING

Estimating the risk related to an anomaly can help a security expert, in the general case, selecting the traffic classes that require to be computed in priority. In our specific case where all traffic classes identified are illegitimate, it can help to distinguish for instance between anomalies that are real attacks from scanning which in general just serves for preparing a possible future attack. In that case, the risk ranking can help superposing relative filtering rules, targeting in priority the most risky one.

The way our risk ranking is done depends on the amount of communication; the more communication in an anomaly or attack, the more risky. In addition, anomalies that appear in many sub-spaces are considered as more dangerous than the ones appearing in a single or on very few sub-spaces. Three features are considered for estimating the risk related to an anomaly. The first one is the number of packets of the anomaly, because, if there are many packets exchanged between the attacking and the victim machines, there exist potentially machines infected by a virus, or a flooding attack attempt. The second feature is the amount of bytes exchanged in the anomaly, because, even if the number of exchanged packets is reduced, a large amount of bytes could have been exchanged. This may correspond to the download of information from a victim machine by an attacker, or the upload of viruses on several victim machines from the attacker machine. Last, the third feature is the communication duration between the attacker and a target machine. The longer the communication duration, the more probable an attacker performing a download, or having an open shell on the victim machine.

The formula used for the risk estimation is: risk = C * (log(nP kts)+log(nBytes)+log(duration + 1)) (6) where C is the number of sub-spaces in which the anomaly appears, nP kts is the number of exchanged packets in the anomaly, nBytes is the number of exchanged bytes, and duration is the duration of the anomaly. +1 appears in log(duration+1) for avoiding some errors as duration ∈ R + whereas (nP kts, nBytes) ∈ R * + . The communication duration can be assimilated as zero because it can be so small that measurement devices can measure it as zero. On the other side, the number of packets or bytes are necessarily greater or equal to 1.

VII. EXPERIMENTAL EVA L UAT I O N I N REAL TRAFFIC

We run the algorithm described in this paper on the honeypot traffic traces gathered at the University of Maryland. For obvious privacy reasons, as well as space limit, we will not present the complete set of attacks evidenced. But we will show on an example how the algorithm behaves, and how it succeed in classifying attacks and anomalies, how it builds the anomaly classes characteristics based on sub-space clustering, evidence accumulation, and anomaly correlation. It also presents the filtering rules that have been autonomously generated and that can serve for automatically configuring security devices as filtering functions of routers, or firewalls. This section starts by showing the complete behavior of the algorithm applied on a single anomaly (designated as anomaly [224]). Figure 2 shows 5 sub-spaces in which anomalies mentioned in this section appear (all sub-spaces in which the anomaly appears are not depicted on figure 2 for space limit reason). On each of these sub-spaces, clusters clearly appear, each corresponding to different anomalies. Figure ?? shows the cliques that have been built for correlating the clusters found in different sub-spaces, and then linking the ones that correspond to the same anomaly or attack. This paper presents an unsupervised algorithm for classifying illicit traffic. This algorithm has several advantages compared to previous work: (i) it works in a completely unsupervised manner, what makes it able to work on top of any monitoring system, and directly usable, without preliminary configuration or knowledge. (ii) It combines robust clustering techniques to avoid classical issues of clustering algorithms, e.g. sensitivity to initial configuration, the required a priori indication of teh number of clusters to be identified, or the sensitivity of results when using less pertinent features. (iii) It automatically builds simple and small signatures fully characterizing attacks; theses signature can then be used in a filtering security device. (iv) It is designed to run in real time by making possible to take advantage of the parallelism of our clustering approach. This algorithm thus opens new perspectives for performing a risk analysis in the Internet -taking advantage of honeypot traffic -and automatically configuring related filtering rules on routers, switches, or firewalls.

Anomaly [224] appears in

Fig. 1 .

 1 Fig. 1. The different filtering rules for sub-spaces (nSyn/nPkts, nDiffDestAddr) What follows gives a real example of the generation of filtering rules based on an anomaly characterization. It is depicted on Figure 1. Clusters 2 and 3 have their value nSyn/nPkts always equal to 1. It then exists for them an absolute rule {nSyn/nP kts == 1}.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Sub-spaces in which anomalies [44], [224], and [327] appear. These sub-spaces correspond to different IP address aggregation levels and different temporal granularities

TABLE I .

 I FEATURE USED FOR THE DETECTION OF DOS, DDOS, NETWORK/PORT SCANS, AND SPREADING WORMS.A NOMALIES OF DISTRIBUTED NATURE 1-TO-N OR N-TO-1 INVOLVE SEVERAL /24 (SOURCE OR DESTINATIONS) ADDRESSES CONTAINED IN A SINGLE /16 ADDRESS.

	Anomaly	Distributed nature	Aggregation type	Clustering result	Impact on traffic features
	DoS (ICMP ∨ SYN)	1-to-1	IPsrc/ * IPdst/ *	Outlier Outlier	nSrcs = nDsts = 1, nPkts/sec > λ1, avgPktsSize < λ2, (nICMP/nPkts > λ3 ∨nSYN/nPkts > λ4).
	DDoS (ICMP ∨ SYN) to several @IP/24	N-to-1	IPsrc/24 (l3) IPsrc/16 (l4) IPdst/ *	Cluster Outlier Outlier	nDsts = 1, nSrcs > α1, nPkts/sec > α2, avgPktsSize < α3, (nICMP/nPkts > α4 ∨ nSYN/nPkts > α5).
	Port scan	1-to-1	IPsrc/ * IPdst/ *	Outlier Outlier	nSrcs = nDsts = 1, nDstPorts > β1, avgPktsSize < β2, nSYN/nPkts > β3.
	Network scan to several @IP/24	1-to-1	IPsrc/ * IPdst/24 (l6) IPdst/16 (l7)	Outlier Cluster Outlier	nSrcs = 1, nDsts > δ1, nDstPorts > δ2, avgPktsSize < δ3, nSYN/nPkts > δ4.
	Spreading worms to several @IP/24	1-to-N	IPsrc/ * IPdst/24 (l6) IPdst/16 (l7)	Outlier Cluster Outlier	nSrcs = 1, nDsts > η1, nDstPorts < η2, avgPktsSize < η3, nSYN/nPkts > η4.

ACKNOWLEDGEMENTS

The author sincerely thanks Michel Cukier and Bertrand Sobesto for providing the traffic traces gathered on the honeypots of the University of Maryland. The author thanks Johan Mazel and Pedro Casas who have been first involved in the research work on the sub-space clustering algorithm applied to the full Internet traffic. The author also thanks Richard Turc who started this work on the analysis of the honeypot traffic during his master internship at LAAS. This work is supported by the ONTIC project, funded by the European commission under grant FP7-ICT-2013-11/619633.