The traditional least-squares fitting or chi-squared minimisation route of fitting a straight line makes the implicit assumption that the errors on the x-axis quantity are negligible. If that is so, then there is no reason why you can't use the uncertainty in the gradient as the uncertainty in $R$.
I guess from your question though, that this is not the case and that $R\delta I$ is in fact comparable with $\delta V$.
One approximate way to proceed would be to do the fit in two steps. After the first step, increase the errors on the y-axis quantity so that they are $((\delta V)^2 + (R\delta I)^2)^{1/2}$. Re-run the fit and this gives you a new value of $R$ and its uncertainty.
Section 15.3 of Numerical Recipes by Press et al. gives one perhaps less approximate route through generalisation of chi-squared fitting.
Cameron Reed (2010) describes a spreadsheet implementation of an iterative least squares method attributable to York (1966).
Gull (1989) describes a Bayesian approach to the problem, where the output would be the posterior probability distribution for the resistance.
Finally, I offer the following: each pair of $V$ and $I$ values gives an independent estimate of $R$ with an uncertainty that could be calculated using standard error propagation formulae. You could then find the weighted mean and uncertainty in the weighted mean of these estimates. Reference.
