2
$\begingroup$

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:

enter image description 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
$\endgroup$
4
  • 1
    $\begingroup$ Please edit the question to show the model that you used, with descriptions of the predictors involved. This type of behavior can occur when there are categorical predictors in a model, but it will be hard to say for sure without more information. Also, note that the p-values aren't that small (none even less than 0.003) and probably aren't significant when you take into account the 20,000 comparisons you have made. $\endgroup$
    – EdM
    Commented Jan 19, 2023 at 15:50
  • $\begingroup$ Thank you! My model does have categorical predictors and I edited my question to include the model used. I appreciate your thoughts. $\endgroup$ Commented Jan 24, 2023 at 0:04
  • $\begingroup$ Try plotting the hazard rate conditional upon the stage; there are roughly seven curves that would match the seven stages. $\endgroup$
    – jbowman
    Commented Jan 24, 2023 at 1:16
  • $\begingroup$ It looks like the criteria for indicator=1, "low gene expression level for the gene of interest and genomic evidence of virus in tumor," are very strict. It's the number of cases with events that primarily determine the standard errors of coefficient estimates in Cox models, and thus the p values for hazard ratios. For example, for the gene TEDC1, how many cases had the combination of the following: its indicator=1, complete data on the other predictors, and an event rather than a censored survival time? $\endgroup$
    – EdM
    Commented Jan 24, 2023 at 13:52

1 Answer 1

2
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ This is a very clear and thoughtful explanation. Looking through my data, you are correct that all the genes in the leftmost curve had only 1 case where indicator = 1 and an event occurred, all the genes in the next curve had 2 such cases and so on. Your point about the importance of correcting for multiple hypothesis testing here is also well taken. Thank you so much for clearing this up for me. $\endgroup$ Commented Jan 25, 2023 at 6:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.