1 Informações prévias
- 1.1 Mapa já construído
2 Descrição dos dados
- 2.1 Gráfico dos tipos de marcadores
3 Filtragem dos dados
- 3.1 Identificar marcadores redundantes
- 3.2 Identificar desvios na segregação
4 Construção do mapa
5 Avaliação individual dos grupos formados
6 Comparação com mapa publicado
7 Avaliação dos marcadores com desvios na segregação esperada
8 Novo mapa

1 Informações prévias

Inicialmente, Fishman et al. (2001) investigaram eventos de adaptação e especiação utilizando uma abordagem de identificação de locos com distorção de taxa de transmissão (TRDLs) com uma população experimental composta por 526 indivíduos \(F_2\). A população \(F_2\) foi sintetizada após a autofecundação de um híbrido \(F_1\) gerado a partir do cruzamento biparental interespecífico entre as linhagens endogâmicas puras IM62 (Mimulus guttatus) e SF5.4 (Mimulus nasutus). A população \(F_2\) foi semeada em vasos em condições de casa de vegetação em um delineamento inteiramente casualizado, e os dados fenotípicos foram coletados para 16 caracteres florais, vegetativos e reprodutivos. O mapa genético foi gerado inicialmente com 255 marcadores (27 microssatélites, 4 marcadores baseados em informação gênica, e 224 AFLP). O mapa genético foi construído com etapas exaustivas de análise gráfica e de reordenamento de marcadores até atingir um padrão biológico esperado com o software MAPMAKER 3.0. O mapa final apresentou um comprimento de 1780 cM (Kosambi) subdividido em 14 grupos de ligação, sendo filtrados 174 marcadores informativos no final do processo - os demais marcadores foram eliminados por inflacionar de forma incoerente o tamanho do mapa. Após a obtenção do mapa genético foram identificados 9 TRDLs com o software ANITA. Adicionalmente, também foi observado que 48% dos marcadores inseridos no mapa apresentavam eventos de distorção (desvio da segregação mendeliana esperada) significativas a um nível de 5% de significância. O padrão de alta frequência de locos com eventos de distorção sugeriu uma forte evidência de interação genômica heteroespecífica entre o genoma das duas espécies avaliadas, o que evidencia uma forte divergência entre espécies do gênero Mimulus.

Neste trabalho serão utilizados dados genotípicos da mesma população \(F_2\) (acesse-os aqui). Um subconjunto da população foi genotipada com novos marcadores incluindo os 27, 4, e 224 locos polimórficos com marcadores microssatélites (marcadores codominantes), marcadores baseados em informação gênica (marcadores codominantes), e marcadores AFLP (marcadores dominantes) já utilizados anteriormente. Detalhes sobre protocolos de genotipagem dos marcadores iniciais podem ser observados em Fishman et al. (2001).

1.1 Mapa já construído

2 Descrição dos dados

Iniciando a análise, podemos visualizar os dados brutos que utilizaremos para a construção do mapa. Para verificarmos qual será o formato de importação, foi realizada uma análise prévia do documento:

#Mostrando apenas as 14 linhas iniciais do arquivo:
head m_feb06.raw -n 14

