Baoling Guo
email: baoling.guo@g2elab.grenoble-inp.fr

Senior member, IEEE Seddik Bacha
email: seddik.bacha@g2elab.grenoble-inp.fr

Member, IEEE Mazen Alamir
email: mazen.alamir@gipsa-lab.grenoble-inp.fr

Member, IEEE Ahmad Hably
email: ahmad.hably@gipsa-lab.grenoble-inp.fr

Cédric Boudinet
email: cedric.boudinet@g2elab.grenoble-inp.fr

Generalized Integrator-Extended State Observer with Applications to Grid-connected Converters in Presence of Disturbances

Keywords: Grid-connected converters, Extended State Observer, generalized integrator, fast-varying sinusoidal disturbances, phase-locked loop, frequency-adaptive mechanism

Large-scale Fast-Varying Sinusoidal Disturbances (FVSD) impact strongly the operation of Grid-connected Converters (GcCs). A conventional Extended State Observer (ESO) is not sufficiently capable to deal with these rapid disturbances. This critical weakness has limited ESO's applications to GcCs. Therefore, a Generalized Integrator-ESO (GI-ESO) is developed from ESO to address such issues. Knowing that grid common disturbances can be expressed as a phase, frequency, or magnitude disturbed sinusoidal signal. Multi-GIs of selective resonant frequencies are hence applied into the disturbance estimation loop, which enables FVSD to be observed with a relatively low bandwidth. Besides, a frequency-adaptive mechanism is introduced in order to mitigate the effects of frequency deviations. The proposed GI-ESO resolves the trade-off between the observer's bandwidth and the noise filtering. A case of threephase phase-locked loop control system of disturbed regimes is in particular studied to assess the performance of GI-ESO and the stability of resulting controlled system. Theoretical analysis and experimental results have proven its effectiveness in applications to GcCs in presence of various disturbances.

I. INTRODUCTION

D ISTURBANCES and uncertainties present in the control systems require the design of robust controllers [START_REF] Chen | Design of unknown input observers and robust fault detection filters[END_REF]- [START_REF] Madoński | Survey on methods of increasing the efficiency of extended state disturbance observers[END_REF]. Thus, various robust approaches based on disturbance rejection have been proposed such as unknown input observer [START_REF] Chen | Design of unknown input observers and robust fault detection filters[END_REF], disturbance observer [START_REF] Chen | Disturbance-Observer-Based Control and Related Methods -An Overview[END_REF], and Extended State Observer (ESO) [START_REF] Han | From PID to active disturbance rejection control[END_REF]- [START_REF] Madoński | Survey on methods of increasing the efficiency of extended state disturbance observers[END_REF], etc. However, the performance of model-based observers and its resulting control systems are sensitive to model parameters deviations [START_REF] Chen | Disturbance-Observer-Based Control and Related Methods -An Overview[END_REF]. ESO has been widely used as a partial-model-based observer, being less dependent on the knowledge of the system [START_REF] Han | From PID to active disturbance rejection control[END_REF]- [START_REF] Guo | A review on ADRC based PMSM control designs[END_REF]. ESO is proposed in the frame of Active Disturbance Rejection Control (ADRC), which can also work independently with other controllers [START_REF] Huang | Active disturbance rejection control: Methodology and theoretical analysis[END_REF]. ESO has been originally proposed with a nonlinear structure [START_REF] Han | From PID to active disturbance rejection control[END_REF], while the linear ESO has been broadly employed as it is easier to design and to implement for practitioners [START_REF] Gao | Scaling and bandwidth-parameterization based controller tuning[END_REF].

The performances of ADRC-or ESO-based controllers primarily depend upon the accuracy of the state and disturbances estimation [START_REF] Han | From PID to active disturbance rejection control[END_REF]. Large efforts have been devoted to the tracking error convergence analysis of nonlinear or linear ESOs [START_REF] Guo | On the convergence of an extended state observer for nonlinear systems with uncertainty[END_REF]- [START_REF] Godbole | Performance Analysis of Generalized Extended State Observer in Tackling Sinusoidal Disturbances[END_REF]. The authors of [START_REF] Guo | On the convergence of an extended state observer for nonlinear systems with uncertainty[END_REF], [START_REF] Yoo | On convergence of the linear extended state observer[END_REF] have respectively proven that ESO can achieve the asymptotic convergence for constant disturbances. In another work [START_REF] Zheng | On stability analysis of active disturbance rejection control for nonlinear time-varying plants with unknown dynamics[END_REF], it has been shown that the ESO estimation error remains bounded if the changing rate of time-varying uncertainties is bounded. Theoretically, the estimation precision of ESO for time-varying disturbances can be improved by increasing its bandwidth. However, the bandwidth is limited by the presence of sensor noises and dynamic uncertainties, otherwise large high-frequency noises from the output side would be transferred to the control signal [START_REF] Sun | A DSP-based active disturbance rejection control design for a 1-kW H-bridge DC-DC power converter[END_REF]. Many attempts have been made to deal with a class of time-varying disturbances in [START_REF] Miklosovic | Discrete implementation and generalization of the extended state observer[END_REF]- [START_REF] Busawon | Disturbance attenuation using proportional integral observers[END_REF]. A Proportional Integral Generalized (GPI) observer-based control approach has been discussed in [START_REF] Cortés-Romero | Robust GPI controller for trajectory tracking for induction motors[END_REF], [START_REF] Cortés-Romero | Generalized proportional integral control for periodic signals under active disturbance rejection approach[END_REF], which is actually in line with the Generalized ESO (GESO) presented in [START_REF] Miklosovic | Discrete implementation and generalization of the extended state observer[END_REF]. A spatial observer-based repetitive controller is then proposed for mechatronic systems in [START_REF] Ramos | Spatial observerbased repetitive controller: An active disturbance rejection approach[END_REF], which is developed from the GPI approach described in [START_REF] Cortés-Romero | Robust GPI controller for trajectory tracking for induction motors[END_REF], [START_REF] Cortés-Romero | Generalized proportional integral control for periodic signals under active disturbance rejection approach[END_REF]. It has been proven that if the rth-order derivative of disturbance becomes zero, the asymptotic convergence for its state estimate can be achieved by employing rth-order GESO or GPI in [START_REF] Miklosovic | Discrete implementation and generalization of the extended state observer[END_REF]- [START_REF] Cortés-Romero | Generalized proportional integral control for periodic signals under active disturbance rejection approach[END_REF]. However, a sinusoid is a type of infinitely differentiable disturbances [START_REF] Miklosovic | Discrete implementation and generalization of the extended state observer[END_REF]. The performance of higher-order ESO in tracking sinusoidal disturbances has been analysed in [START_REF] Godbole | Performance Analysis of Generalized Extended State Observer in Tackling Sinusoidal Disturbances[END_REF]. It has been shown that a higher-order GESO offers better estimation accuracy in comparison with a conventional ESO, provided that the observer bandwidth is chosen sufficiently larger than the frequency of the disturbance. Consequently, the higher-order GESO or GPI can only ensure the precision for limitedly low-frequency sinusoidal disturbances. The GESO/GPI is not sufficiently capable to deal with very Fast-Varying Sinusoidal Disturbances (FVSD) such as harmonics-related disturbances that are multiples of grid fundamental frequency. Moreover, the authors of [START_REF] Busawon | Disturbance attenuation using proportional integral observers[END_REF] propose a proportional integral observer in order to mitigate the effects of high-frequency measurement noises or modeling errors, however, the performance regarding FVSD are not dedicatedly discussed.

Grid-connected Converters (GcCs) have various types of industrial applications such as grid-integrated renewable energy generation systems [START_REF] Andreica | Cross-flow water turbines control under grid disturbances[END_REF], micro-grid systems [START_REF] Wu | Disturbance rejection control of a fuel cell power plant in a grid-connected system[END_REF], uninterruptible power supplies [START_REF] Santos Filho | Comparison of Three Single-Phase PLL Algorithms for UPS Applications[END_REF], etc. In practice, various disturbed regimes have to be taken into account regarding the control design of GcCs. Firstly, transient voltage sags and long-term voltage unbalances commonly happen in power systems [START_REF] Bollen | Voltage Recovery after Unbalanced and Balanced Voltage Dip in Three-Phase Systems[END_REF], the unbalanced voltages then cause second-order harmonics to appear in the control systems. Besides, the power-electronics devices and industrial nonlinear loads are largely penetrated into the power grid, which induce large harmonics in the control system of GcCs [START_REF] Andreica | Cross-flow water turbines control under grid disturbances[END_REF], [START_REF] Etxeberria-Otadui | A Single Synchronous Frame Hybrid (SSFH) Multifrequency Controller for Power Active Filters[END_REF]. The voltage offset introduced by the offset of signal measurements or conversion circuits can not be neglected as well [START_REF] Chung | A phase tracking system for three phase utility interface inverters[END_REF]. Also, the power fluctuations can lead to grid frequency deviations from the nominal value (50Hz or 60Hz) [START_REF] Bevrani | Renewable energy sources and frequency regulation: survey and new perspectives[END_REF]. Note that the grid frequency deviations are limited by international standards [START_REF] Bevrani | Renewable energy sources and frequency regulation: survey and new perspectives[END_REF], [START_REF] Tsili | A review of grid code technical requirements for wind farms[END_REF]. Key features of the disturbances present in GcCs are summarized as follow: 1) They are very fast-varying sinusoidal disturbances; 2) Their frequencies are roughly known, being positive integer multiples of grid fundamental frequency.

Regarding the conventional ESO or the higher-order GESO, neither of them can efficiently identify such FVSD. Then, the performances of resulting GcCs control systems would be seriously affected. Take a three-phase Phase-Locked Loop (PLL) control system for example, it is commonly employed to achieve the grid synchronization for GcCs [START_REF] Chung | A phase tracking system for three phase utility interface inverters[END_REF], [START_REF] Golestan | Three-Phase PLLs: A Review of Recent Advances[END_REF]. The ESO-based PLL can achieve satisfying performance if no harmonics pollute the control system [START_REF] Guo | An enhanced phase-locked loop with extended state observer[END_REF]. In practice, grid voltages are not pure balanced sinusoids but disturbed by unbalances, harmonics, and voltage offsets [START_REF] Andreica | Cross-flow water turbines control under grid disturbances[END_REF], [START_REF] Bollen | Voltage Recovery after Unbalanced and Balanced Voltage Dip in Three-Phase Systems[END_REF]- [START_REF] Chung | A phase tracking system for three phase utility interface inverters[END_REF]. The harmonics components then appear in the estimated output phase angle. Finally, the control performance of GcCs is degraded due to the bad quality output of a PLL. In addition, the SRF-PLLs are commonly equipped with additional filters under disturbed regimes [START_REF] Guo | An enhanced phase-locked loop with extended state observer[END_REF], such as Decoupled PLL (DPLL) [START_REF] Rodriguez | Decoupled Double Synchronous Reference Frame PLL for Power Converters Control[END_REF] and Moving Average Filter PLL (MAF-PLL) [START_REF] Golestan | Moving Average Filter Based Phase-Locked Loops: Performance Analysis and Design Guidelines[END_REF]. The DPLL can deal with the unbalanced grid perfectly, and the results for harmonics polluted grids are practically accepted but not satisfactory [START_REF] Rodriguez | Decoupled Double Synchronous Reference Frame PLL for Power Converters Control[END_REF]. The MAF-PLL, however, presents a slow dynamic response [START_REF] Golestan | Moving Average Filter Based Phase-Locked Loops: Performance Analysis and Design Guidelines[END_REF]. More critically, neither of them can effectively deal with the voltage offset. A Generalized Integrator-Extended State Observer (GI-ESO) is thus proposed in this work, which then can be employed to address such issues for its applications into grid-connected converters.

