BF theory secretly has another name - $Z_n$ gauge theory in the deconfined limit. The parameter $n$ is what appears in front of the $BF$ action. $Z_n$ gauge theory can be defined on any manifold you want - just introduce a lattice approximation of the manifold and compute the lattice gauge theory partition function. By taking the extreme deconfined limit of the lattice gauge theory one can verify that this construction is independent of the way you approximated the manifold.
Basic examples
In the deconfined limit the $Z_n$ flux through every plaquette of the lattice gauge theory is zero. (This is like the constraint imposed by integrating out $B$.) We have a residual freedom to specify the holonomy around all non-contractible loops. Hence $Z_n(M) = Z_n(S^4) n^{b_1(M)}$ where $Z_n(S^4)$ is a normalization constant. Requiring that $Z_n(S^3 \times S^1) =1$ gives $Z_n(S^4) = 1/n$. This condition, $Z_n(S^3 \times S^1) = 1 $, is the statement that the theory has one unique ground state on $S^3$. In general, $Z_n(\Sigma^3 \times S^1)$ is $\text{tr}(e^{-\beta H})$, but since $H=0$ we are simply counting ground states.
As a quick check, $Z_n(S^2 \times S^1 \times S^1) = n$, which is the number of ground states on $S^2 \times S^1$, and $Z_n(T^3\times S^1) = n^3$, which is the number of ground states on $T^3 = (S^1)^3$.
Further relations
Another check on the value of $Z_n(S^4)$: this is the renormalized topological entanglement entropy of a ball in the $3+1$d topological phase described by deconfined $Z_n$ gauge theory. 976
More precisely, the topological entanglement entropy of a ball is $-\ln{Z_n(S^4)}$ which gives $-\log{n}$ in agreement with explicit wavefunction calculations.
We can also consider defects. The $BF$ action is $\frac{n}{2\pi} \int B \wedge d A$. Pointlike particles (spacetime worldlines) that minimally couple to $A$ carry $Z_n$ charge of $1$. Similarly, string excitations (spacetime worldsheets) that minimally couple to $B$ act like flux tubes carry $Z_n$ flux of $2\pi/n$. Now when a minimal particle encircles a minimal flux, we get a phase of $2 \pi/n$ (AB phase), but this also follows from the $BF$ action. Without getting into two many details, the term in the action like $B_{12} \partial_t A_3$ fixes the commutator of $A_3$ and $B_{12}$ to be $[A_3(x),B_{12}(y)] = \frac{2\pi i}{n} \delta^3(x-y)$ (flat space). The Wilson-line like operators given by $W_A = e^{i \int dx^3 A_3}$ and $W_B = e^{i \int dx^1 dx^2 B_{12}}$ thus satisfy $W_A W_B = W_B W_A e^{2 \pi i /n}$ which is an expression of the braiding flux above that arises since the particle worldline pierced the flux string worldsheet.
Comments on the comments
Conservative thoughts
If I understand you correctly, what you want to do is sort of argue directly from the continuum path integral and show how the asymmetry you mentioned arises. I haven't done this calculation, so I can't be of direct help on this point right now. I'll ponder it though.
That being said, it's not at all clear to me that treating the action as $\int A dB$ leads one to just counting 2-cycles. Of course, I agree that 2-cycles are how you get non-trivial configurations of $B$, but in a conservative spirit, it's not at all clear to me, after gauge fixing and adding ghosts and whatever else you need to do, that the path integral simply counts these. 976