Show that $(X+Y+Z)^n = X^n + Y^n + Z^n$ if $XY = qYX$ and $YZ = q ZY$ and $ZX = qX Z$ Let $XY = qYX$ and $YZ = q ZY$ and $ZX = qX Z$ for $q^n = 1$ a root of unity.  Does it follow that:
$$(X+Y+Z)^n = X^n + Y^n + Z^n $$
In other words, all the middle terms cancel out.  If $q^m \neq 1$ for $0 < m < n$ (that is a primtive root of unity) at least we could group into monomial terms.
$$ \sum_{a+b+c = n} \Big[ \sum q^{\dots }\Big] X^a Y^b Z^b $$
After I'm not 100% sure the middle terms all cancel out.  Since neither the left and right sides change under the map $q \mapsto q^k$ (for $k \neq 1$) the only possible invariants are constants... so these coefficients certainly do not depend on $q$ but they might be zero.
If $q=1$ the middle numberes are the multinomial coefficients:
$$ (X+Y+Z)^n = \sum_{a+b+c = n} \binom{n}{a,b,c} X^a Y^b Z^b  $$
Let $n = p$ a  prime number if it makes this easier, so that $q = e^{\frac{2\pi i }{p}}$.  If $p = 5$  we'd have something like:
$$ (X+Y+Z)^5 = X^5 + Y^5 + Z^5 $$
 A: I assume that $ZX = qXZ$ should be $XZ = qZX$ instead.
