Not yet reviewed by Pith; the record is open.

(Abridged) Numerical simulation of strange quark stars (QSs) is challenging due to the strong density discontinuity at the stellar surface. In this paper, we report successful simulations of rapidly rotating QSs and study their oscillation modes in full general relativity. Building on top of the numerical relativity code \texttt{Einstein Toolkit}, we implement a positivity-preserving Riemann solver and a dust-like atmosphere to handle the density discontinuity at the surface. We demonstrate the robustness of our numerical method by performing stable evolutions of rotating QSs close to the Keplerian limit and extracting their oscillation modes. We focus on the quadrupolar $l=|m|=2$ $f$-mode and study whether they can still satisfy the universal relations recently proposed for rotating neutron stars (NSs). We find that two of the three proposed relations can still be satisfied by rotating QSs. For the remaining broken relation, we propose a new relation to unify the NS and QS data by invoking the dimensionless spin parameter $j$. The onsets of secular instabilities for rotating QSs are also studied by analyzing the $f$-mode frequencies. Same as the result found previously for NSs, we find that QSs become unstable to the Chandrasekhar-Friedman-Schutz instability when the angular velocity of the star $\Omega \approx 3.4 \sigma_0$ for sequences of constant central energy density, where $\sigma_0$ is the mode frequency of the corresponding nonrotating configurations. For the viscosity-driven instability, we find that QSs become unstable when $j\approx 0.881$ for both sequences of constant central energy density and constant baryon mass. Such a high value of $j$ cannot be achieved by realistic rotating NSs before reaching the Keplerian limit.