多性状联合统计方法的关联分析(源码)
许多生物遗传的重要性状（身高、体重等）都是数量性状，要揭示遗传学本质就需要将生物个体的性状与基因型相联系。生物的某个性状有时由多个基因共同决定，不同的基因间有高度的连锁，在影响性状的基因中可能会有基因影响多个性状，因此，性状之间也有着极高的相关性，考虑性状间的相关是很有必要的。关联分析广泛应用于人类遗传学研究中，可将个体的性状与基因型相联系，进行基因定位。本文介绍了统计方法在生物多性状方面的应用，最小二乘法、偏最小二乘法和LASSO是常用的统计方法，可以用来研究性状之间的关联，同时与单性状分析结果进行比较。关键字多性状分析；关联分析；偏最小二乘Association analysis based on multi shape joint statistical methodStudent majoring in Information and Computing Science Feng DaliTutor Zhang JinAbstract: Many important traits (height, weight, etc.) of biological inheritance are quantitative traits. To reveal the nature of genetics, we need to associate the traits of the individual with the genotype. A certain character of a creature is sometimes determined by a number of genes, which have a high degree of linkage between different genes, which may have a genetic influence on multiple traits, therefore, there is a very high correlation between traits, and it is necessary to consider the correlation among traits. Association analysis is widely used in the study of *好棒文|www.hbsrm.com +Q: *351916072*
human genetics. Carry on gene mapping can be linked to the genotype and individual traits. This paper introduces the application of statistical methods in the field of biological multi character. Least squares, partial least squares and LASSO are commonly used statistical methods, can be used to study the association between traits. At the same time, the results were compared with single character analysis.1 绪论1.1研究意义有很多生物遗传方面的重要性状都是数量性状，例如人类疾病的治愈、植物的产量等等，数量性状对人类生存和发展的重要性不言而喻。一般说来，数量性状的变异呈间断离散型或者连续型，遗传机理也相对的会比较复杂，而且容易受到环境因素的影响，表现形式为连续变异，控制数量性状的基因在特定的条件下表达，表型与基因型之间关系不明确，因此对于数量性状遗传研究进展并不深入。多个相关性状的联合定位，不但能提高检测的功效和精度，而且能区分紧密连锁与一因多效 [1]。在实际应用上，把植物种子遗传发育的分子机理与其分子育种相结合，将有助于提高产量和外观品质。然而随着生物技术的进步，大量基因标记的出现和数量性状的上位性检测，使得遗传模型中样本容量远大于变量的个数，这就使得常用的参数估计方法不再有效，急需一种新的参数估计方法。而目前的人类遗传学中的关联分析方法均为单标记分析，即一次在模型中只分析一个基因标记，而忽略了基因之间的联系，所以，多性状联合分析方法的研究非常重要。1.2多性状联合统计方法的关联分析简介基于连锁不平衡（Linkage disequilibrium，LD）的全基因组关联分析(genome-wide association study, GWAS)[2][3]是现今解析植物数量性状基因型的主要方法（Price等，2006；Mackay和Powell，2007；Zhu等，2008）。关联分析，将生物个体的表现型与基因型相联系，利用统计方法在整个基因组上搜索控制数量性状的基因位点，以揭示复杂性状的遗传基础。多性状联合分析方法具有许多其他统计方法不具备的优点，不仅可以增加效应估计的准度和精度、提高统计效应估计功效、显著增强基因的发现能力、降低科研资金和显著减少计算时间，还对一因多效基因的检测能力和效应估计的准确度也有着显著的进步。全基因组关联分析是指在生物全基因组范围内找出海量变异序列，即单核苷酸多态性，然后从中筛选出与性状相关的单核苷酸多态性。再根据其在基因组中的位置，和连锁不平衡推测可能的易感基因，并对这些基因和影响的性状的关系进行验证及相关的功能研究[4]。1.3国内外研究背景Jiang C和Zeng Z-B [5]通过研究发现，通过多个相关性状的联合数量性状基因座定位，可以提高数量性状基因座检测的功效和精度，也可以由此区分紧密连锁和一因多效。Henshall等提出logistic用回归的分析方法，用于预先选择的基因型和相关的表型数据定位分析，该研究表明这个方法对数量性状基因座位置和效应的估计渐进无偏的，在多性状联合分析时统计功效有所增加。Xu等对多个不连续的性状进行联合数量性状基因座分析，该方法将多变量的易感性作为缺失值，使之可用最大期望算法法计算，结果表明该方法具有较高的统计功效。Banerjee S和Yandell BS [6]等在2007年提出用贝叶斯方法进行相关性状的联合定位，该方法通过模拟相关表型和不相关的表型2个模型，以马尔可夫链蒙特卡罗算法进行运算。Xu等应用贝叶斯方法分析多性状联合数量性状基因座定位，该方法对模型的选择采用参数压缩的方法估计所有标记区间内的遗传效应。肖静和胡治球 [7]等提出多个相关数量性状的主基因的联合分析方法，得到了不管是控制多个性状还是一个性状的表达，多性状联合分析方法均要比单性状分析更优秀的结论，联合分析不仅可以提高检测效率，还可以增加准确度和精确度。随着多性状联合方法关联分析在国内外的发展，已经取得了比较显著的成绩，成为一种比较常用的生物统计方法，具有广泛的应用前景。2 理论及方法2.1似然比检验[8]似然比(likelihood ratio, LR) 由奈曼（Neyman）和皮尔逊（E.Pearson）在1928年提出，是一种应用较广的检验方法，在假设检验中的地位有如MLE在点估计中的地位。似然比检验应该看作是Neymann-Pearson检验的推广，和极大似然估计没有直接联系，似然比检验比一般的假设检验效果要好，但遗憾的是，该似然比检验统计量在一般场合至今没有统一的精确分布形式但在一般的条件下有一个统一的渐进分布。定义1.1设x1,x2,,xn为来自密度函数为p（x；θ），θ∈Θ的样本，考虑如下检验问题H0θ∈Θ0 vs H1θ∈Θ1=Θ-Θ0 (1.1)令Λ（x1,,xn）=????????????????????Θ????(????1,,????????;????)????????????????????Θ0????(????1,,????????;????), (1.2)则我们称统计量Λ（x1,,xn）为假设（1.1）的似然比，有时也称广义似然比。(1.2)式的Λ（x1,,xn）=????(????1,,????????;????)????(????1,,????????;????0)， (1.3)其中????表示在全参数空间Θ上θ的最大似然估计，????0表示在子参数空间Θ0上θ的最大似然估计。也就是说，Λ（x1,,xn）的分子表示没有假设时的似然函数最大值，分母表示在原假设成立条件下的似然函数最大值，不难看出，如果Λ（x1,,xn）的值很大，则说明θ∈Θ0的可能性要比θ∈Θ1的可能性小，于是，我们有理由认为H0不成立。这样，我们有如下的似然比检验。定义1.2当采用（1.3）式的似然比统计量Λ（x1,,xn）作为检验问题（1.1）的检验统计量，且取其拒绝域为W={Λ（x1,,xn）≥c}，其中临界值c满足Pθ（Λ（x1,,xn）≥c）≤α，?θ?Θ0， (1.4)则称此检验为显著性水平α的似然比检验（likelihood ratio test），简记为LRT。特点计算似然比只涉及特异度和灵敏度，不会受到患病率的影响。2.2多性状定位多性状关联分析，它的研究对象就很有可能是全基因组关联分析或连锁分析的群体。假设一个群体具有m个标记，s个性状，有n个个体。那么联合分析的线性模型能够表示为:????1????=????=1????????????????????1????+????=1????????????????????1????+????1????????????????=????=1????????????????????????????+????=1????????????????????????????+???????????? （2.1）式中yli是n×1的向量，表示第l个性状的表型；zij表示第j个系统环境第i个体的指示变量；rlj表示第j个系统环境第l个性状效应；xij表示第j标记处的第i个体指示变量；βlj表示第j个标记的第l个性状的遗传效应；eli表示服从eil~N（0，????）的多元正态分布的剩余误差；????表示剩余误差协方差矩阵。一般在估计之前都会对数据进行预处理，可以去除系统环境效应，所以可以将方程简化为????1????=????10+????=1????????????????????1????+????1????????2????=????20+????=1????????????????????2????+????2????????????????=????????0+????=1????????????????????????????+???????????? (2.2)矩阵形式表示为Y=Xβ+e (2.3)式中Y表示为多个性状的表型数据矩阵。Y具体表示为Y=????11?????????1???????1?????????????????=????(1)????(2)????(????) (2.4)β表示为标记效应系数矩阵。β具体表示为β=????10?????????0???????1?????????????????=????(1)????(2)????(????) (2.5)X表示为标记指示变量矩阵。X具体表示为X=????10?????1???????????????0????????????? (2.6)ε表示为每个变量的残差矩阵。ε具体表示为e=????11?????????1???????1?????????????????=????(1)????(2)????(????)=????(1)????(2)????(????) (2.7)残差矩阵的数学期望和方差分别是Ee=E????1????2????????=0 (2.8) Vare=???? (2.9)2.3 LASSO算法[9]LASSO（Least absolute shrinkage and selection operator）由1996年Robert Tibshirani首次提出。对于一般的回归模型，形式如下y????=????=1????????????????????????+???????? (3.1)式中 Yi表示因变量m表示自变量个数Xij表示自变量βj表示回归系数ei表示随机误差它的矩阵形式表示为y=Xβ+e (3.2)式中Y是n×1矩阵，表示因变量；X是n×m矩阵，表示自变量矩阵；β是m×1矩阵，表示回归系数向量；e是n×1矩阵，表示随机误差向量。通过最小化(3.3)求解方程Q=????=1????（?????????????=1????????????????????????）2+λ（????=1????????????) (3.3)式中λ表示调节参数，可通过交叉验证获得为了偏优化βu（u=1，2，，m），在βu=????????的梯度通过下式计算?Q?????????|????=????=?1????????=1?????????????????????????????=1????????????????????????+???????????? (3.4)当????????>0，????????<0或????????=0分别的处理（3.4）这个方程。坐标更新有一个简单的算法，遵循以下的方式[10][11]:????????=11+????????（1ni=1nxiuyi?yiu，λ） (3.5)其中????????=????=1,????≠????????????????????????????是因变量的拟合值。b表示1????????=1????????????????????（yi?yi（u））为在u点的最简单的最小二乘估计系数。S（b，β）的值按下式计算signb?????β+=????????? ???????? ????>0 ???????????? ????<????????+???? ???????? ????<0 ???????????? ????<????0 ???????? ????≥???? (3.6)最终可根据方程(3.6)估计出????????。优点最小角回归算法在运算步骤上更简便，其计算的消耗与通常的最小二乘估计是一样的。它能更有效地解决Lasso的计算问题。2.4 多元最小二乘法设给定数据点（x11,x12,,x1n,y1）, （x21,x22,,x2n,y2）,.., （xn1,xn2,,xnn,yn）大致在一直线上，设回归函数为y=b0+b1x1+b2x2++bnxn，则多元线性回归最小二乘问题即求Q????0,..,????????=????=1????[?????????(????0+????1????1????++????????????????????)]2 (4.1)达到最小。对上式的所有参数求偏导数并令偏导数为0，可以得到多元函数的极值?Q?????????=?2????=1?????????????????0+????1????1????++????????????????????????????????=0 (4.2)求解方程，记????????????=????=1?????????????????????????????????????，Sky=????=1????(?????????????????????)(?????????????)，可得????0=?????????1????1??????????????????????11????1+????12????2+????1????????????=????1????????????1????1+????????2????2+????????????????????=???????????? (4.3)即为求解多元线性回归最小二乘解的正规方程组。另一方面，这个式子也能看作一个矩阵Y=BX+ε (4.4)X′Y=X′XB+X′ε (4.5)B=(X′X)-1X′Y (4.6)估计量B称为最小二乘估计。多元最小二乘法的优点第一，最小二乘法同时适用于线性和非线性的数学模型；第二，最小二乘法估计只与算术平均值这个统计量有关，这个统计量较其他统计量重要；第三，最小二乘法的运用范围很广，这也使得它成为获得唯一估值的标准方法。2.5 偏最小二乘（partial square method）偏最小二乘法[12]是一种最新的多元统计方法，它的目标是把误差的平方和最小化，主要用于研究多个因变量对多个自变量的回归模型，尤其当各变量内部具有高度线性相关时，偏最小二乘回归法会更加的有效。此外，偏最小二乘回归也能更好地解决了变量个数多于样本个数等问题。它是由伍德（S.Wold）和阿巴诺（C.Albano）等人在1983年首次提出。在最近几十年内，它在理论、方法和应用方面都得到迅速的发展，被称为第二代回归方法，是多元统计数据分析中的一个巨大的飞跃。偏最小二乘回归相当于是一种把多元线性回归分析、典型相关分析和主成分分析结合起来的一种回归分析，偏最小二乘法同时有这3种分析方法的优点。偏最小二乘回归虽然和主成分分析法一样，都会提取数据大部分信息，但主成分分析法只有自变量矩阵，没有偏最小二乘所具有的“响应”矩阵，因此也不具备偏最小二乘所具备的预测功能。以yi=(y1i y2i yni)T代表第i（i=1,2,q）个性状的观测值，则Y=（y1 y2 yq）n*q为表型性状观测矩阵。以bi=（bi0 bi1 bim）T表示第i个性状的m+1个效应值，其中βi0表示截距项，B=（b1 b2 bq）（m+1）*q效应矩阵。ε=（εij）n*q=（ε1 ε2 εq）为随机误差效应矩阵。则线性模型表示为Y=XB+ε (5.1)其中X表示n行、（m+1）列的系数矩阵，除第一列均取1之外。该模型在n?m且q≥2时为多变量的超饱和模型，因为X矩阵是奇异的，XTX为不可逆矩阵，无法进行模型参数的最小二乘估计，因此，B=（XTX）-1XTY无解，于是采用非线性迭代偏最小二乘法。特点（1）适用于自变量之间具有严重的多重相关性的条件；（2）适用于在变量个数大于样本点个数的条件；（3）在回归模型中，每一个自变量的回归系数将变得更容易解释。（4）在最终模型中可以包含原有的所有自变量；2.6 R语言简介R语言[13]是用于统计分析、绘图的实用软件。它是一个免费、自由、源代码开放的软件，基于S语言。R最早由新西兰奥克兰大学的Ross和Robert开发，如今由“R开发核心团队”担任开发工作。R是开源程序，能够自由地下载使用，也有许多已经编译好的可执行文件可以直接下载使用，还能够在多个平台下运行，操作方式是用命令行，有人开发一些比较好用的图形用户界面，如Rstudio等。R是一套具完整的计算、制图和数据处理于一体的软件系统。它的功能如下大量数据的存储和处理；优秀的矩阵、向量运算工具；统计分析完整而又使用方便；统计制图功能十分强大。其特点为完善、简洁而效率极高的程序设计语言；面向对象的统计编程语言。R/chemometrics是基于R语言的多元统计分析的软件包，有很好的实用性。本文所用到的pls1_nipals和pls2_nipals这两个程序是由K. Varmuza和P. Filzmoser编写的，适用于偏最小二乘的统计计算，分别是单性状独立分析和多性状联合分析。3 仿真实验及结果分析3.1 实验结果本文所使用的样本总量为200。在长度为2000cm的染色体片段上，每隔10cm设置一个标记，总共有201个用于仿真的标记均匀分布在该染色体片段上。在每个样本中需要估计的数量性状效应共有3个，他们都位于标记位置上。我们设定总体均值为10，残差方差为10，用于生成数量形状表型观测值。在本次实验中，我们用相同的参数设定来生成200组数据集，然后再借助于Rstudio中的程序包R/chemometrics中的pls2_nipals 和pls1_nipals，分别对这组数据进行偏最小二乘，再对得到的一组参数估计值通过编写lod程序进行假设检验，最后对所得结果进行处理，剔除p大于0.05的数据，整理剩余的数据，可以得到200个效应的P值均值与对应的标准差、显著的p的计数、lod及对应的标准差、效应值。表1 多性状联合分析结果1geneP值（标准差）显著p计数Lod（标准差）效应值V302.730E-3（7.890E-3）1456.405(8.661)1.774E-1V1004.746E-3（1.087E-2）1325.302(7.572)4.146E-3V1802.890E-3 (9.853E-3)1376.374(9.098)6.285E-1表2 多性状联合分析结果2geneP值（标准差）显著p计数Lod（标准差）效应值V602.760E-3（8.560E-3）1609.369(12.786)3.114E-1V1303.430E-3（8.697E-3）14910.318(14.824)2.744E-2V1802.455E-3（6.925E-3）1529.788(12.950)8.906E-1表3 单性状独立分析结果1geneP值（标准差）显著p计数Lod（标准差）效应值V305.304E-3（1.121E-2）1244.532(7.140)1.815E-1V1004.240E-3（9.571E-3）1354.667(6.980)1.718E-2V1804.065E-3（8.770E-3）1444.499(6.293)6.398E-1表4 单性状独立分析结果2geneP值（标准差）显著p计数Lod（标准差）效应值V603.290E-3（8.540E-3）1477.516(9.716)3.167E-1V1304.118E-3（1.087E-2）1518.153(12.074)2.832E-2V1802.660E-3（7.884E-3）1478.972(11.448)8.840E-13.2结果分析3.2.1计算时间使用pls2_nipals 和pls1_nipals这两个程序分别对40000个基因数据进行参数估计所需时间都为10秒左右，使用lod程序进行假设检验所需时间都为20分钟左右，所以，多性状联合分析与单性状独立分析所需时间相差不大，所消耗的时间都比较短。3.2.2显著性检验p值是一个用于判断假设检验的参数，由R·A·Fisher创立。是体现结果是否可信的一个重要指标。p值越小，就能体现结果具有更好的显著性。反之，p值越大，我们就越不能拒绝原假设，也就表示犯第一类错误的可能性越大。第一类错误也就是拒绝了实际成立的假设，即弃真错误。对比分析表1和表3的p值，V30分别为2.730E-3和5.304E-3；V100分别为4.746E-3和4.240E-3；V180分别为2.890E-3和4.065E-3。对比分析表2和表4的p值，V60分别为2.760E-3和3.290E-3；V130分别为3.430E-3和4.118E-3；V180分别为2.455E-3和2.660E-3。在P值这项指标上，除分析结果1中V100外，其余5个分析结果均是多性状联合分析小于单性状独立分析，所以，总的来说，我们可以得出多性状联合分析比单性状独立分析更显著的结论。3.2.3 连锁关系Lod值是为了确定两个基因座位置是否紧密，是否可能一起遗传。为了判断两对基因之间是否存在着连锁关系，一般都会要求或然比大于10001，即lod>3。也就是说，lod越大，表示基因座的位置越紧密，连锁的概率越高。对比分析表1和表3的lod值，V30分别为6.405和4.532；V100分别为5.302和4.667；V180分别为6.374和4.499。对比分析表2和表4的lod值，V60分别为9.369和7.516；V130分别为10.318和8.153；V180分别为9.788和8.972。很明显，在lod值这项指标上，所有的分析结果都是多性状联合分析均要大于单性状独立分析，所以，我们可以得出结论多性状联合分析比单性状独立分析连锁率更高，基因座的位置也较之更加紧密。3.3讨论此次研究的结果与预期的结果基本大致相符，在对方法的选择和使用分析方法计算得出的结果也比较满意，计算速度也很快。此次研究通过对p值和lod值得比较，可以发现多性状联合分析比单性状分析具有更好的可信度和更高的连锁性，这个结果与预计结果相同。在耗时方面原以为多性状分析会多于单性状，结果也表明，它们所需要的时间相差不大。但是多性状联合分析的检测功效等会比单性状分析差，这个结果与预计结果出入较大，可能原因是偏最小二乘法会想主成分分析一样提取大部分数据信息进行分析，从而会损失部分（一般低于15%）信息，这损失的部分信息影响到了检测功效。但整体来说，多性状联合分析的结果要优于单性状独立分析的结果。致谢参考文献[1] Jiang C, Zeng Z-B, Multiple trait analysis of genetic mapping for quantitative trait loci [J] .Genetics,1995, 140 (3):1111-1127.[2] 严卫丽，复杂疾病全基因组关联研究进展——遗传统计分析[J]. 遗传, 2008, 30 (5): 543-549.[3] Spielman R S, Mc Ginnis R E, Ewens W J, Transmission test for linkage disequilibrium: the insulin gene region and insulin-dependent diabetes mellitus (IDDM) [J]. American journal of human genetics, 1993, 52 (3): 506.[4] Risch N, Merikangas K. The future of genetic studies of complex human diseases. Science, 1996,273:1516-1517[5] Jiang C, Zeng Z-B, Multiple trait analysis of genetic mapping for quantitative trait loci [J] .Genetics,1995, 140 (3):1111-1127[6]Yandell BS,Mehta T,Banerjee S,Shriner D,Venkataraman R,Moon JY,Neely WW,Wu H,Von Smith R,Yi NJ,R/qtlbim:QTL with Bayesian interval mapping in experimental crosses. Bioinformatics,2007,23(5):641-3[7]肖静，胡治球，汤在祥，隋炯明，李欣，徐辰武.多个相关数量性状主基因的联合分析方法[J].中国农业科学，2005，38(9):171724[8] 茆诗松 程依明 濮晓龙.概率论与数理统计教程[M].高等教育出版社.2011.2[9]李红旺. 基于LASSO算法的多性状QTL条件定位方法[D].上海上海交通大学，201423.[10] Donoho D L, Johnstone I M, Ideal spatial adaptation by wavelet shrinkage[J]. Biometrika, 1994, 81: 425-455. [11] Friedman J,Hastie T, Hoefling H,et al, Pathwise coordinate optimization[J]. The Annals of Applied Statistics,2007 2(1): 302-332[12]胡文明. 作物复杂性状QTL定位相关的几个问题的探讨[D].扬州扬州大学，201462.[13]Team RDC.R:A language and environment for statistical computing.ISBN 3-900051-07-0.R Foundation for Statistical Computing.Vienna,Austria,2013.url://www.R-project.org.2005
目录
摘要3
关键词3
Abstract3
Key words3
1、绪论3
1.1研究意义 3
1.2多性状联合统计方法的关联分析简介3
1.3国内外研究背景 4
2.2理论及方法 4
2.1似然比检验4
2.2多性状定位 5
2.3 LASSO算法 6
2.4 多元最小二乘法 7
2.5 偏最小二乘 7
2.6 R语言简介 8
3、仿真实验及结果分析 8
3.1 实验结果 8
3.2结果分析9 3.3讨论 10
致谢10
参考文献10
基于多性状联合统计方法的关联分析
引言
human genetics. Carry on gene mapping can be linked to the genotype and individual traits. This paper introduces the application of statistical methods in the field of biological multi character. Least squares, partial least squares and LASSO are commonly used statistical methods, can be used to study the association between traits. At the same time, the results were compared with single character analysis.1 绪论1.1研究意义有很多生物遗传方面的重要性状都是数量性状，例如人类疾病的治愈、植物的产量等等，数量性状对人类生存和发展的重要性不言而喻。一般说来，数量性状的变异呈间断离散型或者连续型，遗传机理也相对的会比较复杂，而且容易受到环境因素的影响，表现形式为连续变异，控制数量性状的基因在特定的条件下表达，表型与基因型之间关系不明确，因此对于数量性状遗传研究进展并不深入。多个相关性状的联合定位，不但能提高检测的功效和精度，而且能区分紧密连锁与一因多效 [1]。在实际应用上，把植物种子遗传发育的分子机理与其分子育种相结合，将有助于提高产量和外观品质。然而随着生物技术的进步，大量基因标记的出现和数量性状的上位性检测，使得遗传模型中样本容量远大于变量的个数，这就使得常用的参数估计方法不再有效，急需一种新的参数估计方法。而目前的人类遗传学中的关联分析方法均为单标记分析，即一次在模型中只分析一个基因标记，而忽略了基因之间的联系，所以，多性状联合分析方法的研究非常重要。1.2多性状联合统计方法的关联分析简介基于连锁不平衡（Linkage disequilibrium，LD）的全基因组关联分析(genome-wide association study, GWAS)[2][3]是现今解析植物数量性状基因型的主要方法（Price等，2006；Mackay和Powell，2007；Zhu等，2008）。关联分析，将生物个体的表现型与基因型相联系，利用统计方法在整个基因组上搜索控制数量性状的基因位点，以揭示复杂性状的遗传基础。多性状联合分析方法具有许多其他统计方法不具备的优点，不仅可以增加效应估计的准度和精度、提高统计效应估计功效、显著增强基因的发现能力、降低科研资金和显著减少计算时间，还对一因多效基因的检测能力和效应估计的准确度也有着显著的进步。全基因组关联分析是指在生物全基因组范围内找出海量变异序列，即单核苷酸多态性，然后从中筛选出与性状相关的单核苷酸多态性。再根据其在基因组中的位置，和连锁不平衡推测可能的易感基因，并对这些基因和影响的性状的关系进行验证及相关的功能研究[4]。1.3国内外研究背景Jiang C和Zeng Z-B [5]通过研究发现，通过多个相关性状的联合数量性状基因座定位，可以提高数量性状基因座检测的功效和精度，也可以由此区分紧密连锁和一因多效。Henshall等提出logistic用回归的分析方法，用于预先选择的基因型和相关的表型数据定位分析，该研究表明这个方法对数量性状基因座位置和效应的估计渐进无偏的，在多性状联合分析时统计功效有所增加。Xu等对多个不连续的性状进行联合数量性状基因座分析，该方法将多变量的易感性作为缺失值，使之可用最大期望算法法计算，结果表明该方法具有较高的统计功效。Banerjee S和Yandell BS [6]等在2007年提出用贝叶斯方法进行相关性状的联合定位，该方法通过模拟相关表型和不相关的表型2个模型，以马尔可夫链蒙特卡罗算法进行运算。Xu等应用贝叶斯方法分析多性状联合数量性状基因座定位，该方法对模型的选择采用参数压缩的方法估计所有标记区间内的遗传效应。肖静和胡治球 [7]等提出多个相关数量性状的主基因的联合分析方法，得到了不管是控制多个性状还是一个性状的表达，多性状联合分析方法均要比单性状分析更优秀的结论，联合分析不仅可以提高检测效率，还可以增加准确度和精确度。随着多性状联合方法关联分析在国内外的发展，已经取得了比较显著的成绩，成为一种比较常用的生物统计方法，具有广泛的应用前景。2 理论及方法2.1似然比检验[8]似然比(likelihood ratio, LR) 由奈曼（Neyman）和皮尔逊（E.Pearson）在1928年提出，是一种应用较广的检验方法，在假设检验中的地位有如MLE在点估计中的地位。似然比检验应该看作是Neymann-Pearson检验的推广，和极大似然估计没有直接联系，似然比检验比一般的假设检验效果要好，但遗憾的是，该似然比检验统计量在一般场合至今没有统一的精确分布形式但在一般的条件下有一个统一的渐进分布。定义1.1设x1,x2,,xn为来自密度函数为p（x；θ），θ∈Θ的样本，考虑如下检验问题H0θ∈Θ0 vs H1θ∈Θ1=Θ-Θ0 (1.1)令Λ（x1,,xn）=????????????????????Θ????(????1,,????????;????)????????????????????Θ0????(????1,,????????;????), (1.2)则我们称统计量Λ（x1,,xn）为假设（1.1）的似然比，有时也称广义似然比。(1.2)式的Λ（x1,,xn）=????(????1,,????????;????)????(????1,,????????;????0)， (1.3)其中????表示在全参数空间Θ上θ的最大似然估计，????0表示在子参数空间Θ0上θ的最大似然估计。也就是说，Λ（x1,,xn）的分子表示没有假设时的似然函数最大值，分母表示在原假设成立条件下的似然函数最大值，不难看出，如果Λ（x1,,xn）的值很大，则说明θ∈Θ0的可能性要比θ∈Θ1的可能性小，于是，我们有理由认为H0不成立。这样，我们有如下的似然比检验。定义1.2当采用（1.3）式的似然比统计量Λ（x1,,xn）作为检验问题（1.1）的检验统计量，且取其拒绝域为W={Λ（x1,,xn）≥c}，其中临界值c满足Pθ（Λ（x1,,xn）≥c）≤α，?θ?Θ0， (1.4)则称此检验为显著性水平α的似然比检验（likelihood ratio test），简记为LRT。特点计算似然比只涉及特异度和灵敏度，不会受到患病率的影响。2.2多性状定位多性状关联分析，它的研究对象就很有可能是全基因组关联分析或连锁分析的群体。假设一个群体具有m个标记，s个性状，有n个个体。那么联合分析的线性模型能够表示为:????1????=????=1????????????????????1????+????=1????????????????????1????+????1????????????????=????=1????????????????????????????+????=1????????????????????????????+???????????? （2.1）式中yli是n×1的向量，表示第l个性状的表型；zij表示第j个系统环境第i个体的指示变量；rlj表示第j个系统环境第l个性状效应；xij表示第j标记处的第i个体指示变量；βlj表示第j个标记的第l个性状的遗传效应；eli表示服从eil~N（0，????）的多元正态分布的剩余误差；????表示剩余误差协方差矩阵。一般在估计之前都会对数据进行预处理，可以去除系统环境效应，所以可以将方程简化为????1????=????10+????=1????????????????????1????+????1????????2????=????20+????=1????????????????????2????+????2????????????????=????????0+????=1????????????????????????????+???????????? (2.2)矩阵形式表示为Y=Xβ+e (2.3)式中Y表示为多个性状的表型数据矩阵。Y具体表示为Y=????11?????????1???????1?????????????????=????(1)????(2)????(????) (2.4)β表示为标记效应系数矩阵。β具体表示为β=????10?????????0???????1?????????????????=????(1)????(2)????(????) (2.5)X表示为标记指示变量矩阵。X具体表示为X=????10?????1???????????????0????????????? (2.6)ε表示为每个变量的残差矩阵。ε具体表示为e=????11?????????1???????1?????????????????=????(1)????(2)????(????)=????(1)????(2)????(????) (2.7)残差矩阵的数学期望和方差分别是Ee=E????1????2????????=0 (2.8) Vare=???? (2.9)2.3 LASSO算法[9]LASSO（Least absolute shrinkage and selection operator）由1996年Robert Tibshirani首次提出。对于一般的回归模型，形式如下y????=????=1????????????????????????+???????? (3.1)式中 Yi表示因变量m表示自变量个数Xij表示自变量βj表示回归系数ei表示随机误差它的矩阵形式表示为y=Xβ+e (3.2)式中Y是n×1矩阵，表示因变量；X是n×m矩阵，表示自变量矩阵；β是m×1矩阵，表示回归系数向量；e是n×1矩阵，表示随机误差向量。通过最小化(3.3)求解方程Q=????=1????（?????????????=1????????????????????????）2+λ（????=1????????????) (3.3)式中λ表示调节参数，可通过交叉验证获得为了偏优化βu（u=1，2，，m），在βu=????????的梯度通过下式计算?Q?????????|????=????=?1????????=1?????????????????????????????=1????????????????????????+???????????? (3.4)当????????>0，????????<0或????????=0分别的处理（3.4）这个方程。坐标更新有一个简单的算法，遵循以下的方式[10][11]:????????=11+????????（1ni=1nxiuyi?yiu，λ） (3.5)其中????????=????=1,????≠????????????????????????????是因变量的拟合值。b表示1????????=1????????????????????（yi?yi（u））为在u点的最简单的最小二乘估计系数。S（b，β）的值按下式计算signb?????β+=????????? ???????? ????>0 ???????????? ????<????????+???? ???????? ????<0 ???????????? ????<????0 ???????? ????≥???? (3.6)最终可根据方程(3.6)估计出????????。优点最小角回归算法在运算步骤上更简便，其计算的消耗与通常的最小二乘估计是一样的。它能更有效地解决Lasso的计算问题。2.4 多元最小二乘法设给定数据点（x11,x12,,x1n,y1）, （x21,x22,,x2n,y2）,.., （xn1,xn2,,xnn,yn）大致在一直线上，设回归函数为y=b0+b1x1+b2x2++bnxn，则多元线性回归最小二乘问题即求Q????0,..,????????=????=1????[?????????(????0+????1????1????++????????????????????)]2 (4.1)达到最小。对上式的所有参数求偏导数并令偏导数为0，可以得到多元函数的极值?Q?????????=?2????=1?????????????????0+????1????1????++????????????????????????????????=0 (4.2)求解方程，记????????????=????=1?????????????????????????????????????，Sky=????=1????(?????????????????????)(?????????????)，可得????0=?????????1????1??????????????????????11????1+????12????2+????1????????????=????1????????????1????1+????????2????2+????????????????????=???????????? (4.3)即为求解多元线性回归最小二乘解的正规方程组。另一方面，这个式子也能看作一个矩阵Y=BX+ε (4.4)X′Y=X′XB+X′ε (4.5)B=(X′X)-1X′Y (4.6)估计量B称为最小二乘估计。多元最小二乘法的优点第一，最小二乘法同时适用于线性和非线性的数学模型；第二，最小二乘法估计只与算术平均值这个统计量有关，这个统计量较其他统计量重要；第三，最小二乘法的运用范围很广，这也使得它成为获得唯一估值的标准方法。2.5 偏最小二乘（partial square method）偏最小二乘法[12]是一种最新的多元统计方法，它的目标是把误差的平方和最小化，主要用于研究多个因变量对多个自变量的回归模型，尤其当各变量内部具有高度线性相关时，偏最小二乘回归法会更加的有效。此外，偏最小二乘回归也能更好地解决了变量个数多于样本个数等问题。它是由伍德（S.Wold）和阿巴诺（C.Albano）等人在1983年首次提出。在最近几十年内，它在理论、方法和应用方面都得到迅速的发展，被称为第二代回归方法，是多元统计数据分析中的一个巨大的飞跃。偏最小二乘回归相当于是一种把多元线性回归分析、典型相关分析和主成分分析结合起来的一种回归分析，偏最小二乘法同时有这3种分析方法的优点。偏最小二乘回归虽然和主成分分析法一样，都会提取数据大部分信息，但主成分分析法只有自变量矩阵，没有偏最小二乘所具有的“响应”矩阵，因此也不具备偏最小二乘所具备的预测功能。以yi=(y1i y2i yni)T代表第i（i=1,2,q）个性状的观测值，则Y=（y1 y2 yq）n*q为表型性状观测矩阵。以bi=（bi0 bi1 bim）T表示第i个性状的m+1个效应值，其中βi0表示截距项，B=（b1 b2 bq）（m+1）*q效应矩阵。ε=（εij）n*q=（ε1 ε2 εq）为随机误差效应矩阵。则线性模型表示为Y=XB+ε (5.1)其中X表示n行、（m+1）列的系数矩阵，除第一列均取1之外。该模型在n?m且q≥2时为多变量的超饱和模型，因为X矩阵是奇异的，XTX为不可逆矩阵，无法进行模型参数的最小二乘估计，因此，B=（XTX）-1XTY无解，于是采用非线性迭代偏最小二乘法。特点（1）适用于自变量之间具有严重的多重相关性的条件；（2）适用于在变量个数大于样本点个数的条件；（3）在回归模型中，每一个自变量的回归系数将变得更容易解释。（4）在最终模型中可以包含原有的所有自变量；2.6 R语言简介R语言[13]是用于统计分析、绘图的实用软件。它是一个免费、自由、源代码开放的软件，基于S语言。R最早由新西兰奥克兰大学的Ross和Robert开发，如今由“R开发核心团队”担任开发工作。R是开源程序，能够自由地下载使用，也有许多已经编译好的可执行文件可以直接下载使用，还能够在多个平台下运行，操作方式是用命令行，有人开发一些比较好用的图形用户界面，如Rstudio等。R是一套具完整的计算、制图和数据处理于一体的软件系统。它的功能如下大量数据的存储和处理；优秀的矩阵、向量运算工具；统计分析完整而又使用方便；统计制图功能十分强大。其特点为完善、简洁而效率极高的程序设计语言；面向对象的统计编程语言。R/chemometrics是基于R语言的多元统计分析的软件包，有很好的实用性。本文所用到的pls1_nipals和pls2_nipals这两个程序是由K. Varmuza和P. Filzmoser编写的，适用于偏最小二乘的统计计算，分别是单性状独立分析和多性状联合分析。3 仿真实验及结果分析3.1 实验结果本文所使用的样本总量为200。在长度为2000cm的染色体片段上，每隔10cm设置一个标记，总共有201个用于仿真的标记均匀分布在该染色体片段上。在每个样本中需要估计的数量性状效应共有3个，他们都位于标记位置上。我们设定总体均值为10，残差方差为10，用于生成数量形状表型观测值。在本次实验中，我们用相同的参数设定来生成200组数据集，然后再借助于Rstudio中的程序包R/chemometrics中的pls2_nipals 和pls1_nipals，分别对这组数据进行偏最小二乘，再对得到的一组参数估计值通过编写lod程序进行假设检验，最后对所得结果进行处理，剔除p大于0.05的数据，整理剩余的数据，可以得到200个效应的P值均值与对应的标准差、显著的p的计数、lod及对应的标准差、效应值。表1 多性状联合分析结果1geneP值（标准差）显著p计数Lod（标准差）效应值V302.730E-3（7.890E-3）1456.405(8.661)1.774E-1V1004.746E-3（1.087E-2）1325.302(7.572)4.146E-3V1802.890E-3 (9.853E-3)1376.374(9.098)6.285E-1表2 多性状联合分析结果2geneP值（标准差）显著p计数Lod（标准差）效应值V602.760E-3（8.560E-3）1609.369(12.786)3.114E-1V1303.430E-3（8.697E-3）14910.318(14.824)2.744E-2V1802.455E-3（6.925E-3）1529.788(12.950)8.906E-1表3 单性状独立分析结果1geneP值（标准差）显著p计数Lod（标准差）效应值V305.304E-3（1.121E-2）1244.532(7.140)1.815E-1V1004.240E-3（9.571E-3）1354.667(6.980)1.718E-2V1804.065E-3（8.770E-3）1444.499(6.293)6.398E-1表4 单性状独立分析结果2geneP值（标准差）显著p计数Lod（标准差）效应值V603.290E-3（8.540E-3）1477.516(9.716)3.167E-1V1304.118E-3（1.087E-2）1518.153(12.074)2.832E-2V1802.660E-3（7.884E-3）1478.972(11.448)8.840E-13.2结果分析3.2.1计算时间使用pls2_nipals 和pls1_nipals这两个程序分别对40000个基因数据进行参数估计所需时间都为10秒左右，使用lod程序进行假设检验所需时间都为20分钟左右，所以，多性状联合分析与单性状独立分析所需时间相差不大，所消耗的时间都比较短。3.2.2显著性检验p值是一个用于判断假设检验的参数，由R·A·Fisher创立。是体现结果是否可信的一个重要指标。p值越小，就能体现结果具有更好的显著性。反之，p值越大，我们就越不能拒绝原假设，也就表示犯第一类错误的可能性越大。第一类错误也就是拒绝了实际成立的假设，即弃真错误。对比分析表1和表3的p值，V30分别为2.730E-3和5.304E-3；V100分别为4.746E-3和4.240E-3；V180分别为2.890E-3和4.065E-3。对比分析表2和表4的p值，V60分别为2.760E-3和3.290E-3；V130分别为3.430E-3和4.118E-3；V180分别为2.455E-3和2.660E-3。在P值这项指标上，除分析结果1中V100外，其余5个分析结果均是多性状联合分析小于单性状独立分析，所以，总的来说，我们可以得出多性状联合分析比单性状独立分析更显著的结论。3.2.3 连锁关系Lod值是为了确定两个基因座位置是否紧密，是否可能一起遗传。为了判断两对基因之间是否存在着连锁关系，一般都会要求或然比大于10001，即lod>3。也就是说，lod越大，表示基因座的位置越紧密，连锁的概率越高。对比分析表1和表3的lod值，V30分别为6.405和4.532；V100分别为5.302和4.667；V180分别为6.374和4.499。对比分析表2和表4的lod值，V60分别为9.369和7.516；V130分别为10.318和8.153；V180分别为9.788和8.972。很明显，在lod值这项指标上，所有的分析结果都是多性状联合分析均要大于单性状独立分析，所以，我们可以得出结论多性状联合分析比单性状独立分析连锁率更高，基因座的位置也较之更加紧密。3.3讨论此次研究的结果与预期的结果基本大致相符，在对方法的选择和使用分析方法计算得出的结果也比较满意，计算速度也很快。此次研究通过对p值和lod值得比较，可以发现多性状联合分析比单性状分析具有更好的可信度和更高的连锁性，这个结果与预计结果相同。在耗时方面原以为多性状分析会多于单性状，结果也表明，它们所需要的时间相差不大。但是多性状联合分析的检测功效等会比单性状分析差，这个结果与预计结果出入较大，可能原因是偏最小二乘法会想主成分分析一样提取大部分数据信息进行分析，从而会损失部分（一般低于15%）信息，这损失的部分信息影响到了检测功效。但整体来说，多性状联合分析的结果要优于单性状独立分析的结果。致谢参考文献[1] Jiang C, Zeng Z-B, Multiple trait analysis of genetic mapping for quantitative trait loci [J] .Genetics,1995, 140 (3):1111-1127.[2] 严卫丽，复杂疾病全基因组关联研究进展——遗传统计分析[J]. 遗传, 2008, 30 (5): 543-549.[3] Spielman R S, Mc Ginnis R E, Ewens W J, Transmission test for linkage disequilibrium: the insulin gene region and insulin-dependent diabetes mellitus (IDDM) [J]. American journal of human genetics, 1993, 52 (3): 506.[4] Risch N, Merikangas K. The future of genetic studies of complex human diseases. Science, 1996,273:1516-1517[5] Jiang C, Zeng Z-B, Multiple trait analysis of genetic mapping for quantitative trait loci [J] .Genetics,1995, 140 (3):1111-1127[6]Yandell BS,Mehta T,Banerjee S,Shriner D,Venkataraman R,Moon JY,Neely WW,Wu H,Von Smith R,Yi NJ,R/qtlbim:QTL with Bayesian interval mapping in experimental crosses. Bioinformatics,2007,23(5):641-3[7]肖静，胡治球，汤在祥，隋炯明，李欣，徐辰武.多个相关数量性状主基因的联合分析方法[J].中国农业科学，2005，38(9):171724[8] 茆诗松 程依明 濮晓龙.概率论与数理统计教程[M].高等教育出版社.2011.2[9]李红旺. 基于LASSO算法的多性状QTL条件定位方法[D].上海上海交通大学，201423.[10] Donoho D L, Johnstone I M, Ideal spatial adaptation by wavelet shrinkage[J]. Biometrika, 1994, 81: 425-455. [11] Friedman J,Hastie T, Hoefling H,et al, Pathwise coordinate optimization[J]. The Annals of Applied Statistics,2007 2(1): 302-332[12]胡文明. 作物复杂性状QTL定位相关的几个问题的探讨[D].扬州扬州大学，201462.[13]Team RDC.R:A language and environment for statistical computing.ISBN 3-900051-07-0.R Foundation for Statistical Computing.Vienna,Austria,2013.url://www.R-project.org.2005
目录
摘要3
关键词3
Abstract3
Key words3
1、绪论3
1.1研究意义 3
1.2多性状联合统计方法的关联分析简介3
1.3国内外研究背景 4
2.2理论及方法 4
2.1似然比检验4
2.2多性状定位 5
2.3 LASSO算法 6
2.4 多元最小二乘法 7
2.5 偏最小二乘 7
2.6 R语言简介 8
3、仿真实验及结果分析 8
3.1 实验结果 8
3.2结果分析9 3.3讨论 10
致谢10
参考文献10
基于多性状联合统计方法的关联分析
引言
版权保护: 本文由 hbsrm.com编辑，转载请保留链接: www.hbsrm.com/jsj/jsjkxyjs/1705.html
