Can interaction estimates be affected by a total number of events I am testing the interaction of two dummy variables, say Var1 and Var2 in a set of cancers. 
For each cancer type I want to model the number of response events Y:
Y ~ NB(Var1 + Var2 + Var3 + Var1:Var2 + offsetVar)

Var2 here stratifies labels of different cancers, for example, when I am looking at ovarian cancer Var2 = "1" for cases where data comes from ovarian cancer patients and Var2 = "0" for non-ovarian. Var3 is needed to control for confounding variable and offsetVar is used to be able to control for the maximal number of events Y that could have happened.
Regression is performed separately for each cancer.
In the end, I would like to answer a question, if there any cancer-type-specific input in the frequencies of Y and what is their effect size. Therefore, I have to look at the coefficients before Var1:Var2 variable.
I have a different amount of data available for different cancers - however, I don't expect the interaction term to be dependent on this but rather on some intrinsic features of different cancers.
However, after I perform a PCA on the set of coefficients where rows are cancer types and columns are coefficient estimates for different genes I have first PC (~ 50% of variance) and the order of cancers along this PC is matching the ordering by a total number of Y events per cancer type.
If I understand everything correctly, this shouldn't be happening. If you have any idea of what is going on, I would be very happy to hear it.
 A: I would put aside your concerns about what PCA seems to be showing and your fear of loss of information about reference levels of categorical covariates. Build a single model that incorporates all your information about cancer types and covariates. That includes all the information you need while giving you more power by combining information from all of your cases.
You don't lose information about the reference level of a categorical variable when you do this type of analysis. That's an easy misunderstanding based on the usual way that results for (generalized) linear models are presented. For a categorical variable in R, for example, the default is to show coefficients for the differences of each of the other levels from the reference level. So it might seem that information about the reference level has been lost.
That information about the reference level is included the intercept, which reports the predicted outcome when all categorical predictors are at their reference levels and continuous predictors have values of 0. You can then test any combinations of predictors that you want by using the information provided by the regression coefficients (including the intercept) and the covariance matrix among those coefficients. The covariance matrix is typically hidden from view unless you ask for it, but is used by other functions to help perform further tests.
You can avoid the difficulties in designing your own contrast matrix by taking advantage of well-vetted software packages that do the work for you. I have found the anova.rms() and contrast.rms() functions in the rms package in R to be very useful (once I spent the time to understand the syntax). The emmeans package is another useful package for examining model predictions for combinations of predictor values. If you are interested in the magnitude of the difference of one cancer type against the average of all cancer types, once you have the object produced by the emmeans() function applied to your combined model, the default eff.emmc setting for contrasts does exactly that.
For example, if you use the rms package you would fit the model and use the anova.rms() function to determine whether there are any significant interaction terms, in a single test that evaluates all interaction coefficients. If not, you could consider removing the interaction terms and then very simply interpreting the Var1, Var2, and Var 3 coefficients. If there are significant Var1:Var2 interactions then you can't easily interpret either of the Var1 or Var2 coefficients. In that case you specify the particular combinations of predictor values you want to compare to the contrast.rms() function.
Finally, I'm a bit worried about the design and ultimate interpretation. Please be on the lookout for survivorship bias, a potential problem in any study involving cases that may be lost, as in cancer studies. As a trivial example, if the "events" you are examining are trips to a clinician then those with more severe disease might have fewer events simply because they didn't live long enough. I suspect that the events of interest in your study are less trivial, but you still need to be on the lookout for survivorship bias. You don't need to have full time-series data, but information about the times to events or the times between events could be very important to incorporate into your model.
