Question: R analysis for read counts using edgeR
gravatar for Liftedkris
5 months ago by
Liftedkris10 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 • 357 views
ADD COMMENTlink modified 5 months ago by h.mon27k • written 5 months ago by Liftedkris10

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 5 months ago by ATpoint23k

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 5 months ago by Friederike5.1k
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: 531 users visited in the last hour