Esta atividade foi proposta na disciplina de Biometria de Marcadores Genéticos da ESALQ, ministrada pelo professor Antonio Augusto Garcia, mais informações sobre a disciplina e outras atividades referentes à ela podem ser encontradas neste link.

Descrição da atividade proposta:

Considere novamente o conjunto de dados de Mimulus e o caráter fenotípico que você escolheu em exercícios anteriores.

Carregando os dados

Neste primeiro passo vamos carregar os dados genotípicos e fenotípicos contidos no arquivo m_feb06.raw e as informações do mapa genético construído para tais dados contidos em mimulus.map. O mapa genético foi construído uma atividade anterior da mesma disciplina, que pode ser acessada neste link. Nele também podem ser encontradas descrições mais detalhadas sobre o conjunto de dados.

#Carregando os pacotes onemap e qtl:
library(onemap)
library(qtl)
library(knitr)

#Realizando a leitura dos dados do experimento de QTL:
mimulus.qtl <- read.cross(format = "mm", file = "m_feb06.raw", 
                        mapfile = "mimulus.map")
##  --Read the following data:
##  Type of cross:          f2 
##  Number of individuals:  287 
##  Number of markers:      418 
##  Number of phenotypes:   16 
##  --Cross type: f2
#Observando o mapa genético:
plot.map(mimulus.qtl)

Também numa atividade prévia, que pode ser acessada neste link, exploramos uns dos fenótipos contidos no conjunto de dados, referente ao número de polens não viáveis. Foi realizado o mapeamento por intervalos e mapeamento por intervalo composto.

Aqui será realizado o mapeamento por múltiplos intervalos para o mesmo fenótipo, que possui a seguinte distribuição, a qual consideraremos próxima à normal:

plot.pheno(mimulus.qtl, 5)

Mapeamento por Intervalo Composto

Aqui consideraremos o perfil de LOD estabelecido com mapeamento por intervalo composto para definir os QTLs num primeiro modelo para o mapeamento por intervalo múltiplo. Portanto, segue uma breve descrição da análise de mapeamento por intervalo composto já realizada em link.

# Carregando análises do mapeamento por intervalo composto
load("out.cim.20.RData")
load("permic_5.rdata")

# Estabelecendo limiar
lim <- summary(permic_5, alpha = 0.05)

# Gráficos de perfil de lod
plot(out.cim.20,  col= "red")
add.cim.covar(out.cim.20, col="green")
abline(h=lim[1], col = "orange")

plot(out.cim.20,  col= "red", chr = c(5,13))
add.cim.covar(out.cim.20, col="green")
abline(h=lim[1], col = "orange")

# Intervalo de confiança considerando a regra de 1 LOD
intervalo13 <- lodint(out.cim.20, chr=13, drop=1)
intervalo6 <- lodint(out.cim.20, chr=5, drop=1)
intervalo <- rbind(intervalo6, intervalo13)
knitr::kable(intervalo)
chr pos lod
BB186 5 0 4.106504
BB186 5 0 4.106504
c5.loc17 5 17 3.015203
c13.loc146 13 146 1.662089
c13.loc149 13 149 12.893522
c13.loc153 13 153 11.807917 
# Efeitos do QTL
QTL <- summary(out.cim.20, threshold = lim[1])
qtls <- makeqtl(mimulus.sim, chr = QTL$chr, pos = QTL$pos)
qtls
##   QTL object containing imputed genotypes, with 16 imputations. 
## 
##        name chr    pos n.gen
## Q1    5@0.0   5   0.00     3
## Q2 13@153.9  13 153.86     3
model <- fitqtl(mimulus.sim, pheno.col = 5,  qtl = qtls, get.ests = TRUE)
resumo <- summary(model)

CIM_QTLs <- data.frame("marcador"= rownames(QTL),"cromossomo"= QTL$chr,"posição" = 
               QTL$pos, "LOD" = QTL$lod, 
             "efeito aditivo" = c(resumo$ests[2,1],resumo$ests[4,1]) , 
             "efeito de dominância"= c(resumo$ests[3,1],resumo$ests[5,1]),
             "R^2" = resumo$result.drop[,4])

knitr::kable(CIM_QTLs)
marcador cromossomo posição LOD efeito.aditivo efeito.de.dominância R.2
5@0.0 BB186 5 0 4.106504 -24.89030 11.96873 5.325632
13@153.9 c13.loc149 13 149 12.893522 -46.77116 -25.84217 18.879967

No intuito de sermos mais abrangentes, vamos considerar um limiar mais baixo, de valor 2.5.

