You will simply have to find an affine mapping $A$ from the reference tetrahedron $K$ to your arbitrary tetrahedron $T$.
If $T$ is given as the convex hull $T = \textrm{conv}\{n_i : 1\leq i \leq 4\}$ of your nodes $n_i$, one mapping could look like this:
$e_1 \mapsto n_1$, $e_2 \mapsto n_2$, $e_3 \mapsto n_3$ and $\textbf{0} \mapsto n_4$.
Or in closed form: $$\mathbf{x}_{T,i} = A(\mathbf{x}_i),$$
with $A:x \mapsto Mx + n_4,$ and matrix $M:=(\begin{matrix} n_1-n_4,\,n_2-n_4,\,n_3-n_4\end{matrix})$.
The quadrature weights will need to be scaled by the volume of the tetrahedron ($\textrm{det}(M)$/6) divided by the volume of your current tetrahedron (1/6), which equals $$\mathbf{w}_T = \textrm{det}(M)\mathbf{w}.$$
The resulting quadrature rule will usually be different depending on your choice of the ordering of nodes $n_i$.
To know why Gauss quadrature works, you should look at the proof.
As for how to do it, you need to do Gram-Schmidt on the standard polynomial basis to get a degree three orthogonal polynomial. The roots of this polynomial will constitute the nodes. Then, test your rule on the polynomials $1,x$, and $x^2$ to get a linear system of equations to find the weights. The proof explains why this is what you want to do, and this is the algorithm that you always follow for such a problem.
EDIT: How it's done-
First, we apply Gram-Schmidt to the standard basis $v_1=1,v_2=x,v_3=x^2,v_4=x^3$ with respect to the inner product $\langle f,g\rangle =\int\limits_{-1}^1 f(x)g(x)\, dx$. We go up to a cubic because the roots will be the nodes, and we want three nodes. In particular, we want to compute
\begin{align*}
u_1&=v_1\\
v_2&=u_2-\frac{\langle u_2,v_1\rangle}{\langle v_1,v_1\rangle}v_1\\
v_3&=u_3-\frac{\langle u_3,v_1\rangle}{\langle v_1,v_1\rangle}v_1-\frac{\langle u_3,v_2\rangle}{\langle v_2,v_2\rangle}v_2\\
v_4&=u_4-\frac{\langle u_4,v_1\rangle}{\langle v_1,v_1\rangle}v_1-\frac{\langle u_4,v_2\rangle}{\langle v_2,v_2\rangle}v_2-\frac{\langle u_4,v_3\rangle}{\langle v_3,v_3\rangle}v_3
\end{align*} You can find that \begin{align*}
u_1&=1\\
u_2&=x\\
u_3&=x^2-\frac{1}{3}\\
u_4&=x^3-\frac{3}{5}x
\end{align*}
The roots of $u_4$ are $x_{1,2,3}=0,\pm\sqrt{\frac{3}{5}},$ so these are the nodes.
So, our rule is $$\int\limits_{-1}^1f(x)\, dx=c_1f(0)+c_2f\left(\sqrt{\frac{3}{5}}\right)+c_3f\left(-\sqrt{\frac{3}{5}}\right),$$ and all that we have to do is find $c_1,c_2,c_3.$ We impose exactness on $f(x)=1,x,x^2.$ That is, we require that the rule is an equality for these $f$. If we compute what this means, we get that
\begin{align*}
\text{ if }f(x)=1&:\ \ \int\limits_{-1}^1 \, dx=2=c_1+c_2+c_3\\
\text{ if }f(x)=x&:\ \ \int\limits_{-1}^1x \, dx=0=c_2\sqrt{\frac{3}{5}}-c_3\sqrt{\frac{3}{5}}\\
\text{ if }f(x)=x^2&:\ \ \int\limits_{-1}^1x^2 \, dx=\frac{2}{3}=c_2\frac{3}{5}+c_3\frac{3}{5}
\end{align*}
This is a system of linear equations, and you can solve it and find that
\begin{align*}
c_1&=\frac{8}{9}\\
c_2&=\frac{5}{9}\\
c_3&=\frac{5}{9}\\
\end{align*}
This is consistent with the rule that you were looking for.
You can find a proof of the general method of Gaussian quadrature here: https://www.johndcook.com/OrthogonalPolynomials.pdf
Best Answer
Let's look first at $\int_{-1}^{1} f(t) \ dt$.
We want the formula to evaluate $\int_{-1}^{1} dt, \int_{-1}^{1} t \ dt, \int_{-1}^{1} t^{2} \ dt, \int_{-1}^{1} t^{3} \ dt, \int_{-1}^{1} t^{4} \ dt, \int_{-1}^{1} t^{5} \ dt, \int_{-1}^{1} t^{6} \ dt $ and $\int_{-1}^{1} t^{7} \ dt$ exactly.
That leads to the following system of equations:
$$ 2 = \omega_{0} + \omega_{1}$$
$$ 0 = \omega_{0} x_{0} + \omega_{1} x_{1}+ \omega_{2}+\omega_{3}$$
$$ \frac{2}{3} = \omega_{0} x_{0}^{2} + \omega_{1} x_{1}^{2} + 2 \omega_{2} x_{2} + 2 \omega_{3} x_{3}$$
$$ 0 = \omega_{0} x_{0}^{3} + \omega_{1} x_{1}^{3}+ 3 \omega_{2} x_{2}^{2} + 3 \omega_{3} x_{3}^{2} $$
$$ \frac{2}{5} = \omega_{0} x_{0}^{4} + \omega_{1} x_{1}^{4}+ 4 \omega_{2} x_{2}^{3} + 4 \omega_{3} x_{3}^{3} $$
$$0 = \omega_{0} x_{0}^{5} + \omega_{1} x_{1}^{5} + 5 \omega_{2} x_{2}^{4} + 5 \omega_{3} x_{3}^{4} $$
$$\frac{2}{7} = \omega_{0} x_{0}^{6} + \omega_{1} x_{1}^{6} + 6 \omega_{2} x_{2}^{5} + 6\omega_{3} x_{3}^{5} $$
$$ 0 = \omega_{0} x_{0}^{7} + \omega_{1} x_{1}^{7} + 7 \omega_{2} x_{2}^{6} + 7\omega_{3} x_{3}^{6}$$
Solve the system using a numerical solver.
Then use the the fact that $$\int_{a}^{b} f(x) \ dx = \frac{b-a}{2} \int_{-1}^{1} f \Big(a+(1+t)\frac{b-a}{2} \Big) \ dt$$