paired test for differential expression on RNAseq data

I need to make a paired wise differential expression test for a metadata like below:

colData<- DataFrame(row.names = c(       "1",
                                         "2",
                                         "3",
                                         "4",
                                         "5",
                                         "6"),
                    Patient=c("b1","c1","d1","b1","c1","d1"),
                    condition=c("ctr","ctr","ctr","tre","tre","tre"))

DataFrame with 6 rows and 2 columns
      Patient   condition
  <character> <character>
1          b1         ctr
2          c1         ctr
3          d1         ctr
4          b1         tre
5          c1         tre
6          d1         tre

I have two conditions and three patients and I need to do this test paired wise betweent each patient.

I first tried to use deseq2 using this command:

dds <- DESeqDataSetFromMatrix(countData = cstm,
                         colData = colData,design = ~ condition
                          +Patient 
                         )

But the outome had weired result meaning some downregulated genes based on their counts had highly upregulated logfoldchanges.
Then I tried to use edger but keep getting errors:

library(edgeR)

meta<-as.data.frame(colData)
meta$Patient<-as.factor(meta$Patient)

design = model.matrix(~ 0+Patient + condition, data = meta)  

# After the design, create a contrast of the pairs u want to compare.
contr.matrix <- makeContrasts(
  pb1vsc1 = Patientb1- Patientc1,
  pd1vsc1 = Patientd1- Patientc1,
  pb1vsd1 = Patientb1- Patientd1,
  levels = colnames(design))

y <- DGEList(counts = cstm)

# Dispersion  
y <- estimateDisp(y, design)

# Fit
fit <- glmQLFit(y, design)

# contrast fit
vfit <- contrasts.fit(fit, contrasts=contr.matrix)
# Empirical Bayes Statistics for Differential Expression
efit <- eBayes(vfit)
# summary
summary(decideTests(efit))

# log fc cutoff
tfit <- treat(vfit, lfc=1)
dt <- decideTests(tfit)
summary(dt)
Error in contrasts.fit(fit, contrasts = contr.matrix) :
fit must contain stdev.unscaled component

Any idea how to proceed this?

Read more here: Source link