dge_analysis.r

#!/usr/bin/env Rscript

library(edgeR)
RawCounts <- read.delim("count_matrix.txt", row.names = "gene_id")
group <- read.table("SampleGroups.txt", header=TRUE, sep="\t", row.names=1)
dgecomplete <- DGEList(RawCounts, group = group$Condition)
logcpm <- cpm(dgecomplete, log=TRUE)
filtData <- filterByExpr(dgecomplete)
dgecomplete <- dgecomplete[filtData, keep.lib.sizes=FALSE]
dgecomplete <- calcNormFactors(dgecomplete)
dgecomplete <- estimateDisp(y = dgecomplete)
fit <- glmQLFit(y=dgecomplete)
qlf <- glmQLFTest(fit, coef = 2)
diff_results <- topTags(qlf, n=Inf)
write.csv(diff_results, file="edgeR_diff_genes.csv")
write.csv(as.data.frame(logcpm), file="edgeR_normcounts.csv")