TCGA数据库的利用(三)—做差异分析的三种方法

简介: 今天更新TCGA数据库的利用系列第三篇文章,在对TCGA数据进行挖掘时,通常会筛选出来一些表达量显著异常的基因,作为后续研究的对象,这个筛选过程叫做差异分析;本篇文章将分为三大模块对差异分析进行介绍

今天更新TCGA数据库的利用系列第三篇文章,在对TCGA数据进行挖掘时,通常会筛选出来一些表达量显著异常的基因,作为后续研究的对象,这个筛选过程叫做差异分析;本篇文章将分为三大模块对差异分析进行介绍


关于差异分析的官方解释:


差异分析就是将一组资料的总变动量,依可能造成变动的因素分解成不同的部份,并且以假设检定的方法来判断这些因素是否确实能解释资料的变动。


我自己的一点理解:差异分析就是对总体样本数据中非正常数据的筛选。对于转录组数据进行差异分析时,limma 、edgeR、DESeq2这三种程序包都可以(limma相对较受欢迎),大部分科研性文章基本上是用其中一种方法取筛选差异表达基因,但为了使得结果更加准确,在做毕业课题时我把三种方法都做了一遍,把它们结果的交集作为筛选出来的差异表达基因。


不管用那一种方法做差异分析,基本上要做的步骤就是:一,创建表达矩阵;二、创建设计矩阵(分组矩阵);三、得到差异表达分析矩阵。


但不同包基于算法、数据模型不同,所用的函数、筛选标准也大不相同,所以用代码实现时结果有很大的差别。


无论用那种包做差异分析,在做之前必须要保证需要用到的包已经安装成功。在R语言中安装程序包的代码(其中的一种方式)如下:

source('https://bioconductor.org/biocLite.R')#加载的网址都一样
biocLite("edgeR")#把双引号的内容换成你所需要程序包的名称即可

limma做差异分析


传入原始样本基因表达矩阵(表达矩阵格式如下图)



接下来就是对基因表达矩阵进行一些处理,让样本名变成数据矩阵的列名,基因名变成数据矩阵的行名,同时把ensembl_symbl那一列去掉(用express_rec <- express_rec[,(-1)]命令即可),变成如下这个格式:



表达矩阵里面的数据太大,但为了使数据呈现正态分布,需要对数据进行标准化,这里我用的是函数log(express_rec,2),在标准化之前,需要把表达矩阵内为0的数据赋值为1,目的是为了防止取log时,数据变为负无穷大。

