limma做RNAseq差异分析

2020-04-01 16:33:31 浏览数 (1)

limma是一个很强大的用于分析芯片的R包,也可以用于RNA-Seq的差异分析 以两个组比较为例:首先输入count表达矩阵,这里也跟其他差异分析R包一样,不要输入已经标准化的数据。 本文主要参考:https://www.bioinfo-scrounger.com/archives/115/

代码语言:javascript复制
library(limma)
library(edge)

counts <- read.csv("raw_counts.csv",row.names = 1)
#Creates a DGEList
dge_counts <- DGEList(counts = counts,remove.zeros = T)
#Calculate normalization factors to scale the raw library sizes.
#用TMM进行标准化
tmm_counts <- calcNormFactors(dge_counts)
#count进行标准化以及转化为log2的值
logCPM_counts <- cpm(tmm_counts , log=TRUE)

制作分组矩阵

代码语言:javascript复制
#设置分组信息
group_list <- factor(c(rep("control",2), rep("treat",2)))
design <- model.matrix(~group_list)
colnames(design) <- levels(group_list)
rownames(design) <- colnames(counts)

为了避免文库大小在样本间变化的影响,可以使用limma的voom方法进行处理

代码语言:javascript复制
y = voom(logCPM_counts, design, plot = T)

对voom的描述

Transform count data to log2-counts per million (logCPM), estimate the mean-variance relationship and use this to compute appropriate observation-level weights. The data are then ready for linear modelling.

voom()作用是原始counts转换为logCPM值,将所有计数加0.5,以避免取对数零。然后,将logCPM值矩阵进行标准化。在运行voom()之前,应对counts矩阵进行过滤以除去counts非常低的基因。为此,可以使用edgeR包中的filterByExpr函数。

image.png

因为我已经通过edgeR的TMM标准化,所以效果还行,不需要处理,如果效果不好可能由于数据没有过滤好,如下图:

image.png

对voom图的解释可以参考这里: https://stats.stackexchange.com/questions/160255/voom-mean-variance-trend-plot-how-to-interpret-the-plot

差异分析:

代码语言:javascript复制
#不需要voom转化时:
fit <- lmFit(logCPM_counts, design)
fit <- eBayes(fit)
DE_genes <- topTable(fit, coef = 2,p.value = 0.5, lfc = log2(1.5), number = Inf,sort.by="logFC")
代码语言:javascript复制
#不进行TMM转化,即不运行calcNormFactors(),直接进行voom转化
y = voom(counts, design, plot = T)
fit <- lmFit(y, design)
fit <- eBayes(fit)
DE_genes <- topTable(fit, coef = 2,p.value = 0.05, lfc = log2(1.5), number = Inf,sort.by="logFC")

欢迎关注~

参考: https://www.bioinfo-scrounger.com/archives/115/ https://www.jianshu.com/p/616de0ee881a https://mp.weixin.qq.com/s?__biz=MzI4NjMxOTA3OA==&mid=2247483987&idx=1&sn=aa2ca81e7fe128edaaedc47479c517c9&chksm=ebdf8adadca803cc31261a1ccabf8a6bdb835ce12b670cc01969fcfde51cfdb991d0619e7695&scene=21#wechat_redirect

0 人点赞