小数据教程-数据分析

使用GBLUP和交叉验证进行预测
使用GFBLUP和交叉验证进行预测
基因集富集分析

一、使用GBLUP和交叉验证进行预测

这里，我们将展示如何使用qgg中的greml函数执行基因组**线性无偏预测(GBLUP)分析和交叉验证。这包括在训练集中使用限制最大似然估计(REML)估计方差成分，在验证集中使用GBLUP进行预测。这将在果蝇遗传参考面板(DGRP)中可用的“抗饥饿”表型和基因型数据中得到说明。
要执行GBLUP分析，以下输入数据是必不可少的。
y: vector of phenotype
X: design matrix for covariables
W: centered and scaled genotype matrix
G: genomic relationship matrix

该脚本包括以下步骤:
1)加载和准备数据
2)用于估计方差成分的限制最大似然(REML)分析
3)使用GBLUP模型的REML分析和交叉验证。

#library(devtools)#install_github("psoerensen/qgg")library(qgg)

加载并准备GBLUP分析数据

加载表型和共变量数据

load(file = "./phenotypes/starv_inv_wo.Rdata")dim(starv)## [1] 406 20head(starv)## L sex y In2Lt In2RNS In2RY1 In2RY2 In2RY3 In2RY4 In2RY5 In2RY6## 1 100 M 49.28000 INV/ST ST ST ST ST ST ST ST## 2 101 M 47.20000 INV/ST ST ST ST ST ST ST ST## 3 105 M 51.04000 ST ST ST ST ST ST ST ST## 4 109 M 44.96000 INV/ST ST ST ST ST ST ST ST## 5 129 M 33.08475 ST ST ST ST ST ST ST ST## 6 136 M 63.04000 ST ST ST ST ST ST ST ST## In2RY7 In3LP In3LM In3LY In3RP In3RK In3RMo In3RC wo## 1 ST ST ST ST ST INV ST ST y## 2 ST ST ST ST ST ST ST ST n## 3 ST ST ST ST ST INV ST ST n## 4 ST ST ST ST ST ST ST ST n## 5 ST ST ST ST ST ST ST ST n## 6 ST INV/ST ST ST ST INV/ST ST ST y

创建一个抗饥饿表型的向量y。

data <- starvy <- data$y

为共变量准备设计矩阵X。fm是用于在模型中包含相关变量以构建设计矩阵的公式。

fm <- y ~ wo +In2Lt + In2RNS + In3RP + In3RK + In3RMoX <- model.matrix(fm, data=data)dim(X)## [1] 406 12X[1:5,]## (Intercept) woy In2LtINV/ST In2LtST In2RNSINV/ST In2RNSST In3RPINV/ST## 1 1 1 1 0 0 1 0## 2 1 0 1 0 0 1 0## 3 1 0 0 1 0 1 0## 4 1 0 1 0 0 1 0## 5 1 0 0 1 0 1 0## In3RPST In3RKINV/ST In3RKST In3RMoINV/ST In3RMoST## 1 1 0 0 0 1## 2 1 0 1 0 1## 3 1 0 0 0 1## 4 1 0 1 0 1## 5 1 0 1 0 1

加载中心和缩放基因型矩阵W

load(file= "./genotypes/dgrp2_W2.Rdata")dim(W)## [1] 205 1725755W[1:5,1:5]## 2L_5317 2L_5372 2L_5390 2L_5403 2L_5465## 21 -0.2486289 -0.6646914 1.2122772 -0.4059216 -0.4007831## 26 -0.2486289 -0.6646914 -0.8200699 -0.4059216 -0.4007831## 28 -0.2486289 -0.6646914 1.2122772 -0.4059216 -0.4007831## 31 4.0006648 -0.6646914 -0.8200699 -0.4059216 -0.4007831## 32 -0.2486289 -0.6646914 -0.8200699 -0.4059216 -0.4007831

加性基因组关系矩阵G(VanRaden PM. 2008. J Dairy Sci. 91:4414-4423) 使用所有遗传标记构建如下:
G=WW′/m 其中W为中心和缩放的基因型矩阵，m为标记总数。
从基因型矩阵W计算累加的基因组关系矩阵G。

L <- data$LG <- grm(W=W)G <- G[L,L] # 这确保了G中的行id对应于数据中的行id，每一行都有雄性和雌性数据

限制最大似然(REML)分析

