Efeitos genéticos e ambientais no tamanho de esquilos

Imagine que um estudo com esquilos tenha mostrado que animais maiores possuem maior aptidão (fitness) — eles têm, em média, maior probabilidade de sobreviver e/ou se reproduzir. Assim, o tamanho dos esquilos está sob seleção positiva, mas será que ele irá evoluir e, em caso afirmativo, quão rápido? A resposta depende da medida em que a variação de tamanho provém de efeitos genéticos vs. efeitos ambientais (por exemplo, variação na capacidade de forrageamento).

Imagine que temos uma amostra de 500 esquilos capturados no campo, medidos e com amostras de DNA coletadas. As relações familiares (ou seja, o pedigree) foram determinadas por meio da genotipagem de um conjunto de marcadores moleculares. O pedigree abrange duas gerações (100 pais, 400 descendentes). Medidas de qualidade do forrageamento no local de coleta de cada indivíduo também foram feitas.

squirels<-read.table("https://raw.githubusercontent.com/diogro/evofencom/refs/heads/main/Tutoriais/squirels2.txt", header=T)
squirels$ID<-as.factor(squirels$ID)
squirels$sex<-as.factor(squirels$sex)
squirels$SIRE<-as.factor(squirels$SIRE)
squirels$DAM<-as.factor(squirels$DAM)

#read in the pedigree file
squirelsPED<-read.table("https://raw.githubusercontent.com/diogro/evofencom/refs/heads/main/Tutoriais/squirelsPED.txt",header=T)
#2. Automatically installs packages in this list that user does not already have

list_pkgs <- c("lme4","lmerTest","MCMCglmm","pedigreemm", "sommer")
new_pkgs <- list_pkgs[!(list_pkgs %in% installed.packages()[,"Package"])]
if(length(new_pkgs) > 0){ install.packages(new_pkgs) }

library(sommer)
library(MCMCglmm)
library(lme4)
library(lmerTest)
library(pedigreemm)

Estrutura de dados

Vamos dar uma olhada na estrutuda dos dados:

head(squirels)
##      ID aggression     size sex SIRE  DAM      forage
## 401 401   3.675307 23.25477   1 <NA> <NA> -0.14936773
## 402 402   9.898771 17.23754   1 <NA> <NA> -1.14433949
## 403 403   9.549621 20.03212   1 <NA> <NA>  0.45325762
## 404 404   5.139761 24.35055   1 <NA> <NA>  0.12827954
## 405 405   3.819785 21.54173   1 <NA> <NA>  0.03174167
## 406 406   3.368013 18.36796   1 <NA> <NA>  0.54561204
head(squirelsPED)
##    ID DAM SIRE
## 1 401  NA   NA
## 2 402  NA   NA
## 3 403  NA   NA
## 4 404  NA   NA
## 5 405  NA   NA
## 6 406  NA   NA

Espero que isso seja autoexplicativo! Você pode incluir informações sobre o pai (SIRE) e a mãe (DAM) no arquivo de dados principal (como foi feito aqui), embora, por diversos motivos, seja frequentemente mais conveniente ter também um arquivo de dados de pedigree separado ao executar modelos animais (animal models).

Note que o arquivo de pedigree separado é organizado de modo que cada esquilo deve aparecer na coluna “indivíduo” (individual) antes de poder aparecer como genitor. Na prática, isso significa que a geração mais antiga estará no topo do arquivo e os pais desses indivíduos serão desconhecidos, portanto, NA.

Também podemos representar o pedigree graficamente, por exemplo, utilizando o pacote Pedantics do R.

G¨RAFICOS

Estimando a \(V_A\) para o tamanho corporal em esquilos

Dado que os dados que utilizaremos foram, na verdade, simulados a partir de um delineamento simples de irmãos-germanos (full-sib) e meio-irmãos (half-sib) para a segunda geração (em oposição a uma estrutura de pedigree multigeracional complexa), e como também temos os fenótipos dos pais, temos opções de como analisar os dados; portanto, tentaremos algumas abordagens diferentes.

ANOVA simples de dados de irmãos (sibship)

