Using CAD software, I have been able to make a construction in which four arbitrary lines in space (infinitely long) have two lines that intersect with all four. There are only ever two solutions that the software finds and degrees of freedom analysis vaguely suggests this is the expected result.
So we have four lines in space:
$$ \vec I = \vec I_0 + \alpha \vec I_1$$ $$\vec J = \vec J_0 + \beta \vec J_1$$ $$ \vec K = \vec K_0 + \gamma \vec K_1$$ $$ \vec L = \vec L_0 + \kappa \vec L_1$$ where $\vec \cdot$ denotes a vector. Two more lines make the intersections: $$ \vec X = \vec X + \eta \vec X_1$$ $$ \vec Y = \vec Y_0 + \zeta \vec Y_1$$ There will be a unique set of $\alpha$, $\beta$, $\gamma$, $\kappa$ and $\eta$ where the line $X = I$, $X = J$, $X = K$ and $X = L$
There will be another unique set of $\alpha$, $\beta$, $\gamma$, $\kappa$ and $\zeta$ where the line $Y = I$, $Y = J$, $Y = K$ and $Y = L$.
I am supposing the six scalars will each have a parabolic function of a common parameter, say $f(t)$, where the solutions to $f(t) = 0$ are the unique values.
This also creates a challenge of identifying which solution corresponds to which intersection.
I can vaguely see a way through the mathematics but the complexity is getting beyond me. Can anyone help?

To find the line, you need to set up the linear system of equations
$$\vec X_0 + \eta_1 \vec X_1 = \vec I_0 + \alpha \vec I_1$$ $$\vec X_0 + \eta_2 \vec X_1 = \vec J_0 + \beta \vec J_1$$ $$\vec X_0 + \eta_3 \vec X_1 = \vec K_0 + \gamma \vec K_1$$ $$\vec X_0 + \eta_4 \vec X_1 = \vec L_0 + \kappa \vec L_1$$ which imposes that the line $\vec X$ has a common point with all the other lines. This is a nonlinear system of $12$ equations (each equality has to hold for $3$ dimensions) in $13$ unknowns ($3$ values for $\vec X_0$, $3$ values for $\vec X_1$, and $7$ parameters $\eta_2,\ldots,\eta_4,\alpha,\beta,\gamma,\kappa$; there is a scaling factor in $\vec X_1$ that is solved by setting $\eta_1 = 1$). This underdetermined system may not have a solution, so i believe that not all sets of 4 lines have a line intersecting them all.
Edit:
You can further reduce the system by solving for $\vec X_0$ in the first equation, $$ \vec X_0 = \vec I_0 + \alpha \vec I_1 - \vec X_1,$$ then replacing this in the remaining equations $$\vec I_0 + \alpha \vec I_1 - \vec X_1 + \eta_2 \vec X_1 = \vec J_0 + \beta \vec J_1$$ $$ \vec I_0 + \alpha \vec I_1 - \vec X_1+ \eta_3 \vec X_1 = \vec K_0 + \gamma \vec K_1$$ $$\vec I_0 + \alpha \vec I_1 - \vec X_1+ \eta_4 \vec X_1 = \vec L_0 + \kappa \vec L_1$$
Further still:
Could this be useful? Six remaining unknowns with ? equations.
$$\frac{\vec J_0 + \beta \vec J_1}{\eta_2} = \frac{\vec K_0 + \gamma \vec K_1}{\eta_3} = \frac{\vec L_0 + \kappa \vec L_1}{\eta_4}$$
Or perhaps rearranging each and using each dimension:
$$\alpha \vec I_1 - \beta \vec J_1 + (1 - \eta_2) \vec X_1 = \vec J_0 - \vec I_0$$
$$ \begin{bmatrix} {I_1}_x & {J_1}_x & {X_1}_x \\ {I_1}_y & {J_1}_y & {X_1}_y \\ {I_1}_z & {J_1}_z & {X_1}_z \\ \end{bmatrix} \begin{Bmatrix} \alpha \\ -\beta \\ 1 - \eta_2 \\ \end{Bmatrix} = \begin{Bmatrix} {J_0}_x - {I_0}_x \\ {J_0}_y - {I_0}_y \\ {J_0}_z - {I_0}_z \\ \end{Bmatrix} $$
...
However this is confusing as it suggests there is a single solution for each value which we know not to be the case. Alternately it suggests there are three solutions for each value instead of the expected two.