关于基因差异化的那些事 edger Deseq2和limma的使用及一些总结

在基因差异分析的过程中常见的几大差异化常用的R包主要是edgeR Deseq2 和limma(其中limma 包主要内置于edgeR)
那么在分析的过程中呢，每个包有他们各自的一些数据处理方式。
在接受的数据格式上和样本数目上呢也存在一些差异。主要如下

1.首先关于limma 包进行差异分析
参考的博客主要有
https://cloud.tencent.com/developer/article/1667505(R包limma作差异基因分析）

关于limma 包的一些相关知识
limma 主要用于分析来自芯片的基因表达数据或来自RNA-Seq 的基因表达数据，其中最主要的功能是通过使用线形模型来评估多因子实验的差异表达，limma 主要用于小样本量的数据的差异化分析，其尤其为多实验条件变量和预测因素的复杂实验而设计，这种线性的模型和差异化分析的方式，使得它广泛的应用于芯片，RNA-seq 和定量PCR的数据，limma 可用于单通道和双色芯片通道芯片的差异化分析。
通过将log2-counts per million(Log CPM)，估计均值-方差(Mean-variance)关系以确定在线形建模之前每次观察的权重，之后将数据应用于线性建模，并通过经验贝叶斯统计差异基因。

R包的安装

BiocManager::install('limma')library(limma)

2.limma 差异分析的数据准备

2.1 数据读入和处理

今天我们准备的数据是一个从GEO上下载下来的FPKM 数据
那

首先通过read.table ()函数，将下载的压缩包的文件进行读取
exprset <- read.table(file = 'GSE130437_genes.fpkm_table.txt.gz',sep = '\t',header = T,quote = '',fill = T, comment.char = "!",check.names = T) #读取表达数据rownames(exprset) = exprset[,1] #将第一列作为行名
exprset <- exprset[,-1] #去掉第一列
names(exprset) <- names(exprset) %>% substr(3,nchar(names(exprset))-1) # 更改列名
View(exprset)

得到的数据如下

2.2 数据的FPKM 向TPM 的转换

那么接下来，我们需要将FPKM 转化为TPM，主要参考的是生信技能树相关的博客帖子(https://cloud.tencent.com/developer/article/1558031
RNAseq数据，下载GEO中的FPKM文件后该怎么下游分析)

这个地方一个很关键的分析主要是对FPKM 到TPM 的转化，因此下面的那个函数也就是根据此得来的(关于FPKM 和TPM 的一个转换，

主要参考博文
(https://cloud.tencent.com/developer/article/1669450
RNA-Seq的Counts和FPKM数据如何转换成TPM？

https://www.jianshu.com/p/9dfb65e405e8
简单小需求：如何将FPKM转换成TPM？)

expMatrix <- exprset##将上述的FPKM 数据赋给新的expMatrix
fpkmToTpm <- function(fpkm)
{exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}###构建一个将FPKM 转化为TPM 的函数
tpms <- apply(expMatrix,2,fpkmToTpm)###利用apply()函数对expMatrix的每一列进行FPKM 到TPM 的转化。
tpms[1:3,]
colSums(tpms)
#输出结果：
tpms[1:3,]MCF7.Palbo.Replicate.1 MCF7.Palbo.Replicate.2 MCF7.Palbo.Replicate.3 MCF7.Parental.Replicate.1 MCF7.Parental.Replicate.2
ENSG00000000003               24.89911                  23.23                  20.71                     10.16                     11.11
ENSG00000000005                0.02938                   0.00                   0.00                      0.00                      0.00
ENSG00000000419              161.86588                 177.54                 146.70                    194.89                    184.97MCF7.Parental.Replicate.3 MDAMB231.Palbo.Replicate.1 MDAMB231.Palbo.Replicate.2 MDAMB231.Palbo.Replicate.3
ENSG00000000003                     11.35                    29.2956                   26.65212                    27.9496
ENSG00000000005                      0.00                     0.1845                    0.05373                     0.1782
ENSG00000000419                    168.17                   147.5768                  121.83891                   123.6252MDAMB231.Parental.Replicate.1 MDAMB231.Parental.Replicate.2 MDAMB231.Parental.Replicate.3
ENSG00000000003                         10.36                         9.635                         10.25
ENSG00000000005                          0.00                         0.000                          0.00
ENSG00000000419                        113.82                       118.264                        114.09colSums(tpms)MCF7.Palbo.Replicate.1        MCF7.Palbo.Replicate.2        MCF7.Palbo.Replicate.3     MCF7.Parental.Replicate.1 1e+06                         1e+06                         1e+06                         1e+06 MCF7.Parental.Replicate.2     MCF7.Parental.Replicate.3    MDAMB231.Palbo.Replicate.1    MDAMB231.Palbo.Replicate.2 1e+06                         1e+06                         1e+06                         1e+06 MDAMB231.Palbo.Replicate.3 MDAMB231.Parental.Replicate.1 MDAMB231.Parental.Replicate.2 MDAMB231.Parental.Replicate.3 1e+06

2.3limma 的差异分析的过程

构建group_list
group_list <- factor(rep(c('control', 'treat'), each = 3), levels = c('control', 'treat'))#表达矩阵数据校正
exprSet <- tpms
boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2) ###通过上述的FPKM 到TPM 的转换过程中，我们要看看数据是否在每个样本里进行了归一化处理，见图1，那我们可以看到 经过转换后的数据在每个样本里并未进行归一化处理，因此需要进行归一化
library(limma)
exprSet=normalizeBetweenArrays(exprSet)  ###对数据进行归一化处理
boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)  ##进行归一化的数据的可视化，见图2通常情况下，我们为了数据能够正常的进行后续的分析 会对数据进行Log 转换并且对原始数据+1
#判断数据是否需要转换
exprSet <- log2(exprSet+1)#差异分析：
dat <- exprSet
design=model.matrix(~factor( group_list ))
fit=lmFit(dat,design)
fit=eBayes(fit)
options(digits = 4)
topTable(fit,coef=2,adjust='BH')
bp=function(g){library(ggpubr)df=data.frame(gene=g,stage=group_list)p <- ggboxplot(df, x = "stage", y = "gene",color = "stage", palette = "jco",add = "jitter")#  Add p-valuep + stat_compare_means()
}###其实没太看懂这一步函数的作用
deg=topTable(fit,coef=2,adjust='BH',number = Inf)
head(deg)