A regressão pai-filho poderia ser feita aqui, mas vamos ignorá-la e tentar a ANOVA como um ponto de partida simples. Isso é apropriado se utilizarmos apenas dados fenotípicos da geração de descendentes. Para esta geração, os dados são balanceados (número igual de descendentes por mãe e de mães por pai), sem fenótipos ausentes.

Note que, de forma confusa, o comando anova no R é usado genericamente para comparar modelos ou testar a significância de termos, enquanto o R utiliza o comando aov para realizar uma ANOVA de mínimos quadrados ordinários.

#6. Run an ANOVA
mod_aov<-aov(size~SIRE, data=squirels, na.action=na.omit) 
                      # this will drop parental generation for us since they have missing SIRE
summary(mod_aov)  
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## SIRE         19  546.5  28.761   6.686 3.43e-15 ***
## Residuals   380 1634.5   4.301                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 100 observations deleted due to missingness

O efeito do pai (sire) é estatisticamente significativo, o que significa que existe variância genética. Poderíamos executar uma ANOVA de dois fatores definindo a fórmula do modelo como tamanho ~ SIRE + DAM, se quiséssemos (tente fazer isso!). Também poderíamos usar esse modelo estatístico para calcular a herdabilidade. No entanto, o resultado é fornecido em termos de soma de quadrados (Sum Sq), e não como variância na escala real dos dados; portanto, para evitar ter que consultar muitas fórmulas e fazer cálculos aritméticos, vamos mudar para uma estrutura de modelos mistos, que é mais simples de interpretar.

Modelo sire-dam ajustado por REML

Primeiro, usaremos o lme4 como um pacote genérico para ajustar modelos mistos simples (ou seja, modelos em que os fatores aleatórios possuem observações repetidas) por REML. Assim, podemos ajustar o pai (sire) como um efeito aleatório sobre o tamanho, consultar o summary do modelo para ver a estimativa de \(V_{sire}\) e testar sua significância usando um teste de razão de verossimilhança (likelihood ratio test, LRT).

#8. Mixed model with SIRE as random
squirels_nosires = na.omit(squirels) # Omitindo observações sem fenótipos
mod_sire<-lmer(size~(1|SIRE), data=squirels_nosires) 

summary(mod_sire)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: size ~ (1 | SIRE)
##    Data: squirels_nosires
## 
## REML criterion at convergence: 1756.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7349 -0.7010  0.0448  0.6429  3.4329 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  SIRE     (Intercept) 1.223    1.106   
##  Residual             4.301    2.074   
## Number of obs: 400, groups:  SIRE, 20
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  21.3906     0.2681 19.0000   79.77   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rand(mod_sire)            #this will test significance of random effects using LRT
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## size ~ 1 + (1 | SIRE)
##            npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>        3 -878.26 1762.5                         
## (1 | SIRE)    2 -908.02 1820.0 59.513  1  1.215e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Expandindo o modelo

No entanto, irmãos completos também compartilham a mesma mãe, portanto, provavelmente é melhor modelar essa não independência. Ignorar esse fator poderia, potencialmente, levar a uma estimativa inflada da variância do pai (SIRE).

Lembrando que este é um estudo de campo no qual os fatores ambientais não puderam ser controlados experimentalmente, podemos querer controlar estatisticamente, também, quaisquer efeitos da variação na qualidade da forragemento sobre o tamanho corporal. Finalmente, podemos desejar controlar qualquer dimorfismo sexual no tamanho do corpo.

É claro que sempre vale a pena examinar os dados para verificar se os efeitos da forragemento ou do sexo são aparentes!

#9. Look at the data
plot(squirels$forage, squirels$size, pch = 19)

plot(squirels$sex, squirels$size)  #sex 1 is female, 2 is male

Nenhuma relação muito clara entre essas duas variáveis e tamanho corporal, mas vamos incluí-las como efeitos fixos e testar formalmente a relação entre essas variáveis preditores e a variável resposta. Vamos também incluir um efeito randômico adicional relacionado às mães. Nós sempre podemos simplificar o modelo depois.

#10. A SIRE-DAM MODEL 
squirels_nosires = na.omit(squirels) #  omitindo observações sem fenótipo. 
mod_sire_dam<-lmer(size~sex+forage+ (1|SIRE)+(1|DAM), data=squirels_nosires) 
rand(mod_sire_dam)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## size ~ sex + forage + (1 | SIRE) + (1 | DAM)
##            npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>        6 -861.68 1735.3                         
## (1 | SIRE)    5 -868.92 1747.8 14.491  1  0.0001408 ***
## (1 | DAM)     5 -870.76 1751.5 18.172  1  2.018e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

