r语言degseq2_R语言DESeq 包介绍 -

R语言DESeq包介绍

从RNA序列数据中识别差异性基因表达的一个主要任务是通过DESeq包采用统计方法进行差异表达分析。基于负二项分布和一种收缩机制的应用模型能够有效地估计统计参数。

1. 包的安装

输入以下命令即可调用biocLite(

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

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

除了安装这个数据包外，并且为了在后面介绍包中的方法时使用其名称为pasilla的数据包，请确保已经完成相应的配置。

2. 输入数据和准备

2.1 计数表

在该数据表中第i行第j列的位置代表对应的样本中的特定基因所拥有的reads计数数量。本文所采用的数据集基于此结构进行构建。

属于pasilla数据包，函数system.file用于指示数据文件存储的位置。 > datafile = system.file( > datafile）

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

将pasillaCountTable变量设置为从datafile中导入的带有标题的表格。
获取pasillaCountTable表格的前几行内容。

2.2 元数据

缺乏元数据的数据无法实现有效的价值提取；而元数据通常被划分为三个关键组成部分：样本维度（即行），特征维度（即列）以及实验信息维度。

原文

原句保留

+ row.names = colnames( pasillaCountTable ),

+ condition = c( + \

+ libType = c( + > pasillaDesign

此处采用简洁的R代码设置流程，请注意以下几点：第一，在一般情况下是从扩展表中获取所需数据；第二，在此基础之上首先对配对端样本进行详细分析；第三，请明确以下变量定义：1）pairedSamples变量用于筛选配对端样本；2）countTable变量存储对应计数结果；3）condition变量记录相关条件信息；4）请确保上述操作基于pasillaData框架完成；5）最后请运行以下命令查看前几行数据输入：head(countTable)

condition

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

condition = function( within the vector c( CountDataSet, DESeq包的 central data structure > library() ) )

cds = newCountDataSet( countTable, condition )

2.3 规范化

该函数用于计算cds数据集的规模大小因子 > cds = estimateSizeFactors( cds ) > sizeFactors( cds )

将统计数据中的每一列数值依次除以该列所代表的数据量（即大小因子），从而使得各列的数据规模达到统一。该函数用于执行这一计算任务。

head( counts( cds, normalized=TRUE ) )

3. 方差估计

DESeq推断基于估算典型数据间的变异性与平均值之间的关系，这相当于计算偏差与均值的关系。偏差可被定义为生物变异系数的平方。通过执行以下操作可进行偏差估计：> cds = estimateDispersions( cds )

该函数分为三个步骤：首先计算每条基因的离差估计；然后拟合一条曲线以匹配数据；最后为每个基因分配一个离差值，并从每条基因的估计值与匹配值中选择合适的数值作为结果存储在fitInfo对象中。

函数plotDispEsts能够展示各个基因的估计值与平均正常统计值之间的联系。

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

无论何种情况, 可作为子序列测试所需的数据量的离差值位于cds特征数据集中的相关位置.

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

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

为了分析在untreated和treated条件下是否存在差异表达的情况，我们可以直接应用nbinomTest函数来计算结果。具体操作如下：执行以下操作：res = nbinomTest(cds, head(res))

id baseMean baseMeanA baseMeanB foldChange log2FoldChange pval padj

id baseMean baseMeanA baseMeanB foldChange log2FoldChange pval padj

feature identier

average normalized counts, averaged across both experimental conditions, with results derived specifically for the A and B experimental conditions' average normalized counts, and the fold difference between these two conditions

the logarithm (to basis 2) of the fold change

pvalue for the statistical signicance of this change

The p-value was corrected for multiple testing using the Benjamini-Hochberg procedure (as implemented in the R function p.adjust), which regulates the false discovery rate (FDR).

在分析中, 我们首先绘制了log2处理后的数据与平均正常值的关系, 并用红色标记表示对应的基因. > plotMA(res)

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

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

hist(respval, breaks=100, col="通过FDR筛选出显著的基因."> resSig = res[ respadj < 0.1, ] 我们生成最有效的差异表达基因.

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