Hi everyone I am have RNA seq data that I am trying to analyze using EdgeR. Here is the codes I used to analyze the data I just want to make sure that the codes are right and that I can rely on the output for the analysis. The summary of the experiments: I have 4 replicates of infected cells with Cryptosporidium and 4 replicates of control (uninfected)

rawdata <- read.csv("celldata.csv", check.names=FALSE, sep=",", header=T, stringsAsFactors=FALSE)
library(edgeR)
library(limma)
rownames( counts ) <- raw.data[ , 1 ]
> y <- DGEList(counts=rawdata[,2:9], genes=rawdata[,1])

filtering 

> o <- order(rowSums(y$counts), decreasing=TRUE)
> y <- y[o,]
> d <- duplicated(y$genes)
> y <- y[!d,]
> nrow(y)
> y$samples$lib.size <- colSums(y$counts)
>  rownames(y$counts) <- rownames(y$genes) <- y$genes$EntrezGene
> y$genes$EntrezGene <- NULL
>  y <- calcNormFactors(y)
>  y$samples
> plotMDS(y)

Design matrix

cells <- factor(c(1,2,3,4,5,6,7,8))
> status<- factor(c("i","i","i","i","c","c","c","c"))
> data.frame(Sample=colnames(y),cells,status)
> design <- model.matrix(~cells)
>  rownames(design) <- colnames(y)
> Estimating the dispersion 
>y d1$common.dispersion
> y <- estimateTrendedDisp(y)
> y <- estimateTrendedDisp(y)
>y <- estimateCommonDisp(y)
> y<- estimateCommonDisp(y)
> y <- estimateTrendedDisp(y)
> y<- estimateTagwiseDisp(y)
> plotBCV(y)
> fit <- glmFit(y, design)
>  lrt <- glmLRT(fit)
> topTags(lrt)

Really appreciate any feedback

There seem to be multiple problems here. Some are harmless, others are less so:

• Your "filtering" step just removes duplicated genes. The real purpose of filtering is to remove low abundance genes, which you have not done. This will lead to suboptimal performance, as the overabundance of genes with near-zero counts will skew the mean-dispersion trend fit; reduce power by increasing the number of tests (and severity of the multiple testing correction) without increasing the number of rejections, as the counts are too low to reject the null; and increase computational time.
• Your design matrix is constructed using the cells factor. This is almost certainly wrong, as you're basically saying that there are no replicates in your experiment (as each sample corresponds to a different level of cells). You should probably be constructing your design matrix from status.
• I don't know why you're mixing classic edgeR (estimateTrendedDisp, estimateCommonDisp) with GLM edgeR (glmFit, glmLRT). While this doesn't cause too many problems for your two-group setup, you should generally stick with one or the other. Also, the line y d1$common.dispersion can't be syntactically valid.

Other things are probably okay but can be improved. In particular, you don't need to reorder by abundance; recalculation of the library sizes can be done with keep.lib.sizes=FALSE during subsetting; and there's no need to repeat the dispersion estimation calls.

Many thanks for your quick response. I realized there is a problem with filtering since I get the same number before and after filtering. This means filtering did not make any changes or remove any reads. How do you recommend to change this codes? Regarding cells factor I understand the reason why it is wrong but I don't know how to fix it. 1,2,3 nd four belong to the infected group and 5,6,7,8 belong to the control group. Thanks, 

Thanks for providing useful information. Here is my finalized codes. I got differentially expressed genes, GO terms and metabolic pathways based on these codes. How can I conform that these are correct findings?

rawdata <- read.csv("celldata.csv", check.names=FALSE, sep=",", header=T, stringsAsFactors=FALSE)

> y <- DGEList(counts=rawdata[,2:9], genes=rawdata[,1])

keep <- aveLogCPM(y) > 2

y <- y[keep,,keep.lib.sizes=FALSE]

y <- calcNormFactors(y)

 design <- model.matrix(~infected)

y <- estimateDisp(y, design)

fit <- glmQLFit(y, design, robust=TRUE)

res <- glmQLFTest(fit, coef='infected1')

topTags(res)

keg <- kegga(qlf)

topKEGG(keg)

go <- goana(qlf)

> topGO(go)