Introdução

A estatística ABBA-BABA (também conhecida como "estatística D") oferece um teste simples e poderoso para desvios de uma história evolutiva estritamente bifurcada. Portanto, é frequentemente utilizada para testar introgressão genômica utilizando dados de SNP em escala genômica (por exemplo, de sequenciamento de genoma completo ou RADseq).

Neste tutorial, combinaremos o uso de software disponível com código originalmente desenvolvido em R para realizar a análise ABBA-BABA. Analisaremos dados genômicos de várias populações de borboletas Heliconius.

Fluxo de Trabalho

Partindo dos dados de genótipo de múltiplos indivíduos, primeiro inferimos a frequência de alelos para cada SNP. Em seguida, calculamos a estatística D, seguido pelo uso do método block jackknife para testar desvios significativos da expectativa nula de D=0. Finalmente, estimamos a proporção de "introgressão" (f).

Conjunto de Dados

Estudaremos múltiplas espécies de três espécies: Heliconius melpomene, Heliconius timareta e Heliconius cydno. As distribuições dessas espécies se sobrepõem parcialmente, e acredita-se que haja hibridação em áreas de co-ocorrência. Nosso conjunto de amostras inclui dois pares de populações simpátricas de H. melpomene e H. cydno do Panamá e da encosta oeste dos Andes na Colômbia. Na encosta leste dos Andes na Colômbia e no Peru, existem dois outros pares de populações simpátricas: H. melpomene e H. timareta. Finalmente, temos duas amostras da espécie de grupo externo Heliconius numata, necessária para a análise ABBA-BABA.

Todas as amostras foram sequenciadas por sequenciamento de genoma completo em profundidade, e genótipos foram obtidos para cada locus do genoma de cada indivíduo usando um pipeline padrão. Os dados foram filtrados para reter apenas polimorfismos de nucleotideo único bi-alélicos (SNP), e esses dados foram ainda mais refinados para reduzir o tamanho do arquivo para este tutorial.

Pressupostos

Assumimos que o hibridação entre espécies simpátricas resultará em compartilhamento de variação genética entre H. cydno e as raças simpátricas de H. melpomene do oeste, e entre H. timareta e as raças simpátricas correspondentes de H. melpomene dos Andes orientais. Também há outra raça de H. melpomene da Guiana Francesa, que é simpátrica com H. timareta e H. cydno, e que não deve ter passado por troca genética recente com essas duas espécies, servindo como controle.

Além de testar a presença de introgressão, também testaremos a hipótese de que certas partes do genoma experimentam mais introgressão do que outras. Especificamente, sabemos que há pelo menos um locus no cromossomo Z que causa esterilidade feminina em híbridos entre essas espécies, sugerindo incompatibilidade entre cromossomos autossômicos de uma espécie e o cromossomo Z de outra. Portanto, esperamos menos introgressão no cromossomo Z em comparação com autossomos.

Teste de Introgressão Genômica Completa

Na sua formulação mais simples, o teste ABBA-BABA depende da contagem de loci que correspondem aos padrões de genótipo ABBA e BABA. Ou seja, dado três grupos internos e um grupo externo com relação (((P1,P2),P3),O), e dado sequências de genoma único representando cada grupo (ou seja, H1, H2 e H3), os loci ABBA são aqueles onde H2 e H3 compartilham alelo derivado ("B"), enquanto H1 tem estado ancestral ("A"), conforme definido pela amostra do grupo externo. Da mesma forma, BABA representa loci onde H1 e H3 compartilham alelo derivado.

Ignorando mutações recorrentes, esses padrões de SNP só podem ocorrer se certas partes do genoma não sigam a "árvore de espécies", mas agrupem H2 com H3 ou H1 com H3. Se as populações se dividiram recentemente, devido a variações na ordenação da linhagem, espera-se que partes do genomo mostrem essa "inconsistência" na linhagem. Na ausência de desvios de uma topologia estritamente bifurcada, esperamos que aproximadamente a mesma proporção do genoma mostre as duas linhagens inconsistentes (((H2,H3),H1),O) e (((H1,H3),H2),O). Portanto, contando SNP ABBA e BABA em todo o genoma (ou em uma grande parte dele), podemos estimar a proporção do genoma representada por essas duas linhagens inconsistentes, o que significa que esperamos uma proporção de 1:1 de SNP ABBA e BABA. Por exemplo, o fluxo gênico entre as populações P3 e P2 pode causar esse desvio, embora também possa indicar outros fenômenos que quebram nossas suposições, como estrutura populacional ancestral ou taxas de substituição variáveis.

