As suggested in other answers and comments, you can use an exact test that takes into account the origin of the data. Under the null hypothesis that the probability of success $\theta$ is the same in both experiments,
$$P \bigl(\begin{smallmatrix}k_1 & k_2 \\ n_1-k_1 & n_2-k_2\end{smallmatrix}\bigr) = \binom{n_1}{k_1}\binom{n_2}{k_2}\theta^{{k_1 + k_2}}\left({1-\theta}\right)^{{\left(n_1-k_1\right)+\left(n_2-k_2\right)}}
$$
Notice that $P$ is not the p value, but the probability of this result under the null hypothesis. To calculate the p value, we need to consider all the cases whose $P$ is not higher than for our result. As noted in the question, the main problem is that we do not know the value of $\theta$. This is why it is called a nuisance parameter.
Fisher's test solves this problem by making the experimental design conditional, meaning that the only contingency tables that are considered for the calculation are those where the sum of the number of successes is the same as in the example ($1556 + 1671 = 3227$). This condition may not be in accordance with the experimental design, but it also means that we do not need to deal with the nuisance parameter.
There are also unconditional exact tests. For instance, Barnard's test estimates the most likely value of the nuisance parameter and directly uses the binomial distribution with that parameter. Obviously, the problem here is how to calculate $\theta$, and there may be more than one answer for that. The original approach is to find the value of $\theta$ that maximizes $P$. Here you can find an explanation of both tests.
I have recently uploaded a preprint that employs a similar strategy to that of Barnard's test. However, instead of estimating $\theta$, this method (tentatively called m-test) considers every possible value of this parameter and integrates all the results. Using the same notation as in the question,
$P \bigl(\begin{smallmatrix}k_1 & k_2 \\ n_1-k_1 & n_2-k_2\end{smallmatrix}\bigr) = \binom{n_1}{k_1}\binom{n_2}{k_2}\int_{0}^{1}\theta^{{k_1 + k_2}}\left({1-\theta}\right)^{{\left(n_1-k_1\right)+\left(n_2-k_2\right)}}d\theta$
The calculation of the p value can be simplified using the properties of the integral, as shown in the article. Preliminary tests with Monte Carlo simulations suggest that the m-test is more powerful than the other extact tests at different significance levels. As a bonus, this test can be easily extended to more than two experiments, and also to more than two outcomes. The only limitation is in the speed, as many cases need to be considered. I have also prepared an R package to use the test (https://github.com/vqf/mtest). In this example,
> library(mtest)
> m <- matrix(c(1556, 2455-1556, 1671, 2730-1671), nrow = 2,
byrow = FALSE)
> m.test(m)
[1] 0.0837938
In my computer, this takes about 20 seconds, whereas Barnard's test takes much longer.
prop.test
:prop.test(c(1556, 1671), c(2455, 2730))
. $\endgroup$