A lattice QCD calculation of s-wave hadron scattering lengths in the channels π-π, π-N, K-N, K¯-N, and N-N is carried out in the quenched approximation at β=6/g2=5.7. A variant of the method of wall source is developed for this purpose, which reduces the computer time by a factor L3 on an L3×T lattice compared to the conventional point source method and avoids the Fierz mixing problem. A version of the method in which gauge configurations are not fixed to any gauge can be extended to calculate disconnected quark loop contributions in hadron two- and three-point functions. An analytical estimate of statistial errors for this method is worked out, and the magnitude of errors without and with gauge fixing is compared for the case of π-π four-point functions calculated with the Kogut-Susskind quark action. For π-π scattering both I=0 and 2 scattering lengths are evaluated ujsing the Kogut-Susskind and Wilson quark actions on a 123×20 lattice. For the same size lattice, π-N, K-N, and K¯-N scattering lengths are calculated with the Wilson quark action. For the π-π and π-N cases simulation results are consistent with the predictions of current algebra and PCAC within one to two standard deviations up to quite heavy quark masses corresponding to mπ/mρ≊0.74, while for the K-N and K¯-N cases the agreement is within a factor of 2.For N-N scattering a phenomenological study with one-boson exchange potentials indicate that the deuteron becomes unbound if the quark mass is increased beyond 30–400f the physical value. Simulations with the Wilson action on a 204 lattice with heavy quarks with mπ/mρ≊0.74–0.95 show that the nucleon-nucleon force is attractive for both spin triplet and singlet channels, and that the scatteirng lengths are substantially larger compared to those for the π-π and π-N cases even for such heavy quarks. The problem of statistical errors, which has to be overcome toward a more realistic calculation of hadron scattering lengths, is discussed.