给一堆数据后怎么用r处理成正态分布_全基因组关联分析项目设计——基因定位研究中的表型处理...

在上期文章《全基因组关联分析项目设计——基因分型策略》中，周老师介绍了四种基因分型技术的特点及应用，并对这四种技术做了简要比较。本期文章将为大家介绍基因定位研究中表型处理的相关知识。 表型-基因型关联分析，是寻找与性状相关基因的重要手段。在基因型检测的手段不断丰富(重测序、GBS、多重PCR、SNP芯片等)以及成本不断降低的时代背景下，表型检测和数据前处理，就显得尤为重要，因为这将直接影响关联(连锁)分析结果的准确性。

一、表型的类型

常见的表型性状，我们可以将其分为三种： 数量性状，质量性状与分类性状 (表1)。 数量性状在遗传育种研究中十分常见。此类性状由多基因控制，且可以用数字量化，例如产量、株高等，所以它们也比较容易量化且适用于大部分线性回归分析模型。质量性状是一种简单的离散型分类性状。严格意义上讲，单基因组控制的性状才可能被定义为质量性状，例如人类单基因家族遗传病，动植物突变体研究中的单基因突变体。 表1 三种类型性状的特点

以上的两种情况都是我们期望的理想情况。但自然界生物的表型比以上两种情况远为复杂，很多表型既不是质量性状又难以简单地用数字量化。例如，某种植物不同品种的花颜色可能包括蓝、黄、红、紫等多种颜色，叶片形状可能三角形、圆形、椭圆形、长条形等。山羊角的数量，可能是无角、正常(二角)、三角、四角等。 对于这种情况，则需要对表型本质进行剖析，然后对性状进行分级，我们将这种 通过人为观察从而可以进行分类的离散型变量统称为分级性状。 但要注意一点， 分级性状最终的定义，还是部分依赖于我们的经验 。例如，对于植物的抗病性，我们既可以按照叶片病斑的面积将其定义为连续型的数值型性状(0~100%)，也可以人为设定阈值将其定义为分级性状(高、中、低)。对人类的血压，既可以按照血压的高低将其定义为连续型的数值型性状，也可以人为设定一个阈值，然后将人群分为高血压组(病例)和健康组(对照)。 而病例-对照研究，正是人类复杂疾病研究的常见模式。 因此在实际应用中，符合简单离散分布的性状(例如分为两类)，可以是单基因控制的性状，也可以是复杂的多基因组控制性状，只是我们通过一定的标准将样本分为了两类。

二、表型值的处理

1. 分布类型的检验

对于符合一定遗传模式的性状，其性状分布模式也应该符合一定的特性。例如单基因控制的隐性性状，理论上应符合3:1的分离比，我们可以使用卡方检验来判断。对于多基因控制的数量性状，理论上其表型应该符合正态分布(又称作高斯分布)。 当我们拿到一组性状的时候，如何判断其是否符合正态分布呢？

图1 基因性状的正态分布 最简单的方法，可以R语言中自带的shapiro.test命令进行检验。如果P value > 5%，则说明数据分布近似正态分布；另外，也可以通过R语言hist命令对表型数据进行可视化(频率直方图)，从更直观的角度观察其是否符合正态分布的特点。

图2 利用R语言绘制数据的分布模式

关于表型数据的正态性判断，也可以登录我们的Omicshare论坛了解一下：

图3 论坛经验文章《数据正态性检验的方法》 文章链接： http://www.omicshare.com/forum/forum.php?mod=viewthread&tid=790&fromuid=12 正态性检验的意义在于可以用作辅助判断一个性状是否为多基因控制的数量性状，还有另一个意义在于，大部分表型-基因型关联分析的模型属于线性模型，线性模型要求数据符合正态分布，从而保证结果可靠性。 在实际项目中，在样本量较大的情况下，我们一般直接忽略表型数据的正态性检验。因为从统计检验的原理上说，样本数越多，统计检验的自由度越高。这种情况下，只要数值分布轻微偏离正态分布，就会得出很显著的P值(就是统计结果判定数值不符合正态分布)，这很显然是不符合现实情况的。所以， 一般情况下，可以直接通过观察频率直方图，判断性状值大体符合正态分布即可 。 相比性状是否符合正态分布，另一种情况更值得我们注意：那就是离群异常值样本。

