  (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.
