count_analysis.r

library(edgeR)
RawCounts <- read.delim("count_matrix_cleaned.txt", row.names = "ENSEMBL_GeneID")
group <- read.table("sample_groups.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")

# Functional analysis
library(org.Hs.eg.db)
library(GO.db)
row.names(qlf) <- keys(org.Hs.eg.db)[1:nrow(qlf)]
annot <- select(org.Hs.eg.db, row.names(qlf), c("ENTREZID"))
go <- goana(qlf, species="Hs")
kegg <- kegga(qlf, species="Hs")
write.csv(go, file="GO_all.csv")
write.csv(kegg, file="kegg_all.csv")