2. 表型异常值(极端值)的处理

如果数量性状符合正态分布，理论上应该符合中间个体多，两端极端个体逐步减少的特点。 但有时候，数据中可能依然会出现离开群体很远的位置，突然孤零零冒出几个数值异常的个体。这种情况就非常值得我们注意了，这样的样本有大概率是表型检测的错误。例如，在mGWAS研究中(以代谢物为表型的关联分析)，在进行代谢物液相色谱检测的时候，软件有可能在部分样本中会将相邻峰与目标峰混淆，导致部分样本目标代谢物的表达量被高估。而这种情况，几乎只能靠人工去检查峰图文件才能发现问题。所以， 对于这种离群的异常值，我们都有必要回头去检查原始数据，判断其是否准确。

图4 群体整体分布中的极端值个体 如果无法判断真伪，则要考虑 将这些极端值个体去除 。因为极端异常值的存在可能引起关联分析结果的异常。 例如下图，就来源我们一个项目的关联分析结果，曼哈顿图中有大量关联的峰且这些峰的显著性相似。检查后，发现是由于群体中存在数个表型值“离奇”大的极端个体。由于这些个体表型值异常高，导致这几个个体特有的基因型无一例外地被软件判定为与性状强相关(P值显著，且显著性相似)。很显然，这样的结果并不符合数量性状的特点。因此，我们将这几个个体剔除重新分析，结果果然就恢复正常了。

图5 极端值导致的关联分析假阳性 那么，异常值如何才能被找出并去除呢？异常值去除的主要方法包括： (1) 排序观察法 ：即排序观察后手动去除我们认为异常的观测值； (2) 3 sigma规则 ：即在均值加减3倍标准差范围内的值为正常值，其他值为异常值。 如下图，在标准正态分布下，理论上落在3倍标准差外的样本数应该低于群体总样本数的0.3%。考虑到中小型GWAS研究的样本一般只有几百例，3倍标准差的样本如果过多，则显然是异常情况。

图6 标准正态分布的概率分布 (3) 箱线图方法 ：绘制箱线图，在触须以外的值均可以判断为异常值。 在实际关联分析中，如果表型种类较少，我们一般使用第一种方法，按照人工经验判断去除异常值即可。但如果表型种类过多(例如，检查了几十种数量性状)，人眼一个一个看自然效率就太低了，那么则可以考虑用第二种方法，自动将异常样本去除。

3. 数据的标准化

数据的标准化是指将数据按照比例缩放,使之落入一个特定的区间。例如，将数据统一映射到[0,1]区间上。在同时对各类表型开展关联分析的时候，不同类表型的数值以及变异范围可能非常巨大，甚至差了几个数量级。例如，A表型变化范围可能是0.2~0.9，B表型的变化范围可能为20,000~80,000。那么，如此巨大的表型值差异就不利于在完成分析后进行不同表型的比较(例如，表型在某个关联SNP位点的遗传方差)。 为了解决这个问题，则可以对数据进行标准化。数据标准化的基本原则就是在不改变一组数值相对大小的情况下(自然也不影响关联分析的结果)，对数据的整体进行调整。 常用的数据标准化方法有两种，一种是 z-score标准化 ，另一种是 min-max标准化 ： (1) z-score标准化 即我们常说的数据中心化方法之一，公式为：z = (x - μ)/σ； z为标准化后的值，x为原始表型，μ为这组表型的平均值，σ为这组表型的标准差。 在z-score标准化后，这组表型将变成均值为0，标准差为1的一组数，低于平均值的表型变为负数，反之则为正数。如此处理后(如下图)，数据的整体分布模式以及个体之间相对大小并没有改变(大哥依然是大哥)，只是数据的变异范围被压缩到以0为中心的一个小区间。如此，将便于对不同类型的性状进行比较。

