C Villeneuve-Faure

L Boudou

K Makasheva

G Teyssedre

Atomic Force Microscopy developments for probing space charge at sub-micrometer scale in thin dielectric films

Keywords: Atomic Force Microscopy, Dielectric materials, Space charge, Kelvin Force Microscopy, Electrostatic force

Atomic force microscopy developments for probing space charge at sub-micrometer scale in thin dielectric films

Christina Villeneuve-Faure, Laurent Boudou, Kremena Makasheva,

INTRODUCTION

Physical dimension of solid dielectrics, involved in industrial applications, range from few centimeters (electrotechnics) to few tens of nanometers (electrets or microsystems). Whatever their main scale, most insulating materials accumulate electrical charges under applied electric field. This phenomenon is useful for some applications as electrets or nonvolatile memories [START_REF] Gerhard-Mulhaupt | Electrets 3rd edn[END_REF] but in many other instances it leads to a failure [START_REF] Witvrouw | Materials issues in the processing, the operation and the reliability of MEMS[END_REF]. Consequently, performing space charge measurements is needed to predict and to analyze these charging effects.

Conventional space charge measurements are based on acoustic or thermal perturbation whose shape influences strongly the ultimate resolution [START_REF] Fleming | Space Charge Profile Measurement Techniques: Recent Advances and Future Directions[END_REF][START_REF] Holé | Resolution of Direct Space Charge Distribution Measurement Methods[END_REF]. Acoustic perturbation as implemented in Pulsed ElectroAcoustic (PEA) technique [START_REF] Dagher | A Preliminary Study of Space Charge Distribution Measurements at Nanometer Spatial Resolution[END_REF] provides up to 3-4µm of in-depth resolution but without lateral resolution. Concerning thermal perturbation, thermal pulse (TP) and Laser Induced Modulation Method (LIMM) manage to reach 1µm of in-depth resolution [START_REF] Arnaout | Study of signal treatment for a pulsed electro-acoustic measurement cell: a way of improving the transfer matrix condition number[END_REF], but without lateral resolution. The development of Focused LIMM (FLIMM) has brought lateral resolution of the order of 5µm [START_REF] Petre | Threedimensional space charge cartographies by FLIMM in electron irradiated polymers[END_REF]. With recent developments, the actual resolution is improved and now limited mainly by instrumental point of view (perturbation shape, attenuation/dispersion for acoustic waves, deconvolution procedure [START_REF] Arnaout | Study of signal treatment for a pulsed electro-acoustic measurement cell: a way of improving the transfer matrix condition number[END_REF]…) Table 1 summarizes the resolution reached and the expected one for thermal and acoustic methods. Though some improvements have been achieved, up to now such techniques are unsuitable for characterizing the charge accumulation in thin films (a few 100 nm thickness) for microelectronics and microsystems applications or for studies of phenomena occurring at metal / dielectric or semiconductor / dielectric interfaces. Therefore, two ways appear to improve resolution: to push the instrumental limits of existing methods (according to Table 1, 100nm-resolution can be reached which would enable investigation of 1µm-thick dielectric films) or to develop new methods with resolution adapted to the expected scale considering both in plane and in-depth charge localization [START_REF] Teyssedre | Challenges in probing space charge at sub-micrometer scale[END_REF].

3-4µm 0.5µm 100nm

Techniques derived from atomic force microscopy (AFM) due to their sub-nanometer resolution and versatility to measure both the surface topography and electrical properties appear to be suited candidates for probing space charge.

Previous studies have shown the possibility to locally inject charges with a conductive AFM tip, and to study charge retention in dielectric thin films by AFM. In the literature, Electrostatic Force Microscopy (EFM) was extensively used to measure charge densities [START_REF] Schönenberger | Charge flow during metal-insulator contact[END_REF][START_REF] Schönenberger | Observation of single charge carriers by force microscopy[END_REF][START_REF] Melin | Electric force microscopy of individually charged nanoparticles on conductors: An analytical model for quantitative charge imaging[END_REF]. Some of studies highlight the limitation of the EFM and particularly its sensitivity to image charges which biases results to measurement conditions [START_REF] Lambert | Relationship between charge distribution and its image by electrostatic force microscopy[END_REF]. To avoid these effects, Kelvin Force Microscopy (KFM) is more and more used because it provides a direct measurement of surface potential modifications induced by the stored charges [START_REF] Rezende | Detection of charge distribution in insulator surfaces[END_REF]. However, charge density calculation remains difficult although some modeling attempts have been made [START_REF] Palleau | Numerical simulations for quantitative analysis of AFM electrostatic nanopatterning on PMMA by Kelvin Force Microscopy[END_REF], mainly due to the strong assumptions on the distribution of charges in the layer.

To overcome these issues and to obtain quantitative information on the charge density with a spatial resolution in the nanometer range, another technique is required. Our interest is focused on interaction forces between the AFM tip and the surface through force distance curves (FDC) [START_REF] Cappella | Force-distance curves by atomic force microscopy[END_REF][START_REF] Gomez-Monivas | Electrostatic forces between sharp tip and metallic and dielectric surface[END_REF]. FDC is sensitive to all interaction forces as molecular adhesion, Pauli interaction (~1nN), Van der Waals (~5-10nN), capillarity (~100nN) and electrostatic forces (from 1pN to fews tens of nN) [START_REF] Patil | Investigation in the electrostatic forces in scanning probe microscopy at low bias voltage[END_REF][START_REF] Kim | Dependence of approaching velocity on the Force-Distance curve in AFM[END_REF][START_REF] Jenke | Two-dimensional electrostatic force field measurement with simultaneous topography measurement on embedded interdigitated nanoelectrodes using FDC[END_REF]. Especially, the sensitivity to electrostatic forces can be exploited to measure space charge as electrostatic force between the tip and the surface and can be modeled [START_REF] Cappella | Force-distance curves by atomic force microscopy[END_REF][START_REF] Colchero | Resolution enhancement and improved data interpretation in electrostatic force microscopy[END_REF][START_REF] Gros | Two-dimensional, electrostatic finite element study of tip-substrate interactions in electric force microscopy of high density interconnect structure[END_REF], with a good correspondence with the experimental results except for short tip-sample distance [START_REF] Gros | Two-dimensional, electrostatic finite element study of tip-substrate interactions in electric force microscopy of high density interconnect structure[END_REF].

This paper highlights how techniques derived from AFM provide promising information about trapped charges localization in thin dielectric films. After a quick presentation of experimental conditions, charges injection using AFM tip is described. The following two parts address different techniques under investigation. The first part is dedicated to KFM technique. The short introduction of KFM method is followed by some results in order to address relevant parameters which can be extracted. The second part is dedicated to the new technique (named EFDC for Electrostatic Force Distance Curve) which is based on FDC measurement. After the method description, experimental results are presents to highlight information provided by EFDC. Finally, the relative merit of KFM and EFDC is investigated.

EXPERIMENTS

To compare EFDC and KFM methods, well-controlled thin dielectric layers are needed. Our interest is focused on silicon oxynitride layers (a-SiO x N y :H, x < 2, y < 1) whose thickness is of the order of 100 nanometers. The thin layers are elaborated by radiofrequency sustained plasma using the mixture of three gases (SiH 4 , N 2 O and He) [START_REF] Makasheva | Dielectric Layers for RF-MEMS Switches: Design and Study of Appropriate Structures Preventing Electrostatic Charging[END_REF]. The ratio of the precursors γ=[N 2 O]/[SiH 4 ] (partial pressures) allows adjusting silicon contents in the dielectric film and consequently modeling their electrical properties. In this case layers synthesized with γ=100 present properties close to thermal SiO 2 whereas in those with γ=5 the volume conductivity is increased by a order of magnitude due to the higher silicon contents [START_REF] Makasheva | Dielectric Layers for RF-MEMS Switches: Design and Study of Appropriate Structures Preventing Electrostatic Charging[END_REF][START_REF] Villeneuve-Faure | Kelvin force microscopy characterization of charging effect in thin a-SiOxNy:H layers deposited in pulsed plasma enhanced chemical vapor deposition process by tuning the Siliconenvironment[END_REF]. These well controlled materials are consequently ideal to identify the relative merit of KFM and EFDC methods.

All AFM measurements were done using a Bruker Multimode 8 equipment. As electrical measurements are very sensitive to humidity all KFM and EFDC results are acquired under N 2 atmosphere after thermal treatment of the samples at 100°C for 15min. Indeed, the thermal treatment avoids capillarity issues or charge injection in water layer adsorbed on the sample surface, which can screen phenomena occurring in the dielectric layer and be a source of non-reproducibility.

KFM measurements were performed in amplitude modulation mode. As sample surface is very flat (average roughness of 2.5nm), a 5nm-lift was used to keep high lateral resolution and potential sensitivity. A Pt-coated silicon tip provided by Bruker (SCM-PIT) was used for all measurements. Charging was achieved by bringing at contact the AFM tip and the dielectric surface and applying a potential in the range 0 to 30V during fixed times (from 10s up to 5min). To control accurately the contact force during charging, the tip spring constant was controlled periodically.