Let $$\binom{n}{k}_q=\frac{\prod_{i=1}^n(1-q^i)}{\prod_{i=1}^k(1-q^i)\cdot\prod_{i=1}^{n-k}(1-q^i)}$$ denote the $q$-analog of the binomial coefficients. These satisfy the $q$-Pascal identities
$$
\binom{n}{k}_q = q^k\binom{n-1}k_q + \binom{n-1}{k-1}_q=\binom{n-1}k_q + q^{n-k}\binom{n-1}{k-1}_q\tag1
$$
which shows that $\binom{n}{k}_q$ is a polynomial in $q$, not obvious from its definition.
You can $(1)$ to prove by induction on $n$ that, whenever $X$ and $Y$ are non-commuting variables satisfying $XY=qYX$ for $q$ commuting with $X$ and $Y$, then 
$$
(X+Y)^n=\sum_{k\ge0} \binom{n}k_qX^kY^{n-k}\tag2
$$
Specializing to the case when $q^n=1$, $(2)$ proves that $(X+Y)^n=X^n+Y^n$, since $\binom{n}k_q$ has a factor of $(q^n-1)$ in the numerator which is only canceled in the denominator if $k=0$ or $k=n$.
Finally, this is entended to three variables via
$$
\begin{align}
(X+Y+Z)^n = X^n+(Y+Z)^n = X^n+Y^n+Z^n,
\end{align}
$$
where the first equality follows since $X(Y+Z)=q(Y+Z)X$. 
A: Just to answer the original question (with $ZX = qXZ$): the answer is "No". Here is a counterexample for $n = 3$, computed using SageMath (containing a lot of my old legacy code):
class ListOfVectors():

    # List of vectors in some given module (not necessarily a concrete
    # vector space). Certain methods, however, will require the ground
    # ring to be a field, or at least certain elements of it to be
    # invertible.
    # 
    # This probably implements some functionality from
    # https://trac.sagemath.org/ticket/11111 , but it does so
    # in a rather clumsy way.

    def __init__(self, basering, module, xs):
        # Create a list of vectors.
        # Syntax: ``ListOfVectors(R, M, xs)``, where
        # ``R`` is the base ring, ``M`` is the module
        # containing the vectors, and ``xs`` is the list
        # of the vectors.
        # NOTE: Make sure ``xs`` is never mutated. The
        # safe way is to use
        # ``ListOfVectors(R, M, xs.copy())``.
        self.vectors = xs
        self.basering = basering
        self.undermodule = module
        return object.__init__(self)

    def __repr__(self):
        return "List of vectors " + str(self.vectors) + " in " + str(self.basering) + "-module " + str(self.undermodule)

    def base_ring(self):
        return self.basering

    def underlying_module(self):
        return self.undermodule

    def list(self):
        return self.vectors[:]

    def show(self):
        print('List of vectors in module:')
        print self.undermodule
        print('over the base ring')
        print self.basering
        print('The vectors are:')
        for i in self.vectors: print i

    def reduction_blind(self, v):
        # Computes the reduction of a vector v modulo self, under the
        # assumption that self is in echelon form already.
        # It is assumed that the module knows division by elements of
        # the base ring and equality checking.
        w = v
        xs = self.vectors # the vectors in the list - in echelon form,
                          # with highest terms strictly decreasing.
        for t in xs:
            wcoeffs = w.monomial_coefficients()
            tcoeffs = t.monomial_coefficients().items()
            tleader = max(tcoeffs)
            if tleader[0] in w.monomial_coefficients():
                w = w - (wcoeffs[tleader[0]] / tleader[1]) * t
        return w

    @cached_method
    def echelon(self):
        # Returns an echelon form of self. This is a list of vectors
        # which spans the same submodule as self, but has its
        # sequence of highest terms strictly decreasing.
        # It is assumed that the module knows division by elements of
        # the base ring and equality checking.
        echeloned = [] # The echeloned variable contains a list of
                       # vectors already brought into echelon form.
                       # This list will grow step by step until it
                       # spans the same submodule as self.
        xs = self.vectors
        for v in xs:
            w = v
            # Now reduce v modulo echeloned:
            for t in echeloned:
                wcoeffs = w.monomial_coefficients()
                tcoeffs = t.monomial_coefficients().items()
                tleader = max(tcoeffs)
                if tleader[0] in wcoeffs:
                    w = w - (wcoeffs[tleader[0]] / tleader[1]) * t
            # Now w is the reduction of v modulo echelon.
            # If w == 0, then v was linearly dependent on echelon,
            #      and we don't have to do anything.
            # If w != 0, then we now add w to echelon.
            if w != 0:
                echeloned.append(w)
                # w might have been added to the wrong place, so
                # let us sort. I know this is not the best way;
                # if the bisect module would allow for keys, then
                # this would be easy to improve.
                echeloned.sort(key=lambda x : max(x.monomial_coefficients().keys()), reverse=True)
        res = ListOfVectors(self.basering, self.undermodule, echeloned)
        res.echelon.set_cache(res)
        return res

    def echelon2(self):
        # Another way to get an echelon form. Possibly sometimes
        # faster than echelon. It saves the result into the cache
        # of echelon.
        from bisect import insort
        echeloned = [] # The echeloned variable contains a list of
                       # pairs (i, v), with i being a monomial and v
                       # being a vector such that the highest monomial
                       # appearing in v with nonzero coefficient is i.
                       # This list will grow step by step. At every
                       # step k, the vectors v appearing in the list
                       # will be a basis of some submodule of the
                       # module spanned by self (namely, of the
                       # submodule spanned by the first k vectors of
                       # self). At the end, the vectors v appearing in
                       # the list will be a basis of the submodule
                       # spanned by self. At every step, the list will
                       # be sorted by increasing i, so the v's form a
                       # basis in row echelon form (whence the name
                       # "echelon").
        xs = self.vectors
        for v in xs:
            w = v
            # Now reduce v modulo the span of the vectors in echeloned:
            for i, t in reversed(echeloned):
                wcoeffs = w.monomial_coefficients()
                if i in wcoeffs:
                    tleader = t.monomial_coefficients()[i]
                    w = w - (wcoeffs[i] / tleader) * t
            # Now w is the reduction of v modulo the span of the
            # v's in echeloned.
            # If w == 0, then v was linearly dependent on the v's in
            #      echeloned, and we don't have to do anything.
            # If w != 0, then we now add w to echeloned.
            if w != 0:
                insort(echeloned, (max(w.monomial_coefficients().keys()), w))
        # Now forget the i's in echeloned. This is probably not very
        # efficient.
        res = ListOfVectors(self.basering, self.undermodule, [t for i, t in reversed(echeloned)])
        res.echelon.set_cache(res)
        return res

    @cached_method
    def echelon_reduced(self):
        # Returns the reduced echelon form of self. This is an
        # echelon form of self such that the highest term of
        # any of its vectors occurs in no other of its vectors.
        # It is assumed that the module knows division by elements of
        # the base ring and equality checking.
        ech = self.echelon().vectors
        k = len(ech)
        for i in range(k-1, -1, -1):
            w = ech[i]
            # Now reduce v modulo ech[i+1:]:
            for t in ech[i+1:]:
                wcoeffs = w.monomial_coefficients()
                tcoeffs = t.monomial_coefficients().items()
                tleader = max(tcoeffs)
                if tleader[0] in wcoeffs:
                    w = w - (wcoeffs[tleader[0]] / tleader[1]) * t
            ech[i] = w
        res = ListOfVectors(self.basering, self.undermodule, ech)
        res.echelon.set_cache(res)
        return res

    @cached_method
    def echelon_reducer(self):
        # Returns the matrix whose i-th row contains the coefficients
        # of the i-th vector of self.echelon() with respect to the
        # spanning set self.
        # The matrix is returned as a list of lists (the inner lists
        # being its rows).
        from itertools import izip
        echeloned = [] # The echeloned variable contains a list of
                       # pairs (v, hs). The v components are vectors
                       # already brought into echelon form. The hs
                       # component of each v contains the
                       # coefficients of v with respect to self.
                       # This list will grow step by step until its
                       # v components span the same submodule as
                       # self.
        xs = self.vectors
        zero = self.basering.zero()
        one = self.basering.one()
        for (i, v) in enumerate(xs):
            w = v
            whs = [zero] * len(self.vectors)
            whs[i] = one
            # Now reduce v modulo echelon:
            for (t, hs) in echeloned:
                wcoeffs = w.monomial_coefficients()
                tcoeffs = t.monomial_coefficients().items()
                tleader = max(tcoeffs)
                if tleader[0] in wcoeffs:
                    coeff = wcoeffs[tleader[0]] / tleader[1]
                    w = w - coeff * t
                    # Subtract coeff * hs from the vector whs:
                    whs = [whs_i - coeff * hs_i
                           for (whs_i, hs_i) in
                           izip(whs, hs)]
            # Now w is the reduction of v modulo echelon.
            # If w == 0, then v was linearly dependent on echelon,
            #      and we don't have to do anything.
            # If w != 0, then we now add w to echelon.
            if w != 0:
                echeloned.append((w, whs))
                # w might have been added to the wrong place, so
                # let us sort. I know this is not the best way;
                # if the bisect module would allow for keys, then
                # this would be easy to improve.
                echeloned.sort(key=lambda x : max(x[0].monomial_coefficients().keys()), reverse=True)
        return [e[1] for e in echeloned]

    @cached_method
    def syzygies(self):
        # Returns a dictionary whose items `i: xs` stand for the
        # syzygies of self. More specifically, `i: xs` means that
        # the `i`-th vector in self equals the linear combination
        # of the first `i-1` vectors with coefficients taken from
        # `xs`. (Note that `xs` is a length-`i-1` list.)
        from itertools import izip
        echeloned = [] # The echeloned variable contains a list of
                       # pairs (v, hs). The v components are vectors
                       # already brought into echelon form. The hs
                       # component of each v contains the
                       # coefficients of v with respect to self.
                       # This list will grow step by step until its
                       # v components span the same submodule as
                       # self.
        syzzies = {}
        xs = self.vectors
        for (i, v) in enumerate(xs):
            w = v
            whs = [self.basering.zero()] * len(self.vectors)
            whs[i] = self.basering.one()
            # Now reduce v modulo echelon:
            for (t, hs) in echeloned:
                wcoeffs = w.monomial_coefficients()
                tcoeffs = t.monomial_coefficients().items()
                tleader = max(tcoeffs)
                if tleader[0] in wcoeffs:
                    coeff = wcoeffs[tleader[0]] / tleader[1]
                    w = w - coeff * t
                    # Subtract coeff * hs from the vector whs:
                    whs = [whs_i - coeff * hs_i
                           for (whs_i, hs_i) in
                           izip(whs, hs)]
            # Now w is the reduction of v modulo echelon.
            # If w == 0, then v was linearly dependent on echelon,
            #      and we record a syzygy.
            if w == 0:
                syzzies[i] = [-s for s in whs[:i]]
            # If w != 0, then we now add w to echelon.
            else:
                echeloned.append((w, whs))
                # w might have been added to the wrong place, so
                # let us sort. I know this is not the best way;
                # if the bisect module would allow for keys, then
                # this would be easy to improve.
                echeloned.sort(key=lambda x : max(x[0].monomial_coefficients().keys()), reverse=True)
        return syzzies

    def coefficients(self, v):
        # Computes a list [a_1, a_2, ..., a_k] of coefficients
        # such that v = a_1 v_1 + a_2 v_2 + ... + a_k v_k, where
        # self == [v_1, v_2, ..., v_k].
        ech = self.echelon()
        l = len(ech.vectors)
        ring = self.basering
        blinds = ech.coefficients_blind(v)
        k = len(self.vectors)
        res = [self.basering.zero()] * k
        red = self.echelon_reducer()
        return [ring.sum((blinds[i] * red[i][j]
                         for i in range(l)))
                for j in range(k)]

    def coefficients_blind(self, v):
        # Computes a list [a_1, a_2, ..., a_k] of coefficients
        # such that v = a_1 v_1 + a_2 v_2 + ... + a_k v_k, where
        # self == [v_1, v_2, ..., v_k], under the assumption that
        # self is in echelon form, and that v is a linear
        # combination of the vectors of self.
        w = v
        coeffs = [self.basering.zero()] * len(self.vectors)
        for (i, t) in enumerate(self.vectors):
            wcoeffs = w.monomial_coefficients()
            tcoeffs = t.monomial_coefficients().items()
            tleader = max(tcoeffs)
            if tleader[0] in wcoeffs:
                coeff = wcoeffs[tleader[0]] / tleader[1]
                w = w - coeff * t
                coeffs[i] = coeff
        return coeffs

    def reduction(self, v):
        # Computes the reduction of a vector v modulo self.
        return self.echelon().reduction_blind(v)

    def contains_blind(self, v):
        # Finds out whether a vector v lies in the submodule generated
        # by self, under the assumption that self is in echelon form
        # already.
        return (self.reduction_blind(v) == 0)

    def contains(self, v):
        # Finds out whether a vector v lies in the submodule generated
        # by self.
        return (self.echelon().reduction_blind(v) == 0)

    def isbiggerthan(self, anotherlist, verbose=False):
        # Finds out whether the submodule spanned by self contains that
        # spanned by anotherlist (another list of vectors).
        xs = self.echelon()
        for y in anotherlist.list():
            if not xs.contains_blind(y):
                if verbose == True:
                    print "The offending vector is: "
                    print y
                return False
                break
        return True

    def issmallerthan(self, anotherlist, verbose=False):
        # Finds out whether the submodule spanned by self is contained
        # in that spanned by anotherlist (another list of vectors).
        return anotherlist.isbiggerthan(self, verbose=verbose)

    def isequivalentto(self, anotherlist, verbose=False):
        # Finds out whether the submodule spanned by self equals
        # that spanned by anotherlist (another list of vectors).
        # (For instance, self.isequivalentto(self.echelon()) should
        # always return True.)
        return (anotherlist.isbiggerthan(self, verbose=verbose)
                and self.isbiggerthan(anotherlist, verbose=verbose))

    def add(self, anotherlist):
        # Gives the disjoint union of self with anotherlist (another list
        # of vectors, which of course should be in the same module).
        # This union spans the sum of the respective submodules.
        us = self.vectors[:]
        us.extend(anotherlist.list())
        return ListOfVectors(self.basering, self.undermodule, us)

    def intersection_blind(self, anotherlist):
        # Gives the intersection of the span of self with the span
        # of anotherlist (another list of vectors, which of course
        # should be in the same module), as an echelonized list of
        # vectors.
        # This assumes that self and anotherlist are in echelon form
        # already.
        vs = self.vectors
        ws = anotherlist.list()
        R = self.basering
        M = self.undermodule
        n = len(vs)
        m = len(ws)
        vsws = vs + ws
        syzzies = ListOfVectors(R, M, vsws).syzygies()
        reslist = []
        for (k, syz) in syzzies.iteritems():
            # Throw the vs-part of syz away.
            k2 = k - n
            syz2 = syz[n:]
            reslist.append(ws[k2] - M.sum(syz2[i] * ws[i] for i in range(k2)))
        return ListOfVectors(R, M, reslist)

    def intersection(self, anotherlist):
        # Gives the intersection of the span of self with the span
        # of anotherlist (another list of vectors, which of course
        # should be in the same module), as an echelonized list of
        # vectors.
        return self.echelon().intersection_blind(anotherlist.echelon())

    def product(self, anotherlist):
        # Gives the list formed by pairwise products of vectors in
        # self with vectors in anotherlist.
        # This, of course, makes sense only when the underlying module is
        # an algebra.
        # This new list spans the product of the respective submodules
        # (in the sense in which, e. g., the product of ideals is
        # defined).
        us = self.vectors
        vs = anotherlist.list()
        ws = [p * q for p in us for q in vs]
        return ListOfVectors(self.basering, self.undermodule, ws)

    def commutator(self, anotherlist):
        # Gives the list formed by pairwise commutators of vectors in
        # self with vectors in anotherlist.
        # This, of course, makes sense only when the underlying module is
        # an algebra.
        # This new list spans the commutator of the respective submodules.
        us = self.vectors
        vs = anotherlist.list()
        ws = [p * q - q * p for p in us for q in vs]
        return ListOfVectors(self.basering, self.undermodule, ws)

    def power(self, n):
        # Returns the n-th power of the list with respect to the
        # above-defined product function.
        if n == 0:
            return ListOfVectors(self.basering, self.undermodule, [self.undermodule.one()])
        elif n == 1:
            return self
        else:
            m = int(n) / int(2)
            M = n - m
            return self.power(m).product(self.power(M))

    def dimension(self):
        # Gives the dimension of the submodule generated by self.
        return len(self.echelon().list())

    def image(self, f):
        # Returns the list of the images of the vectors under a morphism f.
        # The module in which they lie is the codomain of f.
        ys = [f(x) for x in self.vectors]
        return ListOfVectors(self.basering, f.codomain(), ys)

    def kernel(self, f):
        # Returns a basis of the kernel of a morphism f (restricted
        # to the span of self).
        img = self.image(f)
        syzzies = img.syzygies()
        vects = self.vectors
        R = self.basering
        M = self.undermodule
        ker = []
        for i, xs in syzzies.iteritems():
            ker.append(vects[i] - M.sum((c * vects[j] for (j, c) in enumerate(xs))))
        return ListOfVectors(R, M, ker)

