Cálculo da probabilidade de sobreposição de lista de genes entre um conjunto de dados RNA seq e um chip ChIP

Espero que alguém nesses fóruns possa me ajudar com esse problema básico nos estudos de expressão gênica.

Fiz sequenciamento profundo de um tecido experimental e de controle. Em seguida, obtive valores de enriquecimento dobrado de genes na amostra experimental sob controle. O genoma de referência possui ~ 15.000 genes. 3.000 de 15.000 genes são enriquecidos acima de um certo ponto de corte na minha amostra de interesse em comparação ao controle.

Então: A = população total de genes = 15.000 B = subpopulação enriquecida por RNA-Seq = 3.000.

Em um experimento anterior com chip ChIP, encontrei 400 genes enriquecidos pelo chip ChIP. Dos 400 genes do chip ChIP, 100 genes estão no grupo de 3.000 transcritos de RNA-Seq enriquecidos.

Então: C = número total de genes enriquecidos em chip ChIP = 400.

Qual é a probabilidade de que meus 100 genes de chips ChIP sejam enriquecidos por RNA-Seq apenas por acaso? Em outras palavras, qual é a maneira mais prudente de calcular se minha sobreposição observada entre B e C (100 genes) é melhor do que a obtida apenas por acaso? Pelo que li até agora, a melhor maneira de testar isso é usando a distribuição hipergeométrica.

Usei uma calculadora on-line (stattrek.com) para configurar um teste de distribuição hipergeométrica com os seguintes parâmetros: - tamanho do pop = 15.000 - número de sucessos na população = 3.000 - tamanho da amostra = 400, - número de sucessos na amostra = 100. Eu obtenho o seguinte para Probabilidade hipergeométrica P (x = 100) = 0,00224050636447747

O número real de genes que se sobrepõem entre B e C = 100. Isso é melhor do que por acaso? Não parece que a chance de um gene ser enriquecido seja de 1: 5 (3.000 em 15.000). É por isso que não entendo como meu P (x = 100) calculado acima é 0,0022. Isso equivale a uma chance de 0,2% da sobreposição ocorrer por acaso. Isso não deveria ser muito maior?

Se eu amostrasse 400 genes aleatórios da grande lista de 15.000, seria esperado que 80 desses genes fossem enriquecidos apenas pelo acaso (1: 5). O número de genes que realmente se sobrepõem é 100, então isso é apenas um pouco melhor do que por acaso.

Eu também tentei encontrar uma solução usando as funções dhyper ou phyper em R (usando o que vi em outro post): A = todos os genes do genoma (15.000) B = genes enriquecidos em RNA-Seq (3.000) C = ChIP genes enriquecidos por chip (400) Aqui está a entrada / saída R (adaptada de uma postagem anterior de stackexchange):

> totalpop <- 15000    
> sample1 <- 3000    
> sample2 <- 400    
> dhyper(0:2, sample1, totalpop-sample1, sample2)    
[1] 4.431784e-40 4.584209e-38 2.364018e-36    
> phyper(-1:2, sample1, totalpop-sample1, sample2)    
[1] 0.000000e+00 4.431784e-40 4.628526e-38 2.410304e-36    

Não sei como interpretar esses números. Acredito que 2.36e-36 é a probabilidade de obter uma sobreposição completa entre B e C apenas por acaso? Mas isso não faz sentido, pois essa probabilidade está muito mais próxima de 1: 5. Se eu começar com 15.000 genes, 3.000 serão enriquecidos. Da mesma forma, se eu começar com 400 genes de chip ChIP, 80 deles devem ser enriquecidos apenas no RNA-Seq, devido às chances de enriquecimento de 1: 5 nesse conjunto de dados.

Qual é a maneira correta de calcular o valor-p, de acordo com a distribuição hipergeométrica, para a sobreposição de B e C?

Você está próximo, com o uso de dhypere phyper, mas não entendo de onde 0:2e de onde -1:2viemos.

O valor p que você deseja é a probabilidade de obter 100 ou mais bolas brancas em uma amostra do tamanho 400 de uma urna com 3000 bolas brancas e 12000 bolas pretas. Aqui estão quatro maneiras de calculá-lo.

sum(dhyper(100:400, 3000, 12000, 400))
1 - sum(dhyper(0:99, 3000, 12000, 400))
phyper(99, 3000, 12000, 400, lower.tail=FALSE)
1-phyper(99, 3000, 12000, 400)

Estes dão 0,0078.

dhyper(x, m, n, k)dá a probabilidade de desenhar exatamente x. Na primeira linha, resumimos as probabilidades para 100 - 400; na segunda linha, tomamos 1 menos a soma das probabilidades de 0 - 99.

phyper(x, m, n, k)dá a probabilidade de obter xou menos, phyper(x, m, n, k)o mesmo acontece com sum(dhyper(0:x, m, n, k)).

O lower.tail=FALSEé um pouco confuso. phyper(x, m, n, k, lower.tail=FALSE)é o mesmo que 1-phyper(x, m, n, k), e também é a probabilidade de x+1ou mais. [Eu nunca me lembro disso e sempre tenho que checar.]

$\ge$

max(dhyper(0:400, 3000, 12000, 400))$\sim$

Aqui está uma imagem da distribuição hipergeométrica neste caso. Você pode ver que ele está centrado em 80 (20% de 400) e que 100 está bem distante na cauda direita.

Veja-o desta maneira. Se você o assumiu como um binômio, o que pode não estar correto, mas deve ser razoavelmente aproximado. Seu sigma ^ 2 é 0,8 * .2 * 400 = 64, então sigma = 8. Então de 80 para 100, você fez 2,5 desvios-padrão. Isso é bastante significativo. Deveria ter um pequeno valor-p.