Generalized Integrator (GI), also known as resonant controller, has been widely used for different purposes [START_REF] Andreica | Cross-flow water turbines control under grid disturbances[END_REF], [START_REF] Etxeberria-Otadui | A Single Synchronous Frame Hybrid (SSFH) Multifrequency Controller for Power Active Filters[END_REF], [START_REF] Xin | Re-Investigation of Generalized Integrator Based Filters From a First-Order-System Perspective[END_REF], [START_REF] Zmood | Stationary frame current regulation of pwm inverters with zero steady state error[END_REF]. In this work, a GI-ESO is developed from the conventional ESO. The pure integrator is maintained to track Direct Current (DC)-type or very low-frequency disturbances. Besides, GIs of different resonant frequencies are inserted into the disturbance estimation loop to deal with the FVSD. Furthermore, a frequency-adaptive mechanism is introduced in order to mitigate the effects of frequency deviations. This design enables FVSD to be observed with a relatively low bandwidth, which resolves the trade-off between the bandwidth of ESO and the high-frequency noises filtering. A threephase PLL control system under disturbed regimes is chosen as an example to assess the effectiveness of the proposed GI-ESO and the stability of resulting control system.

The paper is organized as follows: the problem formulation is described in Section II, a case of three-phase PLL system is in particular chosen as an example to formulate the control issues. The proposed GI-ESO is designed and then applied to the three-phase PLL system in Section III, meanwhile, the control stabilities and the tuning guidelines are discussed in details. Then, experimental results under various disturbed regimes are provided to assess its effectiveness in Section IV. Concluding remarks are provided in Section V.

II. PROBLEM FORMULATION

A disturbed GcC with its commonly used control structure is illustrated in Fig. 1. The control of GcCs primarily has three functions: the grid synchronization [START_REF] Guo | An enhanced phase-locked loop with extended state observer[END_REF], the grid current control [START_REF] Andreica | Cross-flow water turbines control under grid disturbances[END_REF], and the DC-link voltage control [START_REF] Guo | A Robust LESO-based DC-Link Voltage Controller for Variable Speed Hydro-Electric Plants[END_REF]. As aforementioned in Section I, large sinusoidal disturbances could be present in the control system of GcCs under disturbed regimes. These disturbances show the following common features:

1) They are very fast-varying sinusoidal disturbances;

2) Their frequencies are roughly known, being positive integer multiples of the grid fundamental frequency. the grid currents, v (abc) the grid voltages, i * d and i * q the d-axis and the q-axis grid current references, V dc the DC-link voltage, C the DC-link capacitance, L and R the filter inductance and resistance, respectively, S (a,b,c) the switching states.

The proper synchronization with grid voltage is important to ensure the whole performances of resulting GcCs control systems. In this work, the case of Synchronous Reference Frame Phase-Locked Loop (SRF-PLL) is chosen as an illustrative example to formulate the control issues and then to assess the effectiveness of the proposed GI-ESO.

A. SRF-PLL modeling and error analysis

A common design diagram of SRF-PLL is described in Fig. 2(a), which is composed of a Phase Detector (PD), a Loop Filter (LF), and a Voltage-Controlled Oscillator (VCO) [START_REF] Golestan | Three-Phase PLLs: A Review of Recent Advances[END_REF]. More design details can be found in [START_REF] Chung | A phase tracking system for three phase utility interface inverters[END_REF], [START_REF] Golestan | Three-Phase PLLs: A Review of Recent Advances[END_REF], [START_REF] Guo | An enhanced phase-locked loop with extended state observer[END_REF], but some key aspects are worth to be highlighted:

-The magnitude V m acts as a gain in the forward path of the linear model as illustrated in Fig. 2(b). This means that the magnitude variations can affect the loop gain and, therefore, the stability margins and dynamic behaviour of SRF-PLL. That is the reason why, a magnitude normalization scheme is usually introduced in order to decouple the effects from magnitude variations [START_REF] Golestan | Three-Phase PLLs: A Review of Recent Advances[END_REF].

-Two saturation blocks are also introduced: the saturation block of d-axis voltage is used to avoid the division by zero in the transient process; the saturation block of the controller is used to ensure that the estimated frequency meets the related international standards. Note that ESO-or GI-ESO-based SRF-PLLs presented later also follow the above principles. The modeling and estimation errors analysis of SRF-PLL for various disturbed regimes have been proposed in [START_REF] Chung | A phase tracking system for three phase utility interface inverters[END_REF]. The unbalanced case is in particular presented in this section. The errors analysis of SRF-PLL in presence of harmonics and voltage offsets are provided in Appendix A. The unbalanced voltages can be represented by:

 m de V v 1 s     ' qe v   g    ( ) LF s a v b v c v v  v  de m v V  qe v 1 s   g  * 0 qe v  e  
  v a v b v c   = V m   cos θ (1 + A) cos(θ -2π 3 ) (1 + B) cos(θ + 2π 3 )   (1)
where A and B represent the unbalanced levels. With Park transform [START_REF] Chung | A phase tracking system for three phase utility interface inverters[END_REF], the q-axis voltage is given by:

v qe = V m sin δ -V m A -B 2 √ 3 (cos θ cos θ -sin θ sin θ) + A + B 6 (sin θ cos θ + cos θ sin θ) (2)
The phase error δ = θ -θ is very small (θ ≈ θ), and the θ + θ ≈ 2θ, the expression can be rewritten by

v qe ≈ V m δ -E 2ω cos(2θ + φ ub ) (3)
where

E 2ω = A -B 2 √ 3 2 + A + B 6 2 φ ub = -tan -1 1 √ 3 A + B A -B .
Under an unbalanced grid condition, the sinusoidal disturbances would be transmitted to the output phase angle [START_REF] Chung | A phase tracking system for three phase utility interface inverters[END_REF]. The phase error caused by the voltage unbalances can be expressed as follows:

0 < |δ| ≤ E 2ω | cos(2θ + φ ub )| (4)

B. ESO-based SRF-PLL design principles and its limitations

An ESO-based controller is used in order to enhance control dynamics of SRF-PLLs in [START_REF] Guo | An enhanced phase-locked loop with extended state observer[END_REF]. First, the design principles are briefly reviewed. Then, the ESO's limitations when used in disturbed grids are analyzed and experimentally revealed.

The control diagram of the ESO-based SRF-PLL is shown in Fig. 3, where an ESO-based corrector is used as the LF. ESO is used to track the unknown disturbances in the PLL system. The estimated disturbances x2 are actively compensated into the signal u 0 in real-time.

   p k 0 u 1 s    ' qe v 0 1 b     g    * 0 qe v      1 s Conventional linear ESO 1 s  2 x 1 x 1  2  1 de v m V qe v  Fig. 3: Conventional ESO-based SRF-PLL control diagram, x1
is the estimation of v qe , x2 is the estimation of total unknown disturbances, β 1 and β 2 are the tuned positive gains of ESO, b 0 is the estimation of control gain

The first step before using an ESO-based controller is to reformulate a practical control issue in the canonical form of cascaded integrators with generalized disturbances [START_REF] Guo | A review on ADRC based PMSM control designs[END_REF]. By differentiating and simply manipulating the expression (3), the canonical form becomes:

1 V de dv qe dt = V m V de d( θ -θ) dt = b 0 ∆ω + f (5)
An appropriate disturbance definition is critical to achieve a high-efficiency ESO [START_REF] Guo | A review on ADRC based PMSM control designs[END_REF]. In the PLL, v de stands for the realtime estimation of V m , therefore, the control gain b = V m /v de approximates 1. It is reasonable to set the estimation of control gain with b 0 = 1. Moreover, the nominal grid frequency ω g is feed-forwarded in the control input. Such compensation improves ESO estimation's precision and response dynamics. Finally, the total disturbances f is formulated by: To highlight, the frequency deviations, the voltage magnitude variations, and phase jumps have a physically limited changing rate in the power system like the majority of related signals. For instance, the frequency jump is related to the speed variations of the synchronous generators which have their own inertia. For the inrush currents and voltages, the phase jump can be considered as ramps even for phenomenon with a very short interval time. In conclusion, the signals we are referring to are differentiable with bounded derivatives.

f = [(b -b 0 )∆ω + (ω g -ω) + d x ] + 2ωE 2ω sin(2θ + φ ub ) = d (0) + d (2ω) (6) where d (0) = (b -b 0 )∆ω + (ω g -ω) + d x ,
In the sequel, the following definitions are used: ∆ω is the control input u, v * qe is the reference r, v qe is the output y, which is defined as a state variable by x 1 = v qe , f is treated as an extended state variable, and expressed by x 2 = f , its time derivative is given by ẋ2 = h [START_REF] Guo | An enhanced phase-locked loop with extended state observer[END_REF]. Now, the state space model can be constructed as follows:

ẋ1 ẋ2 = 0 1 0 0 x 1 x 2 + b 0 0 u + 0 1 h (7)
Based on [START_REF] Guo | LADRC applied to variable speed micro-hydro plants: Experimental validation[END_REF], the ESO can be formulated as follows:

ẋ1 ẋ2 = 0 1 0 0 x1 x2 + b 0 0 u + β 1 β 2 [x 1 -x1 ] ( 8
)
where x1 is the estimation of v qe , x2 is the estimation of total disturbances, β 1 and β 2 are positive gains of ESO. By subtracting ( 8) from ( 7), the error dynamics are defined by:

ė1 ė2 = -β 1 1 -β 2 0 Ae e 1 e 2 + 0 1 h (9)
where e 1 = x 1 -x1 is the state estimation error, and e 2 = x 2 -x2 is the disturbances estimation error. The positive gains β 1 = 2ω o and β 2 = ω 2 o are chosen such that the characteristic polynomial ( 10) is Hurwitz stable [START_REF] Gao | Scaling and bandwidth-parameterization based controller tuning[END_REF].

