Lista 10 - Mapeamento de QTLs IV
Cristiane Taniguti, Fernando Correr, Fernando Espolador, Jhonathan dos Santos
June 2, 2017
- Considere o conjunto de dados de Mimulus, já usado em aulas anteriores. Usando o mapa que você construiu anteriormente com o OneMap, faça mapeamento por intervalo e por intervalo composto para localizar QTL’s para uma das características presentes no conjunto de dados (escolha a que julgar mais interessante).
- Determine o limiar para rejeitar H0 e apresente os resultados de forma gráfica e em tabelas de forma semelhante aos artigos publicados na área.
- Use o software R/qtl, verificando no tutorial do OneMap como exportar o mapa construído para o formato do R/qtl. Use o Fórum para resolver suas dúvidas.
- O exercício pode ser feito em grupo, pela mesma equipe que fez o mapa genético.
- Apresente os resultados em um arquivo html
- Faça um vídeo no Youtube apresentando e interpretando os resultados. Coloque o link do vídeo no Fórum.
A construção do mapa genético do Mimulus foi realizada em uma atividade anterior da mesma disciplina de Biometria de Marcadores Genéticos da Esalq. Mais informações podem ser encontradas acessando o link.
O primeiro passo desta atividade é exportar o mapa através da função write_map do pacote onemap, de forma que possa ser um dos arquivos de entrada do R/qtl.
#Carregando os pacotes onemap e qtl:
library(onemap)
library(qtl)
library(knitr)
# Carregando informações do mapa construído
load("mapa_datas.RData")
# Exportando informações para arquivo de entrada do R/qtl
write_map(mapa, "mimulus.map")Mapeamento por intervalo
No mapeamento de QTLs através do R/qtl é necessário, além do mapa genético, o arquivo bruto do formato MAPMAKER. Isso é necessário para que ele possa reconhecer o tipo de cruzamento, os marcadores e fenótipos envolvidos no experimento. Será verificado, como no exercício de mapeamento genético, que se trata de uma população \(F_2\) onde foram genotipados 418 marcadores em 287 indivíduos e foram avaliados 16 fenótipos que consistem em características florais, reprodutivas e vegetativas (melhor descritas em Fishman et al. 2002).
#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#Observar o mapa genético:
plot.map(mimulus.qtl)
#Reestimar o mapa - checagem do mapa:
newmap <- est.map(mimulus.qtl, tol = 1e-6, map.function = "kosambi")
#Nota-se que os mapas são praticamente idênticos:
plot.map(mimulus.qtl, newmap)
A função est.map faz a reestimativa do mapa genético de acordo com o algoritmo Lander-Green. Comparando o mapa construído (à esquerda) com o mapa reestimado, nota-se que são muito similares, sendo notada uma diminuição mais aparente dos grupos de ligação 1, 5, 8 e 10. Nos grupos 1, 4, 5 e 8 é possível observar ou a separação entre dois marcadores muito próximos ou a aproximação de dois marcadores.
Agora, será realizado o mapeamento por intervalo com o algoritmo EM e a regressão Haley-Knott. Inicialmente é essencial calcular as probabilidades condicionais de cada genótipo através da função calc.genoprob. Serão considerados “passos” de 1cM e será assumida uma taxa de erro de genotipagem de 1%. Em seguida, é feito o escaneamento do genoma para buscar QTL. Aqui consideraremos apenas um dos fenótipos avaliados, o número de pólens não viáveis (nv).
Os objetos out.em e out.hk indicam os métodos EM e Haley-Knott utilizados para a busca de QTLs para essa característica.
#Realizando o mapeamento por intervalo - Com os algoritmo EM e regressão Haley-Knott:
mimulus.qtl.prob <- calc.genoprob(mimulus.qtl, step=1, error.prob=0.01)
##EM:
out.em <- scanone(mimulus.qtl.prob, method = "em", pheno.col = 5)
##Haley-Knott:
out.hk <- scanone(mimulus.qtl.prob, method="hk", pheno.col = 5)
# Maiores valores de LOD para os dois casos
#EM:
kable(summary(out.em))| chr | pos | lod | |
|---|---|---|---|
| CA389C | 1 | 80.17989 | 2.0840674 |
| c2.loc27 | 2 | 27.00000 | 2.2556121 |
| AG19 | 3 | 134.01698 | 2.8268231 |
| c4.loc147 | 4 | 147.00000 | 1.3534722 |
| BB186 | 5 | 0.00000 | 2.3170831 |
| c6.loc36 | 6 | 36.00000 | 1.6093390 |
| c7.loc129 | 7 | 129.00000 | 1.4340664 |
| c8.loc134 | 8 | 134.00000 | 0.9824759 |
| MgSTS536 | 9 | 28.78428 | 2.0666186 |
| c10.loc211 | 10 | 211.00000 | 1.6002020 |
| MgSTS561 | 11 | 113.14222 | 0.9843905 |
| CB216 | 12 | 13.60661 | 0.6481818 |
| c13.loc149 | 13 | 149.00000 | 12.8935216 |
| c14.loc214 | 14 | 214.00000 | 1.5905898 |
#Haley-Knott
kable(summary(out.hk))| chr | pos | lod | |
|---|---|---|---|
| CA389C | 1 | 80.17989 | 2.0480840 |
| c2.loc30 | 2 | 30.00000 | 2.4707325 |
| AG19 | 3 | 134.01698 | 2.8463443 |
| c4.loc147 | 4 | 147.00000 | 1.3454197 |
| c5.loc5 | 5 | 5.00000 | 2.1486857 |
| c6.loc36 | 6 | 36.00000 | 1.6529353 |
| c7.loc130 | 7 | 130.00000 | 1.4944047 |
| c8.loc134 | 8 | 134.00000 | 0.9979106 |
| MgSTS536 | 9 | 28.78428 | 2.1166201 |
| c10.loc210 | 10 | 210.00000 | 1.6173384 |
| c11.loc112 | 11 | 112.00000 | 1.0283502 |
| CB216 | 12 | 13.60661 | 0.6430561 |
| c13.loc151 | 13 | 151.00000 | 11.6791235 |
| c14.loc214 | 14 | 214.00000 | 1.6113803 |
As duas tabelas apresentam os maiores valores de LOD para cada um dos grupos de ligação. Há uma correspondência entre as posições apontadas pelos dois métodos, se nos intervalos a posição não foi a mesma entre os métodos, é possível verificar que são bem aproximadas.
Para detectar que o LOD é significativo, ou seja, maior do que um determinado limiar, será utilizada a estratégia de permutações. Utiliza-se a função scanone com 5000 permutações com a regressão Haley-Knott considerando o fenótipo que já estávamos utilizando. Nas permutações quebra-se a associação entre genótipo e fenótipo, de modo que obteremos uma distribuição empírica dos dados e definiremos o limiar como aquele que indica o percentil 95%.
# Calculando o limiar de LOD por permutações pelo modelo de hk
operm.hk <- scanone(mimulus.qtl.prob, method="hk", n.perm=5000, pheno.col = 5)## Doing permutation in batch mode ...thr <- summary(operm.hk, alpha=0.05)
thr <- as.numeric(thr)
#Fazendo o gráfico do perfil de LOD score
## Linha azul é por EM e verde por hk
plot(out.em, out.hk, col=c("blue", "green"))
#Utilizando o limiar
abline(h = thr, col = "red")
# Intervalo de confiança, considerando regra de 1 LOD
intervalo <- lodint(out.em, chr=13, drop=1)
knitr::kable(intervalo)| chr | pos | lod | |
|---|---|---|---|
| c13.loc129 | 13 | 129 | 11.56160 |
| c13.loc149 | 13 | 149 | 12.89352 |
| c13.loc153 | 13 | 153 | 11.80792 |
# Localização dos QTLs no mapa
QTLs <- summary(out.em, threshold = thr)
mimulus.sim <- sim.geno(mimulus.qtl.prob)
qtlsim <- makeqtl(mimulus.sim, chr = QTLs$chr, pos = QTLs$pos)
plot(qtlsim)
Por meio do gráfico do perfil de LOD score no genoma pelo método com o uso do EM, verifica-se um pico superior ao limiar, localizado no grupo de ligação 13. Com a função effectplot é possível obter estimativas da média fenotípica específica de cada genótipo para o marcador vinculado ao QTL. Com tais valores obtemos os efeitos genéticos aditivos e de dominância do QTL significativo:
\(a = (\mu_{BB} - \mu_{AA})/2\)
\(d = \mu_{AB} - (\mu_{AA} + \mu_{BB})/2\)
# Localizando marcadores relacionados com os maiores picos
picos <- find.marker(mimulus.qtl.prob, QTLs$chr, QTLs$pos)
picos## [1] "MgSTS104"effect_a <- effect_d <- vector()
for (i in 1:length(picos)){
# gráfico que mostra a diferença entre as médias
#effectplot(mimulus.qtl.prob, mname1=picos[i])
# plotando mais detalhes
#plot.pxg(mimulus.qtl.prob, marker=picos[i])
# Obtendo médias para cada genótipo
eff <- effectplot(mimulus.sim, mname1=picos[i], draw=TRUE, pheno.col = 5)
# Calculando efeitos aditivo e de dominância
effect_a <- c(effect_a, (eff$Means[3] - eff$Means[1])/2)
effect_d <- c(effect_d, eff$Means[2] - (eff$Means[1] + eff$Means[3])/2)
}
# Tais valores também podem ser obtidos com a função `fitqtl`
model1 <- fitqtl(mimulus.sim, pheno.col = 5, qtl = qtlsim, get.ests = TRUE)
resumo <- summary(model1)
IM_QTLs <- data.frame("marcador"= rownames(QTLs),"cromossomo"= QTLs$chr,
"posição" = QTLs$pos, "LOD" = QTLs$lod,
"efeito aditivo" = resumo$ests[2,1] ,
"efeito de dominância"= resumo$ests[3,1],
"R^2" = resumo$result.full[1,5])
knitr::kable(IM_QTLs) | marcador | cromossomo | posição | LOD | efeito.aditivo | efeito.de.dominância | R.2 |
|---|---|---|---|---|---|---|
| c13.loc149 | 13 | 149 | 12.89352 | -43.23673 | -24.51985 | 16.35856 |
A próxima etapa é verificar as possíveis alterações neste perfil ao adicionarmos covariáveis no modelo.
Mapeamento por intervalo composto
No mapeamento por intervalo composto serão utilizados marcadores como covariáveis durante a estimativa do efeito nas posições. É importante ressaltar que essas covariáveis são utilizadas ao longo de todo o genoma e bem como no mesmo grupo de ligação onde está sendo escaneado o QTL. Nesse último caso, deve ser definido um tamanho de janelas para “isolar” a região próxima à posição testada.
A função cim realiza a seleção forward para fixar um número de marcadores (especificados pelo argumento n.marcovar), seguido por mapeamento por intervalo, omitindo os marcadores definidos como covariáveis presentes a uma distância fixa (determinada pelo argumento window) do ponto testado.
Para definir o número de covariáveis, foram investigados os resultados de regressão nas marcas e também o número de picos obtidos através do mapeamento por intervalo. Sobre o número de picos na análise anterior foi escolhido abranger apenas o pico considerado significativo localizado no cromossomo 13.
Sobre a regressão nas marcas, definimos um LOD de 3 na tentativa de abranger mais marcadores ao longo do genoma. Esse valor de LOD foi também evidenciado apenas o pico no cromossomo 13.
A decisão tomada foi utilizar apenas um marcador como covariável. Ressalta-se, mais uma vez, que essas investigações foram apenas para definir um número de covariáveis para ser empregado no mapeamento por intervalo composto.
############################################
#Utilizando a regressão nas marcas:
out.mr <- scanone(mimulus.qtl, method="mr", pheno.col = 5)
#Verificando aquelas superior a um LOD de 2.5:
kable(summary(out.mr, threshold=3))| chr | pos | lod | |
|---|---|---|---|
| MgSTS104 | 13 | 153.8571 | 10.06273 |
############## CIM ##################
out.cim.20 <- cim(mimulus.qtl.prob, n.marcovar=1, window=15, pheno.col = 5)###################### LIMIAR
permic_5 <- cim(mimulus.qtl.prob, n.marcovar = 1, window = 15, pheno.col = 5, n.perm = 5000, method = "hk") # Processamento demorado############### Gráfico do perfil de LOD
lim <- summary(permic_5, alpha = 0.05)
plot(out.em,col="blue", ylim = c(0,15));
plot(out.cim.20, col= "red", add = TRUE);
add.cim.covar(out.cim.20, col="green")
abline(h=lim[1], col = "orange")
plot(out.em,col="blue", chr = c(5,13), ylim = c(0,15));
plot(out.cim.20, col= "red", chr = c(5,13), add = TRUE);
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.511880 |
| BB186 | 5 | 0 | 4.511880 |
| c5.loc17 | 5 | 17 | 3.497126 |
| c13.loc130 | 13 | 130 | 2.729978 |
| c13.loc135 | 13 | 135 | 12.860446 |
| c13.loc146 | 13 | 146 | 2.580033 |
Para determinar os efeitos dos QTLs encontrados, é necessário considerar os efeitos da covariável utilizada no mapeamento por intervalo composto:
covar <- attr(out.cim.20,"marker.covar.pos")
covar## chr pos
## CB55 13 137.7798Também consideraremos os QTLs significativos encontrado.
QTL <- summary(out.cim.20, threshold = lim[1])
QTL## chr pos lod
## BB186 5 0 4.51
## c13.loc135 13 135 12.86Com função fitqtl montamos um modelo considerando os efeitos da covariável e do QTL (que no caso coincidem). Dessa forma podemos identificar o efeito do QTL.
qtls <- makeqtl(mimulus.sim, chr = QTL$chr, pos = QTL$pos)
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.51188 | -26.44937 | 2.834529 | 4.897425 |
| 13@137.8 | c13.loc135 | 13 | 135 | 12.86045 | -46.21984 | -7.214839 | 16.155726 |
Pelo modelo, foram identificados os efeitos de aditividade, de dominância e porcentagem de variação do único QTL significativo pelo mapeamento por intervalo composto.