Different methods for solving a linear first order partial differential equation

I currently run into a linear first order partial differential equation $$\frac{\partial}{\partial t}f(q,t)=g(q)f(q,t)-a q\frac{\partial}{\partial q}f(q,t)+b h(q,t).$$

I have noticed that this kind of PDE is usually solved through the method of characteristics. Yet, due to the initial condition issues, it is not easy to find $$f(q,t)$$ explicitly. For example, you may find the general solution by typing this PDE into Mathematica or Maple without initial conditions. However, once the initial conditions are provided, softwares may not work.

Therefore, I am wondering whether it is possible to solve this through transforms. I tried the Laplace and Mellin, while they do not work out well. Any other methods are also welcome.

You may assume whatever assumptions, e.g., initial condtions, $$f(q,t)>0$$, and differentiability of $$g(q)$$, you require to find the solution.

Thanks for the beautiful proof provided by @EditPiAf. The general solution provided by Mathematica is $$f(q,t)=e^{\int _1^q\frac{g(y)}{a y}dy} \left(c_1\left[\frac{a t-\text{Log}[q]}{a}\right]+\int _1^q\frac{b e^{-\int _1^x\frac{g(y)}{a y}dy} h\left(x,\frac{a t-\text{Log}[q]+\text{Log}[x]}{a}\right)}{a x}dx\right).$$ Is it because that Mathematica take $$t$$ as a function of $$q$$ for the characteristics instead of $$q$$ as a function of $$t$$?

Let us write the PDE as follows $$f_t + aqf_q = g(q)f + bh(q,t) ,$$ where $$a$$, $$b$$ are constant parameters. The trajectory of the characteristic curves is governed by the differential equations $$\dot q = aq, \qquad \dot f = g(q)f + bh(q,t)$$ where we have introduced the parametrization $$f(q(t),t)$$ of the unknown. Solving for the coordinate $$q$$, we get the family of curves $$q(t) = q_0 e^{at}$$ starting at $$q_0$$. The ODE for $$f$$ can be solved by integrating factors. In fact, if we introduce the integrating factor $$M = e^{-\int_0^t g(q(\tau))\, d\tau}$$, then we can show that $$(Mf)^{\boldsymbol{\cdot}} = bMh$$ along the characteristic curves. Assuming the initial condition $$f(q,0) = \varphi(q)$$, direct integration in time leads to the formula $$Mf = \varphi(q_0) + b\int_0^t M(\theta)\, h(q(\theta), \theta)\, d\theta$$ at all time $$t$$. Finally, division by $$M$$ and substitution of $$q_0 = q e^{-at}$$ gives the general solution $$f(q,t) = M(t)^{-1} \left( \varphi(q e^{-at}) + b\int_0^t M(\theta)\, h(q e^{a(\theta-t)}, \theta)\, d\theta \right)$$ with $$M(t) = \exp\!\big( {-}\int_0^t g(q e^{a(\tau-t)})\, d\tau\big)$$.