用limma包进行多组差异表达分析

写在前面：最近在使用limma包进行差异表达分析，参考了网上许多教程都觉得说的云里雾里，很不清楚。经过我自己一段时间非常痛苦的钻研，弄明白了，解决了我的实际需求。于是决定将我的分析经验写下来，分享给需要的人。

首先加载前期预处理好的表达矩阵。（我的原始表达矩阵文件已附在文后，大家有需要可以下载实践）

setwd("E://Rworkplace")
data<-read.table("GSE98793.txt")

表达矩阵格式如下图所示，行名为基因名，列名为样本名。

创建样本分类信息表。

所谓样本分类信息表，实际上就是值你的样本名一一对应的属性是什么（譬如case还是control）。

我自己根据GEO数据库中该芯片的分类信息先自己在excel中整理了一张表。因为这张芯片原作者在上传数据的时候，样本分类数据很乱，实验组和对照组的样本有时候交错分布，譬如ID为GSM2612287的样本是第192个样本，是对照组即健康人群。而ID为GSM2612128的样本是第33个样本，取样自重症抑郁症患者。我们在对样本进行差异表达分析之前，一定要做的工作是弄明白每一个样本的属性，我做这个工作，实际上是对样本的属性进行了手工注释。

然后我将上图中的第二第三列，复制单独保存在一个txt文档中。（即gene_list.txt）

第一列（1,2，3，……）只是为了辅助排序的东西，因为这个1，2,3,就是样本在总的样本中的位置，这个表格就是告诉我们第几个样本，它是对照组，第几个样本，取样自重症抑郁症患者。

1	HC
2	HC
3	HC
4	HC
……	……

加载上述分类表gene_list.txt。并按照数据框的第一列升序排序。

取第二列，这样处理得到的就是每一个位置所对应的样本属性信息。我们最终想要得到的也就是向量group。（实际情况大家可能不会像我这样麻烦，具体问题具体分析，只要得到一个类如group的向量）

group_list<-read.table("gene_list.txt")
new_group<- group_list[order(group_list[,1]),]
group<-new_group[,2]

接着，加载limma包，构造如下矩阵。（这一段可以照抄。根据自己的实际情况替换变量）

suppressMessages(library(limma))
design <- model.matrix(~0+factor(group))
colnames(design)=levels(factor(group))
rownames(design)=colnames(data)

然后，我们规定哪一组数据与哪一组数据比较。

这里也是原来我比较困惑的地方，因为我的数据是有MDD，anxiety两个实验组，HC对照组。而我实际只需要比较MDD与HC的差异表达情况，而一般的教程中并未涉及或者介绍模糊。

下面这段的代码的意思就是比较design矩阵中HC与MDD level。（关于多组分析的其它疑惑，可以参考知乎文章）

contrast.matrix<-makeContrasts("HC-MDD",levels=design)

完成以上两个表格之后，就可以直接跑下面的代码，得到差异表达数据了。（有时间具体分析什么意思）

##step1
fit <- lmFit(data,design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
##step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput)

分析得到的nrDEG如下图所示：

保存分析好的差异表达数据。

write.table(nrDEG,"limma_GSE98793.txt")

全部代码：

setwd("E://Rworkplace")
data<-read.table("GSE98793.txt")
group_list<-read.table("gene_list.txt")
new_group<- group_list[order(group_list[,1]),]
group<-new_group[,2]suppressMessages(library(limma))design <- model.matrix(~0+factor(group))
colnames(design)=levels(factor(group))
rownames(design)=colnames(data)contrast.matrix<-makeContrasts("HC-MDD",levels=design)##step1
fit <- lmFit(data,design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
##step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput)
write.table(nrDEG,"limma_GSE98793.txt")

熬夜写到这儿，一定会有些疏漏和错误，欢迎指正，共同进步！

【数据链接】https://me.csdn.net/download/weixin_40640700

重新修改了上传数据所需的积分，设置为零。可是系统还是自动调为了1。

如果参考数据不能顺利下载的话，可以邮箱2456392738@qq.com,我邮箱发送。

百度网盘链接：

链接：https://pan.baidu.com/s/1AEqq2t9jo8bitJcQuesI3Q
提取码：nqos

用limma包进行多组差异表达分析相关推荐

python 分析两组数据的差异_R语言limma包差异基因分析（两组或两组以上）
使用limma包进行差异基因分析时,做最多的是两分类的,例如control组和disease组,但也会碰到按照序列进行的分组.这时,如果逐一使用两两比较求差异基因则略显复杂.其实开发limma包的大神 ...
limma包分析差异表达基因
limma利用广义线性模型进行差异表达基因分析,主要用于分析微阵列数据. Data analysis, linear models and differential expression for mi ...
R语言limma包差异表达分析
目录一.数据准备 1.数据加载 2.做分组信息数据 3.表达数据样本ID顺序与样本信息数据匹配二.数据预处理 (1)缺失值处理 (2)离群值处理 (3)数据归一化三.数据探索 (1)查看数据是否 ...
edger多组差异性分析_R语言利用edgeR package进行基因差异表达分析举例
R语言利用edgeR package进行基因差异表达分析举例实验数据: 同一组织,分为两组,control vs treat,每组7例sample.数据第一列为基因名,后14列为对应的count. ...
ballgown包进行基因差异表达分析
ballgown包可以读入Stringtie 的转录组组装及定量数据,进行基因差异表达分析. 1. 数据读入 # if (!requireNamespace("BiocManager&quo ...
edgeR：一个数字基因表达数据差异表达分析Bioconductor程序包
edgeR:一个数字基因表达数据差异表达分析Bioconductor程序包人们希望在不久的将来,对于许多功能基因组学应用,新兴的数字基因表达(digital gene expression,DGE) ...
edgeR、limma、DESeq2三种差异表达包比较（RNA-seq数据）
文章目录 1. 加载R包和输入数据 2. 表达数据整理 3.edgeR包做差异表达 4.limma包做差异表达 5.DESeq2包做差异表达 6.比较三种包差异表达基因筛选结果总结: 1. 加载R包 ...
edger多组差异性分析_使用edgeR进行无重复差异表达分析
写这篇文章一部分原因是填2年前的一个坑转录组入门(7):差异表达分析. 另一部分原因是GQ最近又在搞一波无重复的差异表达分析, 所以专门去学了edgeR 我个人是不太推荐没有重复的差异表达分析,毕竟 ...
edger多组差异性分析_【step by step】菜鸟学TCGA（4）-用edgeR做差异表达分析
大家好,工作太忙,太久没有更新了,哎,泪-- 有的同学问我要代码,有的发了,后面的还没有发,一个一个发好累啊,大家有建议吗? 感觉某宝的这个课程也不贵,300多,有经济能力的小伙伴可以自己买,学得快些 ...

用limma包进行多组差异表达分析

用limma包进行多组差异表达分析相关推荐

最新文章

热门文章