目前主流的差异分析流程，用的基本都是DESeq2或者edgeR。DESeq2之前说过了。这次来写edgeR。
首先在R中安装
source("https://bioconductor.org/biocLite.R")
biocLite("edgeR")
# 跑的时候提醒我少了个statmod，一起装上
install.packages("statmod")
然后，这次依然是用featureCounts的结果来跑。 可以查看上次的文章。
拿到了all_feature.txt这个结果之后，就可以用在R里面跑了。
library(edgeR)
library(statmod)
# 设置样本的名字
sampleNames <- c("siCtrl_1", "siCtrl_2", "siSUZ12_1", "siSUZ12_2")
# 导入counts结果
data <- read.table("all_feature.txt", header=TRUE, quote="\t", skip=1)
# 把counts数提取处理并且命名
names(data)[7:10] <- sampleNames
# 转换成matrix
countData <- as.matrix(data[7:10])
# 把行名改成基因名
rownames(countData) <- data$Geneid
# 设置condition，也就是对照的内容
condition <- factor(c("siCtrl", "siCtrl", "siSUZ12", "siSUZ12"))
接下来是创建edgeR的文件并且过滤
genelist <- DGEList(counts=countData, group=condition)
# 这里的0.4过滤阈值需要自行调整，一般是用最小的read count总数除以1000000
keep <- rowSums(cpm(genelist) > 0.4) >= 2
table(keep)
genelist.filted <- genelist[keep, , keep.lib.sizes=FALSE]
genelist.norm <- calcNormFactors(genelist.filted)
然后是弄一个condition表格，和DESeq2的condition类似
design <- model.matrix(~0+condition)
colnames(design) <- levels(condition)
design
估计离散值顺便画个图
genelist.disp <- estimateDisp(genelist.norm, design, robust=TRUE)
plotBCV(genelist.disp)
fit <- glmQLFit(genelist.disp, design, robust=TRUE)
计算各种值，用的是和DESeq2一样的BH算法
edg <- makeContrasts(siCtrl-siSUZ12, levels=design)
res <- glmQLFTest(fit, contrast=edg)
# 我在这里我padj值也放进去了
res$table$padj <- p.adjust(res$table$PValue, method="BH")
write.csv(res$table, "edg_output.csv")
write.csv(merge(res$table, countData, by="row.names", sort=FALSE), "all_edg_output.csv", row.names=FALSE)
结果这样就是出来了，然后可以根据结果画各种图或者做后续的通路分析了。