def subalgcomp(A, U, n):
    # INPUT:
    # A: a graded algebra over a base ring.
    # U: a list of listsofvectors, where each vector in U[i] lies in the
    # (i+1)-th graded component of A.
    # n: an integer.
    # OUTPUT:
    # a list of listsofvectors, the i-th of which spans the (i+1)-th
    # graded component of the subalgebra of A generated by U. The list
    # has length n.
    # (As should be clear from the above description, generators in the
    # 0-th graded component are not supported. This is because there is
    # no algorithm for finding the subalgebra generated by them in the
    # general case.)
    BR = A.base_ring()
    A = [0]*n
    l = len(U)
    if l > n: l = n    # forget about the generators that are too high for us to use
    for i in range(l):
        A[i] = U[i]
    for i in range(l, n):
        A[i] = ListOfVectors(BR, A, [])
    for i in range(n):
        for j in range(1,i+1):
            if i - j < l:
                A[i] = A[i].add(A[j-1].product(U[i-j])).echelon()
    return A

def gradedideal(A, U, V, n):
    # INPUT:
    # A: a graded algebra over a base ring.
    # U: a list of listsofvectors, where each vector in U[i] lies in the
    # (i+1)-th graded component of A.
    # V: a list of listsofvectors, where each vector in V[i] lies in the
    # vector subspace of A generated by U[i].
    # n: an integer.
    # OUTPUT:
    # a list of listsofvectors, the i-th of which spans the (i+1)-th
    # graded component of the two-sided ideal generated by V inside the
    # subalgebra of A generated by U. The list has length n.
    # (As should be clear from the above description, generators in the
    # 0-th graded component are not supported. This is because there is
    # no algorithm for finding the subalgebra generated by them in the
    # general case.)
    BR = A.base_ring()
    B = subalgcomp(A, U, n)
    A = [0]*n
    l = len(V)
    if l > n: l = n    # forget about the generators that are too high for us to use
    for i in range(l):
        A[i] = V[i]
    for i in range(l, n):
        A[i] = ListOfVectors(BR, A, [])
    for i in range(n):
        for j in range(1,i+1):
            A[i] = A[i].add(B[j-1].product(A[i-j]).add(A[i-j].product(B[j-1]))).echelon()
    return A