Para quantificar o desvio da proporção esperada, calculamos D, que é a diferença entre as somas dos padrões ABBA e BABA no genoma inteiro dividida pela sua soma:

D = [sum(ABBA) - sum(BABA)] / [sum(ABBA) + sum(BABA)]

Portanto, D varia de -1 a 1 e, sob a hipótese nula, deve ser igual a 0. D>0 indica excesso de ABBA, enquanto D<0 indica excesso de BABA.

Se temos múltiplas amostras de cada grupo, a contagem de loci ABBA e BABA não é tão simples. Uma opção é considerar apenas loci onde todas as amostras do mesmo grupo compartilham o mesmo alelo, mas isso descartaria uma grande quantidade de dados úteis. Uma melhor opção é usar as frequências de alelo de cada locus para quantificar o grau de tendência da linhagem em direção aos padrões ABBA ou BABA. Isso equivale essencialmente a usar todos os possíveis grupos de quatro haplótipos genômicos de cada locus para calcular SNP ABBA e BABA. Portanto, ABBA e BABA não são mais estados binários, mas números entre 0 e 1 que representam a frequência de combinações de alelos que correspondem a cada linhagem. Eles são calculados com base nas frequências do alelo derivado (p) e ancestral (1-p) em cada população, da seguinte forma:

ABBA = (1-p1) x p2 x p3 x 1-pO

BABA = p1 x (1-p2) x p3 x 1-pO

Prática

Preparação

Abra uma janela de terminal e navegue até a pasta onde o exercício será executado e todos os arquivos de dados de entrada e saída serão armazenados. Agora crie um subdiretório chamado "dados" e baixe os arquivos de dados necessários para este tutorial

mkdir dados

cd dados

wget https://github.com/simonhmartin/tutorials/raw/master/ABBA_BABA_whole_genome/data/hel92.DP8MP4BIMAC2HET75dist250.geno.gz

wget https://github.com/simonhmartin/tutorials/raw/master/ABBA_BABA_whole_genome/data/hel92.pop.txt

cd ..

Em seguida, baixe o conjunto de scripts Python necessários para este tutorial no GitHub

wget https://github.com/simonhmartin/genomics_general/archive/master.zip
unzip master.zip

Frequências de Alelo Genômica Completa

Para calcular esses valores com base nos dados genômicos populacionais, primeiro precisamos determinar a frequência do alelo derivado em cada polimorfismo em cada população no genoma. Calcularemos esses valores usando os dados genotípicos de Heliconius fornecidos com scripts Python. O arquivo de entrada foi filtrado para conter apenas loci bi-alélicos. O script de frequência exige que definamos as populações. Essas são definidas no arquivo hel92.pop.txt.

python genomics_general-master/freq.py -g dados/hel92.DP8MP4BIMAC2HET75dist250.geno.gz \
-p mel_mel -p mel_ros -p mel_vul -p mel_mal -p mel_ama \
-p cyd_chi -p cyd_zel -p tim_flo -p tim_txn -p num \
--popsFile dados/hel92.pop.txt --target derived \
-o dados/hel92.DP8MP4BIMAC2HET75dist250.derFreq.tsv.gz

Ao definir --target derivado, obtemos a frequência do alelo derivado em cada locus para cada população. Isso é baseado na população final especificada (H. numata silvana, ou "slv") como grupo externo. Populações com estados ancestrais não fixos em um locus serão descartadas.

Análise ABBA-BABA Genômica Completa

Para entender como o teste ABBA-BABA funciona, desenvolveremos código do zero para realizar o teste. Inicie um novo script R. Isso facilitará a reexecução de toda a análise usando diferentes populações.

• Função R para análise ABBA-BABA

Primeiro, definimos uma função para calcular a proporção ABBA e BABA para cada locus e usá-los para calcular a estatística D. A entrada será a frequência do alelo derivado nas populações P1, P2 e P3 (ou seja, p1, p2 e p3). (A frequência do alelo ancestral no grupo externo é 1 em todos os loci, pois usamos o grupo externo para identificar o alelo ancestral, portanto pode ser ignorada).

estatistica_D <- function(freq_p1, freq_p2, freq_p3) {
    ABBA <- (1 - freq_p1) * freq_p2 * freq_p3
    BABA <- freq_p1 * (1 - freq_p2) * freq_p3
    (sum(ABBA, na.rm=TRUE) - sum(BABA, na.rm=TRUE)) / (sum(ABBA, na.rm=TRUE) + sum(BABA, na.rm=TRUE))
    }