greml函数的幕后工作
REML分析用于估计方差分量，σ2e σ2e是线性混合模型中的随机效应:
y=Xb+Zg+e
其中y是表型观察的向量，X和Z为固定和随机效应的设计矩阵，b是固定效应向量，g是所有遗传标记捕获的基因组值的载体，e是残差向量。随机基因组值和残差假设为独立的正态分布值，如下所述：g∼N(0,Gσ2g) and e∼N(0,Iσ2e). 因此，我们假设观察到的表现型 y∼N(Xb,V) ；V=ZGZ′σ2g+Iσ2e .

在此，我们根据在整个研究人群中观察到的表型，使用GBLUP预测遗传价值:
g=σ2g GZ′V−1(y−Xb)

表型预测为:
y=Xb^+Zg^
greml函数在收敛之前要经过多次迭代(即，连续轮之间的参数变化小于指定的阈值，请参阅greml帮助页面中的“tol”参数)。在本例中，每次迭代都返回方差组件(第3列)和(第4列)的值。
greml函数返回一个列表结构，其中包括对固定效果(b)、随机效果(g)和剩余效果(e)的估计。列表中的其他值在greml帮助页面中描述。将verbose = FALSE更改为verbose = TRUE，可以查看迭代。

fitG <- greml(y=y, X=X, GRM=list(G=G), verbose = FALSE)

REML使用GBLUP模型进行分析和交叉验证

根据训练和验证群体之间的基因组关系来预测遗传价值:

基因组的关系矩阵:

根据训练(t)数据Gtt中的个体之间、验证(v)数据Gvv中的个体之间以及训练和验证数据Gvt中的个体之间的关系进行划分。因此，在训练数据中使用估计的方差成分(檩子2g和檩子2e)来预测总基因组易感性。最右边的术语(yt−Xtbt)构成了训练数据中个体为固定效应而校正的表型。反项V^ - 1t本质上是校正后的表型的方差-协方差结构。这两个术语相乘就是训练数据中个体的标准化和校正的表型，这些表型被投射到训练数据和验证数据之间的总遗传协方差结构上。
估计训练集中观察到的表型的方差成分。根据估计的方差分量，预测验证集的表型为:

通过从1 - 406 (y的长度)中随机取样30个值并重复此取样50次，可以创建50个验证集。验证集保存在验证矩阵中。这个矩阵指定GREML分析中使用的数据行。

n <- length(y)validate <- replicate(50, sample(1:n, 30))cvG <- greml(y = y, X = X, GRM = list(G=G), validate = validate)str(cvG)## List of 4## $ accuracy:'data.frame': 50 obs. of 7 variables:## ..$ Corr : num [1:50] 0.087 -0.049 -0.058 0.388 -0.026 0.264 0.277 -0.004 0.161 0.202 ...## ..$ R2 : num [1:50] 0.008 0.002 0.003 0.15 0.001 0.07 0.077 0 0.026 0.041 ...## ..$ Nagel R2 : num [1:50] NA NA NA NA NA NA NA NA NA NA ...## ..$ AUC : num [1:50] NA NA NA NA NA NA NA NA NA NA ...## ..$ intercept: num [1:50] 50.8 52.7 54.5 49.1 53.7 ...## ..$ slope : num [1:50] 0.017 -0.012 -0.015 0.087 -0.006 0.05 0.076 -0.001 0.056 0.042 ...## ..$ MSPE : num [1:50] 319 153 203 194 277 ...## $ theta :'data.frame': 50 obs. of 2 variables:## ..$ G: num [1:50] 27.3 29.2 30.8 24.6 29.5 ...## ..$ E: num [1:50] 138 150 145 150 140 ...## $ yobs : num [1:1500] 26.6 37.9 40.6 50.1 47.8 ...## $ ypred : num [1:1500] 53.2 55.8 44.9 44.9 58.7 ...library(knitr)library(kableExtra)kable(cvG$accuracy[1:5,], caption = "Cross Validation Predictive Ability") %>%kable_styling(full_width = F, position = "left")

交叉验证预测能力
CorrR2Nagel R2AUCinterceptslopeMSPE
0.0870.008NANA50.8400.017319.024
-0.0490.002NANA52.728-0.012152.879
-0.0580.003NANA54.540-0.015202.904
0.3880.150NANA49.0770.087194.324
-0.0260.001NANA53.673-0.006276.643

kable(cvG$theta[1:5,], caption = "Cross Validation Parameter Estimates") %>%kable_styling(full_width = F, position = "left")

交叉验证参数估计
G-------- E
27.278138.422
29.239150.142
30.825144.888
24.614150.377
29.532140.285

GREML交叉验证分析的输出
输出包括量化模型预测能力的统计数据，通过回归验证数据集观察到的表型与预测的表型进行评估:
y=intercept+y^slope+e
cvG$accuracy中值的解释:

ValueDescription

Corr	预测表型值与观察表型值的相关性。平均列出的50个相关关系产生了预测能力
R2	R2, GBLUP模型所能解释的总方差的比例
–	–
Nagel R2	Nagelkerke’s R
AUC	ROC曲线下面积
–	–
intercept	y回归到y^的y轴截距
slope	y回归到y^的斜率
–	–
MSPE	均数平方预测误差= 1/nv∑nvi=1 (yi−yi^)2, nv =验证集观测数

cvG$theta值的解释:

准备结果的数据框架。

accuMeans <- colMeans(cvG$accuracy)thetaMeans <- colMeans(cvG$theta)accuSEM <- apply(cvG$accuracy, 2, function(x){sd(x)/sqrt(50)})thetaSEM <- apply(cvG$theta, 2, function(x){sd(x)/sqrt(50)})results1 <- round( rbind(accuMeans, accuSEM), digits=3) results2 <- round( rbind(thetaMeans, thetaSEM), digits=3) results <- cbind(c("mean","sem"),results1,results2)rownames(results) <- NULLkable(results, caption = "Results Summary")%>%kable_styling(position = "left")

结果摘要

遗传率(h 2)

h2fm <- cvG$theta$G/(cvG$theta$G+cvG$theta$E)boxplot(h2fm, main = "Genomic Heritability")kable(cvG$theta[1:5,], caption = "Cross Validation Parameter Estimates") %>%kable_styling(full_width = F, position = "left")

交叉验证参数估计

GREML交叉验证分析的输出

输出包括量化模型的预测能力的统计数据，通过回归观察到的表型对验证数据集的预测表型进行评估
y=intercept+y^slope+e

准备结果的数据框架。

accuMeans <- colMeans(cvG$accuracy)thetaMeans <- colMeans(cvG$theta)accuSEM <- apply(cvG$accuracy, 2, function(x){sd(x)/sqrt(50)})thetaSEM <- apply(cvG$theta, 2, function(x){sd(x)/sqrt(50)})results1 <- round( rbind(accuMeans, accuSEM), digits=3) results2 <- round( rbind(thetaMeans, thetaSEM), digits=3) results <- cbind(c("mean","sem"),results1,results2)rownames(results) <- NULLkable(results, caption = "Results Summary")%>%kable_styling(position = "left")

结果摘要

基因组遗传力

theta <- cvG$thetah2f2 <- theta["2R"]/sum(theta)h2f4 <- theta["3R"]/sum(theta)layout(matrix(1:2,ncol=2))boxplot(h2f2, main = "Genomic Heritability 2R", cex.main=0.8)boxplot(h2f4, main = "Genomic Heritability 3R", cex.main=0.8)

计算得出。单标记效果计算如下:

计算出单标记效应的(co)方差为:

个体标记与表型的联系是基于单一标记测试统计，如t统计和该统计的阈值。

其中Var(s^j)是由单标记效应的(co)方差矩阵对角的第j个元素得到的s^的第j个元素的方差估计。
零假设下，s^j=0,假设ts^j服从dfe剩余自由度的分布.(for further information, see Sørensen IF. et al., 2017. Scientific Reports. doi:10.1038/s41598-017-02281-3).

基因集合富集分析

在本例中，与GO项相关的遗传标记被认为是一个标记集(基因组特征)。因此，在这种情况下，基因组特征是由GO项目定义的。这里我们描述了四种不同的集合测试方法，它们都来自于GBLUP模型。
这包括
1)协方差协会组测试(CVAT),
2)基于评分的方法,
3)测试统计基于遗传标记的总和属于同一基因的功能,和
4)基于计算检验统计量的数量遗传标记的功能与特征相关的表型(hyperG)。

在这个例子中，我们考虑了一个竞争性零假设，也就是说，相关的标记是从被测试遗传标记的总数中随机挑选出来的。

协方差关联检验统计量(Tcvat)

在这里，我们考虑所有标记的总基因组效应和基因组特征中遗传标记集的基因组效应之间的协方差


在这个表达式,g^r=∑mr i=1 wis^i ,是剩余标记的基因组效应，例如，不包括在特性中的标记。mf和mr分别为feature中标记的数量和剩余的标记集。TCVAT也与解释的遗传标记的平方和有关,该检验统计量在原假设下的分布很难用精确或近似的分布来描述，需要经验分布。通过对W中的mf列随机抽样，可以得到竞争零假设的经验分布。
gsea函数中的fit参数是使用greml函数对线性混合模型进行拟合得到的fit对象。
gsea函数返回一个数据框架(mma，多标记关联对象)，其中包括协方差测试统计量(setT)、集合中的标记数(g)和p值§。

