Question: R analysis for read counts using edgeR
gravatar for Liftedkris
21 months ago by
Liftedkris20 wrote:

Hi fellows,

I performed RNAseq analysis of human neutrophils infected by Aspergillus fumigatus. Everything seemed to work fine but during the differential expression analysis using EdgeR, I used the standard p.value cut off of 0.05 but surprisingly no genes were found for the conditions at that p. Value .

My question is, what could be the problem?

I have 24 samples, in three groups, each group consisting of 8 samples 4 each for control and test i did alignment with Hisat2 and sorted with samtools and used HTseq to generate read counts. at this point everything seemed to work fine.


data <- setwd("C:/Users/HP Gamer/Desktop/RNA_seq_read_counts")
group <- factor(c(rep("host.30hpi_control", 4), rep("host.30hpi_infected", 4), 
rep("host.90hpi_control", 4), rep("host.90hpi_infected", 4), rep("host.180hpi_control", 4), rep("host.180hpi_infected", 4) )) 

library(edgeR) <- readDGE(list.files(pattern = ".count"), data, columns = c(1,2))$counts <-$counts[1:(nrow($counts)-5), ]$counts <-$counts[rowSums($counts > 3) > 2, ]     
head($counts, 20) 
dim($counts) <- calcNormFactors( 
design <- model.matrix(~0 + group)     
rownames(design) <- colnames($counts)
contrasts <- makeContrasts("Host_30hpi" = grouphost.30hpi_infected - grouphost.30hpi_control, "Host_90hpi" = grouphost.90hpi_infected - grouphost.90hpi_control, "Host_180hpi" = grouphost.180hpi_infected - grouphost.180hpi_control, levels = design)


y <- voomWithQualityWeights(counts =, design = design, plot = TRUE) 

plot.colors <- c(rep("blue", 4), rep("red", 4), rep("orange", 4), rep("black", 4)) 
plotMDS(, main = "MDS Plot for Count Data", labels = colnames($counts), col = plot.colors, cex = 0.9, xlim=c(-2,5)) <- t(cpm( 
dist <- dist( 
plot(hclust(dist), main="Hierarchical Clustering Dendrogram") 

fit <- lmFit(y, design) 
fit <-, contrasts) 
fit <- eBayes(fit) 
top_30hpi <- topTable(fit, coef = "Host_30hpi", adjust = "fdr", number = "Inf", p.value = 0.05, = "P") 
top_90hpi <- topTable(fit, coef = "Host_90hpi", adjust = "fdr", number = "Inf", p.value = 0.05, = "P") 
top_180hpi <- topTable(fit, coef = "Host_180hpi", adjust = "fdr", number = "Inf", p.value = 0.05, = "P") 
write.table(top_30hpi, file = "Host_DEG_annotated30.csv", sep = ",", col.names = NA)
write.table(top_90hpi, file = "Host_DEG_annotated90.csv", sep = ",", col.names = NA)
write.table(top_180hpi, file = "Host_DEG_annotated180.csv", sep = ",", col.names = NA)
limma edger rna-seq • 1.1k views
ADD COMMENTlink modified 21 months ago by h.mon32k • written 21 months ago by Liftedkris20

Do you realize that you are not using edgeR at all for differential analysis but limma? I am not adept with limma so I cannot contribute too much but why don't you simply follow exactly the edgeR manual and then see what comes out?

ADD REPLYlink written 21 months ago by ATpoint44k

Just a couple of comments:

  • lmFit is a function from limma, which is presumably why you chose to use the voom-normalized values for the DE analysis
  • what is this line supposed to do:$counts <-$counts[1:(nrow($counts)-5), ], i.e. why are you removing the last five rows?
  • what do the MDS and Hclust plots actually look like?
  • I'm not sure what your group and design objects look like, do you mind sharing them?
ADD REPLYlink written 21 months ago by Friederike6.7k
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 934 users visited in the last hour