plot(out.cim.20,  col= "red")
add.cim.covar(out.cim.20, col="green")
abline(h=2.5, col = "orange")

Avaliando dessa forma temos a sugestão de 5 QTLs para considerarmos num primeiro modelo a ser testado para o mapeamento por intervalo múltiplo.

Escaneamento bi-dimensional

Antes de montarmos os modelos, podemos também ter uma visão inicial dos efeitos e da interação dos QTLs com a função scantwo, que realiza uma análise bi-dimensional dos QTLs, ou seja, são considerados dois QTLs simultaneamente e a interação entre eles.

mimulus.scan2 <- scantwo(mimulus.qtl.prob, pheno.col=5, method="hk", verbose=FALSE)
save(mimulus.scan2, file = "mimulus.scan2.RData")
plot(mimulus.scan2)

O gráfico mostra valores de \(LOD_i\) no triângulo superior esquerdo e valores de \(LOD_f\) é mostrado no triângulo inferior direito. Na escala de cores à direita, o número à esquerda e à direita correspondem aos valores de \(LOD_i\) e \(LOD_f\), respectivamente.

O \(LOD_f\) mede a melhora no ajuste de um modelo considerando dois locos em relação ao modelo de hipótese nula, o que indica evidências para efeitos dos QTLs e da interação entre eles.

Já o \(LOD_i\) mede a melhora do ajuste de um modelo considerando dois locos em relação ao modelo aditivo, o qual não inclui as interações entre os QTLs. Portanto o \(LOD_i\) destaca efeitos das possíveis interações, uma vez que o modelo aditivo funciona como uma hipótese \(H_0\).

Portanto, o gráfico mostra evidências de interação do QTL no cromossomo 13 com QTLs nos cromossomos 3, 5, 6 e 14.

plot(mimulus.scan2, chr = c(3,5,6,13, 14))

Com o argumento lower podemos alterar o modelo utilizao para traçar perfil de LOD do triângulo inferiro do gráfico. Aqui alteramos para o modelo aditivo para visualizar os efeitos aditivos dos QTLs, sem considerar as interações.

plot(mimulus.scan2, lower="add")

Pelo gráfico gerado podemos verificar no triângulo inferior do gráfico, uma fraca evidência (em relação à demonstrada no cromossomo 13) de efeito aditivo no final do cromossomo 3 e no início do 5, assim como apontado no gráfico de perfil de LOD do mapeamento por intervalo composto.

Além da obsevação dos gráficos, também podemos definir limiares para cada um dos modelos testados automaticamente pela função scantwoe verificar a significância dos efeitos. Uma forma robusta de obter esses limiares é com permutações utilizando o argumento n.permna função scantwo.

load("out.cim.20.RData")
set.seed(85842518)
operm2a <- scantwo(mimulus.qtl.prob, n.perm=200, pheno.col = 5, method = "hk")
save(operm2a, file="perm2a.RData")
set.seed(85842519)
operm2b <- scantwo(mimulus.qtl.prob, n.perm=200, pheno.col = 5, method = "hk")
save(operm2b, file="perm2b.RData")
set.seed(85842520)
operm2c <- scantwo(mimulus.qtl.prob, n.perm=200, pheno.col = 5, method = "hk")
save(operm2c, file="perm2c.RData")
set.seed(85842521)
operm2d <- scantwo(mimulus.qtl.prob, n.perm=200, pheno.col = 5, method = "hk")
save(operm2d, file="perm2d.RData")
set.seed(85842522)
operm2e <- scantwo(mimulus.qtl.prob, n.perm=200, pheno.col = 5, method = "hk")
save(operm2e, file="perm2e.RData")
operm2 <- c(operm2a, operm2b,operm2c, operm2d, operm2e)
head(operm2)
## $full
##             nv
## [1,]  7.752810
## [2,] 10.003507
## [3,]  8.746593
## [4,]  7.826855
## [5,]  7.714861
## [6,]  8.077097
## 
## $fv1
##            nv
## [1,] 7.089636
## [2,] 8.986823
## [3,] 7.325358
## [4,] 6.594705
## [5,] 5.804937
## [6,] 7.344686
## 
## $int
##            nv
## [1,] 6.298393
## [2,] 8.346437
## [3,] 5.286384
## [4,] 5.892345
## [5,] 4.635325
## [6,] 5.762379
## 
## $add
##            nv
## [1,] 4.344257
## [2,] 4.735960
## [3,] 6.089052
## [4,] 3.363841
## [5,] 5.069428
## [6,] 4.564338
## 
## $av1
##            nv
## [1,] 2.499270
## [2,] 2.386185
## [3,] 4.592953
## [4,] 1.969422
## [5,] 3.263803
## [6,] 3.007831
## 
## $one
##            nv
## [1,] 2.393313
## [2,] 2.349775
## [3,] 2.567883
## [4,] 1.721342
## [5,] 2.392513
## [6,] 2.238330
## 
## attr(,"class")
## [1] "scantwoperm" "list"       
## attr(,"method")
## [1] "hk"
summary(operm2)
## nv (1000 permutations)
##     full  fv1  int  add  av1  one
## 5%  11.8 9.97 8.49 7.23 5.66 3.70
## 10% 11.2 9.39 7.92 6.65 4.86 3.37
summary(mimulus.scan2, perms=operm2, alphas=c(0.4, 0.4, 0, 0.4, 0.4),
          pvalues=TRUE)
