How can we define BF theory on a general 4-manifold? (I have rewritten the question some, with new understanding)
4d BF theory is classically presented as the TFT arising from the Lagrangian
$B\wedge F$,
where $B$ is an abelian 2-connection (locally a real 2-form) and $F$ is the curvature of a connection $A$. People often throw about this theory on general 4-manifolds, but it seems the naive definition from the Lagrangian above is anomalous for most situations. If one restricts to the simpler case of $A$ and $B$ being globally defined forms (not considering them as connections) the phase space is $H^1(\Sigma_4,\mathbb{R}) + H^2(\Sigma_4,\mathbb{R})$ with $A$ as the first coordinate and $B$ as the second. $A$ and $B$ are conjugate from the point of view of the above Lagrangian, but this has no chance of being symplectic if the first and second Betti numbers are not equal. If $A$ and $B$ are connections, however, the only solution to the equations of motion (which not only set $dA=dB=0$ but also their holonomies) is the trivial one, and it is not a problem
However, what I'd really like to consider is $nB\wedge F$. With $A$ and $B$ full-fledged connections. Let us first integrate out $B$. $B$ is a (2-)connection, so $dB$ can be any integral 3-form. We can thus write it as $d\beta + \sum_k m^k \lambda_k$, where $\beta$ is an ordinary real 2-form, and $\lambda_k$ form a basis of integral harmonic 3-forms. Then the action becomes $dA \wedge \beta + \sum_k n m^k A\wedge\lambda_k$. The first term sets $dA=0$ after integrating over $\beta$ and the second term sets the holonomies of $A$ to be nth roots of unity. If our 4-manifold doesn't have torsion, there are precisely $n^{b_1}$ such $A$s.
If we do the same thing, integrating out $A$ instead, we get Dirac deltas setting $dB=0$ and the holonomies of $B$ to be nth roots of unity. Again, if we have no torsion, there are $n^{b_2}$ such $B$s.
It looks like the determinants for these Dirac deltas are all 1, so from the naive calculation above, $Z_n = n^{b_1}$ but also $Z_n = n^{b_2}$. This is a problem on a general 4-manifold, but like the simpler situation above, there is no issue in quantizing on something of the form $\Sigma_3 \times \mathbb{R}$. What I think must be happening is that when I am integrating out $A$ or $B$, I am really only integrating out the non-flat parts. There should still be some more factors that correct the discrepancy above. I think I should be more careful in defining the above integral. Perhaps it is possible to define it sensibly on 4-manifolds which bound a 5-manifold a la Chern-Simons theory.
How does one define this theory on general 4-manifolds?
Thanks.
 A: 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.  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.
The way that I know the path integral counts 1-cycles is that I have an alternative formulation, the lattice gauge theory, that is easy to define and unambiguously does the counting.  However, I don't know how it works for the other formulation.  I suppose a naive guess would be to look at lattice 2-gauge theory for $B$, but this doesn't seem to work.
Ground states
One thing I can address is the issue of ground states.  In fact, I favor this approach because one doesn't have to worry about so many path integral subtleties.   All you really needs is the commutator.
Take the example you raise, $S^2 \times S^1$.  There are two non-trivial operators in this case, a $W_B$ for the $S^2$ and a $W_A$ for the $S^1$.  Furthermore, these two operators don't commute, which is why there are only $n$ ground states.  If we define $|0\rangle$ by $W_A |0 \rangle = |0\rangle$, then the $n$ states $\{|0\rangle, W_B |0\rangle, ..., W_B^{n-1} | 0 \rangle\}$ span the ground state space and are distinguished by their $W_A$ eigenvalue.  Importantly, you can also label states by $W_B$ eigenvalue, but you still only get $n$ states.  Poincare duality assures you that this counting of states will always work out no matter how you do it.
Furthermore, the operator $W_B$ has a beautiful interpretation in terms of tunneling flux into the non-contractible loop measured by $W_A$.  It's easier to visualize the analogous process in 3d, but it still works here.
You can also see the difference between the theory in 4d and 4d.  Since both $B$ and $A$ are 1-forms you have different possibilities.  The analogue of $S^2 \times S^1$ might be $S^1 \times S^1$ and this space does have $n^2$ states.  However, that is because both $A$ and $B$ can wrap both cycles.
The remarkable thing is that the $Z_n$ gauge theory formulation always gets it right.  You simply count the number of one-cycles and that's it.
Getting at other spaces
The state counting approach gets you everything of the form $Z_n(\Sigma^3 \times S^1)$, but even spaces like $S^2 \times S^2$ can be accessed.  You can either do the direct euclidean gauge theory computation, or you can regard $Z_n(S^2 \times S^2) = |Z(S^2 \times D^2)|^2$ i.e. the inner product of a state on $S^2 \times S^1 = \partial (S^2 \times D^2)$ generated by imaginary time evolution.
A: I figured this out a little while ago, so as not to leave any of you hanging, here is a (almost complete) resolution.
First, there is a lattice formulation of the 2-gauge theory that agrees with the 1-gauge theory up to the topological factor $n^{\chi(\Sigma_4)}$. One has a $\mathbb{Z}_n$ variable for each 2-cell, gauge transformations on 1-cells, and importantly, 2-gauge transformations (gauge transformations between gauge transformations) on 0-cells. The number of configurations is $n^{b_2}$. Then we divide by the volume of the gauge group, which is the number of gauge transformations divided by the number of 2-gauge transformations. So the partition function of this lattice theory is $n^{b_2-b_1+b_0}$. This gives the same answer as the $A$-theory on any 4-manifold fibered over $S^1$, in particular the 4-torus @Physics Monkey and I were discussing.
To get this answer from the path integral, I gauge fix $B$ using the Lorentz gauge $d*B=0$, introducing a Lagrange multiplier term $\langle \pi_1, \delta B\rangle$ as well as some ghost terms not involving $B$ or $\pi_1$. $\pi_1$ has gauge symmetries, being a 1-form, which I also fix, introducing a term $\langle E_0, \delta \pi_1\rangle$ and some more ghost terms which do not involve $B$, $\pi_1$, or $E_0$. Now I scale $B\rightarrow B/n$, $\pi_1 \rightarrow n\pi_1$, and $E_0\rightarrow E_0/n$ to reduce the integrand to the case with $n=1$ (that integral can be shown to be 1).
I just need to find how the path integral measure scales under these transformations. $DB$ will scale by $n$ to the power of the dimension of the $B$s being integrated over. It's important to notice that the zero modes of the integrand are precisely the harmonic 2-forms, a space of dimension $b_2$. Then if $B_2$ is the (infinite) dimension of all 2-forms, $DB$ scales by $n^{b_2-B_2}$. Terms like $n^{B_k}$ can all be removed by some suitable regularization a la Witten's paper on abelian S-duality. What we are left with is the answer above, $n^{b_2-b_1+b_0}$.
Maybe there is some way to understand the topological factor that appears?