Concerning FDC measurements, two configurations were used. The first one consists of applying DC voltage on AFM tip, whereas the sample back-side is grounded. By acquisition of the resulting electrostatic force the sensitivity of EFDC to AFM probe characteristics and the force itself can be evaluated. In the same way, the influence of tip radius and spring constant were investigated to identify the best experimental conditions using silicon AFM tip with different coating materials (Pt, Co and diamond) provided by Bruker with references SCM-PIT/SCM-PIC, MESP and DDESP, respectively. For each tip used, the curvature radius (R c ) is determined using Scanning Electron Microscopy. The second configuration consists in measuring electrostatic forces induced on the AFM tip by the localized injected charges. These charges were previously injected in the same conditions as for KFM measurements. As FDC measurement is sensitive to piezoelectric system creeping [START_REF] Cappella | Force-distance curves by atomic force microscopy[END_REF], verification was done to insure EFDC is not affected by tip-sample distance range.

TRAPPED CHARGES PROBING BY KFM

KFM PRINCIPLE

Surface potential measurement by KFM is divided in two steps: a first scan to acquire surface topography followed by a second scan, at fixed lift distance, to measure electrostatic force in order to deduce surface potential [START_REF] Mélin | Scanning Probe Microscopy in Nanoscience and Nanotechnology[END_REF]. During the second step a voltage of the following form is applied to the AFM tip:

.

(

) 1
This applied voltage induces an electrostatic force on the tip at the same pulsation as the excitation: ,

where is the surface potential difference between the tip and the surface, to be measured, and C the tip-surface capacitance.

The surface potential difference between the tip and the surface is now determined by tuning the DC voltage to cancel the electrostatic force F. As the force depends on the first derivative of capacitance, the potential measurement is sensitive to parasitic capacitances which limit the lateral spatial resolution [START_REF] Jacobs | Resolution and contract in Kelvin Probe force microscopy[END_REF]. So, the lift distance should be as short as possible to limit parasitic capacitance and improve lateral resolution [START_REF] Zaghloul | A systematic reliability investigation of the dielectric charging process in electrostatically actuated MEMS based on Kelvin probe force microscopy[END_REF]. Moreover, in case of high surface potential difference, a high electrostatic force is induced, which introduces topography error (apparent height) and might screen potential measurement [START_REF] Ziegler | Feed-forward compensation of surface potential in atomic force Microscopy[END_REF]. This introduces limitation for space charge probing in case of large charge density. To improve the lateral resolution other variants of the KFM method based on Frequency Modulation (FM) instead of Amplitude Modulation (AM) can be used, minimizing sensitivity to parasitic capacitance [START_REF] Ziegler | Force gradient sensitive detection in liftmode Kelvin probe force microscopy[END_REF].

For space charge measurements, our objective is not only to record pictures with good spatial resolution based on electrostatic properties contrast, but also to provide procedures to estimate quantitatively and reliably the charge density. The main limitation of KFM is to achieve lateral resolution without in-depth resolution. This represents serious limitation, which can be overcome to some extend for number of applications, as example for interfaces investigation. Indeed, strategies have been implemented for obtaining a potential profile along the inter-electrode distance by using imbedded electrodes [START_REF] Silveira | Microscopic View of Charge Injection in an Organic Semiconductor[END_REF] with resolution well beyond the µm scale. This has been applied to semiconducting materials, and works are ongoing to extend this kind of study to dielectric materials.

CHARGE BUILD-UP INVESTIGATION

Charges were generated with the AFM tip in contact mode as described in the experimental part. During KFM measurement two maps were acquired. The first one consists in surface topography which is not influenced by charges distribution and the second one represents the resulting potential spot after charge injection. As the potential distribution presents a two dimensions axisymmetric shape, our interest is focused on spot cross-section. Figure 1 compares the potential spot crosssections obtained after charging the different SiO x N y layers (charging conditions are adjusted to have the same maximum potential). Even if the KFM surface potential emphasizes different behavior for each layer, the profiles exhibit the same Gaussian shape and three parameters can be extracted: (i) the Full-Width at Half Maximum (FWHM) which characterizes lateral spreading of charges, (ii) the maximum voltage with reference to surface potential without charges and (iii) the area under potential profile which represents in first approximation the charge density (this assumption is valid only for thin films). Comparison of potential maximum to the area under spot provides information on charge release processes. In particular, it is used to distinguish charges spreading on the layer surface from their drift in the dielectric bulk [START_REF] Morita | Microscale Contact Charging on Silicon Oxide in Characterization[END_REF]. Figure 1 emphasizes that each layer exhibits proper behavior. For γ=100 and γ=5 charges are stored close to the injection point (FWHM=0.9 and 4.6µm, respectively) whereas γ=10 layer exhibits strong lateral charges spreading (FWHM=37µm). Figure 2a shows the influence of injection time on KFM profile for γ=100 layer. One should recall here that γ=100 layer presents highly insulating properties, close to those of thermal SiO 2 layer, which actually limit the charge transport. No real shape modification is induced when increasing the charging time (Fig. 2a). Figure 2b represents the evolution of the spot profile parameters with charging time. A saturation of injected charge is observed with increasing charging time (results not shown), whereas the FWHM is quite constant. Quite the same behavior is observed for γ=5 whereas, in γ=10 layer, the charge spot broadens and grows without saturation effect with charging time [START_REF] Villeneuve-Faure | Kelvin force microscopy characterization of charging effect in thin a-SiOxNy:H layers deposited in pulsed plasma enhanced chemical vapor deposition process by tuning the Siliconenvironment[END_REF].

Thorough analysis of the obtained results provides information on the evolution of the charge cloud behavior and infers possible decay mechanisms and trapping processes [START_REF] Tzeng | Charge trapping properties at silicon nitride/silicon oxide interface studied by variable-temperature electrostatic force microscopy[END_REF]. However, more work needs to be done to further understand contrasts observed by KFM on these dielectric layers. In addition, quantitative determination of charges using KFM measurement is difficult due to complications inherent to KFM modeling [START_REF] Palleau | Numerical simulations for quantitative analysis of AFM electrostatic nanopatterning on PMMA by Kelvin Force Microscopy[END_REF]. In-depth localization of the injected charges is impossible. In spite of these limitations, surface potential measurements by KFM have become a standard for exploring injection, trapping and transport of charges in insulating thin films.

ELECTROSTATIC FORCE DISTANCE CURVE: A NEW WAY TO PROBE SPACE CHARGE

EFDC PRINCIPLE

Electrostatic Force Distance Curve method is based on FDC technique which consists in measuring the evolution of the interaction force between tip and sample as a function of their separation distance. As shown on Fig. 3, the procedure consists of three phases, an approach phase, a jump to contact due to the Van der Waals forces and a contact step. An applied bias on the AFM tip modifies the curve shape and more particularly the approach phase reflecting the fact that the electrostatic force is a long-range force. The electrostatic force can be extracted from these measurements by taking the difference in the force curves with and without potential application on the tip (inset in Fig. 3). It must be stressed that the electrostatic force does not vary as the square of the distance as one might expect of an electrostatic force. This is mainly due to the fact that the electrostatic configuration is not plane-plane but tip-plane (typical tip radius being between 5 nm and 120 nm). The EFDC is actually sensitive to tip characteristics, contrary to potential measurement by using KFM. Figure 4 shows the influence of the AFM probe spring constant and radius of curvature on the EFDC shape for 15V applied on tip. It is observed that the curve shape is greatly influenced by the curvature radius while the baseline (long range force) varies a little. Indeed, when the curvature radius decreases, the maximum force decreases, and the jump to contact becomes more abrupt. Moreover, the baseline is more influenced by cantilever spring constant. For the same curvature radius, the shape of the curve is the same as the baseline is more important for low spring constants.

The cantilever spring constant (k) depends on the geometrical parameters, such as its length and width. In particular, k decreases with cantilever surface increasing. Therefore, the base line represents the cantilever contribution to the electrostatic force at a distance far from the surface (typically the height of the tip, of the order of 10µm). It therefore appears that the shape of the curve is influenced by the curvature radius while the base line is more sensitive to the spring constant.

To characterize EFDC curve two parameters can be extracted: baseline which corresponds to the long-range force contribution and maximum force which is the force at jump to contact. Figure 5 shows the evolution of the maximum force (occurring at contact point) as a function of the square of the applied voltage. First of all, the maximum force increases approximately linearly with the square of the potential. The slope of the line provides a trend regarding the sensitivity of the technique as a function of tip characteristics. Obtained results highlight that sensitivity is improved in tips with large curvature radius. Moreover, cantilever spring constant influences the detection threshold (minimum force that can be detected). Indeed, the electrostatic force F e is related to cantilever deflection δ by the spring constant k according to the relation: .

