小数据教程-数据分析

使用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  20
head(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 <- starv
y <- data$y

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

fm <- y ~  wo +In2Lt + In2RNS + In3RP + In3RK + In3RMo
X <- model.matrix(fm, data=data)
dim(X)
## [1] 406  12
X[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 1725755
W[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$L
G <- 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")

交叉验证预测能力
Corr R2 Nagel R2 AUC intercept slope MSPE
0.087 0.008 NA NA 50.840 0.017 319.024
-0.049 0.002 NA NA 52.728 -0.012 152.879
-0.058 0.003 NA NA 54.540 -0.015 202.904
0.388 0.150 NA NA 49.077 0.087 194.324
-0.026 0.001 NA NA 53.673 -0.006 276.643

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

交叉验证参数估计
G -------- E
27.278 138.422
29.239 150.142
30.825 144.888
24.614 150.377
29.532 140.285

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

Value Description
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")


遗传值直方图

hist(fitG$g, xlab = "genetic value", cex.lab=0.8, cex.axis=0.8, main = "REML", cex.main=0.8)


表型的柱状图

hist(y, xlab = "hours", cex.lab=0.8, cex.axis=0.8, main = "Starvation Resistance", cex.main=0.8)

二、使用GFBLUP和交叉验证进行预测

这里，我们将展示如何使用qgg中的greml函数执行基因组特征最佳线性无偏预测(GFBLUP)分析和交叉验证。GFBLUP模型包括一个额外的基因组效应，量化了一组标记的集体行为，即一个基因组特征。在这个例子中，与染色体相关的标记被认为是一个特征(标记集)。GFBLUP的执行涉及使用限制最大似然估计(REML)估计每个特征集的方差分量。GFBLUP模型预测表型的能力通过交叉验证程序进行评估。这将通过Drosophil提供的“抗饥饿”表型和基因型数据来说明
要执行GFBLUP分析，以下输入数据是必不可少的。
y : vector of phenotype
X: design matrix for covariables
Genomic feature marker sets
W: centered and scaled genotype matrix
G: genomic relationship matrix for each feature

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

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

加载并准备用于GFBLUP分析的数据

加载抗饥饿数据的已清理和编辑的数据框。

load(file = "./phenotypes/starv_inv_wo.Rdata")
dim(starv)
load(file = "./phenotypes/starv_inv_wo.Rdata")
dim(starv)

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

data <- starv
y <- data$y
length(y)
## [1] 406

为共变量准备设计矩阵X。fm是设计矩阵中包含相关协变量的公式。

fm <- y ~  wo + In2Lt + In2RNS + In3RP + In3RK + In3RMo
X <- model.matrix(fm, data=data)
X[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 1725755
W[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

加载染色体标记集。
该数据集中的标记分布在6条染色体臂上，即2L, 2R, 3L, 3R, 4和x。在本例中，每条染色体及其相关标记被认为是一个标记集(genomic feature, GF)。

load(file = "./annotation/chrSets2.Rdata")
str(chrSets)
## List of 6
##  $ 2L: chr [1:406577] "2L_5317" "2L_5372" "2L_5390" "2L_5403" ...
##  $ 2R: chr [1:327967] "2R_10037" "2R_10468" "2R_10959" "2R_12079" ...
##  $ 3L: chr [1:390711] "3L_39998" "3L_40145" "3L_40202" "3L_40635" ...
##  $ 3R: chr [1:368096] "3R_1339" "3R_1651" "3R_2158" "3R_2318" ...
##  $ 4 : chr [1:2686] "4_61790" "4_62622" "4_62905" "4_62908" ...
##  $ X : chr [1:229718] "X_19380" "X_19797" "X_20390" "X_20491" ...

只保留W矩阵中的snp。

setsGF <- lapply(chrSets,function(x){ x[x%in%colnames(W)] })

创建对象nsets:集合的数量(染色体的数量)和nsnps:每条染色体的SNPs数量。

nsets <- length(setsGF)
nsnps <- sapply(setsGF,length)

加性基因组关系矩阵G (VanRaden PM)。2008. 利用所有遗传标记G=WW ’ /m，其中W为中心和缩放的基因型矩阵，m为标记总数。基因组特征关系矩阵(GF) Gk=WkW 'k /mk，为第k个遗传标记集的加性基因组关系矩阵。

从基因型矩阵W计算累加的基因组关系矩阵G。

L <- data$L
GF <- lapply(setsGF, function(x) {grm(W = W[L, x])})

测定了每个系雄性和雌性的抗饥饿能力。因此GF每条染色体有406行和406列的维度。

str(GF)
## List of 6
##  $ 2L: num [1:406, 1:406] 0.50583 0.09563 -0.00765 0.11707 -0.03977 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:406] "100" "101" "105" "109" ...
##   .. ..$ : chr [1:406] "100" "101" "105" "109" ...
##  $ 2R: num [1:406, 1:406] 0.996919 -0.000345 -0.000574 -0.005291 -0.005337 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:406] "100" "101" "105" "109" ...
##   .. ..$ : chr [1:406] "100" "101" "105" "109" ...
##  $ 3L: num [1:406, 1:406] 1.0584 -0.0089 0.0159 0.0062 0.0072 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:406] "100" "101" "105" "109" ...
##   .. ..$ : chr [1:406] "100" "101" "105" "109" ...
##  $ 3R: num [1:406, 1:406] 1.4228 -0.0146 0.3536 -0.0101 -0.0248 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:406] "100" "101" "105" "109" ...
##   .. ..$ : chr [1:406] "100" "101" "105" "109" ...
##  $ 4 : num [1:406, 1:406] 0.599 -0.254 -0.319 0.239 0.108 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:406] "100" "101" "105" "109" ...
##   .. ..$ : chr [1:406] "100" "101" "105" "109" ...
##  $ X : num [1:406, 1:406] 1.01855 0.00487 0.02768 0.0022 -0.01247 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:406] "100" "101" "105" "109" ...
##   .. ..$ : chr [1:406] "100" "101" "105" "109" ...

限制最大似然(REML)分析

greml函数返回一个列表结构，其中包括对固定效应(b)、随机效应(yk)、剩余效应(e^)的估计，以及方差成分(xyc2gk)。列表中的其他值在greml帮助页面中描述。将verbose = FALSE更改为verbose = TRUE，可以查看迭代。略…

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

从训练数据(t,365行)中估计表型的基因组参数，从验证数据(v, 41行)中预测表型。根据训练和验证群体之间的基因组关系来预测遗传价值:

k^th的基因组关系矩阵

根据训练数据中的个体之间的关系Gttk，验证数据中的个体之间的关系Gvvk，以及训练数据和验证数据中的个体之间的关系Gvtk进行划分。因此，使用估计的方差组成部分来预测总基因组易感性(utils2g1，…在训练数据中，如:eurch2nf和eurch2e)。最右边的术语(yt−Xtbt)构成了训练数据中个体为固定效应而校正的表型。反项V - 1t本质上是校正后的表型的方差-协方差结构。这两个术语相乘就是训练数据中个体的标准化和校正的表型，这些表型被投射到训练数据和验证数据之间的总遗传协方差结构上。

通过估计训练集的方差分量，预测验证集的表型为:

5验证集是通过从1 - 406 (y的长度)中随机抽样41个值创建的，并重复这个抽样5次。验证集保存在验证矩阵中。这个矩阵指定GREML分析中要使用的数据行。
由于染色体2L和3L的方差分量为0，因此只有染色体2R和3R作为y预测的特征。

n <- length(y)
validate <- replicate(5, sample(1:n, 41))
cvG <- greml(y = y, X = X, GRM = GF[c(2, 4)], validate = validate)
str(cvG)
library(knitr)
library(kableExtra)kable(cvG$accuracy[1:5,], caption = "Cross Validation Predictive Ability") %>%
kable_styling(full_width = F, position = "left")

交叉验证预测能力

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$theta
h2f2 <- 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)


遗传值直方图

layout(matrix(1:2,ncol=2))
hist(fitGF$g[,2], xlab = "genetic value", cex.lab=0.8, cex.axis=0.8,main = "REML, chromosome 2R", cex.main=0.8)
hist(fitGF$g[,4], xlab = "genetic value", cex.lab=0.8, cex.axis=0.8,main = "REML, chromosome 3R", cex.main=0.8)


表型的柱状图

hist(y, xlab = "hours", cex.lab=0.8, cex.axis=0.8, main = "Starvation Resistance", cex.main=0.8)

三、基因集富集分析

在本文中，一系列不同的基因集富集分析方法将通过果蝇遗传参考面板(DGRP)提供的“抗饥饿”表型和基因型数据来说明。
qgg中的gsea功能可以进行四种不同的基因集富集分析。在这里，我们展示如何执行我们自己开发的协方差关联测试(CVAT)，基于分数的基于求和和基于计数的(hyperG)基因集富集分析。有关试验统计数据的详细信息和比较，请查阅doi:10.1534/genetics.116.189498。
一般的程序是获得单一的标记统计(例如总结统计)，从它可以计算和评估一组遗传标记的测试统计，以衡量标记集和表现型之间的联合程度。标记集由基因、生物途径、基因相互作用、基因表达谱等基因组特征来定义。

要进行不同类型的基因集富集分析，以下输入数据是必不可少的。
y : vector of phenotype
X: design matrix for covariables
Genomic feature marker sets
W: centered and scaled genotype matrix
G: genomic relationship matrix for each feature

gblup衍生的标记集检验方法基于三个步骤:首先拟合一个线性混合模型，然后计算单个标记的效果，最后计算标记集(基因组特征)的检验统计量。

该脚本包括以下步骤:
1)加载和准备数据，
2)GBLUP线性混合模型，限制性最大似然分析，
3)CVAT,
4)基于分数的，
5)基于累加的，
6)基于计数的基因集富集分析。

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

