  The Poland--Scheraga model, introduced in the 1970's, is a reference model to
describe the denaturation transition of DNA. More recently, it has been
generalized in order to allow for asymmetry in the strands lengths and in the
formation of loops: the mathematical representation is based on a bivariate
renewal process, that describes the pairs of bases that bond together. In this
paper, we consider a disordered version of the model, in which the two strands
interact via a potential $\beta V(\hat\omega_i,\bar\omega_j)+h$ when the $i$-th
monomer of the first strand and the $j$-th monomer of the second strand meet.
Here, $h\in\mathbb R$ is a homogeneous pinning parameter,
$(\hat\omega_i)_{i\geq 1}$ and $(\bar\omega_j)_{j\geq 1}$ are two sequences of
i.i.d.~random variables attached to each DNA strand, $V(\cdot,\cdot)$ is an
interaction function and $\beta>0$ is the disorder intensity. Our main result
finds some condition on the underlying bivariate renewal so that, if one takes
$\beta,h\downarrow0$ at some appropriate (explicit) rate as the length of the
strands go to infinity, the partition function of the model admits a
non-trivial, i.e. disordered, scaling limit. This is known as an
\textit{intermediate disorder} regime and is linked to the question of disorder
relevance for the denaturation transition. Interestingly and surprisingly, the
rate at which one has to take $\beta\downarrow0$ depends on the interaction
function $V(\cdot,\cdot)$ and on the distribution of $(\hat\omega_i)_{i\geq
1}$, $(\bar\omega_j)_{j\geq 1}$. On the other hand, the intermediate disorder
limit of the partition function, when it exists, is universal: it is expressed
as a chaos expansion of iterated integrals against a Gaussian
process~$\mathcal{M}$, which arises as the scaling limit of the field
$(e^{\beta V(\hat\omega_i,\bar\omega_j)})_{i,j\geq 0}$ and exhibits strong
correlations on lines and columns.
