Informations

Quelles sont les principales raisons de l'échec de la localisation/d'ancrage des séquences dans les assemblages de génomes ?


Ma question concerne l'incorporation de lectures de séquences individuelles dans les chromosomes lors de projets de séquençage de gènes, en particulier ceux avec des génomes plus gros tels que Drosophila melanogaster ou Homo sapiens.

Il existe certaines séquences ou contigs courts qu'il n'a pas encore été possible d'incorporer dans les assemblages chromosomiques plus importants.

Quelles sont les raisons de cet échec ?


Je vais réinterpréter légèrement cette question pour faciliter la réponse :

« Pourquoi les assemblages de génomes se composent-ils fréquemment d'un grand nombre de contigs courts plutôt que d'un nombre relativement petit de longs chromosomes (ou de réplicons complets d'autres types) ? Et comment pourrais-je améliorer mon assemblage ? »

Les assemblages de séquences de novo (généralement) consistent en un ensemble de fragments de séquences contigus (« contigs ») qui ont été dérivés du séquençage par une certaine heuristique. Si l'assemblage de séquences faisait un travail parfait, nous nous attendrions à ce que ces contigs soient des chromosomes complets. Au lieu de cela, ce que nous observons, ce sont de nombreux fragments de séquences avec des relations incertaines les uns avec les autres. Dans ce cas courant, nous pensons que certains des contigs pourraient être cousus ensemble ("localisés", "ancrés", "échafaudés") de manière significative en chromosomes, mais nous ne savons pas comment.

Pour l'assemblage de séquences de novo, il n'y a pas de vérité fondamentale à laquelle se comparer. Ainsi, toutes les informations pour générer des assemblages proviennent des séquences de départ (qui sont généralement des "lectures", de petits lambeaux de séquence qui sortent d'un instrument de séquençage). Les heuristiques qui génèrent des contigs à partir des lectures consistent à trouver des chevauchements entre différentes lectures de départ et à essayer de les étendre dans les plus grands groupes de lectures partiellement chevauchants possibles, en utilisant certaines hypothèses sur les taux d'erreur dans la séquence des lectures (et d'autres variables). Voir wikipedia ou ce document (encadré de la figure) pour un peu plus d'informations à ce sujet.

Il s'avère qu'il existe un corps de théorie statistique qui régit l'ampleur de ces chevauchements, et donc la taille/complétude de vos contigs. Ce corps de théorie dit que la taille et le nombre de contigs sont prévisibles dans une certaine mesure à partir de la taille du génome que vous essayez d'assembler, de la longueur des lectures et du nombre de lectures que vous avez.

Certains génomes, comme le maïs, le blé ou le génome humain, sont également plus difficiles à assembler car ils sont « répétitifs ». Cela signifie qu'il existe de nombreuses séquences très similaires à différents endroits du génome. Ainsi, lorsque vous lisez l'un de ces endroits, vous ne savez pas de quelle région très similaire il provient. Ces régions répétitives sont difficiles à résoudre à moins que vous n'ayez de très longues lectures ou des informations parallèles, elles entraînent donc également des lacunes dans la séquence. De telles régions auront également tendance à être "effondrées", ce qui signifie que vous n'avez qu'un seul contig représentant une répétition qui devrait en fait être présent sous forme de plusieurs contigs avec des emplacements différents dans l'assemblage.

Les sources d'informations parallèles utiles sont la ligature de proximité (méthodes telles que Hi-C), les cartes génétiques issues de croisements ou la cartographie optique. Ces sources d'information orthogonales expriment les relations entre contigs en fonction de leur proximité sur un chromosome. Voici un article avec plus d'informations sur ces sujets. Récemment, les technologies de lecture longue telles que PacBio ou le séquençage des nanopores ont commencé à lire à travers de grandes répétitions, mais elles ne sont généralement pas encore assez fiables pour assembler complètement même des génomes plus petits "télomère à télomère".

Donc, s'il existe une recette pour localiser/ancrer/échafauder des assemblages fragmentés, voici en gros ce qu'il est possible de faire :

  1. obtenir plus de lectures
  2. obtenir des lectures plus longues
  3. obtenir des informations orthogonales sur le placement des contigs (par exemple, ligature de proximité, cartes génétiques, assemblage étroitement lié, etc.)

J'espère que cela pourra aider.


L'assemblage et l'analyse de lectures de séquences de génomes non cartographiées révèlent une nouvelle séquence et une variation chez les chiens

Les chiens sont d'excellents modèles animaux pour les maladies humaines. Ils ont des antécédents vétérinaires étendus, des pedigrees et un système génétique unique en raison des pratiques d'élevage. Malgré ces avantages, un facteur limitant leur utilité est la référence du génome canin (CGR) qui a été assemblée à l'aide d'un seul Boxer de race pure. Bien qu'il s'agisse d'une pratique courante, de nombreuses lectures de haute qualité restent non mappées. Pour traiter ces données de séquence du génome entier de trois races, Border Collie (n = 26), Bearded Collie (n = 7) et Entlebucher Sennenhund (n = 8), ont été analysés pour identifier de nouveaux contigs génomiques non-CGR en utilisant le pseudo- préalablement validéde novo canalisation d'assemblage. Nous avons identifié 256 957 nouveaux contigs et relations appariées ainsi que les scores BLAT ont fourni 126 555 (49 %) contigs de haute qualité avec des coordonnées génomiques contenant 4,6 Mb de nouvelle séquence absente du CGR. Ces contigs comblent 12 503 lacunes connues, dont 2,4 Mb contenant des séquences partiellement manquantes pour 11,5 % d'Ensembl, 16,4 % de RefSeq et 12,2 % des gènes annotés canFam3.1+ CGR et 1 748 contigs non cartographiés contenant 2 366 nouvelles variantes de gènes. Exemples de six gènes associés à la maladie (ÉCHARPE2, RD3, COL9A3, FAM161A, RASGRP1 et DLX6) contenant des lacunes ou des variantes d'épissage alternatif manquantes dans le CGR sont également présentées. Ces résultats provenant de races non-référencées soutiennent la nécessité d'améliorer le CGR actuel réservé aux Boxers pour éviter de manquer des informations biologiques importantes. L'inclusion des séquences de gènes manquantes dans le CGR facilitera l'identification des mutations pathologiques putatives dans diverses races et phénotypes.


Introduction

Hevea brasiliensis (Willd.) Müll.-Arg. (l'hévéa Para) est une culture arborescente pérenne à feuilles caduques indigène des forêts tropicales du bassin amazonien en Amérique du Sud. C'est une espèce monoïque, allogame appartenant à la famille des Euphorbiaceae et a un nombre de chromosomes de 2m = 2X = 36. L'hévéa produit des polymères isoprénoïdes de haute qualité (cis-1,4-polyisoprène) avec des propriétés physiques uniques telles que l'élasticité, la résilience et une dispersion efficace de la chaleur qui sont inégalées par les substituts synthétiques à base de pétrole 1 . Sur les 2500 espèces végétales productrices de latex, H. brasiliensis est la seule espèce qui produit des quantités commercialement viables de caoutchouc naturel de haute qualité, représentant plus de 98 % de la production totale dans le monde 2,3 . Les pays d'Asie du Sud-Est dominent actuellement la production et le commerce du caoutchouc, représentant plus de 76 % des 11,96 tonnes produites dans le monde, la Thaïlande et l'Indonésie étant les deux plus grands producteurs et exportateurs de caoutchouc naturel (http://faostat3.fao.org/ Consulté en septembre 14, 2016). Bien qu'elle soit originaire du bassin amazonien, la production de caoutchouc en Amérique du Sud ne représente que 2% de la production totale dans le monde en raison de la propagation dévastatrice de la brûlure des feuilles en Amérique du Sud (SALB) causée par l'ascomycète. Microcyclus ulei dans les années 1930 4 . L'infestation a entraîné l'effondrement de la plupart des plantations d'hévéas et l'arrêt de la production commerciale de caoutchouc au Brésil et dans d'autres pays d'Amérique du Sud 5 . Bien que Hévéa la culture en Asie du Sud-Est n'a pas été affectée par le SALB, d'autres champignons pathogènes indigènes sont toujours des menaces majeures pour la production de caoutchouc dans cette région. Un certain nombre de clones à haut rendement cultivés commercialement sont sensibles à la maladie anormale de la chute des feuilles causée par divers Phytophtora espèces et jusqu'à 40 % de perte de rendement a été observée lors de l'infection de clones hautement sensibles 6 . Un autre agent pathogène répandu, Corynespora cassiicola, est un champignon nécrotrophe qui cause la maladie de la chute des feuilles à Corynespora, entraînant des pertes significatives de rendement en caoutchouc naturel 7 .

Au cours des deux dernières décennies, les phytogénéticiens ont utilisé des approches conventionnelles de sélections récurrentes pour obtenir des clones avec un rendement de latex amélioré et une résistance aux agents pathogènes fongiques 8 . Cependant, la sélection traditionnelle est intrinsèquement longue en raison du long cycle de sélection de cette espèce. La capacité de sélectionner des caractères commercialement valables à un stade précoce aura un impact considérable sur la réduction du temps et des ressources nécessaires pour développer des clones supérieurs. L'avènement des marqueurs génétiques basés sur l'ADN a fourni une stratégie alternative qui a permis aux sélectionneurs de plantes d'accélérer le processus de développement des cultivars. Les progrès récents de la technologie de séquençage ont simplifié et accéléré la découverte de variantes de séquence, permettant le passage de marqueurs anonymes tels que les microsatellites à des polymorphismes nucléotidiques uniques (SNP) 9,10 plus omniprésents. La disponibilité de marqueurs SNP facilite la construction de cartes de liaison génétique à haute densité, qui sont déterminantes pour les analyses de loci de caractères quantitatifs (QTL) et les études d'association avec des caractères agronomiques précieux 11 .

Plusieurs ressources génomiques pour H. brasiliensis ont été développés récemment, y compris le séquençage du transcriptome de divers tissus et le séquençage du génome RRIM600 cultivé commercialement 10,12,13. Alors que la première ébauche de la séquence du génome 13 a fourni une excellente source d'informations génomiques, l'assemblage est très fragmenté, contenant plus d'un million de contigs avec une longueur N50 de 2,9 kb. L'assemblage d'un génome complexe est difficile en partie en raison de la présence de séquences d'ADN hautement répétitives, qui introduisent une ambiguïté lors de la reconstruction du génome. Il a été estimé que l'ADN répétitif représente 78% du génome de l'hévéa 13 , ce qui pose une difficulté majeure dans le de novo assemblage en particulier lorsqu'il est construit exclusivement à partir de données à lecture courte.

Bien qu'il soit possible de construire un assemblage de haute qualité pour les espèces végétales avec de grands génomes grâce au séquençage hiérarchique, l'utilisation d'une plate-forme de séquençage capable de générer de longues lectures pouvant couvrir des régions de répétitions complexes est susceptible d'être plus rentable et moins à forte intensité de main-d'œuvre/de temps. Le développement récent et les applications des technologies de séquençage à lecture longue ont montré une amélioration substantielle des assemblages de génomes, principalement en joignant les contigs et les échafaudages et en couvrant les espaces autour des régions répétitives 14,15. La technologie de séquençage Pacific Biosciences (PacBio) offre des lectures de la taille d'un kilobase sans biais GC ni erreurs systématiques. Cependant, les données brutes générées à partir de ce système sont intrinsèquement sujettes aux erreurs, avec une précision moyenne de

82 % 16 . La disponibilité d'un logiciel d'assemblage hybride permet aux biologistes de combiner les données de lecture longue de PacBio avec des données de lecture courte complémentaires et de haute précision provenant d'autres plates-formes 17 . Nowak et al. 18 a démontré que l'incorporation d'une couverture supplémentaire de 8 × de données PacBio à l'assemblage existant en lecture seule du Primula veris génome a entraîné la fermeture de 21 % des lacunes et le comblement de 38 % des positions ambiguës dans les lacunes, en plus de l'amélioration de 40 % de la taille du contig N50 18 . De même, un ensemble de longues lectures 16,5 × PacBio a été utilisé pour combler 68% des lacunes existant dans l'assemblage Illumina uniquement du cichlidé africain 14 .

Récemment, le séquençage à lecture courte à haut débit a été utilisé pour explorer la liaison à l'ADN jusqu'à plusieurs centaines de kilobases dans des bibliothèques de ligature de proximité construites à partir de in vitro chromatine reconstituée 19 . Le propriétaire in vitro La méthode d'assemblage du génome à longue distance basée sur la capture de la conformation de la chromatine appelée "Chicago" développée par Dovetail Genomics a été utilisée en combinaison avec une approche standard de séquençage du génome entier pour générer un de novo assemblage du génome humain avec une précision et une contiguïté comparables aux assemblages construits à l'aide de méthodes plus coûteuses. La technique d'échafaudage à longue portée de Chicago a également été utilisée pour améliorer les assemblages de génomes existants chez l'alligator américain 19 , la grenouille à griffes africaine 20 et le manioc 21 .

Aller au-delà de l'assemblage fragmenté du génome nécessite une carte de liaison de haute qualité et à haute densité qui peut aider à ordonner et à orienter de novo échafaudages assemblés en séquences à l'échelle pseudo-chromosomique 22 . Des informations complémentaires dans un de novo un assemblage de fusil de chasse et une carte génétique peuvent être intégrés pour générer une séquence de référence de haute qualité. Des cartes de liaison à haute densité ont été utilisées avec succès pour ancrer des échafaudages à partir d'assemblages de fusils de chasse à génome entier dans des espèces végétales allant d'un haricot commun diploïde 23 à un blé tendre hexaploïde 24 . L'intégration d'une carte génétique et d'un assemblage de génomes nous permet d'étudier les duplications de génomes entiers anciens et récents et d'effectuer une génomique comparative pour étudier l'architecture et l'évolution du génome à travers les espèces 25 .

Ici, nous avons utilisé des technologies complémentaires pour générer des lectures courtes 454/Illumina à couverture profonde et des données de lecture longue PacBio à couverture moyenne pour un de novo assemblage hybride du H. brasiliensis génome du clone BPM24 (Bank Pertanian Malaysia 24). Nous avons ensuite appliqué une technique d'assemblage de Chicago à longue portée 19 pour échafauder notre assemblage préliminaire afin d'obtenir 1,26 Go de séquences d'hévéas assemblées. L'accession BPM24 présente un degré élevé de résistance à deux agents pathogènes fongiques, Phytophthora et Corynespora, que l'on trouve couramment en Asie du Sud-Est 26 , et il est actuellement exploité comme source génétique de résistance fongique dans les programmes de sélection d'hévéas en Thaïlande. La disponibilité de la carte de liaison consensus haute densité construite à partir de deux populations dérivées de BPM24 11 nous a permis d'ancrer et d'orienter un grand nombre d'échafaudages dans ce nouvel assemblage. Ces séquences génétiquement ancrées ont fourni la première preuve concrète de la présence de paléotétraploïdie chez l'hévéa et nous ont permis de réaliser des analyses comparatives parmi les Euphorbiacées.


II. Complexité des génomes végétaux

Les premières cibles du séquençage du génome, en dehors du génome humain, étaient des génomes relativement petits et, par conséquent, plus faciles et moins coûteux à séquencer, tels que Arabidopsis– un cinq chromosomes, c. Plante à génome haploïde de 120 Mbp ( Meinke et al., 1998 ). En réalité, Arabidopsis est beaucoup plus petit que la plupart des génomes végétaux, par ordre de grandeur. Le génome végétal moyen est > 6000 Mbp par génome haploïde pour les angiospermes ( Gregory et al., 2007 ), environ deux fois la taille du génome humain. De nombreux génomes de plantes économiquement importants sont encore plus gros. Le blé, par exemple, est c. 15 Gbp par génome haploïde et le pin a au moins un génome de 26 Gbp ( Valkonen et al., 1993 ). La taille du génome est un contributeur à la complexité du génome végétal. D'autres contributeurs incluent la polyploïdie et les séquences d'ADN répétitives, et, en particulier, les éléments transposables. Ensemble, ces attributs des génomes végétaux augmentent le coût du séquençage et ont un impact négatif sur la qualité de la séquence résultante, d'autant plus que le champ migre du séquençage basé sur une carte (qui sera décrit plus tard) au séquençage à lecture courte du génome entier (WGS) .

Deux facteurs principaux qui contribuent à la taille et à la complexité du génome végétal sont la polypoïdie et les séquences d'ADN répétitives (exemples de cultures illustrés à la figure 1). La polyploïdie est l'accumulation d'ensembles supplémentaires de chromosomes par autopolyploïdie, doublement du même génome, ou allopolyploïdie, deux génomes divergents dans le même noyau. L'augmentation du nombre de chromosomes et de la teneur en ADN sont des conséquences immédiates de la polyploïdie, mais selon le moment où l'événement de polyploïdie s'est produit, l'augmentation du nombre de chromosomes peut ne pas être immédiatement apparente car les anciens événements de polyploïdie sont susceptibles d'être partagés par les taxons frères et/ou la diploïdisation du nombre de chromosomes peut avoir (réduction du nombre de chromosomes par perte et réarrangements). La plupart des plantes terrestres, sinon toutes, ont subi des événements de polyploïdie à divers moments de leur évolution (revue dans Soltis et al., 2004 ). Par exemple, le soja (Glycine max) a subi au moins trois événements polypoïdes qui, du fait d'une séquence génomique de haute qualité ( Schmutz et al., 2010 ), peuvent maintenant être examinés. Le premier événement, et le plus difficile, à détecter était un événement précoce dans l'évolution des plantes partagé par de nombreuses plantes terrestres ( Bowers et al., 2003 ). Le deuxième événement a été c. il y a 45 à 55 millions d'années (Mya) et devrait être partagé avec les légumineuses qui ont divergé après cet événement, telles que Medicago (Canon et al., 2006 ). L'événement le plus récent, c. 5 Ma, était très probablement un événement allopolyploïde ( Gill et al., 2009 ) qui a coïncidé avec l'émergence de la Glycine genre ( Innes et al., 2008 ). Ainsi, le génome du soja de 1,1 Gbp contient des reliques d'au moins trois événements de polyploïdie qui ont abouti à un génome qui est une mosaïque de segments dupliqués ( Schmutz et al., 2010 ). Dans le Glycine genre, cependant, il y a un événement allopolyploïde encore plus récent qui s'est produit dans les espèces vivaces trouvées en Australie ( Doyle et al., 2002 ). Ainsi, la polyploïdie est un processus récurrent qui façonne et façonne les génomes des plantes au cours de l'évolution.

Diagramme de la relation évolutive de plusieurs espèces de cultures majeures montrant des événements de polyploïdie (octogones rouges), la taille relative du génome par rapport au riz (taille des cercles) et le pourcentage d'éléments transposables (couleur du cercle). Le temps estimé de divergence et les événements polyploïdes sont indiqués sur la base de la littérature ( Gaut & Doebley, 1997 Huang et al., 2002 Blanc & Wolfe, 2004 Paterson et al., 2004, 2009 Schlüter et al., 2004 Swigonova et al., 2004 IRGSP, 2005 Jaillon et al., 2007 Zhu et al., 2008 Schnable et al., 2009 Osier et al., 2009 Abrouk et al., 2010 Choulet et al., 2010 Schmutz et al., 2010 ). Mya, il y a des millions d'années.

En plus de la polyploïdie, les séquences d'ADN répétitives et, en particulier, les éléments transposables (ET) composent de grandes fractions de la plupart des génomes végétaux et sont des obstacles à un séquençage efficace du génome. Les TE ont été revus en profondeur ( Bennetzen et al., 2005 ) ici, nous nous concentrerons uniquement sur la contribution à l'obésité génomique chez les plantes et à l'organisation des génomes végétaux car elle contribue à l'obtention de séquences génomiques précises. Il existe plusieurs cas d'amplification rapide de quelques familles TE qui ont entraîné une augmentation de la taille du génome. Oryza australiensis, par exemple, est environ le double de la taille de ses plus proches parents en raison de l'amplification de trois familles TE ( Piegu et al., 2006 ). Le maïs est l'exemple le plus frappant d'obésité génomique résultant de l'amplification TE ( SanMiguel et al., 1996 ).

Le facteur de complication de l'amplification TE sur le séquençage du génome n'est pas principalement l'augmentation de la quantité d'ADN à séquencer, mais plutôt l'effet de nombreuses copies de la même séquence dans tout le génome qui rendent la cartographie et l'assemblage difficiles.Si une famille TE s'est amplifiée récemment, elle peut avoir des milliers de copies dispersées dans tout le génome, toutes avec une identité de séquence très élevée. Si un génome est séquencé via une approche shotgun, comme la plupart des génomes le sont de nos jours, ces TE très similaires compliqueront l'assemblage à moins qu'il n'y ait suffisamment de lectures de paires de partenaires qui couvrent les répétitions (Fig. 2).

Schéma du séquençage du génome entier au fusil de chasse (WGS). Le segment chromosomique est représenté en bas en noir avec des éléments transposables (TE) en orange et des séquences de lecture courte au-dessus. Les lectures courtes sont idéalisées à la fois dans la couverture et dans l'arrangement. (a) L'assemblage de lectures courtes donne trois contigs de séquences non ordonnées, car les lectures des TE ne peuvent pas être placées sans ambiguïté. (b) En utilisant une variété de morceaux d'ADN plus gros (2 à 4, 6 à 10 et 20 Kpb illustrés) avec des séquences de paires de partenaires, les TE peuvent être étendus et des contigs ordonnés en résultent qui peuvent ensuite être placés dans un contexte chromosomique. , séquences d'ADN répétitives, > 98 % d'identité de séquence , séquences de fusil de chasse.


La prévention

L'échange de codes-barres observé sur le HiSeq 4000 était probablement le résultat de la nouvelle conception utilisant des cellules d'écoulement à motifs, et devrait donc probablement se produire également sur le HiSeq 3000 et le HiSeq X Ten et la machine phare NovaSeq. L'utilisation de Flowcells à motifs nécessite le besoin d'un nouveau type de chimie de séquençage, appelé amplification d'exclusion (ExAmp), remplaçant la génération de clusters par une amplification en pont utilisée sur la Flow Cell conventionnelle sans motifs.

Les détails techniques décrivant la nouvelle chimie de séquençage n'ont pas été rendus publics, mais l'examen des brevets associés a permis de comprendre le processus global et il est donc raisonnable de supposer que le saut d'index se produit avant la génération de cluster, puisque tous les réactifs nécessaires pour le cluster génération sont présents dans ce mélange réactionnel. Les amorces à index libres dans ce mélange ont le potentiel d'amorcer des fragments de bibliothèque et d'être prolongées par l'ADN polymérase. Ces molécules, dont les codes-barres sont incorrectement attribués, sont libres de générer des amas d'ADN. Ce n'est pas le cas pour la conception de cellule d'écoulement conventionnelle dans laquelle les codes-barres libres sont éliminés après l'hybridation de l'ADN à la cellule d'écoulement. Ce n'est qu'après cette étape que l'ADN polymérase est ajoutée et que l'extension de séquence peut commencer.

Il n'est pas clair à l'heure actuelle comment mettre en œuvre une solution bioinformatique à usage général pour supprimer de tels artefacts après le séquençage, et en effet, des solutions expérimentales ne sont peut-être possibles. Pour lutter contre le problème, lorsqu'ils travaillent avec un petit nombre d'échantillons, les chercheurs doivent cesser d'utiliser une stratégie de code-barres unique, mais devraient plutôt utiliser des échantillons de code-barres doubles. Après le séquençage, les lectures doivent ensuite être filtrées pour n'autoriser que celles dans lesquelles les deux codes-barres sont identiques et ont une séquence attendue. Étant donné que l'échange de codes à barres doubles devrait être un événement rare, cette étape devrait réduire considérablement les erreurs de classification des échantillons.

Lors du multiplexage de nombreux échantillons (nécessitant déjà un double code-barres, quel que soit le problème de saut de code-barres), les chercheurs de Stanford ont conseillé d'utiliser des paires de codes-barres dans lesquelles chaque code-barres individuel n'était utilisé qu'une seule fois. Bien que cette stratégie réduise le nombre de combinaisons de codes-barres disponibles, elle permet tout de même de traiter de nombreux échantillons sur la même piste. En outre, ils ont proposé des stratégies de purification pour éliminer les amorces libres du pool de la bibliothèque.

Illumina a récemment publié un article confirmant l'apparition de sauts d'index. Ils ont déclaré qu'avec le nettoyage approprié des amorces et la mise en œuvre d'autres techniques expérimentales, il devrait être possible de minimiser l'échange de codes-barres à des niveaux négligeables pour la plupart des applications. Le rapport conseillait également de ne multiplexer des conditions similaires que sur une seule voie. Un exemple d'ARN-Seq cerveau vs foie a été donné dans lequel un transcrit est présent dans le foie à un niveau élevé, mais pas du tout dans le cerveau. En raison du saut d'index, ce transcrit peut en fait sembler être exprimé à un faible niveau dans le cerveau. Illumina suggère qu'en combinant uniquement des échantillons similaires (c'est-à-dire ou foie) dans une voie, ce problème sera évité. C'est loin d'être une solution idéale cependant, car les chercheurs devront savoir à l'avance à quoi s'attendre en ce qui concerne les profils d'expression de leurs échantillons et de plus cette stratégie aggraverait le problème des effets de lots.


Méthodes

Pneumocystis les organismes des rongeurs et des humains

P. murina-des échantillons de poumons infectés ont été obtenus à partir de souris femelles knock-out pour le ligand CD40 26 . P. carinii-des échantillons de poumons infectés ont été obtenus de rats mâles Sprague-Dawley immunodéprimés 24 . P. jirovecii-les échantillons pulmonaires d'autopsie infectés provenaient d'un patient atteint du SIDA, qui n'a été infecté qu'avec un seul P. jirovecii souche basée sur le génotypage à cinq loci génomiques différents 53 . Les directives d'expérimentation animale et humaine des National Institutes of Health ont été suivies dans la conduite de ces études.

Préparation d'échantillons d'ADN génomique

Pneumocystis-les tissus pulmonaires infectés ont été découpés en petits morceaux, homogénéisés dans Qiagen Tissuelyser (5 fois pendant 20 s à une fréquence de 1/30), puis centrifugés à 15 000g pendant 6 minutes. Le culot a été remis en suspension dans une solution Trypsine/EDTA (Lonza), incubé à 37°C pendant 30 min et centrifugé à 15 000g pendant 6 minutes. Le culot a été lavé une fois dans du PBS, remis en suspension dans un tampon de lyse contenant 2,9 % (p/v) de collagénase de type I (Gibco) et 0,1 % (p/v) de DNase I (Sigma), incubé à 37 °C pendant 30 min, et centrifugé à 15 000g pendant 6 minutes. Cette digestion enzymatique séquentielle devait éliminer la plupart des cellules hôtes et de l'ADN ainsi que potentiellement une partie des Pneumocystis formes trophiques. Le culot a été lavé trois fois dans du PBS, remis en suspension dans 100 pi de solution de Zymolyase à 0,5 % (p/v) (Zymo Research) et incubé à 37 °C pendant 1 h. L'ADN génomique a été extrait à l'aide du kit de purification d'ADN de levure MasterPure (Epicentre). L'ARN a été éliminé en utilisant une ARNase sans DNase (Epicentre).

Tous les échantillons d'ADN ont été analysés par PCR quantitative en temps réel (qPCR) à l'aide de sondes de transfert d'énergie par résonance de fluorescence 54 pour mesurer la fraction de Pneumocystis ADN comparé à l'ADN de l'hôte dans chaque échantillon. En qPCR pour les deux P. murina 20 et P. carinii 54, la cible était la dihydrofolate réductase à copie unique (dhfr) gène de Pneumocystis et la cible pour l'hôte était une région hautement conservée de la maladie polykystique des reins à copie unique 1 (pkd1) gène 22 . Les cibles pour P. jirovecii et la qPCR humaine étaient les message la famille de gènes 55 et le gène de la -globine à copie unique 56, respectivement. La pureté de chaque échantillon d'ADN a été estimée par le rapport de Pneumocystis pour héberger le nombre de copies du génome et sur la base d'une taille de génome estimée à 8 Mo pour Pneumocystis et 2,7 à 3,3 Go pour l'hôte. Les résultats de la qPCR ont été validés par séquençage 454 ou Illumina MiSeq d'échantillons sélectionnés. échantillons d'ADN extraits de P. murina ou P. carinii préparations contenant jusqu'à 90 % Pneumocystis ADN, et de enrichi P. jirovecii préparations, 10-25% P. jirovecii ADN, tandis que ceux extraits directement de tissus pulmonaires fortement infectés contenaient <0,4 à 1 % Pneumocystis ADN.

P. murina séquençage et assemblage du génome

Pour évaluer la qualité de la P. murina Les préparations d'ADN pour le séquençage du génome entier, le séquençage à petite échelle ont d'abord été effectués à l'aide d'un séquenceur en titane 454 GS FLX (Roche Applied Science) à Leidos, Inc. (Frederick, MD) selon les 454 protocoles de séquençage standard. Un total de 4 millions de lectures (avec une longueur moyenne de 344 bases) ont été générés, avec ∼ 40 % d'entre eux ayant des coups de souffle avec P. carinii 7 . Par la suite, un séquençage profond a été réalisé à l'aide de la technologie Illumina HiSeq au Broad Institute (Cambridge, MA). Sept différents P. murina Préparations d'ADN (chacune avec une pureté de 60 à 90 % pour P. murina par qPCR) ont été utilisées pour construire de petites bibliothèques d'inserts chacune a été séquencée, et la bibliothèque (préparation B123 à partir d'une seule souris, projet central G11228) avec le pourcentage le plus faible d'ADN de souris hôte contaminant a été utilisée pour l'assemblage. À partir de cette bibliothèque, un total de 34 millions de lectures appariées de 101 bases avec une taille d'insertion moyenne de 153 bases ont été générées. Une deuxième bibliothèque d'inserts plus grande (préparation C2 à partir d'une autre souris, projet central G11230) avec une taille d'insert moyenne de 1 247 bases a été préparée, et un total de 83 millions de lectures appariées à 101 bases ont été générées. Après suppression des séquences de souris (version mm9) et P. murina séquences mitochondriales 22 , ces lectures ont été assemblées avec Allpaths (version R37380) avec les paramètres par défaut 57 l'assemblage résultant a été criblé pour les séquences contaminantes (souris et bactéries) et pour éliminer P. murina séquences mitochondriales 22 en s'alignant sur la base de données de nucléotides non redondante de NCBI et en supprimant les échafaudages avec une couverture élevée qui correspondaient à des organismes non fongiques. Le projet d'assemblage de 24 échafaudages a été encore amélioré en utilisant 454 lectures longues et la PCR pour prendre en charge les jointures d'échafaudage. L'assemblage final contenait 17 échafaudages. Tous les écarts de contig internes ont été comblés par PCR et séquençage de Sanger. Les extrémités des échafaudages sans message les gènes ou les répétitions des télomères ont été étendus par PCR et/ou séquençage PacBio (méthodes supplémentaires).

P. carinii séquençage et assemblage du génome

Après que le séquençage préliminaire de 454 ait démontré une pureté de 63 %, 1 P. carinii Préparation d'ADN (n° B80 de 2 rats fortement infectés, pureté 80% pour P. carinii par qPCR) a été utilisé pour construire 2 bibliothèques avec une taille moyenne d'insert de 180 bases et 793 bases, respectivement. Chacun a été séquencé à l'aide de lectures appariées de 101 bases sur la plate-forme Illumina MiSeq du Broad Institute et, après élimination des séquences de rat hôte et P. carinii séquences mitochondriales 22 , assemblées en utilisant AllPaths-LG (version R47825) avec des paramètres par défaut, résultant en un total de 53 échafaudages. Après avoir rejoint les échafaudages à l'aide de 454 lectures longues et confirmé par PCR, l'assemblage final du génome contenait 17 échafaudages. Tous les écarts de contig internes ont été comblés par PCR et séquençage de Sanger. Les extrémités de six échafaudages se chevauchaient avec les sept séquences de télomères signalées ailleurs 58, elles ont été fusionnées et confirmées par PCR et/ou en alignant les échafaudages fusionnés avec Illumina et 454 lectures brutes. Les extrémités de plusieurs échafaudages sans message les gènes, les gènes de la kexine ou les répétitions des télomères ont été étendus par séquençage PCR et/ou PacBio (méthodes supplémentaires).

P. jirovecii séquençage et assemblage du génome

Trois échantillons d'ADN enrichis (RU7, RU12 et RU817) provenant d'un seul échantillon pulmonaire d'autopsie RU (1,4 à 2 g chacun, avec 10 à 25 % P. jirovecii ADN par qPCR) ont été utilisées pour construire trois banques (taille d'insert moyenne de 156 bases) chacune a été séquencée sur la plate-forme Illumina HiSeq du Broad Institute. Pour améliorer les assemblages initiaux avec un nombre élevé de contigs, une séquence supplémentaire a été générée à partir d'échantillons d'ADN sélectionnés par hybride, comme décrit précédemment. les sondes ont été conçues pour : les assemblages existants 8 comprenant une densité de sondes plus élevée aux extrémités contig, les lectures non assemblées 8 qui ne correspondaient pas au génome humain, les transcrits RNA-Seq assemblés 8 qui ne correspondaient pas aux assemblages existants ou au génome humain et non ancrés message séquences 25 . Les séquences hôtes ont été supprimées par alignement sur l'assemblage du génome humain 19 (GRCh37/hg19), et P. jirovecii les séquences mitochondriales 22 ont été supprimées. Les 101 lectures appariées de bases restantes ont été assemblées séparément avec Spades 2.5.1 (au Broad Institute), Spades 3.0 (à Leidos) et Abyss (Biowulf au NIH), ce qui a donné respectivement 400, 149 et 312 contigs. Tous ces contigs ont été fusionnés avec Sequencher (Gene Codes Co., Ann Arbour, MI), résultant en un total de 20 échafaudages. Toutes les lacunes internes ont été comblées par PCR et séquençage de Sanger. Les échafaudages ont été validés en examinant les lectures brutes alignées par Seqman Pro (DNAStar, Madison, WI) pour s'assurer qu'il n'y avait pas de fausses jointures. Pour améliorer la couverture de message gènes dans l'assemblage, le séquençage PacBio a été effectué comme décrit dans les méthodes supplémentaires.

Électrophorèse chromosomique et hybridation Southern

Pour comparer le P. murina assemblage avec le nombre et la longueur des chromosomes, nous avons utilisé une électrophorèse à champ électrique homogène (CHEF) pour séparer les chromosomes. P. murina les organismes ont été partiellement purifiés à partir de poumons frais de souris par centrifugation en gradient de densité Ficoll-Hypaque 60, puis traités pour CHEF 61 en utilisant le kit CHEF Yeast Genomic DNA Plug (Bio-Rad). En bref, partiellement purifié P. murina les organismes ont été incorporés dans de l'agarose CleanCut à 0,8 %, traités avec 24 U ml -1 de protéinase K pendant une nuit à 50 °C, puis lavés quatre fois dans 1 × Wash Buffer et conservés à 4 °C. L'électrophorèse a été réalisée dans des gels d'agarose à 1% (14 × 21 cm) en utilisant l'appareil CHEF DR II (Bio-Rad). Les gels ont été traités dans du tampon 0,5 × TBE pendant 144 h à 135 V et 12,5 °C, avec une impulsion initiale de 50 s qui a été progressivement augmentée jusqu'à 100 s. Le gel a été coloré avec du bromure d'éthidium (Sigma), puis l'ADN a été transféré sur des membranes Nytran (Schleicher & Schuell) dans des conditions neutres 61 . Nous avons utilisé deux transferts préparés à partir du même tampon d'ADN, et chaque transfert a été hybridé consécutivement à différentes sondes, avec stripping du transfert entre les hybridations. Toutes les sondes étaient des fragments d'ADN double brin amplifiés par PCR marqués à l'aide du kit PCR DIG Probe Synthesis (Roche Applied Science), à ​​l'exception de la sonde pour les répétitions des télomères, qui était un oligonucléotide synthétisé marqué à l'aide du kit DIG Oligonucleotide Tailing (Roche Applied Science). Toutes les séquences d'amorces et de sondes sont fournies dans les données supplémentaires 27. L'hybridation et la détection du signal ont été effectuées en utilisant le système d'hybridation de sonde DIG (Roche Applied Science) comme décrit précédemment 22 .

RNA-Seq, niveaux d'expression et prédiction de gènes

Trois échantillons indépendants de P. murina et P. carinii organismes ont été partiellement purifiés à partir de poumons fortement infectés de trois souris et rats, respectivement, par centrifugation en gradient de densité Ficoll-Hypaque 60 . Les pastilles partiellement purifiées devaient contenir à la fois des formes kystiques et trophiques. L'ARN total a été isolé à l'aide du RNeasy Mini Kit (Qiagen). Le rapport de Pneumocystis à l'ARN hôte a été estimée entre 1:4 et 8:1 sur la base de la densité des bandes d'ARN dans les gels d'agarose pour l'ARNr 28S de l'hôte et Pneumocystis ARNr 26S. Des bibliothèques spécifiques au brin ont été construites en utilisant la méthode de marquage du deuxième brin dUTP 62,63, sauf comme indiqué ci-dessous. Pour tous les échantillons, l'ARN poly(A) a été purifié à l'aide du kit de purification d'ARNm Dynabeads (Life Technologies), avec deux cycles de sélection (avec régénération des billes), résultant en <5% d'ARNr tel que mesuré par la puce Bioanalyzer RNA 6000 Pico (Agilent) programme. L'ARNm résultant (200 ng par échantillon) a été traité avec le kit Turbo DNA-free (Ambion), vérifié par qPCR pour l'absence d'ADN génomique détectable (tableau supplémentaire 5), et fragmenté dans 1 × tampon de fragmentation d'ARN (Affymetrix) pour 4 min à 80°C. Après la synthèse de l'ADN complémentaire du premier brin (ADNc), un volume de 2,0 × volume de billes RNAClean SPRI (Beckman Coulter Genomics) a été utilisé à la place de l'extraction au phénol/chloroforme/alcool isoamylique (25:24:1) et de la précipitation à l'éthanol. La sélection de la taille après la ligature de l'adaptateur a été effectuée avec deux étapes de nettoyage de billes 0,7 × Ampure (Beckman Coulter Genomics) et huit cycles de PCR ont été utilisés pour générer la bibliothèque. Les bibliothèques ont été séquencées sur la plate-forme HiSeq2000 générant en moyenne 26,1 millions de lectures de 76 bases appariées pour chacun des 6 échantillons.

Pour examiner les niveaux d'expression, les lectures d'ARN-Seq de chacun des trois échantillons pour chaque organisme ont été alignées sur les séquences d'ADN codantes (CDS) extraites à l'aide d'un nœud papillon 65 . Les fichiers bam d'alignement ont ensuite été utilisés pour quantifier les abondances de transcrits par RSEM 66 . Nous avons sélectionné le top 5% des niveaux d'expression pour examiner l'enrichissement fonctionnel par GSEA 45 . Pour chaque organisme, nous avons examiné la GSEA des ensembles de gènes classés sur la base des annotations fonctionnelles Pfam, TIGRFAM, KEGG et SignalP. Nous avons également exécuté GSEA sur certains ensembles de gènes sélectionnés manuellement tels que message gènes.

Le niveau d'expression relatif pour chaque gène a été exprimé en fragments par kb d'exon par million de fragments cartographiés (FPKM). Pour les deux P. murina et P. carinii, 99% de tous les gènes annotés ont été exprimés (FPKM>2), avec une couverture médiane de 151 FPKM chez les deux espèces. De plus, les gènes exprimés ont montré une bonne profondeur d'alignement à travers le transcrit, avec 98% de tous les gènes de chaque espèce ayant une profondeur d'ARN-Seq au moins cinq fois supérieure. Tous les gènes sans expression détectable (26 dans P. murina et 23 en P. carinii) contiennent des CDS courts (taille moyenne 250 pb, potentiellement des pseudogènes) ou incomplets.

Les P. murina et P. carinii les génomes ont été annotés à l'aide d'une combinaison de données d'expression (RNA-Seq), d'informations d'homologie et ab initio méthodes de recherche de gènes (méthodes supplémentaires) telles que décrites précédemment 67 . Ces méthodes ont également été utilisées pour annoter l'assemblage RU de P. jirovecii, à l'exception du fait que les données RNA-Seq n'ont pas été utilisées (méthodes supplémentaires). Les ensembles de gènes ont été comparés sur la base d'annotations fonctionnelles et utilisés comme base pour l'analyse synténique et phylogénétique (méthodes supplémentaires).

Analyse de la variation de déformation dans P. jirovecii

Pour examiner les polymorphismes de nucléotide unique (SNP) à l'échelle du génome entre les isolats séquencés, l'assemblage RU7 généré dans cette étude a été utilisé comme référence par rapport aux lectures brutes SE8 8 . Les lectures SE8 ont été récupérées de NCBI (ERR135854 à ERR135863) et converties au format fastq. Les fichiers fastq individuels ont été concaténés et l'ensemble de lecture complet aligné sur l'assemblage RU7 à l'aide de BWA-MEM 68 et converti en bam trié à l'aide de SAMtools 64. Il en a résulté une profondeur d'alignement très élevée sur toutes les bases d'assemblage, avec une profondeur d'alignement médiane de 1 046.

La boîte à outils d'analyse du génome 69 (GATK) v2.7-4-g6f46d11 a été utilisée pour appeler les bases de variantes et de référence à partir des alignements. En bref, les outils Picard (http://picard.sourceforge.net) AddOrReplaceReadGroups, MarkDuplicates, CreateSequenceDictionary et ReorderSam ont été utilisés pour prétraiter les alignements, suivis de GATK RealignerTargetCreator et IndelRealigner pour résoudre les lectures désalignées proches des indels. Ensuite, GATK UnifiedGenotyper (avec le modèle de vraisemblance du génotype haploïde (GLM)) a été exécuté avec SNP et INDEL GLM. Nous avons également exécuté BaseRecalibrator et PrintReads pour le recalibrage du score de qualité de base sur les sites appelés à l'aide de GLM SNP et des variantes réappelées avec UnifiedGenotyper émettant tous les sites. Une dernière étape de filtrage a été utilisée pour supprimer toute position appelée par les deux GLM (c'est-à-dire les indels et les SNP incompatibles) et pour supprimer les positions avec un support de faible profondeur de lecture (<10). Les 24 902 SNP ont ensuite été mappés sur l'ensemble de gènes prédit sur l'assemblage SE8, un total de 14 135 SNP sont situés dans les régions CDS, avec des substitutions affectant les régions codant pour les protéines comme suit : 7 410 non synonymes, 6 690 synonymes, 26 non-sens et 9 extensions.La densité de SNP à travers l'assemblage SE8 a été calculée pour des fenêtres de 5 kb.

Détection de la chitine par marquage par fluorescence

La séquence d'ADNc codant pour le CBD de Bacillus circulans Le gène 70 de la chitinase A1 (2 183 à 2 338 pb, code d'accession GenBank M57601.1) a été optimisé pour l'expression bactérienne (GenScript USA Inc.) et modifié par l'ajout de balises ATG et 3 HA. La séquence modifiée a été synthétisée (GenScript USA Inc., Fig. 15 supplémentaire) et clonée dans le vecteur pET-28 (EMD Biosciences). La construction d'ADNc a été transformée en Escherichia coli souche BL21 (DE3) RIL (Agilent technologies) et exprimée en tant que protéine de fusion His tag. La protéine exprimée a été purifiée en deux étapes, d'abord à l'aide d'un anticorps anti-CBD (New England Biolabs) qui a été immobilisé à l'aide du kit d'immobilisation AminoLink Plus (Thermo Scientific), puis du kit de purification au cobalt Hispur (Thermo Scientific). La protéine purifiée a été biotinylée à l'aide du kit de conjugaison Lightning-Link Biotin (Type A) (Innova Biosciences Ltd) et utilisée pour le marquage par fluorescence (Histoserv, Inc., Germantown, MD) de P. murina-coupes de tissus pulmonaires infectés fixées dans Histochoice (Amresco, Inc.). La streptavidine conjuguée à l'Alexafluor 488 a été utilisée pour la détection du CBD. Pneumocystis les organismes ont été détectés par un anticorps anti-Msg 26 et les kystes par une protéine recombinante dectine-Fc (aimablement fournie par le Dr Chad Steele, Université d'Alabama à Birmingham, Alabama) 71 . Sections de rein de souris infectées par Candidose (gracieusement fourni par le Dr Michail Lionakis, National Institute of Allergy and Infectious Diseases, Bethesda, Maryland) 72 ont été utilisés comme contrôle positif pour la coloration de la chitine (Fig. 5).

Analyse des liaisons chitine glycosyle

Pour préparer la fraction de paroi cellulaire, P. carinii les organismes ont été partiellement purifiés à partir de poumons de rats infectés par centrifugation en gradient de densité Ficoll-Hypaque 60 S. cerevisiae les organismes ont été obtenus à partir de Stratagene (souche YPH499) et cultivés dans une culture de levure standard. Les cellules ont été remises en suspension dans un tampon dénaturant (2% SDS, 1 mM EDTA dans 50 mM Tris pH 8,0) et chauffées à 100°C pendant 10 min. Le culot de paroi cellulaire a été récupéré par centrifugation à 5000g pendant 5 minutes, lavé deux fois avec du PBS et séché dans un concentrateur sous vide rapide.

Pour l'analyse de la chitine, les culots de paroi cellulaire ont été incubés avec 0,1 % (p/v) de chitinase (Sigma) dans un tampon d'acétate de sodium 50 mM, pH 5,6, à 37 °C pendant 72 h. Les matières insolubles après les digestions initiales ont été davantage traitées séquentiellement avec de la pronase (Roche), de la lyticase (Sigma) et de la chitinase pour rechercher la chitine résiduelle. La digestion des protéines a été réalisée en utilisant 0,01 % (p/v) de pronase dans un tampon Tris-HCl pH 8,2 avec 10 mM de CaCl2 à 37°C pendant 48h. La digestion du glucane a été réalisée avec 0,04 % (p/v) de lyticase dans du tampon phosphate de sodium 50 mM, pH 7,5 à 37 °C pendant 24 h. La digestion par la chitinase a été effectuée comme ci-dessus. Après chaque digestion, l'enzyme a été inactivée par chauffage à 100°C pendant 5 min.

Les composants de la paroi cellulaire digérés ont été analysés par MS temps de vol à ionisation laser-désorption assistée par matrice (MALDI/TOF-MS) et MS à ionisation nanospray (NSI-MSn). Une aliquote de surnageant après chaque digestion enzymatique a été analysée par MS pour déterminer la teneur en oligosaccharides. Avant l'analyse MS, les enzymes ont été éliminées par passage à travers une cartouche C18 Sep-PAK. L'écoulement a été récupéré dans de l'acide acétique à 5 % et lyophilisé. Les échantillons séchés ont été perméthylés en utilisant des méthodes déjà publiées 73 et profilés par MS. MALDI/TOF-MS a été réalisé en mode réflecteur d'ions positifs en utilisant du DHBA (acide α-dihyroxybenzoïque, solution à 20 mg ml -1 dans 50 % de méthanol/eau) comme matrice. Le spectre a été obtenu en utilisant un AB SCIEX TOF/TOF 5800 (AB SCIEX). L'analyse NSI-MSn a été réalisée sur un spectromètre de masse LTQ Orbitrap XL (Thermo Fisher) équipé d'une source d'ions nanospray. L'échantillon perméthylé a été dissous dans 1 mM de NaOH dans 50 % de méthanol et infusé directement dans l'instrument à un débit constant de 0,5 l min -1 . Un spectre MS complet de transformée de Fourier a été collecté à une résolution de 30 000. La température capillaire a été fixée à 210 °C et l'analyse MS a été réalisée en mode ion positif. Pour la cartographie ionique totale (analyse MS/MS automatisée), un m/z une plage de 200 à 2 000 a été balayée avec le mode MS de piège à ions dans des fenêtres successives de 2,8 unités de masse qui chevauchaient la fenêtre précédente de 2 UM.

Pour la détermination des liaisons glycosylées, des acétates d'alditol partiellement méthylés (PMAA) ont été préparés à partir de glycanes entièrement perméthylés. Les glycanes perméthylés ont été hydrolysés avec HCl/eau/acide acétique (0,5:1,5:8, en vol.) à 80 °C pendant 18 h, suivi d'une réduction avec 1% de NaBD4 dans 30 mM NaOH pendant une nuit, puis acétylation avec anhydride acétique/pyridine (1:1, v/v) à 100 °C pendant 15 min et analyse par chromatographie en phase gazeuse-MS (GC-MS) sur un GC Agilent 7890A interfacé à un 5975C MSD (détecteur sélectif de masse, mode ionisation par impact électronique). La séparation du PMAA a été réalisée sur une colonne capillaire de 30 m de silice fondue SP2331 (Supelco) pour les dérivés de sucres neutres, et une colonne DB-1 (Agilent) pour les dérivés de sucres aminés.

Identification de la glycosylation des protéines Msg

Purifier P. carinii Protéines Msg, partiellement purifiées P. carinii les organismes 60 ont été remis en suspension dans du SDS à 2 % dans du tampon Tris-HCl 50 mM, pH 8,0, bouillis pendant 10 min et centrifugés. Le surnageant a été récupéré. La procédure d'extraction a été répétée deux fois de plus, les surnageants ont été regroupés et le SDS a été éliminé en utilisant le kit de précipitation SDS-Out SDS (Thermo Scientific). Le Msg solubilisé a été purifié par affinité à l'aide d'un anticorps monoclonal anti-Msg (RA-E7, don des Drs Peter Walzer et Michael Linke 74, Université de Cincinnati, Cincinnati, Ohio) immobilisé sur une colonne AminoLink plus (kit d'immobilisation AminoLink plus) suivant les instructions du fabricant. instructions. L'extrait protéique a été dilué avec un volume égal de solution saline tamponnée Tris et incubé avec l'anticorps immobilisé dans la colonne par mélange bout à bout pendant une nuit à 4°C. La colonne a ensuite été essorée et lavée avec du tampon phosphate de sodium 0,1 M pH 6,9. Le Msg a été élué avec 150 mM d'hydroxyde d'ammonium et les fractions ont été neutralisées avec 1 M de NaH2Bon de commande4. Les fractions contenant Msg ont été regroupées et concentrées en utilisant une unité de filtration centrifuge Microcon-30 kDa (Millipore) et le tampon a été remplacé par du tampon phosphate de sodium 0,1 M pH 6,9. Les aliquotes ont été conservées à -80°C.

Les protéines Msg purifiées ont été soumises à une détection de glycopeptides par LC-MS. L'analyse glycoprotéomique a été réalisée comme décrit précédemment 75 . En bref, les peptides digérés par tryptique des protéines Msg purifiées ont été séparés avec une colonne Magic C18 (15 cm × 75 m, Bruker-Michrom, CA) et analysés avec un spectromètre de masse Orbitrap Fusion (Thermo Scientific) après que Msg ait été dénaturé, réduit, alkylé et dessalé. L'élution par gradient a été réalisée avec un gradient linéaire de 30 min de 5 à 35 % d'acétonitrile dans 0,1 % d'acide formique à un débit de 300 nl min -1 . Les données MS ont été traitées automatiquement par un pipeline de balayage MS dépendant des données, qui intégrait trois méthodes de dissociation, notamment la dissociation collisionnelle à haute énergie (HCD), la dissociation par transfert d'électrons (ETD) et la dissociation induite par collision (CID). Dans ce pipeline, le spectre complet MS a d'abord été acquis à partir des ions les plus abondants entre m/z 350 et m/z 1 550 avec un temps de cycle de 3 secondes à 120 000 résolutions en mode FT. Par la suite, les données MS acquises ont été évaluées par MS/MS avec le mode ETD dépendant du produit HCD ou le mode ETD/CID dépendant du produit HCD avec les ions les plus abondants dans un temps de cycle de 3 s. Pour le déclencheur d'ions HCD MS/MS, les ions glycane oxonium à m/z 204.0867 (HexNAc), m/z 138.0545 (fragment HexNAc) ou m/z 366.1396 (Hexose-HexNAc) dans les spectres HCD ont été utilisés pour déclencher l'acquisition ETD ou CID. Les MS/MS ont été mesurés à une résolution de 15 000.

Pour la cartographie des glycopeptides de P. carinii Les résumés Msg, les données LC-MS ont été analysés à l'aide du logiciel Byonic (Protein Metrics) et l'annotation des données par le logiciel a été validée manuellement. Les paramètres byonic ont été réglés pour autoriser 30 h 00. de tolérance de masse des ions précurseurs et 3,0 p.p.m. de tolérance aux ions fragments avec une masse monoisotopique. Les peptides digérés ont été autorisés avec jusqu'à trois sites de clivage internes manqués, et les modifications différentielles de 57,02146 et 15,9949 Da ont été autorisées pour la cystéine alkylée et l'oxydation de la méthionine, respectivement. La base de données de protéines utilisée pour le blast protéique comprenait l'ensemble de protéines du génome du rat au NCBI (version 104) et du P. carinii génome B80 généré dans cette étude. La base de données de glycanes utilisée était la base de données de N-glycanes de mammifères 309 (base de données par défaut dans le logiciel Byonic). L'analyse des glycanes N-liés libérés de Msg a été réalisée avant la cartographie des glycopeptides.

L'annotation des données par le logiciel a été en outre filtrée comme suit. Toute identification de protéine avec FDR >1%, probabilité de log <4 ou meilleur score Byonic <500 a été exclue. De plus, toutes les annotations de spectres glycopeptidiques avec un score de Byonic <400 ont été exclues. Les identifications de glycopeptides qui restaient après ces filtres stricts ont été validées manuellement en examinant les ions fragments correspondant à la fragmentation théorique des glycopeptides, la présence d'une série d'ions glycane oxonium attendus et la perte neutre de la fraction glycane pour les données MS/MS-HCD.


Méthodes

Identification des UCE

Nous avons identifié les UCE en criblant les alignements du génome entier du poulet (Gallus gallus) et Carolina anole (Anolis carolinensis) préparé par le groupe de bioinformatique du génome de l'UCSC à l'aide d'un script Python personnalisé (http://www.python.org/) pour identifier des séquences d'au moins 60 bases ayant 100 % d'identité de séquence. Nous avons ensuite aligné chaque région conservée des alignements poulet-lézard au génome du diamant mandarin (UCSC taeGut1) en utilisant un programme Python personnalisé et BLAST ( Altschul et al. 1997), et nous avons stocké des métadonnées pour les correspondances ayant une valeur e ≤ 1 × 1 0 − 15 dans une base de données relationnelle (RDB) avec les premiers résultats de dépistage. Nous avons supprimé les doublons du groupe de correspondances contenant des données de poulet, de lézard et de diamant mandarin, et nous avons défini l'ensemble restant de 5599 séquences uniques comme UCE. Nous avons estimé la distance moyenne (± 95% CI) entre chacun de ces UCE en utilisant des positions dans le génome du poulet (UCSC galgal3) car le génome du poulet est actuellement le génome aviaire ou reptile le plus complet et le mieux assemblé.

Conception de sondes à partir d'UCE

Nous avons conçu des sondes d'enrichissement cible en sélectionnant des UCE dans la RDB, en ajoutant une séquence à ces UCE de longueur inférieure à 120 pb en sélectionnant des quantités égales de séquences flanquantes 5' et 3' à partir d'un assemblage de génome de poulet à répétition masquée, et en enregistrant la longueur de flanquement séquence, le cas échéant, ajoutée à chacun. Nous avons masqué tous les UCE tamponnés contenant des régions de type répétition à l'aide de RepeatMasker (http://www.repeatmasker.org/) avant la conception de la sonde. Si les UCE étaient >180 pb, nous avons carrelé les sondes de 120 pb dans les régions cibles à une densité de 2 × (c'est-à-dire que les sondes se chevauchaient de 60 pb). Si les UCE avaient une longueur totale de <180 pb, nous avons sélectionné une seule sonde au centre de l'UCE. Nous avons utilisé LASTZ (disponible sur http://www.bx.psu.edu/miller_lab/) pour aligner les sondes sur elles-mêmes et pour identifier et supprimer les doublons résultant de la conception des sondes. Nous avons inséré ces 5561 sondes dans le RBD et nous avons mis à jour chaque enregistrement de sonde avec des données supplémentaires indiquant si les sondes contenaient des bases ambiguës (N), le contenu Tm et GC de la sonde, le nombre de bases ajoutées pour mettre en mémoire tampon un UCE particulier à la longueur de la sonde. (120 pb), le nombre de bases masquées dans les sondes conçues et les types de mésappariements que nous avons observés pour l'UCE parent de chaque sonde lors du BLASTing d'UCE de poulet-anole contre le diamant mandarin.

Alignement de sondes conçues sur dix génomes amniotes

Nous avons aligné 5561 sondes sur dix génomes amniotes à l'aide d'un programme wrapper Python autour de LASTZ pour faciliter le traitement parallèle des données. Nous n'avons retenu que les correspondances ayant ≥ 92,5% d'identité sur ≥ 100 pb de la séquence de sonde de 120 pb (83%). Nous avons utilisé un programme Python personnalisé pour filtrer les correspondances LASTZ pour les doublons réciproques et non réciproques, et nous avons également exclu les correspondances où le nombre de correspondances observé était inférieur au nombre de sondes conçues. Par exemple, si nous avons superposé deux sondes sur un locus UCE, mais que LASTZ n'a fait correspondre qu'une seule sonde à la séquence du génome, nous avons laissé tomber le locus UCE parent de tout examen plus approfondi.

Enrichissement cible in silico de loci UCE à partir de primates

Pour tester notre flux de travail in vitro putatif, nous avons aligné 5561 sondes sur neuf génomes de primates et un groupe externe de souris à l'aide de LASTZ. Comme ci-dessus, nous n'avons retenu que les correspondances ayant ≥92,5% d'identité sur ≥100 pb de la séquence de sonde de 120 pb (83%), et nous avons ignoré les doublons réciproques et non réciproques pour filtrer les paralogues potentiels. Nous avons également exclu les loci UCE où le nombre de sondes correspondant au génome était inférieur à celui attendu. À travers cet ensemble réduit de correspondances et chaque taxon de primate, nous avons tranché l'emplacement d'alignement des sondes restantes des génomes de référence plus 200 pb de séquence flanquante en amont (5′) et en aval (3′) de l'emplacement d'alignement pour donner une tranche totale de environ 520 pb. Nous avons choisi cette longueur de séquence flanquante car les calculs préliminaires ont révélé que cette longueur était probablement proche de la taille de contig maximale à laquelle nous pouvions nous attendre en utilisant les techniques de préparation de la bibliothèque Illumina Nextera. Nous avons assemblé les sondes + flanc arrière dans les UCE parents qu'ils représentaient à l'aide d'un programme Python personnalisé qui intégrait LASTZ - pour faire correspondre les sondes à leur UCE - et MUSCLE (Edgar 2004) pour assembler plusieurs sondes conçues à partir du même UCE parent. Après assemblage, nous nous référons à chaque UCE plus séquence flanquante comme un locus. Pour chaque locus, nous avons aligné les données sur toutes les espèces de primates à l'aide d'un programme Python personnalisé et de MUSCLE. Nous avons utilisé une moyenne mobile sur une fenêtre de 20 pb pour couper les extrémités de tous les alignements, garantir que les extrémités contenaient au moins 50 % d'identité de séquence et supprimer la séquence mal alignée. Nous avons autorisé les insertions aux extrémités de l'alignement tant que les insertions étaient présentes dans moins de 30 % des taxons individuels au sein d'un alignement donné. Nous avons exclu les loci qui n'étaient pas trouvés dans toutes les espèces de primates et nous avons abandonné les alignements qui manquaient d'espèces de primates. Cela a abouti à une matrice de données complète sans aucun loci manquant.

Enrichissement ciblé in vitro de loci UCE à partir d'oiseaux

A partir de l'ensemble de 5561 sondes, nous avons sélectionné un sous-ensemble de 2560 sondes pour la synthèse où les sondes avaient <10 bases masquées et <50 bases ajoutées (25 de chaque côté). Ces sondes représentaient 2386 UCE conservées chez le poulet, le lézard et le diamant mandarin. Nous avons fait synthétiser ces sondes dans le commerce dans un kit MySelect personnalisé (Mycroarray, Inc.). Nous avons effectué des extractions d'ADN au phénol-chloroforme sur des tissus d'oiseaux à partir d'échantillons de musée justificatifs (tableau supplémentaire 1, disponible auprès de doi: 10.5061/dryad.64dv0tg1), et nous avons préparé des bibliothèques pour le séquençage d'Illumina à l'aide de kits de préparation de bibliothèques Illumina Nextera (Epicentre Biotechnologies).

Pour enrichir les UCE ciblées, nous avons suivi le workflow de base pour l'enrichissement de cibles basé sur des solutions ( Gnirke et al. 2009) avec plusieurs modifications pour les bibliothèques préparées par Nextera. Tout d'abord, nous avons utilisé 1,8X AMPure XP à la place du nettoyage de la colonne après la réaction de tagmentation Nextera, car le nettoyage de la colonne Zymo recommandé a donné des concentrations finales d'ADN inférieures à celles de l'AMPure. Deuxièmement, nous avons amplifié des bibliothèques étiquetées à l'aide d'amorces PCR de longueur réduite et purifiées par HPLC, complémentaires des adaptateurs attachés à chaque fragment d'ADN lors de l'étiquetage. Nous n'avons pas attaché d'étiquettes de séquence aux bibliothèques à ce stade pour réduire les complications potentielles lors de l'enrichissement introduites par des adaptateurs plus longs et variables individuellement. Nous avons augmenté le nombre de cycles de PCR pendant la PCR post-tagmentation à 16 ou 19 pour produire un modèle suffisant (∼ 500 ng) pour l'enrichissement. Après la préparation de la bibliothèque, nous avons substitué un mélange bloquant de 500 M (chacun) d'oligos composés des compléments directs et inverses des adaptateurs Nextera pour le mélange bloquant d'adaptateur fourni par le kit (#3), et nous avons individuellement enrichi des bibliothèques spécifiques à l'espèce en utilisant le sondes d'ARN synthétique décrites ci-dessus. Nous avons incubé les réactions d'hybridation à 65°C pendant 24 h sur un thermocycleur. Après l'hybridation, nous avons enrichi les échantillons en liant l'ADN hybride à des billes recouvertes de streptavidine et nous avons élué l'ADN des billes recouvertes de streptavidine en utilisant 50 L de NaOH, que nous avons neutralisé avec 50 L supplémentaires de Tris-HCl. Nous avons nettoyé l'ADN élué en nous liant à des billes AMPure XP 1,8X (v/v) et nous avons remis en suspension l'ADN enrichi propre dans 30 L d'eau sans nucléase. Nous avons amplifié 14 L d'ADN propre et enrichi dans une réaction PCR de 50 L combinant des amorces directes, inverses et d'indexation avec Nextera Taq ou Phusion DNA Polymerase (New England Biolabs) pour ajouter un ensemble personnalisé de 24 adaptateurs d'indexation Nextera, robustes aux insertions , suppressions et substitutions à chaque échantillon. Après la PCR, nous avons nettoyé les réactions en liant les amplicons à 1,8X (v/v) AMPure XP. Nous avons quantifié les bibliothèques indexées enrichies à l'aide d'un bioanalyseur et nous avons combiné les bibliothèques pour le séquençage à des concentrations équimolaires en supposant que tous les fragments liés par adaptateur étaient de longueur moyenne dans toutes les bibliothèques.

Nous avons séquencé des bibliothèques à l'aide des amorces de séquençage Nextera et d'un analyseur Illumina Genome Analyzer IIx de 100 pb mené par le LSU Genomics Facility. Après le séquençage, la LSU Genomics Facility a démultiplexé les bibliothèques à l'aide du logiciel fourni par Illumina, et nous avons utilisé un pipeline (https://github.com/faircloth-lab/illumiprocessor) implémenté en Python pour éliminer la contamination de l'adaptateur des lectures à l'aide de SCYTHE (https:/ /github.com/vsbuffalo/scythe), les lectures de qualité adaptative à l'aide de SICKLE (https://github.com/najoshi/sickle), excluent les lectures contenant des bases ambiguës (N) et collectent des métadonnées pour chaque groupe de lectures analysé.

Après le pré-traitement des lectures, nous avons assemblé les lectures spécifiques aux espèces dans des contigs à l'aide de VELVET ( Zerbino et Birney 2008) via VelvetOptimiser (http://bioinformatics.net.au/software.velvetoptimiser.shtml), que nous avons utilisé avec les paramètres par défaut ( krange = 59 –79) pour optimiser la longueur, la couverture et la coupure du kmer pour l'assemblage des données de chaque espèce. Velvet résout les sites potentiellement variables à l'allèle majoritaire ( Zerbino et Birney 2008). Nous avons converti les contigs sortis par VELVET au format de banque AMOS, et nous avons calculé la couverture dans les contigs en utilisant un code Python personnalisé pour analyser la sortie des programmes cvgStat et analyser-read-depth fournis dans le progiciel AMOS 3.0.0 (http://sourceforge. net/projets/amos). Nous avons utilisé un programme Python personnalisé intégrant LASTZ (match_contigs_to_probes.py) pour aligner les contigs assemblés sur leur région de sonde/UCE respective tout en supprimant les doublons réciproques et non réciproques et les sondes ayant moins de correspondances que prévu, comme décrit ci-dessus.

Ce programme (match_contigs_to_probes.py) crée une base de données relationnelle des correspondances avec les loci UCE par taxon, et nous avons utilisé un deuxième programme (get_match_counts.py) pour interroger cette base de données et produire un fichier fasta contenant uniquement les contigs construits à partir des UCE présents dans chaque taxon .Ce programme a également la capacité d'inclure des loci UCE tirés de séquences génomiques existantes, dans le but principal d'ajouter des données d'exogroupe de haute qualité à partir de taxons génomiques. Nous avons utilisé cette fonctionnalité pour inclure les loci UCE identifiés dans l'anole de Caroline (UCSC anoCar2) en tant que données d'exogroupe pour la phylogénie des oiseaux. Nous avons aligné et rogné les lectures comme décrit ci-dessus. Nous avons utilisé l'inférence multimodèle et la moyenne des modèles ( Burnham et Anderson 2002) de modèles linéaires généralisés binomiaux (Calcagno et de Mazancount 2010 R Core Development Team 2011) pour évaluer l'effet de combinaisons raisonnables des paramètres suivants sur l'enrichissement des loci UCE (détecté , non détecté) : contenu GC UCE, longueur UCE, sonde TM, nombre de sondes, bases masquées incluses dans les sondes, bases ajoutées aux sondes tampons et taxon.

Estimation des modèles de substitution

Nous avons utilisé des programmes Python personnalisés (run_mraic.py) enveloppant un MrAIC 1.4.4 modifié ( Nylander 2004) pour estimer, en parallèle, les modèles de substitution de sites finis les plus probables pour chacun des alignements générés pour les primates (2030 loci) et les oiseaux ( 854 loci). Nous avons sélectionné le modèle de substitution approprié pour tous les loci à l'aide de l'AIC (Akaike 1974).

Analyse bayésienne des données concaténées

Nous avons regroupé les gènes avec le même modèle de substitution dans différentes partitions, nous avons attribué un modèle de substitution approprié à chaque partition, et nous avons concaténé les partitions et analysé ces données à l'aide de MrBayes 3.1 ( Huelsenbeck et Ronquist 2001). Nous avons effectué toutes les analyses de MrBayes à l'aide de deux exécutions indépendantes (quatre chaînes chacune) de 5 000 000 d'itérations chacune, en échantillonnant des arbres toutes les 100 itérations pour produire un total de 50 000 arbres. Nous avons échantillonné les 25 000 derniers arbres après avoir vérifié la convergence des résultats en visualisant le log de la probabilité postérieure dans et entre les exécutions indépendantes pour chaque analyse, en veillant à ce que l'écart type moyen des fréquences fractionnées soit <0,0001, et en garantissant le facteur de réduction d'échelle potentiel pour l'estimation. paramètres était d'environ 1,0.

Analyse des arbres de gènes et des arbres d'espèces

Nous avons estimé les arbres génétiques sous le maximum de vraisemblance avec PhyML 3.0 (Guindon et al. 2010) en utilisant le modèle de substitution le plus probable pour chaque arbre, que nous avons estimé comme décrit ci-dessus. Nous avons estimé les arbres d'espèces à partir de ces arbres génétiques à l'aide des méthodes STAR (Species Trees based on Average Ranks of coalescences) et STEAC (Species Trees Estimated from Average Coalescent times) mises en œuvre dans le package R Phybase ( Liu et Yu 2010). STAR et STEAC calculent une topologie d'arbre d'espèces analytiquement basée sur les rangs moyens ou les temps d'événements de coalescence dans les collections d'arbres génétiques ( Liu et al. 2009). STAR et STEAC fonctionnent de manière similaire aux méthodes d'arbres d'espèces basées sur la coalescence probabiliste (par exemple, BEST), qui ne sont pas adaptées, d'un point de vue pratique, à la taille des ensembles de données utilisés ici. STAR fonctionne également bien lorsque les arbres de gènes s'écartent de taux d'évolution égaux, probablement le cas dans les phylogénies profondes et taxonomiquement diverses que nous avons étudiées. biaisé ( Liu et al. 2009). Après avoir généré un arbre à espèce unique, nous avons utilisé un programme Python personnalisé sur 250 nœuds d'un cluster Hadoop (http://hadoop.apache.org/) (Amazon Elastic Map Reduce) pour effectuer 1000 réplications bootstrap non paramétriques en rééchantillonnant les nucléotides dans les loci comme ainsi que le rééchantillonnage des loci dans l'ensemble de données ( Seo 2008).


Matériaux et méthodes

Définition des catégories d'inversion

Dans notre analyse, nous définissons trois classes de fréquence de population d'inversion. Travail antérieur dans D. melanogaster a généralement fait référence à quatre catégories d'inversion, « cosmopolite commune », « cosmopolite rare », « endémique récurrent » et « endémique unique » (Mettler et al. 1977 Krimbas et Powell 1992). La seconde moitié de chacun de ces termes fait référence à la distribution géographique de l'inversion. Tant qu'une inversion a atteint une fréquence élevée dans une population, elle n'a pas été fortement impactée par la sélection négative. Nous appelons ces inversions à haute fréquence des inversions « communes ». Nous utilisons « rare » pour désigner des inversions qui n'ont été trouvées que dans des échantillons uniques (à l'exception de In(2R)Mal, qui est présent dans trois échantillons étudiés ici). La distribution d'inversions rares, tout en contenant peut-être des inversions de haute aptitude qui pourraient éventuellement se propager à des fréquences élevées, est susceptible de refléter principalement des biais mutationnels dans leur distribution globale des points de rupture. Pour résumer, "commune cosmopolite", "rare cosmopolite" et "récurrente endémique" tomberont toutes sous notre étiquette "commune", alors que nous nous référons à "unique endémique" comme inversions "rares", de la même manière que l'analyse de Corbett-Detig. (2016).

La troisième classe dans notre cadre, les inversions « fixes », sont les inversions qui sont passées à la fixation au sein d'une lignée pendant la divergence de la D. melanogaster sous-groupe ( Ranz et al. 2007). A l'origine, toutes les inversions corrigées se produisaient comme des événements uniques dans un Drosophile ancêtre. Ils se sont ensuite propagés jusqu'à ce qu'ils atteignent la fixation dans les populations ancestrales des espèces contemporaines dans le melanogaster sous-groupe. Ces inversions fixes ont été découvertes en comparant les emplacements de séquences homologues dans les génomes entre D. melanogaster et ses apparentés ( Lemeunier et Ashburner 1976) et ont déjà été caractérisés moléculairement ( Ranz et al. 2007). Il est important de noter que la grande majorité de ces inversions fixes se sont produites sur le Drosophile yakuba branche et non en direct D. melanogaster ancêtre ( Krimbas et Powell 1992 Ranz et al. 2007). Le génome de référence de D. melanogaster devrait donc généralement refléter l'état ancestral et le fond génétique sur lesquels ces inversions ont pris naissance plutôt qu'un état dérivé évolué après fixation. Les inversions courantes et rares annotées ici se sont produites dans les D. melanogaster populations et donc en l'absence de modifications supplémentaires non liées à la structure du génome, sur un fond génétique similaire à celui sur lequel le D. yakuba les inversions ont été corrigées. Les annotations fonctionnelles utilisées ici sont également basées sur les D. melanogaster arrangement standard, ce qui signifie que ces annotations doivent représenter le fond génétique des trois catégories de fréquence d'inversion.

Alignement de lecture courte

Nous avons obtenu des données à lecture courte sous forme de fichiers fastq à partir de l'archive de lecture de séquence. Toutes les données à lecture courte sont décrites dans Lack et al. (2016) et a été initialement produit dans Pool et al. (2012), Lack et al. (2015), Mackay et al. (2012), Kao et al. (2015) et Grenier et al. (2015). Nous avons aligné les données de lecture courte à l'aide de bwa v0.7.15 à l'aide de la fonction « mem » et des paramètres par défaut ( Li 2013). Tous les post-traitements (tri, conversion au format BAM et filtrage) ont été effectués dans SAMtools v1.3.1 ( Li et al. 2009). Nous avons filtré ces fichiers BAM pour n'inclure que les alignements avec une qualité de mappage minimale de 20 ou plus.

Identification des points d'arrêt rares

Comme dans les travaux précédents qui caractérisaient la variation structurelle à l'aide de bibliothèques Illumina à insertion courte (Cridland et Thornton 2010 Rogers et al. 2014 Corbett-Detig et al. 2019), nous avons d'abord identifié des « clusters » de lecture mappés de manière aberrante. En bref, ici, un cluster est défini comme trois paires de lecture ou plus qui s'alignent dans la même orientation (pour les inversions, il s'agit à la fois d'un mappage vers l'avant ou à la fois d'un mappage inverse) et pour lesquelles toutes les lectures à un bord du cluster correspondent à à moins de 1 ko de toutes les autres lectures dans le cluster. Nous n'avons considéré que les clusters aberrants où les deux extrémités étaient mappées sur le même bras chromosomique que la grande majorité des inversions dans Drosophile sont paracentriques ( Krimbas et Powell 1992). Nous avons exigé que toutes les paires de lecture incluses dans une carte de cluster soient espacées d'au moins 500 ko. Nous n'avons ensuite retenu que les inversions potentielles pour lesquelles nous avons récupéré à la fois des clusters de cartographie directe et inverse qui se trouvaient à moins de 100 kb les uns des autres. Le choix d'une distance maximale entre les coordonnées possibles des points d'arrêt a été inclus pour réduire les taux possibles de faux positifs et parce qu'aucune des inversions connues dont les points d'arrêt ont déjà été caractérisés n'incluait une région dupliquée de 100 kb ou plus ( Ranz et al. 2007 Corbett -Detig et Hartl 2012). Lorsque des assemblages de points d'arrêt existaient à proximité immédiate ou semblaient supprimer de courtes séquences, nous avons défini la taille de duplication sur 1 base. Nous avons en outre filtré tous les assemblages de points d'arrêt qui chevauchaient des éléments transposables d'annotation, car ils sont la principale source de mappage aberrant des clusters de lecture dans des travaux antérieurs ( Corbett-Detig et Hartl 2012).

Comme vérification supplémentaire de l'exactitude de nos points d'arrêt nouvellement découverts, nous avons comparé notre distribution de points d'arrêt rares à la distribution cytogénétique connue et n'avons trouvé aucune différence chromosomique ou par région (P = 0,7, χ 2 testent les données cytogénétiques de Corbett-Detig (2016) qui a résumé Krimbas et Powell (1992)). La taille de l'insert court des expériences de séquençage précédentes variait de ∼200 à ∼600 pb, ce qui peut avoir conduit à un taux de faux négatifs non négligeable de découverte de points d'arrêt, en particulier si les points d'arrêt contiennent des éléments répétitifs ou d'autres grandes insertions d'ADN. Cependant, nous ne nous attendons pas à ce que ces faux négatifs potentiels biaisent nos analyses en aval, et tous les points de rupture d'inversion précédemment caractérisés dans le Mélanogastre complexe d'espèces se sont produits dans des séquences uniques ( Ranz et al. 2007 Corbett-Detig et Hartl 2012). Tous les logiciels utilisés pour effectuer ces analyses sont disponibles depuis les dépôts github associés à ce projet. Plus précisément, les scripts utilisés pour la détection et l'assemblage des points d'arrêt se trouvent sur https://github.com/dliang5/breakpoint-assembly (dernière consultation le 26 mai 2020).

Assemblage de points d'arrêt rares De Novo

Pour chaque inversion putative, nous avons ensuite extrait toutes les lectures pour lesquelles l'une ou l'autre paire correspondait à moins de 5 ko de la position prédite du point d'arrêt. Nous avons converti tous les fichiers de lecture fastq en fichiers fasta et qual comme requis par Phrap, et nous avons assemblé chacun en utilisant des paramètres par défaut, mais en incluant les options de ligne de commande « -vector_bound 0 -forcelevel 10 » ( Corbett-Detig et Hartl 2012 Rogers et al. 2014). Nous avons ensuite utilisé BLAST pour aligner les contigs assemblés de novo résultants sur le D. melanogaster génome de référence pour identifier le contig qui chevauchait le point d'arrêt prévu à l'aide de l'outil flybase BLAST (https://flybase.org/blast/, consulté pour la dernière fois le 26 mai 2020). Nous n'avons retenu que les inversions pour lesquelles nous pouvions de novo assembler des contigs chevauchant les deux points de rupture, et nous avons en outre rejeté tous les contigs où la séquence intervenant dans deux régions génomiques distantes contenait une séquence avec une homologie avec des éléments transposables connus. Toutes les séquences de points d'arrêt assemblées sont disponibles dans le fichier supplémentaire S1 , Supplementary Material en ligne. Les scripts d'assemblage sont disponibles sur https://github.com/dliang5/breakpoint-assembly (dernière consultation le 26 mai 2020).

Chevauchement des inversions et In(2R)Mal

Nous avons également tenté de trouver des ensembles d'inversions qui se chevauchent. En bref, pour les inversions qui se chevauchent, où une inversion apparaît sur un arrière-plan qui contient une autre inversion avec un point d'arrêt à l'intérieur et un à l'extérieur de la région inversée, les clusters de lecture couvrant les points d'arrêt devraient être en grande partie les mêmes que les inversions apparues sur un chromosome d'arrangement standard. Cependant, la principale différence est que plutôt que des paires de clusters de lecture à mappage direct et inverse, nous nous attendons à observer deux clusters de lecture à mappage distant dans les arrangements inverse-forward et forward-reverse. Nous avons appliqué cette approche pour les 17 inversions rares que nous avons initialement découvertes ainsi qu'à tous les échantillons qui contenaient des inversions communes connues de travaux antérieurs ( Corbett-Detig et Hartl 2012 Lack et al. 2015). Nous n'avons trouvé qu'une seule inversion rare qui se chevauche, ce qui est cohérent avec l'inversion chromosomique connue associée à la distorsion de ségrégation. In(2R)Mal, qui est composé de deux inversions superposées ( Presgraves et al. 2009). Dans notre analyse ici, nous traitons ces inversions qui se chevauchent comme indépendantes, mais nos résultats ne sont qualitativement pas affectés si nous excluons simplement la deuxième inversion.

Version du génome, isolant et annotations génétiques

Toutes nos analyses sont basées sur des alignements sur D. melanogaster génome version 6.26 ( Hoskins et al. 2015). Nous avons obtenu des données d'annotation du génome, y compris la localisation des gènes à partir de flybase. Nous avons traité les ARN longs non codants comme des gènes pour nos besoins, car ils remplissent des fonctions essentielles et peuvent être perturbés de la même manière que les gènes codant pour des protéines. Nous avons obtenu les positions des sites de liaison de l'isolant de Nègre et al. (2010, accès GSE16245). Si nécessaire, nous avons converti les coordonnées des caractéristiques génomiques de la version 5 du génome à la version 6 à l'aide de l'outil de conversion par lots de coordonnées flybase (https://flybase.org/convert/coordinates, consulté pour la dernière fois le 26 mai 2020).

Sélection d'ensembles de données publiques pour les domaines topologiques et les marques de chromatine

Nous avons obtenu des données TAD comprenant des annotations de l'état de la chromatine de Sexton et al. (2012). Cet ensemble de données est composé de domaines détectés par séquençage de capture de conformation chromosomique à l'échelle du génome, HiC, sur des embryons à un stade précoce, et annotés avec un état épigénétique à l'aide d'une méthode de regroupement appliquée à une autre source de données épigénomiques linéaires (Sexton et al. 2012). Leurs annotations comprennent quatre catégories : « actif », « nulle », « PcG » (polycomb) et « HP1 » (hétérochromatine centromérique). Par souci de cohérence, nous appelons les domaines « null » de Sexton et al. « inactifs ». Les embryons à un stade précoce sont susceptibles d'être l'environnement dans lequel toute perturbation de la régulation induite par les inversions est la plus délétère étant donné la nature sensible du développement, ce qui en fait une source prometteuse de contexte pour notre analyse de la fréquence d'inversion. Cet ensemble de données nous permet également d'analyser séparément l'occurrence des points de rupture dans les TAD et les états de la chromatine en tandem, car ils sont dérivés de la même source. Il convient toutefois de noter que les annotations de ces TAD sont relativement grossières et peuvent ne pas refléter l'environnement plus local d'un point de rupture d'inversion.

Nous avons donc effectué une deuxième analyse à des échelles plus fines en utilisant le jeu de données de Kharchenko et al. (2011, accès GSE25321). Cet ensemble de données sous sa forme brute se compose de courtes portées marquées par l'un d'un ensemble de marqueurs de chromatine, à la fois dans un modèle à neuf états et dans un modèle à 30 états. Comme nous souhaitions une représentation de l'environnement chromatinien local autour des points de rupture d'inversion, nous avons choisi de regrouper la représentation à neuf états en nombres totaux de bases affectées à un état du type donné sur des fenêtres de 10 kb. Environ 10 ko ont été sélectionnés en fonction de l'hétérogénéité moyenne des fenêtres. Nous voulions que notre taille de fenêtre soit aussi petite que possible, mais que la plupart des fenêtres contiennent au moins une région avec un état de chromatine annoté. Cela a donné une distribution de valeurs pour chaque fenêtre qui représentait l'enrichissement global de chaque état dans chaque tranche de 10 ko. Comme nous manquions de puissance statistique pour évaluer ces types de marques individuellement avec nos ensembles de données de points de rupture d'inversion relativement petits, nous avons en outre attribué à chaque fenêtre de 10 ko un état d'activité basé sur la majorité des marques présentes. Les fenêtres dans lesquelles la grande majorité des sites ont reçu les états un à cinq, annotés par Kharchenko et al. (2011) comme étant divers composants des gènes, y compris les promoteurs, les exons et les introns, ont été désignés « actifs ». Les fenêtres où les états six à neuf, qui incluent PcG, HP1 et d'autres marques hétérochromatiques, étaient les plus importants, ont été désignées « inactives ». Les fenêtres dans lesquelles les deux groupes constituaient chacun au moins 5 % de toutes les marques ont été désignées « mixtes ». Cela donne une représentation alternative des environnements de chromatine entourant les points de rupture d'inversion qui est beaucoup plus fine que les annotations de Sexton et al. (2012).

Nous avons comparé cette représentation aux états de chromatine annotés de Sexton et al. comme vérification supplémentaire de la validité de notre approche. Nous avons constaté que les fenêtres de 10 ko situées dans chaque TAD annoté étaient généralement alignées sur l'annotation de ce TAD, mais qu'une hétérogénéité substantielle des marques de chromatine existe dans chaque étendue de TAD (fig. supplémentaire S3, matériel supplémentaire en ligne). Par exemple, ∼19 % des fenêtres dans les TAD annotées comme « actives » sont enrichies pour l'état de chromatine 9, qui est associé à des régions silencieuses étendues, et inversement, 26 % des fenêtres dans les TAD annotées comme inactives sont enrichies pour l'état de chromatine 2, associé à la transcription active. Cela indique que l'un ne peut pas être traité comme un substitut direct de l'autre.

En guise de vérification finale de la validité des domaines obtenus de Sexton et al. (2012), nous avons obtenu des données sur le domaine polytène d'Eagen et al. (2015), ont répété notre analyse et les ont trouvés globalement cohérents avec nos conclusions. Ces résultats peuvent être trouvés dans le texte supplémentaire S1 , Matériel supplémentaire en ligne.

Permutations et tests statistiques

Pour comparer les positions des points d'arrêt d'inversion à une distribution aléatoire, des permutations pour toutes les catégories d'inversions (rares, communes et fixes) ont été effectuées avec 1 000 itérations d'un groupe de points d'arrêt situés au hasard, en maintenant le nombre d'inversion, les longueurs de duplication et les bras chromosomiques constants. Plus précisément, pour chaque point de rupture d'inversion, 1 000 positions de départ ont été choisies à partir d'une distribution uniforme entre le début de ce bras chromosomique et la fin moins la longueur de la duplication, c'est-à-dire à partir de l'ensemble des points possibles pour cette taille de point de rupture. Des points de rupture aléatoires ont été localisés indépendamment pour la plupart des tests, car la plupart des valeurs ont été calculées pour chaque point de rupture individuellement plutôt que pour l'inversion dans son ensemble. L'exception est le test de mélange de chromatine, dans lequel nous avons en outre contrôlé les longueurs d'inversion pour tenir compte du rôle de la longueur d'inversion dans la polarisation des paires d'environnements de chromatine. Les caractéristiques du génome à chacun de ces points d'arrêt ont été enregistrées comme notre valeur attendue pour la distribution aléatoire des points d'arrêt.

Les tests ont été divisés selon la nature du facteur. Pour les facteurs qui sont une valeur numérique discrète pour chaque rupture, comme la distance à un élément ou la longueur d'une duplication, P les valeurs ont été calculées en tant que centiles des valeurs réelles au sein d'un large ensemble de distributions aléatoires. Des tests entre les catégories de facteurs basés sur la distance et le test de longueur de duplication ont été effectués distribution à distribution avec des tests de somme de rang Mann-Whitney par paires.

Pour les valeurs catégorielles, telles que la perturbation ou non d'une étendue de gènes, les taux d'occurrence de catégorie ont été calculés pour 1 000 permutations. Nous définissons les perturbations des gènes et d'autres éléments comme des ruptures simple brin directes et inverses se produisant au sein d'un élément fonctionnel à annotation unique. Il est important de noter que notre méthode de définition de la perturbation est susceptible de surestimer la proportion de points de rupture d'inversion fixes qui perturbent réellement les séquences géniques.La méthode de Ranz et al. (2007) pour identifier les séquences dupliquées par la rupture d'origine repose sur l'homologie de séquence, et dans les inversions fixes, la divergence des séquences non codantes peut interférer avec l'identification précise des régions des points de rupture. Par exemple, si la région dupliquée d'origine comprend une étendue de gène codant et quelques bases non codantes, une copie complète du gène sera produite avec une duplication partielle. Au fil du temps, la région non codante aura tendance à accumuler plus de mutations que la copie intacte du gène. Dans ce cas, les coordonnées obtenues à partir des alignements BLAST peuvent ne pas détecter l'homologie entre les régions non codantes et à la place ne donner qu'une homologie apparente à partir de la duplication dans l'étendue du gène conservé. Cela serait compté comme un événement de perturbation génétique par notre analyse. Ce biais aura tendance à rendre notre analyse prudente en ce qui concerne l'identification des impacts de la sélection naturelle, car les points de rupture sont plus susceptibles d'être identifiés dans les régions codantes et parce que nous devrions avoir tendance à sous-estimer les tailles des régions dupliquées adjacentes aux points de rupture après que l'homologie de séquence a diminué. . Tous les scripts utilisés pour produire les résultats des tests de permutation décrits ci-dessus sont disponibles dans le référentiel github associé à ce projet https://github.com/jmcbroome/breakpoint_analysis (dernière consultation le 26 mai 2020).

Analyse de phénotype létal et stérile

De plus, nous avons obtenu des données de phénotype de Flybase à l'aide du générateur de requêtes (https://flybase.org/cgi-bin/qb.pl, consulté pour la dernière fois le 26 mai 2020) pour obtenir les identifiants de tous les gènes qui ont des phénotypes mortels et des phénotypes stériles. . Ces données ont été incorporées dans l'analyse de perturbation des gènes et nous avons recherché des preuves de différence dans les taux de perturbation entre les gènes annotés avec ces phénotypes et l'ensemble global de gènes annotés. Le tableau supplémentaire S2, Supplementary Material en ligne, contient l'ensemble des points de rupture d'inversion qui semblent perturber ces gènes.


RÉSULTATS

Détection d'épissage alternatif

Notre analyse de l'épissage alternatif est basée strictement sur des données expérimentales, et non sur des modèles théoriques. Plutôt que de chercher à prédire des épissures alternatives, nous les détectons directement comme de grands inserts dans les données EST des bases de données dbEST ( 20) et UNIGENE (18) accessibles au public. Nous mesurons la preuve d'une véritable épissure alternative via une série de critères (Fig. 1). Tout d'abord, un ensemble d'EST doit correspondre sur toute leur longueur, des deux côtés d'une épissure alternative putative (permettant une erreur de séquence). Un grand insert au milieu d'une telle correspondance parfaite est une épissure alternative candidate. Contrairement à de nombreux autres types de résultats génomiques tels que les SNP et les variations du niveau d'expression, l'épissage alternatif ne ressemble pas au bruit expérimental commun (comme l'erreur de séquençage).

Ensuite, la séquence consensus EST est mappée sur la séquence de projet du génome humain par recherche d'homologie. Étant donné que les gènes humains sont divisés en exons courts, un hit génomique consiste généralement en de nombreuses correspondances courtes. Pour être valides, ces correspondances doivent être parfaites (en ne permettant qu'une erreur de séquençage), doivent toutes être dans la même orientation (brin) et former une marche complète et correctement ordonnée à travers la séquence consensus EST. Nous exigeons que chaque région de correspondance génomique-EST (exon putatif) soit délimitée par des séquences consensus de site donneur d'épissage et de site accepteur dans la séquence génomique voisine (intron). Nos résultats donnent une taille moyenne d'exon interne de 144 pb, avec seulement 4% d'exons internes >300 pb de longueur, similaire aux résultats obtenus pour des gènes connus (21). Seuls 0,2 % (79/39 862) de nos introns étaient inférieurs à 60 pb, et la longueur médiane des introns était de 935 pb. Le modèle génique typique des exons internes courts se terminant par un seul exon 3' long peut généralement être vérifié car les séquences d'extrémité 3' sont fortement représentées dans les données EST, et parce que les EST 3' peuvent être identifiés par leur poly(A) bien visible. queues, qui indiquent directement la fin de l'exon 3′.

Pour évaluer l'exactitude de notre cartographie génétique et de notre structure exon/intron, nous avons comparé les données complètement indépendantes produites par Acembly du NCBI, un effort d'annotation de gènes organisés par l'homme (données téléchargées à partir de ftp://ncbi.nlm.nih.gov/genomes /H_sapiens). LocusLink fournit un lien indépendant entre les gènes individuels RefSeq et les clusters UNIGENE (22). Pour les gènes cartographiés indépendamment sur la séquence génomique par RefSeq et notre procédure, 97,3 % cartographiés sur le même contig génomique. De plus, parmi ces gènes, 95% ont été mappés sur les mêmes nucléotides du contig. Bien que la cartographie d'Acembly ne doive pas être considérée comme parfaite, ce niveau élevé d'accord entre les efforts indépendants est encourageant. Nos détails d'exons (dérivés dans notre procédure de notre détection d'épissage) correspondent aux exons NCBI Acembly dans 97% des cas au site d'épissage 5' et 96% au site d'épissage 3' (dans l'ensemble, 94% des exons étaient identiques) . Nos détails d'épissage correspondaient aux introns NCBI Acembly dans 93% des cas à l'extrémité 5' et 92% à l'extrémité 3' (86 % correspondaient exactement aux deux extrémités). En raison de l'épissage alternatif, une correspondance à 100 % n'est pas attendue.

Un insert d'épissure alternatif candidat (de l'EST) doit réussir une série de tests. Premièrement, il doit également être trouvé dans la séquence génomique, correspondant à une région exonique dans la séquence génomique dont les limites correspondent à des séquences de sites d'épissage connues. Étant donné que ces séquences de sites d'épissage sont pour la plupart introniques, cela fournit une validation indépendante de l'épissage alternatif. Il convient de souligner que les différences d'endroit où les EST commencent et se terminent dans un gène (par exemple, un EST plus court pourrait donner l'apparence d'un produit de gène tronqué) ne seront jamais interprétées comme une épissure alternative par notre procédure. Nous nous concentrons exclusivement sur la détection de l'épissage, c'est-à-dire une région contiguë du transcrit qui a été supprimée lors du traitement de l'ARNm. La détection d'une épissure dans un EST nécessite des correspondances étendues avec les exons en amont et en aval. Notre analyse a identifié 39 862 épissures dans 8429 grappes. Notre analyse ne rapporte que des épissures alternatives, c'est-à-dire des paires d'épissures validées qui s'excluent mutuellement. Ainsi, les introns non épissés ou d'autres contaminants génomiques ne seront jamais signalés, car ils entraînent l'absence d'épissage, et non la création d'un nouvel épissage mutuellement exclusif. Pour appeler une épissure alternative, notre procédure nécessite une paire d'épissures qui correspondent exactement à un site d'épissure et qui diffèrent au niveau du second site d'épissure. Cette procédure peut détecter le saut d'exon, les sites donneurs alternatifs 5' et les sites accepteurs alternatifs 3' (Fig. 1B). 6201 de telles relations d'épissage alternatives ont été identifiées dans 2272 grappes. Ces diverses formes de preuves produisent des scores logarithmiques élevés pour chaque épissure alternative détectée. Une analyse statistique détaillée de ces preuves sera présentée ailleurs (D.Miller, J.Aten, C.Grasso, B.Modrek et C.Lee, manuscrit en préparation).

Comme exemple de validation typique de notre base de données, nous illustrons le gène de la protéine kinase de dystrophie myotonica (DMPK) (Fig. 2), dont l'épissage alternatif a déjà été largement étudié. Dans DMPK, nous avons identifié trois épissures alternatives dans les données EST, qui sont toutes vérifiées par des résultats expérimentaux indépendants dans la littérature existante (23). Sur les trois épissures alternatives, l'une supprime les 15 derniers pb de l'exon 8, une autre saute l'exon 12 et l'exon 13, et la dernière supprime seulement 4 pb de l'exon 14. La figure 2 montre l'une de ces formes d'épissage alternatives, y compris la jonction et la qualité de la correspondance. de la preuve EST par rapport à la séquence génomique.

Nouvelles formes d'épissage alternatif d'un gène connu

La figure 3 montre plusieurs nouvelles épissures alternatives détectées dans un gène bien étudié, HLA-DM . Quatre-vingts EST du cluster UNIGENE Hs.1162 s'alignent pour former une séquence consensus, qui à son tour correspond à une série ordonnée de segments sur un brin du chromosome 6. Les séquences EST correspondent étroitement à la séquence génomique, ce qui correspond à une erreur de séquençage. Les séquences EST délimitent un exon 3' long (359 pb) plus une série de cinq exons courts, dont les tailles (36-288 pb) correspondent à la plage attendue pour les exons internes. Cela correspond à l'emplacement et à la structure connus du gène pour HLA-DM ( 24, 25). Huit épissures sont observées dans ces EST, où la séquence correspondant à une région exonique passe directement à une région exonique en aval, comme indiqué sur la figure 3A. Les 16 limites d'exons putatives impliquées par les EST correspondent précisément à des sites accepteurs et donneurs d'épissage consensuels forts dans la séquence génomique (Fig. 3C).

Quatre formes différentes d'épissage alternatif de HLA-DM β sont observées : épissures 3+4+5 (incluant les exons IV et V dans le produit d'ARNm) épissures 6+5 (en sautant l'exon IV) épissures 3+7 (excluant l'exon V) épissure 8 (en sautant les exons IV et V). Les exons IV et V ont une longueur de 117 et 36 pb, et donc ces épissures alternatives sont toutes dans le cadre. La région codant pour la protéine commence dans l'exon I et se termine dans l'exon VI, de sorte que ces épissures produisent quatre formes différentes de la chaîne β HLA-DM qui diffèrent à leur extrémité C-terminale.

L'analyse de ces formes révèle un effet fonctionnel remarquablement simple et intrigant. HLA-DM est essentiel pour le chargement des molécules du CMH de classe II avec des antigènes peptidiques exogènes, une étape clé dans la présentation de l'antigène et l'activation de la réponse immunitaire humorale. On pense que cela se produit dans les premiers compartiments lysosomal. HLA-DM est normalement ciblé sur les lysosomes et sa chaîne contient un domaine transmembranaire ancrant son extrémité C (26, 27). L'exon IV est court et correspond précisément au domaine transmembranaire. L'exon V est très court et code pour le signal de ciblage lysosomal YTPL, dont le premier résidu commence au début de l'exon. Ainsi, l'épissure alternative régule le ciblage de HLA-DM vers les compartiments endosomaux (en incluant ou excluant le signal YTPL), ainsi que son ancrage à la membrane. Compte tenu de l'importance de HLA-DM dans le traitement et la présentation des antigènes par le CMH de classe II, cette régulation est fonctionnellement intéressante. La suppression de son signal de ciblage redirigerait probablement le HLA-DM d'abord vers la membrane plasmique, de sorte qu'il se déplacerait vers les lysosomes via des voies endocytiques, modifiant la cinétique et les conditions dans lesquelles il rencontre pour la première fois le CMH de classe II. Il semble que la structure des gènes du HLA-DM Le gène β a été soigneusement « conçu » pour permettre le contrôle de la fonction HLA-DM, en extrayant à la fois l'hélice transmembranaire et le signal de ciblage lysosomal dans des exons courts séparés (IV, V) qui peuvent être épissés alternativement dans le cadre (l'exon VI fournit les 4 derniers acides aminés de la protéine, identiques sous toutes leurs formes). Les formes épissées alternativement ont été détectées dans l'utérus (deux EST), le placenta, la lymphe, l'estomac et le côlon. Malgré le fait que le HLA-DM fasse l'objet de recherches intenses, nous n'avons pu trouver aucun rapport sur un tel épissage alternatif dans la littérature publiée, et il est considéré comme nouveau par un expert du HLA-DM (E.Mellins , communication personnelle).

Portée de l'épissage alternatif dans les gènes humains

Notre analyse à l'échelle du génome a détecté des milliers d'épissages alternatifs dans les données actuelles du génome humain accessibles au public (tableau 1). 6201 relations d'épissage alternatif ont été détectées dans lesquelles deux épissures partageaient un site donneur ou accepteur commun, mais épissées à un site différent à leur autre extrémité (c. . Nous avons trouvé des épissures alternatives dans 27% des gènes pour lesquels nous avions suffisamment de séquences exprimées pour couvrir plus d'un seul exon. Cependant, cette estimation, basée sur l'analyse de tous les clusters EST, sous-estime probablement l'occurrence réelle de l'épissage alternatif, car les données EST disponibles ne couvrent généralement qu'une petite partie du gène complet. Pour tester cette hypothèse, nous avons analysé le taux d'épissage alternatif dans les gènes pour lesquels la séquence d'ARNm était disponible (représentant tout ou partie du gène complet). Nous avons détecté une ou plusieurs formes d'épissage alternatives dans 42 % de ces gènes, ce qui est significativement plus élevé que le taux observé dans les clusters EST uniquement. Ceci est en accord étroit avec une étude précédente de clusters de séquences exprimées à base d'ARNm (8). Étant donné que la fragmentation de la séquence génomique peut également bloquer la couverture complète d'un gène, nous avons évalué le taux de détection d'épissage alternatif dans les gènes mappés sur le chromosome 22. Parmi ceux-ci, 43 % contenaient des épissures alternatives, y compris des clusters d'ARNm et d'EST uniquement.

Les données EST actuelles semblent incomplètes. Notre procédure a identifié des épissures (c'est-à-dire des exons multiples) dans seulement 18% des clusters EST cartographiés. Cependant, pour les clusters que nous avons cartographiés sur le chromosome 22 (génomique complet) qui avaient également une séquence d'ARNm, 88 % contenaient au moins une épissure. Une variété de facteurs tels que la fragmentation du projet de séquence du génome humain, la grande taille des introns et la tendance des EST à se regrouper à l'extrémité 3' faussent l'ensemble de données actuel contre la recherche de gènes de pleine longueur, et sous-estiment probablement le vrai niveau d'épissage alternatif. De plus, étant donné que les données EST actuelles pour chaque gène ne représentent qu'un sous-ensemble des tissus et des types cellulaires dans lesquels ce gène est exprimé, il est probable que l'occurrence totale d'épissage alternatif soit bien supérieure à ce que nos analyses peuvent détecter. Une grande fraction des formes d'épissage alternatives EST a été observée plusieurs fois (à partir de différents clones et différentes bibliothèques), indiquant qu'elles constituent une fraction relativement élevée de l'ARNm total. Parmi nos épissures alternatives, 2892 (47 %) ont été observées dans deux séquences EST ou plus. Ces données représentent un sous-ensemble « de confiance élevée » des épissures alternatives détectées.

Notre analyse indique que la grande majorité de notre base de données représente de nouveaux résultats (Fig. 5D). Seulement 13% de nos épissures alternatives ont été détectées dans des séquences d'ARNm de GenBank, qui ont probablement été étudiées de manière approfondie. Les 87 % restants n'ont pu être détectés qu'avec les EST. Notre procédure a également détecté un grand nombre d'événements d'épissage alternatif dans des gènes complètement nouveaux. Environ 1200 épissures alternatives ont été détectées dans des grappes contenant uniquement des EST.

Épissage alternatif dans un nouveau gène humain

La figure 4 illustre un exemple de détection d'épissage alternatif dans un nouveau gène cartographié dans le génome humain par notre procédure. Ce gène a 33 % d'identité avec la chaîne β du récepteur FCε de rat et 25 % d'identité avec CD20, et a un modèle de quatre domaines transmembranaires prédits caractéristiques des deux protéines. Au moins sept formes différentes sont détectables, qui affectent toutes le produit protéique. Dans un modèle qui rappelle de façon frappante HLA-DM β, la région transmembranaire C-terminale et la queue cytoplasmique de la forme principale (forme 1) sont placées sur un seul exon court (exon VI), qui peut être inclus ou exclu pour créer différentes formes. Une forme particulièrement intéressante est créée en ignorant l'épissage normal de l'exon V à l'exon VI, en étendant la région codante de l'exon Va pour 142 pb (que nous avons désigné exon Vb). Un site de polyadénylation est prédit à la fin de cette séquence, et on observe que les EST se terminent en poly(A) à ce stade. Cette terminaison alternative remplace la région codante des exons VI et VII par 40 acides aminés codés par l'exon Vb [terminé par un codon STOP 23 pb avant le site poly(A)]. Curieusement, cette séquence C-terminale de remplacement contient également une séquence transmembranaire prédite, et remplace ainsi parfaitement un nouveau domaine transmembranaire C-terminal et une queue cytoplasmique. La queue cytoplasmique dans les chaînes de récepteur FC équivalentes joue un rôle clé dans l'activation des molécules de transduction du signal cytoplasmique (28, 29), de sorte que cette forme alternative module probablement l'activité de transduction du signal de ce récepteur. Cette forme est détectée dans le placenta et les reins, tandis que la forme majoritaire a été détectée dans de nombreuses bibliothèques différentes.


Résumé

Le séquençage à molécule unique en temps réel développé par Pacific BioSciences offre des longueurs de lecture plus longues que les technologies de séquençage de deuxième génération (SGS), ce qui le rend bien adapté aux problèmes non résolus dans la recherche sur le génome, le transcriptome et l'épigénétique. Le très contigu de novo assemblages l'utilisation du séquençage PacBio peut combler les lacunes dans les assemblages de référence actuels et caractériser la variation structurelle (VS) dans les génomes personnels. Avec des lectures plus longues, nous pouvons séquencer à travers des régions répétitives étendues et détecter des mutations, dont beaucoup sont associées à des maladies. De plus, le séquençage du transcriptome PacBio est avantageux pour l'identification des isoformes de gènes et facilite des découvertes fiables de nouveaux gènes et de nouvelles isoformes de gènes annotés, en raison de sa capacité à séquencer des transcrits complets ou des fragments de longueurs significatives. De plus, la technique de séquençage de PacBio fournit des informations utiles pour la détection directe des modifications de base, telles que méthylation. En plus d'utiliser le séquençage PacBio seul, de nombreux séquençage hybride des stratégies ont été développées pour utiliser des lectures courtes plus précises en conjonction avec des lectures longues PacBio. En général, séquençage hybride les stratégies sont plus abordables et évolutives, en particulier pour les laboratoires de petite taille, que d'utiliser le séquençage PacBio seul. L'avènement du séquençage PacBio a rendu disponible de nombreuses informations qui ne pouvaient pas être obtenues via SGS seul.


Voir la vidéo: Quelles Sont Les Principales Religions Du monde?-Zakir Naik (Janvier 2022).