Lets do this as explicitly as possible (I would like to find the constant and error term)
1 The Original Problem
First, consider the identity $\log\left(1-\frac{2}{p}+\frac{1}{p^2}\right) +\log\left(1- \frac{1}{(p-1)^2} \right)=\log\left(1- \frac{2}{p} \right)$. From this, it follows that
$$\sum_{2<p\leq n}\log\left(1-\frac{2}{p}+\frac{1}{p^2}\right) +\sum_{2<p\leq n} \log\left(1- \frac{1}{(p-1)^2} \right) = \sum_{2<p\leq n}\log\left(1-\frac{2}{p}\right)$$
Multiplying by negative one and exponentiating both sides yields
$$\left(\prod_{2<p\leq n}1-\frac{1}{p}\right)^{-2} \cdot \prod_{2<n<p} \left(1- \frac{1}{(p-1)^2} \right)^{-1} = \left(\prod_{2<p\leq n}1-\frac{2}{p}\right)^{-1} $$
Recall $ \Pi_2=\prod_{2<p} \left(1- \frac{1}{(p-1)^2} \right)$ is the Twin Prime Constant, and that the one product on the left hand side converges to the reciprocal of this. It is then by one of Mertens formulas that we know $$\left(\prod_{2<p\leq n}1-\frac{1}{p}\right)^{-1}=\frac{1}{2}e^\gamma \log n + O(1)$$ where $\gamma$ is the Euler Mascheroni Constant. (Specifically this is Theorem 2.7 (e) in Montgomery Multiplicative Number Theory I. Classical Theory) Upon squaring this asymptotic result, we are able to conclude:
$$\left(\prod_{2<p\leq n}1-\frac{2}{p}\right)^{-1} = \frac{1}{4}e^{2\gamma}\Pi_2^{-1} \log^2n + O(\log n)$$
Hope that helps,
Note: The reason I substituted the twin prime constant in for that other product is because it converges very fast in comparison with the error term. I can give more details if desired, but I'll leave it as an exercise.
2 What is best possible?
Can the error term be made better? Yes. It turns out we can make that error term a lot better. By using the Prime Number Theorem we find
$$\prod_{2<p\leq n} \left( 1-\frac{1}{p}\right)^{-1}=\frac{1}{2}e^\gamma \log n + e^{-c\sqrt{\log n}}$$ where $c$ is the constant used in the proof of the Zero Free Region. Since $e^{-\sqrt{\log n}}$ decreases faster than any power of $\log$, we obtain a much better result upon squaring this estimate. Precisely we have:
$$\left(\prod_{2<p\leq n}1-\frac{2}{p}\right)^{-1} = \frac{1}{4}e^{2\gamma}\Pi_2^{-1} \log^2n + O\left( e^{-c\sqrt{\log n}} \right)$$
(again the convergence to $\Pi_2$ is much to rapid to interfere with the error term)
I would be willing to bet that this is the best we can do, and that any better would imply stronger results regarding the error term for $\pi(x)$, the prime counting function.
3 Numerics
Just for fun, the exact constant in front of the $\log^2n$ is: 1.201303. How close is this? Well for:
$n=10$ we get an error of $0.630811$
$n=50$ we get an error of $1.22144$
$n=100$ we get an error of $0.63493$
$n=1000$ we get an error of $0.438602$
$n=10^4$ we get an error of $0.250181$
$n=10^5$ we get an error of $0.096783$
$n=10^6$ we get an error of $0.017807$
Where each time the error is positive. That is the product seems to be slightly larger than the asymptotic (but converging fairly rapidly). However, my intuition tells me it is almost certain (I will not prove it here) that the error term oscillates between negative and positive infinitely often.
4 Under Riemann Hypothesis
If we assume the Riemann Hypothesis, the error term is bounded by $$\frac{C\log^2 x}{\sqrt{x}}$$ for some constant C. By analyzing the above data with numerical methods, the error seems to be best fitted by $\frac{C\log x}{\sqrt {x}}$.
Leaving some primes out of the Euler product won't affect the location of the zeroes, since you will end up with the Zeta function multiplied by a non-zero analytic function (which won't produce any more zeros), and the formula for the prime counting function depends only on the location of the zeroes. So as far as I understand your algorithm, yes, it will "regenerate" any (finite number of) primes that were initially missing.
EDIT: Corrected my erroneous description of the missing term as "constant".
EDIT: To show that the analytic continuation of the product is the same as the product of the analytic continuation, use the fact that the analytic continuation is unique: "Let $f_1$ and $f_2$ be analytic functions on domains $\Omega_1$ and $\Omega_2$, respectively, and suppose that the intersection $\Omega_1 \cap \Omega_2$ is not empty and that $f_1 = f_2$ on $\Omega_1 \cap \Omega_2$. Then $f_2$ is called an analytic continuation of $f_1$ to $\Omega_2$, and vice versa (Flanigan 1983, p. 234). Moreover, if it exists, the analytic continuation of $f_1$ to $\Omega_1$ is unique." We will also need the fact that $\zeta(s)$ and $1-p^{-s}$ are analytic (and the more basic fact that the product of two analytic functions is analytic).
COMMENT: I think your idea here is pretty interesting! I suspect that it may work even in some cases where an infinite number of primes are discarded.
Best Answer
You can use the Taylor expansion of $ 1 - \frac{x^9}{(1+ x^4)(1+x^5)}$ to pull out as many zeta factors as you need
\begin{eqnarray} && \prod_{p} \left(1 - \frac{1}{(p^4 + 1)(p^5 + 1)} \right) \\ &=& \prod_{p} \left( 1 - \frac{1}{p^9} + \frac{1}{p^{13}} + \frac{1}{p^{14}} - \frac{1}{p^{17}} - \frac{1}{p^{18}} + O(\frac{1}{p^{{19}}}) \right) \\ &=& \prod_{p} \left(1 - \frac{1}{p^9}\right) \left(1 + \frac{1}{p^{13}} + \frac{1}{p^{14}} - \frac{1}{p^{17}} - \frac{1}{p^{18}} + O(\frac{1}{p^{{19}}})\right) \\ && \dots \\ &=& \prod_{p} \left(1 - \frac{1}{p^9}\right)\left(1 + \frac{1}{p^{13}}\right)\left(1 + \frac{1}{p^{14}}\right)\left(1 - \frac{1}{p^{17}}\right)\left(1 - \frac{1}{p^{18}}\right) \left( 1 + O(\frac{1}{p^{{19}}})\right) \\ && \frac{\zeta(13)\zeta(14)}{\zeta(9)\zeta(26)\zeta(28)\zeta(17)\zeta(18)}\prod_{p} \left(1 + O(\frac{1}{p^{{19}}})\right) \\ \end{eqnarray}
In general you can do this for euler products of the form $\prod_p f(p)$ where $f$ has a convergent expansion in $p^{-1}$. The convergence of the above formula can be slow, but a solution is to compute the product explicitly for primes less than some $N$, and then use the formula for the remainder, being sure to adjust the zeta functions for the missing factors.