基于MTAR的罕见变量关联汇总统计量多性状分析
王羡金
2024-09-12
#安装MTAR包
#install.packages("MTAR")
library(MTAR)MTAR包具有以下功能：
1)计算不同类型数据输入的汇总统计量；
2)在必要时调整样本重叠的汇总统计量；
3)使用汇总统计量进行多性状罕见变异关联检验。
这张图展示了MTAR的工作流程，分为三个主要步骤，用于多性状的遗传效应分析。下面是对图中每个部分的解释：
Step 1: 准备多性状的汇总统计数据
功能：Get_UV_from_data()、Get_UV_from_var()和Get_UV_from_beta()
这些函数用于从不同来源获取性状的汇总统计数据。
输出是 \(U\) 和 \(V\)，代表不同性状的估计效应和方差，作为后续分析的输入。
Step 2: 计算性状之间 Z 分数的协方差
功能：
使用 Get_zeta() 函数，从重叠样本中计算不同性状之间的 Z
分数协方差，输出 \(zeta\) 值。
同时还需要考虑几种因素来计算不同协方差矩阵：
MAF（次要等位基因频率）
rho.SNP（单核苷酸多态性相关性）
Genetic correlation（遗传相关性）
rho.trait（性状相关性）
MAF和rho.SNP构建 SNPs 间效应的协方差矩阵 \(B_1\) ，Genetic correlation和rho.trait构建 性状间效应的协方差矩阵 \(B_2\)。
\(B_1\) 和\(B_2\) 合并得到总的效应协方差矩阵 \(B\)。
中间步骤：效应估计
效应估计：
通过第一步输出的 \(U\) 和 \(V\) 值，计算效应估计 \(\hat{\beta} = (\hat{\beta}_1, ..., \hat{\beta}_K)\)，其中 \(\hat{\beta}_K\) 表示第 \(K\) 个性状的效应估计。
协方差矩阵
估计值协方差矩阵 \(\Sigma\)：
将 \(zeta\) 和效应估计 \(\hat{\beta}\) 一起用于计算估计值协方差矩阵 \(\Sigma\)。
Step 3: 进行 MTAR 分析
MTAR 分析：
最后一步是调用 MTAR() 函数，输出最终的 p 值结果。
最终输出：
Final p-values：通过多性状的联合分析得出最终的 p 值，用于判定基因变异是否与多性状有显著关联。
这个流程将多性状和多SNP的效应及其关联整合在一起，以更全面地评估遗传风险。
Step1：为多个特征准备汇总统计信息
在 MTAR中，U 统计量 和它们的 协方差矩阵 V 是必需的汇总统计数据。本文将展示如何使用不同的实用函数来准备这些汇总统计数据。这些实用函数的输出结果可以直接作为 MTAR 函数的输入，用于进行遗传关联分析。
1.给定个体层面的数据
（1）数据和模型：
研究目的是分析m个稀有变异（Minor Allele Frequency, MAF < 5%）的效应，这些稀有变异在某个基因或一组SNP上对K个性状产生影响。
对于第k个性状，当个体水平的基因型\(G_{ik}\)、协变量\(X_{ik}\)和性状数据\(Y_{ik}\)可以得到时，可以通过广义线性模型获得得分统计量\(U_{k}\)及其协方差估计\(V_{k}\)。
广义线性模型：
广义线性模型的似然函数：
\(\exp\left\{\frac{Y_{ik}\left(\beta_{k}^{\mathrm{T}}\mathbf{G}_{ik}+\gamma_{k}^{\mathrm{T}}\mathbf{X}_{ik}\right)-b(\beta_{k}^{\mathrm{T}}\mathbf{G}_{ik}+\gamma_{k}^{\mathrm{T}}\mathbf{X}_{ik})}{a(\phi_{k})}+c(Y_{ik},\phi_{k})\right\}\)
其中\(\beta_{k}\)和\(\gamma_{k}\)是回归参数，\(\phi_{k}\)是离散参数，a、b和c是特定函数，一阶和二阶导数分别用\(b^{\prime}\)和\(b^{\prime\prime}\)表示。
（2）得分统计量和协方差估计：
得分统计量\(U_{k}\)和协方差估计\(V_{k}\)的计算公式也给出：
\(\mathbf{U}_k=a(\widehat{\phi}_k)^{-1}\sum_{i=1}^n\left\{Y_{ik}-b^{\prime}(\widehat{\gamma}_k^T\mathbf{X}_{ik})\right\}\mathbf{G}_{ik}\)
\(\mathbf{V}_{k}=a(\widehat{\phi}_{k})^{-1}\Big[\sum_{i=1}^{n}b^{\prime\prime}(\widehat{\gamma}_{k}^{T}\mathbf{X}_{ik})\mathbf{G}_{ik}\mathbf{G}_{ik}^{T}-\Big\{\sum_{i=1}^{n}b^{\prime\prime}(\widehat{\gamma}_{k}^{T}\mathbf{X}_{ik})\mathbf{G}_{ik}\mathbf{X}_{ik}^{T}\Big\}\Big\{\sum_{i=1}^{n}b^{\prime\prime}(\widehat{\gamma}_{k}^{T}\mathbf{X}_{ik})\mathbf{X}_{ik}\mathbf{X}_{ik}^{T}\Big\}^{-1}\Big\{\sum_{i=1}^{n}b^{\prime\prime}(\widehat{\gamma}_{k}^{T}\mathbf{X}_{ik})\mathbf{X}_{ik}\mathbf{G}_{ik}^{T}\Big\}\Big]\)
其中\(\widehat{\gamma}_{k}\)和\(\widehat{\phi}_{k}\)是在\(H_{0}:\beta_{k}=0\)下的\(\gamma_{k}\)和\(\phi_{k}\)的约束最大似然估计量。
（3）线性回归和逻辑回归模型：
对于线性回归模型和逻辑回归模型，给出了\(a(\widehat{\phi}_{k})\)、\(b^{\prime}(z)\)和\(b^{\prime\prime}(z)\)的具体形式。
线性回归模型：
\(\begin{aligned}&a(\widehat{\phi}_{k})=n^{-1}\sum_{i=1}^{n}(Y_{ik}-\widehat{\gamma}_{k}^{T}\mathbf{X}_{ik})^{2}\\\\&b^{\prime}(z)=z\\\\&b^{\prime\prime}(z)=1\end{aligned}\)
逻辑回归模型：
\(\begin{aligned}&a(\phi)=1\\&b^{\prime}(z)=e^z/(1+e^z)\\&b^{\prime\prime}(z)=e^z/(1+e^z)^2\end{aligned}\)
从个体水平数据计算汇总统计量：
总结数据的个体水平信息后，可以使用函数 Get_UV_from_data
来计算模型中需要的统计量，这个函数需要需要输入性状、协变量和基因型信息。
每个研究的性状、协变量和基因型信息应该是\(n \times K\)、\(n \times D\)和\(n \times m\)矩阵，其中n是样本量，K是性状数，D是协变量数，m是变异数。
每项研究的性状数量可以不同，但必须提供性状名称。与基因型中的 SNP 名称相同。性状或基因型中的缺失值应标注为 NA。性状、协变量和基因型中受试者 ID 的顺序应一致。
模拟数据集：
有三个研究，每个研究都有3个连续性状和10个罕见变异。
研究1有1500名受试者，但每个受试者只有一个性状测量。
研究2和研究3的样本量为500，每个受试者有两个或三个性状测量。
下面用一个模拟数据集 rawdata 来展示如何应用Get_UV_from_data 函数。在 rawdata 中，有 3 项研究、3 个连续性状和 10 个罕见变异。图 2 显示了三项研究中性状间样本重叠的示意图。具体来说，研究 1 有 1500 个受试者，但每个受试者只有一个性状测量值。在研究 2 和研究 3 中，样本量为 500，每个受试者有两个（低密度脂蛋白LDL和高密度脂蛋白HDL）或三个性状（低密度脂蛋白LDL、高密度脂蛋白HDL和总胆固醇TG）测量值。
data("rawdata")
#attach函数将rawdata数据集附加到当前的R环境中。这样做可以直接使用数据集中的变量名来引用数据，而不需要每次都使用数据集的名称作为前缀。
attach(rawdata)
#detach函数可以将附加的数据从搜索路径上移除。
#detach(rawdata)血脂水平研究(LDL、HDL、TG是血液中的脂质成分，通常在血脂检查时测量。)
LDL:低密度脂蛋白（Low-Density Lipoprotein）
HDL:高密度脂蛋白（High-Density Lipoprotein）
TG:和甘油三酯（Triglycerides）
LDL 被认为是“坏”胆固醇，因为它在血管壁上形成斑块，可能导致动脉粥样硬化。HDL 被认为是“好”胆固醇，因为它有助于从血管壁上清除胆固醇，减少心血管疾病的风险。TG 是血液中的脂肪，高甘油三酯水平也可能增加心血管疾病的风险。
每个“SUBJ”代表一个受试者，后面的数值代表他们的LDL、HDL和TG的测量值。HDL和TG的值被标记为“NA”，这通常意味着这些值不可用或未被测量。
str(rawdata)## List of 3
## $ traits.dat:List of 3
## ..$ Study1: num [1:1500, 1:3] -2.617 0.852 1.244 -0.533 0.801 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:1500] "SUBJ1" "SUBJ2" "SUBJ3" "SUBJ4" ...
## .. .. ..$ : chr [1:3] "LDL" "HDL" "TG"
## ..$ Study2: num [1:500, 1:2] -0.479 2.151 -1.679 -1.531 0.62 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:500] "SUBJ1" "SUBJ2" "SUBJ3" "SUBJ4" ...
## .. .. ..$ : chr [1:2] "LDL" "HDL"
## ..$ Study3: num [1:500, 1:3] 0.0945 0.2999 0.7054 -0.5818 -0.8408 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:500] "SUBJ1" "SUBJ2" "SUBJ3" "SUBJ4" ...
## .. .. ..$ : chr [1:3] "LDL" "HDL" "TG"
## $ geno.dat :List of 3
## ..$ Study1: num [1:1500, 1:10] 0 0 0 0 0 0 0 0 0 0 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:1500] "SUBJ1" "SUBJ2" "SUBJ3" "SUBJ4" ...
## .. .. ..$ : chr [1:10] "SNP2148" "SNP2149" "SNP2150" "SNP2151" ...
## ..$ Study2: num [1:500, 1:10] 0 0 0 0 0 0 0 0 0 0 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:500] "SUBJ1" "SUBJ2" "SUBJ3" "SUBJ4" ...
## .. .. ..$ : chr [1:10] "SNP2148" "SNP2149" "SNP2150" "SNP2151" ...
## ..$ Study3: num [1:500, 1:10] 0 0 0 0 0 0 0 0 0 0 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:500] "SUBJ1" "SUBJ2" "SUBJ3" "SUBJ4" ...
## .. .. ..$ : chr [1:10] "SNP2148" "SNP2149" "SNP2150" "SNP2151" ...
## $ cov.dat :List of 3
## ..$ Study1: num [1:1500, 1:2] 1 0 0 0 1 1 0 0 1 1 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:1500] "SUBJ1" "SUBJ2" "SUBJ3" "SUBJ4" ...
## .. .. ..$ : chr [1:2] "cov1" "cov2"
## ..$ Study2: num [1:500, 1:2] 1 0 0 0 0 1 0 0 0 0 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:500] "SUBJ1" "SUBJ2" "SUBJ3" "SUBJ4" ...
## .. .. ..$ : chr [1:2] "cov1" "cov2"
## ..$ Study3: num [1:500, 1:2] 1 1 1 0 1 0 0 0 0 1 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:500] "SUBJ1" "SUBJ2" "SUBJ3" "SUBJ4" ...
## .. .. ..$ : chr [1:2] "cov1" "cov2"表型数据（traits.dat）：一个数字列表，每个子列表包含每项研究的性状信息。在每项研究的 n × K 矩阵，每一行代表一个独立个体，每一列代表一个独立性状。如果受试者 i 没有性状 k，则相应的值设为 NA。每项研究的性状数量可以不同，但性状名称是必需的。
head(traits.dat$Study1)## LDL HDL TG
## SUBJ1 -2.6166166 NA NA
## SUBJ2 0.8519045 NA NA
## SUBJ3 1.2438145 NA NA
## SUBJ4 -0.5327296 NA NA
## SUBJ5 0.8006589 NA NA
## SUBJ6 0.1560920 NA NA协变量数据（cov.dat）：数字列表，每个子列表包含每项研究的协变量信息。每项研究的
n × D 矩阵，每行代表一个独立个体，每列代表一个协变量。
协变量是指不同干预措施前，研究者预计的、会对主要变量分析产生重要影响的因素。这类变量可以是定性，也可以是定量的，抑或是等级资料。其可以是人口统计学指标如年龄、体重、种族/民族等；也可以是一些疾病预后因素如疾病分型、病程或病情严重程度；当然还有一些其它因素如研究中心或研究者等。上述随机对照研究的基线信息均可以看作协变量，是随机对照研究设计和分析时必须要考虑的问题之一。
head(cov.dat$Study1)## cov1 cov2
## SUBJ1 1 -0.29713414
## SUBJ2 0 -0.08074092
## SUBJ3 0 0.55975760
## SUBJ4 0 0.49851956
## SUBJ5 1 -0.07097778
## SUBJ6 1 -0.28825662基因型数据（geno.dat）：数字列表，每个子列表包含每个研究的基因型信息。在每项研究中，数字n × m矩阵，每一行为独立个体，每列为SNP。每个基因型应编码为0、1、2和NA，分别代表AA、AA、AA和缺失，其中A和a分别代表主等位基因和次等位基因。每个研究中的snp数量可以不同，但snp的名称是必需的。此外，在基因型、协变量和性状列表中的研究数量必须相同。在每项研究中，受试者ID在性状、协变量和基因型之间的顺序必须相同。
head(geno.dat$Study1)## SNP2148 SNP2149 SNP2150 SNP2151 SNP2152 SNP2153 SNP2154 SNP2155 SNP2156
## SUBJ1 0 0 0 0 0 0 0 0 0
## SUBJ2 0 0 0 0 0 0 0 0 0
## SUBJ3 0 0 0 0 0 0 0 0 0
## SUBJ4 0 0 0 0 0 0 0 0 0
## SUBJ5 0 0 0 0 0 0 0 0 0
## SUBJ6 0 0 0 0 0 0 0 0 0
## SNP2157
## SUBJ1 0
## SUBJ2 0
## SUBJ3 0
## SUBJ4 0
## SUBJ5 0
## SUBJ6 0obs.stat <- Get_UV_from_data(traits = traits.dat,
covariates = cov.dat,
genotype = geno.dat,
covariance = TRUE)
obs.stat$U## $LDL
## SNP2148 SNP2149 SNP2150 SNP2151 SNP2152 SNP2153 SNP2154
## 0.3438038 1.5488022 1.3252516 0.0000000 2.0216781 -1.8258464 0.6500477
## SNP2155 SNP2156 SNP2157
## 0.0000000 0.1757972 1.5709795
##
## $HDL
## SNP2148 SNP2149 SNP2150 SNP2151 SNP2152 SNP2153
## 0.00000000 0.70525844 1.25507643 0.00000000 -0.03924912 -0.28306623
## SNP2154 SNP2155 SNP2156 SNP2157
## 2.06966868 0.00000000 -12.48457057 0.00000000
##
## $TG
## SNP2148 SNP2149 SNP2150 SNP2151 SNP2152 SNP2153 SNP2154
## 0.0000000 1.3061250 2.2535334 0.0000000 -2.6934948 -0.2769778 1.5616663
## SNP2155 SNP2156 SNP2157
## 0.0000000 -0.3394679 0.0000000detach(rawdata)2.给定分数统计及其方差
当得分统计\(\mathbf{U}_{k}\)及其方差\((\mathrm{diag}(\mathbf{V}_k))\)可用时，协方差矩阵\(\mathbf{V}_{k}\)可以近似为：
\(\mathbf{V}_{k}\approx\mathrm{diag}(\mathbf{V}_{k})^{1/2}\mathbf{R}\mathrm{diag}(\mathbf{V}_{k})^{1/2}\)
其中\(\mathbf{R}=\{R_{j\ell}\}_{j,\ell=1}^{m}\)是
SNP 相关矩阵，估计方法可以基于实际基因型数据或者外部参考（如 Hu
等人，2013 年的研究）。这种变换已经在函数Get_UV_from_varU
中实现。
varU.example：汇总统计量及其方差
给定汇总统计量及其方差，计算 MTAR 的 U 和 V 的示例。
data("varU.example")U：汇总统计列表，每个子列表是每个研究的 m × K 矩阵
varU.example[["U"]]## [[1]]
## LDL HDL TG
## 1:150551327 -268.5912258 229.0655499 -363.297279
## 1:150551576 31.2170161 21.4970739 -11.531591
## 1:150551649 -12.7574302 -18.7336436 -1.618765
## 1:150550965 4.2199502 -0.1519353 -2.521649
## 1:150551447 0.8814878 -3.7007466 1.127816
## 1:150550761 3.2751963 -12.8477062 6.706268varU：汇总统计方差列表，每个子列表是每个研究的 m × K 矩阵
varU.example[["varU"]]## [[1]]
## LDL HDL TG
## 1:150551327 5744.43131 6220.567347 5898.995047
## 1:150551576 490.20387 523.477481 505.804794
## 1:150551649 235.93806 249.270743 247.407921
## 1:150550965 76.26820 80.808048 77.161844
## 1:150551447 4.99751 5.997364 5.997364
## 1:150550761 174.29786 196.580663 181.299473R：所有研究中出现的 SNPs 的 SNP 关联矩阵
varU.example[["R"]]## 1:150551327 1:150551576 1:150551649 1:150550965 1:150551447
## 1:150551327 1.000000000 0 0 -0.007443092 0
## 1:150551576 0.000000000 1 0 0.000000000 0
## 1:150551649 0.000000000 0 1 0.000000000 0
## 1:150550965 -0.007443092 0 0 1.000000000 0
## 1:150551447 0.000000000 0 0 0.000000000 1
## 1:150550761 -0.004294482 0 0 -0.001124344 0
## 1:150550761
## 1:150551327 -0.004294482
## 1:150551576 0.000000000
## 1:150551649 0.000000000
## 1:150550965 -0.001124344
## 1:150551447 0.000000000
## 1:150550761 1.000000000data("varU.example")
attach(varU.example)
obs.stat <- Get_UV_from_varU(U = U, varU = varU, R= R)#输出的内容显示了每个性状对应的 SNP 的得分统计值
#1:150551327表示染色体1的位置150551327
#-268.5912258表示对应的得分统计值
obs.stat$U## $LDL
## 1:150551327 1:150551576 1:150551649 1:150550965 1:150551447 1:150550761
## -268.5912258 31.2170161 -12.7574302 4.2199502 0.8814878 3.2751963
##
## $HDL
## 1:150551327 1:150551576 1:150551649 1:150550965 1:150551447 1:150550761
## 229.0655499 21.4970739 -18.7336436 -0.1519353 -3.7007466 -12.8477062
##
## $TG
## 1:150551327 1:150551576 1:150551649 1:150550965 1:150551447 1:150550761
## -363.297279 -11.531591 -1.618765 -2.521649 1.127816 6.706268detach(varU.example)3.给定遗传效应估计值及其标准误差
当遗传效应估计值 \(\hat{\beta}_k = \{\hat{\beta}_{kj}\}_{j=1}^{m}\) 和它们的标准误差 \(se_{k} = \{se_{kj}\}_{j=1}^{m}\) 可用时，我们可以通过以下方式近似得分统计量 \(U_k = \{U_{kj}\}_{j=1}^{m}\) 和协方差矩阵 \(V_k = \{V_{kj\ell}\}_{j,\ell=1}^{m}\)
得分统计量 \(U_k\) 的近似：
\(U_{kj}\approx\widehat{\beta}_{kj}/se_{kj}^2\)
这里，\(U_{kj}\) 是第k个性状的第j个SNP的得分统计量。这个近似是通过将每个SNP的遗传效应估计值除以其标准误差的平方得到的。
协方差矩阵 \(V_k\) 的近似：
\(V_{kj\ell}\approx R_{j\ell}/(se_{kj}se_{k\ell})\)
这里，\(V_{kj\ell}\) 是第k个性状的第j个和第\(\ell\)个SNP之间的协方差。这个近似是通过将SNP相关性矩阵 \(R\) 中的元素 \(R_{j\ell}\) 除以两个SNP的标准误差的乘积得到的。
data("beta.example")beta.example[["Beta"]]## [[1]]
## LDL HDL TG
## 1:150551327 -0.0467568 0.0368239 -0.0615863
## 1:150551576 0.0636817 0.0410659 -0.0227985
## 1:150551649 -0.0540711 -0.0751538 -0.0065429
## 1:150550965 0.0553304 -0.0018802 -0.0326800
## 1:150551447 0.1763854 -0.6170622 0.1880520
## 1:150550761 0.0187908 -0.0653559 0.0369900beta.example[["Beta.se"]]## [[1]]
## LDL HDL TG
## 1:150551327 0.013194 0.012679 0.013020
## 1:150551576 0.045166 0.043707 0.044464
## 1:150551649 0.065103 0.063338 0.063576
## 1:150550965 0.114506 0.111243 0.113841
## 1:150551447 0.447325 0.408338 0.408338
## 1:150550761 0.075745 0.071323 0.074268beta.example[["R"]]## 1:150551327 1:150551576 1:150551649 1:150550965 1:150551447
## 1:150551327 1.000000000 0 0 -0.007443092 0
## 1:150551576 0.000000000 1 0 0.000000000 0
## 1:150551649 0.000000000 0 1 0.000000000 0
## 1:150550965 -0.007443092 0 0 1.000000000 0
## 1:150551447 0.000000000 0 0 0.000000000 1
## 1:150550761 -0.004294482 0 0 -0.001124344 0
## 1:150550761
## 1:150551327 -0.004294482
## 1:150551576 0.000000000
## 1:150551649 0.000000000
## 1:150550965 -0.001124344
## 1:150551447 0.000000000
## 1:150550761 1.000000000data("beta.example")
attach(beta.example)
obs.stat <- Get_UV_from_beta(Beta = Beta, Beta.se = Beta.se, R = R)obs.stat## $U
## $U$LDL
## 1:150551327 1:150551576 1:150551649 1:150550965 1:150551447 1:150550761
## -268.5912258 31.2170161 -12.7574302 4.2199502 0.8814878 3.2751963
##
## $U$HDL
## 1:150551327 1:150551576 1:150551649 1:150550965 1:150551447 1:150550761
## 229.0655499 21.4970739 -18.7336436 -0.1519353 -3.7007466 -12.8477062
##
## $U$TG
## 1:150551327 1:150551576 1:150551649 1:150550965 1:150551447 1:150550761
## -363.297279 -11.531591 -1.618765 -2.521649 1.127816 6.706268
##
##
## $V
## $V$LDL
## 1:150551327 1:150551576 1:150551649 1:150550965 1:150551447
## 1:150551327 5744.431308 0.0000 0.0000 -4.9266153 0.00000
## 1:150551576 0.000000 490.2039 0.0000 0.0000000 0.00000
## 1:150551649 0.000000 0.0000 235.9381 0.0000000 0.00000
## 1:150550965 -4.926615 0.0000 0.0000 76.2682027 0.00000
## 1:150551447 0.000000 0.0000 0.0000 0.0000000 4.99751
## 1:150550761 -4.297149 0.0000 0.0000 -0.1296334 0.00000
## 1:150550761
## 1:150551327 -4.2971487
## 1:150551576 0.0000000
## 1:150551649 0.0000000
## 1:150550965 -0.1296334
## 1:150551447 0.0000000
## 1:150550761 174.2978621
##
## $V$HDL
## 1:150551327 1:150551576 1:150551649 1:150550965 1:150551447
## 1:150551327 6220.567347 0.0000 0.0000 -5.2771044 0.000000
## 1:150551576 0.000000 523.4775 0.0000 0.0000000 0.000000
## 1:150551649 0.000000 0.0000 249.2707 0.0000000 0.000000
## 1:150550965 -5.277104 0.0000 0.0000 80.8080478 0.000000
## 1:150551447 0.000000 0.0000 0.0000 0.0000000 5.997364
## 1:150550761 -4.748935 0.0000 0.0000 -0.1417088 0.000000
## 1:150550761
## 1:150551327 -4.7489350
## 1:150551576 0.0000000
## 1:150551649 0.0000000
## 1:150550965 -0.1417088
## 1:150551447 0.0000000
## 1:150550761 196.5806630
##
## $V$TG
## 1:150551327 1:150551576 1:150551649 1:150550965 1:150551447
## 1:150551327 5898.995047 0.0000 0.0000 -5.0216182 0.000000
## 1:150551576 0.000000 505.8048 0.0000 0.0000000 0.000000
## 1:150551649 0.000000 0.0000 247.4079 0.0000000 0.000000
## 1:150550965 -5.021618 0.0000 0.0000 77.1618437 0.000000
## 1:150551447 0.000000 0.0000 0.0000 0.0000000 5.997364
## 1:150550761 -4.441177 0.0000 0.0000 -0.1329838 0.000000
## 1:150550761
## 1:150551327 -4.4411774
## 1:150551576 0.0000000
## 1:150551649 0.0000000
## 1:150550965 -0.1329838
## 1:150551447 0.0000000
## 1:150550761 181.2994733detach(beta.example)Step2：计算重叠样本中性状间 Z 值的协方差
如果所有的性状都来自一个研究或来自多个研究重叠的受试者，则性状之间的遗传效应估计是相关的。为了计算这些遗传效应的协方差，我们首先需要计算zeta与性状不相关的SNP的Z分数的性状间协方差，在Get_zeta函数中，应该输入一个Z值的向量，这些Z值来自常见变体和每个性状之间的全基因组关联测试。该函数将进行LD修剪（通过使用1000Genome中的参考独立SNP列表），删除p值<0.05（默认p值截止阈值）的SNP，并计算zeta。
下面展示一个估计矩阵zeta的简化示例，从全球脂质遗传学联合会的1号染色体中提取了737个无效和常见SNPs子集。过滤掉这些非独立的SNP后，只剩下137个SNP。
data("zeta.example")
attach(zeta.example)head(zeta.example[["Zscore"]][["LDL"]])## 1:114680310 1:169533266 1:219743719 1:71439590 1:55738663 1:193038254
## -1.959682 -1.955302 1.953341 1.952537 -1.948679 1.948633head(zeta.example[["Zscore"]][["HDL"]])## 1:151315287 1:159284385 1:155162067 1:203207520 1:172431386 1:170934394
## 1.956903 -1.955491 1.950986 -1.949963 -1.949367 1.947483head(zeta.example[["Zscore"]][["TG"]])## 1:155014223 1:24406647 1:165664552 1:156106185 1:49589847 1:64615775
## 1.959510 1.959106 1.958727 -1.958624 1.956160 1.954811#从1000Genome数据集中下载独立的常见 SNPs
#复制https://raw.githubusercontent.com/lan/MTAR/master/indp_snps.1KG.rda到浏览器下载，放在工作目录
load("indp_snps.1KG.rda")zeta1 <- Get_zeta(Zscore = Zscore, Indp_common_snp = indp_snps.1KG)zeta1## LDL HDL TG
## LDL 1.00000000 -0.02027713 0.06516755
## HDL -0.02027713 1.00000000 0.14705052
## TG 0.06516755 0.14705052 1.00000000detach(zeta.example)Step3. 运行多性状罕见变异关联测试（MTAR）
要运行 MTAR 函数，需要输入得分统计量U、其协方差矩阵V和最小等位基因频率MAF。其他输入为可选项。这里的示例数据集（MTAR.example）提供了 MTAR 函数所需的所有信息。
data("MTAR.example")
names(MTAR.example)## [1] "annotation" "U" "V"
## [4] "MAF" "genetic_cor.trait" "zeta"MTAR.example[["annotation"]]head(MTAR.example[["U"]][["U1"]])## 11:823586 11:821676 11:824584 11:821766 11:824674 11:821778
## 182.289584 -103.571988 18.087092 -5.068688 -4.630540 -11.702198head(MTAR.example[["U"]][["U2"]])## 11:823586 11:821676 11:824584 11:821766 11:824674 11:821778
## -433.2736814 75.5062260 6.3578692 5.6986937 5.5695234 -0.9200333head(MTAR.example[["U"]][["U3"]])## 11:823586 11:821676 11:824584 11:821766 11:824674 11:821778
## 247.668896 -45.542607 6.998428 -7.301066 -2.261961 -5.250828head(MTAR.example[["V"]][["V1"]])## 11:823586 11:821676 11:824584 11:821766 11:824674 11:821778 11:824021
## 11:823586 7086.638 0.000 0.0000 0.00000 0.0000 0.0000 0
## 11:821676 0.000 3045.339 0.0000 0.00000 0.0000 0.0000 0
## 11:824584 0.000 0.000 282.3805 0.00000 0.0000 0.0000 0
## 11:821766 0.000 0.000 0.0000 31.89011 0.0000 0.0000 0
## 11:824674 0.000 0.000 0.0000 0.00000 32.8819 0.0000 0
## 11:821778 0.000 0.000 0.0000 0.00000 0.0000 219.3357 0
## 11:824711 11:824624 11:823729 11:823799 11:823742 11:822565 11:824042
## 11:823586 0 0 0 0 0 0 0
## 11:821676 0 0 0 0 0 0 0
## 11:824584 0 0 0 0 0 0 0
## 11:821766 0 0 0 0 0 0 0
## 11:824674 0 0 0 0 0 0 0
## 11:821778 0 0 0 0 0 0 0head(MTAR.example[["V"]][["V2"]])## 11:823586 11:821676 11:824584 11:821766 11:824674 11:821778 11:824021
## 11:823586 7540.44 0.000 0.0000 0.00000 0.00000 0.000 0
## 11:821676 0.00 3192.652 0.0000 0.00000 0.00000 0.000 0
## 11:824584 0.00 0.000 285.4903 0.00000 0.00000 0.000 0
## 11:821766 0.00 0.000 0.0000 33.22853 0.00000 0.000 0
## 11:824674 0.00 0.000 0.0000 0.00000 32.89736 0.000 0
## 11:821778 0.00 0.000 0.0000 0.00000 0.00000 242.433 0
## 11:824711 11:824624 11:823729 11:823799 11:823742 11:822565 11:824042
## 11:823586 0 0 0 0 0 0 0
## 11:821676 0 0 0 0 0 0 0
## 11:824584 0 0 0 0 0 0 0
## 11:821766 0 0 0 0 0 0 0
## 11:824674 0 0 0 0 0 0 0
## 11:821778 0 0 0 0 0 0 0head(MTAR.example[["V"]][["V3"]])## 11:823586 11:821676 11:824584 11:821766 11:824674 11:821778 11:824021
## 11:823586 7306.384 0.000 0.0000 0.00000 0.00000 0.0000 0
## 11:821676 0.000 3145.555 0.0000 0.00000 0.00000 0.0000 0
## 11:824584 0.000 0.000 285.4614 0.00000 0.00000 0.0000 0
## 11:821766 0.000 0.000 0.0000 32.81863 0.00000 0.0000 0
## 11:824674 0.000 0.000 0.0000 0.00000 34.89334 0.0000 0
## 11:821778 0.000 0.000 0.0000 0.00000 0.00000 224.3848 0
## 11:824711 11:824624 11:823729 11:823799 11:823742 11:822565 11:824042
## 11:823586 0 0 0 0 0 0 0
## 11:821676 0 0 0 0 0 0 0
## 11:824584 0 0 0 0 0 0 0
## 11:821766 0 0 0 0 0 0 0
## 11:824674 0 0 0 0 0 0 0
## 11:821778 0 0 0 0 0 0 0head(MTAR.example[["MAF"]])## 11:823586 11:821676 11:824584 11:821766 11:824674 11:821778
## 0.01267800 0.01715900 0.00071809 0.00937990 0.00019078 0.00041837head(MTAR.example[["genetic_cor.trait"]])## LDL HDL TG
## LDL 1.000 0.091 0.353
## HDL 0.091 1.000 -0.612
## TG 0.353 -0.612 1.000head(MTAR.example[["zeta"]])## LDL HDL TG
## LDL 1.00000000 -0.01893048 0.03298261
## HDL -0.01893048 1.00000000 -0.01622033
## TG 0.03298261 -0.01622033 1.00000000具体来说，genetic_cor.trait表示遗传相关性信息。遗传相关性可通过ldsc软件估算。一些网站还将许多复杂性状之间的遗传相关性存档。如果没有遗传相关性，MTAR将使用性状间的可交换相关结构。zeta可以在步骤2中通过Get_zeta估算。如果没有提供zeta，MTAR将假定性状间没有样本重叠。将返回MTAR-O、cMTAR、iMTAR和cctP的p值及辅助信息。
attach(MTAR.example)
pval <- MTAR(U = U, V = V, MAF = MAF, genetic_cor.trait = genetic_cor.trait, zeta = zeta)
pval## $MTARO
## [1] 5.226324e-08
##
## $cMTAR
## $cMTAR$p
## [1] 5.368409e-08
##
## $cMTAR$rho1.min
## [1] 0
##
## $cMTAR$rho2.min
## [1] 1
##
##
## $iMTAR
## $iMTAR$p
## [1] 2.599168e-08
##
## $iMTAR$rho1.min
## [1] 0
##
## $iMTAR$rho2.min
## [1] 1
##
##
## $cctP
## [1] 3.329077e-06MTAR-O、cMTAR、iMTAR和cctP的p值是多表型关联分析中不同统计测试的结果。每个p值代表了不同模型或方法下遗传变异与表型之间关联的显著性程度。
1. MTAR-O（MTAR Original）
- p值：
5.226324e-08 - 意义：MTAR-O是多表型关联测试的原始方法，用于检测稀有变异在多个表型上的加性效应。p值表示该遗传变异在多个表型上是否与这些表型显著关联。这里的p值非常小，表示该变异与表型的关联具有很高的显著性，表明稀有变异可能在不同的表型上共同起作用。
- 作用：MTAR-O模型适用于检测在所有表型上表现出一致遗传效应的变异。
2. cMTAR（Combined MTAR 或 Conditional MTAR）
- p值：
5.368409e-08 - 意义：cMTAR是一种结合或条件化的MTAR模型，允许稀有变异在不同表型上产生不同强度的效应。这个模型考虑了表型之间的差异，并在不同表型之间更灵活地捕捉遗传效应的变异。
- 作用：cMTAR能处理遗传效应在某些表型上强，而在其他表型上弱或不存在的情况。这个p值同样非常小，表明该稀有变异与至少部分表型之间存在显著关联。
3. iMTAR（Independent MTAR）
- p值：
2.599168e-08 - 意义：iMTAR模型假设稀有变异对不同表型产生独立效应，不依赖于其他表型的结果。p值反映该变异与表型的独立关联。
- 作用：iMTAR适合用来检测在不同表型上独立作用的变异。iMTAR的p值也非常显著，表明该变异可能与多个表型之间具有独立的关联效应。
4. cctP（Cauchy Combination Test）
- p值：
3.329077e-06 - 意义：cctP是使用柯西组合测试（Cauchy Combination Test, CCT）得到的p值，它是一种结合多个测试结果的统计方法，尤其在处理不同的p值时有良好表现。CCT方法通过整合多个p值来得出最终的显著性。
- 作用：cctP可以整合来自多个模型的p值，形成一个更稳健的结果，适用于不同表型中遗传效应复杂多变的情况。尽管p值稍高于其他测试，但依然在统计学上非常显著，表明该稀有变异与表型之间的关联可能是多个模型共同作用的结果。
总结：
- MTAR-O：适合检测在所有表型上表现出一致效应的变异。
- cMTAR：适合检测在部分表型上有不同效应的变异。
- iMTAR：适合检测表型间独立效应的变异。
- cctP：是一种组合测试，整合多个模型的p值，确保结果更加稳健。
这些不同测试方法提供了多角度分析遗传变异与多表型之间的关联，有助于捕捉复杂的遗传效应模式。在本文的结果中，这些p值都非常显著，表明这些稀有变异可能在多个表型上起着重要的作用。