First, in view of Legrende's duplication formula,
$$S=\sum_{n=1}^\infty\frac{\Gamma\left(n+\frac{1}{2}\right)}{(2n+1)^4\,4^n\,n!}=2\sqrt{\pi}\sum_{n=1}^\infty\frac{\Gamma(2n)}{\Gamma(n)\,n!\,(2n+1)^4\, 16^n}
\\=-\frac{\sqrt{\pi}}{3}\int_0^1 \ln^3(x)\sum_{n=1}^{\infty}\frac{\Gamma(2n)}{\Gamma(n)\,n!}\left(\frac{x^2}{16}\right)^ndx\\
=-\sqrt{\pi}-\frac{\sqrt{\pi}}{6}\int_0^1\frac{\ln^3(x)}{\sqrt{1-x^2/4}}dx=-\sqrt{\pi}-\frac{\sqrt{\pi}}{3}\int_0^{\frac{\pi}{6}}\ln^3(2\sin x)dx$$
Claim: for $0<a\leq \frac{\pi}{2}$,
$$\int_0^a \ln^3\left(\frac{\sin x}{\sin a}\right)dx\tag{0}=\frac{4a-3\pi}{2}a^2\ln(2\sin a)-\frac{3\pi}{4}\zeta(3)+3\left(\frac{\pi}{2}-a\right)\Re\left(\frac12 \operatorname{Li}_3(e^{2ia})+\operatorname{Li}_3(1-e^{2ia})\right)+3\Im\left(\frac14\operatorname{Li}_4(e^{2ia})+\operatorname{Li}_4(1-e^{2ia})\right)
$$
Proof.
The idea is exactly identical to the proof displayed in this question. The proof is rather tedious (and obviously inefficient), and ends with a somewhat of a cancellation (implying the existence of a shortcut) , so I omit the boring algebra and outline the main ideas, which can be repeated systematically to obtain closed forms for even higher powers of logsine.
things to know:
$$\ln(2\sin x)=\ln(1-e^{2ix})+i\left(\frac{\pi}{2}-x\right) \tag{1}$$
$$\small\int\frac{\ln^3(1-x)}{x}dx=\ln^3(1-x)\ln(x)+3\ln^2(1-x)\text{Li}_2(1-x)-6\ln(1-x)\text{Li}_3(1-x)+6\text{Li}_4(1-x) \tag{2}$$
$$\int_0^a x\ln(2\sin x)dx=-\frac{a}{2}\text{Cl}_2(2a)-\frac14\Re\text{Li}_3(e^{2ia})+\frac{\zeta(3)}{4}\tag{3}$$
$$\int_0^a x^2\ln(2\sin x)dx=-\frac{a^2}{2}\text{Cl}_2(2a)-\frac{a}{2}\Re\text{Li}_3(e^{2ia})+\frac14\Im\text{Li}_4(e^{2ia})\tag{4}$$
$$\int_0^a \ln(\sin x)dx=-a\ln2-\frac12 \text{Cl}_2(2a)\tag{5}$$
$$\int_0^a \ln^2(\sin x)dx=\frac{a^3}{3}+a\ln^2 2-a\ln^2(2\sin a)-\ln(\sin a)\text{Cl}_2(2a)-\Im\text{Li}_3(1-e^{2ia})\tag{6}$$
$(1)$ is trivial, $(2)$ is not too hard to find, $(5)$ and $(6)$ are shown in the linked answer, and $(3)$&$(4)$ are easily found using $\,\,\ln(2\sin x)=-\sum_{n\geq1}\frac{\cos(2xn)}{n}$.
It is obvious that since we have $(5)$ and $(6)$, the claim $(0)$ depends on a closed form for $\displaystyle\int_0^a \ln^3(\sin x)dx$, and the latter may be evaluated in terms of $\displaystyle\int_0^a \ln^3(2\sin x)dx$.
But, with the help of $(1)$,
$$\int_0^a \ln^3(2\sin x)dx=\Re\int_0^a \ln^3(1-e^{2ix})dx+3\int_0^a \ln(2\sin x)\left(\frac{\pi}{2}-x\right)^2dx\\
=\frac12\Im\int_1^{e^{2ia}}\frac{\ln^3(1-x)}{x} dx+3\int_0^a \ln(2\sin x)\left(\frac{\pi}{2}-x\right)^2dx$$
(Same idea @RandomVariable had in this answer.)
Now we employ $(2),(3),(4),$ and $(5)$. Some expressions cancel and claim follows.$\square $
This result, together with the fact that $e^{i\pi/3}$ and $1-e^{i\pi/3}$ are conjugates, yields $\displaystyle \int_0^{\frac{\pi}{6}} \ln^3(2\sin x)dx=-\frac{\pi}{4}\zeta(3)-\frac94\Im\text{Li}_4(e^{i\pi/3})$,
and
$$S=\sqrt{\pi}\left(\frac{\pi}{12}\zeta(3)+\frac{9}{12}\Im\text{Li}_4(e^{i\pi/3})-1\right)$$
This form is equivalent to @user153012's form, as
$$\frac{2}{\sqrt{3}}\Im\text{Li}_4(e^{i\pi/3})=\sum_{n\geq 0}\frac{(-1)^n}{(3n+1)^4}+\sum_{n\geq 0}\frac{(-1)^n}{(3n+2)^4} \\=\frac{\psi^{(3)}\left(\frac13\right)}{216}-\frac{\pi^4}{81}$$
Also, as noted in the comments in the linked question, this may be used to write a closed form for a certain hypergeometric function.
This serves as a generalisation for the series, because $\displaystyle \sum_{n=1}^{\infty} \frac{\Gamma(n+1/2)}{(2n+1)^4 n!}a^{2n}=-\sqrt{\pi}\left(1+\frac1{6a}\int_0^{\sin^{-1} a}\ln^3\left(\frac{\sin x}{a}\right)dx\right)$
As an example, using closed forms for trilogarithms displayed in this post, we have
$$\int_0^{\frac{\pi}{4}}\ln^3(\sqrt{2}\sin x)dx=-\frac{\pi^3}{128}\ln2-\frac{3\pi}{8}\zeta(3)+\frac34\beta(4)+3\Im\text{Li}_4(1-i)$$
where $\beta(4)=\Im\text{Li}_4(i)$ is a value of Dirichlet's beta function.
Or equivalently,
$$\sum_{n=1}^{\infty} \frac{\Gamma\left(n+\frac12\right)}{(2n+1)^4\,2^n\,n!}=-\sqrt{\pi}-\frac{\sqrt{2\pi}}{6}\left(-\frac{\pi^3}{128}\ln2-\frac{3\pi}{8}\zeta(3)+\frac34\beta(4)+3\Im\text{Li}_4(1-i)\right)$$
Best Answer
Our goal is to see why hypergeometric functions end up in the solution to the integral $$ I=\int \frac{1}{(x+q_1)^{s_1}}\log_2\left(\frac{\left(x+q_{2}\right)^{s_{2}}}{\left(x+q_{1}\right)^{s_{1}}}\right)dx. $$ Letting $\log z$ denote the natural logarithm we have by the change of base formula $\log_2z=\log z/\log2$. Then utilizing properties of logarithms, we are able to write $I$ as $$ I=\frac{s_2}{\log 2}\int \frac{\log(x+q_2)}{(x+q_1)^{s_1}}\,\mathrm dx-\frac{s_1}{\log 2}\int\frac{\log(x+q_{1})}{(x+q_1)^{s_1}}dx. $$ The second integral is easily evaluated in terms of elementrary functions yielding $$ I=\frac{s_2}{\log 2}\underbrace{\int \frac{\log(x+q_2)}{(x+q_1)^{s_1}}\,\mathrm dx}_J-\frac{s_1}{\log 2}\frac{\left(q_1+x\right){}^{1-s_1} \left(-s_1 \log \left(q_1+x\right)+\log \left(q_1+x\right)-1\right)}{\left(s_1-1\right){}^2}. $$ All that is left is to evaluate $J$. Integrating by parts with $u=\log(x+q_2)$ and $\mathrm dv=(x+q_1)^{-s_1}\,\mathrm dx$ we find $$ J=\frac{1}{1-s_1}\frac{\log(x+q_2)}{(x+q_1)^{s_1-1}}-\frac{1}{1-s_1}\underbrace{\int \frac{x+q_1}{x+q_2}\frac{1}{(x+q_1)^{s_1}}\,\mathrm dx}_{J^\prime}. $$ Note that $J^\prime$ is the integral found in the title of this post. Evaluating $J^\prime$ is done by means of a substitution. Substituting $u=(x+q_1)/(x+q_2)$ we obtain $$ J^\prime=\int((q_2-q_1)u)^{1-s_1}(1-u)^{s_1-2}\,\mathrm du. $$ Performing a second substitution of $t=(q_2-q_1)u$ further gives $$ J^\prime=\frac{1}{q_2-q_1}\int \frac{t^{1-s_1}}{(1-zt)^{2-s_1}}\,\mathrm dt, $$ where $z=(q_2-q_1)^{-1}$. By the fundamental theorem of calculus this is equivalent to $$ J^\prime=\frac{1}{q_2-q_1}\int_0^t \frac{v^{1-s_1}}{(1-zv)^{2-s_1}}\,\mathrm dv, $$ which upon substituting $w=v/t$ gives $$ J^\prime=\frac{t^{2-s_1}}{q_2-q_1}\int_0^1 \frac{w^{(2-s_1)-1} (1-w)^{(3-s_1)-(2-s_1)-1}}{(1-ztw)^{2-s_1}}\,\mathrm dw. $$ Comparing this result to DLMF 15.6.1 then allows us to evaluate $J^\prime$ in terms of the hypergeometric function $_2F_1(2-s_1,2-s_1;3-s_1;\cdots)$. The form given by WA can then be obtained using transformation of variables for the hypergeometric function.