The following (hopefully) settles the question. Fix $n$, and let $V(x)$ be the polynomial from the question, for which we want to show that it is
positive on the open interval $(0,1)$. Set
\begin{align*}
a &= 64n^2 - 256/3n + 20\\
b &= 768n^3 - 1024n^2 + 272n - 24\\
c &= 14336/3n^4 - 306304/45n^3 + 35792/15n^2 - 6104/15n + 32
\end{align*}
Then
\begin{equation}
\frac{V(x)}{(1-x^2)^4x^{2n-2}}=(1-x)^3H(x)+(n+1)^2(2n+1)^2x^{22n-4}
(ax^2+bx(1-x)+c(1-x)^2)
\end{equation}
for a polynomial $H(x)$.

Here $a,b,c$ are positive for $n\ge2$. We show that *all*
coefficients of $H(x)$ are positive whenever $n\ge6$.

I'm not sure about the easiest kind to prove the assertion about the coefficients of $H(x)$. But the following works: Multiplying the equation by $(1-x)^4(1+x)^7$ yields
\begin{equation}
(1-x^2)^7H(x)=W(x)
\end{equation}
where $W(x)$ is a polynomial for which we have an explicit (though messy) expression. Write $H(x)=H_0(x^2)+xH_1(x^2)$ and $W(x)=W_0(x^2)+xW_1(x^2)$. Then $H_0(x)(1-x)^7=W_0(x)$ and $H_1(x)(1-x)^7=W_1(x)$. Thus we are reduced to show that $H_0(x)$ and $H_1(x)$ have nonnegative coefficients. We obtain these coefficients by multiplying $W_\star(x)$ with the formal power series $\frac{1}{(1-x)^7}=\sum\binom{k+6}{6}x^k$.

All what follows is assisted/confirmed by the Sage code below.

There is a finite set of pairs $(i,j)$ of nonnegative integers, *independent* of $n$, such that the exponents of $W(x)$ have the form $2ni+j$. Let $a_{i,j}$ be the corresponding coefficient. It is a polynomial in $n$. It turns out that $0\le i\le 11$ and $-4\le j\le 11$.

We have
\begin{equation} 
W = \sum_{i,j}a_{i,j}(n)x^{2ni+j}.
\end{equation}
Thus we obtain $W_0$ or $W_1$ according as to whether $j$ is even or odd.

The exponents of $W_0(x)$ now are $ni+j/2$ for even $j$, and those of $W_1(x)$ are $ni+(j-1)/2$ for odd $j$. For better readability, we call the modified $j$'s again $j$. So the exponents are $ni+j$ for certain pairs $(i,j)$, and the associated coefficients will be called $a_{i,j}$. After this modification, we have $-2\le j\le 5$.
 
If $n\ge8$, then $ni+j>ni'+j'$ if and only if $i>i'$ or $i=i'$ and $j>j'$. This gives a total order on the pairs $(i,j)$. Every exponent of $x$ in $W_\star(x)/(1-x)^7$ which matters has the form $ni+j+y$ where $0\le y<n(i'-i)+j'-j$, where $(i',j')$ is the successor of $(i,j)$ for our order.

The coefficient of $x^{ni_1+j_1+y}$ in $W_\star(x)\frac{1}{(1-x)^7}$
is the sum of all $a_{i_1,j_1}(n)\cdot\binom{ni+j-ni_1-j_1+y+6}{6}$ for all $(i_1,j_1)\le(i,j)$.

Thus we obtain *finitely* many (actually $126$) explicit polynomials $(n,y)$ in $n$ and $y$ (of small degrees in $n$ and $y$) for which we have to prove that they are positive for all $n\ge8$ and all $y$ in the specified range.

We distinguish two cases: If $(n'i+j')-(ni+j)=1$ ($(i',j')$ the successor of $(i,j)$ again), then only $y=0$ is possible. In all cases, we see that $f(n+6,0)$ has nonnegative coefficients. Thus $f(n,0)\ge0$ whenever $n\ge6$.

Now assume that $(n'i+j')-(ni+j)>1$. In these cases it turns out that $i'=i+1$, so $0\le y< n-e$ for some integer. Thus $y=n-e-1-o$ for some non-negative $o$. Thus $f(n,o)=f(y+e+1-o,o)=g(y,o)$. This latter polynomial, in $y$ and $o$, has to be non-negative for all non-negative $o$ and $y$. It turns out that except for a single example actually oll the coefficients are nonnegative. For this single exception, $g(y,o+1)$ has nonnegative coefficients. Thus $g(y,o)\ge0$ except possibly for $o=0$. But $g(y,0)$ has nonnegative coefficients, and things are settled.

Our argument used $n\ge8$. For $n=6$ and $7$, one can directly compute $H(x)$ and check that the coefficients are nonnegative. The cases $2\le n\le5$ are checked directly for $V(x)$.

