paired test for differential expression on RNAseq data

Question

[minoo]

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 outcome had weird 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?

Answer 1

[pilargmarch]

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

dds$condition <- relevel(dds$condition, ref = "ctr")
dds <- estimateSizeFactors(dds)
dds <- DESeq(dds)
results(dds, name="condition_tre_vs_ctr")

If you do this, then a positive log-fold change will mean the gene is upregulated in the treatment sample compared to the control one. Check that you're specifying the right contrast when you get the results (I believe the code above is the one you need).

Answer 2

[ATpoint]

These packages are well established and work. Weird results are usually a user problem. Please show how you extracted counts and the results from deseq2. The edgeR code is wrong, you're mixing limma and edgeR commands. Follow the manual for guidance.