Enriquecimento com teste hipergeométrico (Hypergeometric test) - lmigueel/Bioinformatica GitHub Wiki
1. Sobre
O propósito de um teste de enriquecimento de genes, por exemplo os diferenciais entre duas condições, é você identificar quais genes estão super-representados a partir de uma lista de genes de entrada, levando em conta um base de dados conhecida (background). Mas isso podemos expandir para qualquer tipo de análise, inclusive para Metagenômica, onde você possui um background com a taxonomia de todas as suas OTUs, e deseja descobrir se uma classe delas está enriquecida.
Uma maneira de realizar este enriquecimento é através do Teste Hipergeométrico. Em um primeiro momento, modelamos a associação entre genes e classe GO usando uma distribuição hipergeométrica. O exemplo clássico para o hipergeométrico é a seleção aleatória de bolas k em uma urna contendo m bolas "marcadas" e n "não marcadas", e a observação de uma amostra que contém x bolas marcadas.
2. Pacote para os cálculos
Para realizar os cálculos iremos utilizar a função phyper do pacote stats no R. O pacote stats já é default e não precisa instalação prévia.
phyper(q, m, n, k, lower.tail = TRUE, log.p = FALSE)
| Argumento | explicação |
|---|---|
| q | vetor de quantis que representa o número de bolas brancas tiradas sem substituição de uma urna que contém bolas pretas e brancas. |
| m | o número de bolas brancas na urna. |
| n | o número de bolas pretas na urna. |
| k | o número de bolas retiradas da urna. |
| lower.tail | lógico; se TRUE (padrão), as probabilidades são P[X ≤ x], caso contrário, P [X> x]. |
Quando definimos o parâmetro lower.tail = FALSE no phyper, a interpretação do valor p é P [X> x]. Mas o que precisamos testar é a hipótese nula P [X ≥ x], então subtraímos x por 1. Quando o valor p <= 0,05, podemos rejeitar a hipótese nula e assumimos um enriquecimento significativo, ou seja, há um pequena probabilidade de ver uma sobreposição igual ou maior do que x entre a sua lista e o total de genes/objetos "marcados".
3. Exemplo 1
Para ilustrar a maneira de realizar o teste hipergeométrico, vamos pensar que desejamos verificar se a classe funcional (cell cycle) (GO:0007049) está enriquecida ou não na minha lista de genes. Supomos que existam 13588 genes com algum processo GO (background). Observando o processo GO:0007049, notamos que ele possui 611 genes relacionados, ou seja, 12977 genes não estão relacionados a este processo. Identificamos 59 genes diferenciais em uma comparação, sendo que 19 deles estão relacionados ao processo GO:0007049 (cell cycle). Será que essa classe está enriquecida na minha lista? Qual o p-value?
Definimos os parâmetros do teste hipergeométrico da seguinte maneira:
| Variável | Descrição |
|---|---|
| N = 13588 | número total de genes com alguma anotação GO (background). |
| m = 611 | número de elementos "marcados", ou seja, o número total de genes anotados para o processo cell cycle (GO:0007049). |
| n = N − m | número de elementos "não marcados", ou seja, o número de genes que têm alguma anotação GO, mas não estão associados ao termo GO selecionado (cell cycle - GO:0007049). |
| k = 59 | Tamanho da sua entrada, ou seja, o número de genes da entrada que possuem alguma anotação e estão associados a pelo menos um "Processo Biológico" GO. |
| x = 19 | número de elementos "marcados" na seleção, ou seja, número de genes da sua lista de entrada E associados a classe funcional cell cycle. |
Utilizando a função phyper :
N = 13588
m = 611
n = N-m
k=59
x=19
phyper(q=x-1, m=m, n=n, k=k, lower.tail=FALSE)
Ela irá retornar o p-value:
[1] 4.989683e-12
ou seja, a classe funcional (cell cycle) (GO:0007049) está enriquecida, segundo a minha lista de entrada com p-value = 4.989683e-12.
Importante mencionar que para ajustar seus p-values no R, você pode usar a função p.adjust(p), onde p é o vetor de p-values. Como estamos calculando para várias classes de GO, os softwares de enriquecimento sempre ajustam os valores de p. Sendo assim, caso tenha executado o cálculo com o phyper algumas vezes para algumas classes, utilize:
p.adjust(p, method="bonferroni") # método de Bonferroni (multiplicação simples)
p.adjust(p, method="holm") # método de Holm-Bonferroni (mais potente que o Bonferroni)
p.adjust(p, method="BH") # Benjamini-Hochberg
4. Exemplo 2
Vamos supor agora que estamos lidando com genes da classe funcional de resposta imune (immune response), e possuímos uma lista de 100 genes. Será que os genes de resposta imune estão super-representados em nossa lista de 100 genes? Vamos fazer uma tabelinha rápida e usar o phyper.
| - | Total | genes em resposta imune | genes em não resposta imune |
|---|---|---|---|
| Background | 30.000 | 600 | 29400 |
| Minha lista | 100 | 10 | 90 |
> phyper(9, 600, 29400, 100, lower.tail= FALSE)
[1] 3.276244e-05
ou até utilizando a matriz para o teste de Fisher:
> fisher.test(matrix(c(10,590,90, 29310), nrow=2),alternative="greater")$p.value
[1] 3.276244e-05
ou seja, a classe funcional immune response está enriquecida.
Citação
https://stat.ethz.ch/R-manual/R-devel/library/stats/html/Hypergeometric.html