Hi everyone, this is my first time posting here so I hope I'm not breaking any rules. I am quite a bioinformatics newbie.

My goal is to do a Cox survival analysis based on TCGA data looking into a specific set of gene expression levels. I'd like to do this based on continuous data, so without stratifying data into high/low. I haven't been able to find a tutorial that covers all of this step-by-step, so I built something myself. However, I'm not sure whether this is "bioinformatics-legal": I'd love to know if there are problems with anything I did. I am planning on turning this into a multivariate analysis if everything's okay, incorporating age at diagnosis and possibly disease stage.

Here is my code (with most data quality checks & cleanup removed for readability):

#-----------------------------------------------
# Prepare data
#-----------------------------------------------
# Load from disk    
tcga_data <- readRDS(
  file.path(data_dir, "data", "tcga_data_processed.rds")
)

# Extract counts
tcga_counts <- assay(tcga_data, "unstranded")
tcga_clin <- as.data.frame(colData(tcga_data))
tcga_genes <- as.data.frame(rowData(tcga_data))

# Normalise (TMM + voom)
dge <- DGEList(counts = tcga_counts)
dge <- calcNormFactors(dge)

design <- model.matrix(~ 1, data = tcga_clin)
v <- voom(dge, design, plot = FALSE)

expr_voom <- v$E   # log2 CPM

# Add gene symbols
rownames(expr_voom) <- gsub("\\.\\d+$", "", rownames(expr_voom))

tcga_genes$gene_id <- gsub("\\.\\d+$", "", rownames(tcga_genes))

symbol_map <- tcga_genes$gene_name
names(symbol_map) <- tcga_genes$gene_id

gene_symbols <- symbol_map[rownames(expr_voom)]

#-----------------------------------------------
# Targets
#-----------------------------------------------
# List targets
genes_of_interest <- c(
  "GENE1","GENE2", "GENE3",
)

# Match target symbols
gene_mask <- gene_symbols %in% genes_of_interest
expr_subset <- expr_voom[gene_mask, ]

rownames(expr_subset) <- gene_symbols[gene_mask]

# Transpose to patient level
expr_df <- as.data.frame(t(expr_subset))
expr_df$sample_id <- rownames(expr_df)

#-----------------------------------------------
# Survival
#-----------------------------------------------
# Create survival parameters
tcga_clin$OS_time <- ifelse(
  is.na(tcga_clin$days_to_death),
  tcga_clin$days_to_last_follow_up,
  tcga_clin$days_to_death
)

# Merge
tcga_clin$sample_id <- rownames(tcga_clin)

# Merge expression + clinical data
data <- inner_join(tcga_clin, expr_df, by = "sample_id")

#-----------------------------------------------
# Survival object + Cox per gene
#-----------------------------------------------
surv_object <- Surv(time = data$OS_time,
                    event = data$OS_status)

results <- lapply(genes_of_interest, function(gene){

  formula <- as.formula(
    paste("surv_object ~", gene)
  )

  model <- coxph(formula, data = data)
  s <- summary(model)

  data.frame(
    gene = gene,
    HR = s$coefficients[1, "exp(coef)"],
    p = s$coefficients[1, "Pr(>|z|)"]
  )
})

results_df <- do.call(rbind, results)
results_df$FDR <- p.adjust(results_df$p, method = "fdr")

I'd be much obliged if anyone can identify any problems with what I'm doing, as it seems to work perfectly but I want to make sure it's really okay. I'd love to improve my knowledge and bioinformatics skills, and am most unsure about the whole Cox analysis of gene expression counts thing