For a fixed force, the deflection amplitude is divided by a factor of 10 if the spring constant is multiplied by a factor of 10. It will therefore be easier to detect small forces with low spring constant. In what follows, an AFM tip with a radius of curvature R c = 60 nm and stiffness coefficient of 0.44N/m was used. These results demonstrate that EFDC measurements can be the base of a highly sensitive technique for measuring electrostatic forces induced by the applied bias, however the method is strongly dependent on the probe characteristics. Recently, the electrostatic force was successfully modeled using finite volume technique [START_REF] Boularas | Multi-dimensional modelling of electrostatic forces over dielectric surface: influence of tip geometry and correlation with experiment[END_REF] which reinforces the suggestion that the EFDC method is a promising method to investigate injected charges.

PROBING OF CHARGE INJECTION BY EFDC

Charging was realized following the same protocol as for KFM studies to achieve comparison between the two techniques. Indeed, charge injection is carried out by bringing the tip into contact with the surface and applying a potential difference for a variable time. The electrostatic force is then measured at different sites on the surface as shown in Fig. 6. First, one can notice that the shape of the EFDC recorded over the charge spot has similar shape to that obtained with a bias on tip. As mentioned previously, the electrostatic force does not scale up with the distance as a square function. However, in this case, no offset of the force at long distance is detectable, indicating that the electrostatic force induced by the injected charge has attenuated action on the cantilever. This is probably due to the localization of the charges, at the origin of the forces, in a small spot. Moreover, EFDC shape and maximum force vary with the distance between the measurement point and the injection point. The maximum and the curvature of the EFDC decrease with the distance from the injection point. Consequently, EFDC contains more information than that provided by the maximum force. One can exploit this property to extract information about charge localization.

Considering the maximum force profile as a function of the position related to the injection point, a direct comparison with the potential profile measured by KFM on the same charge spot can be done as shown in the inset in Fig. 6. The profile FWHM is almost the same for the two techniques. This tends to show that the estimated width of the charge spot is not impeded by the resolution of the two methods and corresponds to the actual spot size. Therefore, it can be considered that the lateral spatial resolution of both methods is smaller than the spot size. Additional information is expected to arise after exploring the shape (curvature) of the EFDC. To extract quantitatively this information on the spatial distribution of charge density, a numerical model is required to describe accurately EFDC induced by hypothetical charge clouds. This 3D Finite Element Model, developed on COMSOL, is in progress having first validation obtained for the applied bias on AFM tip [START_REF] Boularas | Multi-dimensional modelling of electrostatic forces over dielectric surface: influence of tip geometry and correlation with experiment[END_REF]. Extension to the computation of forces induced by 3D-charge clouds stored in the dielectric layer is currently on the way. In the following, the influence of charge density and localization are investigated experimentally to highlight EFDC sensitivity. To reach this goal, various charge spots were formed in SiO x N y thin layers with different electrical properties (i.e. different values of the γ ratio). Charge localization was determined by KFM potential measurement (cf. previous part) to determine the maximum potential and FWHM. Different aspects are under investigation, notably the influence of maximum potential and charge lateral spreading.

To extract information from EFDC curve, electrostatic force as function of distance is fitted by a logistic law:

(2

)
Where z is the tip-sample distance, is the force when the tip is far from the surface, is the maximum force (tipsurface contact), characterizes the abruptness of the jump to contact ( =0.01 in our case) and p is a parameter describing the EFDC curvature.

Figure 7 represents EFDC obtained for the same injection conditions as for the results shown in Fig. 2 (different charging times under 20V in γ=100 layer). It is observed that the shape of the curve remains essentially unchanged with the increase of the injection time. Only the maximum force increases. The KFM measurements made on these charge spots show that the maximum potential increases whereas FWHM remains nearly constant, between 0.6 and 0.8µm (Fig. 2). As shown in inset in Fig. 7, the maximum force increases with the area under potential profile which proves the direct link between this parameter and amount of injected charge. The decrease in slope at long charging time is attributed to a sensible spot broadening.

Figure 8 compares EFDC measurements made on three charge spots having different width and the same potential maximum. The maximum and FWHM were determined on KFM profiles shown in Fig. 1. It is observed that the curve shape is modified when the spot size increases. It appears that the curvature is smoother when the spot size is increased. Indeed, parameter p is equal to 1.2, 0.65 and 0.45 for FWHM of KFM spot of 0.6µm, 6µm and 16.5µm respectively. Furthermore, the maximum force increases with the FWHM of KFM profile, which is due to an increase of the amount of stored charges. This is confirmed in a previous study [START_REF] Villeneuve-Faure | Towards 3D charge localization by a method derived from atomic force microscopy : the electrostatic force distance curve[END_REF] highlighting that curvature, represented by p, is linked to charges lateral spreading and maximum force is related to charge amount (potential maximum and width of charge distribution).

To address EFDC in depth sensitivity, electrostatic force measurements were performed on buried electrodes of different width. Table 2 compares maximum force and p for three different electrodes (width and applied potential) buried at a depth of 10nm, 50nm or 100nm-in SiOxNy. These results confirm that curvature is influenced by the lateral expansion of the charged region. Moreover, p has the same value for charge spot or electrode of 6µm-width.

Concerning maximum force, results emphasize that F m is sensitive to applied bias and electrode vertical localization. Indeed, electrostatic force decreases with electrode depth increasing. So these results demonstrate that EFDC is sensitive to in depth charges localization. Even if EFDC shape sensitivity to charges localization is demonstrated experimentally, modelling of electrostatic force appears essential in order to extract the various contributions. According to superposition principle, extraction of space charges profile from EFDC measurement is theoretically possible but not trivial. The first step is to develop accurate modelling of electrostatic force taking into account tip shape and to realize a parametric study on the influence of trapped charge features (density, distribution shape, in-depth localization…) on the EFDC. The second step will consist in the development of inverse method permitting charges density extraction from EFDC curves over charges spot (as on figure 4). The unicity of the solution is not guaranteed at present; however results provided by this method can be coupled e.g. to KFM to level off indetermination.

EFDC AND KFM RELATIVE MERIT FOR SPACE CHARGE MEASUREMENT

To compare KFM and EFDC methods, three aspects are compared: reproducibility, in depth sensitivity and charges extraction strategy.

Concerning measurement reproducibility, the protocol for charges injection was repeated five times at different location and EFDC and KFM were performed each time. Concerning KFM measurement a standard deviation of 8% for maximum potential and of 4% for FWHM was extracted. Concerning EFDC a standard deviation of 4% for maximum force and 7% for p was extracted [START_REF] Villeneuve-Faure | Towards 3D charge localization by a method derived from atomic force microscopy : the electrostatic force distance curve[END_REF]. So these two methods seem to have comparable reproducibility. This conclusion could be limited by the fact that EFDC is more sensitive to tip shape than KFM. But, as charges injection is strongly influenced by tip shape (the curvature radius controls the electric field at injection point) [START_REF] Makasheva | Dielectric charging by AFM in tip-to-sample space mode: overview and challenges in revealing the appropriate mechanisms[END_REF], the tip needs to be well characterized whatever measurement technique is used when localized charges injection is investigated.

Results presented above highlight that EFDC should provide more information than KFM in terms of space charge localizationthrough its sensitivity to in-depth charges localization. Concerning charges density extraction from KFM measurement, two kinds of models have been developed, taking into account or not the fact that KFM is not sensitive to image charges [START_REF] Palleau | Numerical simulations for quantitative analysis of AFM electrostatic nanopatterning on PMMA by Kelvin Force Microscopy[END_REF][START_REF] Charrier | Real versus measured surface potentials in scanning Kelvin probe microscopy[END_REF]. To extract charge densities from these models strong hypotheses are needed concerning charge distribution in the vertical direction, as e.g. homogeneous distribution or Gaussian distribution, which induce a 50% of variation on charges density value [START_REF] Charrier | Real versus measured surface potentials in scanning Kelvin probe microscopy[END_REF]. Concerning EFDC, charges density extraction from curve is not demonstrated, but results present here are promising. Even if the solution unicity is not guaranteed a comparison between results provided by KFM and EFDC could approach charges density in 3D.

Regarding practical implementation, the EFDC measurement is more time-consuming than the KFM one because the force distance curve must be acquired point-bypoint before data processing. This feature is not restrictive because dedicated software for data acquisition (which can be implemented in Peak-Force mode, for example) can solve this issue. The most restrictive aspect of EFDC method is the intrusiveness of measurement. Indeed, during electrostatic force measurement the grounded AFM tip is in contact with the sample surface which alters spot charges by inducing charges quantity decrease. Reducing contact force and time in contact permits to limit this effect. Moreover, the time between two consecutive measurements (typically 5min) should be increased to avoid discharging effect by the tip. This is restrictive for charges decay investigation [START_REF] Villeneuve-Faure | Kelvin force microscopy characterization of charging effect in thin a-SiOxNy:H layers deposited in pulsed plasma enhanced chemical vapor deposition process by tuning the Siliconenvironment[END_REF], but not crippling for the slow discharging dynamic observed classically in dielectric.

Even if the EFDC technique is somewhat heavier to implement than the KFM one, especially for drawing profiles from different positions on a charge spot, EFDC appears very promising for the study of the density of space charge due to its sensitivity to charges localization.