## data type f2 intercross
## 287 418 16
## *AA461 
## DDDDBBBDDDD-D-DDDDBDBBBBDBBDBDBDBDBBDBBDDDDDDDBDDDDBBBDBDDBDBDDBDDBBDDDDDBBDDDBBBBD--DDBDDDDDBBDDDBDBBBDBDDDBDDDDDBDDDDDBDBBDBDDDDDDDDDDBDBBBD-BBBDBDBDDDDBDDDDDDDBD-DDDDD-BBB-DBBDDDBBDBDBBBDDDBDDBDDDBDBDDBDDBDDDDDDDD-DDBDBDDDDBBDDDDDDDDBBDDDDBDDD-DDDDDBBDBDDDDDBDDBDBDBDDD-DDDDDDDDBDDBD-
## *AA420 
## CCCCACCACCA-C-CCCACACAACCCCCCCACCCAACCCCCCCCCCCAAACCCACCCCCACCCCCACCACCCCCCCCCCACCC--CACAACACACCCCCCCCCCACCCCCCCCCCCCCCCCAACCCCCCCACCACAACCACC-ACCCCACCCCCCCCCCAACCC-CCCCC-CAC-CACACCCCCACCAACCCCCACCCCAACACCCCACCCCCCCA-CACCCCCACCCACACCCCCCCACCCCCCC-CCCCCAACCAACCCACACCACCCCC-CCACCCCCCACAC-
## *AA404 
## DDDDDDDDDDD-D-DBBDDDDBDDDBDDBDBBDBDDDDDDDDDDBDDDDBDDDDDDDDDDDDDBDDDDDDDDBDDDDDBDDDD--DDDBDDDDDDDBDDDDDBDBDDBDDDDDDDDDDBDDBDDBDDDDBDDDDDDDBDDDD-DDDBDDDDDBDDDDDBDDDDB-BDBDB-DDB-BDDBDDDBDDDBDDBBDDBBDDDBDDDDDBDDDDDDBBDDD-DDDBDBDDDDBDDDDDBBDDDDDDDDBDB-BDDDDDBDDDDDDBDBDBDDBDDDD-DBBDDBDDDDDDD-
## *AA384 
## DBDDDDBDBDB-D-BDDDBBDDDDDDDDDDDDDBBDBBDDBDDBDDDBDDDDDDBDBDDBDDBDDDDDDDBDBDDDDBBDDDD--DDBBDDDDDDDDDDDBDBDDDBDDBDBDDDBDBDDDDBDBDBDDDDDDBDDBBBBDD-DDBDDBBBDDBBBBDBDBDBD-DDDBD-BBD-DDDDBDDDDBDDDDBBBBDDDDDDDDDBDBDDBDDDDDDDD-DBDBDBBDDDBBDDBBDDDBBDDDDDDBB-DBDBDDDDDBDBDDDBBDDDBDDDB-DDDDDDDDDDBDD-
## *AA378 
## DDDBDBBDBDB-D-DBDDBBDBDDBBDBBBBDDDBBBDDDBBBBDDBDDDDBBDDBBDDDBDDDDDBBDDDDBBBDDBBDDDD--BDDDDDDDBBBDBDDDBBDBBDBBDBDDDDDDDDDDDDBDDDDDDDBDBDDDBBDDD-DDDDBBDDDDDDBDDDDDBDB-DDDBD-DBB-DDDBDBDDDBDDDDDBDDBDDDDDBDDDDDDBBDDBDBDDB-DBDDDDBDBDDBDBDBDDBDDDDBDBDDD-DDBBDDDDBDDBDBDDDDDBBBDBD-DDBDDBDBBDBDD-
## *AA371C 
## AHHHHBHBBAB-B-HBBBHBHAHBBBHHHBHAHBBAHHBHHBBHAAAAHHBHBHAHAABHHHBHBHHHBHBBBHHHBABBAHB--HHBHHBAHHHBHABHBBAHHBHBHHHHHHHBHBHHHHBBBAHAHAHAHHHHHBBAHH-HHBBHHBHABHAHHABBHHBB-BHHHA-HBH-AHAHHAABBHBHAABHBAHAHHHAABBHAHHHHHHABHHHA-HAHHHBHAAHHBHHHHHHHHHAHAHHABB-BHHABAHBBBAHBBABHHHAHBHHH-BHHBAHAABHABH-

Foi verificado que se trata de uma população F2 com 287 indivíduos, genotipados para 418 marcadores, com dados para 16 fenótipos. Como os dados pertencem a uma população F2 podemos encontrar seis classes genotípicas: AA, BB, AB, não BB, não AA, e dados perdidos. Nota-se que será possível encontrar dois tipos de marcadores, dominantes e codominantes. O primeiro é a codificação para marcadores codominantes que indicam os indivíduos homozigotos para um alelo (A), os heterozigotos (H) e os homozigotos para o outro alelo (B). A segunda é a classe dos marcadores dominantes, que pode ocorrer de duas formas: (i) a presença de homozigotos para o alelo (A) e de indivíduos não homozigotos para esse alelo (C); (ii) indivíduos homozigotos para o alelo (B) e de não homozigotos para esse alelo (D).

2.1 Gráfico dos tipos de marcadores