##        pos1f pos2f lod.full pval lod.fv1  pval lod.int  pval     pos1a
## c3:c13   134   139     20.3    0    8.63 0.233   5.353 0.762       134
## c5:c13     6   150     16.6    0    4.92 0.998   0.633 1.000         5
## c6:c13    69   151     18.2    0    6.56 0.810   3.341 1.000        32
##        pos2a lod.add pval lod.av1  pval
## c3:c13   150    15.0    0    3.27 0.285
## c5:c13   151    16.0    0    4.28 0.138
## c6:c13   151    14.9    0    3.22 0.301

Assim, evidencia-se também interações com o QTL no cromossomo 14. Podemos observar com mais detalhe os efeitos das interações utilizando as funções plot.pxg e effectplot.

# effectplot exige dados imputados, portanto utilizamos o mimulus.sim

mar <- find.marker(mimulus.sim, c("6","13"), c(32,151))

plot.pxg(mimulus.sim, marker=mar, pheno.col = 5)

effectplot(mimulus.sim, mname1="6@32", mname2="13@151", pheno.col = 5)

Para esses dois marcadores existem poucos dados imputados (em vermelho) mostrados no primeiro gráfico, pois trata-se de dois marcadores codominantes.

$ grep -A 1 "MgSTS504A" m_feb06.raw 
*MgSTS504A 
HHHAHBHAHHABHHBH-AHHHAHBBBABBHHAHBHHBHHHBHHHAHAHBAHHHHABAA-HBBBHHBA-AH-B-HHBBHB-HHBBAHABHHBABHBBBABA-HBA-BBBBBHBHBBHHHBHHHBHH--HA-HBAAAHHBHHHBH-B-AAHBBAHHAHHABAHHHBAHHHBHHA-HA-HHHHHABBBBAHABHBHAABHBHAAHAA-AHHHHHHH-HA-HAHHHBBHHAHBHAHABBBABAAHBBAHHHBBHAHHABHHABHBHHBHHHHHAHHBHHHBABBBBAHHHB
$ grep -A 1 "MgSTS316" m_feb06.raw 
*MgSTS316 
B--ABHAAH-AHHHBH-AHAHAAH-BHHBBAHHB-AHHHBHHBAHBHAAABH-ABH-BHHH------------AH---B-AHHHHHHHAAB-BHBHHHB-BHBHAHHHHH-BB-H--BHAHA-BH-B-B-A-HABAAHH-A---BBHHA--HHBHHHBHAHBBH-H-HHBBH-HH-A--H-HHBABBAHBHHH-A--HHAABAHHBAAHHBHAHBAAAHHHH-HA-BBAHAHAHBHHBAHHBHHBHHBHHHHAA-BAHB--HHHBBAHBHBHBHH-BBBHHBHHAHB

O segundo gráfico evidencia a interação entre eles já que as linhas se cruzam.

mar <- find.marker(mimulus.sim, c("3","13"), c(134,139))

plot.pxg(mimulus.sim, marker=mar, pheno.col = 5)

effectplot(mimulus.sim, mname1="3@134", mname2="13@139", pheno.col = 5)

Neste caso o marcador AG19 localizado no cromossomo 3 é codominante e o marcador CB55 localizado no 13 é dominante, por isso uma parte maior dos fenótipos são imputados (em vermelho). Portanto, aqui podemos destacar a importância de termos marcadores codominantes no nosso conjunto de dados, pois eles auxiliam na imputação.