CONCLUSION

In this study KFM and EFDC methods were presented and compared as tools to probe space charges in thin dielectric layers at nanometer scale. The information gained is not as direct as that resulting from acoustic or thermal methods, based on perturbation of the equilibrium between electrostatic and elastic forces. However, the nanoscale resolution is not at reach at the present time with these classical methods. KFM provides quick and interesting results about charges localization at the surface of the sample; however the lack of sensitivity on in-depth localization is limiting. Contrary to KFM, EFDC is a point-by-point technique which probes directly the electrostatic forces induced by localized charges. The main advantages of this method compared to KFM are its sensitivity to charge distribution and the fact that EFDC can be easily modelled. These studies are currently underway aiming at quantifying the in-depth resolution of the EFDC. Indeed, the purpose is not only to image charges, but to obtain quantitative information on their density. For this, the development of a 3D model of the electrostatic force is required. Our recent works show a good match between the simulated and measured forces on biased tip, and the next and indispensable step will be to extract quantitatively the charges density.

Figure 1 .

 1 Figure 1. KFM profile over localized injected charges in different SiOxNy layers. Charging conditions have been chosen to provide the same maximum potential.

Figure 2 .

 2 Figure 2. (a) KFM profile after charging under 20V during different times in SiOxNy γ=100 layer. (b) Evolution of maximum potential Vm, area under peak and FWHM as function of injection time.

Figure 3 .Figure 4 .

 34 Figure 3. Influence of tip bias on the EFDC shape. In inset, resulting electrostatic force for 25V applied on tip.

Figure 5 .

 5 Figure 5. Evolution of the maximum force as a function of the square of applied bias V².

Figure 6 .

 6 Figure 6. EFDC measurements at different localization over injected charge spot (charging for 1min under 20V) in γ=100 SiOxNy layer. In inset, KFM and maximum electrostatic force profiles.

Figure 7 .Figure 8 .

 78 Figure 7. EFDC recorded after charging γ=100 SiOxNy layer with AFM tip under 20V for times indicated in the legend. In inset, evolution of the maximum force as a function of the area under KFM spot which represents in first approximation the charge density. 0.0 0.5 1.0 1.5 -4 -3 -2 -1 0 =100 -FWHM = 0.5µm =5 -FWHM = 5µm =10 -FWHM = 15.8µm Tip-sample distance (µm) Electrostatic force (nN)

Table 1 .

 1 Achieved, accessible and ultimate resolutions of classical space charge measurement techniques.

	Resolution	Achieved	Accessible	Ultimate
	Thermal	LIMM: 1µm	0.5µm	LIMM: 100nm
	Methods	TP: 0.6µm		TP: 60nm
	Acoustic
	Methods

Table 1 .

 1 Evolution of p and maximum force Fm of EFDC measured over electrode buried at different depths as function of electrode width and bias applied.

	Depth		Width 6µm 10V	Width 20µm	Width 6µm 4V
	(nm)		applied	10V applied	applied
	100	p	0.62	0.42	0.61
			-24	-31	-6
	50	p	0.6	0.42	0.59
			-26	-35	-7.8
	10	p	0.69	0.46	0.65
			-63	-65	-12.4

ACKNOWLEDGMENT

This work has been supported by Région Midi-Pyrénées under contract Reference 10051265.


Bishnu Sarker

David W Ritchie

Sabeur Aridhi
email: sabeur.aridhi@loria.fr

GrAPFI: predicting enzymatic function of proteins from domain similarity graphs

Keywords: Protein function annotation, Protein network, EC annotation, Label propagation, Domain similarity graph, GrAPFI, K-nearest neighbor

Background: Thanks to recent developments in genomic sequencing technologies, the number of protein sequences in public databases is growing enormously. To enrich and exploit this immensely valuable data, it is essential to annotate these sequences with functional properties such as Enzyme Commission (EC) numbers, for example. The January 2019 release of the Uniprot Knowledge base (UniprotKB) contains around 140 million protein sequences. However, only about half of a million of these (UniprotKB/SwissProt) have been reviewed and functionally annotated by expert curators using data extracted from the literature and computational analyses. To reduce the gap between the annotated and unannotated protein sequences, it is essential to develop accurate automatic protein function annotation techniques. Results: In this work, we present GrAPFI (Graph-based Automatic Protein Function Inference) for automatically annotating proteins with EC number functional descriptors from a protein domain similarity graph. We validated the performance of GrAPFI using six reference proteomes in UniprotKB/SwissProt, namely Human, Mouse, Rat, Yeast, E. Coli and Arabidopsis thaliana. We also compared GrAPFI with existing EC prediction approaches such as ECPred, DEEPre, and SVMProt. This shows that GrAPFI achieves better accuracy and comparable or better coverage with respect to these earlier approaches. Conclusions: GrAPFI is a novel protein function annotation tool that performs automatic inference on a network of proteins that are related according to their domain composition. Our evaluation of GrAPFI shows that it gives better performance than other state of the art methods. GrAPFI is available at https://gitlab.inria.fr/bsarker/bmc_grapfi.git as a stand alone tool written in Python.

Background

Proteins are long sequences of amino acids that form the basis of life and plays vital role in all living organism through out the entire life-cycle. Proteins perform various functions in our body that needs to be understood to understand life, disease processes and guiding drug discovery efforts to combat the diseases. Due to the tremendous advancement in amino-acid sequencing technologies, it is now possible to sequence bulk amount of proteins in a rapid and affordable manner.

Therefore, the number of protein sequences accumulating in public databases is rising at an unprecedented rate. This huge quantity of data calls for further exploitation and enrichment and it presents many challenges for biologists as well as computer scientists in annotating the functional properties of protein sequences. The UniProt knowledge base (UniProtKB) [START_REF] Aridhi | UniProt: a hub for protein information[END_REF] is one of the most comprehensible protein databases as well as the largest public sequence database currently. It is divided into two main components: (i) the UniProtKB/Swiss-Prot database which contains protein sequences with reliable functional annotation of the protein sequences that has been reviewed by expert bio-curators, and (ii) the UniPro-tKB/TrEMBL database that stores unannotated and unreviewed sequences. Thus, for all proteins in UniProtKB, we have the primary amino acid sequence as well as some further information such as InterPro domain definitions which may have been identified from families of similar sequences or 3D protein structures.

The UniProt curators annotate UniProtKB/TrEMBL sequences using two complementary systems. The first, called UniRule, uses a large list of "if-then" rules. These rules have been generated manually, which is both a laborious and time consuming process. The rules in UniRule are generally very reliable but their coverage is low [START_REF] Gattiker | Automated annotation of microbial proteomes in SWISS-PROT[END_REF]. The second system is called Statistical Automatic Annotation System (SAAS), and was developed to support the labourintensive UniRule system [START_REF] Kretschmann | Automatic rule generation for protein annotation with the C4.5 data mining algorithm applied on SWISS-PROT[END_REF]. Automatic annotation rules are generated in SAAS using the annotations of the Swis-sProt sequences and a decision tree algorithm [START_REF] Quinlan | Induction of decision trees[END_REF]. Other approaches exist for automatic protein function annotation. In particular, a number of techniques for predicting Enzyme Commission (EC) numbers that exploit protein structural similarities have been discussed in [START_REF] Dobson | Predicting enzyme class from protein structure without alignments[END_REF][START_REF] Yang | The i-tasser suite: protein structure and function prediction[END_REF][START_REF] Chioko | Prediction of detailed enzyme functions and identification of specificity determining residues by random forests[END_REF]. Many approaches based on sequence similarity have also been discussed in [START_REF] Rahman | EC-BLAST: a tool to automatically search and compare enzyme reactions[END_REF][START_REF] Kumar | Eficaz2. 5: application of a high-precision enzyme function predictor to 396 proteomes[END_REF][START_REF] Quester | EnzymeDetector: an integrated enzyme function prediction tool and database[END_REF][START_REF] Yu | Genome-wide enzyme annotation with precision control: Catalytic families (CatFam) databases[END_REF]. Several machine learning methods have also been studied extensively in [START_REF] Chioko | Prediction of detailed enzyme functions and identification of specificity determining residues by random forests[END_REF][START_REF] Des Jardins | Prediction of enzyme classification from protein sequence without the use of sequence similarity[END_REF][START_REF] Li | Svm-prot 2016: a web-server for machine learning prediction of protein functional families from sequence irrespective of similarity[END_REF][START_REF] Huang | Accurate prediction of enzyme subfamily class using an adaptive fuzzy k-nearest neighbor method[END_REF][START_REF] Lu | ECS: an automatic enzyme classifier based on functional domain composition[END_REF][START_REF] Nasibov | Efficiency analysis of KNN and minimum distance-based classifiers in enzyme family prediction[END_REF][START_REF] Li | DEEPre: sequence-based enzyme EC number prediction by deep learning[END_REF][START_REF] Sarker | Functional Annotation of Proteins Using Domain Embedding Based Sequence Classification[END_REF][START_REF] Shen | Ezypred: a top-down approach for predicting enzyme functional classes and subclasses[END_REF][START_REF] Volpato | Accurate prediction of protein enzymatic class by n-to-1 neural networks[END_REF].