λ(s) = (sI -A e ) -1 = s 2 + β 1 s + β 1 2 = (s + ω o ) 2 (10)
where ω o is referred to as the bandwidth of ESO [START_REF] Gao | Scaling and bandwidth-parameterization based controller tuning[END_REF].

Based on Fig. 3, the control law is formulated by

∆ω = k p (v * qe -v qe ) -x2 b 0 (11)
If the total disturbances were correctly estimated and compensated, the original plant model would be reduced to be a pure integrator [START_REF] Tatsumi | On the enhanced ADRC design with a low observer bandwidth[END_REF]. Then, the closed-loop transfer function can be approximated by

H c (s) ≈ k p 1 s 1 + k p 1 s = 1 1 + s kp (12)
where ω c ≈ k p is referred to as the closed-loop bandwidth.

More details on convergence analysis, stability analysis, and tuning methods can be found in [START_REF] Gao | Scaling and bandwidth-parameterization based controller tuning[END_REF], [START_REF] Guo | LADRC applied to variable speed micro-hydro plants: Experimental validation[END_REF], [START_REF] Tatsumi | On the enhanced ADRC design with a low observer bandwidth[END_REF], [START_REF] Tian | Frequency response analysis of active disturbance rejection based control system[END_REF].

C. ESO's limitations in tracking the FVSD

Based on (9), the transfer function of disturbances estimation error becomes:

f e2 (s) = e 2 (s) f (s) = s(s + 2ω o ) (s + ω o ) 2 (13)
The bode plots of f e2 (s) are shown in Fig. 4, with ω o ∈ {300, 400, 500} (rad/s). It can be observed that ESO behaves like a low-pass filer with respect to the input disturbances [START_REF] Tatsumi | On the enhanced ADRC design with a low observer bandwidth[END_REF]. The conventional ESO can estimate the unknown disturbance very well only when the disturbances are composed of lowfrequency components. The estimation errors of disturbances can be reduced by increasing the bandwidth ω 0 , however, the bandwidth would be limited by the practical noises. Following the tuning rules as [START_REF] Guo | An enhanced phase-locked loop with extended state observer[END_REF], the control parameters are configured with ω c = 100 rad/s and ω o = 400 rad/s. The results of ESO-based SRF-PLL are shown in Fig. 6. The grid phase angle is correctly estimated for the ideal grid; while the phase estimation errors are not negligible in presence of grid unbalances. The reason is that the total disturbances defined by ( 6) are partially filtered by ESO, the second-order harmonics are then introduced into the estimated phase angle. The ESO-based PLL can also be equipped with additional filters under disturbed regimes such as DPLL and MAF-PLL, but neither of them can efficiently deal with the voltage offsets. A GI-ESO is therefore proposed to address such issues for its applications into grid-connected converters.

III. GI-ESO DEIGN AND ITS APPLICATION TO SRF-PLL

In this section, the main properties of GI-ESO are firstly highlighted. Then, a GI-ESO-based SRF-PLL is designed in order to assess the performance of the proposed GI-ESO and the stability of its resulting control system.

A. GI-ESO design

The GI-ESO is developed from the previous linear ESO. It inherits the ESO's original advantages: high robustness, good disturbance rejection abilities within the bandwidth of ESO, and simple tuning methods [START_REF] Gao | Scaling and bandwidth-parameterization based controller tuning[END_REF]. More critically, the capability in dealing with FVSD are enhanced for the proposed GI-ESO.

1) Extended structure and mathematical model

The proposed GI-ESO and the conventional ESO are comparatively presented in Fig. 7. The pure integrator in the disturbance loop is kept to track DC-type or very low-frequency disturbances, owing to the infinite gain of the pure integrator at zero-frequency [START_REF] Xin | Re-Investigation of Generalized Integrator Based Filters From a First-Order-System Perspective[END_REF]. Moreover, the GIs of selective resonant frequencies are inserted in parallel in order to observe the FVSD of different frequencies.

  2 2 2 2 r r k s s   2(0) x 1 2( ) ˆr x  2 x ' 2  ' 1  y 1 x 2 r  0 b u  2 2( ) ˆr x  rad/s  2 2 ri ri k s s   2( ) ˆri x  ri   rad/s 1 2 2 1 r r k s s   1 r  Fig. 7: Proposed GI-ESO design diagram
The inserted GI can be generally expressed by [START_REF] Miklosovic | Discrete implementation and generalization of the extended state observer[END_REF], wherein the infinite gain of the integrator is shifted to ω ri [START_REF] Xin | Re-Investigation of Generalized Integrator Based Filters From a First-Order-System Perspective[END_REF].

R i (s) = k ri s s 2 + ω 2 ri , i = 1, 2, . . . ( 14
)
where k ri is the proportional gain, and ω ri is the selective resonant frequency of each GI.

A GI-ESO is extended from a conventional ESO (8) by

ẋ1 ẋ2 = 0 1 0 0 x1 x2 + b 0 0 u + β 1 β 2 (1 + R(s)s) [x 1 -x1 ] (15) where x2 = x2(0) + x2(ωr1) + x2(ωr2) +• • •+ x2(ωri) is the total estimated disturbances, R(s) is a general expression that is R(s) = 0 for ESO; R(s) = R 1 (s)+R 2 (s)+• • •+R i (s) for GI- ESO.
The observer gains are expressed by β 1 = ξω 0 , β 2 = ω 2 0 , and ξ is the tuning coefficient. Detailed theoretical analysis and selection guidelines are provided later in the sequel.

2) Estimation error assessment

Similarly, by subtracting ( 15) from ( 7), the estimation errors can be expressed by

ė1 ė2 = -β 1 1 -β 2 (1 + R(s)s) 0 e 1 e 2 + 0 1 h (16)
By replacing h with ḟ (s), the transfer function between e 2 (s) and f (s) becomes:

F e2 (s) = s(s + ξω o ) s 2 + ξω o s + ω 2 o (1 + R(s)s) (17)
The bode plots of F e2 (s) are shown in Fig. 8, in which the disturbance estimation errors of ESO and GI-ESO are comparatively presented. For these two observers, the bandwidth is set to ω o ∈ {300, 400, 500} (rad/s), and ξ = 5. The inserted GIs are respectively set to k r1 , = π, ω r1 = ω g for voltage offset, k r2 = 5π, ω r2 = 2ω g for unbalanced grids, k r3 = 10π, ω r3 = 6ω g for 5th-order and 7th-order harmonics (proven in Appendix A), and ω g = 100π rad/s. The bode plots indicate that GI-ESO and ESO performs closely at either the very low-frequency or the high-frequency range. Moreover, being different from ESO, the proposed GI-ESO can also correctly track very FVSD. The reason is that the inserted GIs enable rather high gains at specific resonant frequencies.

B. GI-ESO-based SRF-PLL design

The primary principles described in Section II-B are followed here, moreover, several particular points are highlighted.

1) Disturbances compensation

It has been proven in Section III-A that the proposed GI-ESO is capable of estimating very FVSDs. But another key question is how to correctly compensate the estimated disturbances. The GI-ESO-based SRF-PLL control diagram is shown in Fig. 9. The phase angle is indirectly estimated by controlling the q-axis voltage. If the total disturbances were compensated by the control signal ∆ω, the secondorder frequency components would remain in the estimated frequency ω and the output phase θ. Since the information of phase and frequency are more important in a PLL system. We can force the harmonics components to stay in the q-axis voltage (see Fig. 9). It comes from (3) that, the integration of x2(2ω) actually represents the v qe fluctuations caused by unbalanced voltages. Finally, the sinusoidal errors can be actively removed from the output estimated phase.

   ' p k 0 u GI-ESO 1 s     ' qe v 0 1 b   * 0 qe v     1 s 2(2 ) x  2(0) x g    2 2( ) r g       b
Fig. 9: GI-ESO-based SRF-PLL under voltage unbalances

2) Frequency-adaptive mechanism

Note that the actual grid frequency possibly derives from its nominal value ω g [START_REF] Bevrani | Renewable energy sources and frequency regulation: survey and new perspectives[END_REF]. If the resonant frequency of GI is set to be constant, then the performances of GI-ESO when used in GcCs can be deteriorated. A frequency-adaptive mechanism is hence considered in this design (see Fig. 9). For example, under voltage unbalances, the value of the resonant frequency of GI is adaptively adjusted by considering the estimated frequency deviation with ω r2 = 2(ω g + ∆ω). As the applied resonant frequencies of GIs are corrected in realtime, the negative effects induced by frequency deviations can be significantly mitigated.

3) Control law formulation

The DC-type or low-frequency disturbances x2(0) are directly compensated in the signal u 0 . The integration of x2(2ω) is added to the q-axis voltage reference. Finally, the modified control law is expressed by

∆ω = k p (v * qe + x2(2ω) dt -v qe ) -x2(0) b 0 (18)
The GI-ESO handles all sinusoidal components that are present in the SRF-PLL as generalized disturbances. The structure of GI-ESO can be flexibly adapted to different disturbed regimes by adjusting the selective frequencies.

4) Stability analysis

The stability index is a crucial factor for evaluating the quality of a controller. Stabilities analysis methods of ESO can be referred for that of GI-ESO [START_REF] Tian | Frequency response analysis of active disturbance rejection based control system[END_REF]- [START_REF] Jin | Linear active disturbance rejection control tuning approach guarantees stability margin[END_REF]. The frequency domain properties are therefore quantitatively analyzed in this section. This is done under the following assumptions:

-A mismatch between b and b 0 is assumed, the control gain is set to be b = 1.2, and b 0 = 1. -The parameters k r1 and ω r1 are defined as the proportional gain and the resonant frequency for the GI dealing with voltage offsets, similarly, ω r2 and k r2 for voltage unbalances, ω r3 and k r3 for harmonics. -Phase margins around 45 • are supposed to be acceptable. From expressions ( 15) and ( 18), the transfer function can be obtained:

u(s) y(s) = - 1 b 0 ω c s 2 + ω o 2 s + ξω o ω c s + ω o 2 ω c s(s + ξω o ) + ω o 2 (s + ω c )R(s) (19)
u(s) r(s) = 1 b 0 ω c s 2 + ω o 2 ω c R(s)s + ξω o ω c s + ω o 2 ω c s(s + ξω o ) + ω o 2 (s + ω c )R(s) (20)
where r(s), u(s), and y(s) are the Laplace transforms of r = v * qe , u = ∆ω, and y = v qe , respectively.

By introducing the following notations:

