Jean-Marie Poumirol

Ioannis Paradisanos

Shivangi Shree

Gonzague Agez

Xavier Marie

Cédric Robert

Nicolas Mallet

Peter R Wiecha

Guilhem Larrieu

Vincent Larrey

Frank Fournel

Kenji Watanabe

Takashi Taniguchi

Aurélien Cuche

Vincent Paillard
email: vincent.paillard@cemes.fr

Bernhard Urbaszek
email: urbaszek@insa-toulouse.fr

Unveiling the optical emission channels of monolayer semiconductors coupled to silicon nanoantennas

Monolayers (MLs) of transition metal dichalcogenides (TMDs) such as WSe2 and MoSe2 can be placed by dry stamping directly on broadband dielectric resonators, which have the ability to enhance the spontaneous emission rate and brightness of solid-state emitters at room temperature. We show strongly enhanced emission and directivity modifications in room temperature photoluminescence mapping experiments. By varying TMD material (WSe2 versus MoSe2) transferred on silicon nanoresonators with various designs (planarized versus non-planarized), we experimentally separate the different physical mechanisms that govern the global light emission enhancement. For WSe2 and MoSe2 we address the effects of Mie Resonances and strain in the monolayer. For WSe2 an important additional contribution comes from out-of-plane exciton dipoles. This paves the way for more targeted designs of TMD -Si nanoresonator structures for room temperature applications.

Introduction.-Integrated efficient quantum light sources need to be developed for different applications including quantum computing and transmitted or stored information for cryptography. Key technologies are plasmonic or dielectric nanostructures (NS) coupled to various quantum emitters positioned in the NS near-field. Both the emitting dipole position and orientation with respect to the nanoantenna are important to optimize the emission rate and the emitted light directivity and waveguiding. Therefore, monolayers of semiconducting transition metal dichalcogenides (TMDs) [START_REF] Mak | Photonics and optoelectronics of 2D semiconductor transition metal dichalcogenides[END_REF][START_REF] Wang | Colloquium: Excitons in atomically thin transition metal dichalcogenides[END_REF] such as WSe 2 and also heterobilayers [START_REF] Nagler | Interlayer Excitons in Transition-Metal Dichalcogenide Heterobilayers[END_REF] are very promising as an active medium. Their atomic-monolayer thickness makes them uniformly sensitive to the near field around a nanoantenna, and they host several kinds of emitting dipoles at different wavelengths, such as the bright exciton (X o , dipole in the monolayer plane), the so-called dark exciton (X D , dipole out-of-plane) and other identified complexes such as charged excitons, localized excitons and biexcitons [START_REF] Zhou | Probing dark excitons in atomically thin semiconductors via near-field coupling to surface plasmon polaritons[END_REF][START_REF] Wang | In-Plane Propagation of Light in Transition Metal Dichalcogenide Monolayers: Optical Selection Rules[END_REF][START_REF] Barbone | Charge-tuneable biexciton complexes in monolayer WSe 2[END_REF]. TMDs coupled to plasmonic NS have shown strongly enhanced emission and tailored directivity [START_REF] Wang | Giant photoluminescence enhancement in tungstendiselenide-gold plasmonic hybrid structures[END_REF][START_REF] Kern | Nanoantenna-enhanced lightmatter interaction in atomically thin WS2[END_REF][START_REF] Najmaei | Plasmonic pumping of excitonic photoluminescence in hybrid MoS2-Au nanostructures[END_REF]. However, plasmonics is limited by high losses in metals, heating and incompatibility with complementary-metal-oxide-semiconductor (CMOS) fabrication processes [START_REF] Won | Into the 'Mie-tronic' era[END_REF]. High refractive index dielectric NS such as silicon-NS or GaP-NS [START_REF] Wiecha | Evolutionary multi-objective optimization of colour pixels based on dielectric nanoantennas[END_REF][START_REF] Cihan | Silicon Mie resonators for highly directional light emission from monolayer MoS 2[END_REF][START_REF] Sortino | Enhanced light-matter interaction in an atomically thin semiconductor coupled with dielectric nano-antennas[END_REF], can also provide reduced mode volumes, but in contrast to plasmonics they show low non-radiative losses and are compatible with large-scale semiconductor fabrication processes [START_REF] Bidault | Dielectric nanoantennas to manipulate solid-state light emission[END_REF][START_REF] Kuznetsov | Optically resonant dielectric nanostructures[END_REF]. The emission properties of excitons in TMDs placed on NS have been investigated for several combinations of NS geometry and TMD materials [START_REF] Cihan | Silicon Mie resonators for highly directional light emission from monolayer MoS 2[END_REF][START_REF] Sortino | Enhanced light-matter interaction in an atomically thin semiconductor coupled with dielectric nano-antennas[END_REF][START_REF] Zhang | Unidirectional doubly enhanced MoS2 emission via photonic Fano resonances[END_REF][START_REF] Chen | Enhanced directional emission from monolayer WSe2 integrated onto a multiresonant silicon-based photonic structure[END_REF][START_REF] Lee | Single type of nanocavity structure enhances light outcouplings from various two-dimensional materials by over 100-fold[END_REF][START_REF] Duong | Enhanced Emission from WSe2 Monolayers Coupled to Circular Bragg Gratings[END_REF][START_REF] Bucher | Tailoring photoluminescence from MoS2 monolayers by Mie-resonant metasurfaces[END_REF][START_REF] Mupparapu | Integration of twodimensional transition metal dichalcogenides with Mieresonant dielectric nanostructures[END_REF]. Despite these promising first results, the physical mechanisms that govern changes in the TMD emission as they are placed in proximity of the nanoresonators need to be understood. In addition to the targeted modification of the local electric field due to the resonator, several other mechanisms such as strong local strain induced by the resonator [START_REF] Sortino | Dielectric Nanoantennas for Strain Engineering in Atomically Thin Two-Dimensional Semiconductors[END_REF] and changing optical dipole orientation with respect to the detection contribute to brightening and change in directivity of the emission.

In this work, we experimentally separate the three main mechanisms that govern the photoluminescence emission properties of the coupled dielectric NS-TMD monolayer system at room temperature, namely (i) interactions with broad-band Mie resonances, (ii) strain in the TMD layer and (iii) orientation of the exciton dipole with respect to the detection. First, we address the impact of the dipole orientation by comparing on identical non-planarised resonator samples (see dimer resonators sketched in Fig. 1a) the emission of WSe 2 (in-plane and out-of-plane dipoles) [START_REF] Zhou | Probing dark excitons in atomically thin semiconductors via near-field coupling to surface plasmon polaritons[END_REF][START_REF] Wang | In-Plane Propagation of Light in Transition Metal Dichalcogenide Monolayers: Optical Selection Rules[END_REF][START_REF] Park | Radiative control of dark excitons at room temperature by nano-optical antenna-tip Purcell effect[END_REF] with MoSe 2 (only in-plane emission) [START_REF] Lu | Magnetic field mixing and splitting of bright and dark excitons in monolayer MoSe2[END_REF][START_REF] Robert | Measurement of the Spin-Forbidden Dark Excitons in MoS2 and MoSe2 monolayers[END_REF]. For WSe 2 we show a brightening of the global emission by one order of magnitude on the NS. This is accompanied by a change in overall spectral shape, due to a variation of the relative contributions to the emission from the in-plane and out-of-plane dipoles, which emit at slightly different energies. In the PL mapping, we find a strong dependence of the spectral shape on detection spot position on the NS (edge/center/suspended region). Second, we address the impact of strain as we com-pare TMD monolayer emission from samples on Si nanostructures (height about 100 nm) with samples on the same Si nanostructures surrounded by silica (planarized samples), where the TMD is strain free, compare Figs. 1c ande. We measure in non-planarized samples a clear redshift of the PL energy (correlated to the TMD phonon shift measured in Raman mapping), showing the contribution to enhanced emission of the tensile strain-induced band gap lowering. Strained MoSe 2 behaves differently compared to WSe 2 : although we see a PL energy redshift due to strain, the spectra consist only of a main, single peak corresponding to the in-plane dipole in MoSe 2 (bright exciton). On planarized samples we observe WSe 2 and MoSe 2 emission with less PL intensity enhancement (50%) as compared to the strained case, but with a clear impact of NS geometry on the PL emission pattern and dependence on the excitation laser polarization. Our experiments show an enhanced brightness for all NS-TMD configurations as compared to the TMD alone, so the Si-NS is beneficial in all cases. Samples and optical spectroscopy techniques.-In this section we describe the two different kinds of nanoresonators used, which show broad resonances covering the wavelength region that corresponds to the MoSe 2 and WSe 2 monolayer emission and excitation.

Si-NS are obtained by a top-down approach in a Silicon on Insulator (SOI) wafer, where a hydrogen silsesquioxane (HSQ) resist is patterned by electron beam lithography followed by reactive ion etching. The thickness of the single crystal overlayer which determines the Si-NS height is H =90 nm. After etching and annealing, a SiO 2 layer (from the HSQ deposition) of about 20 nm is left on top of the Si-NS. We will focus on dimers consisting of two nanorods of rectangular cross section having dimensions of L × W × H and separated by a gap G (Fig. 1c). Values of L and G are varying across the resonator array from 50 to 300 nm, for constant W =150nm and H =90 nm. An array of fabricated structures can be seen in Fig. 1b. Before TMD monolayer transfer, we perform dark-field scattering experiments to probe the optical Mie resonance bandwidth of the different dimers [START_REF] Won | Into the 'Mie-tronic' era[END_REF][START_REF] Bidault | Dielectric nanoantennas to manipulate solid-state light emission[END_REF][START_REF] Kuznetsov | Optically resonant dielectric nanostructures[END_REF]. A typical dark-field spectrum is shown in Fig. 1h, where the measured resonances extend up to the relevant exciton transitions in the TMDs. As it can already be seen in Fig. 1g by the different colors, the characteristics of the Mie resonances can be varied by changing L, W and G [START_REF] Wiecha | Evolutionary multi-objective optimization of colour pixels based on dielectric nanoantennas[END_REF].

