The motion of two rigid spherical particles in an arbitrary configuration in an infinite viscous fluid at low Reynolds numbers is considered. The fluid is allowed to slip at the surfaces of the spheres and the particles may differ in radius. The resistance and mobility functions that completely characterize the linear relations between the forces and torques and the translational and rotational velocities of the particles are analytically calculated in the quasi-steady limit using a method of twin multipole expansions. For each function, an expression of power series in r(-1) is obtained, where r is the distance between the particle centers. The agreement between these expressions and the relevant results in the literature is quite good. Based on a microscopic model, the analytical results for two-sphere hydrodynamic interactions are used to find the effect of the volume fraction of particles of each type on the average settling velocities in a bounded suspension of slip spheres. Our results, presented in simple closed forms, agree very well with the existing solutions for the limiting cases of no slip and perfect slip at the particles surfaces. In general, the particle-interaction effects are found to be more significant when the slip coefficients at the particle surfaces become smaller. Also, the influence of the interactions on the smaller particles is stronger than on the larger ones.