方法来源:https://www.jianshu.com/p/84bfdf91118c
Q:该如何选择limma, DESeq2, edgeR
A:各有各自应用的场景
如果是芯片数据,一般选择limma,毕竟它是处理芯片数据之王。 不过edgeR也可以。芯片数据默认符合正态分布,而limma正是基于正态分布。
如果是二代测序的原始count值,一般选择DESeq或edgeR。注意这两者只能处理count,不能处理FPKM等矫正后的数据。二代测序数据符合柏松分布,理论上不能用T检验,只能用非参数检验(秩和),但是统计力度不够,所以还是得用经过矫正后的参数检验。
如果是FPKM等矫正后的表达量,可以用cuffdiff
Sum: 基于以上,对于二代测序数据,先拿到原始count值进行DESeq2差异分析,再转换成TPM进行下游分析。不建议用edgeR和cuffdiff。
edgeR的说明书,包括芯片数据
DESeq找差异基因+火山图 传送门:https://www.jianshu.com/p/7c6a15eddfb6
今天写的是limma进行差异分析。PS:limma的说明书实在太太太长了,不过却是入门生信必读的。
Step1: 准备数据
group数据
group分组信息
表达谱数据
Step2: 构建design 样本分组矩阵
library(limma) design <- model.matrix(~0+group) rownames(design) = colnames(expr1) colnames(design) <- levels(group)
Step3: 构建进行差异分析的对象
%构建DGElist对象 DGElist <- DGEList(counts = expr1, group = group) % 去除表达量过低的基因 % Compute counts per million (CPM) keep_gene <- rowSums(cpm(DGElist) > 1) >= 2 table(keep_gene) DGElist <- DGElist[keep_gene, ,keep.lib.sizes =FALSE]
Step4:构建线性模型
% 计算列的矫正因子 DGElist <- calcNormFactors( DGElist ) % 将count值转化成log2-counts per million (logCPM),准备进行线性回归 v <- voom(DGElist, design, plot = TRUE, normalize = "quantile") % 对每一个基因进行线性模型构建 fit <- lmFit(v, design) % 构建比较矩阵 cont.matrix <- makeContrasts(contrasts = c('cancer-normal'), levels = design) % 构建芯片数据的线性模型,计算估计的相关系数和标准差 fit2 <- contrasts.fit(fit, cont.matrix) % 基于贝叶斯计算T值,F值和log-odds fit2 <- eBayes(fit2)
Step5:得出结果
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"
评论区