Hello everyone. I'm performing a single-cell RNA-seq on a MusMusculus dataset following the DESeq2 instructions in https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html

In my experiment there are two variables of interest: disease/no-disease and receptor1/receptor2, and what i'm attempting to do is to have a comparison disease_vs_nodisease for each of the two receptors.

In this code i'm running, i am wondering if using the contrast for the receptors is the right solution, or if i must consider other ways to perform the analysis.

setwd("Desktop/Normalize")

# we create the beginning dataset
countRawData <- read.table("allSamples_rawData.txt",header = TRUE, sep = "\t")
rownames(countRawData) <- countRawData[,1]
countRawData <- countRawData[,-1] #the first column in original file is set as rownames 
coldata <- read.table("Metadata.txt", header = TRUE, sep ="\t")
rownames(coldata) <- coldata[,1]
coldata <- coldata[,-1]
coldata$Receptor <- as.factor(coldata$Receptor)
coldata$Mouse <- as.factor(coldata$Mouse)  #we need all factors, use str(coldata) to check 
coldata$PolyQ <- relevel(coldata$PolyQ, ref = "Q20")  # we set Q20 as "control" (non-huntington condition)
# make sampleIDs coherent in matrix and metadata. 
colnames(countRawData) <- gsub("\\.","-", colnames(countRawData))
#SampleIDs ordered same way in the two datasets
countRawData <- countRawData[,rownames(coldata)]

# filter out low reads genes
countRawData <- countRawData[rowSums(countRawData) >= 5,]

#We use ZINB-Wave to model our data, to manage zero-inflation
#but, first, make a summarizedExperiment object 
library("zinbwave")
rawCountsSE <- SummarizedExperiment(assays = list(counts=as.matrix(countRawData)), colData=coldata)

#setting as recommended by the ZINB-WaVE integration paper
zinb_counts <- zinbwave(rawCountsSE, K=0, BPPARAM=SerialParam(), epsilon=1e12)

#now, Deseq2, following the workflow @ https://github.com/mikelove/zinbwave-deseq2/blob/master/zinbwave-deseq2.knit.md
library("DESeq2")
dds_counts<- DESeqDataSet(zinb_counts, design=~PolyQ+Receptor)
dds_counts <- DESeq(dds_counts, test="LRT", reduced=~1, sfType="poscounts", minmu=1e-6, minRep=Inf)

#extract results using IHW for more power 
library("IHW")
DESResults <- results(dds_counts, contrast = c("Receptor", "1", "2"), filterFun = ihw)



> DESResults
log2 fold change (MLE): Receptor 1 vs 2 
LRT p-value: '~ PolyQ + Receptor' vs '~ 1' 
DataFrame with 23970 rows and 7 columns
                            baseMean      log2FoldChange             lfcSE
                           <numeric>           <numeric>         <numeric>
ENSMUSG00000051951  1.20874783783011  -0.486272008164474 0.264799182097863
ENSMUSG00000102851  0.17692628048305   0.184089986262215  1.03275481212339
ENSMUSG00000103377  0.45483105797926   0.364930394500295 0.464969162738308
ENSMUSG00000104017 0.256058234805924  -0.368233429339006 0.660452353614404
ENSMUSG00000103025 0.245374687264191  -0.110470841313906 0.626788064771599
...                              ...                 ...               ...
ENSMUSG00000064368  2.74872956747472 -0.0637947150790935  0.17126290563671
ENSMUSG00000064369 0.237388941203093    1.03934379136416 0.532851024735219
ENSMUSG00000064370  3.35773837543451 0.00781897550751081 0.149162326168115
ENSMUSG00000064371 0.200180852346864   0.811713031483214 0.792373303960726
ENSMUSG00000064372  1.76942098675239  0.0567967034292243 0.209702611492413
                                 stat               pvalue               padj
                            <numeric>            <numeric>          <numeric>
ENSMUSG00000051951   3.49645563810861    0.174082175362577  0.892704290739862
ENSMUSG00000102851 0.0638155400003413    0.968595911560463                  1
ENSMUSG00000103377   1.66785646184704    0.434339743967727                  1
ENSMUSG00000104017  0.631714784159911    0.729163427233201                  1
ENSMUSG00000103025   1.40320272869138    0.495790726165263                  1
...                               ...                  ...                ...
ENSMUSG00000064368  0.770642035229798     0.68023223462155  0.953292547184548
ENSMUSG00000064369   14.2439914976536 0.000807154265422844  0.163707387056423
ENSMUSG00000064370   11.8457607335204  0.00267747694383133 0.0957254477349814
ENSMUSG00000064371   5.91910585474172   0.0518420891716626                  1
ENSMUSG00000064372  0.995816817125558    0.607800601683803                  1
                               weight
                            <numeric>
ENSMUSG00000051951   1.20666010946513
ENSMUSG00000102851  0.637350033274622
ENSMUSG00000103377   0.45922123873472
ENSMUSG00000104017  0.514183842111013
ENSMUSG00000103025  0.375903805678487
...                               ...
ENSMUSG00000064368   2.86177930272658
ENSMUSG00000064369  0.514183842111013
ENSMUSG00000064370   4.11903005355909
ENSMUSG00000064371 0.0482764678116727
ENSMUSG00000064372 0.0167393501223644