Recently, the notion of network science [START_REF] Barabási | Linked: The new science of networks[END_REF] has attracted great attention across many scientific communities. Network science has become a multidisciplinary area of research due to its ability to describe complex systems. It has found applications in many real-world scenarios from banking and the internet to modeling the human brain. Several approaches for annotating protein function have used network science and neighborhood based techniques to extract functional information from proteinprotein interaction (PPI) networks and Gene Ontology terms [START_REF] Schwikowski | A network of protein-protein interactions in yeast[END_REF][START_REF] Zhao | An efficient method for protein function annotation based on multilayer protein networks[END_REF][START_REF] Hishigaki | Assessment of prediction accuracy of protein function from protein-protein interaction data[END_REF][START_REF] Chua | Exploiting indirect neighbours and topological weight to predict protein function from protein-protein interactions[END_REF][START_REF] Nabieva | Whole-proteome prediction of protein function via graph-theoretic analysis of interaction maps[END_REF]. A particular feature of biological networks is that they often require specialist biological knowledge to fully understand and exploit the network.

The following methods are widely used for predicting Enzyme Commission (EC) numbers using a variety of approaches based on machine learning, sequence encoding, functional domain similarity, and structural similarity. A deep learning based approach called DEEPre [START_REF] Li | DEEPre: sequence-based enzyme EC number prediction by deep learning[END_REF] predicts EC numbers putting together multiple tools and techniques including PSI-Blast [START_REF] Altschul | Gapped blast and psi-blast: a new generation of protein database search programs[END_REF], HMMER [START_REF] Finn | HMMER web server: interactive sequence similarity searching[END_REF], Convolutional and Recurrent Neural Networks, and sequence encoding using position specific scoring matrix (PSSM) to perform dimensionality uniformization, feature selection, and classification model training. In recent years, deep learning has been applied in many computational biology and healthcare prediction tasks and achieved state-of-theart performance. However, deep learning approaches can suffer from interpretability issues which is necessary in medical research and clinical decision-making [START_REF] Che | Interpretable deep models for ICU outcome prediction[END_REF].

EzyPred [START_REF] Shen | Ezypred: a top-down approach for predicting enzyme functional classes and subclasses[END_REF] is a k-nearest neighbor based method that adopts a top-down approach for predicting main class and sub-class of EC number. EzyPred works based on protein sequences only to perform the annotation task. It starts by predicting whether or not an input protein sequence is an enzyme. Then, EzyPred proceeds by predicting its main EC class and subclass. EzyPred uses pseudo amino acid composition [START_REF] Chou | Pseudo amino acid composition and its applications in bioinformatics, proteomics and system biology[END_REF] and functional encoding by exploiting functional and evolutionary information of proteins. Based on two features, EzyPred proposes a modified K-nearest neighbor classifier called OET-KNN (Optimized Evidence-Theoretic K-Nearest Neighbour). Although EzyPred performs well in terms of accuracy, it predicts only the first two digits of a four-digit EC number. Thus, its predictions are not very specific.

A machine learning based approach called SVM-Prot that uses support vector machine (SVM) for classification is proposed in [START_REF] Cai | Svm-prot: web-based support vector machine software for functional classification of a protein from its primary sequence[END_REF][START_REF] Cai | Enzyme family classification by support vector machines[END_REF][START_REF] Cai | Predicting enzyme subclass by functional domain composition and pseudo amino acid composition[END_REF]. And in 2016 [START_REF] Li | Svm-prot 2016: a web-server for machine learning prediction of protein functional families from sequence irrespective of similarity[END_REF], the performance of SVMProt is improved by adding two more classifiers: 1) K-Nearest Neighbor (KNN) and 2) Probabilistic Neural Networks (PNN). This approach uses important physicochemical properties such as molecular weight, polarity, hydrophobicity, surface tension, charge, normalized van der Waals volume, polarizability, secondary structure, solvent accessibility, solubility, and the numbers of hydrogen bond donors and acceptors in side chain atoms to transform protein sequences into numerical feature representations.

A structure-based protein function annotation approach called COFACTOR is described in [START_REF] Yang | The i-tasser suite: protein structure and function prediction[END_REF][START_REF] Roy | Cofactor: an accurate comparative algorithm for structure-based protein function annotation[END_REF]. The updated version of COFACTOR [START_REF] Zhang | Cofactor: improved protein function prediction by combining structure, sequence and protein-protein interaction information[END_REF] combines information about protein structure and sequence homologs along with Protein-Protein Interaction (PPI) networks to form a hybrid model for jointly predicting GO terms, EC numbers, and ligand-binding .

EFICAz [START_REF] Kumar | Eficaz2. 5: application of a high-precision enzyme function predictor to 396 proteomes[END_REF][START_REF] Tian | Eficaz: a comprehensive approach for accurate genome-scale enzyme function inference[END_REF][START_REF] Arakaki | Eficaz 2: enzyme function inference by a combined approach enhanced by machine learning[END_REF] presents a method for Enzyme Function Inference by Combined Approach. EFICAz combines predictions from four different methods using (i) recognition of functionally discriminating residues (FDRs) in enzyme families obtained by the authors' "CHIEFc" procedure (Conservation-controlled HMM Iterative procedure for Enzyme Family classification), (ii) pairwise sequence comparison using a family-specific sequence identity threshold, (iii) recognition of FDRs in Multiple Pfam enzyme families, and (iv) recognition of multiple Prosite patterns of high specificity.

In ECPred [START_REF] Dalkiran | ECPred: a tool for the prediction of the enzymatic functions of protein sequences based on the EC nomenclature[END_REF], the authors describe a hierarchical prediction model. The model starts by predicting if a query sequence is an enzyme or non-enzyme. Once the query sequence is found to be an enzyme, ECPred predicts the main class to which the query sequence belongs. In the similar fashion, it follows the hierarchy of the EC numbering system to find the sub-class, sub-sub-class and sub-sub-sub-class. ECPred learns independent classifiers for 858 EC classes including 6 main classes, 55 subclass classes, 163 sub-subclass classes and 634 sub-sub-sub classes. The independent predictors that make up ECPred are SPMap, BLAST-kNN and Pepstats-SVM which are based on sub-sequences, sequence similarities, and the physico-chemical features of amino acids, respectively.

In this paper, we give a complete description of our novel graph based annotation approach(GrAPFI) [START_REF] Sarker | Exploiting Complex Protein Domain Networks for Protein Function Annotation[END_REF], and we present an extended experimental analysis using test data from six popular reference proteomes from UniPro-tKB/SwissProt. GrAPFI builds network of proteins using domain and family information and performs neighborhood based label propagation for function annotation.

Although, similar to EZYPred [START_REF] Shen | Ezypred: a top-down approach for predicting enzyme functional classes and subclasses[END_REF] and SVM-Prot-KNN [START_REF] Li | Svm-prot 2016: a web-server for machine learning prediction of protein functional families from sequence irrespective of similarity[END_REF], GrAPFI is a neighborhood based classification technique, it uses different features and different inference mechanism that explores the network. GrAPFI uses Inter-Pro signatures as domain information and label propagation over a weighted undirected graph built on proteins using their domain composition. Unlike COFACTOR [START_REF] Zhang | Cofactor: improved protein function prediction by combining structure, sequence and protein-protein interaction information[END_REF], GrAPFI explores the functional domain architectures extracted from protein sequence instead of protein secondary structure and direct sequence homology. COFAC-TOR includes network information using PPI whereas GrAPFI builds network using jaccard similarity index between the domain composition of proteins. In contrast to ECPred and DEEPre [START_REF] Li | DEEPre: sequence-based enzyme EC number prediction by deep learning[END_REF] which are deep learning based method and learns class specific models for different classes, GrAPFI performs label propagation over weighted protein network to select the best EC annotation. ECPred learns 858 independent classifiers where as for GrAPFI, once the protein network is build, it's ready for the inference of EC number using domain composition of the query protein.

We compare the performance of GrAPFI with the recently published ECPred approach. Along with ECPred, we also present the accuracy for DEEPre and SVMProt as representative examples of other state of the art EC number prediction approaches. Our analysis shows that GrAPFI gives better annotation performance than these earlier approaches.

Results

Data preparation

We have collected 262,564 proteins from the March 2018 release of Uniprot-KB/SwissProt [START_REF] Aridhi | UniProt: a hub for protein information[END_REF] database using the following rules: (i) each protein must contain at least one InterPro signature and (ii) must be annotated with at least one EC annotation. After getting the protein data from each of the proteins, we have extracted InterPro domain composition and EC annotations. Then, we built the protein network as described in "EC annotation performance analysis" section. Each protein forms its own vertex. we did not preprocess training data to remove redundancy. Rather, while performing annotation, it ignores the same protein if it appears in the neighborhood. For example, for a query protein q, GrAPFI will collect the neighbors satisfying a maximum jaccard similarity score. When the maximum jaccard similarity is set to less than 1.0, GrAPFI omits the neighbors with exact same domain composition.