Flakes from TMD bulk crystals are obtained by multiple folding in adhesive tape. The adhesive tape with the flakes of various thicknesses is pressed against a transparent polymer stamp, here PDMS, to transfer the flakes. Monolayers on the PDMS stamp are identified by optical contrast (confirmed by AFM). As a last step [START_REF] Castellanos-Gomez | Deterministic transfer of two-dimensional materials by all-dry viscoelastic stamping[END_REF], the desired WSe 2 or MoSe 2 monolayer is placed above the targeted nano-resonator array and then transferred on the patterned substrate, covering 30 Si-NS with the same flake in Fig. 1b. An atomic force microscope (AFM) profile in Fig. 1d shows that the WSe 2 ML is acting as a nearly conformal coating, following the Si-NS morphology. These samples are referred to as non-planarized and will be used to investigate the impact of strain and dipole orientation on the PL emission (Fig. 1c).

In order to eliminate strain effects and focus solely on the impact of the Si antenna on the TMD emission, another set of samples was fabricated. We followed the same fabrication procedure as for the non-planarized Si-NS and added a planarization step involving HSQ deposition and subsequent annealing and chemical etching [START_REF] Guerfi | Thin-dielectriclayer engineering for 3D nanostructure integration using an innovative planarization approach[END_REF]. This leads to Si nanoantennas buried in silica and we refer to this as planarized samples, see Fig. 1e. The distance between the top Si-NS and the silica smooth surface (see AFM profile Fig. 1f) is about 30 nm over 1 µm. Although the surface is not totally flat, the height variations are smooth enough to avoid strain in the TMD, as confirmed in PL and Raman mappings (see below). The active region is here protected by a 10 nm-thick hexagonal BN layer (high crystal quality, atomically flat [START_REF] Taniguchi | Synthesis of high-purity boron nitride single crystals under high pressure by using Ba-BN solvent[END_REF]) on top of the TMD ML.

Room temperature photoluminescence and Raman experiments are performed using either a Horiba JY Labram HR or a Horiba JY XPlora MV2000 set-up. Hyper-spectral maps are obtained by raster scanning the sample under a laser spot (532 nm excitation wavelength, power of tens of nW) tightly focused through a microscope objective (0.9 numerical aperture -NA). The spectra are obtained in backscattering using the same objective. The high NA allows us to collect mainly emission from in-plane but also to a small extent from out-of plane optical dipoles [START_REF] Wang | In-Plane Propagation of Light in Transition Metal Dichalcogenide Monolayers: Optical Selection Rules[END_REF][START_REF] Brotons-Gisbert | Out-of-plane orientation of luminescent excitons in two-dimensional indium selenide[END_REF]. Although the diffraction limited spot has a diameter comparable to the dimer, the recorded PL emission map of the sample reveals fine details in lateral resolution due to the high signal to noise ratio. The laser polarization is kept constant and the so-called x (linear polarization along the dimer axis) and y polarization (linear polarization perpendicular the dimer axis) are obtained by rotating the sample by 90 • . Here we define as the dimer axis a straight line across the dimer gap that connects the centers of the two Si rectangles.

WSe 2 monolayer on non-planarized nanoresonators.-In Fig. 2a we plot a normalized, spectrally integrated PL map of non-planarized nanoresonators covered by one large WSe 2 monolayer. The nanoresonator positions can be seen directly on the map, PL away from the resonators is very weak (uniform dark blue area). To allow for a direct comparison of the results in Figs. 2, 3, 4 we show spectra from the Si dimer nanoresonator with the same dimensions in all cases of gap G = 300 nm and length L = 300 nm (white circle in Figs. 2a andg, 3a). This dimer has broad optical Mie resonances at the emission as well as the laser excitation energy (See Fig. 1h). Strain will clearly play an important role for non-planarized nanoresonators and similar structures are used to generate localized excitons (0D) at low temperature for quantum optics [START_REF] Branny | Deter-ministic strain-induced arrays of quantum emitters in a two-dimensional semiconductor[END_REF][START_REF] Palacios-Berraquero | Large-scale quantum-emitter arrays in atomically thin semiconductors[END_REF]. Note that we avoid exciton localization in our room temperature measurements as the thermal energy of excitons is larger than the localization energy.

In our experiment in Fig. 2a the excitation laser is polarized along the main axis of the dimer (x), Fig. 2g shows identical measurements, but carried out with laser polarization perpendicular to the dimer (y) (white arrows indicating the electric field in the figures).

We first comment on x-polarized excitation. The WSe 2 PL intensity is substantially enhanced by one order of magnitude on top of the nanoresonators. In addition to intensity changes, we clearly observe that the spectra recorded on top and next to the nanoresonator shown in Fig. 2b are very different, compare panels Fig. 2c,d ande. The PL emission is fitted using three peaks, in agreement with the different transitions identified in the literature for WSe 2 ML even at room temperature : the free exciton at about 1.66 eV (745 nm) (dashed blue line), the charged excitons (trions) (dashed green line) at about 1.64 eV (756 nm), and the dark (out-of-plane) exciton near 1.61 eV (770 nm) (dashed red line) [START_REF] Zhou | Probing dark excitons in atomically thin semiconductors via near-field coupling to surface plasmon polaritons[END_REF][START_REF] Wang | In-Plane Propagation of Light in Transition Metal Dichalcogenide Monolayers: Optical Selection Rules[END_REF][START_REF] Park | Radiative control of dark excitons at room temperature by nano-optical antenna-tip Purcell effect[END_REF]. Please note that we have identified the respective dipole orientation of the dark and bright excitons previously by imaging the Fourier plane of the microscope objective in experiments on ML WSe 2 [START_REF] Wang | In-Plane Propagation of Light in Transition Metal Dichalcogenide Monolayers: Optical Selection Rules[END_REF]. The relative weight of the different peaks strongly varies as a function of the position along the map, especially around and above the Si-NS (see Fig. 2f top panel).

The impact of strain on the optical transitions is revealed in Fig. 2f bottom panel, where we plot the emission energy of the three PL components as a function of position. A shift of the PL emission to lower energy for all three transitions is a clear confirmation for tensile strain in the layer [START_REF] Schmidt | Reversible uniaxial strain tuning in atomically thin WSe2[END_REF]. The same conclusion can be independently drawn from the Raman maps, presented in the supplement. The enhancement observed for emission on the nanoresonators benefits therefore from the combined effect of strain (local lowering of the band gap that results in funneling of excitons) and optical resonances, which we aim to separate in the next section. Interestingly, as one can see in Fig. 2 (f) top panel, the emission linked to the dark exciton (red line), which has an optical dipole out-of the monolayer plane, shows strongest enhancement on the sample edges. Here we suggest that, as the monolayer is folded with respect to the microscope objective axis (by 30 • to 40 • , as given by AFM profile), more light emitted by the out-of-plane dipole is detected. These arguments are developed in more detail in the discussion section. In addition, enhanced emission from excitons with out-of-plane dipoles might be due to lowering of the crystal symmetry [START_REF] Branny | Deter-ministic strain-induced arrays of quantum emitters in a two-dimensional semiconductor[END_REF][START_REF] Luo | Exciton dipole orientation of strain-induced quantum emitters in WSe2[END_REF]. The impact of strain for WSe 2 MLs on GaP dimers is discussed in references [START_REF] Sortino | Enhanced light-matter interaction in an atomically thin semiconductor coupled with dielectric nano-antennas[END_REF][START_REF] Sortino | Dielectric Nanoantennas for Strain Engineering in Atomically Thin Two-Dimensional Semiconductors[END_REF].

In Fig. 2g we plot the same intensity map as in Fig. 2a, but the excitation laser polarization is now perpendicular to the dimer axis. We observe clear differences: the global enhancement on the NS is now a factor of 3, as compared to roughly 8 for the other laser polarization, and the spatial region around the NS position that shows enhancement is smaller. The maximum enhancement is localized in the gap between the two rectangular rods constituting the dimer, as can be seen in Fig. 2h. In the dimer gap the monolayer is suspended, and therefore presumably less strained than on the dimer edges. As a result, the PL emission energy does not change as a function of position as we scan across the gap, see Fig. 2l bottom panel.

Our measurements also show, in contrast to laser excitation parallel to the dimer axis, the global shape of the PL spectra is not affected in the perpendicular case (Fig. 2i,j,k). The relative amplitude of the different PL components remains constant when scanning over an antenna (Fig. 2l top panel). We suggest in the discussion section below, that this striking dependence on laser excitation polarization (compare Figs. 2a,b with Figs. 2g,h) is linked to an optical antenna effect.

MoSe 2 monolayer on non-planarized nanoresonators.-For WSe 2 we observe on several Si-NS dimers a strong PL intensity enhancement above each element of the dimer, showing amplification as high as 8 as compared to the bare monolayer. We now investigate emission from a MoSe 2 monolayer deposited on an identical set of resonators, see Fig. 3a. We choose this material for comparison as the emission wavelength is close to WSe 2 and can be covered by the broadband resonances of the Si-NS as shown in Fig. 1. An additional advantage for this material is the comparatively simple emission spectrum as the dark, z-polarized states (optical dipole out-of-plane) are in energy above the bright states [START_REF] Lu | Magnetic field mixing and splitting of bright and dark excitons in monolayer MoSe2[END_REF][START_REF] Robert | Measurement of the Spin-Forbidden Dark Excitons in MoS2 and MoSe2 monolayers[END_REF] and therefore have only very small impact on PL emission here, see spectra in Fig. 3d and scheme in Fig. 5a.

