We have $$ I=\int_{0}^{1}\frac{\log\left(x\right)\log^{3}\left(1-x^{2}\right)}{1-x}dx=\int_{0}^{1}\frac{\log\left(x\right)\log^{3}\left(1-x^{2}\right)}{1-x^{2}}dx+\int_{0}^{1}\frac{x\log\left(x\right)\log^{3}\left(1-x^{2}\right)}{1-x^{2}}dx
$$ and so if we put $x=\sqrt{y}$ we get $$I=\frac{1}{4}\int_{0}^{1}\frac{y^{-1/2}\log\left(y\right)\log^{3}\left(1-y\right)}{1-y}dy+\frac{1}{4}\int_{0}^{1}\frac{\log\left(y\right)\log^{3}\left(1-y\right)}{1-y}dy$$ and recalling the definition of beta function $$B\left(a,b\right)=\int_{0}^{1}x^{a-1}\left(1-x\right)^{b-1}dx
$$ we have $$\frac{\partial^{h+k}B}{\partial a^{h}\partial b^{k}}\left(a,b\right)=\int_{0}^{1}x^{a-1}\log^{h}\left(x\right)\left(1-x\right)^{b-1}\log^{k}\left(1-x\right)dx
$$ hence $$I=\frac{1}{4}\frac{\partial^{4}B}{\partial a\partial b^{3}}\left(\frac{1}{2},0^{+}\right)+\frac{1}{4}\frac{\partial^{4}B}{\partial a\partial b^{3}}\left(1,0^{+}\right).$$ For the computation of the limit, we can use the asymptotic $\Gamma(x)=\frac{1}{x}+O(1)$ when $x\rightarrow 0$ and the relations between the polygamma terms and zeta.