How to write down the collision integral for any interaction in the Boltzmann equation?
The collision term $$ {\cal I}_{coll} \simeq \frac{f-f_0}{\tau} $$ known as the relaxation time (or BGK, Bhatnagar-Gross-Krook) approximation, is a simple phenomenological model that can be motivated, but not derived, from more microscsopic theories (of course, if you compute only one number, there is always a value of $\tau$ that will give the right answer, and this is sometimes used to provide explicit formulas for $\tau$).
In general, the problem of deriving the Boltzmann equation and the collision term from microscopic theories is very complex. The simplest situation, dealt with in standard text books, arises when the system is dominated by elastic 2-body scattering. Then $$ \left.\frac{\partial f}{\partial t}\right|_{coll} = \int dv_2 \int d\Omega\, \frac{d\sigma}{d\Omega} \, |v_1-v_2| \left( f(v_1')f(v_2') - f(v_1)f(v_2)\right) $$ where $d\sigma/d\Omega$ is the cross section for $(v_1,v_2) \to (v_1',v_2')$. This cross section can be either classical (derived from classical trajectories in a potential) or quantum (the square of an amplitude, calculated, for example, using Feynman diagrams). In the classical case this follows from the Liouville equation and the BBGKY hierachy. In the quantum case this was first studied by Baym and Kadanoff, and is described in their text book. This result is for classical particles (Boltzmann statistics), but it is easily generalized to Bose or Fermi statistics, by adding Pauli blocking or Bose enhancement factors.
The literature contains a variety of extensions, for example for multiple species, inelatic collisions, three-body collisions, radiation (as in QED), etc. Deriving these collision terms from microscopic theory becomes increasingly difficult.
Finally: How to motivate the BGK appproximation? The collision term is a non-linear functional of the distribution function, $$ \left. \frac{\partial f}{\partial t}\right|_{coll}={\cal I}_{coll}[f] . $$ We can define a linearized collison kernel by expanding around the equilibrium distribution $f_0$ $$ {\cal I}_{coll}[f]={\cal I}_{coll}[f_0+\delta f]\simeq L[\delta f]. $$ One can show that $L$ is a hermitean, negative semi-definite, operator. This means that we can diagonalize $L$ $$ L = -\sum_i \frac{|\chi_i\rangle\langle \chi_i|}{\tau_i} $$ where the eigenvalues $\tau_i$ are a set of relaxation times, and the $\chi_i(v)$ are the corresponding eigenfunctions of $L$ (and I haved used bra-ket notation as in QM). If I consider the Boltzmann equation in the long-time (hydrodynamic) limit, the longest relaxation time dominates and $$ L[\delta f] \simeq -\frac{|\chi_0\rangle \langle \chi_0| \delta f\rangle}{\tau_0} \simeq -\frac{|\delta f\rangle}{\tau_0} $$ which is the BGK form.
The "quantum Boltzmann equation"
There is a semiclassical extension of the Boltzmann equation by Uehling and Uhlenbeck [1] (sometimes called the "quantum Boltzmann equation"), meant to account for how Fermi-Dirac or Bose-Einstein statistics affect the likelihood of a particular collision occurring. Their collision operator reads $$ \mathcal{I}_{\mathrm{coll}}(\vec v_1) = \int |\vec v_1-\vec v_2| \frac{d\sigma}{d\Omega} [ f_1' f_2'(1 + \eta f_1)(1 + \eta f_2) - f_1 f_2(1 + \eta f_1')(1 + \eta f_2') ] d^3v_2 d^2\Omega $$ where "1" and "2" label the velocities of the two colliding particles, primes denote the post-collision velocity, and decorations on the distribution functions $f$ are abbreviations for the velocity argument, e.g., $f_1' = f(\vec v_1')$. The quantity $\eta$ depends on the statistics of the particles: $$ \eta = \begin{cases} 1 &\text{bosons} \\ -1 &\text{fermions} \\ 0 &\text{classical} \end{cases} $$ In particular, one can see that in the classical case, one recovers the ordinary Boltzmann equation. There are two ways in which quantum mechanics appears in this collision operator.
The new factors involving $\eta$ either increase the probability of collisions between bosons or decrease the probability between fermions. In the ordinary Boltzmann equation ($\eta=0$), the quantities like $f_1 f_2$ have the interpretation $$ f_1 f_2 = \text{the probability that, prior to the collision, particle 1 has velocity $\vec v_1$ and particle 2 has velocity $\vec v_2$} $$ The new quantum-statistical factors modify this as appropriate for bosons or fermions. For concreteness, consider Fermions: $$ f_1 (1 - f_1')= \text{ the probability that, prior to the collision, particle 1 has velocity $\vec v_1$ *but not* $\vec v_1'$ } $$ (and repeat for particle 2). That second factor enforces the Pauli Exclusion Principle. It prevents collisions which would cause a colliding particle to end up in a state that is already occupied by some other particle. (In literature, this is often called "Pauli blocking".) For bosons, $\eta = +1$ casuses the opposite effect, making it preferential for particles to collide so that they end up already-occupied states.
One important property of the Uehling-Uhlenbeck equation is that the equilibrium distribution function is $$ f_0(\vec v) = \frac{1}{e^{(\frac{1}{2}m|\vec v|^2 - \mu)/k_BT} + \eta} $$ which is the Bose-Einstein, Fermi-Dirac, or Maxwell-Boltzmann distribution, according to the value of $\eta$.
[1] E. A. Uehling and G. E. Uhlenbeck, Phys. Rev. 43, 552 (1933).
The quantum cross-section
Regardless of whether one uses the quantum or classical Boltzmann equation, one should in principle always use a quantum-mechanically calculated collision cross-section. This means two things:
The cross-section should come from the quantum theory of two-particle scattering. This could mean a phase-shift calculation, a Born approximation, or whatever. Any introductory QM textbook should cover these topics. If you are working with material conditions that are hot and/or dilute, you can get away with calculating the cross-section classically.
For collisions between particles of the same species, the cross-section should account for the fact that identical particles are indistinguishable. This will either reduce or enhance the cross-section depending on whether that particles are fermions or bosons. Usually, this topic is treated in graduate-level QM textbooks. Note that there is no classical limit for this phenomenon, so even if you are treating the binary-collision problem classically, that cross-section should still be corrected to account for the identity of particles.
The quantum relaxation time
If your application is couched in the relaxation-time approximation, then QM shows up in two ways:
The equilibrium distribution function, $f_0$, is either the Fermi-Dirac, Bose-Einstein, or Maxwell-Boltzmann distribution as described in the first section.
The relaxation time should be given by a quantum cross-section, as described in the second section. That is, one may still take $$ \tau(\vec v) = \frac{1}{n|\vec v|\sigma^{(1)}(|\vec v|)} $$ provided that the momentum-transfer cross-section $$ \sigma^{(1)}(|\vec v|) = \int_0^{2\pi}\int_0^\pi (1 - \cos\theta) \frac{d\sigma}{d\Omega} \sin\theta\,d\theta\,d\phi $$ is calculated using a quantum-mechanical differential cross-section.