Hello, I'm very sorry for posting again a question regarding RNA-seq with no replicates. I've read most of the posts and still don't know what's the best way to proceed. My experimental design was: two industrial fermenters, one as control, the other as the treatment and we took samples every 24 hours for a total of 4 samples. Due to a fail on the fermenter I have replicates only for the first time points. So My first approach was to use the last two time points as one (replicates). I've used edge R, and make contrast for each pairwise comparison I want to test. My BCV is really high (2.06) and I've used glmQLFit. So I got the FC for each of the contrasts. But my main problem is that the replicates are very very different from each other so my supervisor suggest to redo the analysis by individual fermenter run (which means without replicates). Does this seems sensible? If so, I've been checking edgeR on how to proceed but as soon as I get to calculate the dispersion it fails, I get this error: "Warning message:

In estimateGLMCommonDisp.default(y = y$counts, design = design,  :
  No residual df: setting dispersion to NA".

My code:

count_cols <- c("Gene","GLU24H","GLU48H","GLU72H","GLU.FA24H","GLU.FA48H","GLU.FA72H")
x<-read.delim('~/Desktop/PhDProject/3rd_Year_RESULTS/Counts/counts_run1.csv',skip=0, sep=",", check.names=FALSE, colClasses='character', na.strings=c())
x[,count_cols] <- apply(x[,count_cols], 2, function(v) as.numeric(v))     
counts <- x[, count_cols]
keepMin <- apply(counts, 1, max) >= 1
keepCpm <- rowSums(cpm(counts)> 1) >= 1         
keep <- keepMin & keepCpm
x <- x[keep,]
counts <- counts[keep,]
counts
group<-factor(c("GLU24H","GLU48H","GLU72H","GLU.FA24H","GLU.FA48H","GLU.FA72H"))
y <- DGEList(counts=counts, group = group)
y <- calcNormFactors(y, method="TMM")
y$samples
genes<-as.numeric(x$GENE) 
typeof(genes)
PGRA_rpkm<-rpkm(counts, log=FALSE, gene.length = genes)
PGRA_rpkm2<-cbind(x$Gene,PGRA_rpkm)
plotMD(cpm(PGRA_rpkm, log=TRUE), column=6) > abline(h=0, col="red", lty=2, lwd=2)

design <- model.matrix(~0+group)
y <- estimateGLMCommonDisp(y,design)

Any help will be really useful as I'm currently stuck at this. 

Thanks in advance!!

You don't need replicates for every time point to estimate the dispersion. If you have replicates at one time point, edgeR will use those to estimate the dispersion, i.e., you only need non-zero residual degrees of freedom in the entire model. For each gene, we model the counts using a negative binomial distribution that has the same dispersion for every sample, so it doesn't matter where the replication occurs.

But anyway, I'll move onto the data set you've presented. As it's currently set up, you don't have any residual d.f. because you have 6 coefficients to estimate from 6 samples. This means that there's no information left to estimate the dispersion after fitting of the GLM, resulting in the error from estimateGLMCommonDisp. There are two solutions; the first is to follow the advice in Section 2.11 of the user's guide, and just fit the GLM with a "sensible" dispersion value for the purpose of computing log-fold changes for your desired contrasts:

fit <- glmFit(d, design, dispersion=0.05)
res <- glmLRT(fit, contrast=makeContrasts(groupGLU.FA24H-groupGLU24H, levels=design))

You could also sacrifice some residual d.f. in your model to obtain that sensible estimate:

time <- factor(rep(c("24", "48", "72"), 2))
fermenter <- factor(rep(c("GLU", "GLU.FA"), each=3))
design0 <- model.matrix(~time+fermenter)
# ... estimate dispersion with design0.

I wouldn't put much faith in the p-values, though. The dispersion estimates will be somewhat inflated as design0 assumes that the two fermenters are mostly similar in their effect on expression over time. More importantly, testing with dispersion estimates obtained with a different model will only work with glmLRT, not glmQLFTest; and you won't be able to model gene-specific dispersion values, which will reduce the accuracy of your results. This is the cost of not having replicates - see point 3 in Section 2.11.

The second strategy is to assume linearity in the responses, which allows you to pull out some residual d.f. from your design:

time.as.numeric <- rep(c(24, 48, 72), 2)
designL <- model.matrix(~0+fermenter+time.as.numeric:fermenter)
# ... run through edgeR with designL

You can then test for differences in the gradient of the time response, which tells you if there's a difference in the time effect. Testing for differences in the intercept also tell you if there's already a difference at zero hours. If you had more samples, you could weaken the assumption of linearity by fitting a fermenter-specific spline instead, which will model a smooth curve with respect to time.

Hello Aaron, thank you very much!!! Ok following your nice explanation I'll proceed with the second strategy. Although I have some questions (I'm new to R and edgeR so I've learned basically using this package and looking up on internet how to do stuff, this said my understanding of how to code/interpret some things is not very good). Hence, my next questions are,

1)If I proceed with edgeR using designL, do you recommend then to use glmLRT in order to get the FC values for each of the contrast I'm looking for? (which basically are a pairwise comparison between fermenters and time).

2)I do have more samples but come from another a fermentation ran 6 months after. Biologically speaking I spect some difference... my question is, statistically speaking do you think I can use those samples on the same model to get a better fit or should I work with them separately? this has been one of my starting questions but I ended up in circles. Also I use them all together (fermentation 1 and fermentation 2) to calculate the rpkms in edgeR, does this seems sensible to you? 

Stefany