Q：該如何選擇limma, DESeq2, edgeR
A：各有各自應(yīng)用的場景
- 如果是芯片數(shù)據(jù)第晰，一般選擇limma颓屑，畢竟它是處理芯片數(shù)據(jù)之王棕叫。 不過edgeR也可以仗岸。芯片數(shù)據(jù)默認(rèn)符合正態(tài)分布允耿，而limma正是基于正態(tài)分布。
- 如果是二代測序的原始count值扒怖，一般選擇DESeq或edgeR较锡。注意這兩者只能處理count，不能處理FPKM等矯正后的數(shù)據(jù)盗痒。二代測序數(shù)據(jù)符合柏松分布蚂蕴，理論上不能用T檢驗(yàn)低散，只能用非參數(shù)檢驗(yàn)（秩和），但是統(tǒng)計(jì)力度不夠骡楼，所以還是得用經(jīng)過矯正后的參數(shù)檢驗(yàn)熔号。
- 如果是FPKM等矯正后的表達(dá)量，可以用cuffdiff
Sum: 基于以上鸟整，對于二代測序數(shù)據(jù)引镊，先拿到原始count值進(jìn)行DESeq2差異分析，再轉(zhuǎn)換成TPM進(jìn)行下游分析篮条。不建議用edgeR和cuffdiff弟头。
DESeq找差異基因+火山圖 傳送門：http://www.reibang.com/p/7c6a15eddfb6
今天寫的是limma進(jìn)行差異分析兑燥。PS：limma的說明書實(shí)在太太太長了亮瓷，不過卻是入門生信必讀的。
Step1: 準(zhǔn)備數(shù)據(jù)
-
group數(shù)據(jù)
group分組信息 -
表達(dá)譜數(shù)據(jù)
Step2: 構(gòu)建design 樣本分組矩陣
library(limma)
design <- model.matrix(~0+group)
rownames(design) = colnames(expr1)
colnames(design) <- levels(group)
Step3: 構(gòu)建進(jìn)行差異分析的對象
%構(gòu)建DGElist對象
DGElist <- DGEList(counts = expr1, group = group)
% 去除表達(dá)量過低的基因
% Compute counts per million (CPM)
keep_gene <- rowSums(cpm(DGElist) > 1) >= 2
table(keep_gene)
DGElist <- DGElist[keep_gene, ,keep.lib.sizes =FALSE]
Step4：構(gòu)建線性模型
% 計(jì)算列的矯正因子
DGElist <- calcNormFactors( DGElist )
% 將count值轉(zhuǎn)化成log2-counts per million (logCPM)降瞳，準(zhǔn)備進(jìn)行線性回歸
v <- voom(DGElist, design, plot = TRUE, normalize = "quantile")
% 對每一個(gè)基因進(jìn)行線性模型構(gòu)建
fit <- lmFit(v, design)
% 構(gòu)建比較矩陣
cont.matrix <- makeContrasts(contrasts = c('cancer-normal'), levels = design)
% 構(gòu)建芯片數(shù)據(jù)的線性模型嘱支，計(jì)算估計(jì)的相關(guān)系數(shù)和標(biāo)準(zhǔn)差
fit2 <- contrasts.fit(fit, cont.matrix)
% 基于貝葉斯計(jì)算T值，F(xiàn)值和log-odds
fit2 <- eBayes(fit2)
Step5：得出結(jié)果
nrDEG_limma_voom = topTable(fit2, coef = 'cancer-normal', n = Inf)
nrDEG_limma_voom = na.omit(nrDEG_limma_voom)
head(nrDEG_limma_voom)
% 篩選出符合要求的差異基因
library(dplyr)
res<-cbind(rownames(nrDEG_limma_voom),nrDEG_limma_voom)
res_1<-res %>% dplyr::filter((logFC>1 | logFC < (-1)) & adj.P.Val < 0.05)
colnames(res_1)[1]<-"Symbol"
