Can we combine varying coefficients in GAM to get them as a whole? I'm mainly interested in an estimation of smoothed yearly trend of rate of patients by a disease type.  Currently I'm using model with varying coefficients in a form:
    fit <- mgcv::gam(no_patients ~ 
                              s(year, by=disease_type) + 
                              s(clinic_ID, year, bs="re")
                              sex + 
                              age + 
                              disease_type + 
                              offset(population), 
                              family="poisson", 
                              method="REML")

Is there a possibility to get at the same time from the model one (collective) smoothed yearly trend for all disease types?

I think you would need to use something like
gam(no_patients ~  s(year) +
                   s(year, by = disease_type, m = 1) + 
                   s(clinic_ID, bs = "re") +
                   sex + 
                   age + 
                   disease_type + 
                   offset(population), 
    family = "poisson", 
    method = "REML")

where we have what we called a "global" smoother s(year) and then specify the factor by smooth to have a 1 derivative penalty (m = 1) in order to allow better identifiability of these repeated smooths of year. This s(year) term is in some sense the smooth of time averaged over all the data (for some definition of "all" - it is still conditional upon all the other terms in the model).
I'm not sure what you were hoping to induce by the s(clinic_ID, year, bs = "re") term? This is going to give you a random linear slope of year for each clinic but not allow each clinic to have it's own intercept. I have removed this term and left only the random intercept per clinic.
If you have more than a few diseases, and you expect the trends per disease are in some way similarly smooth (same amount of wiggliness) then you could use an alternative factor smooth interaction, the fs basis:
gam(no_patients ~  s(year) +
                   s(year, disease_type, bs = "fs") + 
                   s(clinic_ID, bs = "re") +
                   sex + 
                   age +
                   offset(population), 
    family = "poisson", 
    method = "REML")

where we remove m = 1 to allow a second derivative penalty as the fs basis is fully penalised. We also remove the disease_type parametric term because the fs basis contains a random intercept per disease_type.
For more on these models, see our Open Access paper on HGAMs in PeerJ (Pedersen et al 2019).
Pedersen, E.J., Miller, D.L., Simpson, G.L., Ross, N. (2019) Hierarchical generalized additive models: an introduction with mgcv. PeerJ 7:e6876.
