Let P-beta((V))(N-I) be the probability that a N x N beta-ensemble of random matrices with confining potential V (x) has N-I eigenvalues inside an interval I = [a, b] on the real line. We introduce a general formalism, based on the Coulomb gas technique and the resolvent method, to compute analytically P-beta((V))(N-I) for large N. We show that this probability scales for large N as P-beta((V))(N-I) approximate to exp [-beta N-2 psi((V))(N-I/N)], where beta is the Dyson index of the ensemble. The rate function psi((V))(k(I)), independent of beta, is computed in terms of single integrals that can be easily evaluated numerically. The general formalism is then applied to the classical beta-Gaussian (I = [-L, L]), beta-Wishart (I = [1, L]), and beta-Cauchy (I = [-L, L]) ensembles. Expanding the rate function around its minimum, we find that generically the number variance var(N-I) exhibits a nonmonotonic behavior as a function of the size of the interval, with a maximum that can be precisely characterized. These analytical results, corroborated by numerical simulations, provide the full counting statistics of many systems where random matrix models apply. In particular, we present results for the full counting statistics of zero-temperature one-dimensional spinless fermions in a harmonic trap.