In Fig. 3a we plot the PL emission map above the resonators and we see strong PL emission enhancemement (up to 14 ×) on several resonators covered by the ML flake. PL emission next to the resonators is weak (uniform dark area in the map). The most striking observations are: First, a strong enhancement of the PL emission intensity above the Si-NS, see Fig. 3b. Second, in contrast to our WSe 2 results, the PL emission for MoSe 2 does not change shape, as can be seen in Fig. 3d. The PL band could be fitted using a single peak attributed to the neutral bright exciton. Thus, the red shift of the PL energy shown in Fig. 3c is deduced from the maximum of the measured PL peak. As in the case of WSe 2 , this is a signature of tensile strain but resulting in a smaller energy shift, in qualitative agreement with previous studies on strain tuning [START_REF] Island | Precise and reversible band gap tuning in single-layer MoSe 2 by uniaxial strain[END_REF]. Although due to the manual stamping method the coupling efficiency of the TMD monolayer might not be uniform across the nanoresonator array, we observe for the investigated MoSe 2 monolayer very similar enhancement for several nanoresonators analyzed, as can be seen in the map in Fig. 3a. Interestingly for dimers with narrower gaps two emission lobes appear, which might be due to an antenna effect that redirects the light at the dimer edge towards our microscope objective. The simple PL shape makes MoSe 2 a suitable material for analyzing more complex effects related to Mie resonances, that might be obscured by the multiple transitions in WSe 2 . Please note that MoSe 2 monolayers have in general comparatively weaker emission than WSe 2 monolayers at room temperature [START_REF] Wang | Spin-orbit engineering in transition metal dichalcogenide alloy monolayers[END_REF][START_REF] Withers | WSe2 light-emitting tunneling transistors with enhanced brightness at room temperature[END_REF][START_REF] Zhang | Experimental evidence for dark excitons in monolayer WSe 2[END_REF] and hence a wider margin for increasing emission efficiency.

Concerning the dependence on laser excitation polarization for MoSe 2 we observe a stronger emission for excitation along the dimer axis than perpendicular to it. The difference is about 15 %, see supplement, much less pronounced than for the WSe 2 case, which suggests that the contributions from out-of-plane dipole emission plays an important role for the strong polarization difference observed in WSe 2 . Note that we analyze for both materials the same nanoresonators with the same optical resonances (as shown in Fig. 1h), so major differences in our comparative studies between MoSe 2 and WSe 2 emission have their origin in the different bandstructure, sketched in Fig. 5a.

Strain free WSe 2 monolayer on planarized nanoresonators.-In the case of the planarized samples, neither strain-induced Raman shifts of the TMD phonon can be detected (see supplement) nor a systematic redshift of the PL energy, as shown in our measurements in Fig. 4c andf. The AFM profile in Fig. 1f shows that the maximum roughness above a Si-NS is limited to about 30 nm extended over 1 µm. For this strain-free WSe 2 monolayer on planarized resonators we make the following observations: We record an increase of the PL of 50% due to the presence of Si-NS, this is less than the factor of about 8 achieved on the same nanoresonator but non-planarized, compare Fig. 2a and Fig. 4a. Interestingly, the overall emission spectrum for the planarized sample of WSe 2 does not change shape but only the amplitude is modified by the antenna, see Fig. 4b ande. The bright exciton contribution in the PL band is the strongest peak according to our fits and the weight of the different contributions (dark and charged excitons) are constant when we scan across the Si-NS. The PL maps in Fig. 4 are dependent on the laser polarization for WSe 2 , compare panels Fig. 4a,c with Fig. 4d,f. For y-polarized excitation we see a maximum in PL signal over the dimer gap, see Fig. 4d,f. Our measurements show that the emission of WSe 2 is governed by an optical antenna effect in the planarized samples, whereas strain effects play a key role in the nonplanar samples, as pointed out in earlier work on strain WSe 2 on GaP dimers [START_REF] Sortino | Enhanced light-matter interaction in an atomically thin semiconductor coupled with dielectric nano-antennas[END_REF]. The measured PL enhancement is up to 50% of the integrated intensity, which can be further improved. Our measurements on planarized MoSe 2 samples gave very similar results.

In future experiments the Si-NS shape and dimensions can be optimized in order to reach higher electric field amplitudes. It is important to take into account that, after HSQ planarization and hBN deposition, the dielectric contrast between Si-NS and the surrounding is lower than for non-planarized resonators. This lowers light confinement and needs to be compensated by optimizing the antenna design.

Numerical simulations and discussion.-Our measurements of the PL emission of ML WSe 2 and MoSe 2 coupled to non-planarized and planarized nanoresonators revealed for all cases an enhancement of the PL emission as compared to the bare monolayer. In the discussion below we suggest different mechanisms at the origin of this enhancement taking into account the polarization-dependent near field at the excitation energy and the interaction of dipole in-plane and out-of the ML plane with the Mie resonances of the dimer structures at the emission energy.

The Si-NS is optically resonant at the excitation wavelength. We first simulate the spatial distribution of the electric near-field intensity in the plane of the TMD monolayer (see Fig. 5c) with the finite-difference time-domain (FDTD) method [START_REF] Taflove | Computational electrodynamics: the finite-difference time-domain method[END_REF], using an open-source software package [START_REF] Oskooi | A flexible free-software package for electromagnetic simulations by the FDTD method[END_REF], see results in Fig. 5b. For laser excitation parallel to the dimer axis our simulation shows a stronger near-field enhancement than in the case where the polarization is perpendicular to the dimer axis. The near-field enhancement is localized above and around the Si blocks in the parallel polarization case. This can be seen in particular along the edges (y-oriented) perpendicular to the polarization. In contrast, there is almost no enhancement for the perpendicular case (even lower field above each block of the Si dimer). This behavior, relevant at the laser excitation wavelength, is in qualitative agreement with the PL maps of the chosen example described previously, compare Fig. 2a and Fig. 2g. This indicates that optical Mie resonances contribute to PL enhancement.

Nevertheless Mie resonances alone cannot explain the differences in the enhancement and the polarization dependence observed in the four configurations (ML WSe 2 and ML MoSe 2 on non-planarized and planarized Si-NS, respectively), otherwise the observations for MoSe 2 and WSe 2 for the polarization dependence would have been identical.

To go further, we study the antenna effect induced by Si-NS on the emission of the two kinds of exciton dipoles (the bright exciton X o and the dark exciton X D ) [START_REF] Zhou | Probing dark excitons in atomically thin semiconductors via near-field coupling to surface plasmon polaritons[END_REF][START_REF] Wang | In-Plane Propagation of Light in Transition Metal Dichalcogenide Monolayers: Optical Selection Rules[END_REF][START_REF] Park | Radiative control of dark excitons at room temperature by nano-optical antenna-tip Purcell effect[END_REF]. In a first approximation we consider a single oscillating dipole with orthogonal orientation depending on the nature of the exciton (in-plane X o and out-of-plane X D -see Fig. 5a). Considering the non-planarized and planarized geometries (Fig. 5c ande), we compute the intensity of the electric field collected in backscattering geometry through the microscope objective NA for several key locations of the emitting dipoles. The electric field radiated by the dipolar source in presence of a Si-NS is obtained by the Green Dyadic Method [START_REF] Girard | Near fields in nanostructures[END_REF] using the pyGDM toolkit for electrodynamics simulations [START_REF] Wiecha | pyGDM-A python toolkit for fullfield electro-dynamical simulations and evolutionary optimization of nanostructures[END_REF][START_REF] Wiecha | Enhancement of electric and magnetic dipole transition of rare-earth-doped thin films tailored by high-index dielectric nanostructures[END_REF]. In Fig. 5d-f are shown the emission patterns of different dipoles emitting at 750 nm (approximately the average of the bright and dark exciton wavelengths). Depending on the position around the Si-NS, the dipole represents either a dark exciton (out-of-plane dipole, black arrow) or a neutral exciton (in-plane, red arrow). In the upper left part of Fig. 5d, we represent a dipole marked by a black arrow which is located at one Si-NS edge, at mid-height (z= 55 nm) and oriented along the x-axis of the dimer. Now we need to take into account that our monolayer changes orientation : as represented in Fig 5c, the TMD ML is nearly conformal to the Si-NS i.e. following its shape. That means, at the outer dimer edges a dipole oriented along the x-axis of the dimer is oriented out of the folded monolayer plane. Hence it can be considered as a dark exciton X D with a dipole out of the TMD ML plane. For this dipole our calculations show that the intensity measured through the microscope objective NA (represented by the dotted circle) is very high. In contrast, the collected intensity originating from the dipole in the TMD plane at the same position (represented by a red arrow in Fig. 5c) is negligible. Indeed, the emission pattern (not shown) is oriented along the x-axis due to the combined effects of the dipole orientation parallel to the microscope axis and the Si nanoantenna. These simulations give an indication why the emission signal detected at the X D energy is particularly strong at the dimer edges.

In contrast, if the dipole is positioned at mid-height of the long side of the nanoblock (upper right part of Fig. 5d), the collected intensity drops. We believe that this antenna effect plays an important role in the increase of the measured signal. It adds to the polarizationdependent enhancement presented in Fig. 5b, which leads to the creation of more excitons (including X D ones) in the regions of high electric field. Thus, the combination of these two fundamental mechanisms leads to the boost in the PL enhancement observed for WSe 2 on nonplanar Si-NS, only for a parallel polarization of the laser with respect to the Si dimer axis.

In the gap of the Si dimer (lower part of Fig. 5d), the integrated intensity corresponding to X D emission is approximately the same as the one due to X o . Both are contributing, thus leading again to the signal increase. Compared to the emission pattern of the dipole in a TMD monolayer alone hosting X D (WSe 2 ) (See Fig. 5a), it is clear that the relative contributions of X D and X o are very different as a function of the emitting dipole position in the case of the non-planar sample.

Our simple model is already sufficient to understand the underlying physics leading to the strong PL enhancement when the TMD hosts out-of-plane excitons. For instance, since this cumulative effect mainly occurs along the y-oriented edges (small side of a nanoblock) for a parallel polarization, it helps explaining the extension of the PL patterns along the x-axis well beyond the dimer limits in the maps of Fig. 2. For the same polarization configuration, it is also in good agreement with the PL band shape modification related to the strong enhancement of the fitted (red) peak associated to X D at the y-oriented (outer) edges of the Si-NS in the spectra of Fig. 2.

