1

My data represents the count of mutations by gene in two groups cases and ctrl. I would like to do comparison gene-wise between the two groups. But in my case I don’t think anova is appropriate since count data follows a Poisson distribution. Here a small subset of my Db, as you can see the samples size (cases and ctrls is different).

gene   Samples_ID value  Type
AB       1    28 Cases
AB     100    22 Cases
AB     101    36 Cases
AB     102    57 Cases
AB     105    29 Cases
AB     106    25 Cases
AB     108    23 Cases
AB 4928    18  Ctrls
AB 4929    18  Ctrl
AB 4930    24  Ctrl
AB 4931    20  Ctrl
AB 4932    25  Ctrl
AB 4933    15  Ctrl
AB 4934    25  Ctrl
AB 4935    22  Ctrl
AB 4936    30  Ctrl
AB 4937    15  Ctrl
AB 4938    18  Ctrl
AB 4939    21  Ctrl
FG       1    21 Cases
FG     100    16 Cases
FG     101    21 Cases
FG     102    34 Cases
FG     105    22 Cases
FG     106    23 Cases
FG     108    23 Cases
FG 4928     8  Ctrl
FG 4929     3  Ctrl
FG 4930     7  Ctrl
FG 4931     6  Ctrl
FG 4932     5  Ctrl
FG 4933    15  Ctrl
FG 4934     8  Ctrl
FG 4935    11  Ctrl
FG 4936     1  Ctrl
FG 4937     7  Ctrl
FG 4938     6  Ctrl
FG 4939     8  Ctrl
SYU       1    27 Cases
SYU     100    23 Cases
SYU     101    35 Cases
SYU     102    39 Cases
SYU     105    24 Cases
SYU     106    25 Cases
SYU     108    30 Cases
SYU 4928     5  Ctrl
SYU 4929     6  Ctrl
SYU 4930     6  Ctrl
SYU 4931    16  Ctrl
SYU 4932     5  Ctrl
SYU 4933    11  Ctrl
SYU 4934    12  Ctrl
SYU 4935    11  Ctrl
SYU 4936    15  Ctrl
SYU 4937     8  Ctrl
SYU 4938    13  Ctrl
SYU 4939    11  Ctrl

Where "gene" is the gene name, "Samples_ID" represent the patient, "value" is the number of mutation for the given gene, "type" is the group.
I have tried a generalized linear model as follows:

fit <- glm(value ~ gene + Type, data = file, family = poisson())

but I'm not convinced it's the right way.
PS
Please, give me a simple explanation, I'm not a statistician and my background about statistic is ~ 0.
Thank you in advance.
Best.

Update
The number of cases is 98
The number oc ctrls is 40
Number of genes is 2780
Mutations are RNA editing events.

Salvatore Di Giorgio
  • 11
  • 2
  • Can you give more information on the data structure. If the left column is the name of the gene, why do they differ between groups? (or is it just the location of your snapshopt?) – stefgehrig Jan 04 '22 at 18:30
  • They differ between groups because they came from different patient (Samples_ID columns). The column "gene" is the gene name, "values" is the number of mutation for the given gene, "type" column indicate the group. – Salvatore Di Giorgio Jan 04 '22 at 18:39
  • Then you could add the gene X Type interactions to the model you showed. You will have many interaction terms with test statistics and p-values then, however, because there are many genes and you get an interaction term with Type for each one, and you might want to correct for multiplicity. As an omnibus Test for all interaction terms, you could use LRT. However, you should make sure that Poisson is reasonable assumption, i.e. check for over or underdispersion of the counts (and have a look at the distribution of the counts in general - are they skewed?) – stefgehrig Jan 04 '22 at 18:45
  • 1
    How many genes and samples? By "mutations" do you mean germ-line variants or something like somatic mutations in a tumor? If the former, what's your reference genome for defining variants? If the latter, how are you ending up with so many mutations in a single gene of an individual? Please provide that information by editing the question, as comments are easy to overlook and can get deleted. – EdM Jan 05 '22 at 14:17
  • The formula `value ~ gene * Type` ought to do what you want. Alternatively, it sounds like you need to run your model `value ~ Type` separately for each gene. It's not perfectly clear, because your statement "there is a significant difference between the number of mutations for the same gene between the two groups" is ambiguous. Could you restate your objective more definitely? – whuber Jan 05 '22 at 20:24

1 Answers1

1

I don’t think anova is appropriate since count data follows a Poisson distribution.

You're correct that count data might not be handled properly by standard ANOVA, but they also aren't necessarily handled correctly by Poisson models. Poisson models assume independence of observations and no confounding variables omitted from the analysis.

If such assumptions don't hold, the standard Poisson equivalence of mean and variance might not hold. This issue is often addressed by negative binomial models that allow for more general associations between mean and variance of counts.

You also have issues of repeated measures on the same individuals that need to be taken into account for a potential lack of independence. Furthermore, with 2780 genes to test, you have a substantial multiple-comparison problem. Without correction for multiple comparisons, at p = 0.05 you would find 139 "significant" differences even if there weren't any true ones at all.

Negative binomial count models that take all those issues into account are used for RNA-seq gene expression data. They would seem to be ideally suited to your situation: count data, a small number of conditions, a large number of genes, and an intermediate number of individuals having data on each gene. Instead of re-inventing an analytical approach, see whether standard tools like those provided by the DESeq2 Bioconductor package can be adapted to your situation. (If you're getting into bioinformatics you should be familiar with such tools and their statistical underpinnings in any event.)

EdM
  • 57,766
  • 7
  • 66
  • 187
  • There are no issues of repeated measures on the same individuals because each gene appears once per sample. – Salvatore Di Giorgio Jan 10 '22 at 22:35
  • @SalvatoreDiGiorgio but each individual has values for 2780 genes. You are likely to have correlations among genes within individuals, maybe resulting from something as simple as a different baseline "mutation" rates among individuals. That should ideally be taken into account. I think that the negative binomial modeling of DESeq2 can handle that type of thing directly, although I haven't tried to apply it to this type of situation. – EdM Jan 10 '22 at 22:55