Franck Tirode

Karine Laud-Duval

Alexandre Prieur

Bruno Delorme

Pierre Charbord

INSERMU830 Olivier Delattre
email: olivier.delattre@curie.fr

Mesenchymal stem cell features of Ewing tumors

The cellular origin of Ewing tumor (ET), a tumor of bone or soft tissues characterized by specific fusions between EWS and ETS genes, is highly debated. Through gene expression analysis comparing ETs with a variety of normal tissues, we show that the profiles of different EWS-FLI1-silenced Ewing cell lines converge toward that of mesenchymal stem cells (MSC). Moreover, upon EWS-FLI1 silencing, two different Ewing cell lines can differentiate along the adipogenic lineage when incubated in appropriate differentiation cocktails. In addition, Ewing cells can also differentiate along the osteogenic lineage upon long-term inhibition of EWS-FLI1. These in silico and experimental data strongly suggest that the inhibition of EWS-FLI1 may allow Ewing cells to recover the phenotype of their MSC progenitor.

Significance

The precise knowledge of the normal counterpart of cancer cells is dramatically lacking for Ewing sarcoma, one of the most frequent primary bone tumors. Since its initial description more than 80 years ago, endothelial, epithelial, neural and mesenchymal origins have been hypothesized. We have developed original cell models to regulate EWS-FLI1 expression. In silico analyses and in vitro differentiation experiments provide strong evidences that Ewing tumors originate from a mesenchymal precursor that may have both connective and vascular differentiation potentials. This origin helps to clarify the issue of the predominant localization of ET in bone or soft tissue, two tissues that constitute major sources of MSC.

Introduction

Ewing tumor (ET), the second most frequent bone tumor in adolescent and young adult, harbors characteristic translocations which fuse the 5' portion of the EWS gene with the 3' region, encoding DNA binding domain, of one of five ETS family genes [START_REF] Arvand | Biology of EWS/ETS fusions in Ewing's family tumors[END_REF][START_REF] Janknecht | EWS-ETS oncoproteins: the linchpins of Ewing tumors[END_REF]. The specificity of these translocations for ET has now been extensively documented. The most frequent fusion occurs with the FLI1 gene as a consequence of the t(11;22)(q24;q12) translocation [START_REF] Delattre | Gene fusion with an ETS DNA-binding domain caused by chromosome translocation in human tumours[END_REF]. The resulting EWS-FLI1 chimeric protein behaves as an aberrant transcription regulator that is thought to induce transformation through induction or repression of specific target genes. Accordingly, a number of genes controlling cell growth, signal transduction or differentiation have been shown to be regulated by EWS-FLI1 (reviewed in [START_REF] Janknecht | EWS-ETS oncoproteins: the linchpins of Ewing tumors[END_REF].

Histologically, ET appears as uniformly undifferentiated small round cells containing a vesicular nuclei and a small cytoplasm within a sparse intercellular stroma. These characteristics are poorly indicative of the cell lineage that gives rise to ET. Since the initial proposal of an endothelial origin by James Ewing in 1921 [START_REF] Ewing | Diffuse endothelioma of bone[END_REF], the issue of the ET progenitor is highly debated. A neural crest origin is hypothesized by the observation that EWS-ETS fusions characterize a spectrum of tumors with undifferentiated ETs at one end and differentiated peripheral primitive neuroectodermal tumors, expressing various neural markers, at the other end. This hypothesis is further strengthened by the observation that ET derived cell lines grown in vitro can undergo neural differentiation upon treatment with various differentiationinducing agents [START_REF] Cavazzana | Experimental evidence for a neural origin of Ewing's sarcoma of bone[END_REF][START_REF] Noguera | Patterns of differentiation in extraosseous Ewing's sarcoma cells. An in vitro study[END_REF]. A recent study suggests that EWS-FLI1 expression itself may be responsible, at least in part, for the observed primitive neuroectodermal phenotype of ET through the regulation of neural specific genes [START_REF] Hu-Lieskovan | EWS-FLI1 fusion protein up-regulates critical genes in neural crest development and is responsible for the observed phenotype of Ewing's family of tumors[END_REF]. This study indicates that the ultimate phenotype of ET might be determined not only by the cell lineage from which it derives but also by intrinsic effects of EWS-FLI1. EWS-FLI1 has also been shown to inhibit neural differentiation of neuroblastoma cells therefore leading to the proposal that EWS-FLI1 may shift the sympathetic differentiation program of neuroblastoma to the parasympathetic one of ET [START_REF] Rorie | The Ews/Fli-1 fusion gene switches the differentiation program of neuroblastomas to Ewing sarcoma/peripheral primitive neuroectodermal tumors[END_REF].

Other experiments indicate that EWS-FLI1 inhibits mesenchymal differentiation. In particular, enforced expression of EWS-FLI1 inhibits adipogenic and osteogenic differentiation in murine marrow stromal cells [START_REF] Torchia | Ewing tumor fusion proteins block the differentiation of pluripotent marrow stromal cells[END_REF]. It also impairs myogenic differentiation in C2C12 myoblasts [START_REF] Eliazer | Alteration of mesodermal cell differentiation by EWS/FLI-1, the oncogene implicated in Ewing's sarcoma[END_REF]. The hypothesis of a mesenchymal origin of ET is further supported by two recent reports indicating that EWS-FLI1 can transform primary bone marrow stromal cells and induce the formation of tumors resembling Ewing sarcomas in SCID mice [START_REF] Castillero-Trejo | Expression of the EWS/FLI-1 oncogene in murine primary bone-derived cells Results in EWS/FLI-1-dependent, ewing sarcoma-like tumors[END_REF][START_REF] Riggi | Development of Ewing's sarcoma from primary bone marrow-derived mesenchymal progenitor cells[END_REF]. This indicates that, unlike other cell types which cannot be transformed by EWS-FLI1 alone and necessitate additional events, mesenchymal stem cells display permissiveness for EWS-FLI1-mediated cell transformation and may therefore constitute an appropriate background for ET development.

Finally, together with the observation of occasional cytokeratin staining, the study of tight junction structural proteins has suggested that ETs exhibit partial epithelial differentiation [START_REF] Schuetz | Intercellular junctions in Ewing sarcoma/primitive neuroectodermal tumor: additional evidence of epithelial differentiation[END_REF].

Presently, all attempts to identify the ET cell of origin have relied on the overexpression of EWS-FLI1 in various heterologous backgrounds. In this study, we have tested the hypothesis that Ewing cells, inhibited for EWS-FLI1 by a specific small hairpin RNA (shRNA), may recover characteristics of their primary origin and responsiveness to differentiation programs.

Experimental Procedures

Inhibition of EWS-FLI1

