Hi Aaron and Ryan,

I am getting an error while designing a model (design <- model.matrix(~individuals). This code is not working.

Could you please go through the following code?

Thanks for your help.

R  code: 

x<-read.delim("test.txt",header=T,sep="\t")
y<-DGEList(counts = x[,2:13]) #using only 12 individuals for test 
y$samples
#Filter lowly expressed genes
keep<-rowSums(cpm(y)>1)>=2
y<-y[keep,keep.lib.sizes=FALSE]
head(y$counts)
#individuals1<-factor(rep(1:12,each=3))
individuals<-factor(y$counts)
design <- model.matrix(~individuals) # it's not working
#design <- model.matrix(object = ~y$counts)
y<-estimateGLMCommonDisp(y, design)
y<-estimateGLMTrendedDisp(y,design)
y<-estimateGLMTagwiseDisp(y.design)
fit<-glmFit(y,design)
res <- glmLRT(fit, coef=2:ncol(design))

Hm... well, if you can re-use particular individuals for new experiments, perhaps individual-specific inferences might be useful. Still seems a bit iffy, though. Anyway, if you want to identify genes that are unique to an individual, there's a couple of things you can do. The least specific approach would be:

plant <- factor(rep(1:50, each=3))
design <- model.matrix(~0 + plant)
# ... edgeR stuff ...
con <- makeContrasts(plant1 - (plant2 + ... + plant50)/49, levels=design)
res <- glmLRT(fit, contrast=con)

... which would identify genes that are up/down-regulated in plant 1 compared to the average of all other plants. This should give you a fair number of DE genes, but many of these may not be specific to plant 1 in the strict sense - for example, if a gene is expressed in both plant 1 and 50, you might still detect it as being DE if all other plants have lower expression.

The most specific approach would be something like this:

con <- makeContrasts(plant1 - plant2, ..., plant1 - plant50, levels=design)
all.pvalues <- list()
all.signs <- list()
for (x in seq_len(ncol(con))) {
    res <- glmLRT(fit, contrast=con[,x])
    all.pvalues[[x]] <- res$table$PValue
    all.signs[[x]] <- sign(res$table$logFC)
}
final.pvalue <- do.call(pmax, all.pvalues)
final.fdr <- p.adjust(final.pvalue, method="BH")
all.signs <- do.call(cbind, all.signs)
same.sign <-  rowSums(all.signs[,1]!=all.signs)==0L
keep <- which(final.fdr <= 0.05 & same.sign)

... which will identify genes that are DE in the same direction between plant 1 and each other plant. This will probably give no genes, as it's very conservative - lack of evidence for DE between any two plants will cause the gene to be discarded.

A compromise might be to do an ANOVA:

res <- glmLRT(fit, contrast=con) # using the same contrast above.
current.fdr <- p.adjust(res$table$PValue, method="BH")
all.signs <- sign(res$table[,grep("logFC", colnames(res$table)),drop=FALSE])
same.sign <-  rowSums(all.signs[,1]!=all.signs)==0L
keep <- which(current.fdr <= 0.05 & same.sign)

... which will identify all genes that have any DE between plant 1 and any other plant, and where plant 1's expression is either higher or lower than that of all other plants.

I should add that you'll have to set up different contrasts for each plant. This can be done without re-writing the long makeContrasts call every time. For the first approach, you can just do:

con <- rep(-1/49, 50)
con[x] <- 1 # where 'x' is the plant of interest.

... while for the second and third approaches, you can do:

con <- -diag(50)
con[x,] <- 1
con <- con[,-x]