Hi,

I am trying to conduct a DE analysis comparing my wildtype single cell data with 4 other genotypes. If there is a change in gene expression in at least 1 of the 4 genotypes compared to WT than that gene is considered DE. We ran the following code taking into account time course however the only output is one padj/p value for each gene, which I assume is the average across all 5 genotypes. However, I want a padj/p value for each gene for each genotypes so I can compare my WT to the others. Any suggestions on how to do this?

Code should be placed in three backticks as shown below

time_dynam_deseq <- function(pt) {
  metadata <- data.frame(timepoint = pt$timepoint, strain = pt$strain, pseudotime = pt$slingPseudotime_1)
  rownames(metadata) <- rownames(colData(pt))

  dds <- DESeqDataSetFromMatrix(countData = as.matrix(counts(pt) + 1),
                                colData = metadata,
                                design = ~strain* pseudotime)

  dds <- DESeq(dds, test = "Wald")
  dds
}

dds_strains <- unique(dds$strain)
isolated_strains <- list()
for (strain in dds_strains) {
  isolated_strains[[strain]] <- subset(dds, strain == strain)
}

#For example, both strain results yield the same values
dmd3_res <- results(dds, contrast=c("strain", "dmd-3(-)", "WT"))
wt_res <- results(dds, contrast=c("strain", "WT", "dmd-3(-)"))

# include your problematic code here with any corresponding output 
# please also include the results of running the following in an R session 

sessionInfo( )

The two contrasts at the bottom are identical but for the sign of the logFC, so you should expect the same results (except for flipped sign of the logFC).