limma包做差异分析之分组矩阵和比较矩阵
原创
limma包做差异分析之分组矩阵和比较矩阵
原创
小辫子
发布于 2025-06-13 18:03:07
发布于 2025-06-13 18:03:07
title: "limma包做差异分析之分组矩阵和比较矩阵"
output: html_document
editor_options:
chunk_output_type: console
##注：pairinfo为配对信息，解决GEO芯片中配对样本如何做差异分析的问题
##方法一
####分组矩阵model.matrix
#分组矩阵model.matrix
rm(list = ls())
library(limma)
group <- c(rep('before',18),rep('after',18))
group
#> [1] "before" "before" "before" "before" "before" "before"
#> [7] "before" "before" "before" "before" "before" "before"
#> [13] "before" "before" "before" "before" "before" "before"
#> [19] "after" "after" "after" "after" "after" "after"
#> [25] "after" "after" "after" "after" "after" "after"
#> [31] "after" "after" "after" "after" "after" "after"
group <- factor(group)
group
#> [1] before before before before before before before before
#> [9] before before before before before before before before
#> [17] before before after after after after after after
#> [25] after after after after after after after after
#> [33] after after after after
#> Levels: after before
pairinfo = factor(rep(1:18,2));pairinfo
#> [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 1 2
#> [21] 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
#> Levels: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
design <- model.matrix(~0 +group + pairinfo)
design2 <- model.matrix(~0 +group)
design3 <- model.matrix(~0 +pairinfo)"group"的levels为"after"和"befor"，所以design2只有两列；"pairinfo"的levels为c(1:18),所以design3有18列;design同design3，只是前两列的列名是"groupafter"和"groupbefor"
####分组矩阵命名
colnames(design)[1:length(levels(group))] <- levels(group)
colnames(design2) <- levels(group)
colnames(design3) <- levels(pairinfo)
#只是修改了列名而已####比较矩阵
contrast.matrix <- makeContrasts(after - before, levels=design)
contrast.matrix2 <- makeContrasts(after - before, levels=design2)"contrast.matrix"为配对样本的比较矩阵，共18行1列；"contrast.matrix2"为常规两组比较，只有2行1列。这也是后续"coef"取1的原因。
####常规步骤
## 线性拟合
#fit <- lmFit(exprSet, design)
#fit <- contrasts.fit(fit, contrast.matrix)
#fit2 <- eBayes(fit)
## 提取差异结果，注意这里的coef是1
#allDiff=topTable(fit2,adjust='fdr',coef=1,number=Inf) ##方法二
rm(list = ls())
library(limma)
pairinfo = factor(rep(1:18,2))
group <- c(rep('before',18),rep('after',18))
group
#> [1] "before" "before" "before" "before" "before" "before"
#> [7] "before" "before" "before" "before" "before" "before"
#> [13] "before" "before" "before" "before" "before" "before"
#> [19] "after" "after" "after" "after" "after" "after"
#> [25] "after" "after" "after" "after" "after" "after"
#> [31] "after" "after" "after" "after" "after" "after"
group <- factor(group,levels = c("before","after"),ordered = F)
group
#> [1] before before before before before before before before
#> [9] before before before before before before before before
#> [17] before before after after after after after after
#> [25] after after after after after after after after
#> [33] after after after after
#> Levels: before after
design <- model.matrix(~group+pairinfo)
#fit <- lmFit(exprSet,design)
#fit2 <- eBayes(fit)
#最终获取差异基因有两种方式：
#allDiff2=topTable(fit2,adjust='fdr',coef=2,number=Inf ,p.value=0.05)
#allDiff21=topTable(fit2,adjust='fdr',coef="groupafter",number=Inf ,p.value=0.05)
注：这里coef取2或者"groupafter",也就是"design"里的第2列。"makeContrasts"函数的作用只是让组间比较看起来更清晰，不容易出错。两组方法均可，没有优劣之分，自己习惯用哪个就用哪个。
引用自果子学生信
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
评论
作者已关闭评论
推荐阅读
腾讯云开发者
