I have a large presence/absence matrix of genes from different strains of S. aureus that I obtained through microbiome shotgun sequencing. I conducted a PCA and observed differences between the strains living in the conditions (healthy and diseased). Now, I would like to find out which genes drive this variation. My presence/absence matrix comes with genes in columns and observations/samples in rows. I have a grouping variable as class factor, which specifies which samples belong to which condition (i. e. healthy or diseased).
> head(gene.matrix)
gene1 gene2 gene3
sample1 1 0 0
sample2 1 1 0
sample3 0 1 0
in total I have >3500 genes and 12 strains.
> grouping.variable
[1] healthy healthy diseased
How can I find out, which genes are (preferentially) present in the strains on diseased sites? Through a loop, I already tried a glm Call:
glm(formula = Gene_x ~ condition, family = binomial(link = "logit"), data = gene.matrix) and received a list with an entry for every single gene that holds the results of the glm call. Is this appropriate? Can I just filter out those genes with a pvalue <0.05? And If yes, how can I do that? Would I need to adjust the pvalue for multiple testing?