On the other hand, for the planarized sample, the dark exciton contribution is always very minor and the PL band is dominated by the bright exciton contribution, compare Fig. 4e and Fig. 2e. The radiation pattern of the X D is weakly redirected by the antenna in the angular window defined by the microscope objective NA, neither when it is positioned along the edges, nor in the gap region (Fig. 5f). As a consequence, no enhancement associated to the X D is expected. It is in good agreement with the experimental observation where this dark exciton is intrinsically eliminated (MoSe 2 samples) or does experience no enhancement (WSe 2 on planarized samples). In both cases, the polarization-dependent enhancement, resulting from a pure optical effect associated to the Mie resonances in these Si-NS, is about 50%, far from the roughly 10× enhancement in PL from WSe 2 on nonplanar resonators.

We unveil here three intertwined mechanisms that contribute to the PL enhancement. The Mie-assisted enhancement of the dark exciton contribution plays therefore an important role in the global PL response of the WSe 2 -resonator coupled systems. Interestingly, its contribution seems to be of the same order of magnitude as the strain one. The amplitude of the latter can be estimated with MoSe 2 on non-planar samples. Indeed, about 15 % variation of the PL is observed as the polarization is rotated, indicating a very strong contribution to the global factor PL enhancement (14 ×) originates from the strain generated by the Si nanoblocks in this sample.

Finally, as expected for non-optimized resonators, the pure optical effect associated to Mie resonances observed in planarized samples in Fig. 4 (no strain, no dark exciton) is relatively weak compared to the two other mechanisms. It should easily be further improved by optimizing the geometry of the nanoresonators.

In conclusion, we have shown strong PL enhancement for WSe 2 and MoSe 2 on non-planarized nanoresonators due to the combined effect of amplifying contributions from out-of plane dipole emission, strain and pure electric field amplification. On planarized nanoresonators we show that the TMD monolayer is strain free and the PL emission is only enhanced by the Mie resonances of the Si dimer. In addition to presence and absence of strain in the active medium for non-planarized and planarized samples, respectively, also the dielectric contrast plays a role. For non-planarized resonators dielectric contrast is high, whereas it is lower for planarized samples, which demands further work on planarization techniques such as the use of low refractive index dielectric around the silicon nanostructures, and a very small thickness to position the TMD as close as possible to the Si nanoantenna. Finally the optimization of the nanoantenna design will be essential to control strain, resonances at both the excitation and emission wavelength, and directivity of the emitted light. WSe 2 is a very promising material as the in-plane and out-ofplane exciton emission can be controlled independently to increase the global emission yield. MoSe 2 shows also strong PL emission enhancement and a simple emission spectrum for targeted interaction with Mie resonances. . Bottom: orientation of the in-plane neutral exciton Xo (red arrow) and the out-of-plane dark exciton XD (blue arrow) in the x-z plane with respect to the TMD layer. The gray shaded areas represent the emission lobes of the dipoles. b) Near field intensity in the TMD layer of the non-planarized system for two polarizations parallel and perpendicular to the dimer axis computed in the x-y plane for an incident plane wave at the excitation wavelength λ = 532 nm. c) Geometry of the TMD on a non-planarized dimer in the x-z plane, direction of Xo (red) and XD (dark blue) are indicated by arrows. d) Emission patterns of a dipole emitting at 750 nm and placed at different positions and heights (z) with respect to the non-planarized dimer shown in (b). The relative positions of the Si-NS and the dipoles are displayed in the x-y plane. The dotted circle corresponds to the 0.9 NA of the microscope objective, we hence detect emission from inside this circle. The dipole can be assimilated to an in-plane exciton (red arrow) or out-of-plane exciton (dark blue arrow) depending on its orientation and position. The surrounding refractive index is n= 1. e) Geometry of the TMD on a planarized dimer in the

FIG. 1 .

 1 FIG. 1. Nanoresonator arrays for non-planarized and planarized samples.a) Schematic of a bare monolayer on a non-planarized Si nanoresonator, which is referred to as a dimer and consists of 2 rectangles parallel to each other, separated by a gap G. b) Bright field optical image of the sample. Array of rectangular Si dimers with constant width W = 150 nm but varying length (L) and gap (G) covered by a single exfoliated WSe2 ML. c) Illustration of the TMD layer configuration when deposited on non-flattened Si antennas. d) Measured AFM profiles for configuration shown in c). e) Illustration of the TMD layer configuration when deposited on planarized Si nanoantennas for the same dimers shown in b). f ) Measured AFM profiles for configuration shown in e). g) Dark-field optical images showing the light scattered by the antennas. h) Normalized dark field intensity, measured for G = 300 nm, L = 300 nm dimer without the TMD MLs (bottom left on figure g -white circle). The two first modes (below 2 eV) are the magnetic (MD) and electric dipole (ED) modes. They cover the range of the exciton emission. The resonances at higher energy correspond to higher order modes (quadrupoles) and cover the excitation laser energy. Normalized PL spectra measured on ML WSe2 (orange) and MoSe2 (red) (far from antennas), laser excitation

3 FIG. 2 .

 32 FIG. 2. WSe2 monolayer on non-planarized nanoresonators.a) Normalized WSe2 PL color maps of 12 dimer resonators (from G = 300 nm, L = 300 nm to G = 50 nm, L = 200 nm), for x-polarized incident light (along the antennas). b) Normalized PL color map at the vicinity of resonator G = 300 nm, L = 300 nm for x-polarized incident light (corresponding to white circle in a)). Dimer position and shape are marked by the black rectangles on the map. c) d) e) show typical PL spectra taken along the white line on figure b) Here we discern three contributions namely the bright exciton X0 (blue curves and data points), the trion (green) and the dark exciton XD (red). f ) top: Normalized total PL amplitude (black line) and individual Lorentzian fitting component amplitude (red, green, blue) measured across the white line on (b) ; bottom: energy of the individual Lorentzian peaks versus position. The dimer position is shown by the blue square, dashed grey lines mark the positions where the PL spectra shown in panel c) d) e) have been taken. g) as a) but for y-polarized incident light (perpendicular to the antennas). h) Normalized PL color map at the vicinity of resonator G = 300 nm, L = 300 nm for y-polarized incident light. Dimer position marked by white rectangles. i) j) k) show typical PL spectra taken along the white line on figure h). Experimental data appear in black continuous line. Blue, green and red dashed lines are the constituting Lorentzian peak used to fit the PL, the corresponding global fit is shown as black dashed lines. l) same as f) but for y-polarized incident light.

X 7. 5 X 2 FIG. 3 .

 523 FIG. 3. MoSe2 monolayer on non-planarized nanoresonators. a) Normalized MoSe2 PL color maps of 12 dimer resonators (same nanoresonators as in Fig. 2a), for x-polarized on non-planarized nanoresonators. b) Normalized PL color map at the vicinity of nanoresonator G = 300 nm, L = 300 nm for x-polarized (corresponding to white circle in a)). Dimer positions and shape is marked by the back rectangles on the map. c) Normalized PL integrated amplitude (black line) and energy of the PL maximum is (blue line) versus position, cross section taken along the white dashed line on (b). The dimer position is marked by the blue squares. d) Typical PL spectra taken at 3 different positions marked by matching color dashed vertical lines in (b). The spectra are normalized, the normalization factor is indicated.

3 1 2 XFIG. 4 .

 24 FIG. 4. Strain free WSe2 monolayer on planarized nanoresonators. a) Normalized PL images of WSe2 monolayer deposited on G = 300 nm, L = 300 nm planarized nanoresonator for x-polarized laser excitation. Dimer position shown as black rectangles. b) Typical PL spectra taken along the white line on panel (a). The PL spectra are fitted using the same three Lorentzian peaks as detailed in Fig. ?? c) Top: Normalized total PL amplitude (black line) and individual Lorentzian fitting component amplitude (red, green, blue) measured across the white line on (a). The dimer position is marked by the blue squares. Bottom: Energy of the individual Lorentzian peaks versus position along the white line on panel (a). The positions where the PL spectra have been taken appear as dashed line, with matching color. d), e), f ) same as a), b) and c) but for y-polarized incident light.

Acknowledgements.-LPCNO acknowledges funding from ANR 2D-vdW-Spin, ANR MagicValley, ITN 4PHOTON No. 721394 and the Institut Universitaire de France. CEMES, LAAS, LETI and LPCNO acknowledge funding from ANR HiLight (ANR-19-CE24-0020-01). This study has been partially supported through the EUR grant NanoX 2D Light (ANR-17-EURE-0009) in the framework of the Programme des Investissements d'Avenir and by the LAAS-CNRS micro and nanotechnologies platform, a member of the French RENATECH network. K.W. and T.T. acknowledge support from the Elemental Strategy Initiative conducted by the MEXT, Japan , Grant Number JPMXP0112101001, JSPS KAKENHI Grant Number JP20H00354, and the CREST, JST Grant Number JPMJCR15F3. CEMES team acknowledges the CALMIP computing facility (grant P12167).

Supplementary Information.-The following file TMD-SiNSv11-Supplement.pdf as Supporting Information is available free of charge via the internet at http://pubs.acs.org. The Supporting Information includes the Raman mapping of WSe 2 monolayers transferred on non planarized (under strain) and planarized (strain free) nanoantennas. They are shown in Fig. S1 and Fig. S3, respectively. Fig. S2 presents the photoluminescence mapping of a MoSe 2 monolayer on non planarizd nanoantenna as function of the excitation laser polarization.


Pacôme Prompsy
email: pacome.prompsy@curie.fr

Pia Kirchmeier

Justine Marsolier

Marc Deloger

Nicolas Servant

Céline Vallot
email: celine.vallot@curie.fr

Interactive analysis of single-cell epigenomic landscapes with ChromSCape