O summary do modelo revela que as variâncias de pais e mães são parecidas, e que os efeitos fixos são pequenos, porém significativos.

summary(mod_sire_dam)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: size ~ sex + forage + (1 | SIRE) + (1 | DAM)
##    Data: squirels_nosires
## 
## REML criterion at convergence: 1723.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.60525 -0.64051  0.00728  0.65391  3.08974 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  DAM      (Intercept) 0.8224   0.9068  
##  SIRE     (Intercept) 1.0155   1.0077  
##  Residual             3.4693   1.8626  
## Number of obs: 400, groups:  DAM, 80; SIRE, 20
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  21.0116     0.2807  24.2002  74.864  < 2e-16 ***
## sex2          0.7482     0.1897 331.1339   3.943 9.81e-05 ***
## forage        0.2883     0.1032 361.7645   2.795  0.00547 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr) sex2  
## sex2   -0.339       
## forage -0.048  0.123

Às vezes, usando um modelo desse tipo, chamado SIRE-DAM Model, você pode achar que a variância de DAMS é muito maior que a variância de SIRE. Isso indica contribuição importante de efeitos de dominância ou de efeitos maternos, ambos contribuindo para a variância do fenótipo.

Calculando a herdabilidade \(h^2\)

Nesse desenho experimental, a variância associada aos pais, a SIRE variance, é um quarto da variância genética aditiva (\(v_{sire} = 0.25 \times V_A\), você consegue imaginar por que?). A variância total é a soma de todas as variâncias. Nós podemos usar esses dois fatos pra calcular uma estimativa da herdabilidade.

#11. Calculate h2
h2_sire_dam=(4*1.0155)/(0.8224+1.0155+3.4693)
h2_sire_dam |> round (2)
## [1] 0.77

Modelos animais ajustados usando REML

Opções de software

Infelizmente, não são muitos os pacotes genéricos de modelos mistos REML que conseguem ajustar modelos animais (animal models). Uma boa opção é o ASReml-R, que é um pacote de modelos mistos muito potente e flexível que permite ajustar modelos animais básicos e estendê-los de diversas formas interessantes para aplicações em genética quantitativa. No entanto, ele não é gratuito e não temos acesso a ele hoje. O WOMBAT é um bom programa REML gratuito (e um pouco complicado), mas precisa ser executado fora do R.

Portanto, por enquanto, usaremos um pacote chamado pedigreemm, que faz interface com o lme4. Esses pacotes não são os ideais para ajustar modelos animais, mas, para fins ilustrativos, podemos usá-los para estimar a \(h^2\) (herdabilidade) em um cenário simples. Mais tarde, ajustaremos o mesmo modelo animal usando o MCMCglmm. Este é um pacote de modelos mistos para R, potente e gratuito, que consegue ajustar modelos animais. No entanto, o MCMCglmm exige que adotemos um pensamento bayesiano (o que tem tanto vantagens quanto desvantagens) e verificações de modelo pós-ajuste mais extensas para garantir que “acreditamos” nos resultados.

Ajustando o modelo no pedigreemm

O pedigreemm precisa que o objeto do pedigree seja separado nos dados fenotípicos, num formato levemente diferente do que nós temos usado. Depois de ajustar o pedigree, podemos ajustar o modelo animal, mantendo os mesmos efeitos fixos que usamos anteriormente, mas substituindo os efeitos aleatórios de SIRE e DAM com um efeito associado ao ID de cada indivíduo e ponderado pelo pedigree.

#13.Animal model in pedigreemm
squirelsPED2<-pedigree(squirelsPED$SIRE,squirelsPED$DAM,squirelsPED$ID)  #restructure ped file

mod_animalREML<-pedigreemm(size~sex+forage+(1|ID),pedigree=list(ID=squirelsPED2), 
                   data=squirels, REML=TRUE, 
                   control=lmerControl(check.nobs.vs.nlev="ignore",
                                       check.nobs.vs.nRE="ignore"))