$ grep -A 1 "AG19" m_feb06.raw 
*AG19 
AHB-HHHAHHHHHBA-AHHAHBHHBHHAHHHHHBABHBBHABHAAAHHHHHHAHHHHHHHAAHHAHHHHHBHB-AAAHBHHHHHBBAHBA-HHHAAAAAHBHHABHHAAHAAHBHHBAABHH-AHHHHHAHHHH-ABAABH---A--HAHHAHHHHAB--AA-HHAA-BHAH-HAH-AHABHBHAAHAH-B-ABBA-AAB-HHABAABAABABHAHAAA-BBAA--AHABBHHB-HHHABHABHHBH-HHHHHBHAHAH-HBHHHHHHHHHHBAA-HHBHHHHABAB
$ grep -A 1 "CB55" m_feb06.raw 
*CB55 
BDBDBDDDDBDDD-DDDB-BDBDBBBDDDBDBDBDBBBDDBDDDDDDDBDDBBDBDBDDDBDDDDBB-BBDBBDBDDBBDBDDDDDDDD-DDBBDDBBDDDDDDBBBDBDDDBDDDDDBBDDDDBDBDDDDD-DDBBDBDBD--BBDBBDDDDDDDDDBDDDDDBDDDBDDDDDDDDDBDBDDDDDDDDDDB-DBBB-DDBBDDDDDBDDBDDDDDDDDDBDDBDBBBB-BDDDDDBDBBDDDBDDDBDDDDBDDDBBDDBDDBBDBBBDDDBDBBDBDBDDDDDDB

Assim como visto no caso anterior, o segundo gráfico também demonstra interação entre os locos.

Com tais análises podemos sugerir QTLs e interações entre eles para serem considerados no modelo de mapeamento múltiplo. Uma forma automática de obter sugestões seria utilizando a função stepwiseqtl.

outsw2 <- stepwiseqtl(mimulus.qtl.prob, max.qtl=8, pheno.col = 5)

save(outsw2, file = "outsw2.RData")
outsw2
##   QTL object containing genotype probabilities. 
## 
##        name chr pos n.gen
## Q1  3@132.0   3 132     3
## Q2    5@5.0   5   5     3
## Q3 13@146.0  13 146     3
## 
##   Formula: y ~ Q1 + Q2 + Q3 + Q1:Q3 
## 
##   pLOD:  11.52

Portanto, obtemos sugestões de QTLs a serem considerados no modelo de mapeamento múltiplo com o mapeamento por intervalo composto, com as análises bi-dimensionais de QTLs e com a função stepwiseqtl.

Mapeamento múltiplo

A partir das avaliações consideraremos os seguintes QTLs e suas interações.

# Selecionando os QTLs

qtl <- makeqtl(mimulus.sim, chr=c(3,5,6,13,13,14),
               pos=c(134, 6, 69, 150, 139,213))

plot(qtl)

# Primeiro modelo

model1 <- fitqtl(mimulus.sim, qtl=qtl, formula=y~Q1+Q2+Q3+Q4+Q5+Q6+Q1*Q5+Q2*Q4+Q3*Q4+Q4*Q5+Q1*Q6, pheno.col = 5)
summary(model1)
## 
##      fitqtl summary
## 
## Method: multiple imputation 
## Model:  normal phenotype
## Number of observations : 275 
## 
## Full model result
## ----------------------------------  
## Model formula: y ~ Q1 + Q2 + Q3 + Q4 + Q5 + Q6 + Q1:Q5 + Q1:Q6 + Q2:Q4 + 
##                     Q3:Q4 + Q4:Q5 
## 
##        df        SS        MS      LOD     %var Pvalue(Chi2) Pvalue(F)
## Model  32  775251.5 24226.609 35.88588 45.17081            0         0
## Error 242  941014.9  3888.491                                         
## Total 274 1716266.4                                                   
## 
## 
## Drop one QTL at a time ANOVA table: 
## ----------------------------------  
##                   df Type III SS    LOD   %var F value Pvalue(Chi2)
## 3@134.0           10      116843 6.9892 6.8080  3.0048        0.000
## 5@0.0              6       96877 5.8514 5.6446  4.1523        0.000
## 6@68.8             6       95066 5.7471 5.5391  4.0747        0.000
## 13@153.9          14      146474 8.6389 8.5344  2.6906        0.000
## 13@137.8          10      140150 8.2906 8.1660  3.6042        0.000
## 14@211.5           6       12032 0.7587 0.7011  0.5157        0.745
## 3@134.0:13@137.8   4       44311 2.7477 2.5818  2.8489        0.013
## 3@134.0:14@211.5   4        4929 0.3120 0.2872  0.3169        0.838
## 5@0.0:13@153.9     4       11558 0.7290 0.6734  0.7431        0.500
## 6@68.8:13@153.9    4       70181 4.2953 4.0892  4.5121        0.001
## 13@153.9:13@137.8  4       26224 1.6414 1.5280  1.6860        0.109
##                   Pvalue(F)    
## 3@134.0            0.001349 ** 
## 5@0.0              0.000543 ***
## 6@68.8             0.000650 ***
## 13@153.9           0.001094 ** 
## 13@137.8           0.000175 ***
## 14@211.5           0.796188    
## 3@134.0:13@137.8   0.024603 *  
## 3@134.0:14@211.5   0.866528    
## 5@0.0:13@153.9     0.563463    
## 6@68.8:13@153.9    0.001560 ** 
## 13@153.9:13@137.8  0.153853    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Neste primeiro modelo os efeitos foram diferentes do esperado conforme as análises anteriores, supomos que poderia ser pela incerteza nas posições dos QTLs, para isso utilizamos a função refineqtl para melhor verificarmos os efeitos.