• Lendo nossos dados de frequência de alelo.

tabela_frequencia = read.table("dados/hel92.DP8MP4BIMAC2HET75dist250.derFreq.tsv.gz", header=TRUE, as.is=TRUE)

Isso cria um objeto chamado tabela_frequencia, que contém a frequência do alelo derivado para cada SNP. Podemos verificar o número de loci nesta tabela e também visualizar as primeiras linhas para entender os dados.

nrow(tabela_frequencia)

head(tabela_frequencia)

Observe que as duas primeiras colunas fornecem o nome do scaffold (ou seja, cromossomo) e a posição de cada locus no cromossomo. As colunas restantes são as frequências de alelo de diferentes subespécies, conforme mostrado acima.

• Estatística D

Agora, para calcular D, precisamos definir as populações P1, P2 e P3. Começaremos com um caso de teste óbvio e publicado anteriormente: perguntaremos se há evidências de introgressão entre H. melpomene rosina (mel_ros) e H. cydno chioneus (cyd_chi). Eles são P2 e P3, respectivamente. P1 será nossa população de espécies diferentes, H. melpomene melpomene da Guiana Francesa (mel_mel).

Definimos essas populações e então calculamos D extraindo as frequências do alelo derivado para todos os SNP dessas três populações.

pop1 <- "mel_mel"
pop2 <- "mel_ros"
pop3 <- "cyd_chi"

D <- estatistica_D(tabela_frequencia[,pop1], tabela_frequencia[,pop2], tabela_frequencia[,pop3])

print(paste("D =", round(D,4)))

Obtemos uma estatística D (lembrando que D varia de -1 a 1) que indica excesso de ABBA sobre BABA. Isso sugere que H. cydno chioneus (cyd_chi) do Panamá compartilha mais variação genética com H. melpomene rosina (mel_ros) simpátrico do Panamá do que com H. melpomene melpomene (mel_mel) alopátrico da Guiana Francesa. Isso é consistente com hibridação e fluxo gênico entre essas duas espécies simpátricas.

No entanto, atualmente não sabemos se esse resultado é estatisticamente confiável. Especificamente, não sabemos se o excesso de ABBA está uniformemente distribuído em todo o genoma. Se for causado por ancestrais estranhos em parte do genoma, não confiaríamos tanto na existência de introgressão significativa.

Para testar um sinal genômico consistente, usamos o procedimento block-jackknife.

• Block Jackknife

Embora os loci não sejam independentes, o Jackknife nos permite calcular a variância de D. Um método mais tradicional, ou seja, reamostrarmos aleatoriamente locos e recalculemos D, não é adequado porque loci próximos no genoma têm ancestrais semelhantes, tornando-os observações não independentes.

O procedimento block jackknife estima o desvio padrão dos chamados "pseudovalores" da D média do genoma inteiro, onde cada pseudovalor é calculado excluindo um bloco definido do genoma, tomando a diferença entre a D média do genoma inteiro e a D calculada quando o bloco é excluído.

Para considerar a não independência entre loci vinculados, o tamanho do bloco precisa exceder a distância na qual ocorre autocorrelação. Em nosso caso, usaremos um tamanho de bloco de 1 Mb, pois sabemos que o desequilíbrio de ligação decai para níveis de fundo a distâncias muito menores que 1 Mb.

O código para executar o processo jackknife é bastante simples, mas não o escreveremos aqui. Em vez disso, uma função R para este propósito é fornecida em um script separado, que agora podemos importar.

source("genomics_general-master/jackknife.R")

O primeiro passo do processo é definir os blocos que serão excluídos do genoma em cada iteração. A função get_block_indices no script jackknife fará isso e retornará os "índices" (ou seja, as linhas na tabela de frequência) correspondentes a cada bloco. Ele nos pede para especificar o tamanho do bloco para cada locus e os cromossomos e posições a serem analisados.

indices_bloco <- get.block.indices(tamanho_bloco=1e6,
                                   posicoes=tabela_frequencia$posicao,
                                   cromossomos=tabela_frequencia$scaffold)

num_blocos <- length(indices_bloco)

print(paste("Genoma dividido em", num_blocos, "blocos."))

Agora podemos executar o processo block jackknife para calcular a média e o erro padrão de D. Fornecemos a função estatistica_D criada anteriormente, que será aplicada em cada iteração. Também fornecemos as frequências de cada locus e os índices do bloco que serão usados para excluir todos os loci de um determinado bloco.