Chromatin modifications orchestrate the dynamic regulation of gene expression during development and in disease. Bulk approaches have characterized the wide repertoire of histone modifications across cell types, detailing their role in shaping cell identity. However, these population-based methods do not capture cell-to-cell heterogeneity of chromatin landscapes, limiting our appreciation of the role of chromatin in dynamic biological processes. Recent technological developments enable the mapping of histone marks at single-cell resolution, opening up perspectives to characterize the heterogeneity of chromatin marks in complex biological systems over time. Yet, existing tools used to analyze bulk histone modifications profiles are not fit for the low coverage and sparsity of single-cell epigenomic datasets. Here, we present ChromSCape, a user-friendly interactive Shiny/R application that processes singlecell epigenomic data to assist the biological interpretation of chromatin landscapes within cell populations. ChromSCape analyses the distribution of repressive and active histone modifications as well as chromatin accessibility landscapes from single-cell datasets. Using ChromSCape, we deconvolve chromatin landscapes within the tumor micro-environment, identifying distinct H3K27me3 landscapes associated to cell identity and breast tumor subtype.

Introduction

Histone modifications are key regulators of gene expression, driving chromatin folding and genes accessibility to transcription machineries. The recent development of single-cell methods to study epigenomes now enables the appreciation of the heterogeneity of chromatin modifications within a population. These experimental methods assess the distribution of histone marks at single-cell resolution by coupling next generation sequencing to highthroughput microfluidics DNA barcoding (scChIP-seq) [START_REF] Grosselin | High-throughput single-cell ChIP-seq identifies heterogeneity of chromatin states in breast cancer[END_REF][START_REF] Rotem | Single-cell ChIP-seq reveals cell subpopulations defined by chromatin state[END_REF] or in situ reactions (scChIL-seq [START_REF] Harada | A chromatin integration labelling method enables epigenomic profiling with lower input[END_REF] , scChIC-seq [START_REF] Ku | Single-cell chromatin immunocleavage sequencing (scChIC-seq) to profile histone modification[END_REF] , scCUT&Tag [START_REF] Kaya-Okur | CUT&Tag for efficient epigenomic profiling of small samples and single cells[END_REF] ). In contrast to scATAC-seq approaches which identify open regions of the chromatin [START_REF] Chen | A rapid and robust method for single cell chromatin accessibility profiling[END_REF][START_REF] Cusanovich | Multiplex single cell profiling of chromatin accessibility by combinatorial cellular indexing[END_REF][START_REF] Lareau | Droplet-based combinatorial indexing for massive-scale single-cell chromatin accessibility[END_REF] , these methods can capture various chromatin states, enriched in repressive or active histone marks (H3K27me3 or H3K4me3 for example). Using these approaches, we can study the heterogeneity of epigenomes within complex biological samples, such as tumors [START_REF] Grosselin | High-throughput single-cell ChIP-seq identifies heterogeneity of chromatin states in breast cancer[END_REF] , and start appreciating the role of epigenomic diversity and the dynamics of chromatin in disease and development.

Existing tools used to analyze bulk ChIP-seq experiments are not fit for the low coverage and sparsity of these single-cell histone modifications datasets, which is due to the inherent low number of copy of DNA molecules per cell -maximum two for a diploid genome. Several computational methods for the analysis of scATAC-seq have been developed to deal with the specificities of single-cell DNA-based datasets. They were recently benchmarked [START_REF] Chen | Assessment of computational methods for the analysis of single-cell ATACseq data[END_REF] , with SnapATAC 10 , CisTopic [START_REF] Bravo González-Blas | cisTopic: cis-regulatory topic modeling on single-cell ATAC-seq data[END_REF] and Cusanovich2018 [START_REF] Cusanovich | The cis-regulatory dynamics of embryonic development at singlecell resolution[END_REF] being the top-three performing methods.

These tools, initially dedicated to scATAC-seq and without graphic interface, require some scripting skills. Biologists with limited computational training can manipulate and analyze scRNA-seq and scATAC-seq datasets using applications such as 'scOrange' [START_REF] Stražar | scOrange-a tool for hands-on training of concepts from single-cell data analytics[END_REF] and 'SCRAT' [START_REF] Ji | Single-cell regulome data analysis by SCRAT[END_REF] .

With ChromSCape (Figure 1), we propose a user-friendly, step-by-step and customizable Shiny/R application to analyze all types of sparse single-cell epigenomic datasets. The user can interactively identify subpopulations with common epigenomes within heterogeneous samples, find differentially enriched regions between subpopulations and interpret epigenomes by linking regions to associated genes and pathways. The pipeline starts from aligned sequences or count tables, and is designed for high-throughput single-cell datasets with samples containing as low as 100 cells with a minimum of 1,000 reads per cell up to 25,000 cells on a standard laptop. ChromSCape accepts multiple samples to allow comparisons of cell populations between and within samples. ChromSCape can determine cell identities from single-cell histone modification profiles, whatever the technology, as well as scATAC-seq datasets. We showcase the use of ChromSCape by deconvolving chromatin landscapes within the tumor microenvironment; we identify distinct H3K27me3 landscapes associated to cell identity and breast tumor subtype.

Results

ChromSCape identifies cell identities from scChIP-seq data

To test the efficiency of ChromSCape (Figure 1) in identifying cell sub-populations based on their epigenome (H3K27me3), we generated an in-silico dataset with known ground truth, mixing 4 different human cell types: Jurkat B cells, Ramos T cells, MDA-MB-468 breast cancer cells and HBCx-22 tumor cells derived from a luminal breast tumor PDX model [START_REF] Grosselin | High-throughput single-cell ChIP-seq identifies heterogeneity of chromatin states in breast cancer[END_REF] . Interestingly, Jurkat and Ramos cells were processed within the same microfluidics experiment, preventing the existence of any batch effect between them (see Grosselin et al., [START_REF] Grosselin | High-throughput single-cell ChIP-seq identifies heterogeneity of chromatin states in breast cancer[END_REF] ). We compared ChromSCape to methods specifically designed for single-cell epigenomic datasets (scATACseq) for their ability to identify cell identities. Based on a recent scATAC-seq benchmark [START_REF] Chen | Assessment of computational methods for the analysis of single-cell ATACseq data[END_REF] , we selected the top-performing methods, namely Cusanovich2018 12 , SnapATAC 10 and CisTopic 11 .

We also benchmarked EpiScanpy [START_REF] Danese | EpiScanpy: integrated single-cell epigenomic analysis[END_REF] , a recent analysis pipeline for various single-cell epigenomic data (scATAC-seq, scDNA methylation, …) developed in Python. We applied hierarchical clustering on the reduced feature space obtained by each method and used an ARI metric to evaluate their ability to identify cell phenotypes. ChromSCape with default parameters manages to separate almost perfectly the 4 cell types, with an ARI of 0.998 (Figure 2a), as Cusanovich2018 and CisTopic (both an ARI of 0.996, Figure 2b), followed closely by EpiScanpy (ARI of 0.940, Figure 2b). ChromSCape, EpiScanpy and SnapATAC were all run on 50kbp bins, but SnapATAC had noisier clusters and a slightly poorer ARI (0.822).

We also compared the agility of ChromSCape to manipulate and interpret scChIP-seq datasets to two applications with graphic interface, using the same reference dataset (Figure 2c), with either default settings or optimizing input and settings. scOrange is a stand-alone platform allowing researchers to create workflows to analyze single-cell datasets, offering a wide variety of analytical modules. While for scRNA-seq many workflows have been developed and are ready-to-use, in the case of epigenomic datasets, users need to have prior computational knowledge to organize a proper workflow. We managed to group cells according to sample of origin only with an optimized workflow (Figure 2c, default vs optimized). SCRAT is a Shiny/R package presented as a user interface to analyze single-cell epigenomic data. The default option for SCRAT is to count reads within 'ENCODE Clusters' corresponding to co-regulatory open chromatin regions obtained from DNAse-seq datasets, not adapted for the analysis of repressive histone marks like H3K27me3. In order to use SCRAT for our reference scChIP-seq datasets, we had to pre-compute counts on pre-defined peaks called on the 'pseudo-bulk' (see Methods and 'Optimized' panel), limiting the usability of SCRAT. In addition, in contrast to ChromSCape, both applications do not propose functionalities to associate genomic regions to gene annotation, limiting the biological interpretation of the results obtained with differential analysis.

Like single-cell transcriptomics approaches, single-cell epigenomic technologies can be influenced by various batch effects, e.g. library preparation, batch of hydrogel beads or efficacy of immuno-precipitation or cleavage. To overcome this, we have implemented in ChromSCape a module for batch correction based on the fastMNN method [START_REF] Haghverdi | Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors[END_REF] . To test this functionality, we used three datasets: two from our previous study (HBCx-95 and HBCx95-CapaR, collected with batch1 beads) and a new dataset, HBCx95 -batch 2, which is a biological replicate of HBCx-95, i.e a PDX tumor from the same PDX model but a different mouse, processed with a new batch of beads. Using ChromSCape, we analyzed together the three H3K27me3 tumor samples. As shown in Supplementary Figure 2, a strong batch effect separates the two biological replicates before batch correction (left panel). After applying batch correction (right panel), cells from the two biological replicates of untreated HBCx-95 tumors successfully mix, but not with cells from the resistant tumor, suggesting that the module corrects for batch effect without overcorrecting biological differences.

We also evaluated the usability of ChromSCape for other types of single-cell histone modification data obtained by other technologies than scChIP-seq. We analyzed two public datasets of scCUT&Tag and scChIC-seq targeting H3K27me3 and H3K4me3 marks respectively. ChromSCape facilitates the analysis of such public dataset as the user can directly upload the GEO single-cell BED files into the application. We recommend here for H3K27me3 mark -accumulating in broad peaks -to aggregate the signal into 50kbp bins, and for H3K4me3 mark -accumulating in sharp peaks -to count within 5kbp bins or around gene TSS (+/-2500bp). As shown in Supplementary Figure 1a, for the scCUT&Tag dataset the two K562 replicates showed no batch effect and were clustered together separately from H1 cells (ARI = 0.976). For the scChIC-seq dataset (Supplementary Figure 1b), 7 clusters are clearly observable on the UMAP representation, as was found by the authors in their study [START_REF] Ku | Single-cell chromatin immunocleavage sequencing (scChIC-seq) to profile histone modification[END_REF] .