# refinando as posições dos QTLs
rqtl <- refineqtl(mimulus.sim, qtl=qtl,
                  formula=y~Q1+Q2+Q3+Q4+Q5+Q6+Q1*Q5+Q2*Q4+Q3*Q4+Q4*Q5+Q1*Q6, 
                  verbose = FALSE, pheno.col = 5)

# Verificando novamente os efeitos com as posições atualizadas
model1.up <- fitqtl(mimulus.sim, qtl=rqtl,
                    formula=y~Q1+Q2+Q3+Q4+Q5+Q6+Q1*Q5+Q2*Q4+Q3*Q4+Q4*Q5+Q1*Q6,
                    pheno.col = 5)
summary(model1.up)
## 
##      fitqtl summary
## 
## Method: multiple imputation 
## Model:  normal phenotype
## Number of observations : 275 
## 
## Full model result
## ----------------------------------  
## Model formula: y ~ Q1 + Q2 + Q3 + Q4 + Q5 + Q6 + Q1:Q5 + Q1:Q6 + Q2:Q4 + 
##                     Q3:Q4 + Q4:Q5 
## 
##        df        SS        MS      LOD     %var Pvalue(Chi2) Pvalue(F)
## Model  32  829548.1 25923.379 39.43487 48.33446            0         0
## Error 242  886718.3  3664.125                                         
## Total 274 1716266.4                                                   
## 
## 
## Drop one QTL at a time ANOVA table: 
## ----------------------------------  
##                   df Type III SS     LOD    %var F value Pvalue(Chi2)
## 3@134.0           10      190151 11.6019 11.0793   5.190        0.000
## 5@12.6             6       99777  6.3675  5.8136   4.538        0.000
## 6@103.3            6       80238  5.1729  4.6751   3.650        0.001
## 13@153.9          14      205127 12.4266 11.9519   3.999        0.000
## 13@137.8          10      266150 15.6742 15.5075   7.264        0.000
## 14@104.5           6       32550  2.1528  1.8965   1.481        0.128
## 3@134.0:13@137.8   4       94926  6.0731  5.5309   6.477        0.000
## 3@134.0:14@104.5   4       14965  0.9994  0.8720   1.021        0.331
## 5@12.6:13@153.9    4       16271  1.0859  0.9481   1.110        0.287
## 6@103.3:13@153.9   4       51940  3.3993  3.0264   3.544        0.004
## 13@153.9:13@137.8  4       81784  5.2683  4.7652   5.580        0.000
##                   Pvalue(F)    
## 3@134.0            6.81e-07 ***
## 5@12.6             0.000221 ***
## 6@103.3            0.001735 ** 
## 13@153.9           3.31e-06 ***
## 13@137.8           4.96e-10 ***
## 14@104.5           0.185386    
## 3@134.0:13@137.8   5.75e-05 ***
## 3@134.0:14@104.5   0.397027    
## 5@12.6:13@153.9    0.352308    
## 6@103.3:13@153.9   0.007850 ** 
## 13@153.9:13@137.8  0.000259 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plotLodProfile(rqtl, ylab="Profile LOD score")

Podemos observar que o QTL no cromossomo 14 e algumas das interações sugeridas não foram significativas ao nível de 5%, portanto podemos desconsiderá-los do modelo.

# Segundo modelo
qtl <- makeqtl(mimulus.sim, chr=c(3,5,6,13,13),
               pos=c(134, 6, 69, 150, 139))

model2 <- fitqtl(mimulus.sim, qtl=qtl,
                 formula=y~Q1+Q2+Q3+Q4+Q5+Q1*Q5+Q3*Q4+Q4*Q5,
                 pheno.col = 5)