(express_rec[express_rec==0<-1)

下面进行分组矩阵的组建,首先提前创建好一个矩阵列表,如下,行名为样本编码,列名为样本类型,如下面这种格式:

而limma包用到的设计矩阵是下面这种格式:

实现代码如下:

rownames(group_text) <- group_text[,1]
group_text <- group_text[c(-1)]
Group <- factor(group_text$group,levels = c('Tumor','Normal'))
design <- model.matrix(~0+Group)
colnames(design) <- c('Tumor','Normal')
rownames(design) <- rownames(group_text)

接下来的步骤依次进行数据拟合、经验贝叶斯检验、筛选差异表达基因

fit <- lmFit(express_rec,design)
#制作比对标准;
contrast.matrix <- makeContrasts(Tumor - Normal,levels=design)
fit2 <- contrasts.fit(fit,contrast.matrix)
#进行经验贝叶斯检验;
fit2 <- eBayes(fit2)
#基于logFC为标准,设置数量上限为30000,调整方法为fdr;
all_diff <- topTable(fit2, adjust.method = 'fdr',coef=1,p.value = 1,lfc <- log(1,2),number = 30000,sort.by = 'logFC')#从高到低排名;


limma包的另一种方法,精确权重法(voom),然后把筛选得到的差异表达基因写入csv文件中。

#limma的另一种方法;
dge <- DGEList(counts = express_rec)
dge <- calcNormFactors(dge)#表达矩阵进行标准化;
v <- voom(dge, design,plot=TRUE)#利用limma_voom方法进行差异分析;
fit <- lmFit(v, design)#对数据进行线性拟合;
fit <- eBayes(fit)#贝叶斯算法组建
all <- topTable(fit, coef=ncol(design),n=Inf)#从高到低排名;
sig.limma <- subset(all_diff,abs(all$logFC)>1.5&all$P.Value<0.05)#进行差异基因筛选;
write.csv(sig.limma,'C:/Users/FREEDOM/Desktop/TCGA_data/limm_diff.csv')#写入csv文件中;

DESeq2做差异分析

第一步跟limma程序包一样,读入表达矩阵,对表达矩阵进行数据处理

 
         
express_rec<- read.csv('C:/Users/FREEDOM/Desktop/TCGA_data/after_note2.csv')#读取数据
group_text <- read.csv('C:/Users/FREEDOM/Desktop/TCGA_data/group_text.csv')
library('DESeq2')#加载包;
install.packages('rpart')#含有这个包可忽略,没有的时候才安装;
express_rec <- express_rec[,-1]
rownames(express_rec) <-express_rec[,1]
express_rec <- express_rec[(-1)]#表达矩阵的处理;
rownames(group_text) <- group_text[,1]

创建分组矩阵

rownames(group_text) <- group_text[,1]
group_text <- group_text[c(-1)]#分组矩阵的数据处理;
all(rownames(group_text)==colnames(express_rec))#确保表达矩阵的列名与分组矩阵行名相一致;

构建 DESeq2 所需的 DESeqDataSet 对象

dds <- DESeqDataSetFromMatrix(countData=express_rec, colData=group_text, design<- ~ group)  #DESeq2的加载
head(dds)
dds <- dds[rowSums(counts(dds)) > 1, ] #过滤一些low count的数据;

使用DESeq进行差异表达分析,返回 results可用的DESeqDataSet对象

> dds <- DESeq(dds)#DESeq进行标准化;
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 2819 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing

可以用summary函数查看表达基因上下调分布基本情况

> summary(res)#查看经过标准化矩阵的基本情况;
out of 24823 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 7590, 31%
LFC < 0 (down)     : 5120, 21%
outliers [1]       : 0, 0%
low counts [2]     : 0, 0%
(mean count < 0)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

最后提取差异表达基因,存入csv文件中。

edgeR进行差异分析

这个包的操作步骤与limma包类似,首先就是读入数据,创建表达、分组矩阵。

path = 'C:/Users/FREEDOM/Desktop/TCGA_data/after_note2.csv'
path1 ='C:/Users/FREEDOM/Desktop/TCGA_data/group_text.csv'
express_rec <- read.csv(path,headers <- T)#读取表达矩阵;
group_text <- read.csv(path1,headers <- T)#读取分组矩阵;
library(edgeR)#加载edgeR包
express_rec <- express_rec[,-1]
rownames(express_rec) <- express_rec[,1]
express_rec <- express_rec[(-1)]#创建表达矩阵;
rownames(group_text) <- group_text[,1]
group_text <- group_text[c(-1)]#加载分组矩阵;
group <-factor(group_text$group)
dge <- DGEList(counts = express_rec,group = group)#构建DEList对象;
y <- calcNormFactors(dge)#利用calcNormFactor函数对DEList对象进行标准化(TMM算法) 
#创建设计矩阵,跟Limma包相似;
Group <- factor(group_text$group,levels = c('Tumor','Normal'))
design <- model.matrix(~0+Group)
colnames(design) <- c('Tumor','Normal')
rownames(design) <- rownames(group_text)#创建分组矩阵;

edgeR包创建的分组矩阵与limma一样,是以factor因子格式展现出来

接下来依次进行构建DGEList对象、利用TMM算法对数据进行标准化、估计离散值、数据拟合、创建对比矩阵、对数据做QL-text检验、差异表达基因写入csv文件中

dge <- DGEList(counts = express_rec,group = group)#构建DEList对象;
y <- calcNormFactors(dge)#利用calcNormFactor函数对DEList对象进行标准化(TMM算法) 
y <- estimateDisp(y,design)#估计离散值(Dispersion)
fit <- glmQLFit(y, design, robust=TRUE)#进一步通过quasi-likelihood (QL)拟合NB模型
TU.vs.NO <- makeContrasts(Tumor-Normal, levels=design)#这一步主要构建比较矩阵;
res <- glmQLFTest(fit, contrast=TU.vs.NO)#用QL F-test进行检验
# ig.edger <- res$table[p.adjust(res$table$PValue, method = "BH") < 0.01, ]#利用‘BH’方法;
result_diff <- res$table#取出最终的差异基因;
write.csv(edge_diff,'C:/Users/FREEDOM/Desktop/TCGA_data/edgeR_diff2.csv')
相关文章
|
1月前
|
SQL Linux 数据库
|
2天前
|
存储 关系型数据库 MySQL
mysql数据库查询时用到的分页方法有哪些
【8月更文挑战第16天】在MySQL中,实现分页的主要方法包括:1)使用`LIMIT`子句,简单直接但随页数增加性能下降;2)通过子查询优化`LIMIT`分页,提高大页码时的查询效率;3)利用存储过程封装分页逻辑,便于复用但需额外维护;4)借助MySQL变量实现,可能提供更好的性能但实现较复杂。这些方法各有优缺点,可根据实际需求选择适用方案。
|
5天前
|
存储 SQL 数据库
触发器的设计、掌握存储过程的基本概念和创建、执行、删除方法。掌握数据库备份的方法和数据库恢复的方法。
这篇文章介绍了数据库中触发器的设计概念,包括创建、修改、删除触发器的方法,并通过实验内容教授如何使用SQL命令创建DML触发器以及如何利用触发器实现数据的完整性和自动化处理。
触发器的设计、掌握存储过程的基本概念和创建、执行、删除方法。掌握数据库备份的方法和数据库恢复的方法。
|
1天前
|
存储 SQL 关系型数据库
关系型数据库和非关系型数据库的区别和选择方法?
【8月更文挑战第17天】关系型数据库和非关系型数据库的区别和选择方法?
4 0
|
1月前
|
Java 数据库连接 应用服务中间件
数据库连接管理与性能优化方法
数据库连接管理与性能优化方法
|
1月前
|
存储 人工智能 分布式数据库
现代数据库技术的发展与应用前景分析
随着信息时代的发展,数据库技术在各行各业中扮演着至关重要的角色。本文探讨了现代数据库技术的最新发展趋势,以及其在未来的应用前景,涵盖了分布式数据库、区块链技术与数据库融合、人工智能驱动的数据管理等领域。
|
1月前
|
Oracle 关系型数据库 数据库
Oracle数据恢复—Oracle数据库误删除的数据恢复方法
相信有很多oracle数据库用户都遇到过在操作Oracle数据库时误删除某些重要数据的情况,这个时候如果数据库没有备份且数据十分重要的,怎么才能恢复误删除的数据呢?北亚企安数据恢复工程师下面简单介绍几个误删除Oracle数据库数据的恢复方法。
|
1月前
|
小程序 JavaScript 安全
【微信小程序-原生开发】转发给好友/群,分享到朋友圈(含单页模式访问云开发数据库的方法)
【微信小程序-原生开发】转发给好友/群,分享到朋友圈(含单页模式访问云开发数据库的方法)
40 0
|
1月前
|
存储 JavaScript 前端开发
文本,三步走构思,富文本点击提交能够存储到数据库当中(下),最快的方法,还是会看资料,因此会整合资料最好,直接看资料最快,因为是JWT的资料,我们要设置好登录的内容,看登录的地方怎样写的
文本,三步走构思,富文本点击提交能够存储到数据库当中(下),最快的方法,还是会看资料,因此会整合资料最好,直接看资料最快,因为是JWT的资料,我们要设置好登录的内容,看登录的地方怎样写的
|
1月前
|
Java 数据库连接 应用服务中间件
数据库连接管理与性能优化方法
数据库连接管理与性能优化方法