Bayesian modeling and FDR correction Question
Assume that we have 2 experiments, each of which yields $N$ observations. We label them $(c_1^i,\ldots,c_N^i)$ with $i=1,2$. Each observation is a positive integer, and can be modeled as a Poisson variate. We wish to know where the two experiments disagree.
Assume further that in the classical context, the appropriate analysis for this experiment would be to


*

*Perform a rate ratio test on all pairs of observations, asking whether $c_i^1$ and $c_i^2$ have the same Poisson rate ($\mathcal{H}_0$) or not ($\mathcal{H}_1$)

*Perform a FDR (Benjamin-Hochberg) correction on the $p$-values to account for multiple testing

*Select those pairs with $q<0.05$


I would like to know what is the Bayesian take on this, in particular at step 2.
Ideas


*

*I have pursued a model comparison approach, but am stuck at step 2


First, I write the probability for two observations and one rate, model $M_1$: 
$$p(M_1|c_i^1,c_i^2) = \int \text{d}\lambda p(c_i^1,c_i^2|\lambda)p(\lambda)$$
Or two rates, model $M_2$:
$$p(M_2|c_i^1,c_i^2) = \int \text{d}\lambda_1\text{d}\lambda_2 p(c_i^1|\lambda_1)p(c_i^2|\lambda_2)p(\lambda_1)p(\lambda_2)$$
I use a proper prior for $p(\lambda)$, which I don't want to discuss here.
Then, I compute the Bayes factor for these two models: $K_i = p(M_2|c_i^j)/p(M_1|c_i^j)$ assuming equal prior odds. The Bayes factor for all observations then factorizes into the individual pairs: $K=K_1\cdots K_N$
Now, the big question, is how do we assign most probable models to each pair?
a) Maximize K, e.g. give 1 rate to points for which $K_i<1$ and 2 rates to the others. To me, this sounds like no $p$-value correction for multiple testing, i.e. a bad idea.
b) Seek a combination to reach, say, $K\simeq 3$. But how do I do that? That sounds much closer to FDR to me.


*

*Write down the posterior probability of one rate being larger than the other, and use this probability in a setting similar to the definition of a $p$-value. I did not try this yet but thought I could write it here.

 A: I am not sure to fully answer your question but think that the following can help.
What you can try is to a design multilevel model including a partial pooling of the $\lambda_i$ by considering for example:
$$ \begin{align} c_i &\sim P(\lambda_i) \\ λ_i &\sim Gamma(a,b) \\ a &\sim \text{Find a good prior here} \\ b &\sim findTheGoodPriorHere \end{align}$$
Or another possibility:
$$\begin{align} c_i &\sim P(\lambda_i). \\ log(\lambda_i) &=\alpha + \theta_i \\ \mbox{with } \theta_i &\sim N(0,\sigma^2) \\ \alpha &\sim u_{R} \\ \sigma &\sim gamma(\epsilon,\epsilon) \end{align}$$
or any other suitable pooling, so that the margin posterior for the $\lambda_i$:
$$
p(\lambda_i|(c_j)_{j=1:N}) \propto \int_{R}[ p(c_i|\lambda_i) p(\lambda_i|\beta) d\beta ] \cdot \prod_{j\ne i} [\int_{R^+} \int_{R} p(c_j|\lambda_j) p(\lambda_j|\beta) d\beta  d\lambda_j]
$$
(calling for generality $\beta$ the hyperparameter of $p(\lambda_j|\beta)$)
are no longer independant each other (because the hyperparameter $\beta$ is common to all the $\lambda_i$) and does no more result in independent inferences/comparisons (while this needs a dedicated discussion e.g. Why don't Bayesian methods require multiple testing corrections?).
The important question is  weither or not the partial pooling model is suitable for your design. 
Here (http://www.stat.columbia.edu/~gelman/research/unpublished/multiple2.pdf) is a reference for such a consideration by Gelman (it does not use Bayes factor but to my knowledge, there is no limitation in using such a pooling model with Bayes factor)
