The inclusive radiative decay of the $B$ meson is known to provide strong constraints on many popular extensions of the Standard Model. Such constraints crucially depend on precision of the Standard Model predictions. One of the main contributions to the theoretical uncertainty is due to certain Next-to-Next-to-Leading Order QCD corrections whose values at the physical charm quark mass $m_c$ have been estimated using interpolation between the $m_c=0$ and $m_c \gg m_b$ limits. A direct determination of such corrections at the physical value of $m_c$ requires calculating hundreds of two-scale four-loop propagator diagrams with unitarity cuts. Applying the integration-by-parts method, we express the corrections in terms of master integrals. Asymptotic expansions of these integrals at $m_c \gg m_b$ serve as boundary conditions for differential equations in $z=m_c^2/m_b^2$ that are being numerically solved. Here, we present our final results for the diagrams involving massless and massive fermion loops on the gluon lines. For the two-body cuts, we confirm the analytical expressions and/or numerical fits that are already present in the literature. In the four-body case, we make the correction complete by including several diagrams that have previously been only estimated using interpolation in $m_c$. We also report the status of the ongoing calculation of the remaining diagrams where no closed fermion loops on the gluon lines are present.