R语言DESeq包介绍

分析RNA序列数据的一个主要任务是探测基因的差异表达，DESeq包提供了测试差异表达的方法，应用负二项分布和收缩的分布方程估计。

1. 包的安装

输入如下命令，DESeq和相关的包就可以自动下载和安装。 > source(\> biocLite(\

相关的包会自动下载安装，安装的包如下：

中间会有个选择需要更新相关的包，选择更新全部，更新的包有：

另外还要安装一个数据包，供下面介绍包中的方法时使用，包名为pasilla.

2. 输入数据和准备

2.1 计数表

数据表的第i行第j列元素表示第j个样本的第i个基因有多少个reads。本文使用的数据来

自于pasilla数据包，函数system.file告诉我们数据文件保存的路径。 > datafile = system.file( \> datafile

[1] \在R中读取这个文件，使用read.table函数。

> pasillaCountTable = read.table( datafile, header=TRUE, row.names=1 ) > head( pasillaCountTable )

2.2 元数据

没有元数据的数据是没有用的，元数据可以分为三组，分别是样本(行)，特征(列)和整个实验的信息。

首先需要样本的描述信息，data.frame的列表示各种信息，行表示7个样本。 > pasillaDesign = data.frame(

+ row.names = colnames( pasillaCountTable ),

+ condition = c( \+ \

+ libType = c( \+ \> pasillaDesign

这边简单地使用R代码进行设定，通常情况是从扩展表中读取这些数据。为了分析这些样本，我们需要解释single-end和paired-end的方法，这边首先简单的分析paired-end样本。 > pairedSamples = pasillaDesign$libType == \> countTable = pasillaCountTable[ , pairedSamples ] > condition = pasillaDesign$condition[ pairedSamples ] 现在，我们有下面的数据输入 > head(countTable)

> condition

对于自己的数据，可以简单的创建因子。 > #not run

> condition = factor( c( \我们现在举例CountDataSet，DESeq包的核心数据结构 > library( \

> cds = newCountDataSet( countTable, condition )

2.3 规范化

函数estimateSizeFactors估计统计数据的大小因子 > cds = estimateSizeFactors( cds ) > sizeFactors( cds )

如果统计数据的每列除以这列的大小因子，这样统计值就变成同一规模，使它们具有可比性。函数counts可以做这个计算。

> head( counts( cds, normalized=TRUE ) )

3. 方差估计

DESeq推断依靠估计典型的数据间的方差和平均关系，或者等效的离差和均值。离差可以理解为生物变异系数的平方。估计离差可以使用以下命令。 > cds = estimateDispersions( cds )

函数estimateDispersions做了三步，首先估计每条基因的离差，然后，通过估计匹配一条曲线，最后，每个基因分配一个离差，从每条基因估计值和匹配值选一个。为了让用户知道中间过程，fitInfo对象被储存下来。 > str( fitInfo(cds) )

函数 plotDispEsts可以画出每条基因的估计值和平均正常统计值的关系。 > plotDispEsts( cds )

图1 经验的(黑点)和匹配的(红线)离差值与平均正常统计值的关系图

在任何情况下，可以被子序列测试使用的离差值被存储在cds的特征数据集里。 > head( fData(cds) )

4. 推断：称之为差异表达

4.1 两个实验条件下的标准比较

为了看在条件”untreated\和”treated\下是否有差异表达，我们简单的使用nbinomTest函数， > res = nbinomTest( cds, \> head(res)

各列表示的是： id baseMean baseMeanA baseMeanB foldChange log2FoldChange pval padj

feature identier

mean normalised counts, averaged over all samples from both conditions mean normalised counts from condition A mean normalised counts from condition B fold change from condition A to B

the logarithm (to basis 2) of the fold change

pvalue for the statistical signicance of this change

pvalue adjusted for multiple testing with the Benjamini-Hochberg procedure (see the R functionp.adjust), which controls false discovery rate (FDR)

我们首先画log2折叠变换和平均正常统计量的关系，红色的点表示在10yR的基因。 > plotMA(res)

图2 log2折叠变换和平均正常统计量的关系图

图3 从nbinomTest的p值统计直方图

> hist(res$pval, breaks=100, col=\我们可以通过FDR过滤有效的基因。 > resSig = res[ res$padj < 0.1, ] 列举最有效的差异表达基因

> head( resSig[ order(resSig$pval), ] )