搜索
第22章数据分析2020/5/20在微阵列实验中，基因表达数据的后续分析显得尤为重要，Bioconductor数据分析类的14个包提供了对基因表达数据的后续处理方法。•daMA包提供了双色因子微阵列实验的设计方法以及对实验结果的分析方法；•edd包和factDesign包应用于微阵列表达数据分析和数据挖掘；•genefilter包提供了一系列根据芯片实验数据筛选基因的工具；•globaltest包检测一组基因是否和感兴趣的临床诊断结果存在显著性；•gpls包和pamr包主要对基因表达数据进行分类；•limma包提供一系列分析基因表达芯片数据工具的库，特别是使用线性模型来分析实验以及评估差异表达，提供了对多个RNA目标基因进行分析比较的工具；•siggenes利用芯片的相关性分析以及对芯片的经典贝叶斯分析确定不同的表达基因并估计检出率；•ROC包是和受试者工作特征曲线(ROC)相关的R语言和函数的集合，这些函数对DNA芯片实验进行ROC分析；•splicegear包提供了对可变剪接起作用的一套工具；MeasurementError.cor包是相关系数的尺度误差模型估计。22.1daMA包daMA包主要应用于设计双色因子微阵列实验，并对相关的实验结果数据进行统计分析。双色微阵列实验应用于检测基因在实验样本与对照样本中的差异表达。一个因子的情况比较简单，只要用Cy3和Cy5两种染料分别标记实验样本和对照样本，然后与同一块芯片杂交，检测信号，分析数据就能实现目的。而对于某些需要考察2个或2个以上的因子的实验，每个因子有2个或2个以上的水平，这时实验设计显得尤为重要，因为好的实验设计对于获取可靠的数据是至关重要的，对这些数据的后续分析和生物学解释也更有意义。daMA包就是用于2因子多水平这类微阵列实验的设计和数据分析的。22.1.1简介22.1.2基本函数基本函数包括designMA和analyseMA，这些函数都涉及一些共同的向量和矩阵，它们也是与微阵列设计相关的一些数据结构，所以在介绍函数之前先给出描述。daMA包中的基本数据结构见表22.1（p352），在描述这些数据结构之前，先介绍设计矩阵X，这是一个N×(K+2)的矩阵，行表示微阵列实验，第1、2列表示样本的标记情况，一般有两种不同的染色（绿和红），后面的K列表示对样本的K种处理（如果某一列的元素值等于1则表示样本经过了对应的处理，若等于0表示样本没有经过对应的处理）。设计矩阵确定了每个实验问题对应的染色和样本的处理情况。表22.1daMA包中的基本数据结构┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┃┃描述实验问题的对照矩阵。这个数值矩阵描述了实验问题，每个实验问题用一个对照向量┃┃（cmat中的一行）或一个对照矩阵（cmat中的某些行）来描述，列的次序依赖于设计矩阵┃Cmat┃X。因此，每行的头两个元素都表示两种染色(Cy3，Cy5)，如对照实验仅使用两种染色，┃┃则我们可以设定对应行向量为（1，-1，0．…，0），-1表示红色．1表示绿色，0表示不标记┣━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┃┃描述对照矩阵cmat中行结构的向量。cinfo中的元素个数对应于实验问题的数目。1表示┃Cinfo┃向量形式的对照，整数n1，表示n次对照，以矩阵形式出现。通常，对于每一个对照矩┃┃阵cmat都必须创建这样一个向量。同时，cinfo各元素的和应当等于cmat矩阵的行数┣━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┃┃描述实验问题的常数对照矩阵。这个数值矩阵描述了实验问题，每一个实验问题用一个对┃┃照向量（cmat中单行）或者一个对照矩阵（cmat中的数行）描述。列的顺序依赖于对应的┃cmatB.AB┃设计矩阵X，因此，每行的头两个元素都保存为两种染色Cy3和Cy5．矩阵cmatB.AB的第一行┃┃描述了主效应B┣━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┃┃描述对照矩阵cmatB.AB中行结构的向量。它的第一个元素B指出第一个实验问题(主效┃cinfoB.AB┃应B)，是用向量形式的单个对照来描述的，第二个元素AB指出第二个实验问题（A和B┃┃间的交互作用），是以对照矩阵的形式给出的。每一个对照矩阵都必须创建对应的向量┣━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┃┃300012X18的微阵列实验数据矩阵。这些数据来自微阵列实验，包括使用药物或不使用药┃┃物处理的三个细胞系的表达谱数据。cDNA微阵列是由300012个人类cDNA探针点样制成┃Data.3×2┃的，数据矩阵由300012行以及18列组成，每一行表示一个探针点，每一列表示一次微阵列实验┣━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━1．微阵列实验设计函数designMA用于因子微阵列实验的设计，它根据用户定义的设计矩阵列表，设计有效的双色因子微阵列实验。每一个设计矩阵描述一个实验问题（对照）。使用形式为designMA(design.list,cmat,cinfo,type=c(d,e,t),tol=1e-06)其中，design.list为设计矩阵的名字列表，每个设计矩阵必须满足的条件是矩阵总行数等于微阵列芯片数，列数等于实验条件数，头两列描述染色标签（绿色或红色），余下的列描述实验条件；type用于表明使用哪种最优性判据；“d”为决定性的；“e”为特征值；“t”为回溯；tol用于指示估计对比检查误差的值。2．实验结果分析函数analyseMA基于实验设计对双色因子微阵列实验结果进行分析。实验设计包括，以对照形式包含实验问题的用户自定义矩阵(cmat)以及区分向量和矩阵对照的向量(cinfo)。对每个探针点单独执行该分析。首先去除那些数据中含有NA值的探针点，然后针对每次的实验问题（对照向量或者对照矩阵）对结果进行线性方程的估计。如果感兴趣的线性方程是可估计的，则计算t检验或F检验。使用形式为analyseMA(data,design,id,cmat,cinfo,padj=c(“none,”bonferroni,fdr),tol=1e-06)其中，data为维数为G×N的矩阵，包含用来分析的归一化或者标准化的数据，G为探针数目，N为微阵列实验数目，该矩阵的每一行数据对应一个基因，矩阵的列数应对应于微阵列实验，这样，每一列包含的数据对应于单个微阵列，矩阵中不能含有任何ID变量，这需要单独输入，缺失的数据被当作NA值来处理；design为大小为N×(K+2)的设计矩阵，K为实验条件的数目，这是源于线性模型理论的设计矩阵X，其元素取值为0，1和-1，0表示相关的联合参数并没有应用到相应的观察中，前两列被保留为两种染色Cy3和Cy5，并且常常被分别填充为1和-1；id为ID向量，长度为G，用于标记每个基因；padj为一个被引用的字符串，用于指示可能用到的多样性调节；“none”为没有多样性调节；“bonferroni为Bonferroni单步调节；“fdr”为Benjamini和Hochberg线性step-up程序；tol为用于指示估计对比检查误差的值。3．其余数据designs.basic是描述双色因子3×2微阵列实验的基本设计矩阵列表。矩阵的行表示微阵列，矩阵的列表示参数。Designs.composite是描述使用18个微阵列进行3×2双色因子微阵列实验的合成设计的矩阵列表，包含10个18X9的矩阵列表。22.2edd包edd是表达密度诊断学(expressiondensitydiagnostics)的缩写，该包提供封装的参考分布函数，计算每个基因表达谱的分布函数，并根据分布函数对基因进行分类。22.2.2基本用法edd包用于微阵列实验表达数据的探索性数据分析。基因表达谱的分布可以是常见的分布或分布的组合，edd把这些常见的分布作为参考分布，对于新的基因表达谱，通过采样、分类等迭代过程确定其分布，并根据分布对基因表达谱进行分类。输入的数据是exprSet类的基因表达谱数据，输出的是每个基因表达谱的分布。22.2.1简介分类步骤如下。首先，对基因特异性表达数据向量进行零均值以及单位化平均绝对偏差(MAD:MeanAbsoluteDeviation)的变换。然后，建立变换分布的参考目录，定义一系列感兴趣的谱分布。每一个变换的表达向量都关联到参考目录中的一个元素（'doubt'或奇异点）。谱分布的估计在数据挖掘中是基本的行为，常常使用盒图、柱状图、密度估计以及其他工具达到。对于微阵列数据的分析，主要问题是如何有效地处理成千上万个基因表达谱的分布。在考虑分布分类时，我们通常使用的方法是从高斯分布或者更普遍的描述分布的特征出发，用简单的图形得到很好的效果。下面通过示例来说明edd包的用法。1．加载edd包和基因表达谱数据基因表达谱数据保存在eset中，这是exprSet类的对象，包含500个基因26个样本的数据。library(edd);data(eset);print(getSlots(eddDist));stubparmsmedianmadtagplotlimlatexTagcharacternumericnumericnumericcharacternumericcharacterstub是一个被预先设定为“p”“d”“r”“q”的字符串，用来得到从分布中计算累计分布函数、密度、采样或者分位数的R函数的名称；向量parms指出分布的参数，一般被指定为数字型；向量的均值和标准方差有时无法通过直接计算得到，因此在模拟仿真时预先将之存储在变量median和mad中以供参考；tag保存一个字符串，用来清楚地指出邻近的分布；在必须使用latexrendering时，latexTag使用数学符号；plotlim是一个数字型二维向量，指定一个分布密度绘图的x限制。print(names(eddDistList))；##得到edd中所有分布的名称，结果为：[1]N01T3LN01CS1B82U01B28MIXN1MIXN2print(eddDistList[[1]])；##输出$N01的内容如果要得到某一项的某个slot的值，则可使用类中的函数：CDFname(eddDistList[[1]])结果为：[1]pnorm我们还可以通过相关函数察看每种分布的图形化，如：plotED(eddDistList[[3]])；我们得到如图22.1所示的LN01分布的图形：2．对基因表达谱数据进行分析导入分析需要的包，之所以载入golubEsets包，是因为要用到其中的数据集golubMerge。该数据集包含了47个急性淋巴癌(ALL)病人样本的数据以及25个急性白血病(AML)病人样本的数据，这些样本的数据都是通过Affy-metrixHgu6800芯片分析得到的，包括7129个基因的表达数据。library(golubEsets);data(Golub_Merge);data(golubMerge)；##导入数据包接下来，从golubMerge数据集中筛选出符合我们需要的基因表达水平数据，本例子是从数据集中筛选出最小值大于300并且平均绝对偏差大于中位数的那些基因表达数据。madvec-apply(exprs(Golub_Merge),1,mad);minvec-apply(exprs(Golub_Merge),1,min);keep-(madvecmedian(madvec))&(minvec300);gmfilt-Golub_Merge[keep==TRUE,]；##取出符合条件的表达数据将数据集分成ALL和AML两类样本，从而生成我们所需要的特异性数据集。ALL-gmfilt$ALL.AML==ALL