Hi,

I am using EdgeR for differential expression. The experiment includes five individuals who receive the vaccine in 4 different time points (pre, 2hour post, 24hour post, 14-day post). I want to find DE genes between pre and post samples, but unfortunately, i couldn't find a way to design my matrix in a way that i get results. i did PlotMD as you saw in the image sample 5-8 show separate batch in contrast to others. I intend to make the identification of DE genes using a log2 fold change and likelihood ratios (LR) Test in edgeR and significantly expressed genes had an FDR adjusted P-value of < 5%. Here is the script i use for my analysis:

library(edgeR)
setwd("~/Desktop/counts/low")
x_low<-read.csv("counts.txt")
time<-factor(c("pre","2hr-post","24hr-post","14d-post",
           "pre","2hr-post","24hr-post","14d-post",
           "pre","2hr-post","24hr-post","14d-post",
           "pre","2hr-post","24hr-post","14d-post",
           "pre","2hr-post","24hr-post","14d-post"))
time <- relevel(time, ref="pre")
y_low<-DGEList(counts=x_low[,2:21], genes = x_low[,1], group = time) #read counts.csv and make table  
keep <- filterByExpr(y_low)
y_low <-y_low[keep, keep.lib.sizes=FALSE] 
y_low <- calcNormFactors(y_low)
plotMDS(y_low)
data.frame(sample=colnames(y_low),time) #data frame
      sample      time
1   Sample_1       pre
2   Sample_2  2hr-post
3   Sample_3 24hr-post
4   Sample_4  14d-post
5   Sample_5       pre
6   Sample_6  2hr-post
7   Sample_7 24hr-post
8   Sample_8  14d-post
9   Sample_9       pre
10 Sample_10  2hr-post
11 Sample_11 24hr-post
12 Sample_12  14d-post
13 Sample_13       pre
14 Sample_14  2hr-post
15 Sample_15 24hr-post
16 Sample_16  14d-post
17 Sample_17       pre
18 Sample_18  2hr-post
19 Sample_19 24hr-post
20 Sample_20  14d-post
design<-model.matrix(~time) 
y_low<-estimateDisp(y_low,design)
fit<-glmFit(y_low,design)
lrt<-glmLRT(fit)
deg = topTags(lrt,p.value = 0.05)$table

Thanks in advance

You have repeated measures on five different subjects, so you should block on subject. See 3.4.1 in the edgeR User's Guide.

Also, you should provide the design matrix to filterByExpr, otherwise you will filter as if you have only one group.

When you use a function in R (or most software, so far as I know), you can use named arguments (like counts = x_low[,2:21]) or you can use positional arguments, where R infers what a given argument is, by its position in the function call. You did this:

y_low<-DGEList(counts=x_low[,2:21], genes = x_low[,1],time)

Where, by mistake or intention you have two named arguments and one positional argument (the third argument, time). An interesting question in this instance is what will happen with the positional argument, since it follows two named arguments? I didn't know what would happen, so I made a stupid test:

> z <- DGEList(counts = matrix(rnbinom(100, 3, 0.5), 10), genes = data.frame(genes = LETTERS[1:10]), letters[1:10])
> z
An object of class "DGEList"
$counts
   Sample1 Sample2 Sample3 Sample4 Sample5 Sample6 Sample7 Sample8 Sample9
1        5       0       1       2       2       6       0       0       2
2        4       7       6       2       5       7       2       3       3
3        2       3       4       2      10       3       4       7       4
4        0       6       1       1       3       2       2       3       6
5        4       2       1       2       2       6       3       5       1
6        8       4       7       5       2       3       4       2       7
7        2       6       5       4       7       5       0       9       4
8        1       0       1       2       5       4       0       0       5
9        6       2       8       3       9       3       2       4       2
10       8       4       2       1       1       6       3       3       2
   Sample10
1         4
2         2
3         1
4         1
5         2
6         1
7         3
8         3
9         3
10        7

$samples
         group lib.size norm.factors
Sample1      1        a            1
Sample2      1        b            1
Sample3      1        c            1
Sample4      1        d            1
Sample5      1        e            1
Sample6      1        f            1
Sample7      1        g            1
Sample8      1        h            1
Sample9      1        i            1
Sample10     1        j            1

$genes
   genes
1      A
2      B
3      C
4      D
5      E
6      F
7      G
8      H
9      I
10     J

With a mix of positional and named arguments, it appears that R will infer the position of the positional argument based on the position of the named arguments (the arguments for DGEList are counts, lib.size, norm.factors, samples, group, genes, remove.zeros). Since you used names for the first and sixth arguments, the positional argument is then inferred to be the lib.size. In other words, you already 'took' the first and sixth position, so the first available 'empty' position is the second one, which is lib.size.

In other words, you told DGEList that your library sizes were your time object, which is a factor! And it is unsurprising that when the offsets are computed and the underlying function expects a numeric vector but gets a factor, you get an error:

> estimateDisp(z, model.matrix(~gl(2,5)))
Error in .compressOffsets(y, lib.size = lib.size, offset = offset) : 
  offsets must be finite values