处理 sommer 包中的不同性别

Dealing with separate sexes in the sommer package

我正在尝试弄清楚如何在 R 的 sommer 包中指定特定的混合模型。

我测量了一群男性个体的两个特征,以及一群女性亲属的两个特征。目的是估计这 4 个性状内部和之间的遗传(协)方差。但是,因为个体不能是两种性别,所以我想拟合模型,以便它估计女性特征 1 和女性特征 2 之间的残差协方差,以及男性特征 1 和男性特征 2 之间的残差协方差,但不是残差协方差在任何男性和女性特征之间(因为数据中不应该有关于这些协方差的信息)。在 MCMCglmm 中,可以使用如下代码完成此操作(假设两个女性特征位于响应变量矩阵的第 1 列和第 2 列,两个男性特征位于第 3 列和第 4 列):

rcov = ~us(at.level(trait, 1:2)):units + us(at.level(trait, 3:4)):units

但在 sommer 中似乎没有等效函数:我收到错误消息

Error: On the meantime the only rcov structures available are:
     'rcov=~units' or 'rcov=~at(.):units'.

我接下来尝试为每个人分配一个数字,并像这样拟合个人级别的随机效应:

library(sommer)

# Generate some fake data: 
# 100 males and 100 females
# Two traits are measured on each male, and two traits on each female
# 20 individuals per sex are measured for each of 5 different genotypes 
df <- data.frame(
  sex = rep(c("female", "male"), each = 100),
  female_trait_1 = c(rnorm(100), rep(NA, 100)),
  female_trait_2 = c(rnorm(100), rep(NA, 100)),
  male_trait_1 = c(rep(NA, 100), rnorm(100)),
  male_trait_2 = c(rep(NA, 100), rnorm(100)),
  genotype = rep(rep(1:5, each = 20), 2),
  individual = 1:200
)
df$genotype <- as.factor(df$genotype)
df$individual <- as.factor(df$individual)

sommer_test <- mmer2(

  # four traits as multivariate response
  cbind(female_trait_1, 
        female_trait_2,
        male_trait_1, 
        male_trait_2) ~ 1,

  # Fit the random effect of genotype (to estimate genetic covariance within and between sexes)
  # Try to fit US covariance matrices to specific levels of 'trait' (does not work)
  random =~ 
    us(trait):genotype +
    us(at.levels(trait, c("female_trait_1", "female_trait_2"))):individual +
    us(at.levels(trait, c("male_trait_1", "male_trait_2"))):individual,

  data = df
)
summary(sommer_test)

然而,后者也不起作用 - 它运行,但美国个人矩阵只适合两次特征的所有组合(包括我试图避免使用的男性 - 女性组合 at.levels).所以,似乎 at.levels 的工作方式与 MCMCglmm 中的 at.level 不同,因为我在这里使用它似乎什么都不做。

我是不是漏掉了一些关键的东西,或者 sommer 目前缺少这个功能?

非常感谢!

为了在 sommer 中处理不同的性别,你需要安装 sommer >=3.7。该模型直接操作 Gtc 参数(在 vs 函数中),该参数处理应用于方差-协方差分量的约束。值 1、2、3 对应于 positive、unconstrained、fixed。

鉴于此,数据的结构清楚地表明所有性状之间的协方差只能针对 "genotype" 随机效应进行估计:

> unsm(4)
     [,1] [,2] [,3] [,4]
[1,]    1    2    2    2
[2,]    0    1    2    2
[3,]    0    0    1    2
[4,]    0    0    0    1

而 "individual" 不是一个完全非结构化的特征,而是看起来像这样:

> mm <- adiag1(unsm(2),unsm(2));mm
     [,1] [,2] [,3] [,4]
[1,]    1    2    0    0
[2,]    0    1    0    0
[3,]    0    0    1    2
[4,]    0    0    0    1

长话短说模型如下:

# Generate some fake data: 
# 100 males and 100 females
# Two traits are measured on each male, and two traits on each female
# 20 individuals per sex are measured for each of 5 different genotypes 
set.seed(3434)
df <- data.frame(
  sex = rep(c("female", "male"), each = 100),
  female_trait_1 = c(rnorm(100), rep(NA, 100)),
  female_trait_2 = c(rnorm(100), rep(NA, 100)),
  male_trait_1 = c(rep(NA, 100), rnorm(100)),
  male_trait_2 = c(rep(NA, 100), rnorm(100)),
  genotype = rep(rep(1:5, each = 20), 2),
  individual = 1:200
)
df$genotype <- as.factor(df$genotype)
df$individual <- as.factor(df$individual)

library(sommer)
mm <- adiag1(unsm(2),unsm(2));mm
mix <- mmer(cbind(female_trait_1, 
                  female_trait_2,
                  male_trait_1,
                  male_trait_2) ~ 1,
     random=~vs(genotype,Gtc=unsm(4)) + vs(individual,Gtc=mm),
     rcov=~vs(units), na.method.Y = "include",
     data=df)

mix$sigma
cov2cor(mix$sigma$genotype)
cov2cor(mix$sigma$individual)

您可能希望使用真实数据来获得一些有意义的结果,但这说明了方法。干杯。