# Atualizando posições
rqtl <- refineqtl(mimulus.sim, qtl=qtl,
                  formula=y~Q1+Q2+Q3+Q4+Q5+Q1*Q5+Q3*Q4+Q4*Q5, 
                  verbose = FALSE, pheno.col = 5)

# Verificando novamente os efeitos com as posições atualizadas
model2.up <- fitqtl(mimulus.sim, qtl=rqtl,
                    formula=y~Q1+Q2+Q3+Q4+Q5+Q1*Q5+Q3*Q4+Q4*Q5, pheno.col = 5)
summary(model2.up)
## 
##      fitqtl summary
## 
## Method: multiple imputation 
## Model:  normal phenotype
## Number of observations : 275 
## 
## Full model result
## ----------------------------------  
## Model formula: y ~ Q1 + Q2 + Q3 + Q4 + Q5 + Q1:Q5 + Q3:Q4 + Q4:Q5 
## 
##        df        SS        MS      LOD     %var Pvalue(Chi2) Pvalue(F)
## Model  22  778592.4 35390.565 36.09827 45.36548            0         0
## Error 252  937674.0  3720.929                                         
## Total 274 1716266.4                                                   
## 
## 
## Drop one QTL at a time ANOVA table: 
## ----------------------------------  
##                   df Type III SS    LOD   %var F value Pvalue(Chi2)
## 3@134.0            6      174715 10.203 10.180   7.826        0.000
## 5@12.6             2       89551  5.447  5.218  12.033        0.000
## 6@103.3            6       78667  4.811  4.584   3.524        0.001
## 13@153.9          10      177473 10.351 10.341   4.770        0.000
## 13@137.8          10      283582 15.779 16.523   7.621        0.000
## 3@134.0:13@137.8   4      112647  6.775  6.564   7.568        0.000
## 6@103.3:13@153.9   4       48751  3.027  2.841   3.275        0.007
## 13@153.9:13@137.8  4       81038  4.950  4.722   5.445        0.000
##                   Pvalue(F)    
## 3@134.0            9.75e-08 ***
## 5@12.6             1.02e-05 ***
## 6@103.3            0.002291 ** 
## 13@153.9           2.83e-06 ***
## 13@137.8           1.28e-10 ***
## 3@134.0:13@137.8   8.99e-06 ***
## 6@103.3:13@153.9   0.012174 *  
## 13@153.9:13@137.8  0.000321 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plotLodProfile(rqtl, ylab="Profile LOD score")

O QTL no cromossomo 6, apesar de significativo, tem um efeito bem menor comparado aos outros. Além disso, a função stepwiseqtl já não indicava sua inclusão no modelo. Por isso decidimos removê-lo.

# Terceiro modelo

qtl <- makeqtl(mimulus.sim, chr=c(3,5,13,13),
               pos=c(134, 6, 150, 139))

model3 <- fitqtl(mimulus.sim, qtl=qtl,
                 formula=y~Q1+Q2+Q3+Q4+Q1*Q4+Q2*Q3+Q3*Q4,
                 pheno.col = 5)

# Atualizando posições
rqtl <- refineqtl(mimulus.sim, qtl=qtl,
                  formula=y~Q1+Q2+Q3+Q4+Q1*Q4+Q2*Q3+Q3*Q4,
                  verbose = FALSE, pheno.col = 5)

# Verificando novamente os efeitos com as posições atualizadas
model3.up <- fitqtl(mimulus.sim, qtl=rqtl,
                    formula=y~Q1+Q2+Q3+Q4+Q1*Q4+Q2*Q3+Q3*Q4, pheno.col = 5)
summary(model3.up)
## 
##      fitqtl summary
## 
## Method: multiple imputation 
## Model:  normal phenotype
## Number of observations : 275 
## 
## Full model result
## ----------------------------------  
## Model formula: y ~ Q1 + Q2 + Q3 + Q4 + Q1:Q4 + Q2:Q3 + Q3:Q4 
## 
##        df        SS       MS      LOD     %var Pvalue(Chi2) Pvalue(F)
## Model  20  717850.9 35892.55 32.35009 41.82631            0         0
## Error 254  998415.5  3930.77                                         
## Total 274 1716266.4                                                  
## 
## 
## Drop one QTL at a time ANOVA table: 
## ----------------------------------  
##                   df Type III SS    LOD   %var F value Pvalue(Chi2)
## 3@134.0            6      204477 11.126 11.914   8.670        0.000
## 5@12.6             6      117051  6.620  6.820   4.963        0.000
## 13@153.9          10      128422  7.226  7.483   3.267        0.000
## 13@137.8          10      246424 13.173 14.358   6.269        0.000
## 3@134.0:13@137.8   4      132934  7.464  7.746   8.455        0.000
## 5@12.6:13@153.9    4       17926  1.063  1.044   1.140        0.298
## 13@153.9:13@137.8  4       49264  2.876  2.870   3.133        0.010
##                   Pvalue(F)    
## 3@134.0            1.36e-08 ***
## 5@12.6             7.98e-05 ***
## 13@153.9           0.000541 ***
## 13@137.8           1.40e-08 ***
## 3@134.0:13@137.8   2.03e-06 ***
## 5@12.6:13@153.9    0.338103    
## 13@153.9:13@137.8  0.015376 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plotLodProfile(rqtl, ylab="Profile LOD score")

