We introduce the method of desingularization of multi-variable multiple zeta-functions (of the generalized Euler-Zagier type), under the motivation of finding a suitable rigorous meaning of the values of multiple zeta-functions at non-positive integer points. We reveal that multiple zeta-functions (which are known to be meromorphic in the whole space with infinitely many singular hyperplanes) turn out to be entire on the whole space after taking the desingularization. The desingularized function is given by a suitable finite "linear" combination of multiple zeta-functions with some arguments shifted. It is shown that specific combinations of Bernoulli numbers attain the special values at their non-positive integers of the desingularized ones. We also discuss twisted multiple zeta-functions, which can be continued to entire functions, and their special values at non-positive integer points can be explicitly calculated.