加载和准备数据

加载抗饥饿数据的已清理和编辑的数据框。

load(file = "./phenotypes/starv_inv_wo.Rdata")
dim(starv)
## [1] 406  20
head(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 <- starv
y <- data$y
length(y)
## [1] 406

为共变量准备设计矩阵X。fm是设计矩阵中包含相关协变量的公式。

fm <- y ~  wo + In2Lt + In2RNS + In3RP + In3RK + In3RMo
X <- model.matrix(fm, data=data)
X[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 1725755
W[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

基于基因本体(GO)加载SNP集、goset。

load(file="./annotation/goSets2.Rdata")
str(goSets[1:5])
## List of 5
##  $ GO:0000022: chr [1:315] "2R_2636861" "2R_2637023" "2R_2637293" "2R_2637711" ...
##  $ GO:0000027: chr [1:644] "X_1771396" "X_1772294" "X_1772337" "X_1772796" ...
##  $ GO:0000028: chr [1:283] "X_1376457" "X_1376829" "X_1376835" "X_9448856" ...
##  $ GO:0000045: chr [1:1602] "X_7208847" "X_7209453" "X_7209470" "X_7209504" ...
##  $ GO:0000060: chr [1:1816] "X_2056263" "X_2056399" "X_2056417" "X_2056810" ...

只保留W矩阵中的snp。

goSets <- lapply(goSets,function(x){ x[x%in%colnames(W)] })

创建对象nset:集合的数量和nsnp(每个GO项的snp数量)。

nsets <- length(goSets)
nsnps <- sapply(goSets,length

加性基因组关系矩阵G (VanRaden PM)。2008. 利用所有遗传标记G=WW ’ /m，其中W为中心和缩放的基因型矩阵，m为标记总数。

从基因型矩阵W计算累加的基因组关系矩阵G。

L <- data$L
G <- grm(W = W)
G <- G[L,L]

测定了每个系雄性和雌性的抗饥饿能力。因此G的维数是406行和406列。

str(G)
##  num [1:406, 1:406] 1.0327 0.006 0.0889 0.0128 -0.0133 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : chr [1:406] "100" "101" "105" "109" ...
##   ..$ : chr [1:406] "100" "101" "105" "109" ...

GBLUP线性混合模型

GBLUP基于线性混合模型，仅包含一种随机基因组效应:
y=Xb+Zg+e
其中y为表型观察向量，X和Z为固定效应和随机效应的设计矩阵，b为固定效应向量，g为所有遗传标记捕获的基因组值向量，e为残差向量。随机基因组值和残差被假设为独立的正态分布值，描述如下:g ~ N(0, g ggg2g)和e ~ N(0,I = bx2e)。因此，我们假设观察到的表型为y ~ N(Xb,V)，其中V=ZGZ ’ g2g +I g2e。

限制最大似然(REML)分析

使用REML分析(greml函数)估计方差分量,对于用于GBLUP分析的线性混合模型中的随机效应，分别为两个型号的2g和2e。
根据在整个研究群体中观察到的表型来预测遗传价值:

表型预测为:

greml函数在收敛之前要经过多次迭代(即，连续轮之间的参数变化小于指定的阈值，请参阅greml帮助页面中的“tol”参数)。在本例中，每次迭代都返回方差组件(第3列)和(第4列)的值.
greml函数返回一个列表结构，其中包括对固定效果(b)、随机效果(g)和剩余效果(e)的估计。列表中的其他值在greml帮助页面中描述。

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

单一标记的效应

单一标记效应的计算建立在gsea函数中。
单标记效应s^由GBLUP模型预测的总基因组值

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

计算出单标记效应的(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.6866151
mma <- 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值。

qgg包（续）-小数据教程-数据分析相关推荐

  1. qgg包-续2-大数据集教程

    qgg包-续2-大数据集教程 目录 qgg包-续2-大数据集教程 使用1000基因组项目数据的教程 下载1000 基因组数据 准备qgg的基因型数据 摸拟数据 计算基因组关系矩阵(GRM) 使用GRM ...

  2. 贝叶斯软件genle教程_手把手教你用R的gemtc包对生存数据进行贝叶斯网状Meta分析...

    大家好,本教程将介绍如何使用R的gemtc包对生存数据(HR为效应量)进行贝叶斯网状Meta分析. 前提条件: 需要下载R软件(推荐使用的R版本为3.5.3),以及RStudio(一个R的友好交互界面 ...

  3. linux tcp 包大小,linux – 通过大量连接和小数据包流量高的千兆网络提高TCP性能...

    我正在尝试通过"具有大量连接和小数据包流量的千兆网络"来提高TCP吞吐量.我的服务器操作系统是Ubuntu 11.10 Server 64bit. 有大约50.000(和不断增长的 ...

  4. zstd优秀的数据压缩算法，大数据小数据包

    zstd是Facebook在2016年开源的新无损压缩算法,优点是压缩率和压缩/解压缩性能都很突出. 在我们测试的文本日志压缩场景中,压缩率比gzip提高一倍,压缩性能与lz4.snappy相当甚至更 ...

  5. 微信小程序缓存获取数据教程

    微信小程序缓存获取数据教程 每个微信小程序都可以有自己的本地缓存,可以通过 wx.setStorage(wx.setStorageSync).wx.getStorage(wx.getStorageSy ...

  6. TCP/IP详解--TCP传输小数据包效率问题

    摘要:当使用TCP传输小型数据包时,程序的设计是相当重要的.如果在设计方案中不对TCP数据包的 延迟应答,Nagle算法,Winsock缓冲作用引起重视,将会严重影响程序的性能.这篇文章讨论了这些 问 ...

  7. R统计绘图-rgbif包下载GBIF数据及绘制分布图

    1 基本信息 博士退学前,做完斑马鱼的Phylogenomics分析,系统进化树冲突.基因流.ILS和种群历史动态等分析了之后,需要看一下Danio属物种的地理分布,希望能跟Phylogenomics ...

  8. 模型 标签数据 神经网络_大型神经网络和小数据的模型选择

    模型 标签数据 神经网络 The title statement is certainly a bold claim, and I suspect many of you are shaking yo ...

  9. pb程序怎么发布到iis_怎么使用抖音小程序第三方平台系统开发制作发布抖音小程序+教程...

    怎么使用抖音小程序第三方平台系统开发制作发布抖音小程序+教程 抖音短视频APP发布<2019年抖音数据报告>显示,其日活跃用户已经于2020年1月达到4亿,抖音APP受到广大用户追捧的同时 ...

  10. python抓取数据包_python抓数据包

    广告关闭 腾讯云11.11云上盛惠 ,精选热门产品助力上云,云服务器首年88元起,买的越多返的越多,最高返5000元! 前言:数据科学越来越火了,网页是数据很大的一个来源. 最近很多人问怎么抓网页数据 ...

最新文章

  1. 如何卸载office201032位_微软官方安装卸载修复工具、恶意软件删除工具，了解下！...
  2. springboot分页展示功能_springboot+vue实现分页功能
  3. 创建支持依赖注入、Serilog 日志和 AppSettings 的 .NET 5 控制台应用
  4. 没有主清单属性_原神：晴知的主C诺艾尔大型进阶攻略初版
  5. C++TCP和UDP属于传输层协议
  6. 70行Python代码，获取中国数据库大会（DTCC）全部PPT
  7. Unity面试题精选（4）
  8. editplus 快捷键及设置tab空白符及删除空格空行
  9. shell编程四剑客之 grep
  10. linux之shell快速入门系列＜8＞ | shell工具cut、sed、awk、sort
  11. 安卓手机屏幕投射电脑能同步声音
  12. MER：高通量测序应用于病原体和害虫诊断——综述与实用性建议
  13. python模型保存与恢复_tensorflow1.0学习之模型的保存与恢复(Saver)
  14. [开发过程]＜项目管理＞TAPD工具
  15. JavaScript高级程序设计读书笔记(第6章面向对象的程序设计之创建对象)
  16. CDbCriteria CArrayDataProvider zii.widgets.grid （1）
  17. php 递归 递归方式与算法
  18. STM32 SWD下载口无法下载的原因和解决办法
  19. 重积分 | 二重积分中 dx x dy = ρ dρ x dθ
  20. Ypai windows部署小白也可以哦

热门文章

  1. 阜阳师范学院java,刘冬冬 - 阜阳师范学院 - 计算机与信息工程学院
  2. 简单实验uwsgi+flask 部署caffe模型
  3. winedit 永久试用的办法
  4. python处理wrf气象数据_气象编程 | Python3之WRF的投影转换
  5. PS3主机今日发售 附官方问答
  6. 理财笔记 - 给朋友的建议
  7. OpenAI 最强对话模型 ChatGPT 注册使用笔记
  8. warmup与余弦退火学习率
  9. html制作一个视频播放器,H5 打造属于自己的视频播放器（HTML 篇）
  10. 疯狂原始人服务器维修,疯狂原始人服务器互通详解 安卓和ios能一起玩么