C(s) = 1 b 0 ω c s 2 + ω o 2 s + ξω o ω c s + ω o 2 ω c s(s + ξω o ) + ω o 2 (s + ω c )R(s) (21)
H(s) = ω c s 2 + ω o 2 ω c R(s)s + ξω o ω c s + ω o 2 ω c ω c s 2 + ω o 2 s + ξω o ω c s + ω o 2 ω c ( 22
)
the ESO-based controller can be formulated to a two-degreeof-freedom (2dof) closed-loop system as shown in Fig. 10. Now, the feedback control loop has a standard structure of the feedback and the controlled plant [START_REF] Jin | Linear active disturbance rejection control tuning approach guarantees stability margin[END_REF]. Firstly, the phase margins of the loop transfer function of the controlled system are discussed. The loop transfer function is given by:

L(s) = C(s)G p (s) = 1 b 0 ω c s 2 + ω o 2 s + ξω o ω c s + ω o 2 ω c s(s + ξω o ) + ω o 2 (s + ω c )R(s) b s (23)
where the plant model is G p (s) = b/s, and b = V m /v de . The unbalanced voltage case is exclusively studied in order to analyze the effects on stability with respect to each parameter (ξ, ω o , ω c , k r ). The bode plots of L(s) for different coefficients ξ ∈ {2, 3, 4} are shown in Fig. 11. The remaining parameters are set to ω o = 400 rad/s, ω c = 100 rad/s, k r2 = 5π, ω r2 = 2ω g . It can be observed that the gain margin is negative infinitive, however, the phase margin cannot sufficiently meet stabilities requirements for ξ = 2 or ξ = 3. While too large ξ enables slow dynamics response [START_REF] Gao | Scaling and bandwidth-parameterization based controller tuning[END_REF]. Therefore, β 1 = 4ω 0 and β 2 = ω 2 0 are chosen in this work as well as in the forthcoming analysis. The bode plots of L(s) for different ω o and ω c are shown in Figs. 12 and13, respectively. The GI is applied with k r2 = 5π, and ω r2 = 2ω g . The ESO's tuning rule can be referred, that is commonly configured by ω o = 3 ∼ 5ω c [START_REF] Gao | Scaling and bandwidth-parameterization based controller tuning[END_REF]. Then, we make the following settings: ω o ∈ {300, 400, 500} (rad/s), ω c = 100 rad/s for Fig. 12; ω o = 400 (rad/s), ω c ∈ {80, 100, 120} (rad/s) for Fig. 13. It can be seen that the phase margin will decrease with larger observer bandwidth ω o , which can deteriorate stabilities of the closed-loop control system. However, the phase margins are insensitive to the changes of ω c . The bode plots of L(s) of different k r2 are shown in Fig. 14. The GIs are set to k r2 ∈ {π, 5π, 10π}, and ω r2 = 2ω g . Larger k r2 enables higher gains of GI at the selective frequencies which helps improving the estimation accuracy of GI-ESO. However, the phase margins are sensitive to the changes of k r2 as shown in Fig. 14. The high gains tend to bring destabilizing effects on the resulting closed-loop system. Too large gains would be risky in respect to stability issues. Consequently, the gains of inserted GI are limited accordingly. Now, the stabilities of GI-ESOs with different selective frequencies are analyzed. The resonant frequencies of applied GIs are respectively configured with ω r1 = ω g for voltage offset, ω r2 = 2ω g for grid unbalance, and ω r3 = 6ω g for 5th-order and 7th-order harmonics (see Appendix A).

The bode plots of L(s) for different selective frequencies are provided in Fig. 15. The same gains of GIs with k r1 = k r2 = k r3 = 5π are used. It can be seen that the phase margin decreases critically for system that applies a GI with lower resonant frequency. Therefore, the gains of GIs must be selected by taking into account its applied resonant frequency. Fig. 15: Bode plots of L(s) for various selective frequencies Fig. 16 shows the bode plots of L(s) that apply different gains k ri , the remaining parameters keep the same configurations. Compared to the ESO, although the stability margins decrease due to the high gains of GIs, the stability condition can be ensured by properly adjusting k ri . The bode plots of L(s) for multi-GIs in parallel are provided in Fig. 17, with the same parameters settings as that applied in Fig. 16. It can be observed that multi-GIs cause a further decrease of phase margin (ξ = 4, P m = 38.3 • ). The coefficient ξ is thus adjusted by using ξ = 5 for multi-GIs to ensure that the stability margin is sufficient. The stability of the pre-filter H(s) [START_REF] Bollen | Voltage Recovery after Unbalanced and Balanced Voltage Dip in Three-Phase Systems[END_REF] is guaranteed, as all coefficients in its denominator are positive [START_REF] Tian | Frequency response analysis of active disturbance rejection based control system[END_REF]. Therefore, the closed-loop stability of the loop that linked the synthesized reference signal r * = x2(2ω) dt + v * qe to the output can be determined by

G cl (s) = C(s)G p (s) 1 + C(s)G p (s) (24)
The bode plots of G cl (s) are shown in Fig. 19, the same parameters are configured as that of Fig. This part provides a general stability analysis methodology, which can be adapted to future implementations of GI-ESOs.

6) Parameter selection guidelines

Firstly, the structure and the selective frequencies of GI-ESOs need to be adapted to the sinusoidal disturbances present in a real applications. Then, proper parameters have to be selected by taking into account several factors: system stability, estimation's precision and response dynamics. The parameter tuning guidelines are summarized as follow:

-Firstly, regarding parameters ω c and ω 0 , similar tuning rules discussed in [START_REF] Gao | Scaling and bandwidth-parameterization based controller tuning[END_REF], [START_REF] Guo | LADRC applied to variable speed micro-hydro plants: Experimental validation[END_REF] can be followed. -Then, ω ri is chosen by adapting to frequencies of the present sinusoidal disturbances, and k ri is tuned to ensure both the estimation precision and the stability. -Moreover, the coefficient ξ have to be corrected for GI-ESO to ensure the stability condition. -Finally, the applied control parameters must be retuned and validated in the experimental tests. Finally, the control parameters are configured as given in Table I, which will be validated in the forthcoming tests.

IV. EXPERIMENTAL RESULTS AND ANALYSIS

The GI-ESO-based SRF-PLL is tested in order to check the effectiveness of the proposed GI-ESO in dealing with FVSD. Various disturbed schemes (unbalances, harmonics, and voltage offset) have been considered. Experimental results are obtained by using the dSPACE controller (DS2005).

First of all, the grid unbalances are set to be the same as the settings in Fig. 5. The performances of GI-ESO-based PLL are shown in Fig. 20. By introducing the GI, the second-order grid frequency disturbance is approximately estimated by the GI-ESO, namely x2(2ω) ≈ d (2ω) . Then, the fluctuations are actively transmitted to the q-axis voltage. Fig. 21 shows that the integration of x2(2ω) enables a tight tracking of the desired theoretical value ∆v qe(2ω) . The sinusoidal components are forced to remain in the q-axis voltage. Finally, the errors are removed from the estimating phase angle. The performances are significantly improved compared with the ESO-based PLL (see Fig. 6 for comparison). x  (2 ) The sinusoidal disturbances induced by voltage offset have the same frequency component as that of grid voltage [START_REF] Chung | A phase tracking system for three phase utility interface inverters[END_REF]. As shown in Fig. 23, the results indicate that the sinusoidal disturbances are tightly tracked by GI-ESO. The PLL enables the performance precision regarding voltage offsets.

d  10    o (a) Phase jump +10 • 2(2 ) x  2(0) x  (2 )
Regarding harmonics regime, the 5th-order harmonic of magnitude 10V and the 7th-order harmonic of magnitude 5V are injected simultaneously into three-phase voltages. The 6thorder harmonic disturbances appear in the SRF-PLL control system (proven in Appendix A). The control parameters are given in Table I. Fig. 24 shows that the harmonic disturbances are approximately estimated by GI-ESO. The high-order harmonic errors in the estimated phase become negligible. For a real application, the PLLs are expected to be robust against unbalances, offsets, and harmonics. Fig. 25 provides the results of GI-ESO-based PLL that deal with various disturbances simultaneously. Multi-GIs of different resonant frequencies are incorporated in parallel, being configured in Table I. The total multi-disturbances can be identified by the GI-ESO. The overshoot of estimated phase angle is less than 1.1%, and the estimated phase quickly follow the phase of grid voltage. The proposed GI-ESO-based PLL achieves practically satisfied performances in presence of various disturbances. In this paper, the presented results are focused on PLL issues. However, the proposed GI-ESO can be used in other applications that might be encountered in GcCs like the DClink control or the grid current control contexts. These are topics of undergoing works.

V. CONCLUSIONS This work is a development/extension of a previously proposed work on linear ESO, therefore, it inherits the ESO's original advantages. Moreover, GIs of different resonant frequencies are applied into the disturbance estimation loop, which enables FVSD to be observed with a relatively low bandwidth. Besides, a frequency-adaptive mechanism is introduced in order to mitigate the effects of frequency deviations. Furthermore, the structure of GI-ESO can be flexibly adjusted regarding different disturbed regimes present on the field.

A case of three-phase SRF-PLL control system is chosen as an example to assess GI-ESO's performance and the stability of its resulting control system. Firstly, the frequency domain properties, particularly the phase margins are quantitatively analyzed with respect to different control parameters. The parameters tuning guidelines are highlighted as well. Moreover, the experimental results indicate that various FVSD can be efficiently observed by introducing GI-ESOs. The sinusoidal components are removed from the estimated phase angle by compensating correctly the estimated disturbances. Besides, in the case of frequency deviations, the frequency-adaptive mechanism mitigates the phase estimation errors. The GI-ESO also achieves high robustness under various internal or external disturbances. Both theoretical analysis and experimental results have proven the effectiveness of the proposed solution when used in GcCs in presence of various disturbances.

The proposed GI-ESO solution breaks the limits of ESObased controllers or ADRC for their applications in various disturbed control schemes of GcCs (grid synchronization, current control, and DC-link voltage control). Moreover, it has a wide applications scope for many systems are impacted by the presence of very FVSD. If the frequencies of sinusoidal disturbances are not known, a simple PLL can be incorporated to obtain the frequency information.

Fig. 1 :

 1 Fig.1: Schematic diagram of a disturbed grid-connected converter with the commonly used control structure, with i (abc) the grid currents, v (abc) the grid voltages, i * d and i * q the d-axis and the q-axis grid current references, V dc the DC-link voltage, C the DC-link capacitance, L and R the filter inductance and resistance, respectively, S (a,b,c) the switching states.

Fig. 2 :

 2 Fig. 2: (a) Schematic diagram of conventional SRF-PLL, (b) linear model, with v qe , v de the voltages in a dq synchronous reference frame, v α , v β the voltages in a αβ rotating reference frame, V m the grid voltage magnitude, ω g the grid nominal angular frequency, ω the estimated angular frequency, ∆ω the estimated angular frequency deviation, θ the estimated phase angle, and LF (s) is the transfer function of LF

  and (b -b 0 )∆ω is due to the small mismatch between b and b 0 , ω g -ω is the grid frequency variations, d x represents extra disturbances such as low-frequency disturbances induced by phase jumps, time-varying parameters, nonlinear dynamics, etc.; d (2ω) = 2ωE 2ω sin(2θ + φ ub ) represents the harmonics disturbances.