ChromSCape classifies cells from scATAC-seq data

In order to assess the capacity of ChromSCape to analyze all types of single-cell epigenomic data, we re-analysed a scATAC-seq dataset (GSE99172) containing 8 cell lines and 4 patientderived cells. This dataset was partly produced and analyzed in a study by chromVar, a dedicated scATAC-seq analytical tool [START_REF] Schep | chromVAR: inferring transcription-factor-associated accessibility from single-cell epigenomic data[END_REF] ; we used the same color code as in this study. This dataset contains various biological samples as well as technical replicates for two cell lines, K562 & GM12878, for which there are 6 and 4 technical replicates respectively. Due to a relatively low number of cells per sample (n=96 per sample), we set the read count threshold to 1,000 for cells to be included in the analysis. We measured the ability of ChromSCape to classify cells according to cell type of origin using assignment scores for each sample X and each cluster Y (number of cells from sample X assigned to cluster Y / total number of cells of sample X).

In the unsupervised analysis, we identified the optimal number of clusters to be k = 5 according to the relative change in area under the CDF curve (Figure 3a and data not shown Supplementary Figure 3a). Cells from different technical replicates of K562 and GM12878 all grouped together in clusters 5 and 2 with assignment scores of 99.7% and 98.3% respectively (Figure 3b). Analyzing all samples together, ChromSCape could robustly identify -i.e as stable separated clusters by consensus clustering -TF1, K562 and GM12878 cells, affecting in average 99.4% of cells to correct cluster. Cluster C1 grouped together AML, Mono and LMPP samples with HL60 cell line, which is also originally derived from leukocytes of a patient with an AML cancer.

