Multiband processing in one pass
The most direct way to make these calculations is to design a set of filters to your specifications, then read the data out of the WAV files and run it through the filters. In Python/SciPy, you can use Python's wave module to read and write WAV files, SciPy's firwin to design a windowed sinc filter like SoX's sinc effect, and SciPy's lfilter to apply each filter. You will have to design each filter separately and call lfilter for each filter.

A simple way to approximate the energy in several bands is to first calculate the Discrete Fourier Transform (DFT) of consecutive frames, and sum up the energy in the bins corresponding to your bands of interest. The DFT approximates a signal with a sum of harmonically related sinusoids. The advantage of this technique is that it's easy to sum up several arbitrarily spaced bands. A big disadvantage is you have to apply a window to control spectral leakage. This makes the energy calculation an approximation. You can use NumPy's rfft function to calculate the Fourier Transforms.

If your bank of filters is related, for instance all the same bandwidth, or centered at harmonically or logarithmically related frequencies, there are some clever techniques that fall under the concept of filterbanks. This is more specialized, and I don't know of much pre-packaged code you could use.