（四）逆正态转换 / 开普饭

逆正态转换(inverse normal transformation，INT)

基本介绍

考虑连续性状的遗传关联分析。如果残差分布呈偏态或峰度分布(容易出现异常值)，标准线性回归可能无法控制中等样本量的I类错误。当残差分布为非正态分布时，即使样本足够大，标准线性回归也不能提供有效的推理。可能呈现非正态分布的特征包括体质指数、循环代谢物、基因表达、多导睡眠图信号和肺活量测量。在这种情况下，基于秩的逆正态变换(inverse normal transformation，INT)被用来抵消偏离的正态性。

逆正态转换的两种方式

1. qx <- qnorm((rank(x)-0.5)/sum(x))

参考文献：Yang, Jian, et al. "FTO genotype is associated with phenotypic variability of body mass index." Nature 490.7419 (2012): 267.

# ifyou have missing data

qx<-qnorm((rank(x,na.last="keep")-0.5)/sum(!is.na(x))

2. RNOmni Package

RNOmni Package进行INT时，首先将样本测量值映射到概率尺度，用秩分数替换观测值，然后使用probit函数将其转换为z-score。在以下示例中，一个大小为n=1000的样本来自卡方分布。经过变换，测量值的经验分布与标准正态分布没有明显区别。

if(!require(RNOmni)) install.packages("RNOmni")

if(!require(ggpubr)) install.packages("ggpubr")

# Sample from the chi-1 square distribution

y = rchisq(n=1000,df=1);

# Rank-normalize

z = rankNorm(y);

A=data.frame(type=rep(c("original","RankNorm"),each=1000),feature=c(y,z))

ggdensity(A, x="feature", add = "mean", rug = TRUE, color = "type", fill = "type", palette = c("#00AFBB", "#E7B800"))

(1) 模拟数据

模拟数据产生1000个样本，在连锁平衡中绘制1000个基因座的基因型，其中MAF=0.25。

在模型矩阵X中，包含一个截距项、4个标准正态协变量Z和前4个遗传主成分。截距项设为1，其余回归系数模拟为随机效应。用协变量解释表型变异的比例为20%，主成分解释的变异占5%。模拟了两种具有加性残基的表型：① yn，具有标准正态残差；② yt，有t3残差，按比例得到单位方差。

## Simulated Data

set.seed(100);

# Sample size

n = 1e3;

## Simulate genotypes

G = replicate(rbinom(n,size=2,prob=0.25),n=1e3);

storage.mode(G) = "numeric";

# Genetic principal components

S = svd(scale(G))$u[,1:4]; #矩阵奇异值分解svd

S = scale(S);

# Covariates

Z = scale(matrix(rnorm(n*4),nrow=n));

# Overall design

X = cbind(1,Z,S);

# Coefficient

b = c(1,rnorm(n=4,sd=1/sqrt(15)),rnorm(n=4,sd=1/sqrt(60)));

# Linear predictor

h = as.numeric(X%*%b);

# Normal phenotype

yn = h+rnorm(n);

# T-3 phenotype

yt = h+rt(n,df=3)/sqrt(3);

A=data.frame(type=rep(c("Normal Pheno","T3 pheno"),each=1000),pheno=c(yn,yt))

ggdensity(A, x="pheno", add = "mean", rug = TRUE, color = "type",palette = c("#00AFBB", "#E7B800"))

(2) 数据格式

结局y：数值向量

基因型G：数值矩阵（行为样本，列为基因型）

模型矩阵X：数值矩阵（调整协变量和群体结构，包含一个截距项）

提前扩展factors和交互：model.matrix()

无论在结局向量y还是模型矩阵X中，都不期望出现缺失值，但是基因型矩阵G可能包含缺失值。

只有在基因型缺失的位点，才会将缺失基因型的观察值从关联检验中排除。

(3) 基本的关联检验

基本的关联检验是(未转换)表型对基因型和协变量的线性回归。BAT()提供了一个使用表现型y、基因型G和模型矩阵X的高效实现。标准输出包括分数统计量（score statistic）、它的标准误、z-score和 p值（对G的每一列都可以得到一行相应的统计量）。设置test="Wald"指定一个Wald检验。Wald检验可能提供更多的power，但通常比较慢。设置simple=T只返回p值。

## Basic Association Test

# Basic Association Test, Normal Phenotype

Results_BAT_Yn = BAT(y=yn,G=G,X=X);

round(head(Results_BAT_Yn),digits=3);

# Basic Association Test, T3 Phenotype

Results_BAT_Yt = BAT(y=yt,G=G,X=X);

round(head(Results_BAT_Yt),digits=3);

(4) 逆正态转换

假设对n个样本中的每一个都观察到一组连续的测量ui，当ui按升序排列时，让rank(ui)表示每个样本的测量ui的秩。基于秩的逆正态转换（INT）定义为：

其中，φ-1,是probit函数；k∈（0,1/2）为可调整的偏移量（adjustable offset）。默认采用，k=3/8的Blom offset。

1) 直接逆正态转换（Direct INT，D-INT）