Fig. 4 :Fig. 5 :

 45 Fig. 4: Bode diagram of the transfer function f e2 (s)

Fig. 6 :

 6 Fig. 6: Performance of conventional ESO-based SRF-PLL, f the total disturbances (50/div), x2 the estimated disturbances by ESO (50/div), δ the estimated phase error (2 • /div)

Fig. 8 :

 8 Fig. 8: Bode diagram of the transfer function F e2 (s)

Fig. 10 :

 10 Fig. 10: 2dof closed-loop diagram of GI-ESO-based SRF-PLL control system, d w represents external disturbances

Fig. 11 :

 11 Fig. 11: Bode plots of L(s) for coefficient ξ ∈ {2, 3, 4}

Fig. 12 :Fig. 13 :

 1213 Fig. 12: Bode plots of L(s) for ω o ∈ {300, 400, 500} (rad/s)

Fig. 14 :

 14 Fig. 14: Bode plots of L(s) for k r2 ∈ {π, 5π, 10π}

16 :

 16 Bode plots of L(s) different k ri

Fig. 17 :

 17 Fig. 17: Bode plots of L(s) for multiple selective frequencies

17

 17

Fig. 18 :Fig. 19 :

 1819 Fig. 18: Bode plots of L(s) for b ∈ {0.5, 1, 1.5}

  Common parameters ωc = 100 rad/s, ωo = 400 rad/s, b0 = 1 Coefficient ξ ξ = 4 for single GI, ξ = 5 for multi-GIs Offsets kr1 = π, ωr1 = (ωg + ∆ω) Unbalances kr2 = 5π, ωr2 = 2(ωg + ∆ω) Harmonics kr3 = 10π, ωr3 = 6(ωg + ∆ω)

Fig. 20 :Fig. 21 :

 2021 Fig. 20: Performance of GI-ESO-based PLL, d (2ω) the sinusoidal disturbances (50/div), x2(2ω) the estimated sinusoidal disturbances (50/div), x2(0) the estimated DC-type or lowfrequency disturbances (50/div), δ the phase error (2 • /div)

Fig. 22 :

 22 Fig. 22: Robustness test of a GI-ESO-based PLL, d (2ω) the sinusoidal disturbances (50/div), x2(2ω) the estimated sinusoidal disturbance (50/div), x2(0) the estimated low-frequency disturbance (50/div), δ the phase error (5 • /div)

  .

Fig. 23 :

 23 Fig. 23: Performance of GI-ESO-based PLL in presence of voltage offset, d (ω) the sinusoidal disturbances (50/div), x2(ω) the estimated sinusoidal disturbances (50/div), x2(0) the estimated DC or low-frequency disturbances (50/div), δ the phase error (2 • /div)

Fig. 24 :

 24 Fig. 24: Performance of GI-ESO-based PLL in presence of harmonics disturbances, d (6ω) the 6th-order harmonic disturbances (50/div), x2(6ω) the estimated 6th-order harmonic disturbances (50/div), x2(0) the estimated DC or low-frequency disturbances (50/div), δ the phase error (2 • /div)

Fig. 25 :Fig. 26 :

 2526 Fig. 25: Performance of GI-ESO-based PLL in presence of multi-harmonic disturbances, d 2(ω,2ω,6ω) the multi-harmonic disturbances (100/div), x2(ω,2ω,6ω) the estimated disturbances (100/div), x2(0) the estimated DC or low-frequency disturbances (50/div), δ the phase error (2 • /div)

TABLE I :

 I Configurations of control parameters

This work is supported by PSPC Innov'hydro project, which brings together General Electric (GE) Renewable Energy, eDF ( Électricité de France), Grenoble INP, and other key players in the hydroelectric sector, in France.

B.

ERROR ANALYSIS AND DISTURBANCES DEFINITION

The error analysis of voltage harmonics and offset are recalled in order to define the disturbances present in the control system of GI-ESO-based SRF-PLL.

A.1. Voltage harmonics

The grid voltages with harmonics can be represented by:

where V 1 , V 5 , V 7 , . . . represent the harmonics magnitudes.

The q-axis voltage can be derived and linearised by:

The phase error caused by voltage harmonics is given by

where

The voltage harmonics enables the phase errors with the frequency components of 6ω, 12ω, . . . , . The total disturbance of ESO-based PLL control system can be represented by:

A.2. Voltage offset

The measurement offset or the signal conversion circuits can sometimes cause voltage offset [START_REF] Chung | A phase tracking system for three phase utility interface inverters[END_REF]. The unbalanced grid voltages can be represented by:

where V ao , V bo , and V co represent the offset voltage levels.

The voltages in the αβ frame are as follows:

where

The q-axis voltage can be derived and linearised by:

where

The phase error caused by the offset has the same frequency component as that of the grid voltage [START_REF] Chung | A phase tracking system for three phase utility interface inverters[END_REF]. The phase error caused by the voltage offset is given by:

The total disturbances present in ESO-based PLL control system can be expressed by:


Sebastien Theil
email: sebastien.theil@inrae.fr

Etienne Rifa
email: etienne.rifa@inrae.fr

Florent Murat

Simon Roux

rANOMALY: AmplicoN wOrkflow for Microbial community AnaLYsis [version 1; peer review: 2 approved]

Keywords: 16S, ITS, amplicon sequencing, metagenomics, microbial community, R package Conceptualization, Methodology, Software, Validation, Writing -Original Draft Preparation, Writing - Conceptualization, Methodology, Software, Validation, Writing -Original

Bioinformatic tools for marker gene sequencing data analysis are continuously and rapidly evolving, thus integrating most recent techniques and tools is challenging. We present an R package for data analysis of 16S and ITS amplicons based sequencing. This workflow is based on several R functions and performs automatic treatments from fastq sequence files to diversity and differential analysis with statistical validation. The main purpose of this package is to automate bioinformatic analysis, ensure reproducibility between projects, and to be flexible enough to quickly integrate new bioinformatic tools or statistical methods. rANOMALY is an easy to install and customizable R package, that uses amplicon sequence variants (ASV) level for microbial community characterization. It integrates all assets of the latest bioinformatics methods, such as better sequence tracking, decontamination from control samples, use of multiple reference databases for taxonomic annotation, all main ecological analysis for which we propose advanced statistical tests, and a cross-validated differential analysis by four different methods. Our package produces ready to publish figures, and all of its outputs are made to be integrated in Rmarkdown code to produce automated reports.

Introduction

Studies of microbial communities tends to become a daily routine analysis for lots of laboratories and the main method to explore microbial diversity is metabarcoding, which is an amplicon targeted sequencing method (16S for bacteria and ITS for fungi). Metabarcoding generates a large amount of data and a lot of applications already exist for their processing (FROGS 1 , qiime [START_REF] Bolyen | Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2[END_REF] ). Methods and software are continuously evolving and the main challenge for bioinformaticians is to implement the most recent and effective ones in their analysis. Here we present rANOMALY, a scalable and lightweight R package which is able to handle every step of a metabarcoding analysis, from read cleaning, contaminant filtering, taxonomic assignment, to advanced statistical analysis. rANOMALY is fully implemented in R language in which each step correspond to one function, allowing to easy implementation of new features or tools while being easy to use and maintain. The package allows the workflow to be executed on any R environment. rANOMALY only needs a CSV table describing the metadata for each sample, and a folder containing the corresponding fastq files as input. It can produce high quality figures for Rmarkdown reports along with statistical tests ready for publication.

The workflow is illustrated in Figure 1.

Implementation

rANOMALY is an R package depending on other CRAN, Bioconductor, and git R packages. It is easy to install via devtools::install_git function.

Raw sequence processing

Samples must have been previously demultiplexed into one file per sample with the file name following this syntax: {sampleid}_R [12].fastq. The denoising process is handled using the dada2 R package [START_REF] Callahan | DADA2: Highresolution sample inference from Illumina amplicon data[END_REF] which produces amplicon sequence variants (ASV) as a taxonomic unit. This improves resolution of the potential presence of microbial organisms by using a prediction model to correct sequencing errors before aggregating similar sequences. rANOMALY handles processing any region of 16S (V1 to V9) and ITS amplicon (ITS1, ITS2) sequences, in which an additional step of cutadapt [START_REF] Martin | Cut adapt removes adapter sequences from highthroughput sequencing reads[END_REF] removes ITS probes left in some short sequences. For 16S amplicon, primers are trimmed based on the primer sequence length. ASV identifiers are sequences translated into MD5 hashes which are unique identifiers based on the DNA sequences offering the possibility to be compared between projects. This step results in an object with representing sequences, a raw ASV counts table and a text file containing statistics from the denoising process.

Taxonomic assignment

Taxonomic assignments of ASVs are carried out by IDTAXA (part of the DECIPHER package), an algorithm based on a machine learning method [START_REF] Murali | IDTAXA: a novel approach for accurate taxonomic classification of microbiome sequences[END_REF] . We implemented a functionality to compute assignment with two reference databases. For example, 16S amplicon sequences can be assigned with an environment specific database (DAIRYdb 6 , HITdb [START_REF] Ritari | Improved taxonomic assignment of human intestinal 16S rRNA sequences by a dedicated reference database[END_REF] , MIDAS [START_REF] Nierychlo | MiDAS 3 : An ecosystem-specific reference database, taxonomy and knowledge platform for activated sludge and anaerobic digesters reveals species-level microbiome composition of activated sludge[END_REF] ) and a general database like SILVA [START_REF] Quast | The SILVA ribosomal RNA gene database project: improved data processing and web-based tools[END_REF] or GreenGenes [START_REF] Desantis | Greengenes, a Chimera-Checked 16S rRNA Gene Database and Workbench Compatible with ARB[END_REF] . The assignment with the best confidence and the lowest rank is kept, thus increasing assignment depth and accuracy. As taxonomy can differ between reference databases, we implemented a taxonomy validation step in rANOMALY to unravel taxonomy inconsistency like taxa with multiple ancestors or empty ranks. As supplementary features, rANOMALY functions allow users to create their own IDTAXA formatted specific reference database. The first step consists of filling the taxonomic table empty fields with the last known rank, and checking for taxonomy incongruencies as in the assignment function. A taxid file as used with RDP classifier [START_REF] Wang | Naïve Bayesian Classifier for Rapid Assignment of rRNA Sequences into the New Bacterial Taxonomy[END_REF] is constructed and then a last function takes as input, the corrected taxonomy table, the fasta file and the taxid file to generate the IDTAXA formatted database.

Phylogenetic tree

