Basic RNA Differential Expression in R

Question

I have two matrices, one for individuals before treatment and one for the same individuals after treatment. Both matrices are raw read counts of RNA expression.

        Treated_1  Treated_2  Treated_3  Treated_4  Treated_5
RNA_1   105        283        64         155        51
RNA_2   359        27         47         348        84
RNA_3   99         10         89         345        77
RNA_4   48         100        77         74         83

        Untreated_1  Untreated_2  Untreated_3  Untreated_4  Untreated_5
RNA_1   100          200          50           130          33
RNA_2   200          10           30           300          75
RNA_3   90           10           60           320          60
RNA_4   66           50           33           60           55

How can I tell which RNA's are significantly differentially expressed in R using edgeR or DESeq2? I'm thinking first normalize using TMM then compare the mean of the normalized values for treated vs. untreated?

Answer 1

See here for a Differential Expression guide which discusses carrying out differential expression using DESeq2:

https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html

It starts with a standard "quick start" workflow, and goes into quite a lot of detail explaining all the different things that can be done:

dds <- DESeqDataSetFromMatrix(countData = cts,
                              colData = coldata,
                              design= ~ batch + condition)
dds <- DESeq(dds)
resultsNames(dds) # lists the coefficients
res <- results(dds, name="condition_trt_vs_untrt")

Answer 2

I'm thinking first normalize using TMM then compare the mean of the normalized values for treated vs. untreated?

That is a bit too simple. You need to analyse your experiment as a paired comparison. You need to combine the treated and untreated count matrices into one and form factors for individual and treatment:

Counts <- cbind(Counts.treated, Counts.untreated)
Individual <- gl(5,1,length=10)
Treatment <- gl(2,5,length=10,labels=c("Treated","Untreated"))
Treatment <- relevel(Treatment, ref="Untreated")
design <- model.matrix(~Individual+Treatment)

Then follow a standard workflow: https://bioconductor.org/packages/release/workflows/vignettes/RnaSeqGeneEdgeRQL/inst/doc/edgeRQL.html

library(edgeR)
y <- DGEList(counts=Counts)
y <- calcNormFactors(y)
y <- estimateDisp(y,design)
fit <- glmQLFit(y,design)
test <- glmQLFTest(fit)
topTags(test)