If you have lots of samples, you can include gene expression directly in the cox regression via e.g., ~ sex + age + stage + gene1_expr + gene2_expr + .... (in survminer or other package) Technically you only need 1 more sample than total number of variables, but for 15 total genes I would hesitate to take this approach for anything less than 100 samples.

Another problem is that many clinicians find it difficult to interpret survival data with continuous independent variables; it's a lot easier to discretize the expressions into 3 bins (high expression, low expression, background) -- usually 20%/60%/20% or 25%/50%/25% to plot survival curves. Taking this approach would expand 15 quantitative variables to 45 binary variables; but in this way you can directly look for expression bins that are prognostic.

The most typical approach is to define an "expression score" from all of the selected genes. This is basically taking the coefficients you got from rbsurv and using them to form a weighted average of expression to generate a single score. You can then cut that score into tertiles or quartiles to compare survival curves.

However, because you have used rbsurv to select these genes, you must ignore all of the statistics generated from these regressions (unless you are using a new set of samples) as they will be miscalibrated. To obtain calibrated statistics from this approach, you will need to use a permutation approach, followed by rbsurv selection and then the post-hoc regression. How to permute appropriately in the presence of clinical covariates is tricky and worth its own post.