A673 and SK-N-MC were obtained from ATCC. EW24 derived from a bone tumour was kindly provided by Pr Gilbert Lenoir [START_REF] Zucman | Combinatorial generation of variable fusion proteins in the Ewing family of tumours[END_REF]. All three cell lines exhibit a type 1 fusion transcript. For siRNA mediated silencing, A673 Ewing cells expressing a type 1 fusion transcript were transfected mainly as previously described [START_REF] Prieur | EWS/FLI-1 silencing and gene profiling of Ewing cells reveal downstream oncogenic pathways and a crucial role for repression of insulin-like growth factor binding protein 3[END_REF], except that transfection medium was removed 6 hours after transfection instead of 24 hours and replaced by fresh medium changed every two days. For doxycycline inducible silencing, the oligonucleotides of the EWS-FLI1 short cassette from pTER/shEF-1 vector was inserted in the pTRIP lentiviral expression vector. Packaging and titration using GFP as a reporter were performed according to the tronolab website (http://tronolab.epfl.ch/).

Isolation of human bone marrow stromal cells (BMSC).

Human primary BMSCs were obtained after informed consent from bone marrow aspirates (iliac crest) of patients undergoing hip replacement surgery. Mononuclear cells were plated at the density of 50 000 cells/cm 2 and cultured in α-MEM without nucleotides (Gibco BRL, Invitrogen, Cergy-Pontoise, France), supplemented with 10% FCS (Hyclone), L-glutamine 2 mM (Gibco BRL) and 1ng/mL bFGF (AbCys).

After two weeks, confluent cells were detached with trypsin/EDTA and replated at 1000 cells/cm 2 (passage 1, P1). Cells were phenotyped as BMSC by flow cytometry (negative for CD34, CD31 and CD45 and positive for CD73, CD90, CD105 and CD44). The three-lineages (adipocyte, osteocyte, and chondrocyte) differentiation potential of these cells was controlled.

Isolation of RNAs and microarray analyses.

Total RNAs were isolated using the Trizol Reagent (Invitrogen). RNAs were isolated 9 days following transfection of A673 Ewing cells with the specific EWS-FLI1 siRNA (siEF1) or with the control siRNA (siCT) and 3 days after infection of EW24 and SK-N-MC with lentiviruses encoding EWS-FLI1-specific or control shRNA. RNAs from 27

ET samples and from freshly isolated (P1) BMSCs were also used for microarray analyses. Experimental procedures for GeneChip microarray were performed according to the Affymetrix GeneChip Expression Analysis Technical Manual (Affymetrix, Santa Clara, CA) using HG-U133A arrays. CEL files from E-AFMX-5 [START_REF] Su | A gene atlas of the mouse and human proteinencoding transcriptomes[END_REF] and from E-MEXP-167 and E-MEXP-168 [START_REF] Boquest | Isolation and transcription profiling of purified uncultured human stromal stem cells: alteration of gene expression after in vitro cell culture[END_REF] datasets were downloaded from the EBI's ArrayExpress repository [START_REF] Parkinson | ArrayExpress--a public repository for microarray gene expression data at the EBI[END_REF].

Additional CEL files for EWS-FLI1-silenced A673 cells were kindly provided by Steve Lessnick [START_REF] Smith | Expression profiling of EWS/FLI identifies NKX2.2 as a critical target gene in Ewing's sarcoma[END_REF]. All microarray data were simultaneously normalized using the gcrma package version 1.1.4 in the R environment (R Development Core Team, 2005). Principal component analyses were performed using the PCA tool of MeV TM4 software [START_REF] Saeed | TM4: a free, open-source system for microarray data management and analysis[END_REF]. SAM analyses were performed with the samr package version 1.22 [START_REF] Tusher | Significance analysis of microarrays applied to the ionizing radiation response[END_REF]. Microarray data are available at Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) public database with the accession number GSE7007.

Real-time quantitative reverse transcription-PCR.

cDNAs were synthesized from 1µg of RNA using the GeneAmp RNA PCR core Kit (Applied Biosystem, Courtaboeuf, France). Quantitative PCR analyses were performed using TaqMan Assays-on-demand Gene expression reagents (Applied Biosystem) with qPCR Mastermix Plus without UNG (Eurogentec, Belgium).

Reactions were run on an ABI/PRISM 7500 (Applied Biosystem) and analyzed using the 7500 system SDS software. The following Assays-on-demand were used: FABP4 RPLP0 expression, which does not vary upon shRNA induction or differentiation conditions, was used for normalization. EWS-FLI1 expression was quantified by SYBR green (Applied Biosystem) using the following primers: forward 5'GAGGCCAGAATTCATGTTATTGC3'; reverse 5'GCCAAGCTCCAAGTCAATATAGC3'.

Cell cycle analyses by flow cytometry

Cell cycle analyses were performed by BrdU incorporation as previously described [START_REF] Dauphinot | Analysis of the expression of cell cycle regulators in Ewing cell lines: EWS-FLI-1 modulates p57KIP2and c-Myc expression[END_REF]. Samples were subjected to FACS analysis (FACScalibur, BD Biosciences, San Diego, USA) and the data were processed using CELL Quest software (BD Biosciences). For cell surface markers detection, 10 5 cells were incubated 30 min at 4°C in 100 µL of PBS with phycoerythrin (PE)-conjugated mAb.

Cells were then washed twice and then proceeded immediately for FACS analysis.

The PE-conjugated antibodies from BD Biosciences used were: IgG1-PE (MOPC-2), CD29-PE (MAR4), CD31-PE (WM59), CD44-PE (G44-26), CD45-PE (HI30), CD54-PE (HA58), CD59-PE (H19), CD90-PE (5E10) and CD166-PE (3A6). CD105-PE (SN6) was from Invitrogen.

Differentiation assays.

For adipogenic differentiation, cells were plated and grown for 2 days at 37°C, 5%

CO 2 in DMEM or RPMI containing 10% FCS, 0.5mM IsoButylMethylXanthine (Sigma), 1µM dexamethasone (Sigma) and 10µg/mL Insulin (Sigma). Medium was then replaced by DMEM or RPMI, 10% FCS, with 10µg/mL Insulin and cells were grown for six additional days with medium changes every two days, then collected for quantitative-RT-PCR or fixed in 4% paraformaldehyde (Sigma) prior to staining with 0.3% w/v Oil-Red-O (Sigma) in 60% isopropanol. For osteogenic differentiation, cells were grown in DMEM, 10% serum, 0.1µM dexamethasone, 0.15mM L-ascorbic acid (Sigma), 2mM β-glycerophosphate (Sigma) and 1mM NaH 2 PO 4 for 21 days.

Paraformaldehyde fixed calcium depositions were stained by the Von-Kossa protocol including incubation with 5% silver nitrate (Sigma) for 30 min, extensive washing with PBS, exposition for 20 min to UV and incubation in 5% sodium thiosulfate for 5 min.

For chondrogenic differentiation, cells were grown in pellet as previously described [START_REF] Johnstone | In vitro chondrogenesis of bone marrow-derived mesenchymal progenitor cells[END_REF]. Briefly, 4x10 5 cells centrifuged 5 min at 1200 rpm in a 15-ml polypropylene tube were grown in pellet in 500µl of DMEM without serum but containing 0.1µM dexamethasone, 0.15mM L-ascorbic acid, 0.35 mM proline, 1mM sodium pyruvate, 1x ITS (Roche, Meylan, France), 10ng/mL TGFβ3 (R&D system, Lille, France), and 10ng/mL IGF1 (Abcys) for 21 days.

Results

Convergence of EWS-FLI1-silenced Ewing cells toward mesenchymal stem cells in silico

To gain insight into the origin of Ewing tumors (ET), we compared the expression profile of 27 ETs with that of various normal human tissues [START_REF] Su | A gene atlas of the mouse and human proteinencoding transcriptomes[END_REF], human Mesenchymal Stem Cells (MSCs) [START_REF] Boquest | Isolation and transcription profiling of purified uncultured human stromal stem cells: alteration of gene expression after in vitro cell culture[END_REF] expression by RNA interference, were also included in this in-silico analysis. These silencing experiments were initially performed with an EWS-FLI1-specific siRNA for the A673 cell line which can be efficiently transfected (Prieur et al.2004). Given the low efficiency of such transfection in other cell lines, we developed lentiviral vectors to express a shRNA corresponding to the junction between EWS and FLI1 genes.

This lentivirus was used to infect EW24 and SK-N-MC cells. In addition, data recently published for EWS-FLI1-silenced A673 cells [START_REF] Smith | Expression profiling of EWS/FLI identifies NKX2.2 as a critical target gene in Ewing's sarcoma[END_REF] were also included in this analysis (a complete description of the whole dataset is provided in Table S1).

All experiments, performed on Affymetrix HG-U133A chips, were simultaneously normalized using gcrma. Hierarchical clustering and principal component analyses (PCA) using all probe sets showed that ETs clustered, together EWS-FLI1-silenced (efRNAi) and control (ctrlRNAi) samples, apart from all other tissues. This indicated that ET form a homogeneous group not evidently related to any normal tissue (Figure 1A). Based on the hypothesis that EWS-FLI1 may modulate differentiation genes, we performed a new PCA analysis with genes strongly and significantly modulated by EWS-FLI1 as identified by SAM analysis comparing efRNAi and ctrlRNAi samples.

As shown in figure 1B, ETs were still clustered with ctrlRNAi apart from other tissues.

In contrast, efRNAi samples shifted away from the ET cluster to join the mesenchymal group containing undifferentiated BMSCs or MSCs but also the more differentiated adipocytes and myocytes. Hierarchical clustering analysis using the same probesets also enabled to visualize the vicinity of efRNAi with tissues derived from the mesenchymal lineage (Figure 1C and Figure S1). These in silico analyses provided a first clue that EWS-FLI1 may regulate a subset of genes involved in mesenchymal commitment or differentiation.

We more precisely investigated similarities between genes modulated by EWS-FLI1 in Ewing cell lines and genes differentially expressed between MSCs and Ewing tumors. SAM analyses enabled the identification of 3238 probe sets ( 2430genes) differently expressed between MSCs and ETs (absolute fold change>2; qvalues < 10 -2 )(Table S2). Comparisons of MSCs/ETs genesets with the 289 genes identified between efRNAi and ctrlRNAi using the same criteria demonstrated a very significant overlap of 110 genes (p< 10 -16 , exact binomial test) (Table S2). Most of these genes behaved in a consistent manner between both comparisons, genes downregulated by EWS-FLI1 being overexpressed in MSCs as compared to ETs and genes upregulated by EWS-FLI1 being underexpressed in MSCs (Figure 1D). We next investigated gene sets specifically modulated by EWS-FLI1 in each individual cell lines. As expected, for each cell line, the overlap between MSCs/ETs-differential genes and EWS-FLI1-regulated genes was highly significant although it was larger for A673. In addition to a common core set of genes regulated by EWS-FLI1 in all three cell lines, this analysis pointed out cell line-specific genesets (Figure 1E and Table S2). Among genes common to the MSC/ET comparison and regulated by EWS-FLI1 in at least 2 of the Ewing cell lines (grey area in Figure 1E), various signaling molecules (DKK1, IGFBPs, CYR61) are key regulators of pathways controlling the commitment of MSCs to different lineages [START_REF] Kiepe | Differential expression of IGF system components in proliferating vs. differentiating growth plate chondrocytes: the functional role of IGFBP-5[END_REF][START_REF] Schutze | Differential expression of CCN-family members in primary human bone marrow-derived mesenchymal stem cells during osteogenic, chondrogenic and adipogenic differentiation[END_REF][START_REF] Van Der Horst | Downregulation of Wnt signaling by increased expression of Dickkopf-1 and -2 is a prerequisite for late-stage osteoblast differentiation of KS483 cells[END_REF]. The overlap also contains numerous genes whose products are involved in the cytoskeleton and matrix remodeling associated with mesenchymal cell biology (Table 1). In agreement with the proposed endothelial potential of MSCs, some genes modulated by EWS-FLI1 and expressed by MSCs have critical role in angiogenesis (Table 1). Finally, the inhibition of EWS-FLI1 also induced the expression of some master genes of the mesenchymal differentiation (including SOX9 and BMP1) that are not expressed in undifferentiated MSCs.

Together, these results strongly suggested that the EWS-FLI1-silenced Ewing cells recover an expression pattern related to that of MSCs. They also showed that a number of neural genes expressed in ETs, but not in MSCs, are regulated by EWS-FLI1 (Table 1).

EWS-FLI1-silenced Ewing cells harbor mesenchymal features.

From these in silico analyses, we postulated that terminal differentiation in the mesenchymal lineages might be achieved upon combined effects of EWS-FLI1

silencing and treatment of cells with appropriate differentiation cocktails. To test this hypothesis, EWS-FLI1 silenced EW24, SK-N-MC and A673 cells were used to conduct in vitro adipogenic differentiation experiments (Figure 2). Upon differentiation conditions, EW24 infected with the control lentivirus and hence harboring a normal EWS-FLI1 expression exhibited a mild induction of PPARγ2 and FABP4 and a strong induction of LPL, three specific markers of the adipocyte lineage (Figure 2A). The silencing of EWS-FLI1 in EW24 cells grown in standard or differentiation media led to a much stronger induction of FABP4 and to a lesser extent of PPARγ2. Oil-Red-O staining, which labels lipidic vesicles, was not significant for EW24 infected with the control virus either in basal or differentiation conditions. The inhibition of EWS-FLI1 by itself induced a light Oil-Red-O staining of a few cells (Figure 3B). However, this staining was much more obvious when EWS-FLI1-silenced cells were incubated in the presence of differentiation medium (Figure 2B). These results indicated that EW24 cells exhibit some degree of adipogenic conversion upon incubation with specific cocktail or upon EWS-FLI1 silencing but that this commitment is much more pronounced when both conditions are combined. A similar observation was made with A673 (data not shown and see below). In contrast, neither induction of markers, nor Oil-Red-O staining was observed in the SK-N-MC cell line.

In addition to adipogenic differentiation, MSC can also commit to the osteogenic or chondrogenic lineages. However, these differentiations necessitate up to several weeks of incubation with specific media, which is hardly compatible with the si or shRNA approaches used in the previous experiments. Indeed, due to a progressive overgrowth of poorly inhibited cells, a long term and stable EWS-FLI1

silencing could not be achieved with these systems. To circumvent this technical constrain we constructed Ewing cell lines expressing a doxycycline (DOX) inducible shRNA targeting EWS-FLI1. Two independent clones (shA673-1C and -2C) were isolated. As shown in figure 3A, the induction of the shRNA led to a dramatic reduction of EWS-FLI1 transcript and protein levels observed by 24h after addition of DOX in the culture medium and remaining stable over time. In agreement with previous observations based on siRNAs [START_REF] Chansky | Targeting of EWS/FLI-1 by RNA interference attenuates the tumor phenotype of Ewing's sarcoma cells in vitro[END_REF][START_REF] Prieur | EWS/FLI-1 silencing and gene profiling of Ewing cells reveal downstream oncogenic pathways and a crucial role for repression of insulin-like growth factor binding protein 3[END_REF], inducible EWS-FLI1 silencing resulted in a considerable decrease of cell proliferation associated with a reduced number of cells in the S-phase of the cell cycle and a simultaneous increase of the amount of G1-phase cells (Figure S2). We first confirmed that adipogenic differentiation conducted with both inducible clones led to results similar to those described above. In the presence of EWS-FLI1 (-DOX), differentiation medium led to a significant induction of the three adipogenic markers (Figure 3B). The knock-down of EWS-FLI1 in standard medium also led to an induction of these transcripts. Nevertheless, the increased expression of these three genes was higher when differentiation was conducted in the presence of DOX, suggesting that inhibition of EWS-FLI1 and differentiation cocktail may provide additive or even synergistic effects (see PPARγ2 and LPL). Finally, Oil-Red-O staining was positive only for cells treated by both DOX and the differentiation cocktail (Figure 3C). Altogether these data show that EWS-FLI1-inhibited Ewing cells can accumulate lipid vesicles and express adipogenic markers, two features highly suggestive of a differentiation along the adipocyte lineage.

Similar experiments were conducted for osteogenic differentiation using quantitative RT-PCR analyses of three markers specific for this lineage (SPP1, ALPL and RUNX2) and Von-Kossa staining. The EWS-FLI1 silencing and subsequent growth arrest observed with DOX treatment of shA673 cells were stable over time and long-term differentiation experiments could thus be conducted. The three aforementioned osteogenic markers were induced by DOX and/or differentiation mix with different patterns. The induction of ALPL was clearly dependent upon the medium but not influenced by EWS-FLI1 (Figure 3D). Contrarily, differentiation medium and DOX exhibited additive or even synergistic effects for RUNX2 and SPP1 inductions.

Even more strikingly, the Von-Kossa staining of shA673 grown for three weeks in differentiation medium and DOX demonstrated the presence of a calcified matrix whereas this was observed neither in cells grown in differentiation medium without DOX nor in standard medium with or without DOX (Figure 3E). Together these results show that EWS-FLI1-silenced Ewing cells can exhibit characteristic features of osteocytes.

As MSCs are also precursors of chondrocytes, we tested whether shA673 cells could differentiate along this lineage. Pellet cultures of shA673 cells were performed either in presence or absence of DOX and chondrogenic differentiation medium. Whereas EWS-FLI1-expressing shA673 cells remained flat and rapidly died, EWS-FLI1-silenced shA673 formed a round pellet, but ultimately died prior to the delay necessary to evaluate terminal chondrocyte differentiation. However, quantitative RT-PCR experiments showed that, in addition to the afore-mentioned RUNX2 induction, COL10A1 and SOX9, two markers of the chondrocyte lineage, were upregulated upon EWS-FLI1 silencing (Figure 3D). These results suggested that Ewing cells can engage the initial steps of chondrocyte differentiation but cannot fully complete this process.

Finally, since MSCs can be characterized by the expression of several surface antigens [START_REF] Covas | Mesenchymal stem cells can be obtained from the human saphena vein[END_REF][START_REF] Izadpanah | Characterization of multipotent mesenchymal stem cells from the bone marrow of rhesus macaques[END_REF][START_REF] Oswald | Mesenchymal stem cells can be differentiated into endothelial cells in vitro[END_REF][START_REF] Pittenger | Multilineage potential of adult human mesenchymal stem cells[END_REF], we investigated the expression of 10 surface markers using FACS analyses in the presence or absence of EWS-FLI1 (Figure 3F). Five of these (CD44, CD59, CD73, CD29 and CD54) were expressed on shA673 cells and induced by EWS-FLI1 inhibition. Three (CD90, CD105 and CD166) were expressed but not modified by EWS-FLI1 modulation. Finally, CD45 and CD31, two markers expressed in hematopoeitic stem cells or endothelial cells, but not in mesenchymal stem cells,

were not detected at the surface of Ewing cells whatever the EWS-FLI1 status.

Altogether these results strengthen the hypothesis of EWS-FLI1-silenced Ewing cells sharing phenotypic characteristics with MSCs.

Discussion

The fusion between EWS and members of the ETS family, usually FLI1, is a critical event of Ewing sarcoma development. EWS-ETS fusions are only observed in tumors of the Ewing family, suggesting that transformation is dependent upon the specific action of this oncogene within a precise cell background. Consequently, a thorough phenotypic analysis of EWS-FLI1-silenced Ewing cells should give insights into this cellular context. We have addressed this issue investigating the transcription profile of EWS-FLI1-silenced Ewing cells. Numerous genes involved in neural development and differentiation are expressed in ET and repressed upon inhibition of EWS-FLI1. This observation is in complete agreement with the occasional presence of neural markers in ETs and with the recently reported shift from myogenic to neural phenotype induced by EWS-FLI1 expression in rhabdomyosarcoma cells [START_REF] Hu-Lieskovan | EWS-FLI1 fusion protein up-regulates critical genes in neural crest development and is responsible for the observed phenotype of Ewing's family of tumors[END_REF]. It confirms that EWS-FLI1 by itself can promote the induction of some neural crest genes. However, our data show that the neural profile of Ewing cells is not sufficient for ETs to cluster among neural tissues (Figure 1A).

The more striking in silico observation is that three different EWS-FLI1inhibited Ewing cell lines cluster with MSCs when the analysis is performed with EWS-FLI1-modulated genes. MSCs were initially described by Friedenstein [START_REF] Friedenstein | Fibroblast precursors in normal and irradiated mouse hematopoietic organs[END_REF] who isolated plastic adherent cells from bone marrow capable of differentiating into adipocytes, chondrocytes and osteocytes. More recently, a number of investigators have shown that multi-potent MSCs can be isolated from a variety of sources including bone marrow, adipose tissue, connective tissues, placenta and umbilical cord blood (for review see [START_REF] Dennis | Origin and differentiation of human and murine stroma[END_REF][START_REF] Gregory | Non-hematopoietic bone marrow stem cells: molecular control of expansion and differentiation[END_REF]. In addition to the afore-mentioned lineages, MSCs can also give rise to stromal cells supporting hematopoiesis and also probably to endothelial and vascular smooth muscle cells. Our study show that markers currently used for the isolation of MSCs may be induced upon EWS-FLI1 inhibition (CD44, CD54, CD59 and CD73). The in silico analysis described above is strongly supported by experiments showing that EWS-FLI1-inhibited Ewing cells, but not parental Ewing cells, exhibit differentiation characteristics of mesenchymal lineages. Indeed, two out of three Ewing cell lines can be clearly differentiated along the adipogenic lineage.

This indicates that at least in these two cell lines, the apparatus responsible for the adipogenic differentiation is present but impaired by EWS-FLI1. Long term inhibition of EWS-FLI1 also shows that terminal osteogenic differentiation can be observed with synthesis of a calcified matrix. These experimental data are fully consistent with the previous description of osteoid matrix in a subgroup of EWS-FLI1-positive tumors called small cell osteosarcoma [START_REF] Oshima | Abdominal small round cell tumor with osteoid and EWS/FLI1[END_REF]. A less evident commitment of EWS-FLI1 silenced cells toward the chondrogenic lineage is also observed. Despite the death of Ewing cells observed upon chondrocyte differentiation conditions, impairing the evaluation of the terminal differentiation, COL10A1, SOX9 and RUNX2, three markers suggestive of a commitment toward the chondrocyte lineage are induced upon EWS-FLI1 silencing. SOX9 can be induced in chondrogenic but also in neural differentiation. To further investigate the hypothesis of SOX9 induction being related to differentiation along the neural lineage, we also analyzed the expression of SNAI2, FOXD3 and PAX7, three markers which are required for neural crest cells differentiation [START_REF] Basch | Specification of the neural crest occurs during gastrulation and requires Pax7[END_REF]. PAX7 was strongly (10X) inhibited by EWS-FLI1

silencing. FOXD3 and SNAI2 were slightly inhibited or induced, respectively.

Interestingly, none of these markers were induced upon treatment of Ewing cells with DMSO/HBA or DMSO/RA, two protocols previously described for neural differentiation [START_REF] Scintu | Differentiation of human bone marrow stem cells into cells with a neural phenotype: diverse effects of two specific treatments[END_REF][START_REF] Woodbury | Adult bone marrow stromal stem cells express germline, ectodermal, endodermal, and mesodermal genes prior to neurogenesis[END_REF][START_REF] Woodbury | Adult rat and human bone marrow stromal cells differentiate into neurons[END_REF].

These results strongly suggest that the SOX9 induction is rather related to a chondrocyte than a neural commitment.

Altogether, these in silico and experimental data strongly suggest that Ewing The MSCs compartment consists of heterogeneous population of cells and presently, no single marker can be used to isolate subpopulations of MSCs committed to specific lineages. Stem cells with an extended differentiation potential, in particular in the angiogenic lineage, can be isolated from adipose tissue and bone marrow. In that respect, it is striking that different genes involved in angiogenesis are modulated by EWS-FLI1. The hypothesis of a vascular differentiation potential of ET cells, which was proposed by James Ewing himself, is further supported by recent microarray data and by the observation of vascular mimicry in ET associated with the ability of some Ewing cell lines to form vascular-like tubes in vitro [START_REF] Staege | DNA microarrays reveal relationship of Ewing family tumors to both endothelial and fetal neural crest-derived cells and define novel targets[END_REF][START_REF] Van Der Schaft | Tumor cell plasticity in Ewing sarcoma, an alternative circulatory system stimulated by hypoxia[END_REF]. In that respect, it is interesting to mention that FLI1, the wild type counterpart of EWS-FLI1, is a key regulator of early blood vessel formation [START_REF] Brown | Insights into early vasculogenesis revealed by expression of the ETS-domain transcription factor Fli-1 in wild-type and mutant zebrafish embryos[END_REF]. (B) Same analysis than in (A) but performed with EWS-FLI1-modulated probe sets (SAM analyses between efRNAi and ctrlRNAi samples, absolute fold change >3; qvalue <0.01). The two first components account for 41.8% of the variation.

(C) Hierarchical clustering analysis (average linkage method, 1-rank correlation distance metric) using EWS-FLI1-modulated probe sets. efRNAi samples cluster with mesenchymal tissues and stem cells. The complete annotated cluster is provided as

(D) A significant proportion of genes differentially expressed between efRNAi and ctrlRNAi also distinguishes MSCs from ET. The blue circle schematizes the set of genes exhibiting significant variation (>2 fold; q-value<0.01) between EWS-FLI1silenced or -expressing Ewing cells (efRNAi/ctrlRNAi) and the red circle represents the set of genes differentially expressed between MSCs and ETs (the complete list of genes is provided in Table S2). The table represents the distribution of the 110 genes found in the overlap. (F) FACS analyses of four cell surface markers (CD59, CD73, CD54 and CD44).

Cells were treated (black lines) or not (grey lines) with DOX during 6 days prior to FACS analyses. Dashed lines represent control IgE.

  antisense) were annealed and then cloned between the BglII and HindIII restriction sites of the pTER vector (van de[START_REF] Van De Wetering | Specific inhibition of gene expression using a stably integrated, inducible small-interfering-RNA vector[END_REF]. Expression plasmids were transfected in A673 Ewing cell line with the effectene transfection reagent (Qiagen, Paris, France) according to the manufacturer's instructions. A Tet repressor expressing clone (A673tetR5BIII) was selected on Blasticidine (20µg/mL) then transfected with the pTER/shEF-1 construct encoding Zeocin resistance (200µg/mL). Two clones (shA673-1C and -2C), exhibiting strong inhibition of EWS-FLI1 upon treatment with 1µg/mL Doxycycline (DOX), were selected. For lentivirusmediated silencing of EW24 and SK-N-MC cell lines, the H1 promoter-tetO-shEF-1

(

  Hs00609791_m1); LPL (Hs00173425_m1); PPARG (Hs00234592_m1); ALPL (Hs00758162_m1); RUNX2 (Hs00231692_m1); SPP1 (Hs00167093_m1); SOX9 (Hs00165814_m1); COL10A1 (Hs00166657_m1), PAX7 (Hs00242962_m1), FOXD3 (Hs00255287_m1), SNAI2 (Hs00161904_m1) and RPLP0 (Hs99999902_m1).

  cells arise from MSC which terminal mesenchymal differentiation is blocked by EWS-FLI1. It also shows that the phenotype of Ewing cells is perverted from the mesenchymal to neural lineage by EWS-FLI1. However, some Ewing cell lines (SK-N-MC), despite a strong convergence toward a MSC gene expression profile upon EWS-FLI1 silencing do not differentiate, at least along the adipogenic lineage. This indicates that the Ewing cell background, although related to MSC, harbor cell linespecific differences which may reflect diverse mesenchymal potentials of the original progenitors in which the EWS-FLI1 fusion occurred. Alternatively, some Ewing cell lines may have acquired secondary genetic changes impairing differentiation processes.

Figure legends Figure 1 .

 legends1 Figure legends

  99 genes behave consistently between MSC/ET and ef/ctrlRNAi: 44 genes upregulated in MSCs as compared to ETs are downregulated by EWS-FLI1 and 55 genes downregulated in MSCs as compared to ETs are upregulated by EWS-FLI1. (E) Venn diagram showing, among the 2540 genes which discriminates MSCs and ETs, the number of genes regulated by EWS-FLI1 in a cell line specific manner. For each cell line, the number of EWS-FLI1-modulated genes (Welch test p-value<0.01, absolute fold change >2) common with the gene set of the MSC/ET comparison is represented. Among the 782, 364 and 382 genes regulated by EWS-FLI1 in A673, EW24 or SK-N-MC, 320, 93 and 82 are common with the MSC/ET comparison, respectively.

Figure 2 .

 2 Figure 2. Differentiation of EWS-FLI1-silenced EW24 cells along the adipogenic

Figure 3 .

 3 Figure 3. Long term silencing of EWS-FLI1 induces mesenchymal features.

  Taken together, in silico and experimental results suggest that ETs are derived from MSCs that may have both mesenchymal and vascular potentials. ET may therefore constitute an appropriate model to further study the characteristics of a common mesenchymal and vascular progenitor.

	The identification of the Ewing primary cell constitutes a key step toward a
	better understanding of ET biology. The MSC origin may account for the predominant
	localization of ET in bones or soft tissues, two major sources for these stem cells.
	Through a precise spatio-temporal targeting of EWS-FLI1 expression in the mouse, it
	should foster the development of a Ewing tumor model, which is presently lacking
	and which may be of considerable help to achieve comprehensive insights into tumor
	initiation and progression processes.

Table 1 .

 1 Main genes modulated by EWS-FLI-1 and differentially expressed between MSC and ET

	Symbol	ET/MSC	A673 a	SKNMC a EW24 a	Symbol	ET/MSC	A673 a	SKNMC a	EW24 a
	Mesenchymal differentiation				Extracellular matrix and cytoskeleton
	CALD1	-13.8	-3.5	-3.4	-3.2	ARHGAP29	-22.1	-4.9	-2.8	-5.2
	CYR61	-21.3	-6.0	-10.7	-11.6	ITGB5	-10.9	-2.7		-2.9
	DKK1	-2.2	-33.7	-11.5	-24.3	MARCKS	-4.1	-2.2		-2.2
	IGFBP5	-7.0	-25.8	-8.7	-8.5	MID1	-29.7		-2.0	-2.6
	IGFBP7	-3.8	-29.8	-17.9	-18.9	PALLD	-11.3	-2.4	-2.7	-3.4
	PPP1R1A	152.9	12.9		4.8	PLAT	-26.5	-11.1		-2.1

a Mean fold change between ctrlRNAi and efRNAi samples. -indicates EWS/FLI-1-repressed genes.

Acknowledgments

We thank Zofia Maciorowski from the department of tumor biology of the Institut Curie for the BrdU analyses, Hans Clevers for the generous gift of the pTER and the Tet repressor expression plasmid (pTet-R) and Pierre Charneau for providing us the pTRIP lentivirus. We also thank Gilbert Lenoir and Steve Lessnick for materials used in this study. This work was supported by grants from the Institut National de la Santé et de la Recherche Médicale, the Institut Curie, the Ligue Nationale contre le Cancer (Equipe labellisée and CIT program), the Réseau National des Génopoles and the Association for International Cancer Research (AICR).


Valentina Boeva

Didier Surdez

Noëlle Guillon

Franck Tirode

Anthony P Fejes

Olivier Delattre

Emmanuel Barillot

Noe ¨lle Guillon

De novo motif identification improves the accuracy of predicting transcription factor binding sites in ChIP-Seq data analysis

De novo motif identification improves the accuracy of predicting transcription factor binding sites in ChIP-Seq data analysis

INTRODUCTION

The appearance of next-generation sequencing technologies has propelled forward the development of new techniques among which ChIP-Seq has become an important method for genome-wide discovery of binding sites for DNA-associated proteins and in particular for TFBSs. ChIP-Seq consists of the immunoprecipitation of protein-DNA complexes followed by massively parallel sequencing of short ends of immunoprecipitated DNA [START_REF] Barski | High-resolution profiling of histone methylations in the human genome[END_REF][START_REF] Johnson | Genome-wide mapping of in vivo protein-DNA interactions[END_REF][START_REF] Robertson | Genome-wide profiles of STAT1 DNA association using chromatin immunoprecipitation and massively parallel sequencing[END_REF]. This technique succeeded the ChIP-on-chip technique (4) and has nearly replaced the latter because of the increased accuracy in identification of TFBSs [START_REF] Johnson | Genome-wide mapping of in vivo protein-DNA interactions[END_REF].

At the completion of a ChIP-Seq experiment, millions of short ($35-50 bp) directional DNA tags are obtained, which can be positioned or aligned to the reference genome for the sample organism (Supplementary Figure S1). Each short tag represents an extremity of a longer DNA fragment ($200-400 bp depending on the experiment) isolated from the immunoprecipitation. Thus, in the analysis of the short representative tags, it is important to take this experimental fact into consideration to identify the full length of the original fragment that gave rise to the tag. By extending each tag, it is then possible to identify areas of overlap, which represent the location of the protein binding event. The density profile of DNA fragment coverage can then be calculated and 'peaks' corresponding to putative binding sites can be extracted. This idea was elegantly implemented in the FindPeaks software [START_REF] Fejes | FindPeaks 3.1: a tool for identifying areas of enrichment from massively parallel short-read sequencing technology[END_REF]. However, the accuracy of peak calling can be considerably improved by incorporating information about genomic sequences of peaks in addition to coverage depth information.

In this article we present an algorithm implemented in the MICSA software (Motif Identification for ChIP-Seq Analysis) that is based on the idea that functional binding sites of transcription factors (TFs) should contain a consensus motif (or a set of motifs). Consensus motifs are the composite sequences of DNA for which a DNA-binding protein, such as a TF or restriction enzyme, has a high affinity. Such motifs can be identified from the small subset of peaks with a high DNA fragment coverage.

The MICSA algorithm is innovative in the context of ChIP-Seq data analysis for simultaneous: (i) de novo TFBS motif identification and (ii) functional binding site prediction using information about motif occurrences in peaks along with coverage depth information. Here, motif identification is not a post-processing step as in other ChIP-Seq analysis pipelines (6) but a key element which allows keeping even low peaks if they have a strong motif occurrence.

Since MICSA checks for motif occurrences in all peaks including those with very low coverage depth, there is no need in the explicit selection of threshold on DNA tag/ fragment coverage. The only parameter that remains to be specified is the maximal number of expected false positive hits among selected peaks or the maximal false discovery rate (FDR).

Using the procedure developed by Kharchenko et al. [START_REF] Kharchenko | Design and analysis of ChIP-seq experiments for DNA-binding proteins[END_REF], we compared the peak identification performance of MICSA and 10 other published tools [START_REF] Fejes | FindPeaks 3.1: a tool for identifying areas of enrichment from massively parallel short-read sequencing technology[END_REF][START_REF] Ji | An integrated software system for analyzing ChIP-chip and ChIP-seq data[END_REF][START_REF] Kharchenko | Design and analysis of ChIP-seq experiments for DNA-binding proteins[END_REF][START_REF] Zhang | Model-based analysis of ChIP-Seq (MACS)[END_REF][START_REF] Rozowsky | PeakSeq enables systematic scoring of ChIP-seq experiments relative to controls[END_REF][START_REF] Valouev | Genome-wide analysis of transcription factor binding sites based on ChIP-Seq data[END_REF][START_REF] Nix | Empirical methods for controlling false positives and estimating confidence in ChIP-Seq peaks[END_REF][START_REF] Boyle | F-Seq: a feature density estimator for high-throughput sequence tags[END_REF][START_REF] Mortazavi | Mapping and quantifying mammalian transcriptomes by RNA-Seq[END_REF][START_REF] Jothi | Genome-wide identification of in vivo protein-DNA binding sites from ChIP-Seq data[END_REF]. The dataset selected for the comparison was generated by Johnson et al. [START_REF] Johnson | Genome-wide mapping of in vivo protein-DNA interactions[END_REF] for the neuron-restrictive silencer factor (NRSF). MICSA showed a considerable increase in the performance over 10 other approaches. To increase the statistical basis we performed the same comparison procedure for selected algorithms on other ChIP-Seq datasets, including those for GA-binding protein (GABP) [START_REF] Valouev | Genome-wide analysis of transcription factor binding sites based on ChIP-Seq data[END_REF], signal transducer and activator of transcription 1 (STAT1) [START_REF] Rozowsky | PeakSeq enables systematic scoring of ChIP-seq experiments relative to controls[END_REF] and CCCTC-binding factor (CTCF) [ENCODE project, the Broad Institute and the Bradley E. Bernstein lab at the Massachusetts General Hospital/ Harvard Medical School [START_REF] Birney | Identification and analysis of functional elements in 1% of the human genome by the ENCODE pilot project[END_REF]]. The results of the comparison indicated that use of MICSA for ChIP-Seq data analysis allows us to significantly reduce the number of false positive predictions for TFBSs.

The MICSA package was also used on our ChIP-Seq data [START_REF] Guillon | The oncogenic EWS-FLI1 protein binds in vivo GGAA microsatellite sequences with potential transcriptional activation function[END_REF]. Immunoprecipitation was performed with a specific antibody directed against the oncogenic TF EWS-FLI1 [START_REF] May | Ewing sarcoma 11;22 translocation produces a chimeric transcription factor that requires the DNA-binding domain encoded by FLI1 for transformation[END_REF] to obtain biological insight into the functioning of this TF, which is known to be the major oncogene in Ewing sarcoma. Using our technique, based on motif identification, we confirmed the existence of two consensus motifs, one representing a (GGAA) n microsatellite, and the second containing the RCAGGAARY consensus sequence (16) (R = A/G, Y = T/C). Further analysis of the EWS-FLI1 data, together with expression arrays, suggested that EWS-FLI1 bound to (GGAA) n microsatellites can activate transcription of neighboring genes; while EWS-FLI1 bound to RCAGGAARY sites may, depending on genes, activate or repress transcription. Our analysis confirmed five known direct target genes of EWS-FLI1 and has also predicted many new genes that are putatively regulated directly by EWS-FLI1.

The algorithm we developed is pioneering in its use of motif information when predicting sites of specific binding for TFs from ChIP-Seq data. It allows identification of several motifs which, as shown at the EWS-FLI1 example, can possibly carry different biological function.

MATERIALS AND METHODS

Candidate peak identification

Candidate peaks are identified by FindPeaks [START_REF] Fejes | FindPeaks 3.1: a tool for identifying areas of enrichment from massively parallel short-read sequencing technology[END_REF] which is included in the MICSA package. The java class DeleteRegions is then used to eliminate peaks in satellite regions which are often a source of noise in ChIP-Seq data. We propose to use blacklist regions masked as centromeric repeats by RepeatMasker (http://www.repeatmasker.org) in UCSC genome browser [START_REF] Karolchik | The UCSC Genome Browser Database: 2008 update[END_REF]. Files are provided in the MICSA package.

Files with candidate peaks are then processed by the java class FilterPeaks in order to filter out false peaks occurring both in ChIP and control datasets because of biases due to PCR errors. MICSA needs to be supplied with a reference genome from which it would extract DNA sequences of peaks.

Classes

The whole set of candidate peaks is decomposed into classes. Class C i contains regions from which i overlapping DNA fragments were immunoprecipitated during ChIP experiment. For peaks from class C i the probability p i to be a false binding site can be evaluated using control data: p i = min{1,(#peaks in control class C' i )/(#peaks in ChIP class C i )}. Peaks highly enriched in mapped DNA tags are likely to be true positives, while peaks with low DNA tag coverage can probably be false since lots of peaks with the same values of coverage could be identified in the control data. Thus, p i will be smaller for higher values of i.

Here our hypothesis is that the normalized number of peaks in the control dataset is an estimate of the number of false binding sites in our ChIP dataset. For the same purpose, other methods use, for example, Monte-Carlo simulations [START_REF] Fejes | FindPeaks 3.1: a tool for identifying areas of enrichment from massively parallel short-read sequencing technology[END_REF] or Poisson approximation [START_REF] Zhang | Model-based analysis of ChIP-Seq (MACS)[END_REF]. We believe that using the control data is more appropriate in this case since it takes into account the sequencing bias and guarantees not to underestimate the number of false binding sites.

Motif identification

We use a subset of the highest peaks from classes with zero p i value to identify over-represented motifs. Moreover, we do not consider whole peak sequences but only sequences at the location of the maximal enrichment in each peak. The MEME motif finding tool [START_REF] Bailey | Fitting a mixture model by expectation maximization to discover motifs in biopolymers[END_REF] is then run automatically on this set of sequences in order to identify the most over-represented motif. Sequences from this set which do not contain the identified motif are subject to the second MEME run. Finally, one uses the top significance motifs, each in the form of a position-specific scoring matrix (PSSM) with the minimal threshold value, which occur in the areas of maximal enrichment of the high peaks.

It is only one motif in the case when the protein does not undergo allosteric modifications which could change its binding motif, e.g. distance between its half-sites; one obtains more than one motif in the case when such change is possible. A large part of genome is scanned to yield the real frequencies of extracted motifs.

Optimization

In what follows below, we consider separately peaks on different chromosomes. For motif M and class C i we will call peaks which have an occurrence of M with a PSSM score above T M,i within the top peak area defined by (h-M,i ), where h is a peak height, i.e. maximal number of overlapping DNA fragments coming from this region. The question is how to choose T M,i and M,i for each particular class. It can be done through the calculation of expected number of false binding sites selected by the peak calling procedure given the null hypothesis of the Markov(0) model for nucleotide distribution in peaks. For each pair (T M,i, M,i ) and the class C i we can evaluate the expectation E i (T M,i, M,i ) of number of peaks called by chance. This value is equal to a sum of the motif P-values of each peak in C i . The motif P-value means a probability to observe a motif by chance in a sequence of given length. We use the Poisson approximation for the P-values: P-value &1 -(1 -MotifProbability) TextLength -MotifLength + 1 . Here, the motif probability, which is the probability of observing a motif occurrence above a given PSSM threshold on a given position, is considered to be equal to the genomic frequency of the motif with the given threshold. Finally, since p i is the estimate of the probability of any given peak in class C i to be a false binding site, the product E i (T M,i, M,i )Áp i estimates the total number of false binding sites which would be called by this selection procedure. The number of selected peaks S i in class C i also depends on (T M,i, M,i ). Our optimization procedure maximizes AES i so that AEE i p i stays below the predefined threshold on the total number of false positives.

Score calculation

A score is reported for each peak selected during the optimization procedure. It is equal to the product of p i of corresponding class and the smallest motif P-value among motif P-values of all motif occurrences in the peak region.

The second term is small if the motif is situated near the location of the maximal enrichment of the peak and its sequence is close to the consensus.

Additional methods

Information about ChIP-Seq library construction and sequencing for EWS-FLI1 data was published in [START_REF] Guillon | The oncogenic EWS-FLI1 protein binds in vivo GGAA microsatellite sequences with potential transcriptional activation function[END_REF]. More detail on the optimization procedure is available in Supplementary Data.

Software availability

The MICSA algorithm is implemented as a Java package with a graphical user interface in Java. It is freely available for nonprofit use at http://bioinfo-out.curie.fr/projects/ micsa/. All data presented in this study (mock control and ChIP-Seq data, and peak call coordinates) are available at the same website.

RESULTS

Theoretical framework

The key idea of MICSA is to use the information about TFBS motif occurrences to predict true binding sites from ChIP-Seq data. The MISCA workflow consists of four main steps (Figure 1): (i) identify all candidate peaks, (ii) identify TFBS motifs from a small subset of peaks, (iii) find motif occurrences in candidate peaks, and (iv) optimize the peak calling output by calculating statistics that take into account both motif occurrence and depth of coverage information. The last step roughly corresponds to filtering out of insignificant peaks without motif occurrences.

The first step of the MICSA algorithm consists in identifying all regions that are enriched in mapped DNA tags (Figure 1). Such areas of enrichment, also called peaks, represent potential binding sites for the protein of interest and can often exceed tens of thousands of locations. To detect these regions we use the previously developed tool FindPeaks [START_REF] Fejes | FindPeaks 3.1: a tool for identifying areas of enrichment from massively parallel short-read sequencing technology[END_REF]. One of the reasons that we use FindPeaks is that it generates UCSC compatible custom 'WIG' track files from aligned-read files. We calculate the false discovery rate (FDR) predictor to estimate the proportion of false peaks to called peaks for each dataset. FDR is estimated using both the ChIP and mock control data. Generally, the FDR estimate for high peaks (peaks with greater depth of DNA tag or fragment coverage) is smaller than for low ones since the former are less likely to occur randomly or in control data.

The second step is de novo motif identification. DNA sequences of peaks are extracted from the reference genome using position information. In order to limit the In the third step, discovered motifs are identified for all regions of enrichment located by FindPeaks (peaks file). The combined information on motif location and position-specific scoring matrix (PSSM) scores is then used to calculate a motif P-value for each peak (see 'Materials and Methods' section). In general, the P-value will be small if the motif occurrence closely resembles the consensus (and thus has a high PSSM score) or if the motif occurs close to the region of maximum enrichment of the peak. In contrast, the motif P-value will be high if the motif occurs in the periphery of the peak and the motif observed has a low PSSM score. The final peak score appearing in the MICSA output is a product of the motif P-value described above and the estimated FDR. Supposing that, for a given threshold on the peak height, the number of false binding sites in the ChIP sample is equal to the normalized number of peaks in the control sample, and nucleotides along such peak sequences are distributed according to the Markov(0) model, then MICSA's peak score estimates the probability of the given peak to be a false binding site and have a motif occurrence just by chance.

In the fourth step we select binding sites out of candidate peak dataset based on motif occurrence and depth of coverage information. The raw peak candidate prediction set contains a large number of false positive predictions. To filter them out, other methods commonly determine a threshold on DNA tag/fragment coverage [START_REF] Ji | An integrated software system for analyzing ChIP-chip and ChIP-seq data[END_REF][START_REF] Valouev | Genome-wide analysis of transcription factor binding sites based on ChIP-Seq data[END_REF][START_REF] Mortazavi | Mapping and quantifying mammalian transcriptomes by RNA-Seq[END_REF]. However, our experience and that of other researchers (C. Wadelius, personal communication) shows that even regions with relatively low DNA tag coverage often contain functional binding sites. Here we propose an approach based on the presence of motif information in peak sequences. We use this information to retain additional peaks containing strong motifs, especially those with low DNA tag coverage. For each subset of peaks with a given depth of DNA fragment coverage, we choose a criterion for peak retention based on motif strength and position within the peak. We then evaluate a number of false positive peaks that we expect to satisfy the criterion (see 'Materials and Methods' section). We use an optimization procedure that applies different criteria in an attempt to maximize the total number of selected peaks for all subsets without having the estimated total number of false positive peaks passing the selection exceed the user defined threshold. MICSA outputs the list of peaks which are selected by the optimization procedure and ranked according to their scores and associated motifs. As a result of the optimization, high peaks are usually kept by MICSA even if there is no strong motif hit. In other words, though the motif based filtering is applied by MICSA to all peaks, it is only effective on lower ones (Supplementary Figure S2).

Performance of MICSA

In order to assess relative performances of MICSA and other existing tools, we used the algorithm comparison method developed by Kharchenko et al. [START_REF] Kharchenko | Design and analysis of ChIP-seq experiments for DNA-binding proteins[END_REF]. The tested programs were MACS (8), PeakSeq (9), QuEST [START_REF] Valouev | Genome-wide analysis of transcription factor binding sites based on ChIP-Seq data[END_REF], wdt [START_REF] Kharchenko | Design and analysis of ChIP-seq experiments for DNA-binding proteins[END_REF], Useq [START_REF] Nix | Empirical methods for controlling false positives and estimating confidence in ChIP-Seq peaks[END_REF], F-Seq [START_REF] Boyle | F-Seq: a feature density estimator for high-throughput sequence tags[END_REF], CisGenome (6), ERANGE/ChIPSeqMini (13), SISSRs [START_REF] Jothi | Genome-wide identification of in vivo protein-DNA binding sites from ChIP-Seq data[END_REF] and FindPeaks [START_REF] Fejes | FindPeaks 3.1: a tool for identifying areas of enrichment from massively parallel short-read sequencing technology[END_REF]. Each application was used with its default parameters, according to the instructions given in the manual (see Supplementary Table S1 for command lines). In cases where only a small number of peaks were extracted, parameters were modified in order to increase the number of identified peaks. ChIP-Seq data for the NRSF (2,20) were selected for testing because they had already been widely used by other groups to validate ChIP-Seq analysis software [START_REF] Johnson | Genome-wide mapping of in vivo protein-DNA interactions[END_REF][START_REF] Ji | An integrated software system for analyzing ChIP-chip and ChIP-seq data[END_REF][START_REF] Kharchenko | Design and analysis of ChIP-seq experiments for DNA-binding proteins[END_REF][START_REF] Zhang | Model-based analysis of ChIP-Seq (MACS)[END_REF][START_REF] Nix | Empirical methods for controlling false positives and estimating confidence in ChIP-Seq peaks[END_REF][START_REF] Jothi | Genome-wide identification of in vivo protein-DNA binding sites from ChIP-Seq data[END_REF]. Additionally, we run MICSA on datasets for GABP [START_REF] Valouev | Genome-wide analysis of transcription factor binding sites based on ChIP-Seq data[END_REF], STAT1 [START_REF] Rozowsky | PeakSeq enables systematic scoring of ChIP-seq experiments relative to controls[END_REF] and CTCF [ENCODE project, the Broad Institute and the Bradley E. Bernstein lab at the Massachusetts General Hospital/Harvard Medical School [START_REF] Birney | Identification and analysis of functional elements in 1% of the human genome by the ENCODE pilot project[END_REF]].

For NRSF we used two positive TFBS sets to assess the sensitivity of each of the different methods. The first one is a list of 3000 high-scoring motif instances designed using canonical sequence binding motifs for NRSF by Kharchenko et al. [START_REF] Kharchenko | Design and analysis of ChIP-seq experiments for DNA-binding proteins[END_REF]. The second is composed of 83 binding sites verified by qPCR (2). To compare the methods' sensitivities, we selected increasing numbers of top peaks returned by each method and analyzed the fraction of peaks containing the motif of interest (Figure 2 and Supplementary Figure S3). MICSA clearly outperforms other methods on both tested datasets, with almost any threshold on a number of called peaks. We compared algorithm performances considering the best 3000 peaks called by each program. For the first positive set of 3000 high-scoring motif instances we found that 1422 of them were identified by MICSA (Figure 2A). According to this test we could rank the other tools as uSeq (1254), wdt (1229), PeakSeq (1227), F-Seq (1217), FindPeaks 3.3 (1216), CisGenome (1203), MACS (1195), SISSRs (1194), QuEST 2.0 (1132), ERANGE 3.1 (1118) (Figure 2A, Supplementary Figure S3a and Supplementary Table S2). However, it is unlikely that all 3000 best NRSF matrix matches are true functional binding sites. Thus, we reduced the positive set to the best 500 motif instances with the highest score. Using this smaller positive set, within 3000 best peaks, MICSA was able to successfully identify 447 out of all 500 instances (Figure 2B). Other tools can be ranked accordingly: uSeq (437), MACS (435), F-Seq (432), CisGenome (431), PeakSeq (431), wdt (430), FindPeaks 3.3 (428), SISSRs (424), ERANGE 3.1 (412) and QuEST 2.0 (402) (Figure 2B, Supplementary Figure S3b and Supplementary Table S3).

When evaluating the tool sensitivity on the set of 83 qPCR verified NRSF-binding sites, comparable results were obtained: MICSA (55 sites), MACS (52), CisGenome (52), FindPeaks 3.3 (51), SISSRs (51), F-Seq (51), uSeq (51), wdt (51), PeakSeq (51), QuEST 2.0 (50), ERANGE 3.1 (50) (Figure 2C, Supplementary Figure S3c and Supplementary Table S4). However, this can be probably explained by the fact that only regions that were highly enriched in mapped reads were tested by qPCR. It is noteworthy that MICSA was able to identify the maximum number of experimentally tested regions within the dataset of 3000 reported peaks.

Motifs identified by MICSA during the analysis procedure nicely corresponded to the known NRSF-binding motifs (Supplementary Figure S4).

In the comparison above MICSA was the only tool in which DNA motif information was utilized. Since the CisGenome package also contains a module for identification of enriched motifs we added to the comparison the subset of peaks identified by CisGenome which also contain strong hits for enriched motifs. Though motifs identified by CisGenome were close to those identified by MICSA, the latter called more peaks from the positive datasets (Supplementary Figure S5).

We run MICSA on datasets for GABP [START_REF] Valouev | Genome-wide analysis of transcription factor binding sites based on ChIP-Seq data[END_REF], STAT1 (9) and CTCF [ENCODE project, the Broad Institute and the Bradley E. Bernstein lab at the Massachusetts General Hospital/Harvard Medical School [START_REF] Birney | Identification and analysis of functional elements in 1% of the human genome by the ENCODE pilot project[END_REF]] and compared the reported sets of peaks with those reported by other peak calling tools: FindPeaks, uSeq, QuEST and PeakSeq. To compare algorithms between each other we applied the same procedure [START_REF] Kharchenko | Design and analysis of ChIP-seq experiments for DNA-binding proteins[END_REF] as previously for NRSF. We used canonical sequence motifs for binding by GABP (Genomatix, http://www.genomatix.de), STAT1 [START_REF] Horvath | A STAT protein domain that determines DNA sequence recognition suggests a novel DNA-binding domain[END_REF] and CTCF [START_REF] Kim | Analysis of the vertebrate insulator protein CTCF-binding sites in the human genome[END_REF] to create positive sets of 3000 peaks. For any considered TF, within a given number of peaks selected by each algorithm there were more peaks from the positive set in the output of MICSA than in the output of any other program (Supplementary Figure S6). For all three TFs, binding motifs identified by the MICSA pipeline were highly similar to the canonical motifs of the same TF (Figure 3).

The results of the test above might appear to be anticipated in advance. Indeed, MICSA is the only algorithm to use motif information in prediction and the motifs identified by MICSA perfectly matched to the known binding motifs for the considered TFs. So we performed an additional test to assess MICSA's performance. We took the latest ChIP-Seq data for NRSF [ENCODE project, Myers Lab at the HudsonAlpha Institute for Biotechnology [START_REF] Birney | Identification and analysis of functional elements in 1% of the human genome by the ENCODE pilot project[END_REF]]. There the sequencing depth was increased to obtain about 13 million of uniquely mapped reads instead of 2 million in the Johnson dataset [START_REF] Johnson | Genome-wide mapping of in vivo protein-DNA interactions[END_REF]. Since the depth of sequencing had increased, we expected peaks selected by MICSA from the old dataset to appear in the new dataset and to have a read coverage greater than peaks selected by FindPeaks from the old dataset but rejected by MICSA. Indeed, there are 4572 peaks among 7780 selected by MICSA which depth increased more than or exactly twice, while there are 11 748 peaks among 22 116 peaks selected by FindPeaks. This constitutes the overall ratio of 59% against 53% (Supplementary Table S5). The difference is higher for low peaks and disappears with the increase of peak heights (Supplementary Figure S7). This test though not being based on the information about a priori known motifs, demonstrated the advantage of using MICSA for filtering low peaks.

Biological application of MICSA

To demonstrate the use of the MICSA package to obtain biological insight, we then applied the package to ChIP-Seq data for the oncogenic TF EWS-FLI1. EWS-FLI1 is a chimeric protein produced by a fusion of the EWS and FLI1 genes. This abnormal TF is a key oncogene in Ewing sarcoma [START_REF] May | Ewing sarcoma 11;22 translocation produces a chimeric transcription factor that requires the DNA-binding domain encoded by FLI1 for transformation[END_REF][START_REF] Riggi | EWS-FLI-1 expression triggers a Ewing's sarcoma initiation program in primary human mesenchymal stem cells[END_REF]. To investigate EWS-FLI1 DNA binding, we reanalyzed ChIP-Seq data previously collected by our team [START_REF] Guillon | The oncogenic EWS-FLI1 protein binds in vivo GGAA microsatellite sequences with potential transcriptional activation function[END_REF]. The DNA fragments retained on complexes immunoprecipitated by an FLI1-specific antibody were processed by the Illumina/ Solexa cluster station and 1G analyzer for the A673 Ewing cancer cell line and aligned with the Maq [START_REF] Li | Mapping short DNA sequencing reads and calling variants using mapping quality scores[END_REF] software with a maximum two mismatches to the unmasked human reference genome (NCBIv36, hg18) [START_REF] Karolchik | The UCSC Genome Browser Database: 2008 update[END_REF][START_REF] Lander | Initial sequencing and analysis of the human genome[END_REF]. A control was obtained using the same anti-FLI1 antibody in a rhabdoid tumor cell line (MON) that does not express EWS-FLI1 or ETS family TFs.

A key characteristic of these data is that the total amount of sequenced DNA is insufficient for straight-forward identification of the majority of binding sites. Indeed, our previous analysis of these data showed a very limited number of regions of EWS-FLI1 specific binding (246). Upon re-analysis with the MICSA tools, we were able to discover 2264 sites with an expectation FDR of 5%. MICSA was also able to identify two known consensus motifs occurring in the most highly enriched regions of called peaks [START_REF] Guillon | The oncogenic EWS-FLI1 protein binds in vivo GGAA microsatellite sequences with potential transcriptional activation function[END_REF]. The first one represents a (GGAA) !6 microsatellite and it is found in 496 peak sequences; the second motif, found in 1768 peak sequences, corresponds to the consensus RCAGGAARY (R = A/G, Y = T/C) (Figure 4A).

Since the single RCAGGAARY consensus motif resembles the known binding motifs for the ETS TF family (26) (Figure 4B), we will refer to it as the ETS motif. Interestingly, the extracted ETS motif, although resembling the consensus motif of FLI1 ( 27) (CCGGAA RY) (Figure 4C), does not completely coincide with the latter, even though EWS-FLI1 shares the same DNA-binding domain as FLI1.

We compared the locations of discovered peaks with gene organization and expression data [START_REF] Tirode | Mesenchymal stem cell features of Ewing tumors[END_REF] for the Ewing cancer cell line A673 in both the presence and absence of EWS-FLI1 using a random set of peaks as a control (Figure 5). To create the random set we randomly selected 2264 locations in the annotated part of the human genome (NCBIv36, hg18) [START_REF] Karolchik | The UCSC Genome Browser Database: 2008 update[END_REF][START_REF] Lander | Initial sequencing and analysis of the human genome[END_REF]. From the expression data, we extracted a list of putative target genes of EWS-FLI1: 557 genes downregulated by EWS-FLI1 and 577 upregulated genes (fold change >|2| with a Welsh P-value <0.01). These are genes modulated by EWS-FLI1 in A673 and SK-N-MC Ewing cell lines.

Our analysis revealed the tendency of sites bearing microsatellites to upregulate neighboring genes [sites found from 150-kb upstream to 50-kb downstream of gene transcription start sites (TSSs)] (Figure 5A and Supplementary Figure S8), while sites with the ETS motif do not seem to have a definite activator function (Figure 5B). ETS sites show some transcriptional inhibitory influences on gene expression when located in the first 50-kb downstream of the TSSs. However, when ETS-sites are found further away from genes (within 1 Mb upstream or downstream but not in the first 50-kb downstream TSS), both activatory and inhibitory influences are observed for EWS-FLI1 transcriptional activity. Among other hypotheses, this could be explained by competitive binding of EWS-FLI1 and native repressor or activator TFs. Our analysis confirmed binding sites for five known direct target genes of EWS-FLI1: C-Myc, CCND1, TGFbRII (29), CAV1 (30) and IGF1 [START_REF] Cironi | IGF1 is a common target gene of Ewing's sarcoma fusion proteins in mesenchymal progenitor cells[END_REF]. However, in two cases out of five, a binding site was identified inside the gene and not in the promoter region. MICSA predicted many new possibly direct targets of EWS-FLI1 (Supplementary Table S6). Among them we find PPP1R1A, LBH, FAS, CAV2 and NBL1. This information will aid in the construction of a detailed and accurate regulation network for this particular type of cancer.

Our motivating idea was that low peaks without motifs are likely to be false positives while low peaks with motifs may indicate real binding events. To test it we performed ChIP-qPCR for 16 interesting genes having a low peak of EWS-FLI1 within 50-kb upstream of gene TSS (Supplementary Table S7). To create a control set we took seven peaks which had been discarded by MICSA as low peaks without strong motif occurrence. The heights of selected peaks vary from 3.9 to 8 for the set of peaks selected by MICSA and from 4 to 8 for the control set.

For tested genes from MICSA's set we got a clear positive response for CCND1, GYG2 and PAPPA, and a positive trend for AKAP7 and SLCO5A1 (Supplementary Figure S9). Interestingly, none of the 7 peaks from the control set was found to be positive in our experiment. Though we performed the experiments on a limited set of peaks, we believe that the results clearly reinforce our idea that MICSA is a useful tool to distinguish between possible binding events in the case of low peaks.

DISCUSSION

To our knowledge, MICSA is the first tool developed for peak identification in ChIP-Seq data that uses an approach combining knowledge about DNA fragment coverage in ChIP and control experiments along with motif discovery. MICSA is able to automatically identify overrepresented motifs in a single run, as well as to use motif occurrence probabilities to enhance the result set returned. MICSA achieves a higher accuracy in . EWS-FLI1 binding to GGAA microsatellites results in significant expression activation of neighboring genes. EWS-FLI1 binding to single ETS sites can produce both negative and positive effects on transcription of neighboring genes. The P-values were directly evaluated by Monte-Carlo simulations of random peaks. Distances from the TSSs of modulated genes to random peaks (iterative trials) and to predicted sites were calculated. The P-values correspond to the probability to get at least the observed number of distances falling within a given 50-kb window, under the hypothesis that peaks are randomly distributed and their coordinates are independent of coordinates of TSSs of EWS-FLI1 modulated genes. Bars above the dashed line correspond to a P-value <0.05.

identifying regions of TF binding in comparison to other methods. For example, no other tested tool was able to identify more than 45% of predicted motif occurrences (Figure 2) within the top 3000 selected peaks. With default parameters, only Useq (11) (with the best 6073 peaks), MACS (8) (with 6450 peaks), FindPeaks (5) (with 7097 peaks) and F-Seq (12) (with 8576 peaks) managed to obtain the same coverage. However, that required a significantly larger set of selected peaks. This example shows that using MICSA helps to avoid inclusion of thousands of low peaks which do not carry TF-binding motifs and thus are likely to be false binding sites.

Two other major advantages of MICSA, in addition to its greater accuracy, are the score calculation for each reported peak and evaluation of the number of false positives in the total output.

MICSA uses two previously published tools: FindPeaks (5) for candidate peak calling and MEME [START_REF] Bailey | Fitting a mixture model by expectation maximization to discover motifs in biopolymers[END_REF] for de novo motif finding. They were chosen as the best performing from the point of view of speed and result quality. Additionally, FindPeaks allows elimination of duplicate tags due to PCR errors and supports many input formats including: both Maq's .map and mapview formats [START_REF] Li | Mapping short DNA sequencing reads and calling variants using mapping quality scores[END_REF], ELAND, ELAND Extended and BED format.

Because of its high sensitivity, MICSA can be used on medium quality datasets with low average DNA fragment coverage. EWS-FLI1 data (16) represents one such example. In spite of non-perfect initial data quality, we were able to get rich biological results and some insights in the function of EWS-FLI1 in Ewing cancer cells.

For the EWS-FLI1 TF we identified two different types of motifs carrying different biological functions, i.e. microsatellites (GGAA) n and single RCAGGAARY motif (R = A/G, Y = T/C). Interestingly, none of them completely coincide with the known binding motif of FLI1, which has the same DNA-binding domain as EWS-FLI1. This may be a result of the presence in EWS-FLI1 of an additional EWS domain. The observed motif difference provides evidence that de novo motif finding is an important issue in ChIP-Seq data analysis. Our results suggest that EWS-FLI1 binding to a site bearing (GGAA) n microsatellites can activate gene expression if the site occurs within 150-kb upstream and 50-kb downstream region from gene TSS. Occasionally, even more distant sites bearing (GGAA) n microsatellites appear to moderate the activator function of EWS-FLI1. EWS-FLI1 binding to a site without microsatellites can, depending on the gene, activate or repress transcription (sites within 1-Mb upstream 1-Mb downstream of gene TSS). This change of regulatory function depending on binding motif provides an insight into the molecular mechanisms of EWS-FLI1 function. One of the hypotheses is a conformation change induced by dimerization of EWS-FLI1 on microsatellites. Verification of this hypothesis was out of the scope of the article, however, if true, this would indicate a potentially new way to target Ewing cancer by disrupting dimerization.

In conclusion, the MICSA package proposes solutions for a great number of problems including peak calling with predefined false positives number, peak score calculation and de novo motif identification and should be a useful tool in ChIP-Seq data analysis of TFBSs.

Figure 1 .

 1 Figure 1. Main steps of the MICSA pipeline.

Figure 2 .

 2 Figure 2. Performance comparison of MICSA with FindPeaks, PeakSeq, QuEST and uSeq. As a positive set of binding sites of NRSF we used (A) 3000 best matches of the canonical NRSF matrix in the human genome, (B) 500 best matches of the canonical NRSF matrix in the human genome, (C) 83 q-PCR verified NRSF-binding sites in the human genome. Peaks extracted by each algorithm were ranked according to in-built scores or P-values. For each number of top peaks the frequency of identified positive sites among them was plotted. 'ToolName^' means that the default parameters of the tool were modified to make it report more peaks.

Figure 3 .

 3 Figure 3. Binding motifs identified by MICSA in ChIP-Seq data for GABP, STAT1 and CTCF resemble canonical motifs. (A) GABP motif logos [Weblogos (32)], canonical motif from (Genomatix, http:// www.genomatix.de), (B) STAT1 motif logos (21), (C) motif logos for CTCF (22).

Figure 4 .

 4 Figure 4. Motifs identified by MICSA in EWS-FLI1 ChIP-Seq data resemble but are not identical to the canonical binding motif of FLI1. (A) Consensus motifs identified by MICSA [Weblogos (32)], (B) canonical motif for ETS family of TFs including the TF FLI1 (26), (C) canonical motif for the TF FLI1 (27).

Figure 5 .

 5 Figure 5. Histogram of distances between predicted/random peaks and genes up/downregulated by EWS-FLI1. (A) Predicted sites containing (GGAA) n microsatellites; (B) ETS sites (site without microsatellites). EWS-FLI1 binding to GGAA microsatellites results in significant expression activation of neighboring genes. EWS-FLI1 binding to single ETS sites can produce both negative and positive effects on transcription of neighboring genes. The P-values were directly evaluated by Monte-Carlo simulations of random peaks. Distances from the TSSs of modulated genes to random peaks (iterative trials) and to predicted sites were calculated. The P-values correspond to the probability to get at least the observed number of distances falling within a given 50-kb window, under the hypothesis that peaks are randomly distributed and their coordinates are independent of coordinates of TSSs of EWS-FLI1 modulated genes. Bars above the dashed line correspond to a P-value <0.05.

  PAGE 3 OF 9 Nucleic Acids Research, 2010, Vol. 38, No. 11 e126 amount of time required, we only process three hundred peaks with high DNA fragment coverage. Moreover, the motif extraction is done only for central area of these peaks to avoid getting false motifs that are over-represented, such as those of cooperating TFs. The MICSA package utilizes the expectation maximization MEME algorithm (19) to detect the strongest motifs in peak sequences.

ACKNOWLEDGEMENTS

The authors thank A. Zinovyev, P. Kharchenko, I.V. Kulakovskiy and N. Rajewsky for valuable discussion and help, M. Hue for the idea of the optimization procedure and K. Bleakley for the proofreading of the manuscript.

FUNDING

Curie Institute, the « Ligue Nationale contre le Cancer » (V.B., E.B., D.S., N.G., F.T. and O.D. are members of labeled team, the project was also supported by the CIT program 'Carte d'Identite´des Tumeurs'); Institut National de la Sante´et de la Recherche Me´dicale and the Agence Nationale de la Recherche (SITCON project). Funding for open access charge: Institut Curie, Paris, France.

SUPPLEMENTARY DATA

Supplementary Data are available at NAR Online.

Conflict of interest statement. None declared.


Florian Dupuy

Thierry Hedde

Pierre Roubin

Pierre Durand

Statistical downscaling of WRF mesoscale simulations using an artificial neural network to forecast winds at local, unresolved scale

WRF mesoscale (3km horizontal resolution) forecasts and local observations are available hourly on a 1 year period Use of an Artificial Neural Network (ANN) to generate a statistical forecast of the wind at 2m

The method is assessed by using the following indexes, in complement to usual statistical indexes:

) Winds lower than 0.50 𝑚/𝑠 are not considered for the calculation of performances because of the lack of precision on their direction.

II. Methodology

The Artificial Neural Network

• is a statistical tool commonly used for atmospheric studies (Gardner and Dorling, 1998), only with observed data up to now • is particularly suitable for non linear problems • searches a relation between some input variables (WRF variables) and the wind components observed at 2m

Several WRF variables (linked with the stability, the daily cycle, the synoptic forcing and the wind) are used as inputs to the ANN  The 4 most suitable variables are:

• the wind components at 10m 𝑢 10 and 𝑣 10 • a difference of potential temperature ∆𝜃 (110𝑚 -2𝑚)

• the wind speed at 110m 𝑈 110 .

III. Results

Motivation

Calculate the atmospheric dispersion of pollutants emitted in an industrial center (Cadarache), using daily forecasts made with the WRF model at a 3 km horizontal resolution.

Context

• Complex topography (situated in a small valley)  has effects on the low level transport • High occurrence of clear skies  favors stable conditions that decreases vertical mixing • The KASCADE experiment (Duine, 2015) highlights the importance of such conditions that generate a Cadarache down-valley wind (CDV) At a 3 km resolution, these effects cannot be simulated, which may invalidate the prediction of atmospheric dispersion in low wind stable conditions.

I. Introduction

• The statistical forecast of a valley wind based on using an ANN in conjunction with mesoscale simulations is successful • The ANN performance remains dependent on the precision of the mesoscale simulations (for the input variables)

• In spite of the large size of the dataset, the ANN failed to reproduce wind directions falling in the red sector (figure 2) • The ANN needs both wind and stability data to calculate the valley wind because they often have a thermal origin

V. Conclusion

Relative role of each input variable

• 𝑢 10 and 𝑣 10 are fundamental for the calculation of outputs • ∆𝜃 (110𝑚 -2𝑚) improves the calculation of down-valley thermal winds

• The role of 𝑈 110 is not well defined The performance is dependent on the speed, increasing with the wind speed Better predictions of the local valley wind are obtained when using the same input variables with their observed values, highlighting the effect of the quality of input variables on the performance of the ANN.

IV. Discussion

European

  Meteorological Society annual meeting -16 September 2016 -Trieste, Italy Duine, G.-J. (2015). Characterization of down-valley winds in stable stratification from KASCADE field campaign and WRF mesoscale simulations. Theses, Université Toulouse III Paul Sabatier Gardner, M. and Dorling, S. (1998). Artificial neural networks (the multilayer perceptron)-a review of applications in the atmospheric sciences. Atmospheric environment, 32(14) :2627-2636.

Fig 5 :

 5 Fig 5: DACC values calculated for different wind speed classes (blue columns). Number of data in each class (black line).

Fig 2 :

 2 Fig 2: Wind rose for the observations at 2m and separation of the wind into three classes (blue for down-valley winds, green for up-valley and red for others).

Fig 1 :

 1 Fig 1: Comparison of the 2m wind observed and forecasted, represented by the probability density function in color, for the direction (a) and the speed (b). The colored frames refer to the wind sectors of figure 2.

Fig 3 :

 3 Fig 3: Same as Fig 1 for the comparison of the 2m wind observed and calculated by the ANN.