Analyse d'enrichissement par niveau de duplication de gènes

Contexte biologique

Au fil du temps, certaines espèces végétales ont tendance à dupliquer leurs génomes entiers, obtenant une copie supplémentaire de chaque gène. En raison de l'instabilité de cette configuration, bon nombre de ces gènes sont ensuite supprimés, et le génome se réorganise et se stabilise, prêt à être dupliqué à nouveau. Ces événements de duplication sont associés à des événements de spéciation et d'invasion, et la théorie est que la duplication aide les plantes à s'adapter plus rapidement à leurs nouveaux environnements.

Lupinus, un genre de plante à fleurs, a envahi les Andes dans l'un des événements de spéciation les plus rapides jamais détectés, et de plus, il semble avoir plus de copies en double dans son génome que le genre le plus proche, Baptisia.

Et maintenant le problème mathématique:

Les génomes d'un membre de Lupinus et d'un membre de Baptisia ont été séquencés, fournissant des données brutes sur environ 25 000 gènes dans chaque espèce. En interrogeant une base de données de gènes de fonction connue, j'ai maintenant une "meilleure estimation" des fonctions que ce gène pourrait remplir - ainsi, par exemple, Gene1298 pourrait être associé au "métabolisme du fructose, réponse au stress salin, réponse au stress froid". Je veux savoir, s'il y a eu un événement de duplication entre Baptisia et Lupinus, si la perte de gènes a eu lieu au hasard, ou si les gènes remplissant des fonctions particulières étaient plus susceptibles d'être conservés ou supprimés.

J'ai un script qui produira un tableau comme celui illustré ci-dessous. L * est un décompte de tous les gènes de Lupinus associés à la fonction. L 1+ est un nombre de gènes de lupinus associés à la fonction où il existe au moins une copie en double. Je peux le faire produire L 2+, L 3+ etc., bien que L 1+ soit un groupe beaucoup plus fiable que L 2+ en raison du processus de séquençage.

Function            | L *  | L 1+ | B *  | B 1+ |
fructose metabolism | 1000 | 994  | 1290 | 876  |
salt stress         | 56   | 45   | 90   | 54   |
etc.

Ce que je voudrais faire, c'est tester, pour chaque fonction génique, s'il y a plus ou moins de gènes avec des doublons que ce à quoi on pourrait s'attendre par pur hasard dans Lupinus et Baptisia, et si Lupinus diffère de Baptisia dans le rapport observé / attendu.

La meilleure chose que j'ai à ce jour

Des études antérieures sur différentes espèces ont utilisé l'analyse d'enrichissement, avec le test exact de Fisher et la correction FDR pour un échantillonnage multiple, pour effectuer un test de contingence sur chaque ligne.

Ce serait bien de l'améliorer; Je ne suis pas sûr que cela semble être la meilleure façon de le faire.

Glen_b a suggéré d'utiliser un GLM pour analyser les données; J'ai joué avec les GLM en JMP8, ce qui a été intéressant, mais j'avoue ne pas vraiment les comprendre.

Cela dit, j'essaie d'utiliser R à la place maintenant.

Pourquoi est-ce que j'utilise cela?

À l'origine, cela devait faire partie d'un court projet de recherche que je fais à l'université, mais s'est maintenant étendu à un énorme projet d'annotation du génome. Pourquoi? Parce que la bioinformatique est cool. Être capable de prendre une chaîne de A, T, C et G et de l'utiliser pour déduire des informations sur des événements qui se sont produits il y a des millions d'années est incroyable.

Inutile de dire que je ne vais pas essayer de soumettre une réponse aimablement fournie comme mon propre travail. Je serais heureux d'inclure une reconnaissance dans le document si j'utilise une méthode suggérée ici dans le travail soumis.

Bien que je convienne que le test de Fisher (ou quelque chose de similaire) peut être l'approche la plus naturelle ici, qu'en est-il de ceci:

1. Pour chaque gène unique, vous déterminez la différence de nombre de duplications en L et B
2. Ordonner les gènes par cette différence. Maintenant, les gènes montrant la plupart des différences entre les espèces seront en haut de votre liste.
3. Appliquer un test d'enrichissement de l'ensemble de gènes à la liste ordonnée de gènes. Par exemple, vous pouvez utiliser une méthode de Fisher modifiée de mon package tmod , pour laquelle vous devrez définir vos ensembles de gènes (cela devrait être assez simple). Notez que la méthode de Fisher n'est pas liée au test de Fisher.

Le test de Fisher modifié (surnommé CERNO par les auteurs qui l'ont décrit pour la première fois dans ce contexte) prend en entrée toute liste ordonnée de gènes, tant que vous pouvez les regrouper en certaines catégories utiles.

L'avantage de cette approche est qu'en dehors d'une valeur de p, vous pouvez facilement calculer la taille d'effet de l'enrichissement et la visualiser (par exemple, sous la forme d'une courbe ROC sur la liste ordonnée de gènes). Cela vous donne une bien meilleure idée de l'importance de ce que vous observez pour la biologie que vous étudiez.

Comme vous le dites, vous posez deux questions distinctes.

Question 1 "le rapport L * / L1 + est-il différent de B * / B1 + pour une fonction génétique donnée"

cela pourrait être mieux répondu par le test exact de Fisher en utilisant des données sur la ligne comme vous l'avez trouvé précédemment.

Question 2 "est le rapport: gènes où il y a une seule copie / gènes où il y a plus d'une copie, différent entre les fonctions des gènes?"

Je pense que cela pourrait aussi être mieux répondu par le test exact de Fisher. Vous testeriez le rapport L * / L1 + pour la fonction génique 1 contre L * / L1 + pour la fonction génique 2. Ensuite, la fonction génique 1 contre la fonction génique 3, etc.

Aucune de ces séries de questions ne permet de savoir si elles sont maintenues / supprimées plus rapidement que prévu par pur hasard, seulement si elles sont supprimées / maintenues à des taux différents les uns des autres. Pour savoir si elles étaient supprimées / maintenues à un taux différent de celui du hasard, vous devez connaître le rapport des copies uniques / des copies multiples pour de nombreuses régions d'ADN qui ne sont affectées que par hasard. Si vous pouviez trouver de telles régions, vous vous retrouveriez avec un "groupe de fonctions" où la fonction est "aucune". Vous devez ensuite comparer cela à vos autres groupes de fonctions génétiques de la même manière que je l'ai décrit à la question 2.