Many integrable stochastic particle systems in one space dimension (such as TASEP - totally asymmetric simple exclusion process - and its various deformations, with a notable exception of ASEP) remain integrable when we equip each particle $x_i$ with its own jump rate parameter $ν_i$. It is a consequence of integrability that the distribution of each particle $x_n(t)$ in a system started from the step initial configuration depends on the parameters $ν_j$, $j\le n$, in a symmetric way. A transposition $ν_n \leftrightarrow ν_{n+1}$ of the parameters thus affects only the distribution of $x_n(t)$. For $q$-Hahn TASEP and its degenerations ($q$-TASEP and directed beta polymer) we realize the transposition $ν_n \leftrightarrow ν_{n+1}$ as an explicit Markov swap operator acting on the single particle $x_n(t)$. For beta polymer, the swap operator can be interpreted as a simple modification of the lattice on which the polymer is considered. Our main tools are Markov duality and contour integral formulas for joint moments. In particular, our constructions lead to a continuous time Markov process $\mathsf{Q}^{(\mathsf{t})}$ preserving the time $\mathsf{t}$ distribution of the $q$-TASEP (with step initial configuration, where $\mathsf{t}\in \mathbb{R}_{>0}$ is fixed). The dual system is a certain transient modification of the stochastic $q$-Boson system. We identify asymptotic survival probabilities of this transient process with $q$-moments of the $q$-TASEP, and use this to show the convergence of the process $\mathsf{Q}^{(\mathsf{t})}$ with arbitrary initial data to its stationary distribution. Setting $q=0$, we recover the results about the usual TASEP established recently in [arXiv:1907.09155] by a different approach based on Gibbs ensembles of interlacing particles in two dimensions.