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 modelandtesting

-- replacing outliersandrefittingfor2819genes

-- DESeq argument'minReplicatesForReplace'=7

-- original counts are preserved in counts(dds)

estimating dispersions

fitting modelandtesting

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

> summary(res)#查看经过标准化矩阵的基本情况；

outof24823withnonzero totalreadcount

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'argumentof?results

[2] see'independentFiltering'argumentof?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')

以上就是做差异分析三种R包的使用方法，关于本篇文章涉及到的完整源码的获取方式：关注公众号：Z先生点记，后台回复关键词：TCGA差异分析即可。