Let $\{a_1,a_2,\ldots,,a_n\}$ be a random sample of a Poisson distribution. Consider the following random variables $X_1=\mathrm{Binomial}(a_1,q), ~X_2=\mathrm{Binomial}(a_2,q),\ldots,~X_n=\mathrm{Binomial}(a_n,q)$.
Is there a method to find the distribution of $X=\sum_{i=1}^nX_i$?
The answer to this particular problem is immediate if you think in terms of how these distributions usually arise.
Consider a homogeneous Poisson process of rate $\mu.$ In the first interval of time you observe $a_1$ events, then in the second interval you observe $a_2 $ events, and so on. The Binomial probabilities reflect the results of selecting randomly, with equal probabilities $q,$ from each set of $a_i$ events. The result is a Poisson process that has been "thinned" by keeping only a fraction $q$ of all events. Thus, its rate must be $q\mu$ and it has run for $n$ time intervals. Consequently, writing $\lambda=n\mu$ for the overall rate during these time intervals, the answer must be the Poisson probability,
$$\Pr(X=k) = e^{-q\lambda} \frac{(q\lambda)^k}{k!}.$$
Let's now address the general question where the $a_i$ arise independently according to any distribution.
Let $A=a_1+a_2+\cdots+a_n$ be the random variable equal to this sum. By conceiving of a Binomial$(a,q)$ distribution as being that of the sum of $a$ independent Bernoulli$(q)$ distributions, it is evident that the sum of the $X_i$ is distributed as the sum of all $a_1+a_2+\cdots + a_n$ Bernoulli variables. Let the probability function for $A$ be $p,$ so that for $a=0, 1, 2, \ldots,$
$$p_a = \Pr(a_1+a_2+\cdots+a_n=a).$$
Conditional on $A=a,$ the Binomial probability law is
$$\Pr(k) = \binom{a}{k}q^k(1-q)^{a-k}.$$
The event $X=k$ is the disjoint union of the events $(A,X)=(0,k),$ $(A,X)=(1,k),$ $(A,X)=(2,k),$ and so on. Therefore its probability is the sum of the probabilities of these events, each of which can be computed in terms of conditional probabilities as
$$\Pr(A,X)=(a,k) = \Pr(X=k\mid A=a)\Pr(A=a) = \binom{a}{k}q^k(1-q)^{a-k}\, p_a.$$
Their sum is
$$\Pr(X=k) = \sum_{a=0}^\infty \Pr(a,k) = \sum_{a=0}^\infty\binom{a}{k}q^k(1-q)^{a-k}\, p_a.$$
Notice that $X=k$ is impossible unless $a \ge k.$ Let us therefore reindex the sum in terms of $a-k=0, 1, 2, \ldots,$ writing $n=a-k:$
$$\Pr(X=k) = \sum_{n=0}^\infty\binom{n+k}{k}q^k(1-q)^{n}\, p_{k+n} = q^k\sum_{n=0}^\infty\binom{n+k}{k}(1-q)^{n}\, p_{k+n}.$$
This the general solution. It was obtained using only the basic axioms of probability (and therefore many explanations can be found in textbooks, if you wish to do further research).
When the $a_i$ have Poisson distributions of rates $\lambda_i,$ their sum has a Poisson distribution of rate $\lambda = \sum_i\lambda_i,$ whence
$$p_a = e^{-\lambda} \frac{\lambda^a}{a!}.$$
Plugging this into the foregoing formula gives
$$\Pr(X=k) = q^k\sum_{n=0}^\infty\binom{n+k}{k}(1-q)^{n}\, e^{-\lambda} \frac{\lambda^{k+n}}{(k+n)!} = e^{-\lambda}\frac{(\lambda q)^k}{k!}\sum_{n=0}^\infty\frac{(\lambda(1-q))^{n}}{n!}.$$
This last sum is the series for $\exp(\lambda(1-q)),$ simplifying the result to
$$\Pr(X=k) = e^{-\lambda}\frac{(\lambda q)^k}{k!}e^{\lambda(1-q)} = e^{-q\lambda}\frac{(\lambda q)^k}{k!},$$
which is the Poisson probability for $k$ with rate $\lambda q$ we obtained at the outset (using no calculation at all).