Pelo perfil da curva e a ausência de marcadores na região entre os dois picos de QTLs do cromossomo 13, é possível que componham, na verdade, apenas um QTL. Portanto, removeremos um dos QTLs do cromossomo 13.

# Quarto modelo

qtl <- makeqtl(mimulus.sim, chr=c(3,5,13),
               pos=c(134, 6, 139))

model4 <- fitqtl(mimulus.sim, qtl=qtl,
                 formula=y~Q1+Q2+Q3+Q1*Q3+Q2*Q3,
                 pheno.col = 5)

# Atualizando posições
rqtl <- refineqtl(mimulus.sim, qtl=qtl,
                  formula=y~Q1+Q2+Q3+Q1*Q3+Q2*Q3,
                  verbose = FALSE, pheno.col = 5)

# Verificando novamente os efeitos com as posições atualizadas
model4.up <- fitqtl(mimulus.sim, qtl=rqtl,
                    formula=y~Q1+Q2+Q3+Q1*Q3+Q2*Q3, pheno.col = 5)
summary(model4.up)
## 
##      fitqtl summary
## 
## Method: multiple imputation 
## Model:  normal phenotype
## Number of observations : 275 
## 
## Full model result
## ----------------------------------  
## Model formula: y ~ Q1 + Q2 + Q3 + Q1:Q3 + Q2:Q3 
## 
##        df        SS        MS      LOD     %var Pvalue(Chi2) Pvalue(F)
## Model  14  595552.7 42539.476 25.44988 34.70048            0         0
## Error 260 1120713.8  4310.438                                         
## Total 274 1716266.4                                                   
## 
## 
## Drop one QTL at a time ANOVA table: 
## ----------------------------------  
##                  df Type III SS     LOD    %var F value Pvalue(Chi2)
## 3@134.0           6      252451 12.1313 14.7093  9.7612        0.000
## 5@12.6            6       91423  4.6828  5.3269  3.5350        0.001
## 13@137.8         10      455626 20.3713 26.5475 10.5703        0.000
## 3@134.0:13@137.8  4      157728  7.8631  9.1902  9.1480        0.000
## 5@12.6:13@137.8   4        6124  0.3254  0.3568  0.3552        0.827
##                  Pvalue(F)    
## 3@134.0           1.06e-09 ***
## 5@12.6             0.00221 ** 
## 13@137.8          5.33e-15 ***
## 3@134.0:13@137.8  6.27e-07 ***
## 5@12.6:13@137.8    0.84030    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plotLodProfile(rqtl, ylab="Profile LOD score")

Podemos observar que a interação do QTL do cromossomo 5 e do 13 passa a ser não significativa, mostrando um modelo semelhante ao sugerido pela função stepwiseqtl. Portanto, utilizaremos o modelo proposto por ela.

# Quinto modelo

qtl <- makeqtl(mimulus.sim, chr=c(3,5,13),
               pos=c(132, 5, 146))

model5 <- fitqtl(mimulus.sim, qtl=qtl,
                 formula=y~Q1+Q2+Q3+Q1*Q3,
                 pheno.col = 5)

# Atualizando posições
rqtl <- refineqtl(mimulus.sim, qtl=qtl,
                  formula=y~Q1+Q2+Q3+Q1*Q3,
                  verbose = FALSE, pheno.col = 5)

# Verificando novamente os efeitos com as posições atualizadas
model5.up <- fitqtl(mimulus.sim, qtl=rqtl,
                    formula=y~Q1+Q2+Q3+Q1*Q3, pheno.col = 5, get.ests = TRUE)