Phylogenetic tree is generated in three steps. First, sequences are aligned with AlignSeqs function from DECI-PHER package [START_REF] Wright | DECIPHER: harnessing local sequence context to improve protein multiple sequence alignment[END_REF] by the guided-tree method. Then, distance matrix is calculated with the dist.ml function from the phangorn package. And finally, neighbour joining and pml function computes the likelihood of the phylogenetic tree.

The abundance, metadata, taxonomic table, reference sequences and phylogenetic tree are merged into a phyloseq object [START_REF] Mcmurdie | phyloseq: An R Package for Reproducible Interactive Analysis and Graph-ics of Microbiome Census Data[END_REF] .

Decontamination

ASVs have the advantage of enabling the distinction between contaminants and the real community. We have integrated the decontam package 14 into the workflow. Indeed using OTU based clustering methods can agglomerate contaminant sequences with real sample sequences, the whole cluster could hence be considered as contaminant by mistake. Working with ASVs allows the use of R package decontam which will sensitively exclude contaminant ASVs. It integrates two main methods, one based on the prevalence of the contaminant in the control samples, and another one based on the DNA concentration of the samples. Moreover, our decontamination step allows users to apply various filters such as low ASVs frequency, low ASVs prevalence in real samples, and the minimum number of reads per sample.

Graphical plots and statistical analyses

Statistical analyses are key features of rANOMALY workflow. Main descriptive analyses are integrated thanks to phyloseq functionalities. In addition, we automatized graphical representations and advanced statistical tests for alpha, beta diversity and composition plots. Above all, we have included the four most up-to-date differential analyses to assess differentially abundant features.

Community composition plots, alpha and beta diversity analyses rANOMALY allows users to explore microbial community composition with three different types of plots : classical interactive bar plot [START_REF] Sievert | Interactive Web-Based Data Visualization with R, plotly, and shiny[END_REF] of raw and relative taxa abundances, rarefaction curves to check sampling effort, and Krona interactive pie charts [START_REF] Ondov | Krona: Interactive Metagenomic Visualization in a Web Browser[END_REF] .

Alpha diversity indices representing the specific richness are calculated (Richness, Simpson, Shannon...) with the vegan R package [START_REF] Oksanen | Community Ecology Package. R package version 2[END_REF] . We added statistical tests such as multi-factors analysis of variance and pairwise Wilcoxon tests to assess significant differences between tested categories. We included repeated measures ANOVA to handle within-subjects variation. For example, it can be used when there are measures taken on the same individuals at different time points. This step outputs graphical representations and tables with results of statistical tests which are both saved in files and returned as list objects for markdown reports.

Beta diversity analysis allows users to estimate the community differences between two samples, it is also based on the vegan package [START_REF] Oksanen | Community Ecology Package. R package version 2[END_REF] . rANOMALY can calculate all different distances such as the one based on ASVs abundance (BrayCurtis), rank based Jaccard indices, and phylogenetic distances as UniFrac, weighted Unifrac. Graphical representations PCoA, NMDS and more are available. This analysis can be processed at different taxonomic levels and categorical factors chosen by the user. The additional statistical test of PERMANOVA uses the distance matrix to determine if microbial communities of sample groups are significantly different from each other. We added a pairwise PERMANOVA [START_REF] Arbizu | pairwiseAdonis: Pairwise multilevel comparison using adonis[END_REF] test to confirm significant differences between specific group of samples.

Differential analysis

Differential analyses are meant to assess potential differentially abundant taxas between tested conditions chosen by users. rANOMALY wraps three methods:

• DESeq2 19 uses negative binomial generalized linear model with a variance stabilizing transformation on abundances.

• metagenomeSeq 20 uses a zero-inflated log-normal model with cumulative sum scaling normalization.

• metacoder [START_REF] Foster | Metacoder: An R package for visualization and manipulation of community taxonomic diversity data[END_REF] applies a total sample sum normalisation and uses a non-parametric Wilcoxon Rank Sum test to compare the log ratio of mean proportions.

The use of multiple methods for differential analysis allows the user to investigate which feature can be considered as differentially abundant between conditions. rANOMALY function using DESeq2 and metagenomeSeq outputs tables and plots with significant features. Metacoder outputs heat tree plots allowing users to infer differentially abundant features at each taxonomy rank and their position on the phylogenetic tree. With basic data management, and according to the identification of all significant differentially features (i.e. ASVs, genus, species) from the three methods, we generate a single table to find significant features in one or more methods and ease the interpretation. Additional information are added to the final table, like mean relative abundance for each feature and condition. A column in which condition features are significantly more abundant, features taxonomies and related sequences are added to the aggregated table. To complete differential analysis, we have included the well recognized PLS-DA from the mixOmics package [START_REF] Rohart | mixOmics: An R package for 'omics feature selection and multiple data integration[END_REF] . It is a supervised classification method allowing users to identify features discriminating the sample groups.

Additional features and export

Functions and procedures are available to help the user to generate additional figures or to export the data to third party software. For instance, shared taxa between conditions are useful to explore. We use a function that can generate Venn diagrams, or for more complex visualisation, we can produce files readable by Cytoscape [START_REF] Shannon | Cytoscape: a software environment for integrated models of biomolecular interaction net-works[END_REF] to produce shared taxa networks. Krona diagrams can be displayed to explore sample microbial composition, where samples can be merged by a specified factor. rANOMALY allows users to generate inputs for STAMP [START_REF] Parks | STAMP: statistical analysis of taxonomic and functional profiles[END_REF] , which is a graphical software that provides statistical hypothesis tests and exploratory plots for analysing taxonomic profiles.

Operation rANOMALY requires R 3.6.3 or upper and can be run on any operating system with common specifications (1Go disk space, 4Go RAM, multicore CPU is recommended).

Use case

For this example we are using a dataset from Fretin et al. study [START_REF] Fretin | Bacterial community assembly from cow teat skin to ripened cheeses is influenced by grazing systems[END_REF] in which samples are from four different environments: cow milk, cow cheese (rind and core) and cow teat skin. This dataset and metadata are available on NCBI-SRA website: BioProject accession PRJNA421256. To ease access to this dataset, fastq files along with pre-formatted metadata are available on this repository. Main outputs of this function are:

Install

• read_tracking.csv summarizes the read number after each filtering step (Table 1).

• raw_otu-table.csv the raw ASV table.

• rep-seqs.fna fasta file with all representative sequences for each ASV.

• We share the latest databases we use in the IDTAXA format in this link. Users can also generate your own IDTAXA formatted database following those instructions and scripts we provide at this page.

# Taxonomy assignment step tax.table = assign_taxo_fun(dada_res = dada_res, id_db = c("PATHTO/SILVA_SSU_r138_2019.RData", "PATHTO/DAIRYdb_v1.2.4_20200604_IDTAXA.rdata"), verbose = 3, output = "./analysis/02_idtaxa_silva_ddb")

Main file outputs:

• robjects.Rdata with taxonomy in phyloseq format in tax.table object.

• final_tax_table.csv the final assignation table (tax.table) outputed in CSV format.

• allDB_tax_table.csv raw assignations from the two databases, mainly for debugging.

Phylogenetic tree.

The phylogenetic tree from the representative sequences is generated using phangorn and DECIPHER packages.

tree = generate_tree_fun(dada_res)

Phyloseq object. To create a phyloseq object, we need to merge four objects and one file:

• the ASV table otu.table and the representative sequences seqtab.nochim in dada_res variable.

• a taxonomy table (tax.table).

• the phylogenetic tree (tree).

• metadata from from csv file.

# Phyloseq object generation data = generate_phyloseq_fun(dada_res = dada_res, tax.table = tax.table, tree = tree, metadata = "./sample-metadata.csv", output = "./analysis/03_phyloseq/") Decontamination. The decontam_fun function uses decontam R package along with control samples (PCR control) to filter out contaminants. The decontam package offers two main methods, frequency and prevalence (users can also combine those methods). For frequency method, it is mandatory to have the DNA concentration of each sample in the phyloseq object (and hence in the metadata.csv). The prevalence method does not need DNA quantification, this method allows to compare presence/absence of ASV between real samples and control samples and then identify contaminants.

Tips: sequencing plateforms often quantify the DNA before sequencing, but do not usually give the information. Just ask for it ;).

Our function integrates the basic ASV frequency

_ _ _ _ nb reads ASV Freq nb total reads =      
and minimum prevalence in overall samples filtering. We have also included an option to filter out ASV based on their taxa names for known laboratory recurrent contaminants.

Main outputs:

• robjects.Rdata with contaminant filtered phyloseq object named data.

• Exclu_out.csv list of filtered ASVs for each filtering step.

• Krona plot before and after filtering.

• raw_asv-table.csv and relative_asv-table.csv.

Here we are going to filter out ASVs representing less than 0.1% (freq) of the reads and that are present in more than 4 samples (prev). Moreover, we are excluding "unassigned" taxas for this use case. Our sample dataset do not contain control samples, this step will be skipped.

data_filtered = decontam_fun(data = data, domain = "Bacteria", skip = TRUE, number = 4000, freq = 0.001, prev = 4, unassigned = TRUE, krona = TRUE, output = "./analysis/04_decontam/")

We obtain the final phyloseq object used for downstream analysis: Plots, diversity and statistic analyses Rarefaction curves. In order to observe the sampling depth of each sample we start by plotting rarefaction curves. Those plots are generated by plotly which makes them interactive. (Figure 2) rareplot = rarefaction(data_filtered, "source_location", 100)

Community composition plot. The bars_fun function allows user to generate interactive community composition plot. Figure 3 presents the composition plot with relative abundances for the top 20 genera existing in our samples. The function allows to plot at different taxonomy rank and to modify the number of taxa to show.

Here two arguments controlling the composition plot aesthetics:

• Ord1 option order the sample along the X axis.

• Fact1 option control labels of the X axis. Fact1="sample.id" if user don't want the sample to be renamed.

source_location factor shows very different bacterial community between milk, cheese and cow teats environments. from others. Figure 4 shows boxplots of diversity indices, cow teats environment has much more ASV than other environments. Shannon index reveals more differences between cheese rinds and cheese cores. Cheese rinds show a higher Shannon index highlighting a more balanced bacterial community.

divAlpha = diversity_alpha_fun(data = data_filtered, output = "./analysis/07_all_res/plot_div_ alpha/", column1 = "source_location", measures = c("Observed", "Shannon"))

Results of the analysis of variance and pairwise wilcoxon test on Shannon index: The alpha_diversity_fun function returns a list which contains:

# > divAlpha$Shannon # $anova # Df
• boxplots comparing conditions with chosen indices. ($plot)

• a table of indices values. ($alphatable) And for each of the computed indices :

• an ANOVA analysis. (${measure}$anova)

