Dear all, and Aaron,

I am using the following R code for the paired design analysis of a 6*2 RNA-seq samples : please could you let me know if there are any changes/modifications you would suggest . thanks !

library("limma")
library("edgeR")

### reading the expression dataset

eset <- read.delim("mm10.strand.both.count.exons.condense.genes.NOADJ",row.names="Symbol")

### reading the library sizes

libsize <- c(28182007, 28599704, 26386209, 28069058, 26361374, 27108450, 27097918, 12155292, 27124214, 26774125, 26241068, 26296757)

### setting up the groups and the subjects

group <- factor(c("csc","csc","csc","csc","csc","csc","non","non","non","non","non","non"))
subject <- factor(c(1,2,3,4,5,6,1,2,3,4,5,6))

### setting up the design and the contrast matrix

design <- model.matrix(~0+group+subject)
contrast.matrix <- makeContrasts(groupcsc-groupnon, levels=design)

### filtering the genes based on CPM :

y <- DGEList(counts=eset,group=group)

keep <- rowSums(cpm(y, lib.size=libsize)>1) >= 3
y <- y[keep,]

y$samples$lib.size <- colSums(y$counts)

### computing the normalization factors :

y <- calcNormFactors(y)

### using the VOOM transformation :

v <- voom(y,design,plot=TRUE)

### doing the LINEAR FIT in LIMMA :

fit <- lmFit(v, design)

fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)

### obtaining and writing the results :

results <- topTable(fit2, coef=1, adjust="fdr", number=24073)

write.table(results, file="results", sep="\t", eol="\n", row.names=TRUE, col.names=TRUE)

This is quite a standard analysis and I guess you already know that your code is basically correct.

Considering that the library sizes are almost equal, I would personally run a slightly simpler, cleaner analysis using limma-trend:

library("edgeR")
eset <- read.delim("mm10.strand.both.count.exons.condense.genes.NOADJ",row.names="Symbol")
group <- factor(c("csc","csc","csc","csc","csc","csc","non","non","non","non","non","non"))
group <- relevel(group, ref="non")
subject <- factor(c(1,2,3,4,5,6,1,2,3,4,5,6))
design <- model.matrix(~group+subject)
y <- DGEList(counts=eset, group=group)
keep <- rowSums(cpm(y) > 0.5) >= 6
y <- y[keep,, keep.lib.sizes=FALSE]
y <- calcNormFactors(y)
logCPM <- cpm(y, log=TRUE, prior.count=3)
fit <- lmFit(logCPM, design)
fit <- eBayes(fit, trend=TRUE, robust=TRUE)
plotSA(fit)
results <- topTable(fit, coef=2, number=Inf)

Thanks a lot Gordon for confirmation. A question though, if I may add : although I have looked through several webpages,

still not very sure about the difference between :

design <- model.matrix(~group+subject)

and

design <- model.matrix(~0+group+subject)

in terms of the practical output, I have not noticed any difference, I am asking just for my knowledge. Thanks again !