mma <- gsea(fit = fit, sets = goSets, W = W, method = "cvat")

基于评分的方法统计(Tscore)

它是由似然的一阶导数推导出来的，就像Rao的分数测试一样(Rao et al. 1948)。比饶的分数测试关键的区别是,只有一阶导数的二次项的基础测试统计(Goeman et al. 2004, Wu et al. 2011, Wang et al. 2013 )从一个论点,这是唯一的一部分,涉及到数据(Huang and Lin, 2013)。因此，这里给出的基于分数的方法相当于序列核关联测试(SKAT, Wu et al. 2011)。因此得分统计量可以写成:

其中固定效应b和表型协方差矩阵V在零模型下估计。空模型的目的是为了调整环境的非遗传因素，遗传因素不是基因组特征的一部分，包括群体结构。对于GBLUP零模型，基因组效应可以定义为g ~ N(0, g ggg2g)，也可以定义为g ~ N(0,Gr Gr Gr Gr 2r) (r = remaining)。在第一种情况下，基因组关系矩阵是使用所有遗传标记计算，因此null模型只需要拟合一次。在后一种情况下，只使用不包括在基因组特征中的遗传标记进行计算，这需要为每个基因组特征拟合不同的模型。
得分法的测试统计量集合可以改写为:


通过对W中的mf列随机抽样，得到了在竞争零假设下得分检验统计量的经验分布。
gsea函数返回一个数据框架(mma)，其中包括协方差测试统计量(stat)，集合中标记数(g)和p值§。

mma <- gsea(fit = fit, sets = goSets, W = W, method = "score")

基于总和(Tsum)设置检验统计量
这个总和集测试是基于位于特征集中的所有标记汇总统计量的测试统计量之和，例如:
其中ti表示第i个单变量检验统计量，如标记效应(s^)或t统计量。可以从线性模型分析(从PLINK或使用qgg lma近似)中获得单一标记汇总统计信息，或者greml函数中的单个或多个组件REML分析(GBLUP或GFBLUP)。如果基因组特征含有许多具有小到中度影响的遗传标记，那么总和检验是强有力的。该检验统计量在竞争零假设下的分布是未知的，需要一个经验分布。s^和t统计量都用来计算Tsum。

基于sum的set测试可以由gsea函数执行，并将方法参数指定为“sum”。

ma <- lma(y=y,X=X,W=W[L,])mma <- gsea(stat = ma[,"stat"]**2, sets = goSets, method = "sum")head(mma)## m stat p## GO:0000022 315 424.3648 0.253## GO:0000027 644 921.7163 0.145## GO:0000028 283 307.9064 0.466## GO:0000045 1602 1659.4118 0.583## GO:0000060 1816 1587.1936 0.881## GO:0000070 3416 3458.3678 0.703

基于hyperG计算测试统计量(Tcount)

该检验统计量是通过对特征中与性状表型相关的遗传标记的数量进行计数，计算方法为:

其中mf为特征中标记数，ti为第i个单标记检验统计量(如t统计量)，t0为单标记检验统计量的任意选择阈值,I是和指示器函数，如果满足参数(abs(ti)>t0)，它将取值1。在零假设(即个体相关的标记随机分布)下，假设Tcount ~ Hyper(m,ma,mf)是一个参数为m的超几何分布(测试的遗传标记总数)的实现，ma(所有标记中相关遗传标记的总数)和mf(特征中遗传标记的总数)。gsea函数可以执行基于sum的set测试，并将方法参数指定为“hyperg”。

ma <- lma(y=y,X=X,W=W[L,])head(ma)## s se stat p## 2L_5317 0.4483611 0.6761037 0.6631544 0.5076099## 2L_5372 0.2898735 0.7143767 0.4057712 0.6851255## 2L_5390 0.6692177 0.7099992 0.9425612 0.3464687## 2L_5403 0.4075746 0.6727596 0.6058251 0.5449711## 2L_5465 0.4404401 0.6781146 0.6495069 0.5163799## 2L_5598 -0.2743689 0.6795625 -0.4037434 0.6866151mma <- gsea(stat = ma[,"p"], sets = goSets[1:10], method = "hyperg", threshold = 0.05)head(mma)## GO:0000022 GO:0000027 GO:0000028 GO:0000045 GO:0000060 ## 5.036082e-02 3.891661e-06 9.190028e-01 9.513077e-01 9.999934e-01 ## GO:0000070 ## 9.990277e-01

结果表明，每组snp与表型之间的显著性相关的p值。