resultado_jackknife_D <- block.jackknife(indices_bloco=indices_bloco,
                                         FUN=estatistica_D,
                                         tabela_frequencia[,pop1], tabela_frequencia[,pop2], tabela_frequencia[,pop3])

print(paste("Média do jackknife de D =", round(resultado_jackknife_D$mean,4)))

Com base na média e na estimativa não viesada do erro padrão de D, podemos calcular a pontuação Z para testar se D se desvia significativamente de zero.

pontuacao_Z_D <- resultado_jackknife_D$mean / resultado_jackknife_D$standard_error

print(paste("Pontuação Z de D = ", round(pontuacao_Z_D,3)))

Normalmente, pontuações Z maiores que 3 ou 4 são consideradas significativas, portanto, uma pontuação Z grande neste caso significa que o desvio de zero é muito significativo.

• Estimação da proporção de introgressão

A estatística D fornece um poderoso teste para introgressão genética, mas não quantifica a proporção do genoma compartilhado. Já foi desenvolvido um método relacionado para estimar a proporção de "introgressão" (f).

A ideia por trás deste método é que cmopararemos o excesso observado de loci ABBA em relação a BABA com o esperado sob introgressão completa. Para aproximar o esperado sob introgressão completa, recalculamos ABBA e BABA, mas substituímos a segunda população da espécie P3 por P2. Se você não tiver uma segunda população, pode simplesmente dividir as amsotras de P3 em duas. Neste caso, temos duas populações para representar cada espécie, portanto, se usarmos H. cydno chioneus (cyd_chi) como P3a, podemos usar H. cydno zelinde (cyd_zel) como P3b.

Precisamos escrever nossa própria função para calcular f. A entrada será a frequência do alelo derivado em cada população, mas agora incluímos P3a e P3b.

estatistica_f <- function(freq_p1, freq_p2, freq_p3a, freq_p3b) {
    numerador_ABBA <- (1 - freq_p1) * freq_p2 * freq_p3a
    numerador_BABA <- freq_p1 * (1 - freq_p2) * freq_p3a

    denominador_ABBA <- (1 - freq_p1) * freq_p3b * freq_p3a
    denominador_BABA <- freq_p1 * (1 - freq_p3b) * freq_p3a

    (sum(numerador_ABBA, na.rm=TRUE) - sum(numerador_BABA, na.rm=TRUE)) /
    (sum(denominador_ABBA, na.rm=TRUE) - sum(denominador_BABA, na.rm=TRUE))
    }

Agora podemos escolher P3a e P3b e estimar f.

pop3a <- "cyd_chi"
pop3b <- "cyd_zel"

f <- estatistica_f(tabela_frequencia[,pop1], tabela_frequencia[,pop2], tabela_frequencia[,pop3a], tabela_frequencia[,pop3b])

print(paste("Proporção de introgressão = ", round(f,4)))

Isso indica que mais de 25% do genoma simpátrico de H. melpomene e H. cydno é compartilhado. A proporção de introgressão pode ser itnerpretada como a proporção média de linhagem estrangeira em qualquer genoma individual. Alternativamente, pode ser interpretada como a frequência esperada de alelos estrangeiros em qualquer locus dado no genoma dessa população.

Podemos novamente estimar o desvio padrão de f e obter um intervalo de confiança. Os índices do bloco jackknife já foram calculados, portanto, podemos simplesmente executar novamente a função Jackknife, desta vez especificando a função f como a função a ser executada em cada iteração.

resultado_jackknife_f <- block.jackknife(indices_bloco=indices_bloco,
                                         FUN=estatistica_f,
                                         tabela_frequencia[,pop1], tabela_frequencia[,pop2], tabela_frequencia[,pop3a], tabela_frequencia[,pop3b])

O intervalo de confiança de 95% é a média +/- ~1.96 erro padrão.

limite_inferior_IC_f <- resultado_jackknife_f$mean - 1.96*resultado_jackknife_f$standard_error
limite_superior_IC_f <- resultado_jackknife_f$mean + 1.96*resultado_jackknife_f$standard_error

print(paste("Intervalo de confiança de 95% de f =", round(limite_inferior_IC_f,4), round(limite_superior_IC_f,4)))

Análise ABBA-BABA por Cromossomo

• Todos os cromossomos mostram evidências de introgressão?

