I would like some guidance when it comes to using DESeq2 to look for DEGs in a bulkRNA seq experiment. I have at my disposal samples that are enriched for endothelial cells from 2 tissues and at 6 developmental stages. I have 3 samples from each condition for a total of 36 samples. An analysis of the endothelial purity of the samples shows that samples from tissue A have more contamination from mural cells that tissue B samples. Uncorrected, a list of DEGs comparing A to B gives many (mostly?) hits that are mural cell genes. Therefore I would like to add a covariate for each sample with the estimated proportion of mural cells. (Here I am hesitant, as I think this has been discouraged in related threads, but I am not entirely sure).
Because I have so few samples as each timepoint, I am reluctant to only use input from a single stage in my design.
So my code now looks something like this:
dds <- DESeqDataSetFromMatrix(countData = cts,
colData = coldata,
design = ~ Mural + Tissue + Stage + Tissue:Stage)
I would like to identify DEGs when comparing Tissue A vs Tissue B at a specific stage. I think this gets down to specifying contrasts, but I am unsure about how to proceed.
It seems similar to the vignette describing a time course analysis but I don't really want to compare my developmental stages back to a reference stage. I think the code below would give me differences at a given stage but after correcting for the reference level of the stage covariate.
ddsTC <- DESeq(dds, test="LRT", reduced = ~ Mural + Tissue + Stage)
res <- results(ddsTC, name="TissueA.Stage1", test="Wald")
Is there a suggestion for how to handle this type of analysis within DESeq2 or alternatively a suggestion for a better solution?
Thank you for any guidance or advice.
-ross