*Remark (answering Jason's question):* It is a known fact that a polynomial $V$ is positive on $(0,1)$ if and only if it is a nonnegative linear combination of polynomials $x^i(1-x)^j$. The problem is that it may involve terms where $i+j$ is bigger than $\deg V$. I used an LP solver to play a little with $V$ and $\frac{V}{(1-x^2)^4x^{2n-2}}$. From that a pattern showed up which lead to a solution.

By the way, $V$ actually seems to be a nonnegative linear combination of $x^i(1-x)^{\deg V-i}$. This can be obtained from the recursion (fixed $n$)
\begin{equation}
V_0=V,\;\; V_{i+1}=\frac{V_i-V_i(1)x^{\deg V_i}}{1-x}.
\end{equation}
I have an explicit expression of the coefficients $a_i$ in $V=\sum a_ix^i(1-x)^{\deg V-i}$. But I don't see an easy argument showing positivity.

Here is the Sage code:

<pre>
# Formally compute W
var('x n')
P = (1+x^(4*n-1))*(1+x^(2*n))*(1-x^(2*n+1))
Q = (1+x^(2*n+1))*(1-x^(2*n+2))
V = P*Q^2*P.diff(x,2)+(P*Q.diff(x))^2-P^2*Q*Q.diff(x,2)-(P.diff(x)*Q)^2
a = 64*n^2 - 256/3*n + 20
b = 768*n^3 - 1024*n^2 + 272*n - 24
c = 14336/3*n^4 - 306304/45*n^3 + 35792/15*n^2 - 6104/15*n + 32
W = V*(1+x)^3*x^(2-2*n)-(1-x)^4*(1+x)^7*(n+1)^2*(2*n+1)^2*x^(22*n-4)*(a*x^2 + b*x*(1-x)+ c*(1-x)^2)

# Compute the pairs (i,j) for the exponents 2ni+j of W, and the
# corresponding coefficients a_ij. Somewhat clumsy, as I don't know how
# to deal with polynomials where exponents are symbolic expressions
#
# l is the list of summands of W
l = [z.canonicalize_radical() for z in W.expand().operands()]
aij = {}
for term in l:
    c = term(x=1) # get coefficient of term
    e = (term.diff(x)/c)(x=1) # get exponent of term
    i,j = ZZ(e.diff(n))//2,ZZ(e(n=0)) # Clumsy method to compute pairs (i,j)
    aij[i,j] = aij[i,j] + c if (i,j) in aij else c

# Check if the pairs (i,j) and coefficients aij[i,j] were correctly computed
Wnew = sum(c*x^(2*n*i+j) for (i,j),c in aij.items())
print (W-Wnew).canonicalize_radical().is_trivial_zero()

# pairs (i,j) and coeffcients for W_0 and W_1, respectively
aij0 = {(z[0],z[1]//2):c for z,c in aij.items() if is_even(z[1])}
aij1 = {(z[0],(z[1]-1)//2):c for z,c in aij.items() if is_odd(z[1])}

# compute and collect coefficients of W_0(x)/(1-x)^7 and W_1(x)/(1-x)^7
var('y o')
hcoeffs = []
for aij in [aij0,aij1]:
    aijk = aij.keys()
    aijk.sort()
    for k in range(len(aijk)-1):
        (i,j) = aijk[k]
        # f(n,y) is the coefficient of x^(ni+j+y) for 0 <= y < yr
        f = sum(aij[z]*binomial(n*i+j+y-n*z[0]-z[1]+6,6) for z in aijk[:k+1])
        (ii,jj) = aijk[k+1]
        yr = n*(ii-i)+(jj-j)
        hcoeffs.append((f,yr))

# check non-negativity of the polynomials in hcoeffs
trouble = [] # trouble list of polynomials, which need further checking
for k,ff in enumerate(hcoeffs):
    if ff[1] == 1:
        # Checks if the coefficient of x^(ni+j), which is a
        # polynomial in n, has nonnegative coefficients upon replacing
        # n with n+6
        f = ff[0](n=n+6,y=0).polynomial(QQ)
        if f == 0:
            continue
        if min(f.coefficients()) < 0:
            print "error"
            break
    else:
        # The coefficient of x^(ni+j+y), which is a polynomial in n and y,
        # has to be nonnegative for 0 <= y < yr=ff[1], so y = yr-1-o for
        # non-negative o. Here yr has the form n-e for an integer e. Thus
        # n = y + n - yr + 1 + o. So upon replacing n with y + n - yr + 1
        # + o, we get a polynomial in y and o which has to be nonnegative
        # for all y and o. In most cases, this holds because the
        # coeffcients are non-negative.
        f = ff[0](n=y+n-ff[1]+1+o).polynomial(QQ)
        if f == 0:
            continue
        if min(f.coefficients()) < 0:
            trouble.append(f)

print "exactly one trouble polynomial? :", len(trouble) == 1
f = trouble[0]
# f has is nonnegative for all o >= 1.
print min(f(o=o+1).polynomial(QQ).coefficients()) >= 0
# So only o=0 remains to check. We see that f(o=0) has nonnegative coeffcients
print min(f(o=0).coefficients()) >= 0
</pre>