The "right" thing to do is to treat the options as derivative contracts. Let's say for simplicity that you are using Monte Carlo to compute VaR. Then you would simulate the equity prices on each iteration, and then apply an option-pricing formula to get the corresponding option prices on that iteration. This lets you obtain an accurate simulated portfolio value.
The VaR then just comes out as the usual tail measure of the simulated distribution.
Technically, what is going on is that $\text{VaR}_b^\tau$ is a quantile of the portfolio value distribution
$$ \Pi^\tau = \sum_{i=1}^N A_i^\tau $$
where some of the instruments $A_i$ may be options. That is,
$$ \text{VaR}_b^\tau = Q_b(\Pi^\tau). $$
The Delta-Gamma approximation is giving you inaccurate values of the $A_i$ (though they are fast to compute and often "good-enough" in real-world situations).
If you look at commercial packages like RiskMetrics, they offer the user an option to use Delta-Gamma, or alternatively to price options as derivatives. In the latter case you can also simulate volatility changes and credit risk changes to get even more precise values.