DefaultBase = CyclotomicField()
q = CyclotomicField().gen(3) # Primitive 3rd root of unity.

def FA(N):
    # Free algebra over the cyclotomic field in `N` generators.
    return FreeAlgebra(DefaultBase, N, names='x')

def FABasis(N, n):
    # Basis of the `n`-th graded component of the free algebra
    # over the cyclotomic field in `N` generators.
    FAN = FA(N)
    gens = FAN.gens()
    import itertools
    xs = [FAN.prod(w) for w in itertools.product(*([gens] * n))]
    return ListOfVectors(DefaultBase, FAN, xs)

A = FA(3) # Free algebra on 3 generators
X, Y, Z = A.gens()

def MSIdeal(n):
    # The `i`-th graded components, for `i = 1, 2, ..., n`,
    # of the ideal of `A` generated by
    # `XY - qYX, YZ - qZY, ZX - qXZ`.
    nothing = ListOfVectors(DefaultBase, A, [])
    gs = [X * Y - q * Y * X, Y * Z - q * Z * Y, Z * X - q * X * Z]
    gslist = ListOfVectors(DefaultBase, A, gs)
    V = [(gslist if h == 1 else nothing) for h in range(n)]
    U = [FABasis(3, h+1) for h in range(n)]
    return gradedideal(A, U, V, n)

MSI3 = MSIdeal(3)

MSI3[2].contains((X+Y+Z)**3 - X**3 - Y**3 - Z**3)

The answer is "False", which signifies that $\left(X+Y+Z\right)^3 - \left(X^3 + Y^3 + Z^3\right)$ is not in the ideal generated by the relations you gave.
Alternatively, you can check this by using the multigrading on the free algebra (in which the degree of a word $w$ is $\left(\text{number of $X$'s in $w$}, \text{number of $Y$'s in $w$}, \text{number of $Z$'s in $w$}\right)$). Clearly, if $\left(X+Y+Z\right)^3 - \left(X^3 + Y^3 + Z^3\right)$ would lie in the ideal, then the multidegree-$\left(1,1,1\right)$ component $XYZ+YZX+ZXY+XZY+YXZ+ZYX$ of $\left(X+Y+Z\right)^3 - \left(X^3 + Y^3 + Z^3\right)$ would also lie in this ideal (because the ideal is homogeneous with respect to this multigrading). But an easy verification shows that this is not the case.