Utilizando o OneMap é possível verificar a distribuição das classes pelos marcadores e indivíduos, conforme o gráfico abaixo:

#Carregar o pacote OneMap
library(onemap)

#Ler os dados do formato do Mapmaker
data <- read_mapmaker(file = "m_feb06.raw")

##  --Read the following data:
##  Type of cross:          f2 
##  Number of individuals:  287 
##  Number of markers:      418 
##  Missing trait values:       
##    fl: 11 
##    fs: 19 
##    ft: 11 
##    ll: 20 
##    nv: 12 
##    pa: 35 
##   pal: 11 
##    pv: 12 
##    sa: 11 
##    sl: 11 
##    ss: 28 
##    tl: 11 
##    tp: 12 
##    tw: 11 
##    vi: 12 
##    ww: 11

#Gráfico da distribuição das classes dos marcadores:
plot.onemap(data)

## No id variables; using all as measure variables

É possível perceber a segregação de marcadores dominantes e codominantes, bem como a ocorrência de dados perdidos.

3 Filtragem dos dados

3.1 Identificar marcadores redundantes

Para identificar a presença de marcadores redundantes utilizou-se a função find_bins, cujo resultado apontou que não haviam redundâncias.

#Exemplo do comando find_bins
data_bins <- find_bins(data)

3.2 Identificar desvios na segregação

data_seg <- test_segregation(data)
plot.onemap_segreg_test(data_seg)

# Selecionando apenas os marcadores com segregação
data_index_distorted <- select_segreg(data_seg, distorted = TRUE)

# Retirando os marcadores com desvio do arquivo inicial
index_all <- colnames(data$geno)
temp <- index_all %in% data_index_distorted
num_distorted <- which(temp == TRUE)
num_ndistorted <- which(temp == FALSE)

Aqui manipularemos o arquivo de entrada para separar os marcadores com desvios na segregação esperada em um arquivo diferente. Para facilitar, disponibilizamos os arquivos depois de manipulados: arquivo com os marcadores sem desvio na segregação e arquivo com os marcadores com desvios na segregação.

Retiramos com editor de texto o cabeçalho e os fenótipos do arquivo. Foi feito um script em bash para inserir um tab entre o nome do marcador e os genótipos e facilitar a manipulação do arquivo:

k=1
for i in $(cat data_whead)
do
   printf -- $i
   printf "\t"
   if(($k % 2 == 0)); then
      printf "\n"
   fi
   ((k++))
done

chmod 777 make_tab
./make_tab.sh > data_whead_t.txt

Dentro do ambiente R, o arquivo ‘data_whead_t.txt’ foi interpretado como uma data.frame e removemos os marcadores com desvios na segregação.

data_whead <- read.table(file = "data_whead_t.txt", sep="\t")
head(data_whead)

data_whead_wdistorted <- data_whead[-num_distorted,-3]
data_whead_cdistorted <- data_whead[-num_ndistorted, -3]

dim(data_whead)
dim(data_whead_wdistorted)
dim(data_whead_cdistorted)

# Produzindo novos arquivos apenas com os marcadores sem desvio na segregação

write.table(x = data_whead_wdistorted, file = "data_head_wdistorted", sep = "\n", quote = FALSE, col.names = FALSE, row.names = FALSE)


# Produzindo novos arquivos apenas com os marcadores com  desvio na segregação

write.table(x = data_whead_cdistorted, file = "data_head_cdistorted", sep = "\n", quote = FALSE, col.names = FALSE, row.names = FALSE)

# Adicionamos manualmente o cabeçalho e os fenótipos no arquivo

data1 <- read_mapmaker(file = "data_with_head.raw")

##  --Read the following data:
##  Type of cross:          f2 
##  Number of individuals:  287 
##  Number of markers:      356 
##  Missing trait values:       
##    fl: 11 
##    fs: 19 
##    ft: 11 
##    ll: 20 
##    nv: 12 
##    pa: 35 
##   pal: 11 
##    pv: 12 
##    sa: 11 
##    sl: 11 
##    ss: 28 
##    tl: 11 
##    tp: 12 
##    tw: 11 
##    vi: 12 
##    ww: 11

# Mantemos essas informações para testarmos o posicionamento dos marcadores com desvios na segregação ao final da construção do mapa