图1

图2 经过归一化处理后，我们发现在每个样本中，数据的表达平齐了

关于edgeR 和Deseq2的差异化分析下次遇到的时候，再进行补充。

###ID 转换
library(org.Hs.eg.db)
list=select(org.Hs.eg.db,keys=deg$ENSEMBL,columns = c("ENTREZID","SYMBOL"), keytype="ENSEMBL")
View(list)
MCF_1_0_A_Deg=merge(deg,list,by="ENSEMBL")
View(MCF_1_0_A_Deg)
write.table(MCF_3_0_A_Deg,file="D:/desk_user/Data_of_lab/FX/GSE99116/MCF_1_0_A_Deg.txt",sep="\t")

关于基因差异化的那些事 edger Deseq2和limma的使用及一些总结相关推荐

公主恋人汉化组的事我出来说明下
这样的: 有这么件事就是汉化组 [宣传]给我看了某帖子反正就是与我有关吧,根本来说在喷公主恋人汉化组 [宣传]就问我该怎么办 (三大妈某娃娃: 怕什么,换马甲上去喷他不就得了) 我这里就不说 ...
【生活】年化收益率、七日年化收益率这些事
缘由: 互联网人使用互联网的工具越来越多,支付宝.微信已是必备的手机App,相信我们猿猿们使用最多的理财就是"财付通"和"理财通"了,然而对其中的一些名词应该不 ...
DESeq2，EdgeR，limma对比
参考: [陈巍学基因]RNA-seq - 知乎关于基因差异化的那些事 edger Deseq2和limma的使用及一些总结_forever luckness 的博客-CSDN博客_des ...
生物学重复吗？还有技术重复？
转录组测序多少生物重复合适?2个?3个?48个? 生信宝典 2020-03-30 21:10:44 664 收藏 1 分类专栏: 生物信息版权 2016年英国邓迪大学的Geoffrey J Ba ...
如果不是没有钱，谁想测3个重复？
这篇文章上次发出后,有朋友留言说到底要测几个重复?其实也没有定论,有钱多多益善.只是需要知道重复少时,发现的差异基因会有不少假阴性,获得不了结果时,可尝试加测一些,可获得更稳定的结果. 重复少时,抽 ...
转录组测序多少生物重复合适？2个？3个？48个？
2016年英国邓迪大学的Geoffrey J Barton教授在RNA发表一篇文章专门评估这一问题.作者对野生型和snf2突变型酵母样品分别测序了48个生物学重复:质控后,野生型样品保留42个生物学重 ...
MOVICS系列教程(三) RUN Module
前言通过学习前面几个模块,我们已经发现了基于多组学数据找出的乳腺癌各亚型间具有非常显著的分子差异,而我们如果想深入挖掘其背后机制,就需要找出差异表达的基因是哪些,以及这些基因具有什么样的功能.而MO ...
edgeR、limma、DESeq2三种差异表达包比较（RNA-seq数据）
文章目录 1. 加载R包和输入数据 2. 表达数据整理 3.edgeR包做差异表达 4.limma包做差异表达 5.DESeq2包做差异表达 6.比较三种包差异表达基因筛选结果总结: 1. 加载R包 ...
关于单细胞批次矫正那些事(二) KBET 用于单细胞批次矫正结果的评估
之前我们根据GB 那篇文章了解到了对于常见的用于单细胞批次矫正分析的常见的方法Harmony LIGER以及Seurat3 等,那么文章也推荐如果是批次矫正后用于下游的分析的话推荐ComBat,MNN ...

关于基因差异化的那些事 edger Deseq2和limma的使用及一些总结

关于基因差异化的那些事 edger Deseq2和limma的使用及一些总结相关推荐

最新文章

热门文章