If the univariate hypergeometric is your only tool you have to get it into something where you have two classes.
One approach (not the only one):
Break the total up as follows --
Draw 2 non-blue + draw 3 non-blue + ... + draw 5 non-blue.
Then work out the probability under each case; e.g. the first one is:
$P(\text{two non-blue balls in 5 draws}) \times$
$\hspace{0.5 cm} P(\text{exactly one red} |$ $\hspace{ 3cm} \text{two balls that are either green or red drawn from the original pool})$
So for the second part, you're essentially drawing two balls from (2 red, 3 green) and working out the probability of exactly 1 red. So it should be the product of two hypergeometric probabilities.
The second term would be
$P(\text{three non-blue balls in 5 draws}) \times$
$\hspace{0.5 cm} [P(\text{exactly one red} |$ $\hspace{ 3cm} \text{three balls that are either green or red drawn from the original pool})$ $\hspace{0.3cm}+P(\text{exactly two red} |$ $\hspace{ 3cm} \text{three balls that are either green or red drawn from the original pool})]$
This is the sum of two hypergeometric probabilities, times a hypergeometric probability; however, you can write it as the difference of two phyper calls (which doesn't save anything in this term, but will on the remaining ones.
So the overall thing would be a sum of a vector of terms of the form dhyper(...)*(phyper(...)-phyper(...)). Note that you should be able to do the whole thing with a single call of dhyper and two phyper calls (since you can pass vector arguments).
where the dhyper call covers the "draw $i$ non-blue balls" and the difference of phyper terms covers the range of how many reds are drawn out of $i$.
If you do have a multivariate hypergeometric pmf, you should be able to write it as a sum of terms.
You could also approach it in terms of
$P(\bar{B}\geq 2)*
[1-P(R=0|\bar{B}\geq 2)-P(G=0|\bar{B}\geq 2)+
P(R=0,G=0|\bar{B}\geq 2)]$
This, too, will involve a sum of terms, but you can generate these using vector arguments as well.
I may come back and try to make this answer more broadly useful.
Best Answer
Yes, it has, but it doesn't seem to be well known. It has the same form as the mgf of the (univariate) hypergeometric as given in https://en.wikipedia.org/wiki/Hypergeometric_distribution, so let us start with that.
The hypergeometric probability mass function is $$ f(x) = \frac{\binom{a}{x}\binom{b}{n-x}}{\binom{a+b}{n}} $$ and then the mgf can be written $$ M(t)=\frac{_2F_1(-n,-a;b-n+1;e^t)}{_2F_1(-n,-a;b-n+1;1)} $$ and other forms can be given. Kendall's Advanced theory of Statistics gives it as the solution of a differential equation, while there is a short paper by Roanld Lessing "An Alternative Expression for the Hypergeometric Moment Generating Function" which gives it as an $n$th derivative of some related function.
Generalization to the multivariate case is given in the paper by K G Janardan & G P Patil "A Unified Approach for a class of Multivariate Hypergeometric Models" in Sankhya. First some notation.
A finite population consists of $s+1$ subpopulation, indexed $0,1,\dotsc,s$. The population size is $N$ and $N_i$ the subpopulation sizes. We are sampling without replacement , sample size $n$. $X$ is the random vector $X=(X_1, X_2, \dotsc, X_s)$. The multivariate hypergeometric (MH) pmf can then be written as $$ mh(x;n,N)= \frac{\binom{N_0}{n-\sum x_i}\prod_1^s\binom{N_i}{x_i}}{\binom{N}{n}} $$ and the multivariate moment generating function can then be written as $$ M_X(t_1,\dotsc,t_s)=\frac{N_0^{(n)}}{N^{(n)}} {}_2F_1(-n,-N_1,\dotsc,-N_s;N_0-n+1;e^{t_1},\dotsc,e^{t_s}) $$ and the symbol $\alpha^{(n)}=\alpha (\alpha-1) \dotsm (\alpha-n+1)$.
This paper gives similar information also for related distributions like mv (multivariate) multinomial, mv inverse hypergeometric, mv negative hypergeometric, mv negative inverse hypergeometric and extensions like mv Polya and mv inverse Polya.