The training network covers 25 level-2, 31 level-3 and 408 level-4 EC classes from 41,618 oxidoreductases, 70,530 transferases, 100,027 hydrolases, 14,677 lyases, 25,551 isomerases, and 29,735 ligases which are linked using 10,866 InterPro signatures.

In the training network, there are 1) 4.3% of the proteins are single-domain proteins i.e. proteins having only one domain in their domain composition (Fig. 1a,c), 2) 5.7% of the proteins have more than one EC numbers assined with them (Fig. 1b,c), and 3) Around 15% of the training nodes have incomplete EC annotations i.e. the EC numbers assigned with these proteins do not have all four digits. In the Fig. 1a, we show the distribution of EC numbers against domain composition. There are 13713 unique domain compositions in the training data. In the X-axis we put the domain compositions and along Y-axis we have the number of different EC annotations found for each domain composition sorted in descending order. It is evident from the figure that few of the domain compositions contain significantly higher number of EC numbers. For example, for some domain composition, there are more than 50 EC numbers found in the training data. We also show the distribution of domain compositions per EC number i.e. the different domain compositions found for each EC annotation shown in Fig. 1b. There are many cases when a higher number of domain compositions are mapped to a single EC. For example, in some cases, it is around 500 distinct domain compositions found against a single EC number. In essence, these two distributions reflect the dominance of many-to-many relationship between domain composition and EC annotation in the training data.

To valdiate GrAPFI, we used six popular reference proteomes from Uniprot-KB/SwissProt to as test set. The reference proteomes are the following: (i) Rattus norvegicus (UP000002494) containing 1,953 proteins, (ii) Mus musculus (UP000000589) containing 3,682 proteins, (ii) Saccharomyces cerevisiae (UP000002311) containing 1,581 proteins, (iv) Homo sapiens (UP000005640) containing 3,843 proteins, and (v) Arabidopsis thaliana (UP000006548) containing 5,352 proteins. (vi) E. Coli (UP000000625) containing 1465 proteins. For each of the reference proteomes, we collected the InterPro domains and EC labels from Uniprot-KB/Swissprot. We kept only the proteins which have at least one InterPro domain and are annotated with a single EC number.

To prepare the COFACTOR benchmark dataset, we used the 318 protein sequences published in [START_REF] Zhang | Cofactor: improved protein function prediction by combining structure, sequence and protein-protein interaction information[END_REF], and we ran InterProScan [START_REF] Jones | Interproscan 5: genome-scale protein function classification[END_REF] on these sequences to get their InterPro domain signatures. We only used IntePro domain signatures for the purpose of EC annotation.

EC annotation performance analysis

To validate the annotation performance of GrAPFI, we computed the accuracy, Macro-precision, Macro-recall, and Macro-F1 score at different levels of EC number. For each query sequence, we picked the top ranked annotation only. The validation method we have used is similar to leave one out cross validation. For each proteomes, when annotating a protein, we have removed that protein from the training set so that a direct mapping is not present. we also present a 10-fold cross validation for enzyme vs. non-enzyme classification (Fig. 2a andb). The following formula (as used in [START_REF] Li | DEEPre: sequence-based enzyme EC number prediction by deep learning[END_REF][START_REF] Sarker | Functional Annotation of Proteins Using Domain Embedding Based Sequence Classification[END_REF]) were used to compute the evaluation scores:

accuraccy(y, y ) = 1 N N-1 i=0 1(y i = y i ),
Here, y and y are the list of ground truths and predicted annotations. The accuracy is computed for each level of EC annotation. As the problem is a multi-class classification problem, we computed class-wise Macro-precision, Macro-recall, and Macro-F1 score as follows:

Macro -precision(y, y ) = 1 |M| l∈M precision(y l , y l ), Macro -recall(y, y ) = 1 |M| l∈M recall(y l , y l ), Macro -F1(y, y ) = 1 |M| l∈M F1 measure(y l , y l ),
Here, y l is the part of y with the label l and y l is the part of y with label l. And M is the set of classes. In general the precision, recall, and F1-Measure are computed as follows when two sets A and P are given:

precision = |A ∩ P| |P| , recall = |A ∪ P| |A| , F1 -measure = 2 × precision × recall precision + recall .
Here, A is the set of ground truths and P is the set of predictions. As EC numbers are hierarchical with 4 levels, we Fig. 2 The precision, recall, F1, accuracy and coverage score for various minimum Jaccard similarity index for the Enzyme vs. Non-enzyme classification. a and b show the performance for upper similarity index of 1 and less than 1 respectively report level-wise precision, recall and F1-measure. Level-1 denotes main class, level-2 denotes sub-class, level-3 denotes sub-sub-class and level-4 denotes sub-sub-sub class. We also report coverage which is calculated according to Coverage = M/T, where T is the total number of proteins in the test set and M is the number of proteins for which at least one EC is predicted. For each query sequence, we consider the top-most annotation. For evaluation purposes, we split the 4-digit EC annotation into its constituent parts. Then, for level-1 we consider first digit, for level-2 we take first 2 digits, for level-3 we take first 3-digits and finally for level-4 we take all four digits together.

For the validation dataset, GrAPFI was run by setting different minimum jaccard similarity index ranging from 0.05 to 0.5, and setting an upper limit of the similarity to 1 or less than 1.

Figure 3a to h show the GrAPFI performance for the reference proteome of A. thaliana for various min. Jaccard similarity indices. Similarly, Figs. 4, 5, 6, 7 and 8 show the performance of GrAPFI for the reference proteomes of Homo sapiens, S. cerevisiae, Mus musculus, Rattus norvegicus and E. Coli, respectively. In Fig. 9, we also show the performance of GrAPFI on COFACTOR benchmark dataset for various Jaccard domain similarity index ranging from 0.05 to 0.5, and setting an upper limit of the similarity to 1 and less than 1.

In these figures, we show the annotation accuracy (Y axis) against different Jaccard similarity thresholds (X axis) for the respective proteomes. We have considered similarity thresholds ranging from 0.05 to 0.5 as the annotation coverage falls significantly after 0.5. For each of the thresholds, we present the accuracy for EC-1, EC-2, EC-3 and EC-4 digit prediction shown in green, blue, orange and black color respectively. Along with accuracy, we also present the coverage of annotation (red curve). For each of the figures, we have two parts. The first part shows the accuracy and the coverage considering only the neighbors who have a Jaccard similarity of less than 1. The second part considers the Jaccard similarity of less than 1. It can be seen from these figures that GrAPFI performs very well for all of the cases with a good coverage.

To compare GrAPFI with other state of the art methods, we considered three machine learning based methods, namely ECPred [START_REF] Dalkiran | ECPred: a tool for the prediction of the enzymatic functions of protein sequences based on the EC nomenclature[END_REF], DEEPre [START_REF] Li | DEEPre: sequence-based enzyme EC number prediction by deep learning[END_REF], and SVMProt [START_REF] Li | Svm-prot 2016: a web-server for machine learning prediction of protein functional families from sequence irrespective of similarity[END_REF]. The performance is compared based on the COFACTOR [START_REF] Zhang | Cofactor: improved protein function prediction by combining structure, sequence and protein-protein interaction information[END_REF] benchmark having 318 sequences. The SVMProt prediction results cover three different conditions: (i) using SVM only, (ii) using KNN only, and (iii) using SVM, KNN and PNN combined. Figure 10a shows the performance analysis for EC level-1 and EC level-2 prediction. The results presented here are achieved using a lower Jaccard similarity index of 0.3 and upper similarity index of 1.0. A much lower similarity threshold brings false positives that significantly reduce the accuracy. Based on the obtained results, a similarity index of 0.3 achieves a good trade-off between accuracy and coverage.

Because not all of the methods can make predictions for all four EC levels, we compared GrAPFI only with ECPred and DEEPre for 4-digit EC numbers as shown Fig. 3 The accuracy and coverage for various minimum Jaccard similarity index for the A. thaliana reference proteome. a to d and e to h show the performance for upper similarity index of 1 and less than 1 respectively in Fig. 10b. In Fig. 10c, we show the annotation coverage of the methods considered here. It shows that ECPred has superior coverage compared to other methods. The reason GrAPFI fails to achieve highest coverage is due to the fact that it is a neighborhood based annotation method. GrAPFI performs label propagation by filtering out weakly linked neighbors determined by a minimum similarity threshold. Due to this filtering action, for few cases, GrAPFI fails to suggest any appropriate annotation for query proteins. This reduces the total annotation Fig. [START_REF] Quinlan | Induction of decision trees[END_REF] The accuracy and coverage for various minimum Jaccard similarity index for the Homo Sapiens reference proteome. a to d and e to h show the performance for upper similarity index of 1 and less than 1 respectively Fig. 5 The accuracy and coverage for various minimum Jaccard similarity index for the S. Cerevisiae reference proteome. a to d and e to h show the performance for upper similarity index of 1 and less than 1 respectively coverage. However, on the other hand, GrAPFI increases the accuracy by considering strongly linked neighbors. As shown in Fig. 10b andc, GrAPFI has better accuracy compared to ECPred and DEEPre, but it gives slightly less coverage than ECPred.