## 'as(<dtTMatrix>, "dtCMatrix")' is deprecated.
## Use 'as(., "CsparseMatrix")' instead.
## See help("Deprecated") and help("Matrix-deprecated").
summary(mod_animalREML)
## Linear mixed model fit by REML ['lmerpedigreemm']
## Formula: size ~ sex + forage + (1 | ID)
##    Data: squirels
## Control: 
## lmerControl(check.nobs.vs.nlev = "ignore", check.nobs.vs.nRE = "ignore")
## 
## REML criterion at convergence: 2150.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.29066 -0.54643 -0.01826  0.56947  2.75949 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 2.252    1.501   
##  Residual             2.615    1.617   
## Number of obs: 500, groups:  ID, 500
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 20.93889    0.18834 111.174
## sex2         0.80969    0.18230   4.441
## forage       0.28035    0.09156   3.062
## 
## Correlation of Fixed Effects:
##        (Intr) sex2  
## sex2   -0.289       
## forage -0.014  0.088

A variância associada ao ID agora é uma estimativa direta da variância genética aditiva, e podemos calcular a herdabilidade simplesmente dividindo a variância aditiva pela variância total.

#15. And calculate h2
h2_animalREML=2.252/(2.252+2.615)
h2_animalREML
## [1] 0.462708

Por que as duas estimativas de herdabilidade são tão diferentes?

O modelo sire-dam (pai-mãe) produz uma estimativa de \(h^2 = 0,765\), que é notavelmente superior à estimativa do modelo animal de \(h^2 = 0,463\). Por que isso ocorre? Na verdade, ambas são estimativas válidas e nenhuma das abordagens é inerentemente enviesada neste caso. Como os dados foram simulados, eu sei que a estimativa mais baixa está muito mais próxima da “verdade” (ou seja, a \(h^2\) que foi definida para a população hipotética da qual os pais foram amostrados). É claro que, em situações reais, a “verdade” é desconhecida — mas espera-se que a estimativa do modelo animal seja melhor no sentido de ser mais precisa porque:

  • Embora o tamanho da amostra seja baixo para um estudo de genética quantitativa, o modelo sire-dam é incapaz de utilizar todos os dados fenotípicos disponíveis aqui. Ele não consegue usar os dados fenotípicos dos pais.

  • O modelo sire-dam não utiliza os dados de forma eficiente. A estimativa de \(V_A\) é, na verdade, feita pela comparação das médias fenotípicas entre meio-irmãos paternos, das quais existem apenas 20. Ao usar um modelo sire-dam (e assumindo apenas herança aditiva), o número de pais — em vez de mães por pai ou descendentes por mãe — é a principal limitação para a precisão.

  • O modelo animal utiliza todos os dados com máxima eficiência. A estimativa de \(V_A\) é informada por todas as comparações pareadas possíveis entre os indivíduos. Isso inclui pares de pais-prole, irmãos completos, meio-irmãos e indivíduos não aparentados.

Caracterizando a incerteza

A combinação pedigreemm/lme4 utilizada aqui não gera prontamente medidas de incerteza (ex: erro padrão ou intervalo de confiança) para \(V_A\) ou \(h^2\) em nenhum dos modelos. Poderíamos provavelmente criar um método (ex: bootstrapping), mas a codificação não é trivial. Quantificar a incerteza nos parâmetros de genética quantitativa é crucial, e essa limitação é um bom motivo para trabalhar com pacotes mais especializados. Por exemplo, usando o ASReml-R (código não mostrado), as estimativas de \(V_A\) e \(h^2\) são (tranquilizadoramente!) as mesmas obtidas acima. No entanto, este programa também estima o erro padrão (EP) que, embora aproximado por razões técnicas, fornece uma imagem clara da precisão relativa das duas estratégias de modelagem. Também é fácil testar a \(V_A\) (e quaisquer outros componentes de variância) usando testes de razão de verossimilhança no ASReml.

  • Modelo sire-dam (ASReml): \(V_{pai} = 1,02\) \((0,457)\); \(V_A = 4,06\) \((1,83)\); \(h^2 = 0,765\) \((0,290)\)
  • Modelo animal (ASReml): \(V_A = 2,25\) \((0,480)\); \(h^2 = 0,463\) \((0,078)\)