I'm trying to analyse RNAseq data using the DESeq package, however i am unsure on what the right design would be for my experimental setup. I have three different treatments, where i am always comparing the effect of the treatment on either KO or WT. I have separated these into three different dataframes and am running the DESeq analysis on them separately.
Here is an example code of what I am using for all three comparisons. Since in all three cases I essentially want to compare the same thing, I am using the same setup for all three. Based on what I understood from reading different vignettes and guides, using just genotype in the design should be enough. I am also setting the WT as baseline using relevel.
dds <- DESeqDataSetFromMatrix(countData=subset_neutral,
colData=metadata_neutral,
design= ~genotype, tidy = FALSE)
dds$genotype<- relevel(dds$genotype, "WT")
dds <- DESeq(dds)
results_neutral<-results(dds_neutral)
Now the issue I have. This code seems to work perfectly well for my first subset of the original dataset. In the results file all values are calculated correctly and I get nice visualization of the differentially expressed genes in the volcanoplot. However, when I do this same analysis on my other treatments something seems to go wrong. I have very little differentially expressed genes. The dataframes I get back have a lot of genes for which no pvalue/adjusted pvalue has been calculated. Now this is the case still after I have removed all the genes with low read counts from my analysis.
The low number of DEGs could of course be actual biology, but to be honest based on the treatment and genotype, as well as a lot of other data we have, that is the last thing I expect. I just want to make sure that I am not doing anything wrong with the analysis itself.