Enzyme vs. non-enzyme classification

GrAPFI can be used in Enzyme vs. Non-enzyme classification task in a similar fashion as described in above section. However, the training graph must include non-enzyme proteins. To experiment with enzyme vs. non-enzyme Fig. 6 The accuracy and coverage for various minimum Jaccard similarity index for Mus Musculus reference proteome. a to d and e to h show the performance for upper similarity index of 1 and less 1 Fig. 7 The accuracy and coverage for various minimum Jaccard similarity index for the Rattus Norvegicus reference proteome. a to d and e to h show the performance for upper similarity index of 1 and less than 1 respectively Fig. 8 The accuracy and coverage for various minimum Jaccard similarity index for the E.Coli reference proteome. a to c and d to f show the performance for upper similarity index of 1 and less than 1 respectively Fig. 9 The accuracy and coverage for various minimum Jaccard similarity index for the COFACTOR benchmark proteins. a and b show the performance for upper similarity index of 1 and less than 1 respectively classification, To evaluate the method, we have used a well defined dataset of enzyme and non-enzyme proteins curated from UniprotKB [START_REF] Aridhi | UniProt: a hub for protein information[END_REF]. This dataset is called "NEW" and was constructed as described in [START_REF] Li | DEEPre: sequence-based enzyme EC number prediction by deep learning[END_REF]:

1. The SWISS-PROT (released on September 7, 2016) database was separated into enzymes and non-enzymes based on their annotation. 2. To guarantee uniqueness and correctness, enzyme sequences with more than one set of EC numbers or incomplete EC number annotation were excluded. 3. To avoid fragment data, enzyme sequences annotated with 'fragment' or with less than 50 amino acids were excluded. Enzyme sequences with more than 5000 amino acids were also excluded. 4. Redundancy bias is removed using CD-HIT [START_REF] Fu | Cd-hit: accelerated for clustering the next-generation sequencing data[END_REF] with 40% similarity threshold to sift the raw dataset, resulting in 22,168 low-homology enzyme sequences. 5. To construct the non-enzyme part, 22,168 non-enzyme protein sequences were randomly collected from the SWISS-PROT (released on September 7, 2016) non-enzyme part, which were also subject to the above (ii-iv) steps. Thus the original dataset contains 22,168 enzymes and an equal number of non-enzymes. The dataset contains the protein sequences along with their respective EC annotations. We have run Inter-ProScan5 [START_REF] Jones | Interproscan 5: genome-scale protein function classification[END_REF] to identify the domains contained in the sequences. Later, with the domain information, we have built the training graph. This graph contains 40040 proteins with 54% enzymes and 46% non-enzymes connected based on their domain composition.

To evaluate the annotation performance, we present 10fold cross validation on the training graph and average macro-precision, macro-recall, macro-F1 scores are computed for various jaccard similarity indices. The result shows performance of enzyme vs. non-enzyme classification only. The experimental outcomes are shown in Fig. 2a andb.

It is evident from the experimental outcome that GrAPFI can distinguish enzyme and non-enzyme proteins with a a good score in all evaluation metrics. However, the coverage goes down as we move towards higher similarity thresholds. One of the things to be noted that considering exact similarity match does not change the performance significantly as can be seen in Fig. 2b.

Discussions

Here, we explore new ways of connecting proteins. The proteins are connected based on domains that are potentially linked to the protein functions. This eventually means that GrAPFI is biologically meaningful approach. One of the major advantages of using GrAPFI to annotate proteins is that it produces explainable high quality annotations with a relatively simple annotation pipeline. The potential is evident from the experimental results. Although GrAPFI performs well, there are few drawbacks of using GrAPFI. For example, GrAPFI works on domain composition that can be achieved using another tool. GrAPFI can not be used with proteins without domain information. And also for the proteins with single domain, in most cases, GrAPFI fails to find appropriate annotation. The reason for this failure is that for a single domain protein, it is highly unlikely that there will be any high quality neighbors that can share annotations that eventually left the protein without any labels or wrong one. In any case, if GrAPFI fails to find an annotation, it is possible to identify the reason behind the failure and it restricts itself from predicting any annotation. This attitude reduces the false positives. However, from the experiment, it is evident that GrAPFI performs with high annotation coverage. Unlike other hierarchical classification models like ECPred [START_REF] Dalkiran | ECPred: a tool for the prediction of the enzymatic functions of protein sequences based on the EC nomenclature[END_REF] and DEEPre [START_REF] Li | DEEPre: sequence-based enzyme EC number prediction by deep learning[END_REF], GrAPFI does not learn model for every class. Instead, it builts a giant network of proteins and apply label propagation for each query proteins. The described approach could easily be distributed in order to handle large protein databases. The method is scalable for larger dataset using big data processing frameworks like Hadoop/Spark. We therefore aim to extent GrAPFI to use a distributed processing framework for the large scale annotation of the entire UniProtKB/TrEMBL database. Moreover, there is still scope of improvement specially for level-3 and level-4 predictions. As a future plan, we envision to improve the method for more precise predictions and also to apply the similar approach for protein function annotation using Gene Ontology Terms.

Conclusion

In this paper, We have extended and validated GrAPFI [START_REF] Sarker | Exploiting Complex Protein Domain Networks for Protein Function Annotation[END_REF], a novel network based approach for automatic protein function annotation using the domain composition of the proteins. Pairs of proteins in the network are linked based on their jaccard similarity coefficient using InterPro domain composition.

Our neighborhood based label propagation algorithm was applied to the network in order to propagate annotations from reviewed proteins to non-reviewed query proteins. This approach was validated using six popular reference proteomes from UniProtKB/SwissProt. We also compared GrAPFI results with those of ECPred, DEEPre, and SVMProt as examples of state of the art EC prediction approaches using the COFACTOR dataset. This comparison shows that GrAPFI achieves better accuracy and comparable or better coverage with respect to these earlier approaches.

Methods

GrAPFI combines the notion of protein domain similarity with a graph neighborhood inference technique for automatic EC number annotation. More specifically, the functional annotations of reviewed proteins in SwissProt are used to predict those of non-reviewed proteins in TrEMBL using label propagation on a complex network representation of protein sequence data. The GrAPFI algorithm first constructs an undirected weighted graph of the proteins using the domain composition of the reviewed proteins. Then, given an non-reviewed protein, a label propagation algorithm is applied to the protein graph in order to infer appropriate annotations.

Notation

In this section, we first present some definitions and notations used in the paper.

Graph: A graph is a collection of objects denoted as G = (V , E), where V is a set of vertices/nodes and E ⊆ V × V is a set of edges.

Weighted Graph: A weighted graph is a graph which is represented as a three tuple G = (V , E, W ) where:

• V is a set of nodes, • E ⊆ V × V is a set of edges,
• W is a weight matrix where each cell W uv represents a numerical weight of the edge (u, v) ⊆ E.

Labeled Graph: A labeled graph is a graph which is represented as a four tuple G = (V , E, L, I) where:

• V is a set of nodes, • E ⊆ V × V is a set of edges, • L is a set of labels, • I : V ∪ E -→ L is a labeling function. Directed Graph: A Directed graph G = (V , E) is a col- lection of
objects where V is a set of vertices/nodes and E ⊆ V × V is a set of edges with ordered pair of vertices

(u, v) such that (u -→ v) ∈ E.
Undirected Graph: An undirected graph is a collection of objects denoted as G = (V , E), where V is a set of vertices/nodes and E ⊆ V × V is a set of edges with unordered vertices

(u, v) such that if (u -→ v) ∈ E exists then (v -→ u) ∈ E must exists.
Neighbors: The neighbors of a node u are defined as

N(u) = {v|(u, v) ∈ E}.
Degree: The degree of a node in a graph is the number of edges which touch it. The degree of a node u in a graph G is denoted deg(u) = N(u).

Average Degree: The average degree of a graph G = (V , E) is a measure of how many edges are in the set E compared to number of vertices in the set V. The average degree of a graph G = (V , E) is defined by Avgdeg = 2|E|/|V |.

Graph construction

We present here a novel way of connecting protein sequences using their associated InterPro domains. Domains may be considered as natural building blocks of proteins. Due to evolution, protein domains may have gone through changes such as duplication, fusion, recombination to produce proteins with distinct structures and functions [START_REF] Kummerfeld | Protein domain organisation: adding order[END_REF]. Here, each node of the graph represents a protein, while a link between two nodes means that the proteins exhibit a given level of domain similarity. Thus, each node u is identified by a set of labels L(u), has a set of neighbours N(u), and for every neighbour v it has an associated weight W u,v . The overall aim is to propagate labels (i.e. annotations) from nodes having labels to similar nodes that lack labels.