Pelo teste de segregação com correção por Bonferroni 15% dos dados ou 62 marcadores apresentaram distorção na segregação. Esses marcadores não serão considerados nas etapas inciais da construção do mapa. Esses marcadores serão considerados após ordenados todos aqueles que não sofreram desvio de segregação, para não prejudicarem a construção do mapa.

4 Construção do mapa

4.1 Teste de dois pontos

Incialmente, são calculadas as frações de recombinação a cada dois marcadores com a função ‘rf2points’. Aqui já estabelecemos um limiar para considerar ligação entre os marcadores, com o valor de LOD recomendado conforme o número de testes realizados (número de marcadores) e da fração de recombinação.

O valor de LOD recomendado por ser obtido pela função ‘suggest_lod’.

LOD <- suggest_lod(data1)

rf2points <- rf_2pts(data1, LOD = LOD, max.rf = 0.5)

## Computing 63190 recombination fractions:
## 
##  0%  ..................................................  89%
##  89% .....................   100%

4.2 Agrupamento

Os marcadores são então grupados conforme a fração de recombinação calculados e os critérios de LOD e fração de recombinação estabelecidos anteriormente. Segundo o tutorial do programa, é possível alterar os critérios sobrepondo os argumentos, mas manteremos os já determinados na função ‘rf2points’.

Seq.Mark <- make_seq(rf2points, "all")

LGs <- group(Seq.Mark)

##    Selecting markers: 
##    group    1 
##     .......
##    group    2 
##     ......................
##    group    3 
##     ............................................................
##     .......................
##    group    4 
##     .......................................
##    group    5 
##     ...............
##    group    6 
##     ...................................................
##    group    7 
##     ................................
##    group    8 
##     ..........................
##    group    9 
##     .............
##    group    10 
##     ......
##    group    11 
##     ........................
##    group    12 
##     .....................
##    group    13 
##     ......
##    group    14 
##     ...

4.3 Rapid Chain delineation

Uma vez formados os grupos, podemos ter uma visão inicial deles aplicando um teste de dois pontos. Dessa forma, poderemos identificar avaliando os Heatmaps gerados se houve agrupamentos incentos.

# Primeira visualização com RCD

LGs.RCD <- vector(l = c(LGs$n.groups), "list")
for(i in 1:LGs$n.groups)
{
  LG.Temp <- make_seq(LGs, i)
  
  png(paste("LG", i, "_RCD_", length(LG.Temp$seq.num), 
            "_Markers_First_Check.png", sep = ""))
  LG.Temp.RCD <- rcd(LG.Temp)
  rf_graph_table(LG.Temp.RCD, scale = 1.5, main = paste("LG", i, "- RCD"),
                 inter = FALSE)
  dev.off()
  
  LGs.RCD[[i]] <- LG.Temp.RCD
}

save(LGs.RCD, file = "LGs_RCD_Markers_First_Check.RData")

Com a observação dos Heatmaps gerados, podemos averiguar que o LG3 provavelmente é composto por vários grupos.

Após essa verificação inicial e mais simples, realizamos a abordagem multiponto para maior exploração.

4.4 Cadeia de Markov

  LGs.HMM.safe <- LGs.HMM.force <- vector(l = c(LGs$n.groups), "list")
  for(i in 1:LGs$n.groups)
  {
    LG.Temp <- make_seq(LGs, i)
    
    png(paste("LG", i, "_HMM_", length(LG.Temp$seq.num), 
              "_Markers_Second_Check.png", sep = ""))
    LG.Temp.HMM <- order_seq(LG.Temp)
    LG.Temp.seq.HMM <- make_seq(LG.Temp.HMM, "force")
    LG.Temp.seq.HMM.safe <- make_seq(LG.Temp.HMM, "safe")
    rf_graph_table(LG.Temp.seq.HMM, scale = 1.5, main = paste("LG_force", i, "- HMM"),
                   inter = FALSE)
    rf_graph_table(LG.Temp.seq.HMM.safe, scale = 1.5, main = paste("LG_safe", i, "- HMM"),
                   inter = FALSE)
    dev.off()
    LGs.HMM.force[[i]] <- LG.Temp.seq.HMM
    LGs.HMM.safe[[i]] <- LG.Temp.seq.HMM.safe
  }
  
  save(LGs.HMM.force, file = "LGs_force_HMM_Markers_First_Check.RData")
  save(LGs.HMM.safe, file = "LGs_safe_HMM_Markers_First_Check.RData")