• a pairwise wilcox test result comparing conditions and giving the pvalue of each comparison tested: ${measure}$wilcox_col1: wilcox test results on the first or unique factor, ${measure} $wilcox_col2_fdr: wilcox test results on the second factor, ${measure}$wilcox_col2_ collapsed: wilcox test results on collapsed factor 1 and factor 2.

• a mixture model if your dataset includes repeated measures, ie. column3 option.

(${measure}$anovarepeat, ${measure}$mixedeffect)

Beta diversity. The diversity_beta_light function allows to generate specific tests and figures ready to publish in rmarkdown report as in the example below. It is based on the vegan package function vegdist for the distance calculation and phyloseq-extended in addition to ordinate funtion for the ordination plot.

We include statistical tests to ease the interpretation of results. A permutational ANOVA is carried out on matrix distance to compare groups by testing if centroids and dispersion are equivalent for all groups. User have to inform col argument and optionally cov (covariable) to assess PERMANOVA to determine significant differences between groups. A pairwise-PERMANOVA is processed to determine which condition is significantly different from another (based on p-value).

As a return, you will get a list that contains:

• An ordination plot ($plot).

• The permANOVA results ($permanova).

• The pairwise permANOVA ($pairwisepermanova)

Here we present results of beta diversity analysis on source_location factor, Figure 5 show ordination plot and it confirms the big differences between community of cheese, milk and cow teats.

divBeta = diversity_beta_light(psobj = data_filtered, output = "./analysis/07_all_res/plot_div_beta/", col = "source_location", cov = NULL, dist0 = "bray", ord0 = "MDS")

The permanova tests on BrayCurtis distance shows significant p-value for the source_location factor. Pairwise permanova test is used to define which level of the factor tested is significantly different from others.

Differential analysis.

We choose three different methods to process differential analysis which is a key step of the workflow. The main advantage of the use of multiple methods is to cross validate deferentially abundant taxa between tested conditions. For this use case, we choose to focus on milk and cow teat environment to compare community at genus level.

Metacoder Metacoder is the most simple differential analysis tool of the three. Counts are normalized by total sum scaling to minimize the sample sequencing depth effect and it uses a Kruskal-Wallis test to determine significant differences between groups. The metacoder_fun function allows the user to choose the taxonomic rank, which factor to the test (column1), and a specific pairwise comparison (comp) to launch the differential analysis.

It produces pretty graphical trees, representing taxas present in both groups and coloring branches depending on which group this taxa is more abundant (Figure 6). Two trees are produced, a raw one, where everything is displayed and a filtered one where only significant features are represented (p-value <= 0.05). # Metacoder mtc = metacoder_fun(data = data_filtered, output = "./analysis/07_all_res/DA/metacoder", column1 = "source_location", rank = "Genus", signif = TRUE, plottrees = TRUE, min = "10", comp = "milk ~teat")

Main output is a list with :

• for each comparison (${comparison}), two heattree, one with all features (${comparison}$raw) and an other one with only significant features (${comparison}$signif).

• a table with all wilcoxon test results ($table).

DESeq2 DESeq2 is a widely used method, primarily for RNAseq applications, for assessing differentially expressed genes between controlled conditions. Its use for metabarcoding datas is sensibly the same and well documented. The deseq2_fun allows to process differential analysis as metacoder_fun, and users can choose the taxonomic rank, the factor to test and which condition to compare. DESeq2 algorithm uses negative binomial generalized linear models with VST normalization (Variance Stabilizing Transformation).

Main output is a list with:

• a plot showing Log2FoldChange value of each significant feature (${comparison}$plot).

• a table with statistics (LogFoldChange, pvalue, adjusted pvalue...) (${comparison}$table).

# DESeq2 deseq = deseq2_fun(data = data_filtered, output = "./analysis/07_all_res/DA/deseq/", column1 = "source_location", verbose = 1, rank = "Genus", comp = "milk ~teat")

MetagenomeSeq MetagenomeSeq uses a normalization method able to control for biases in measurements across taxonomic features and a mixture model that implements a zero-inflated Gaussian distribution to account for varying depths of coverage. As deseq2_fun, metagenomeseq_fun returns a table with statistics and a plot with significant features for each comparison. # MetagenomeSeq mgSeq = metagenomeseq_fun(data = data_filtered, output = "./analysis/07_all_res/DA/metagenomeseq/", column1 = "source_location", verbose = 1, rank = "Genus", comp = "milk ~teat")

Results aggregation step

The aggregate_fun function allows to merge the results from the three differential analyses methods computed previously to obtain one unique table with all informations of significant differentially abundant features. Figure 7 shows the most significant and differential abundant Genera between the two environments.

# Aggregate results resF = aggregate_fun(data = data_filtered, metacoder = "./analysis/07_all_res/DA/metacoder/metacoder_signif_Family.csv", deseq = "./analysis/07_all_res/DA/deseq/", mgseq = "./analysis/07_all_res/DA/metagenomeseq/", output = "./analysis/07_all_res/DA/aggregate_diff/", column1 = "source_location", column2 = NULL, verbose = 1, rank = "Genus", comp = "milk ~teat")

The generated table include the following fields:

• seqid: ASV ID.

• Comparaison: Tested comparison.

• Deseq / metagenomeSeq / metacoder: differentially abundant with this method (0 no or 1 yes).

• sumMethods: sum of methods in which feature is significant.

• DESeqLFC: Log Fold Change value as calculated in DESeq2.

• absDESeqLFC: absolute value of Log Fold Change value as calculated in DESeq2. • MeanRelAbcond1 / MeanRelAbcond2: Means relative abundance in condition 1 and 2.

• Condition: in which the mean feature relative abundance is higher.

• Taxonomy and representative sequence.

Here is an overview of the aggregate_fun table informing in which methods each feature is significant, their DESeq2 LogFoldChange value, taxonomy and representative sequences:

Miscellaneous function

Heatmap. User can generate an interactive heatmap with heatmap_fun function to explore relative abundance of top taxa through samples as showed in Figure 8.

# heatmap heatmap_plot = heatmap_fun(data = data_filtered, column1 = "source_location", top = 100, output = "./analysis/07_all_res/plot_heatmap/", rank = "Genus") Venn diagram of shared taxa. Figure 9 shows Venn diagram comparing the three environment of this study, this method is useful to determine shared taxa between group of samples. This function uses venn::venn function which handles up to 7 groups comparison.

# venn outvenn = ASVenn_fun(data = data_filtered, output = "./analysis/07_all_res/ASVenn/", rank = "Genus", column1 = "source_location", shared = TRUE)

ASVenn_fun generate a table allowing user to find which taxa are shared between conditions. for (i in c("ASV", "Species", "Genus", "Family")) { phy2tsv_fun(data_prune, output = glue::glue("./analysis/08_tsv_table/{i}_ table/"), rank = i) } export_to_stamp_fun function creates two files that can be imported into the STAMP software.

export_to_stamp_fun(data = data, output = "./analysis/stamp/", correc = TRUE) csv2phyloseq_fun function allows user to import data from other bioinformatic pipeline like FROGS, Qiime2. This function needs tabulated ASV, taxonomy, metadata and DNA sequence table as inputs. It can generate phylogenetic tree if missing and output a phyloseq object ready for downstream analyses.

Conclusions

rANOMALY allows users to handle metagenomic data from raw sequences quality control to final differential analysis with ready to publish results in an easy and reproducible manner. Users have access to all sources of the rANOMALY package that can be deployed on any operating system or server allowing them to analyse anything from a few samples to several thousand. This workflow combines all of the latest developments in the field: the use of high resolution amplicon sequence variant, contaminant filtering, double automatic taxonomic assignation, integrated statistical analyses and four differential analyses with cross validation. rANOMALY help users those who don't have big programming skills. And since rANOMALY uses phyloseq objects and standards, original functions from the phyloseq package can still be used to split, filter and select specific samples previously to be visualized and/or tested with rANOMALY functions. For exploratory analysis and interactive experience, a shiny application based on this workflow is under active development and will be shared with the community.

Software availability

Up-to-date source code, and tutorials are available at: https://forgemia.inra.fr/umrf/ranomaly. Package documentation is also provided at: https://umrf.pages.mia.inra.fr/ranomaly/use_case.html Archived source code as at time of publication are available from: https://doi.org/10.5281/zenodo.4338833 26 License: Creative Commons Attribution 4.0 International

