$\int_0^{\infty} \frac{x^{k_1}}{\left(x^{n_1} + a_1 \right)^{m_1}} \cdot \frac{x^{k_2}}{\left(x^{n_2} + a_2 \right)^{m_2}}\:dx$ Recently I was able to find a result to a common definite integral:
\begin{equation}
 J(n_1, k_1, m_1) = \int_0^{\infty} \frac{x^{k_1}}{\left(x^{n_1} + a_1 \right)^{m_1}}\:dx = \frac{a^{\frac{k_1 + 1}{n_1} - m_1}}{n_1}\Gamma\left(m_1 - \frac{k_1 + 1}{n_1}\right)\Gamma\left(\frac{k_1 + 1}{n_1}\right)
\end{equation}
The advantage with this integral is that if you can introduce a free parameter into your integral (ala Feynman's Trick) then under a given transformation using that parameter (i.e. derivatives, Laplace Transform, Fourier Transforms) that you can isolate out the parameter. This makes taking the inverse transform quite easy in most cases. 
I'm hoping to expand this method to cater to more difficult integrals where multiple parameters are introduced. In doing so, I've come this very similar but much more complicated integral.  
\begin{equation}
 H\left(a_1, a_2,n_1,n_2,k_1,k_2,m_1,m_2\right)= \int_0^{\infty} \frac{x^{k_1}}{\left(x^{n_1} + a_1 \right)^{m_1}} \cdot \frac{x^{k_2}}{\left(x^{n_2} + a_2 \right)^{m_2}}\:dx
\end{equation}
Where $a_1, a_2,n_1,n_2,k_1,k_2,m_1,m_2 \in \mathbb{R}^{+}$ 
Unfortunately I'm stuck with this one and am interested to find out if anyone has encountered this form before and if so, if they know a method to find a solution expressed in terms of elementary or non-elementary functions. 
Any starting points would be greatly appreciated. 
Also, if the solution can be represented with certain restrictions on the parameters, please post up. The only absolute conditions are $a_{1},a_{2},m_{1},m_{2} \gt 0$. Keen on all solutions restricted or not.
Edit - Just thinking now that the one thing that can be done to simplify the integral is to let $u = x^{n_1}$ or $u = x^{n_2}$. Here I will use the later to render the integral as 
\begin{align}
 &H\left(a_1, a_2,n_1,n_2,k_1,k_2,m_1,m_2\right)= \int_0^{\infty} \frac{\left( u^{\frac{1}{n_2}}\right)^{k_1}}{\left(\left( u^{\frac{1}{n_2}}\right)^{n_1} + a_1 \right)^{m_1}} \cdot \frac{\left( u^{\frac{1}{n_2}}\right)^{k_2}}{\left(u + a_2 \right)^{m_2}}\frac{1}{n_2}u^{\frac{1 - n_2}{n_2}}\:du\\
\quad& = \frac{1}{n_2}\int_0^{\infty} \frac{u^{\frac{k_1 + k_2 + 1 - n_2}{n_2}}}{\left( u^{\frac{n_1}{n_2}} + a_1 \right)^{m_1}} \frac{1}{\left(u + a_2\right)^{m_1}}\:du = \frac{1}{n_2}\int_0^{\infty} \frac{u^{k_3}}{\left(u^{n_3} + a_1\right)^{m_1}} \frac{1}{\left(u + a_2\right)^{m_2}}\:du
\end{align}
 A: This is an integral of two Meijer G-functions which gives a Fox H-function. Let $\nu_1 = 1/n_1$, $\,\nu_2 = 1/n_2$, $\,\sigma = k_1 + k_2 + 1$. The Mellin transform is
$$\mathcal M_{x \to p}[(x + a)^{-m}] =
\frac {\Gamma(p) \Gamma(m - p)} {\Gamma(m)} a^{p -m}, \\
\mathcal M[(x^{n_1} + a_1)^{-m_1}] =
\nu_1 \mathcal M[(x + a_1)^{-m_1}](\nu_1 p), \\
\mathcal M[x^{-\sigma} (x^{-n_2} + a_2)^{-m_2}] =
\nu_2 \mathcal M[(x + a_2)^{-m_2}](\nu_2 (\sigma - p)).$$
The convolution is mapped to the product:
$$\mathcal M[f * g] =
\mathcal M {\left[\int_0^\infty f(x)
 \,g {\left( \frac \omega x \right)} \,\frac {dx} x \right]} =
\mathcal M[f] \mathcal M[g].$$
The product is again a rational combination of linear gamma functions times $z^p$, the inverse transform of which is, by definition, the H-function. Evaluating it at $1$ ($\omega = 1$ above) gives
$$\int_0^\infty
 \frac {x^{k_1}} {(x^{n_1} + a_1)^{m_1}}
 \frac {x^{k_2}} {(x^{n_2} + a_2)^{m_2}} dx = 
\frac
  {\nu_1 \nu_2 a_1^{-m_1} a_2^{\nu_2 \sigma - m_2}}
  {\Gamma(m_1) \Gamma(m_2)} \times \\
 \mathcal M^{-1}[
  \Gamma(m_1 - \nu_1 p) \Gamma(\nu_2 \sigma - \nu_2 p)
  \Gamma(\nu_1 p) \Gamma(m_2 - \nu_2 \sigma + \nu_2 p)
  (a_1^{-\nu_1} a_2^{\nu_2})^{-p}](1) = \\
\frac
  {\nu_1 \nu_2 a_1^{-m_1} a_2^{\nu_2 \sigma - m_2}}
  {\Gamma(m_1) \Gamma(m_2)}
 H_{2, 2}^{2, 2} {\left( a_1^{-\nu_1} a_2^{\nu_2} \middle|
  {(1 - m_1, \nu_1), (1 - \nu_2 \sigma, \nu_2) \atop
  (0, \nu_1), (m_2 - \nu_2 \sigma, \nu_2)} \right)}.$$
The integral exists if $\sigma < n_1 m_1 + n_2 m_2$.
A: To long for the comment.
Maybe it's too trite, but 
In the case $n_1=n_2$ this integral can be rewrite through the hypergeometric function $_2F_1$.
I do not want to pay attention to the coefficient. I this case the integral can be rewrite in the following form
$$
I= \alpha\int_0^{\infty} \frac{u^{k_3}}{\left(u + 1\right)^{m_1}} \frac{1}{\left(u + c\right)^{m_2}}\:du
$$
like it was done by DavidG.
After that one can do the change of variables $y=u/(1+u)$
$$
I=\frac{\alpha}{c^{m_2}}\int_0^{1} \frac{y^{k_3}}{\left(1 - y\right)^{k_3+2-m_1-m_2} (1-y(1-1/c))^{m_2}}\:du=\frac{\alpha}{c^{m_2}}Г(k_3+1) Г(m_1+m_2-k_3-1)  _2F_1(m_2,k_3+1,m_1+m_2,1-1/c)
$$
I am not sure about the existance of the answer in the general case.
A: HINT
Let us start from the OP result
\begin{align}
&H_1\left(a_1, a_2, m_1, m_2, n_3, k_3\right) 
= \frac{1}{n_2}\int\limits_0^{\infty} \frac{u^{k_3}}{\left(u^{n_3} + a_1\right)^{m_1}} \frac{1}{\left(u + a_2\right)^{m_2}}\:du.\tag1
\end{align}
Taking in account, that identity

\begin{align} &\frac{x^n}{\left(x^n + a\right)^m}
 =\frac{x^n+a-a}{\left(u^n + a\right)^m}
 =\frac{1}{\left(x^n + a\right)^{m-1}}
 -\frac{a}{\left(x^n + a\right)^m}\tag2 \end{align}

allows to manage the degree in the numerator both decreasing and increasing,
assume WLOG $\mathbf{k_3<m_1,\ k<m_2}.$
Also, due to the scaling substitution, let WLOG $\mathbf{a_1=1}.$
Under these conditions, enough to consider the integral
\begin{align}
&I(a,k,n,p,q) 
= \int\limits_0^{\infty} \frac{u^k}{\left(u^n + 1\right)^p} \frac{1}{\left(u + a\right)^{q}}\:du,\quad k<p,\quad k<q, \quad a>1.\tag3
\end{align}
Is known the binomial decomposition

\begin{align}
 &(1-z)^{-d-1}=\sum\limits_{j=d}^\infty\binom{j}{d}z^{j-d} 
 = \sum\limits_{j=d}^\infty\dfrac{\Gamma(j+1)}{\Gamma(d+1)\Gamma(j-d+1)}z^{j-d},\quad|z|<1.\tag4
 \end{align}

If $|u|<1$ then from $(4)$ should Maclaurin series
\begin{align}
&(1+u^n)^{-q} = \sum\limits_{j=q-1}^\infty(-1)^{j-q+1}\dfrac{\Gamma(j+1)}{\Gamma(q)\Gamma(j-q+2)}u^{n(j-q+1)}.\tag5
\end{align}
If $|u|>1$ then from $(4)$ can be obtained Laurent series
\begin{align}
&(1+u^n)^{-q} = u^{-nq}(1+u^n)^{-q} = \sum\limits_{j=q-1}^\infty(-1)^{j-q+1} \dfrac{\Gamma(j+1)}{\Gamma(q)\Gamma(j-q+2)}u^{-n(j-q+1)}.\tag6
\end{align}
At the same time, obtained integrals on the both intervals can be expressed via hypergeometric function

$$\int_0^1\dfrac{x^k}{(x+a)^m}\:dx=\dfrac{a^{-m}}{k+1}\,{_2F_1\left(k+1,m,k+2,-\dfrac1a\right)}\tag7$$

(see also Wolfram Alpha),

$$\int_1^\infty\dfrac{x^k}{(x+a)^m}\:dx=\dfrac{_2F_1(m,-k+m-1,m-k,-a)}{-k+m-1}\tag8$$

(see also Wolfram Alpha).
Splitting the issue integral to the intervals $(0,1)$ and $(1,\infty)$ and applying solutions $(5)-(8),$ one can present the result $I(a,k,n,p,q)$ as the sum of the pair of series, using gamma function and hypergeometric functions.
Thanks to the other members of the discussion for their interesting ideas.
A: This is not a solution. It is an approach that shows where some of the potential roadblocks lay. Of course, this is not to say they are not insurmountable, but $\ldots$
The approach to be use will make use of the so-called Schwinger parametrisation that makes use of the well-known observation of
$$\frac{1}{\beta^p (x)} = \frac{1}{\Gamma (p)} \int_0^\infty u^{p - 1} e^{-u \beta (x)} \, du,$$
for suitable conditions on the function $\beta (x)$.
Let us use this parametrisation on your first integral $J(a_1, n_1, m_1, k_1)$. From Schwinger parametrisation we observe that
$$\frac{1}{(x^{n_1} + a_1)^{m_1}} = \frac{1}{\Gamma (m_1)} \int_0^\infty u^{m_1 - 1} e^{-u (x^{n_1} + a_1)} \, du,$$
(here $\beta (x) = x^{n_1} + a_1$ and $p = m_1$). Your integral for $J$ then becomes
$$J(a_1, n_1, m_1, k_1) = \frac{1}{\Gamma (m_1)} \int_0^\infty u^{m_1 - 1} e^{-u a_1} \, du \int_0^\infty x^{k_1} e^{-u x^{n_1}} \, dx,$$
where a change in the order of integration has been made. Now if we let $t = u x^{n_1}$ we have
\begin{align}
J(a_1, n_1, m_1, k_1) &= \frac{1}{\Gamma (m_1)} \int_0^\infty u^{m_1 - 1} e^{-u a_1} \, du \frac{1}{n_1} u^{-\frac{(k_1 + 1)}{n_1}} \int_0^\infty t^{\frac{k_1 + 1}{n_1} - 1} e^{-t} \, dt\\
&= \frac{1}{n_1 \Gamma (m_1)} \Gamma \left (\frac{k_1 + 1}{n_1} \right ) \int_0^\infty u^{\frac{m_1 n_1 - k_1 - 1}{n_1} - 1} e^{-u a_1} \, du.
\end{align}
Enforcing a substitution of $u \mapsto u/a_1$ leads to
\begin{align}
J(a_1, n_1, m_1, k_1) &= \frac{a_1^{\frac{k_1 +1}{n_1} - m_1}}{n_1 \Gamma (m_1)} \Gamma \left (\frac{k_1 + 1}{n_1} \right ) \int_0^\infty  u^{\frac{n_1 m_1 - k_1 - 1}{n_1} - 1} e^{-u} \, du\\[2ex]
&= \frac{a_1^{\frac{k_1 +1}{n_1} - m_1}}{n_1 \Gamma (m_1)} \Gamma \left (m_1 - \frac{k_1 + 1}{n_1} \right ) \Gamma \left (\frac{k_1 + 1}{n_1} \right ).
\end{align}

Now attempting such an approach on the integral $H$. Here the following two Schwinger parametrisations are used:
$$\frac{1}{(x^{n_1} + a_1)^{m_1}} = \frac{1}{\Gamma (m_1)} \int_0^\infty y^{m_1 - 1} e^{-y (x^{n_1} + a_1)} \, dy,$$
and
$$\frac{1}{(x^{n_2} + a_2)^{m_2}} = \frac{1}{\Gamma (m_2)} \int_0^\infty z^{m_2 - 1} e^{-z (x^{n_2} + a_2)} \, dz.$$
So the integral for $H$, after a change in the order of integration has been made, becomes
$$H = \frac{1}{\Gamma (m_1) \Gamma (m_2)} \int_0^\infty y^{m_1 - 1} e^{-y a_1} \, dy \int_0^\infty z^{m_2 - 1} e^{-z a_2} \, dz \int_0^\infty x^{k_1 + k_2} e^{-(y x^{n_1} + z x^{n_2})} \, dx.$$
The roadblock is what is to be done with that inner $x$-integral if $n_1 \neq n_2$?  