在D-INT中，INT转化的表型在基因型和协变量上回归。当表型可能是由一个潜在的正态性状的秩保留转换（rank-preserving transformation）而产生时，D-INT是最强大（powerful）的。DINT()使用表型y、基因型G和模型矩阵X实现关联检验。标准输出包括检验统计量、标准误、z-score和 p值（对G的每一列都可以得到一行相应的统计量）。可用Wald和score检验。设置simple=T只返回p值。

###Direct INT

# Direct INT Test, Normal Phenotype

Results_DINT_Yn = DINT(y=yn,G=G,X=X);

round(head(Results_DINT_Yn),digits=3);

# Direct INT Test, T3 Phenotype

Results_DINT_Yt = DINT(y=yt,G=G,X=X);

round(head(Results_DINT_Yt),digits=3);

2) 间接逆正态转换（indirect INT，I-INT）

在I-INT中，首先，表型和基因型在协变量上回归得到残差。表型残差按秩归一化处理。然后，在基因型残差上对逆正态转化的表型残基进行回归。当表型在协变量上是线性时，I-INT最强大（powerful），但是残差分布呈偏态或峰度分布。IINT()使用表型y、基因型G和模型矩阵X实现关联检验。标准输出包括检验统计量、标准误、z-score和p值。设置simple=T只返回p值。

### Indirect INT

# Indirect INT Test, Normal Phenotype

Results_IINT_Yn = IINT(y=yn,G=G,X=X);

round(head(Results_IINT_Yn),digits=3);

# Indirect INT Test, T3 Phenotype

Results_IINT_Yt = IINT(y=yt,G=G,X=X);

round(head(Results_IINT_Yt),digits=3);

3) 综合逆正态转换（Omnibus INT，O-INT）

由于D-INT和I-INT都不是最强大的（powerful），因此O-INT自适应地将它们组合成一个健壮的、统计上强大的检验。在方法内部，OINT()同时应用D-INT和I-INT，然后选择为关联提供更多证据的检验。具体地说，综合统计量（omnibus statistic）是D-INT和I-INT中p值的最小值。为该统计量分配p值需要估计D-INT和I-INT检验之间的相关性。当考虑多个基因座（loci）时，一种计算效率高的策略是通过对相邻基因座进行平均来估计相关性。如果考虑的位点较少，或者需要对特定位点进行估计，则可以通过bootstrap估计相关性。最后，用户可以手动提供他们自己的相关性估计。

OINT使用表型y、基因型G和模型矩阵X实现了综合检验。标准输出包括由D-INT、I-INT和O-INT分别估计的p值。设置keep.stats=T保留D-INT和I-INT 的Z-statistics。设置keep.rho=T保留D-INT和I-INT之间的估计相关性。

## Omnibus INT

cat("Omnibus Test, Normal Phenotype, Average Correlation Method\n");

Results_OINT_Avg_Yn = OINT(y=yn,G=G,X=X,method="AvgCorr");

round(head(Results_OINT_Avg_Yn),digits=3);

cat("\n");

cat("Omnibus Test, Normal Phenotype, Bootstrap Correlation Method\n");

Results_OINT_Boot_Yn = OINT(y=yn,G=G[,1:10],X=X,method="Bootstrap",B=100);

round(head(Results_OINT_Boot_Yn),digits=3);

cat("\n");

cat("Omnibus Test, T3 Phenotype, Average Correlation Method\n");

Results_OINT_Avg_Yt = OINT(y=yt,G=G,X=X,method="AvgCorr");

round(head(Results_OINT_Avg_Yt),digits=3);

cat("\n");

cat("Omnibus Test, T3 Phenotype, Bootstrap Correlation Method\n");

Results_OINT_Boot_Yt = OINT(y=yt,G=G[,1:10],X=X,method="Bootstrap",keep.rho=T,B=100);

round(head(Results_OINT_Boot_Yt),digits=3);

cat("\n");

cat("Replicate the Omnibus Test on the T3 Phenotype, Manually Specifying Correlation\n");

Results_OINT_Boot_Yt = OINT(y=yt,G=G,X=X,method="Manual",set.rho=Results_OINT_Boot_Yt[,"Corr"],keep.rho=T);

round(head(Results_OINT_Boot_Yt),digits=3);

cat("\n ");

注意：并行计算

所有关联检验都可以并行运行。为此，请注册一个并行后端，

## Parallelization

install.packages("doMC", repos="http://R-Forge.R-project.org")

# register a parallel backend, e.g.

doMC::registerDoMC(cores=4) #then specify the parallel=T option.

?registerDoMC

注：考虑NA的处理方式，是否会对分布造成影响

参考阅读：

https://www.biostars.org/p/312945/

https://cran.r-project.org/web/packages/RNOmni/vignettes/RNOmni.html

（四）逆正态转换

相关推荐