Why should the modes of the linearized metric perturbation be "wavefunctions" of gravitons (in the Randall-Sundrum model)?

This is a slight abuse of terminology, related to talking about 'second quantization.' The word 'wave function' in this case really refers to the 'one particle wave function,' which happens to correspond to the solutions of the (linear) classical equations of motion. It does not refer to the 'wave functional,' ie the Schrodigner representation of the full quantum field theory which is of course a rather complicated object.

Let's work with a scalar field $\phi(x)$ living on some fixed geometry described by $g_{\mu\nu}$ (it's not hard to generalize this discussion to include spin). We will discuss the free theory, which can be thought of as the first step to setting up a perturbative treatment of an interacting theory.

For a free scalar obeying the operator equation \begin{equation} \square \phi = \frac{1}{\sqrt{-g}}\partial_\mu\left(g^{\mu\nu}\partial_\nu \phi\right) = 0, \end{equation} you can expand in mode functions \begin{equation} \phi(x) = \sum_{n} a_n u_n(x) + a_n^\dagger u_n^{\star}(x), \end{equation} where $u_n(x)$ are the mode functions associated with the wave operator, appropriately normalized. As an example, for Minkowski spacetime, we can replace the sum over $n$ with an integral $\sum_n \rightarrow \int d^3 p / (2\pi)^3$ and the appropriately normalized mode functions are $u_n(x)\rightarrow u_p(x) = e^{i p\cdot x}/\sqrt{2 E_p}$.

You can then define a position eigenstate (at least formally) by \begin{equation} |x\rangle = \phi(x) |0\rangle = \sum_n u_n^\star a_n^\dagger |0\rangle = \sum_n u_n^\star |n\rangle. \end{equation}

In an interacting relativistic quantum field theory the notion of a one particle position space eigenstate isn't obviously well defined because you can't localize particles arbitrarily precisely without having a probability to create particle/anti-particle pairs. But, at the level of the free theory it is ok and if we work perturbatively and squint we can imagine that the notion of position eigenstate will be approximately correct so long as we keep in mind that this is an idealization and we don't try to compute physical quantities that rely on a particle being localized to within its compton wavelength (or you could be more rigorous and define everything in terms of wavepackets).

Carrying on ignoring that subtlety, we can define the wave function of the one particle state $|n\rangle$ by its projection onto the position eigenstate \begin{equation} \psi_n(x) \equiv \langle x | n \rangle = u_n(x), \end{equation} where $|n\rangle = a_n^\dagger |0 \rangle$ (note that $n$ here is not an occupation number, it is a label for the modes). Intuitively a particle is a quanta of the mode $n$, roughly corresponding to a small vibration in a given classical mode of the field which gives some sense of where you might be likely to find that particle.

Thus, there is a direct relationship between the one particle position eigenstate (the 'wavefunction'), and the modes of the classical equations of motion.

What you typically want in extra dimensional models like Randall-Sundrum is for the one particle eigenstates to be localized in the extra dimensions. Let's say that spacetime is $4+1$ dimensional, and our universe lives on a brane with 3 spatial dimensions living in this larger space. The idea is that if an observer on the brane creates a graviton, that graviton should 'look like' a graviton appropriate for 3+1 dimensions which in turn will guarantee that scattering amplitudes computed by an observer on the brane will 'look like' scattering amplitudes for an observer in a 3+1 dimensional spacetime. (The exact definition of 'looks like' depends on your setup and on the experimental constraints on extra dimensions). If this condition is not satisfied, the model would be ruled out, for example by tests of the inverse square law. We can satisfy the condition if the graviton wavefunction (=the classical mode associated with the graviton wave operator = the 1 graviton state projected into the position basis) is sufficiently localized on the brane. In that case, acting on the (4+1) vacuum with a one particle graviton creation operator creates a state that 'looks like' the state you would get in 3+1 dimensions by acting with a graviton creation operator.

In physics, it is not out of the ordinary to quantize a composite system in stages: First quantize some of its components and then proceed to the rest. We do this all the time, sometimes unconsciously. Sometimes it is just a way of thinking. Taking for example, the simplest case of the Kaluza-Klein expansion when the internal space is $S^1$: $$g_{\mu\nu}(x, \theta) = \sum_{n\in \mathbb{Z}} g_{\mu\nu}^{(n)} e^{in\theta}$$ This expression is on one hand an expansion in terms of square integrable functions on the circle. The exponential functions are just an orthonormal basis of wave-functions on the circle.
In the Kaluza-Klein theory, the fields $$g_{\mu\nu}^{(n)} $$ are charged massive gravitons with $n$ units of charge. (The charge appears as the coupling constant to the U(1) field $A_{\mu})$. We may say that these gravitons spin in the fifth dimension and that is how they get their charge, or equivalently they have an internal wave function $e^{in\theta}$

This description is not limited to Kaluza Klein theories. A more realistic example would be the wave function of a spinning particle which we usually write as: $$\Psi = \psi(x) \otimes \psi_{s, m_s}$$ Where $x\in \mathcal{M}$ the configuration manifold. We usually think of the spin wave function $\psi_{s, m_s}$ as a spinor with a finite number of components. But, we can also describe the spin wave function as a section of a line bundle on the sphere $S^2$. Thus in this case we can think of the total wave function (locally) as a function $\Psi(x, \theta, \phi) $of the $\mathcal{M} \otimes S^2$ (more generally an $S^2$ bundle over $\mathcal{M}$). And instead of saying that the wave function describes a spinning particle we can view it as a wave function of a scalar particle moving on an internal space $S^2$.

The whole subject of expressing mechanical system in this fashion is known as sub-Riemannian geometry.