Apparent artifact in R volcano plot of hazard ratio vs. p-value

Question

I've used R and the survival package's coxph() function to calculate hazard ratios and p-values for around 20,000 genes using data from The Cancer Genome Atlas. When I plot the hazard ratios and -log10(P) values in a volcano plot, there appears to be an artifact in the data - the data points all seem to fall along one of several curves, as shown here:

I'm wondering why this artifact is present. My best guess is that it has something to do with how HRs and p-values are calculated, either in the coxph() function or intrinsic to how Cox regression models work, but I'm not sure. I'd greatly appreciate any insights or thoughts about this. Thank you!

1/23/23: Edited to add the model I am using, with descriptions of the terms:

summaryCoxph <- 
  summary(coxph(
     Surv(times, died) ~ indicator + grade + stage + sex + age
  )
)

• indicator: An indicator variable which is 1 if the patient has a low gene expression level for the gene of interest and genomic evidence of virus in tumor, and 0 otherwise
• grade: Tumor grade (1, 2, 3 or 4)
• stage: Tumor stage (I, II, IIIA, IIIB, IIIC, IVA or IVB)
• sex: Male or female
• age: Numeric age

Answer 1

These are the hazard ratios and corresponding p-values for the indicator = 1 predictor for each of the 20,000 genes analyzed separately. This type of display can happen with categorical predictors.

I suspect that you have very small numbers of cases that both meet the indicator = 1 binary criterion for a gene (low gene expression and evidence of virus in the tumor) and experienced events. In a Cox model, it's the number of cases with events that provides power, not the total number of cases themselves. The separate curves might then be associated with each of those small numbers of cases (1, 2, 3,...) that meet that criterion.

The green dot at the lower left is a hazard ratio of 0 with a p value of 1. I suspect that includes all the genes for which no cases with events had indicator = 1. That's what you get when there aren't events for a level of a binary predictor.

I suspect that the leftmost curve represents genes for which exactly 1 case has both an event and indicator = 1. Depending on the clinical characteristics of the case meeting those criteria for a given gene, you get a range of hazard ratios and p-values.

That curve passes the horizontal dashed line (presumably representing p = 0.05) at a hazard ratio of about 0.2. It's easier to think about this in the original log-hazard coefficient scale, where that is equivalent to a coefficient of -1.6. To pass significance at p < 0.05, the standard error of that coefficient estimate needs to be about half of that in magnitude, or about 0.8 standard error.

Now look at the next two curves toward the right and see the hazard ratios at which they pass the same horizontal line. I estimate HRs of about 0.3 (coefficient of -1.2) and 0.375 (coefficient of -0.98). The corresponding standard errors at significance (where they cross that line) are thus about 0.6 and 0.49, respectively.

As a first approximation in a Cox model, the standard error of a coefficient estimate decreases with the square root of the number of corresponding observations. The horizontal spacing of the curves at p = 0.05 is very close to that relationship. Based on my estimate for 1 case and a standard error of 0.8, for 2 cases you would have a standard error about $0.8/\sqrt{2}\approx 0.57$, and for 3 cases $0.8/\sqrt{3}=0.46$.

Also, note that none of your individual p-values is less than about 0.003. That's unlikely to pass a significance test once there's a correction for multiple comparisons. If my suspicion is correct, you simply do not have enough cases to draw any reliable conclusions about the association of indicator = 1 with outcome.

This type of single-gene association with clinical outcome is tricky. You admirably tried to account for some known clinical associations with outcome. You might have made things harder, however, by using an all-or-none criterion for low gene expression. It's typically better not to categorize a continuous predictor like gene expression, but rather to model it flexibly.