To illustrate the construction of the protein graph, let us consider five proteins with symbolic names P1, P2, P3, P4, and P5. Let us assume that these proteins are composed of domains D1 = (d1, d2, d3, d4), D2 = (d1, d3, d5), D3 = (d1, d2, d10), D4 = (d5, d6, d1), and D5 = (d4, d1, d10, d40, d7, d9, d12, d52, d100), respectively. Domain composition of a protein is the set of domains found in a protein sequence and considered irrespective of order of appearance in the sequence. For example the domain information in D1 = (d1, d2, d3, d4) can be used in any other order D1 = (d1, d4, d3, d2). Therefore, the composition is not strictly linear. The overlapping of domains are not considered as long as the overlapped domains has a new domain identification.

It is then evident that proteins P1 and P2 contain two domains d1 and d3 in common. Therefore, proteins P1 and P2 may be linked and the number of shared domains may serve as a link weight given by

W P1,P2 = |(d1, d2, d3, d4) ∩ (d1, d3, d5)| = |(d1, d3)| = 2.
In a similar way, proteins P1 and P5 may be linked with a link weight of |(d1, d2, d3, d4) ∩ (d4, d1, d10, d40, d7, d9, d12, d52, d100)| = |(d1, d4)| = 2. In both cases, the link weight is 2. However, the link weight computed in this way does not reflect the relative strength of the relationship among the proteins. More specifically, in the first case the two proteins have |(d1, d2, d3, d4) ∪ (d1, d3, d5)| = |(d1, d2, d3, d4, d5)| = 5 different domains, of which two are shared. In the second case, there are |(d1, d2, d3, d4) ∪ (d4, d1, d10, d40, d7, d9, d12, d52, d100)| = 11 different domains of which again two are shared. Although two domains are shared in each case, P1 is intuitively more aligned with P2 than P5. Therefore, instead of using the above raw similarity score, we used the Jaccard similarity index, or Jaccard similarity coefficient, to reflect better the similarity in composition. This is calculated as |A∩B| |A∪B| , where A and B are the two sets of constituent domains. Using the Jaccard similarity index, the link weights for P1 and P2 are calculated as

W P1,P2 = |(d1, d2, d3, d4) ∩ (d1, d3, d5)| |(d1, d2, d3, d4) ∪ (d1, d3, d5)| = |(d1, d3)| |(d1, d2, d3, d4, d5)| = 2 5 = 0.4.
Similarly, for P1 and P5, the link weight is calculated as

W P1,P5 = |(d1, d2, d3, d4) ∩ (d4, d1, d10, d40, d7, d9, d12, d52, d100)| |(d1, d2, d3, d4) ∪ (d4, d1, d10, d40, d7, d9, d12, d52, d100)| = 2 11 = 0.18.
Using the Jaccard similarity index, the final graph is built in two simple steps. In the first step, the data files that contain reviewed protein information are parsed to collect the constituent domains of each protein. If the training data contains only sequences, InterProScan [START_REF] Jones | Interproscan 5: genome-scale protein function classification[END_REF][START_REF] Quevillon | InterProScan: protein domains identifier[END_REF] is used to find the domains associated with each of the protein sequences. Then the graph is built using the domain composition of the proteins.

It is worth mentioning that the order of the domains is not maintained while computing jaccard similarity index.

Enzyme commission numbers

Enzymes are usually labelled following the Enzyme Commission (EC) system [START_REF] Cornish-Bowden | Current IUBMB recommendations on enzyme nomenclature and kinetics[END_REF], the widely used numerical enzyme classification scheme. The EC System assigns each enzyme a four digits number. This classification system has a hierarchical structure. The first level consists of the six main enzyme classes: (i) oxidoreductases, (ii) transferases, (iii) hydrolases, (iv) lyases, (v) isomerases and (vi) ligases, represented by the first digit. Each main class node further extends out several subclass nodes, specifying subclasses of the enzymes, represented by the second digit. Similarly, the third digit indicates the sub-subclass and the fourth digit denotes the sub-sub-subclasses. Let us consider as an example a Type II restriction enzyme, which is annotated as EC 3.1.21.4. The first digit, 3, denotes that it is a hydrolase. The second digit, 1, indicates that it acts on ester bonds. The third digit, 21, shows that it is an endodeoxyribonuclease producing 5phosphomonoesters. The last digit, 4, specifies that it is a Type II site-specific deoxyribonuclease.

Label propagation for protein function annotation

After building the graph from the reviewed proteins, the graph is ready to be used for the function annotation of new protein sequences. A neighborhood based label propagation algorithm is designed to perform the annotation task. Given the constituent domains of an input protein sequence, all of its neighboring proteins and their annotations are retrieved from the graph. Once the neighbors have been obtained, the weighted frequency of the labels are computed using the following formula:

f i u = v∈N(u) W u,v δ(v i , i) v∈N(u) W u,v
, where f i u is the weighted score of the candidate function i for the query protein u. And δ(v i , i) is 1 if the function v i of the protein v is same as function i, otherwise, 0. The details of the label propagation algorithm is described in Algorithm 1. Overall, for a given input sequence, the annotation algorithm works according to the flow diagram shown in Fig. 11.

Let us consider a query protein P with a set of domains D = (d5, d6, d101). Our aim is to annotate this protein with an EC Number following the label propagation algorithm, as illustrated in Fig. 12. Based on the domain similarity, protein P will have connection with proteins P2 and P4 in the running example graph. The dotted lines show the links from P to P2 and P4 in the graph along with the associated weights. Therefore, the protein P will have P2 and P4 as it's neighbors. After finding the neighbors, Algorithm 1 Label Propagation in a protein graph 1: Input: A weighted undirected protein graph ("EC annotation performance analysis" section), G = (V , E), Minimum Jaccard Similarity Index, θ, and a query protein u with domain composition d 2: Output: Weighted EC Annotations 3: procedure LABEL PROPAGATION for each i ∈ ECs do 8:

f i u = v∈N W u,v δ(v i ,i) v∈N W u,v 9:
Annotations ← Annotations ∪ {f i u } 10:

Rank the Annotations 11:

Select the top ranked annotations and assign it to the protein u return N the functional annotations of all the neighbors are propagated along with the corresponding weights. All of the functional annotations are ranked based on their cumulative weights. The top ranked function is selected as the best functional annotation for protein P. In this example, the weighted annotations for P are EC3, EC5, EC6, EC1, EC2 with cumulative weights of 0.70, 0.50, 0.50, 0.20, and 0.20, respectively. Therefore, the functional annotation for the protein P is EC3 as it has the highest weight among the propagated labels. Clearly, it is possible to select more than one high scoring functional annotations if we wish to propose more than one candidate annotation. Furthermore, node neighbours could be selected in other ways to reflect the requirements of the problem at hand.

Fig. 1 a

 1 Fig. 1 a: Distribution of EC numbers per domain composition, b: Distribution of domain compositions per EC number and c: training set statistics like a proportion of single-domain proteins, b single-EC proteins and c proteins with incomplete EC number

Fig. 10

 10 Fig. 10 Part a shows a performance comparison of GrAPFI with SVMProt (SVM, KNN, and Mixed), DEEPre, and ECPred for 2-digit EC number predictions. Part b shows the accuracy comparison of GrAPFI with DEEPre and ECPred for all 4 level of EC prediction. Part c shows the coverage of the considered methods

Fig. 11

 11 Fig.11The annotation workflow used in GrAPFI. The right-hand portion of the workflow depicts the graph construction using reviewed proteins from the UniprotKB/Swissprot. The left part shows the annotation flow

N 6 :

 6 ← FilterNeighbors(N(u), θ) ECs ← list of distinct ECs present among the neighbors N 7:

14 :

 14 function FilterNeighbors (N(u), θ))

Acknowledgements

This work was partially supported by the CNRS-INRIA/FAPs project "TempoGraphs" (PRC2243). This work is dedicated to the memory of David W. Ritchie, who recently passed away.

Funding

Not Applicable

Availability of data and materials

The data and materials can be accessed here: https://gitlab.inria.fr/bsarker/bmc_grapfi.git

Abbreviations

GrAPFI: Graph based automatic protein function inference EC: Enzyme commission GO: ontology HMM: Hidden Markov model BLAST: Basic local alignment search tool PSI-BLAST: Position specific iterative basic local alignment search tool SAAS: Statistical automatic annotation system PPI: Protein protein interaction PSSM: Position specific scoring matrix OET-KNN: Optimized evidence theoretic K-nearest neighbor SVM: Support vector machine PNN: Probabilistic neural network KNN: K-nearest neighbor FDR: Functionally discriminating residues CHIEFc: Conservation-controlled HMM Iterative procedure for enzyme family classification Authors' contributions BS and SA conceived the approach. BS developed the software and performed the analysis of the results. BS, DWR, and SA wrote the manuscript. The author(s) read and approved the final manuscript.

Authors' information

• Bishnu Sarker is a PhD student at Inria Grand-Est center and University of Lorraine in Nancy, France. • David W. Ritchie is a senior researcher at Inria Grand-Est center in Nancy, France. • Sabeur Aridhi is an associate professor at University of Lorraine, France.

Ethics approval and consent to participate Not Applicable

Consent for publication Not Applicable

Competing interests

The authors declare that they have no competing interests.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.