R语言学习转录组转录组流程

limma差异分析

2020-01-07  本文已影响0人  PriscillaBai
Q:该如何选择limma, DESeq2, edgeR
A:各有各自应用的场景

Sum: 基于以上,对于二代测序数据,先拿到原始count值进行DESeq2差异分析,再转换成TPM进行下游分析。不建议用edgeR和cuffdiff。

edgeR的说明书,包括芯片数据

DESeq找差异基因+火山图 传送门:https://www.jianshu.com/p/7c6a15eddfb6
今天写的是limma进行差异分析。PS:limma的说明书实在太太太长了,不过却是入门生信必读的。

Step1: 准备数据
  1. group数据


    group分组信息
  2. 表达谱数据


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"

上一篇 下一篇

猜你喜欢

热点阅读