In order to get more insight into cell identities within cluster C1, we re-analyzed cells from C1 with ChromSCape. In contrast to the first round of analysis, ChromSCape was able to distinguish cell identities within samples, and detect individual clusters of cells, with high assignment scores for the two normal samples (LMPP & Monocytes) and HL-60 (average assignment score 98.3%, Figure 3c &d). Additionally, AML blasts from patient SU070 show a greater monocyte signature than patient SU353 (Figure 3d, p-value = 0.0025, Fisher's exact test, respectively 85.0 % and 25.7 % of SU070 and SU353 blasts cells cluster with monocytes), as previously described for these cells in [START_REF] Corces | Lineage-specific and single-cell chromatin accessibility charts human hematopoiesis and leukemia evolution[END_REF] . ChromSCape identifies distinct populations within the normal immune cell environment based on their chromatin accessibility. Within AML patient samples, ChromSCape matches each cancer cell to the closest resembling healthy population.

ChromSCape deconvolves epigenomes of the tumor micro-environment

To further showcase the use of ChromSCape, we interrogated the heterogeneity of chromatin states within the tumor micro-environment of two breast tumor subtypes: luminal and triplenegative (TNBC) breast tumors. Tumor micro-environment is a key player in tumor evolution processes, and can vary between tumor types and with response to cancer therapy. Here our goal was to compare H3K27me3 landscapes of cells from the tumor micro-environment of luminal and TNBC subtypes, resistant or not to cancer treatment. The HBCx-22 and HBCx-22-TamR datasets correspond to mouse cells from a pair of luminal ER + breast PDXs 1 : HBCx-22, responsive to Tamoxifen and HBCx-22-TamR, resistant to Tamoxifen. The HBCx-95 and HBCx-95-CapaR correspond to a triple-negative breast cancer (TNBC) tumor model of acquired resistance to chemotherapy [START_REF] Grosselin | High-throughput single-cell ChIP-seq identifies heterogeneity of chromatin states in breast cancer[END_REF] . We analyzed together these four H3K27me3 mouse scChIP-seq datasets, two of which had not been analyzed in our previous study [START_REF] Grosselin | High-throughput single-cell ChIP-seq identifies heterogeneity of chromatin states in breast cancer[END_REF] . Using ChromSCape, we propose a comprehensive view of cell populations based on their chromatin profiles, and show the identification of tumor-type and treatment-specific cell populations and respective chromatin features. All plots in Figure 4 were automatically generated by the application and are downloadable from the interface. In the quality filtering step, a threshold of 2,000 reads per cell was set due to a relatively high initial number of cells (n = 5,516 cells).

After the dimensionality reduction step (Figure 4a & b), we applied our consensus clustering approach on the filtered dataset with k = 2 to k = 10 clusters. We chose to partition the data into k = 4 clusters based on the knee method, as a plateau in the relative change in area under the CDF curve was observed between k = 4 and k = 5 clusters (Figure 4c-d & Supplementary Figure 34a-b). Consensus score matrix in Figure 4d shows that most of the cells were stably assigned to four chromatin-based populations (mean consensus score for elected clusters of 0.91 is significantly higher than mean consensus score for other clusters, 0.17, p-value = 2.2e-16, twosided Student's t-test). Assignment of cells to cluster C2 and C4 is significantly less stable than C1 and C3 (p-value < 2.2e-16, Student's two-sided t-test, mean consensus scores are respectively 0.84 and 0.90 for C1-C3 and 0.70 and 0.71 for C2-C4, see Supplementary Figure 3a), suggesting that cells from C2 and C4 might share H3K27me3 features, whereas cells from C1 and C3 have distinct H3K27me3 landscapes. Clusters C1, C2 and C4 contain cells from all four samples, with a significantly higher proportion of HBCx-22-TamR for C1 (p-value = 3e-05, Pearson's Chi-squared test) (Figure 4e). On the other hand, cluster C3 is almost exclusively composed of cells from model HBCx-95 (Figure 4c&e), revealing a stromal cell population specific to the triple negative breast cancer model (HBCx-95).

To further identify the specific features of each chromatin-based population, we proceeded to peak calling, differential analysis and gene set enrichment analysis using default parameters (see Methods). As H3K27me3 is a repressive histone mark, we focused our analysis on loci depleted in H3K27me3, where transcription of genes can occur, in cells from each cluster versus all other cells. The differential analysis identified respectively 189, 210, 83 and 9 significantly depleted regions for clusters C1 to C4 (Figure 4g, logFC < 1, adjusted p-value < 0.01). We found loci devoid of H3K27me3 specific to cluster C2, enriched for genes involved in apical junction such as Bcar1 (Figure 4f) and Ptk2, which are characteristic of genes expressed in fibroblasts. We found a depletion of H3K27me3 specific to cluster C3 over the genes Nrros (Figure 4f) and Il10ra, two genes characteristic of immune expression programs. Depletion of H3K27me3 over the transcription start site of Rap1gap2, a gene expressed in endothelial cells, was a key feature of cluster C4 (Figure 4f). For cluster C1 and C2, we found a depletion of H3K27me3 over Eln, a gene expressed in fibroblasts.

Gene set enrichment analysis for genes located in regions depleted of H3K27me3 enrichment only revealed very few enriched gene lists, mostly for cluster C2 (q-value < 0.1, Figure 4h, multiple gene sets related to stem and cancer cells) and one list for C1 ("LPS_VS_CONTROL_MONOCYTE_UP"). Linking H3K27me3 enrichment to transcription is indeed indirect, we envisage such enrichment analysis more appropriate for H3K4me3 scChIP-seq in which enriched regions are directly associated to gene transcription.

Overall, these results are consistent with our previous analysis of HBCx-95 scRNA-seq datasets where subpopulations were differentially expressing markers of fibroblasts, endothelial and macrophage cells [START_REF] Grosselin | High-throughput single-cell ChIP-seq identifies heterogeneity of chromatin states in breast cancer[END_REF] . This new analysis comprising the HBCx-22 dataset allowed us to identify the H3K27me3 signature of potential endothelial cells (cluster C4). These cells are present in each model, but might not have been previously detected in the previous scChIP-seq analysis due to low cell representation. In addition, the H3K27me3 signature of potential immune cells is restricted to cells from the TNBC model (cluster C3), suggesting that these immune cells are absent from the luminal tumor.

Discussion

ChromSCape is a Shiny/R application designed for both biologists and bioinformaticians to analyze complex chromatin profiling datasets such as scChIP-seq datasets. The comprehensive application is quick to take over plus the direct visualization of cells clusters combined to configurable parameters and incremental saving of intermediary R objects eases bench-marking of parameters. We show that ChromSCape performs as well or better than state of the art singlecell epigenomic analytic tools to identify cell identities from an in-silico mix of H3K27me3 scChIP-seq datasets. It also manages to identify sub-populations within a complex scATACseq benchmarking dataset, showing its wide range of application for epigenomic analysis. In addition, using ChromSCape to study the epigenome of mouse stromal cells in breast tumors, we can identify the various epigenomes within the tumor micro-environment. Overall, we see ChromSCape as a useful tool to probe heterogeneity and dynamics of chromatin profiles in various biological settings, not only in cancer development but also in cell development and cellular differentiation.

Methods

Implementation

ChromSCape is an R package developed in Shiny/R. It uses various Shiny related packages (shinyjs, shinydashboard, shinyDirectoryInput) for the user interface. The application takes advantage of public R libraries for data vizualisation (RcolorBrewer, colorRamps, Rtsne, umap, colourpicker, kableExtra, knitr, viridis, ggplot2, gplots, png, grid, gridExtra, DT) as well as for data manipulation (Matrix, dplyr, tidyr, stringr, irlba, rlist, qualV, stringdistr). ChromSCape uses Bioconductor packages (i) for the manipulation of single-cell data with SingleCellExperiment, scater [START_REF] Mccarthy | Scater: pre-processing, quality control, normalization and visualization of single-cell RNA-seq data in R[END_REF] , scran [START_REF] Lun | A step-by-step workflow for low-level analysis of single-cell RNA-seq data with Bioconductor[END_REF] , (ii) for the manipulation of genomic regions with IRanges and GenomicRanges 17 , (iii) for the manipulation of genomic files with Rsamtools and BiocParallel, (iv) for the correction of batch effects with batchelor 18 and (v) to determine the optimal number of clusters with ConsensusClusterPlus [START_REF] Wilkerson | ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking[END_REF] . In addition, ChromSCape makes use of custom R functions which serve for both manipulation and visualization of datasets. Brief command lines enable users without any bioinformatics skills to install all R dependencies and run the application in a web browser.

Demonstration application

A demonstration of ChromSCape is freely available at https://vallotlab.shinyapps.io/ChromSCape/.

Input datasets, quality control and pre-processing

Input files for ChromSCape are either one or multiple count matrices with genomic regions in rows and cells in columns or single-cell BAM or BED files. In this case, a directory containing single-cell BAM or BED files must be specified and the count matrix created by aggregating the signal into successive genomic bins, peaks (BED file must be provided by user) or into regions around genes Transcription Start Sites (TSS). For H3K27me3 scChIP-seq datasets, with a distribution in broad peaks, we recommend using bins of 50kbp, while for H3K4me3 scChIPseq or scATAC-seq datasets we recommend using smaller bins (e.g 5kb), knowns peaks or regions around TSS. The 'condition' or 'label' of each cell is then heuristically determined using file names and the number of conditions specified by the user. Guidelines and link towards datasets are given in the user guide (https://vallotlab.github.io/ChromSCape/ChromSCape_guide.html).

In order to efficiently remove outlier cells from the analysis, e.g. cells with excessively high or low coverage, the user sets a threshold on a minimum read count per cell and the upper percentile of cells to remove. The latter could correspond to doublets, e.g. two cells in one droplet, while lowly covered cells are not informative enough or may correspond to barcodes ligated to contaminant DNA or library artifacts. Regions not supported by a minimum userdefined percentage of cells that have a coverage greater than 1,000 reads are filtered out.

Defaults parameters were chosen based on the analysis of multiple scChIP-seq datasets from our previous study 1 : a minimum coverage of 1,600 unique reads per cell, filtering out the cells with the top 5% coverage and keeping regions detected in at least 1 % of cells. Post quality control filtering, the matrices are normalized by total read count and region size. At this step, the user can provide a list of genomic regions, in BED format, to exclude from the subsequent analysis, in cases of known copy number variation regions between cells for example.

To reduce the dimensions of the normalized matrix for further analysis, principle component analysis (PCA) is applied to the matrix, with centering, and the 50 first PCs are kept for further analysis. The user can visualize scChIP-seq data after quality control in the PCs dimensional space. The t-distributed stochastic neighbor embedding (t-SNE) algorithm [START_REF] Maaten | Visualizing Data using t-SNE[END_REF] and UMAP [START_REF] Mcinnes | UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction[END_REF] is applied on the PCA to visualize the data in two dimensions. The PCA and t-SNE plots are a convenient way to check if cells form clusters in a way that was expected before any clustering method is applied. For instance, the user should verify whether the QC filtering steps and normalization procedure were efficient by checking the distribution of cells in PC1 and PC2

space. Cells should group independently of normalized coverage. In our hands, for our scChIPseq H3K27me3 datasets, minimum coverage of 1,600 unique reads per cell was required to separate cells independently of coverage post normalization [START_REF] Grosselin | High-throughput single-cell ChIP-seq identifies heterogeneity of chromatin states in breast cancer[END_REF] . A batch correction option using mutual nearest neighbors 'FastMNN' function from 'batchelor' package [START_REF] Haghverdi | Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors[END_REF] is implemented to remove any known batch effect in the reduced feature space.

Hierarchical clustering, filtering and consensus clustering

Using the first 50 first PCs of computed PCA as input, hierarchical clustering is performed, taking 1-Pearson's correlation score as distance metric. To improve the stability of our clustering approaches and to remove from the analysis isolated cells that do not belong to any subgroup, cells displaying a Pearson's pairwise correlation score below a threshold t with at least p % of cells are filtered out (p is set at 1 % by default). The correlation threshold t is calculated as a user-defined percentile of Pearson's pairwise correlation scores for a randomized dataset (percentile is recommended to be set as the 99th percentile). Correlation heatmaps before and after correlation filtering and the number of remaining cells are displayed to inform users on the filtering process.

ChromSCape uses Bioconductor ConsensusClusterPlus package [START_REF] Wilkerson | ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking[END_REF] to determine what is the appropriate k-partition of the filtered dataset into k clusters. To do so, it evaluates the stability of the clusters and computes item consensus score for each cell for each possible partition from k=2 to 10. For each k, consensus partitions of the dataset are done on the basis of 1,000 resampling iterations (80% of cells sampled at each iteration) of hierarchical clustering, with Pearson's dissimilarity as the distance metric and Ward's method for linkage analysis. The optimal number of clusters is then chosen by the user; one option is to maximize intra-cluster correlation scores based on the graphics displayed on the 'Consensus Clustering' tab after processing. Clustering memberships can be visualized in two dimensions with the t-SNE or UMAP plot.

Peak calling for genomic region annotation

This step of the analysis is optional, but recommended in order to refine the peak annotation prior to enrichment analysis. To be able to run this module, some additional command line tools are required such as SAMtools [START_REF] Li | The Sequence Alignment/Map format and SAMtools[END_REF] , BEDTools 23 and MACS2 [START_REF] Zhang | Model-based analysis of ChIP-Seq (MACS)[END_REF] . The user needs to input BAM files for the samples (one separate BAM file per sample), with each read being labeled with the barcode ID. ChromSCape merges all files and splits them again according to the previously determined clusters of cells (one separate BAM file per cluster). Customizable significance threshold for peak detection and merging distance for peaks (defaults to p-value=0.05 and peak merge distance to 5,000) allows to identify peaks in close proximity (<1,000bp) to a gene transcription start site (TSS); these genes will be later used as input for the enrichment analysis.

For the annotation, ChromSCape uses the reference human transcriptome Gencode_hg38_v26, limited to protein coding, antisense and lncRNA genes.

Differential Analysis and pathway enrichment analysis

To identify differentially enriched regions across single-cells for a given cluster, ChromSCape can perform (i) a non-parametric two-sided Wilcoxon rank sum test comparing normalized counts from individual cells from one cluster versus all other cells, or cluster of choice, or (ii) a parametric test comparing raw counts from individual cells, using edgeR [START_REF] Robinson | a Bioconductor package for differential expression analysis of digital gene expression data[END_REF] , based on the assumption that the data follows a negative-binomial distribution. We test for the null hypothesis that the distribution of normalized counts from the two compared groups have the same median, with a confidence interval 0.95. The calculated p-values are then corrected by the Benjamini-Hocheberg procedure [START_REF] Benjamini | Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing[END_REF] . The user can set a log2 fold-change threshold and corrected p-value threshold for regions to be considered as significantly differentially enriched (default settings are a p-value and log2 fold-change thresholds respectively of 0.01 and 1). If users have specified batches, the differential analysis is done using the 'pairwiseWilcox' function from the scran package [START_REF] Lun | A step-by-step workflow for low-level analysis of single-cell RNA-seq data with Bioconductor[END_REF] , setting the batch of origin as a 'blocking level' for each cell.

For the top 100 most significant differential regions, single-cell H3K27me3 enrichment levels can be visualized overlaying H3K27me3 counts for each cell at selected genes onto a t-SNE plot. Using the refined annotation of peaks done in previous step, the final step is to look for enriched gene sets of the MSigDB v5 database [START_REF] Subramanian | Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles[END_REF] scATAC-seq datasets: The scATAC-seq dataset is composed of two cell types derived from two acute myeloid leukaemia (AML) (patients SU070 and SU353 blastocytes (blast) and leukemic stem cells (LSC) from [START_REF] Corces | Lineage-specific and single-cell chromatin accessibility charts human hematopoiesis and leukemia evolution[END_REF] ) as well as multiple cell lines : GM12878 (4 replicates), TF1, BJ, H1, HL60, K562 (3 replicates) from [START_REF] Buenrostro | Single-cell chromatin accessibility reveals principles of regulatory variation[END_REF] , K562 (3 replicates) from [START_REF] Schep | chromVAR: inferring transcription-factor-associated accessibility from single-cell epigenomic data[END_REF] ; monocytes (Mono) and lymphoid primed multipotent progenitor (LMPP) from [START_REF] Corces | Lineage-specific and single-cell chromatin accessibility charts human hematopoiesis and leukemia evolution[END_REF] . The count matrix of reads in peaks was downloaded from GEO accession number GSE99172

[https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE99172], split into distinct matrices for each sample and formatted to be accepted as input by ChromSCape.

H3K27me3 scChIP-seq mouse datasets: The samples correspond to mouse cells from patientderived xenograft (PDX) originating from two different human donors [START_REF] Grosselin | High-throughput single-cell ChIP-seq identifies heterogeneity of chromatin states in breast cancer[END_REF] . Raw FASTQ reads were processed using the latest version of our scChIP-seq data engineering pipeline (see above)

to produce 50kbp binned count matrices given as input to ChromSCape (matrices available at https://figshare.com/projects/Single-Cell_ChIP-seq_of_Mouse_Stromal_Cells_in_PDX_tumour_models_of_resistance/66419).

Cell line

MDA-MB-468 cells, bought at ATCC (HTB-132™), were cultured in DMEM 1640 (Gibco-BRL) and supplemented with 10% heat-inactivated fetal calf serum. Cell numbers, as judged by Trypan Blue exclusion test, were determined by counting cells using a Countess automated cell counter (Invitrogen). Cells were cultured at 37 °C in a humidified 5% CO2 atmosphere.

The cell line was mycoplasma negative. The MDA-MB-468 cells were trypsinized (Trypsin, Gibco-BRL). Prior to single-cell ChIP-seq, cells were then re-suspended in PBS/0.04% BSA (ThermoFisher Scientific, # AM2616).

Patient-derived xenograft (PDX)

Female Swiss nude mice were purchased from Charles River Laboratories and were maintained under specific pathogen-free conditions. Their care and housing were in accordance with institutional guidelines and the rules of the French Ethics Committee (project authorization no.

02163.02). A PDX from a residual triple-negative breast cancer post neo-adjuvant chemotherapy (HBCx-95) was previously established at Institut Curie with informed consent from the patient [START_REF] Grosselin | High-throughput single-cell ChIP-seq identifies heterogeneity of chromatin states in breast cancer[END_REF] . Prior to single-cell ChIP-seq, PDX were digested at 37°C for 2h with a cocktail of Collagenase I (Roche, # 11088793001) and Hyaluronidase (Sigma, # H3506). Cells were then individualized at 37°C using a cocktail of 0.25% trypsin/Versen (ThermoFisher

  H3K4me3 scChIC-seq human white blood cells dataset 4 : n=285 white blood cells from a human male donor were downloaded as gzipped single-cell BED files from GSE105012[https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE105012], inputted directly into

	ChromSCape and aggregated into 50kbp bins (default).
	H3K27me3 scCUT&Tag human H1 and K562 cells 5 : A replicate of K562 cell line comprising
	of	n=908		cells		from		GSE124680
	[https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE124680], another replicate of
	n=479	K562	cells	and	n=486	H1	cells	from	GSE124690
	[https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE124690] were downloaded as
	gzipped single-cell BED files, inputted directly into ChromSCape and aggregated around gene
	TSS (+/-2500bp).
	Datasets
	H3K27me3 scChIP-seq human in-silico mix of 4 cell types: The samples correspond to n=326
	human tumor cells from untreated PDX (HBCx-22), n=201 human T cells (Jurkat) and n=306
	B cells (Ramos) taken from 1 and n=454 cells from the MDA-MB-468 triple-negative breast
	cancer cell line (HBCx-22, Jurkat and Ramos data are from GSE117309

within differentially enriched regions (either enriched or depleted regions in the studied histone mark). We apply hypergeometric tests to identify gene sets from the MSigDB v5 database over-represented within differentially enriched regions, correcting for multiple testing with the Benjamini-Hochberg procedure. Users can then visualize most significantly enriched or depleted gene sets corresponding to the epigenetic signatures of each cluster and download gene sets enrichment tables. [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE117309]), MDA-MB-468 is available at GSE152502 [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE152502]). H3K27me3 scChIP-seq human datasets: The samples correspond to human cells from patientderived xenograft (PDX) originating from two different human donors 1 . For this study, we added a new scChIP-seq dataset, corresponding to a biological replicate of HBCx-95 (GSE152502 [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE152502]), processed with a novel batch of hydrogel beads.

Funding

This research project was supported by the ATIP Avenir program, by Plan Cancer and by the SiRIC-Curie program SiRIC Grants #INCa-DGOS-4654 and #INCa-DGOS-Inserm_12554.

Code availability

Source code, guidelines for installation and use of the application are provided at https://github.com/vallotlab/ChromSCape. A docker container containing the application and it's dependecies is available on DockerHub (pacomito/chromscape:v0.0.9001), instructions on how to launch it are available on the github page. Codes for the benchmark of 'SnapATAC', 'CisTopic', 'Cusanovich2018' and 'EpiScanpy' are available at https://github.com/vallotlab/ChromSCape_benchmarking. Data availability In this study, in addition of using publicly available single-cell datasets, we produced H3K27me3 scChIP-seq data for MDA-MB-468 sample and a new replicate of untreated HBCx-95. The sequencing data that support the findings of this study have been deposited in the National Center for Biotechnology Information Gene Expression Omnibus (GEO) and are accessible through the GEO Series accession number GSE152502 [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE152502]. All other relevant data are available from the corresponding author on request.

Scientific, #15040-033), Dispase II (Sigma, #D4693) and Dnase I (Roche, # 11284932001).

Red Blood Cell lysis buffer (ThermoFisher Scientific, # 00-4333-57) was then added to degrade red blood cells. In order to increase the viability of the cell suspension, dead cells were removed using the Dead Cell Removal kit (Miltenyi Biotec). Cells were re-suspended in PBS/0.04% BSA (ThermoFisher Scientific, # AM2616).

Single-cell ChIP-seq

The protocol for scChIP-seq was rigorously the same as in Grosselin et al., [START_REF] Grosselin | High-throughput single-cell ChIP-seq identifies heterogeneity of chromatin states in breast cancer[END_REF] and can be resumed by the main following steps. Cells were first compartmentalized into droplets containing Mnase in a microfluidics chip, then fused with barcoded hydrogel beads. After fusion of cellcontaining droplets and bead-containing droplets, Fast-link DNA ligase [Lucigen, # LK0750H] was used to ligate segmented DNA to barcodes. Droplets were pooled and used for chromatin immuno-precipitation against anti-H3K27me3 antibody ([Cell Signaling Technology, # 9733]).

After treatment with RNAse A (ThermoFisher Scientific, #EN0531) and Proteinase K (ThermoFisher Scientific, # EO0491), barcoded-nucleosomes were then amplified by in-vitro transcription using the T7 MegaScript kit (ThermoFisher Scientific, # AM1334) and reversetranscribed. After RNA digestion, DNA was amplified by PCR. The final product was sizeselected by gel electrophoresis. Single-cell ChIP-seq libraries were finally sequenced on an Illumina NextSeq 500 MidOutput 150 cycles.

Demultiplexing and alignment of H3K27me3 scChIP-seq datasets

Raw FASTQ reads were processed using the latest version of our scChIP-seq data engineering pipeline that allowed a more precise removal of PCR and RT duplicates ( barcodes and m genomic bins used in downstream analysis.

Benchmark of tools for scChIP-seq data analysis

Three methods dedicated to the analysis of scATAC-seq with the best performance according to Chen et al., 2019 [START_REF] Chen | Assessment of computational methods for the analysis of single-cell ATACseq data[END_REF] , were tested on a mixture of H3K27me3 scChIP-seq datasets (see Datasets below), namely 'SnapATAC', 'CisTopic' and 'Cusanovich2018'. The scripts were taken from the GitHub repository of the benchmark paper (https://github.com/pinellolab/scATACbenchmarking). For 'CisTopic' and 'Cusanovich2018', peaks were called using MACS2 with options '--nomodel --extsize 300 --keep-dup all --broad'. Peaks closer to 5000bp were merged together using BEDTools. For 'SnapATAC', 50kbp bins were counted from BAM files using 'SnapTools'. In addition, we also tested a recent method for single-cell epigenomic analysis, 'EpiScanpy', following the basic steps described in the tutorial for scATAC-seq (https://github.com/colomemaria/epiScanpy/blob/master/docs/tutorials/Tutorial_Hackathon_B uenrostro_2.html) with the same 50kbp matrices used for ChromSCape. We extracted from each method the matrix of reduced feature space, and used hierarchical clustering with Pearson's dissimilarity as the distance metric and Ward's method for linkage. The adjusted Rand's index (ARI), a widely-used measure to quantify clustering accuracy, was calculated for each method using R package 'mclust' [START_REF] Scrucca | mclust 5: Clustering, Classification and Density Estimation Using Gaussian Finite Mixture Models[END_REF] , taking samples of origin as 'true' clusters.

In addition, two softwares with graphic interface, dedicated to the analysis of single-cell data, 'scOrange' [START_REF] Stražar | scOrange-a tool for hands-on training of concepts from single-cell data analytics[END_REF] and 'SCRAT' [START_REF] Ji | Single-cell regulome data analysis by SCRAT[END_REF] , were also tested on the same set of cells both with 'default' parameters and manually 'optimized' parameters. 'default' for scOrange corresponds to using the template called 'Loading data from 10X protocols', a workflow meant for analyzing scRNA-seq of bone marrow cells, replacing the input by our matrices of selected cells in 50kbp bins. The 'optimized' workflow is available at www.github.com/vallotlab/ChromSCape_benchmarking and can be opened with the 'scOrange' software. For 'SCRAT', we found that the 'optimized' counting method corresponded to counting signal within peaks called on the 'pseudo-bulk' (see above).

In order to be able to compare the distinct methods, ChromSCape was first run on the raw count matrices and a set of 1287 cells passing the quality control thresholds were selected to be used as input for all methods. As the number of cells in each sample was unbalanced (e.g. the raw MDA-MB-468 containing n=3,382 cells while others have a maximum of n=456 cells), 500 cells from MDA-MB-468 were randomly sub-sampled using ChromSCape 'Perform Subsampling' option. We removed from the analysis the segments corresponding to known amplifications and homozygous loss of DNA of the Triple Negative Breast Cancer cell line MDA-MB-468, corresponding to a total of 77Mbp, previously found by analyzing the input of bulk ChIP-seq of the same cells (see Supplementary Note 2).

End Notes

Authors contribution

PP, PK and CV wrote code for ChromSCape and performed data analysis. JM performed scChIP-seq experiments on PDX and breast cancer cell line. PP, MD and NS processed raw scChIP-seq datasets into count matrices. PP and CV wrote the manuscript. CV conceptualized and supervised this work.

Competing interests

Licenses have been filed on some aspects of this work by Institut Curie and CNRS; contributors may receive payments related to exploitation of these licenses under their employer's rewards to inventor scheme.

Legends Figure 1. Representation of ChromSCape workflow

Users upload single-cell epigenomic data formatted as count matrices, single-cell BAM or single-cell BED files to start the analysis. The application includes Quality Control (QC),

Classification and Interpretation tools. The user can save plots and tables in png, pdf or csv formats, and R analysis objects in RData format.