Acima, estudamos o grau de introgressão em todo o genoma. Supondo que cada cromossomo tenha número suficiente de SNPs, podemos realizar uma análise semelhante ao nível do cromossomo para avaliar a introgressão em cromossomos individuais.

O primeiro passo para executar esta tarefa é identificar as linhas na tabela de frequência que correspondem a cada um dos 21 cromossomos de Heliconius.

Primeiro, usamos a função unique para identificar todos os nomes de cromossomo existentes no conjunto de dados. Em seguida, precisamos identificar as linhas da tabela que correspondem a cada cromossomo. Para isso, usamos a função lapply, que aplica uma função simples várias vezes para criar uma saída em formato de tabela R. Neste caso, usaremos o nome do cromossomo para aplicar a função, e a função que aplicaremos simplesmente perguntará quais valores na coluna scaffold da tabela correspondem a esse cromossomo, utilizando a função which do R.

nomes_cromossomos <- unique(tabela_frequencia$scaffold)
indices_cromossomos <- lapply(nomes_cromossomos, function(crom) which(tabela_frequencia$scaffold == crom))
names(indices_cromossomos) <- nomes_cromossomos

Isso criará uma lista com 21 elementos - um para cada cromossomo. Cada elemento é um vetor de todos os loci da tabela que vêm desse cromossomo. Podemos verificar quantos SNP cada cromossomo tem aplicando a função length à lista que acabamos de criar.

sapply(indices_cromossomos, length)

(sapply é semelhante a lapply, mas tenta simplificar a saída, portanto, aqui retorna um vetor em vez de uma lista de vetores).

Agora podemos usar esses índices para calcular o valor D para cada cromossomo. Usamos sapply novamente, aplicando a função estatistica_D e indexando apenas as linhas da tabela de um cromossomo específico para cada caso.

valores_D_por_cromossomo <- sapply(nomes_cromossomos,
                                    function(crom) estatistica_D(tabela_frequencia[indices_cromossomos[[crom]], pop1],
                                                               tabela_frequencia[indices_cromossomos[[crom]], pop2],
                                                               tabela_frequencia[indices_cromossomos[[crom]], pop3]))

Também precisamos determinar se o D de cada cromossomo é significativamente diferente de zero. Primeiro, definiremos os blocos para cada cromossomo.

indices_bloco_por_cromossomo <- sapply(nomes_cromossomos,
                                        function(crom) get.block.indices(tamanho_bloco=1e6,
                                                                        posicoes=tabela_frequencia$posicao[tabela_frequencia$scaffold==crom]),
                                                                        simplify=FALSE)

Este comando retorna uma lista de listas. É uma lista com 21 elementos - um para cada cromossomo. Cada elemento é uma lista que fornece o índice de cada bloco dentro desse cromossomo.

Podemos verificar o número de blocos por cromossomo e o número de SNP por bloco para cada cromossomo.

sapply(indices_bloco_por_cromossomo, length)

lapply(indices_bloco_por_cromossomo, sapply, length)

Agora usamos o jackknife para calcular a pontuação Z para D de cada cromossomo.

resultados_jackknife_D_por_cromossomo <- sapply(nomes_cromossomos,
                                                function(crom) block.jackknife(indices_bloco=indices_bloco_por_cromossomo[[crom]],
                                                                               FUN=estatistica_D,
                                                                               tabela_frequencia[indices_cromossomos[[crom]], pop1],
                                                                               tabela_frequencia[indices_cromossomos[[crom]], pop2],
                                                                               tabela_frequencia[indices_cromossomos[[crom]], pop3]))

resultados_jackknife_D_por_cromossomo <- as.data.frame(t(resultados_jackknife_D_por_cromossomo))

resultados_jackknife_D_por_cromossomo$Z <- as.numeric(resultados_jackknife_D_por_cromossomo$mean) / as.numeric(resultados_jackknife_D_por_cromossomo$standard_error)

resultados_jackknife_D_por_cromossomo

Vemos que os cromossomos 1-20 mostram todas evidências significativas de introgressão (Z > 4), enquanto o cromossomo 21 (cromossomo Z) não. Na verdade, o D do chr21 é negativo, indicando que a população alopátrica de H. melpomene tem mais variação com H. cydno do que a população simpátrica de H. melpomene com H. cydno, embora a diferença não seja significativa. Isso indica uma redução significativa da introgressão no cromossomo sexual em comparação com o resto do genoma, consistente com seleção forte contra alelos de introgressão em loci no cromossomo Z. Se existirem uma ou mais incompatibilidades que afetem loci no cromossomo Z, é isso que esperaríamos.