R语言ESeq 包介绍.docxVIP

R语言DESeq包介绍分析RNA序列数据的一个主要任务是探测基因的差异表达，DESeq包提供了测试差异表达的方法，应用负二项分布和收缩的分布方程估计。包的安装输入如下命令，DESeq和相关的包就可以自动下载和安装。source(/biocLite.R)biocLite(DESeq)相关的包会自动下载安装，安装的包如下：中间会有个选择需要更新相关的包，选择更新全部，更新的包有：另外还要安装一个数据包，供下面介绍包中的方法时使用，包名为pasilla.输入数据和准备2.1 计数表数据表的第i行第j列元素表示第j个样本的第i个基因有多少个reads。本文使用的数据来自于pasilla数据包，函数system.file告诉我们数据文件保存的路径。 datafile = system.file( extdata/pasilla_gene_counts.tsv, package=pasilla ) datafile[1] D:/Program Files/R/R-2.15.3/library/pasilla/extdata/pasilla_gene_counts.tsv在R中读取这个文件，使用read.table函数。 pasillaCountTable = read.table( datafile, header=TRUE, s=1 ) head( pasillaCountTable )2.2 元数据没有元数据的数据是没有用的，元数据可以分为三组，分别是样本(行)，特征(列)和整个实验的信息。首先需要样本的描述信息，data.frame的列表示各种信息，行表示7个样本。 pasillaDesign = data.frame(+ s = colnames( pasillaCountTable ),+ condition = c( untreated, untreated, untreated,+ untreated, treated, treated, treated ),+ libType = c( single-end, single-end, paired-end,+ paired-end, single-end, paired-end, paired-end ) ) pasillaDesign这边简单地使用R代码进行设定，通常情况是从扩展表中读取这些数据。为了分析这些样本，我们需要解释single-end和paired-end的方法，这边首先简单的分析paired-end样本。 pairedSamples = pasillaDesign$libType == paired-end countTable = pasillaCountTable[ , pairedSamples ] condition = pasillaDesign$condition[ pairedSamples ]现在，我们有下面的数据输入 head(countTable) condition对于自己的数据，可以简单的创建因子。 #not run condition = factor( c( untreated, untreated, treated, treated ) )我们现在举例CountDataSet，DESeq包的核心数据结构 library( DESeq ) cds = newCountDataSet( countTable, condition )2.3 规范化函数estimateSizeFactors估计统计数据的大小因子 cds = estimateSizeFactors( cds ) sizeFactors( cds )如果统计数据的每列除以这列的大小因子，这样统计值就变成同一规模，使它们具有可比性。函数counts可以做这个计算。 head( counts( cds, normalized=TRUE ) )方差估计DESeq推断依靠估计典型的数据间的方差和平均关系，或者等效的离差和均值。离差可以理解为生物变异系数的平方。估计离差可以使用以下命令。 cds = estimateDispersions( cds )函数estimateDispersions做了三步，首先估计每条基因的离差，然后，通过估计匹配一条曲线，最后，每个基因分配一个离差，从每条基因估计值和匹配值选一个。为了让用户知道中间过程，fitInfo对象被储存下来。 str( fitInfo(cds) )函数plotDispEsts可以画出每条基因的估计值和平均正常统计值的关系。 plotDispEsts( cds )图1 经验的(黑点)和匹配的(红线)离差值与平均正常统计值的关系图在任何情况下，可以被子序列测试使用的离差值被存储在cds的特征数据集里。 head( fData(c