summary(model5.up)
## 
##      fitqtl summary
## 
## Method: multiple imputation 
## Model:  normal phenotype
## Number of observations : 275 
## 
## Full model result
## ----------------------------------  
## Model formula: y ~ Q1 + Q2 + Q3 + Q1:Q3 
## 
##        df        SS        MS      LOD     %var Pvalue(Chi2) Pvalue(F)
## Model  10  589428.5 58942.850 25.12445 34.34365            0         0
## Error 264 1126837.9  4268.325                                         
## Total 274 1716266.4                                                   
## 
## 
## Drop one QTL at a time ANOVA table: 
## ----------------------------------  
##                  df Type III SS    LOD  %var F value Pvalue(Chi2)
## 3@134.0           6      265436 12.631 15.47  10.365            0
## 5@12.6            2       85299  4.357  4.97   9.992            0
## 13@137.8          6      449502 20.046 26.19  17.552            0
## 3@134.0:13@137.8  4      167336  8.268  9.75   9.801            0
##                  Pvalue(F)    
## 3@134.0           2.58e-10 ***
## 5@12.6            6.56e-05 ***
## 13@137.8           < 2e-16 ***
## 3@134.0:13@137.8  2.09e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Estimated effects:
## -----------------
##                        est      SE      t
## Intercept          108.393   4.429 24.472
## 3@134.0a             2.849   6.548  0.435
## 3@134.0d           -37.711   8.764 -4.303
## 5@12.6a            -24.260   6.046 -4.012
## 5@12.6d             -2.372   8.230 -0.288
## 13@137.8a          -41.132   6.197 -6.637
## 13@137.8d           -2.071   9.209 -0.225
## 3@134.0a:13@137.8a -15.031   8.019 -1.874
## 3@134.0d:13@137.8a  52.723  11.582  4.552
## 3@134.0a:13@137.8d -24.757  13.409 -1.846
## 3@134.0d:13@137.8d  55.650  19.343  2.877
plot(qtl)

plotLodProfile(rqtl, ylab="Profile LOD score")

Consideramos o modelo 5 como o mais adequado, a partir dele calculamos os efeitos dos QTLs e interações (mostrados a cima).

Verificando outras possibilidades automaticamente.

# Busca automatica por mais QTL:
out.ap <- addpair(mimulus.sim, qtl=rqtl, formula=y~Q1+Q2+Q3+Q1*Q3,verbose=FALSE,
                  pheno.col = 5)
save(out.ap, file = "out.ap.RData")
summary(out.ap, perms=operm2, alphas=c(0.4, 0.4, 0, 0.4, 0.4),
          pvalues=TRUE)
##        pos1f pos2f lod.full  pval lod.fv1  pval lod.int  pval     pos1a
## c2:c13  19.0   154     8.52 0.713    5.41 0.987    1.53 1.000        19
## c6:c6   12.6    14    10.14 0.261    8.53 0.257    5.49 0.717       194
##        pos2a lod.add  pval lod.av1  pval
## c2:c13   154    6.98 0.068    3.88 0.180
## c6:c6    195    4.64 0.623    3.03 0.356

O QTL no cromossomo 6 já foi testado no modelo e se demonstrou pouco significativo e não foi incluído no modelo. O QTL apontado no cromossomo 2 não foi evidenciado em nenhum das análises anteriores e também demonstra baixo valor de LOD.full, portanto também não o incluimos no modelo.

# Testando automaticamente todas as interações
inte <- addint(mimulus.sim, qtl=rqtl, formula=y~Q1+Q2+Q3+Q1*Q3, pvalues=FALSE,
       pheno.col = 5)

As outras interações possíveis para o modelo não são significativas.

Calculo do intervalo de confiança conforme regra do 1.5LOD.

# Intervalo de confiança
lodint(rqtl, qtl.index=1)
##       chr      pos       lod
## CA220   3 115.0992  2.829366
## AG19    3 134.0170 12.631190
## AG19    3 134.0170 12.631190
lodint(rqtl, qtl.index=2)
##         chr      pos      lod
## BB186     5  0.00000 3.956317
## CC124     5 12.59292 4.357411
## MgSTS40   5 22.78345 2.479228
lodint(rqtl, qtl.index=3)
##          chr      pos       lod
## CA315     13 114.3072  7.219894
## CB55      13 137.7798 20.045915
## MgSTS104  13 153.8571 16.690566

Comparação com outros métodos de mapeamento

load("out.em.RData")

plotLodProfile(rqtl, col = "green", ylab="Profile LOD score")
plot(out.em,col="blue", chr = c(3,5,13), ylim = c(0,20), add=TRUE)
plot(out.cim.20,  col= "red", chr = c(3,5,13), add = TRUE)

O perfil de LOD em verde se refere ao mapemento por intervalo múltiplo, em vermelho ao mapeamento por intervalo composto e em azul ao mapeamento por intervalo.