图7 Z score中心化前后的比较 类似的，在eGWAS分析中，是以基因表达量为表型。由于不同基因的表达量差异极大，也往往需要进行Z score矫正后再开始后续的分析。 还有另外一个地方，也会用到Z score标准化，那就是绘制热图，相关内容可以查看以下帖子。

图8 论坛经验文章《Omishare Tools 中“热图”工具使用教程》 文章链接： http://www.omicshare.com/forum/forum.php?mod=viewthread&tid=407&fromuid=12 (2)min-max标准化 min-max标准化是一种将数据转换到0~1之间的方法，公式为y = (x – min(x))/(max(x) – min(x))；y为标准化后的值，x为原始值。 如果矫正后，最小值将转化为0，最大值将转化为1，其他数值根据相对大小，被重新分布在0~1之间。这种矫正效果也类似Z score矫正，但唯一的不足是个别极端值对矫正后的总体分布影响比较大，所以对于基因表达量等这类数据效果可能不大好。

4. 利用哑变量定义难以数值化的性状

一些性状属于描述型的多分类性状，不好直接数字化。例如，群体花色有红、黄、蓝3种颜色。由于红、黄、蓝没有明显的线性梯度关系，那么就不能将它们简单赋值为1、2、3，而是需要将它们按照合理的逻辑重新进行归类。这就需要引入哑变量对它们进行归类。 哑变量，又称为虚拟变量，是一种人为定义的变量类型，通常取值0和1。如花色的例子，将可以类似下图，针对每种颜色进行分组并赋予0,1变量，那么原来的一组颜色变量就被拆分为3组重新定义的变量。

图9 颜色分类的哑变量 然后就可以对3种方式赋值后的结果分别进行关联分析，获得与不同类型颜色关联的结果。 赋值时需要注意1和0比例不要太悬殊，否则可能会降低检验的功效 。   当然，哑变量的本质是对性状的重新定义。还是花色的例子，如果我们对颜色的理解更深刻，还可以用其他更好的方式定义它们。例如，我们如果确定了不同花色形成的本质原因是某些代谢物的含量或组成，那么我们就可以用代谢物含量或组成来代表花的颜色。  

图10 利用代谢物进行花色的关联分析

5. 多年多点表型的处理

为获得可靠的关联结果，我们通常会对同一个性状观测多次，多次观测可能是相同年份不同地点，也可能是不同年份相同(或不同)地点。 对于此类数据，我们可以根据性状的遗传机制选择不同的处理方式。 如果性状遗传力高，受环境影响不大，我们可以根据多年多点的结果取均值或BLUP值作为该性状的代表值进行分析。 如果性状遗传力低，受环境影响大，我们可以每年每点单独分析后综合评判结果，在获得定位结果后(例如获得了10个关联位点)，那么可以利用多元回归模型开展基因-环境互作的分析。

图11 受环境因素影响大的表型可以多次结果进行独立的分析和比较 关于BLUP分析(全称最佳线性无偏估计)可能是很多小伙伴关心的。这个算法可以利用模型来估算年份、地点等效应对表型影响的效应值，从而我们可以计算剔除这些环境因素后剩下的基因型效应是多少，最后可以将多年多点的表型数据整合为一致的基因型效应数据(剔除环境效应后的表型值)来开展关联分析。 关于BLUP分析的R语言操作，也可以收看我们的在线操作视频：

课程链接： https://www.omicshare.com/class/home/index/series?id=41

视频观看方式

电脑端： 登录Omicshare课堂 www.omicshare.com/class 观看学习 手机端： 通过点击基迪奥微信公众号底部菜单栏【视频教程】观看学习

实用科研工具推荐      
详实生信软件教程分享
前沿创新组学文章解读
独家生信视频教程发布