Here we show all packages used in this workflow and their version numbers: The manuscript itself is very clear, nicely outlines how the different components of the pipeline are connected, and also clearly distinguish between which functions are provided by this pipeline and which functions are provided by third-party tools or packages. My only major comment is about user input (either user questions or bug reports), for which there is currently no process or web page mentioned in the manuscript. For instance, should users report bugs in the gitlab repo ( https://forgemia.inra.fr/umrf/ranomaly/-/issues)? (Incidentally, I did encounter an unexpected error with DeSeq: Error in `assays<-`(`*tmp*`, value = `*vtmp*`) : please use 'assay(x, withDimnames=FALSE)) <-value' or 'assays(x, withDimnames=FALSE)) <-value' when the dimnames on the supplied assay(s) are not identical to the dimnames on DESeqDataSet object 'x' In addition: Warning message: In DESeqDataSet(se, design = design, ignoreRank) : some variables in design formula are characters, converting to factors).

sessionInfo() # R version 3.6.
If the authors do intend users to submit gitlab issues, this should be specified in the conclusion and/or the software availability sections. If not, the alternative(s) path to provide feedback and report issues should be clearly described.

Other than that, I only have minor comments on the manuscript which are listed below: Abstract: "uses amplicon sequence variants (ASV) level for". Shouldn't this read "uses amplicon sequence variants (ASV) for"?

In the keywords, I would argue "metagenomics" should be "metabarcoding" ○ p. 3 Introduction:

○ "Studies of microbial communities tends to" should be "Studies of microbial communities tend to" "(FROGS 1 , qiime 2 )": it should be made clear that these are examples of softwares but many more exist, e.g. by changing to "(e.g. FROGS 1 , qiime 2 )", or explicitly stating "such as FROGS 1 or qiime 2 ".

"allowing to easy implementation of" should be "allowing for easy implementation of" or "allowing to easily implement".

p. 3 Raw sequence processing: "This improves" should be clarified as "This approach improves" or "This tool improves" (to make it obvious that "this" refers to dada2).

"an additional step of cutadapt" should be "an additional step using cutadapt".

○ p. 4 Taxonomic assignment: "with representing sequences" should probably be clarified as "with representative sequences for each ASV" (since the concept of "representative sequences" has not been described earlier)

"IDTAXA formatted" should be "IDTAXA-formatted" ○ p. 4 Phylogenetic tree: "with AlignSeqs function from DECIPHER package" should be "with the AlignSeq function from the DECIPHER package"

"computes the likelihood of the phylogenetic tree" This may suggest that a tree exists and that this function evaluates the likelihood of this existing tree, but I believe the first function actually builds a tree? So there may be a clause missing here, e.g. "And finally, the neighbour joining function is used to build a phylogenetic tree based on the distance matrix, and the pml function computes the likelihood of the phylogenetic tree."? ○ p. 5 Install:

The authors may want to add a "library(devtools)" line after the install of the package devtools and before installing the package ranomaly, in case some users already have the devtools package installed.

The authors should also specify that a user will need to replace "PATHTOWORKINGDIRECTORY" with the path to their working directory on their own computer. (The same is true in the Taxonomic assignment section when users have to specify a path to the taxonomic databases).

○ p. 6 Processing of raw sequences: "files name match samples name" should be "file names match sample names" ○ p. 6 Taxonomic assignment: "givenby" should be "given by" "different genus" should be "different genera" "correction step" should be removed, or clarified as e.g. "as part of a correction step" "your own" should be "their own" p. 8 Rarefaction curves: On my version of R, I have to call "rareplot" after the line "rareplot = rarefaction(data_filtered, "source_location", 100)" otherwise no plot appears. If this is the case in most versions (which I believe it is), this line ("rareplot") should be included in the authors' instructions so that a user would not be confused to not see any plot after calling the line "rareplot = rarefaction(data_filtered, "source_location", 100)" (this is also true for the other interactive plots, e.g. community composition, alpha and beta diversity, etc).

○ p. 8 Community composition plot: "source_ location": there seems to be an extra space between "source_" and "location" here ○ p. 10 Beta diversity: "funtion" should be "function"

"on matrix distance" should be "on the distance matrix" "A pairwise-PERMANOVA is processed" should be "A pairwise-PERMANOVA is performed" ○ p. 11 Differential analysis: "of the use of multiple" would be more clear as "of using multiple" "cross validate" should be "cross-validate"

"which factor to the test" should be "which factor to test" or "which factor to perform the test on" "which group this taxa is more abundant" should be "which group this taxa is more abundant in" ○ p. 12 DESeq2: "and well documented": Could the authors add 1 or 2 references here? ○ p. 15 Venn diagram: Shouldn't "the three environment" be "the four environments". ○ p. 15 Export: "function allows user" should be "function allows users" (also for csv2phyloseq_fun).

○ p. 15 Conclusion: "metagenomic data" should be "metabarcoding data"

"help users those" should be "help users"

"big programming skills" should be "advanced programming skills" "previously to be visualized and/or tested" should be "prior to visualization/test" or "before being visualized and/or tested"

Is the rationale for developing the new software tool clearly explained? Yes

Is the description of the software tool technically sound? Yes

Are sufficient details of the code, methods and analysis (if applicable) provided to allow replication of the software development and its use by others? Yes

Is sufficient information provided to allow interpretation of the expected output datasets and any results generated using the tool? Yes

Are the conclusions about the tool and its performance adequately supported by the findings presented in the article? Yes

Competing Interests: No competing interests were disclosed.

Reviewer Expertise: metagenomics, microbial ecology, phage genomics, viral ecology I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.

from fastq treatment (contaminant filtering, taxonomic assignment, alpha and beta diversity indexes) to statistical analysis (using three main differential analysis methods, DESeq2, Metacoder, metagenomeSeq) in an automatic, step by step, and reproducible way. Thus, this is a suitable workflow for exploratory and educational purposes. I was pleased to be able to reproduce all analyses and figures following the use case of the paper.

Hence, I only have a few minor issues that should be discussed:

The authors mentioned that the developed package is customizable. This is an important aspect given that projects are different from one to the other, and bioinformatic tools evolve tremendously fast so it will be crucial to adapt/update the workflow accordingly. The workflow actually allows this, however, it would be useful (particularly for beginners) to provide a vignette/tutorial on the gitlab page (and/or in a section of the paper) showing how this can be done with simple examples. The use case is based on Illumina sequences, a use case dealing with Ion torrent sequences could be relevant in this context.

1.

The use case is based on a dataset from Fretin et al. 2018, it would be relevant to interpret the results of the workflow, and to describe to what extent they match or complement the results of the previous analysis. This would be relevant to benchmark the workflow and highlight the novelties brought by this analysis.

2.

The access to the use case fastq sequences and to the IDTAXA formatted databases is given via links in the text. It would be convenient for users working with remote servers to download it with command lines. Please provide command lines to do so (via wget for instance), either in the text or on the gitlab page.

3.

There is a typo on page 6, the IDTAXA formatted database that can be downloaded from the link is DAIRYdb_v1.2.4_20200603_IDTAXA.rdata and not DAIRYdb_v1.2.4_20200604_IDTAXA.rdata.

4.

Page 8: In the "bars_fun" function, additional space has been introduced between "source_" and "location" leading to an error if one runs this piece of code.

5.

Page 9: In the "diversity_alpha_fun" function, the font has been modified for the closing double quotes after "Observed" and "Shannon". Thus, this leads to an error if one runs this piece of code.

6.

Page 15: In the "phy2tsv" function, the phyloseq object "data_prune" was never created upstream. Would it be better to test the function "phy2tsv" with the phyloseq object "data" previously created in the use case?

7.

Is the rationale for developing the new software tool clearly explained? Yes

Is the description of the software tool technically sound? Yes

Are sufficient details of the code, methods and analysis (if applicable) provided to allow replication of the software development and its use by others?

Figure 1 .

 1 Figure 1. rANOMALY workflow.

#

  Relative abondance prel = bars_fun(data_filtered, "Genus", top = 20, Ord1 = "source_location", Fact1 = "source_ location", relative = TRUE, outfile = "plot_compo_rel.html") Alpha diversity. The alpha_diversity_fun function can computes various alpha diversity indexes. It uses the estimate_richness function from phyloseq (Available measures : Observed, Chao1, ACE, Shannon, Simpson, InvSimpson, Fisher). Here we calculate ASV richness and Shannon index and carry out an analysis of variance on the source_location factor. Sequencing depth is automatically taken into account in this test. A pairwise wilcoxon test is added to ANOVA to define which group might be significantly different

Figure 2 .

 2 Figure 2. Rarefaction plot for all samples, separated by source_location factor.

Figure 3 .

 3 Figure 3. Composition plot of relative abundance of 20 most abundant genus.

Figure 4 .

 4 Figure 4. Boxplot of diversity alpha indexes.

Figure 5 .

 5 Figure 5. Ordination plot: MDS on Braycurtis distance.

Figure 6 .

 6 Figure 6. Heattree of significant features generated by metacoder package.

Figure 7 .

 7 Figure 7. Top significant and differentially abundant genera between Milk and Teat samples.

Figure 8 .

 8 Figure 8. Heatmap of top genuses relative abundance.

Figure 9 .

 9 Figure 9. Venn plot of shared genuses.

ASV definition with DADA2. The

  Up to date code is hosted by the INRAE gitlab, users can simply download and install the package in R console with following command lines. first step is to define ASVs thanks to the dada2 package. In rANOMALY, only one function is needed to compute all the different steps require from this package. Here sample names are extracted from the file name, thus be sure that files name match samples name in the metadata file.

	Processing of raw sequences
	# ASV resolution step, raw ASV table and representative sequence generation dada_res = dada2_fun(amplicon = "16S", path = "./fastq_files", outpath = "./analysis/01_results_dada/", plot = FALSE, compress = FALSE, verbose = 1, paired = FALSE, torrent_single = FALSE, returnval = TRUE)
	# Install ranomaly install.packages("devtools") devtools::install_git("https://forgemia.inra.fr/umrf/ranomaly")
	# Call ranomaly package library(ranomaly) setwd("PATHTOWORKINGDIRECTORY")
	# Create the analysis directory dir.create("./analysis")

  robjects.Rdata with saved dada_res list containing raw ASV table and representative sequences in objects otu.table, seqtab.export, seqtab.nochim.

# Adjust sample names to match with metadata colnames(dada_res$otu.table) = tools::file_path_sans_ext(colnames(dada_res$otu.table

)

)

Taxonomic assignment. assign_taxo_fun function uses IDTAXA function from DECIPHER package, and allows to use two different databases. It keeps the best assignment on two criteria, resolution (depth in taxonomy assignment) and confidence (value givenby IDTAXA). The final taxonomy is validated by checking for multiple ancestors (i.e. same species assigned to different genus) and incongruities (i.e. empty fields or incomplete lineage) correction step.

Table 1 . Read tracking table overview. sample.id input filtered denoisedF nonchim

 1

	SRR6365127.fastq 35690 33511	33402	33402
	SRR6365128.fastq 38871 38183	38111	38046
	SRR6365129.fastq 50970 49765	49631	49092
	SRR6365130.fastq 44207 42657	42524	42524
	SRR6365131.fastq 68008 66517	66257	66087
	SRR6365132.fastq 44014 42158	41945	41848
	SRR6365133.fastq 36704 35884	35799	35444
	...	...	...	...	...

  Export and compatibility phy2cyto_fun allows to generate input files (SIF format) for Cytoscape which is useful to visualise shared taxa and easily modify each nodes and arrows position and aesthetic.

	outvenn$TABf #table with shared taxa, taxonomies and sequences
	as.data.frame(head(apply(outvenn$TABf, 2, stringr::str_sub, 1, 10))) # teat milk cheese_rind cheese_core ASV taxonomy seq # 1 1 1 1 1 a0bc0ee9c9 k__Bacteri ATCTTCGGCA # 2 1 1 1 1 13b622bf4d k__Bacteri ATATTGCACA # 3 1 1 1 1 15875b7f67 k__Bacteri ATATTGCACA # 4 0 0 1 1 fe8453c0a4 k__Bacteri ATATTGCACA # 5 1 1 1 1 907e61bc24 k__Bacteri ATATTGCACA # 6 0 0 1 1 fe1733ab7e k__Bacteri ATATTGCACA

phy2tsv_fun function allows user to generate tabulated format abundance table at different taxonomic rank with following commnands:

F1000Research 2021, 10:7 Last updated: 09 AUG 2021

Acknowledgements

Authors deeply thank Philippe Ruiz (INRAE, MEDIS) and Marco Meola (Agroscope) for testing and feedback on the rANOMALY package.

The author(s) declared that no grants were involved in supporting this work.

Yes

Is sufficient information provided to allow interpretation of the expected output datasets and any results generated using the tool? Yes

Are the conclusions about the tool and its performance adequately supported by the findings presented in the article? Yes Competing Interests: No competing interests were disclosed.

Reviewer Expertise: Genomics, Bioinformatics, Evolution I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.

The benefits of publishing with F1000Research:

Your article is published within days, with no editorial bias • You can publish traditional articles, null/negative results, case reports, data notes and more •

The peer review process is transparent and collaborative • Your article is indexed in PubMed after passing peer review • Dedicated customer support at every stage • For pre-submission enquiries, contact research@f1000.com