5 Avaliação individual dos grupos formados

Avaliando esse segundo ordenamento é possível identificar marcadores mal posicionados em cada um dos grupos. Os grupos foram avaliados separamente de forma criteriosa. De forma geral, foram retirados os marcadores mal posicionados e averiguada a correspondência com o mapa já construído. A descrição mais detalhada desta etapa para cada um dos grupos conforme o novo conjunto de dados pode ser obtida nos links abaixo:

GL1
GL2
GL3
GL4
GL5
GL6
GL7
GL8
GL9
GL10 - Este grupo não foi modificado nesta etapa por apresentar poucos marcadores
LG11
GL12
GL13
GL14

6 Comparação com mapa publicado

Em seguida, verificamos a correspondência dos novos grupos gerados com o mapa já publicado.Foi observado que nove (2,4,5,6,7,8, 10, 11, 12) dos grupos formados tinham grande correspondência, dois pares de grupos (9 e 14; 1 e 13) formariam apenas um grupo de ligação no mapa publicado. Um dos grupos (3) correponde a três outros grupos de ligação do mapa publicado. A partir de agora iremos nos referir a grupos de ligação com a numeração já estabelecida para o mapa publicado, conforme tabela a seguir.

Grupos no mapa publicado	Grupos correspondentes no novo conjunto de dados
1	1 e 13
2	2
3	3
4	4
5	5
6	6
7	3
8	7
9	8
10	3
11	10
12	11
13	12
14	9 e 14

Nas descrições presentes nos links a seguir foram feitas as separações/junções conforme a tabela anterior.

7 Avaliação dos marcadores com desvios na segregação esperada

Uma vez esbelecidos os grupos conforme o mapa publicado, realizamos novamente o agrupamento dos marcadores, mas dessa vez considerando os marcadores que apresentaram desvios na segregação. Esta etapa esta descrita no seguinte link:

Avaliando marcadores com desvios na segregação esperada

8 Novo mapa

Assim foi gerado um novo mapa contendo 346 marcadores, ou seja, 172 marcadores a mais que o mapa já publicado. Portanto, do total de 418 marcadores desse conjunto de dados 82.78% foram informativos, sendo uma taxa superior do que foi obtida no primeiro trabalho, em que do total de 255 marcadores apenas 68.24% foram informativos. As distâncias em centimorgans foram estabelecidas pela função de Kosambi.

Abaixo apresentamos os Heatmaps da versão final do mapa.

load("mapa_datas.RData")
data1 <- read_mapmaker(file = "data_with_head.raw")

##  --Read the following data:
##  Type of cross:          f2 
##  Number of individuals:  287 
##  Number of markers:      356 
##  Missing trait values:       
##    fl: 11 
##    fs: 19 
##    ft: 11 
##    ll: 20 
##    nv: 12 
##    pa: 35 
##   pal: 11 
##    pv: 12 
##    sa: 11 
##    sl: 11 
##    ss: 28 
##    tl: 11 
##    tp: 12 
##    tw: 11 
##    vi: 12 
##    ww: 11

data_seg2 <- read_mapmaker(file = "data_with_distorted_2.raw")

##  --Read the following data:
##  Type of cross:          f2 
##  Number of individuals:  287 
##  Number of markers:      62

data <- combine_onemap(data1, data_seg2) 

LOD1 <- suggest_lod(data)
rf2points1 <- rf_2pts(data, LOD = LOD1, max.rf = 0.5)

## Computing 87153 recombination fractions:
## 
##  0%  ........................................... 100%

for (i in 1:14){
rf_graph_table(mapa[[i]], scale = 1.5, main = paste("LG", i, "- HMM"),
               inter = FALSE, axis.cex = 0.75)
}

O software MapChart foi utilizado para montar as figuras abaixo:

Mapa Mimulus

Cristiane Taniguti, Fernando Correr, Fernando Espolador, Jhonathan Pedroso

Maio de 2017