Informations

Trop peu de transcrits de l'assembleur de transcriptome Oases


J'essaie d'exécuter Oases pour l'assemblage du transcriptome. Le résultat est loin d'être attendu, alors je voudrais demander si je l'exécute correctement ? Merci.

Voici ma commande en cours d'exécution :

python scripts/oases_pipeline.py -m 25 -M 29 -o output -d " -strand_specific -shortPaired data/reads.fa" -p " -min_trans_lgth 100 -ins_length 300"

Ma bibliothèque est spécifique au brin et se termine par une paire avec une longueur de 67 pb. Les lectures sont mélangées comme :

>0(left_mate_forwarded) ACTC… >1(right_mate_reverse_complemented) TATA…

J'ai quelques transcriptions, mais je suis loin des transcriptions annotées, loin aussi du résultat de Trinity. Le contig le plus long d'Oases est d'environ 2500 pb (contre ~ 10000 pb pour les boutons de manchette et ~ 6000 pb pour Trinity). La valeur N50 est également faible. Il ne rapporte que 20 contigs qui couvrent toute la longueur de certaines transcriptions de Cufflinks (totalement ~ 4000), tandis que Trinity rapporte ~ 650.

L'ensemble de données que j'utilise est un sous-ensemble de S. pombe. Est-ce que ça importe?

Quelqu'un pourrait-il m'aider à indiquer si quelque chose ne va pas ici?


TraRECo : une approche gourmande basée sur un assembleur de transcriptome de novo avec correction d'erreur de lecture à l'aide d'une matrice de consensus

Les défis lors du développement d'un bon assembleur de transcriptome de novo comprennent la façon de traiter les erreurs de lecture et les répétitions de séquence. Presque tous les assembleurs de novo utilisent un graphe de Bruijn, avec lequel la complexité augmente linéairement avec la taille des données tout en souffrant d'erreurs et de répétitions. Bien que l'on puisse corriger les erreurs en inspectant la structure topologique du graphe, ce n'est pas une tâche facile lorsqu'il y a trop de branches. Deux directions de recherche visent à améliorer soit la fiabilité du graphe, soit la précision de la recherche de chemin, et dans cette étude, nous nous sommes concentrés sur la première.

Résultats

Nous présentons TraRECo, une approche gourmande de l'assemblage de novo utilisant la construction de graphes tenant compte des erreurs. Dans l'approche proposée, nous avons construit des contigs par alignement de lecture directe dans une marge de distance et effectué une recherche de jonction pour construire des graphes d'épissage. Ce faisant, un contig de longueur je était représenté par un 4 × je matrice (appelée matrice de consensus), dans laquelle chaque élément était le nombre de base des lectures alignées jusqu'à présent. Une séquence représentative a été obtenue en prenant la majorité dans chaque colonne de la matrice consensus à utiliser pour un alignement de lecture supplémentaire. Une fois les graphes d'épissage obtenus, nous avons utilisé IsoLasso pour trouver des chemins avec une profondeur de lecture notable. Les expériences utilisant des lectures réelles et simulées montrent que la méthode a permis une amélioration considérable de la sensibilité et des performances modérément meilleures lors de la comparaison de la sensibilité et de la précision. Ceci a été réalisé par la construction de graphe consciente des erreurs à l'aide de la matrice de consensus, avec laquelle les lectures comportant des erreurs ont été rendues utilisables pour la construction de graphe (sinon, elles auraient pu être finalement rejetées). Cela a amélioré la qualité des informations de profondeur de couverture utilisées dans l'étape de recherche de chemin suivante et enfin la fiabilité du graphe.

Conclusion

L'assemblage de novo est principalement utilisé pour explorer des isoformes non découvertes et doit être capable de représenter autant de lectures que possible de manière efficace. En ce sens, TraRECo nous offre une alternative potentielle pour améliorer la fiabilité des graphes même si la charge de calcul est beaucoup plus élevée que le seul k-mer dans l'approche du graphe de de Bruijn.


Assemblage du transcriptome de novo

Résumé
Contexte Le kiwi [Actinidia deliciosa (A Chev) Liang et Ferguson] est une vigne subtropicale de la famille des Actinidiaceae originaire de Chine. Cette espèce possède un génome allohexaploïde (issu d'un parent diploïde et autotétraploïde) contenu dans 174 chromosomes produisant un fruit climatérique et charnu appelé kiwi. Actuellement, il n'y a pas trop d'informations génomiques et transcriptomiques sur cette espèce. Dans ce contexte de faible connaissance moléculaire, l'objectif principal de ce travail est de construire un assemblage de transcriptome de novo tissu-spécifique générant une analyse d'expression différentielle entre ces tissus spécifiques pour obtenir une nouvelle base de données utile pour une meilleure connaissance de la croissance végétative, florale et fruitière dans différents états phénologiques d'Actinidia deliciosa cv. 'Hayward'.

Résultats Dans la présente étude, nous avons analysé différents transcriptomes entiers de pousses, feuilles, boutons floraux, fleurs et fruits à 4 stades de développement (7, 50, 120 et 160 jours après la floraison DAF) dans le kiwi en utilisant RNA-seq. Nous avons séquencé vingt-quatre bibliothèques, obtenant 604 735 364 lectures qui ont été assemblées à l'aide du logiciel Trinity. La première version du transcriptome d'Actinidia deliciosa de novo contenait 142 025 contigs (x̅ = 1 044 pb, N50 = 1 133 pb). CEGMA et BUSCO ont été utilisés pour l'évaluation de la qualité des assemblages, obtenant près de 90,0% (35,1% partiel) et plus de 85,0% (18,3% partiel) des gènes ultraconservés pour les eucaryotes et les plantes, respectivement. L'annotation a été effectuée avec BLASTx par rapport à la base de données de protéines TAIR10 et nous avons trouvé une proportion d'annotation de 35,6% (50 508), laissant 64,4% (91 517) de l'assemblage de contigs sans annotation.

Conclusions Ces résultats représentent un transcriptome de référence pour les kiwis allohexaploïdes générant une base de données de gènes d'Actinidia deliciosa liés au développement des feuilles, des fleurs et des fruits. Ainsi, la présente étude fournit des informations précieuses, identifiant plus de 20 000 gènes exclusifs, y compris toutes les comparaisons de tissus, qui sont associés aux protéines impliquées dans différents processus biologiques et fonctions moléculaires. L'assemblage et le raffinage du transcriptome ainsi que l'évaluation de la métrique d'assemblage ont impliqué une qualité suffisante pour constituer une base de données putative de cette espèce et un grand nombre de protéines ultra-conservées ont été trouvées. En ce qui concerne le transcriptome, près de 65 % des contigs ne correspondaient à aucune protéine. Par conséquent, une future annotation fonctionnelle sera nécessaire afin d'obtenir une meilleure connaissance du développement tissu-spécifique.


CHOIX DU SÉQUENÇAGE ET DE L'ANALYSE

Le choix de la technologie de séquençage et l'approche d'analyse des données sont essentiels à la réussite d'une expérience. Les trois technologies de séquençage mentionnées produisent un énorme volume de données de haute qualité, mais chacune a des applications pratiques spécifiques. Le séquençage Illumina et SOLiD produit des ensembles de données courts mais de grande profondeur. Pour le séquençage Illumina, l'utilisateur peut actuellement sélectionner la longueur des lectures dans la plage de 36 nt à 150 nt qui peuvent être séquencées soit à partir d'une extrémité d'un fragment d'ADN (lectures à une seule extrémité) ou à partir des deux extrémités d'un fragment d'ADN ( lectures appariées). Les lectures plus longues et les lectures appariées sont généralement sélectionnées dans les projets d'assemblage de novo, mais des lectures plus courtes sont parfois choisies pour l'alignement sur un génome de référence. Le score de confiance pour une base donnée dans une séquence diminue à mesure que la longueur de lecture augmente, ce qui peut entraver l'alignement et l'analyse en aval. Les données des lectures de séquençage Illumina sont représentées sous forme de séquence nucléotidique réelle et l'analyse peut procéder directement à l'alignement sur un génome de référence ou à un assemblage de novo.

Dans le système SOLiD, l'utilisateur peut actuellement choisir des longueurs de lecture de 35 nt à 75 nt en format single-end ou paired-end. Le système SOLiD séquence deux bases à la fois (il y a donc 16 combinaisons possibles à interroger), et toute base unique doit être séquencée deux fois pour identifier la vraie séquence à une seule position. On pense que cette méthode améliore l'identification des erreurs de séquençage dans l'analyse post-séquençage. Cependant, pour les chercheurs sans génome de référence, ce système de codage à 2 bases est un inconvénient, car la séquence résultante est codée numériquement et ne sera pas immédiatement reconnaissable par un biologiste. Ce n'est que par une analyse ultérieure que la pertinence biologique d'une lecture de séquençage SOLiD est restaurée. Habituellement, les lectures SOLiD sont alignées dans leur format codé sur 2 bases (appelé format « espace colorimétrique ») sur un génome codant sur 2 bases pour reconvertir la séquence en espace nucléotidique, mais sans génome de référence, cela peut nécessiter des informations informatiques supplémentaires pour avoir un sens du séquençage. La conversion directe des lectures de séquençage est possible mais non recommandée car toutes les bases qui suivent une seule erreur dans l'espace couleur créeront des erreurs dans toutes les bases suivantes d'une lecture. Le lecteur est renvoyé au site Web du fabricant pour une explication plus détaillée du double codage. Si un chercheur dans un système non modèle choisit d'utiliser le système SOLiD, alors le génome d'un parent proche peut être l'option la plus directe pour l'analyse en aval.

La distribution des longueurs de lecture des systèmes Illumina et SOLiD est très uniforme et la plupart des lectures correspondent exactement à la longueur demandée par le chercheur. Dans le séquençage 454 de Roche, les lectures ont une distribution de longueur de séquence plus large et 454 lectures sont également codées dans l'espace nucléotidique normal. La plupart des 454 lectures sont maintenant plus longues que 500 nt, avec un mode autour de 700 nt et une longueur maximale supérieure à 1000 nt. Les longues lectures du séquenceur 454 aboutissent généralement à des assemblages de transcriptome de haute qualité, mais ces ensembles de données ont une profondeur beaucoup plus faible par dollar de séquençage dépensé. L'analyse des données de lecture courte à grande profondeur est fondamentalement différente de l'analyse des lectures longues à faible profondeur, et par conséquent, les ressources de calcul et les approches d'analyse diffèrent considérablement selon le choix de la plate-forme. En partie, ces différences sont enracinées dans la préparation de la bibliothèque.


Méthodes

Matériel végétal

Safran (C. sativus L.) ont été récoltées sur une terre agricole ouverte d'un village situé dans la ville de Pampore du district de Pulwama, Jammu-et-Cachemire, Inde. Différents tissus, y compris le corme, le tépale, la feuille, le stigmate et l'étamine ont été récoltés sur les plantes et immédiatement congelés dans de l'azote liquide et stockés à -80 °C jusqu'à une utilisation ultérieure.

Isolement d'ARN et séquençage du transcriptome

L'ARN total des tissus ci-dessus a été isolé dans trois réplicats biologiques en utilisant le réactif TRI (Sigma Life Science, USA). La quantité et la qualité de l'ARN total ont été déterminées par le spectrophotomètre Nanodrop (Thermo Fisher Scientific) et le Bioanalyzer (Agilent technologies, Singapour). La pureté de l'ARN total a été vérifiée en estimant le rapport d'absorbance à 260/280 et 260/230 et le nombre d'intégrité de l'ARN (RIN). La qualité de l'ARN total isolé des tissus du corme et du stigmate ne répondait pas à la norme minimale pour le séquençage d'Illumina. Par conséquent, nous avons modifié le protocole standard pour obtenir une meilleure qualité d'ARN, qui comprenait le lavage du culot d'ARN avec du NaCl 5 M (2 à 3 fois) avant de le dissoudre dans de l'eau sans RNase. L'ARN total de haute qualité (260/280, 1,8-2,0 260/230 > 2,0 RIN > 7,5) regroupé en quantité égale à partir des trois réplicats biologiques pour chaque échantillon a été utilisé pour le séquençage du transcriptome à l'aide de la plate-forme Illumina pour générer des paires de 100 nt de long. fin de lecture. Pour obtenir des données propres de haute qualité pour de novo l'assemblage, un contrôle de qualité rigoureux a été effectué pour supprimer les lectures de mauvaise qualité et le rognage de l'adaptateur à l'aide de notre kit d'outils NGS QC (v2.3) interne 42 .

De novo assemblage du transcriptome

Des lectures de haute qualité ont été assemblées en contigs à l'aide de divers assembleurs de lecture courte couramment utilisés, tels que Velvet (v1.2.01) 43 , Oases (v0.2.04) 44 , ABySS (v1.2.6) 45 , SOAPdenovo (v1.04) 46 , CLC Genomics Workbench (v4.7.2) et Trinity (v2012-05-18) 47 . L'assemblage du transcriptome a été réalisé en utilisant deux approches différentes, comme décrit précédemment 48 . Dans la première approche (meilleure k-mer), des lectures de haute qualité ont été assemblées à divers k-mer longueur 39-99 en utilisant Velvet, Oases, ABySS et SOAPdenovo, tandis que les logiciels CLC Genomics Workbench et Trinity ont été utilisés avec les paramètres par défaut. Dans la seconde approche (additif k-mer suivi de TGICL), une stratégie en deux étapes a été utilisée pour l'assemblage. Premièrement, les contigs générés pour tous k-mers par assembleur respectif ont été fusionnés et la redondance a été supprimée à l'aide de l'outil CD-HIT. Ensuite, l'ensemble non redondant de contigs a été assemblé à l'aide de la suite TGICL (v2.0) 49 avec une longueur de chevauchement minimale de 40 et une identité maximale de 90. Analyse du contenu GC de C. sativus transcriptome a été réalisé à l'aide d'un script perl interne.

Annotation fonctionnelle

Attribuer la fonction putative à chaque transcription de C. sativus, une recherche de similarité utilisant BLASTX 50 a été effectuée contre Arabidopsis et les protéomes de riz suivis des bases de données NCBI non redondantes et UniRef90 avec un E-valeur seuil de ≤10 -5 pour trouver la meilleure correspondance significative pour chaque transcrit. Les termes GOSlim ont été attribués à chaque C. sativus transcrit sous fonction moléculaire, processus biologique et catégories de composants cellulaires en comparant la séquence avec les protéines d'Arabidopsis. De même, le classement des C. sativus transcriptions dans différentes catégories fonctionnelles a été réalisée en utilisant la base de données KOG. Identification des familles de TF dans C. sativus Le transcriptome a été effectué sur la base d'une recherche de profil de modèle de Markov caché (HMM) (obtenu à partir de la base de données PFAM ou généré à partir des alignements de domaines conservés) en utilisant des critères donnés dans la base de données des facteurs de transcription des plantes (http://plntfdb.bio.uni-potsdam .de/v3.0/) comme décrit précédemment 21 .

Identification des SSR

C. sativus transcriptome a été scanné pour la présence de séquences répétées simples (SSR) en utilisant MISA (MicroSAtellite) aux paramètres par défaut 51. Le nombre minimum d'unités répétées pour le di-nucléotide était de six, alors que pour les tri-, tétra-, penta- et hexa-nucléotides, le nombre minimum d'unités répétées était supérieur à cinq dans les critères de recherche.

Analyse de l'expression génique différentielle

Pour estimer le modèle d'expression de chaque transcrit dans différents échantillons de tissus, des lectures de haute qualité de chaque échantillon ont été cartographiées sur l'assemblage final du transcriptome à l'aide de CLC Genomics Workbench. Un maximum de deux discordances a été autorisé pour le mappage. Le nombre de lectures a été normalisé en calculant le nombre de lectures par kilobase par million (RPKM) pour chaque transcrit dans un tissu individuel. L'analyse différentielle de l'expression génique a été réalisée à l'aide du logiciel DESeq (v1.10.1) 52 sur la base d'une distribution binomiale négative. UNE P-valeur seuil de ≤ 0,05 avec au moins deux fois un changement a été utilisé pour identifier une expression différentielle significative des transcrits. La carte thermique montrant les modèles d'expression spécifiques aux tissus (log2 fold change) pour les transcrits impliqués dans diverses voies ont été générés via TIGR MultiExperiment Viewer (MeV, v4.8).

Analyse PCR en temps réel

Pour l'analyse PCR en temps réel, les amorces spécifiques du gène (tableau S5) ont été conçues à l'aide du logiciel Primer Express (v3.0) (Applied Biosystems, USA). Des PCR en temps réel ont été réalisées dans trois réplicats biologiques indépendants et trois réplicats techniques pour chaque réplicat biologique de chaque échantillon de tissu, comme indiqué précédemment 53 . Ubiquitine a été utilisé comme gène de contrôle interne pour la normalisation.

Disponibilité des données

Les données de séquence générées dans cette étude ont été déposées dans Gene Expression Omnibus sous le numéro d'accès GSE65103. L'assemblage du transcriptome, l'annotation fonctionnelle, les SSR et les données d'expression génique ont été mis à disposition sur la page Web Saffron Transcriptome (http://nipgr.res.in/mjain.html?page=saffron).


Fond

Les transcriptomes peuvent désormais être étudiés par séquençage. Cependant, en l'absence d'un génome de référence, l'assemblage de novo reste une tâche difficile. La principale difficulté vient certainement du fait que les lectures de séquençage sont courtes, et les séquences répétées au sein des transcriptomes pourraient être plus longues que les lectures. Ce problème de lecture courte/répétition longue n'est bien entendu pas spécifique au séquençage du transcriptome. C'est un vieux problème qui existe depuis les premiers algorithmes d'assemblage de génomes. Même si les problèmes que causent les répétitions dans les deux contextes sont similaires, ils ont aussi des caractéristiques propres à chacun. Dans l'assemblage du génome, les répétitions ont tendance à être plus longues et présentes en plus de copies. Dans l'assemblage du transcriptome, les répétitions sont situées dans les gènes et ont tendance à être plus courtes et en moins de copies. Cependant, dans ce dernier cas, la couverture ne peut pas être appliquée pour discriminer les contigs qui correspondent à des répétitions, comme cela peut être le cas en génomique en utilisant par ex. La statistique A de Myers [6, 7], car la couverture d'un gène ne reflète pas seulement son nombre de copies dans le génome, mais aussi et surtout son niveau d'expression. Certains gènes sont fortement exprimés et donc fortement couverts, alors que la plupart des gènes sont faiblement exprimés et donc mal couverts. De telles spécificités compliquent l'application d'une stratégie de résolution de répétition génomique au contexte transcriptomique.

Initialement, on pensait que les répétitions ne seraient pas un problème majeur dans l'ARN-seq, car elles se trouvent principalement dans les introns et les régions intergéniques. Cependant, la vérité est que de nombreuses régions que l'on pense être intergéniques sont transcrites [8] et les introns ne sont pas toujours déjà épissés lorsque l'ARNm est collecté pour être séquencé [9]. Les répétitions, en particulier les éléments transposables, sont donc très présentes dans les échantillons réels et causent des problèmes majeurs dans l'assemblage du transcriptome, si elles ne sont pas traitées correctement.

La plupart, sinon tous les assembleurs de transcriptome à lecture courte actuels sont basés sur des graphes de Bruijn. Parmi les plus connues figurent Oases [3], Trinity [2], et dans une moindre mesure Trans-Abyss [10] et IDBA-tran [11]. Le point commun à tous est l'absence d'un modèle clair et explicite pour les répétitions dans les données RNA-seq. Les heuristiques sont ainsi utilisées pour essayer de gérer efficacement les répétitions. Par exemple, dans les oasis, on pense que les sommets courts correspondent à des répétitions et ne sont donc pas utilisés pour assembler des gènes. Ils sont ajoutés dans une deuxième étape, ce qui, espérons-le, empêche les gènes partageant des répétitions de s'assembler. Dans Trinity, il n'y a aucune tentative de traiter les répétitions en les modélisant explicitement. Le premier module de Trinity, Inchworm, tentera d'assembler le contig le plus couvert qui, espérons-le, correspond au transcrit alternatif le plus abondant. Ensuite, des exons alternatifs sont collés à ce transcrit majeur pour former un graphe d'épissage. La dernière étape consiste à énumérer toutes les transcriptions alternatives. Si des répétitions sont présentes, leur couverture élevée peut être interprétée comme un lien fortement exprimé entre deux transcrits non apparentés. Dans l'ensemble, les transcrits assemblés peuvent être chimériques ou épissés en de nombreux sous-transcriptions.

Dans la méthode que nous avions développée précédemment, KisSplice, qui est un assembleur local de transcriptome [12], les répétitions sont moins problématiques puisque le but n'est pas d'assembler des transcrits complets. KisSplice vise plutôt à trouver des variantes dans les transcriptomes (SNP, indels et épissages alternatifs). Cependant, comme nous l'avons signalé dans [12], KisSplice n'était pas capable de traiter de grandes portions d'un graphe de de Bruijn contenant des sous-graphes associés à des séquences très répétées, par ex. éléments transposables, les Composants Biconnectés complexes.

Ici, nous essayons d'atteindre trois objectifs : (1) donner une formalisation claire de la notion de répétitions avec un nombre de copies élevé dans les données RNA-seq, (2) l'appliquer à l'assemblage local du transcriptome en donnant un moyen pratique d'énumérer les bulles qui sont perdues à cause de telles répétitions, et (3) l'appliquent à l'assemblage global du transcriptome en montrant que la topologie du sous-graphe autour d'un transcrit peut donner quelques indications sur son niveau de confiance. Rappelons que nous sommes dans un contexte de novo, nous supposons donc que ni un génome/transcriptome de référence ni une base de données de répétitions connues, par ex. RepBase [13], sont disponibles.

Premièrement, nous introduisons formellement un modèle pour représenter les répétitions à grand nombre de copies et exploitons ses propriétés pour déduire que les sous-graphes associés aux répétitions dans un graphe de Bruijn contiennent peu d'arcs compressibles. Cependant, nous montrons que le problème d'identification, dans un graphe de de Bruijn, d'un sous-graphe correspondant à des répétitions selon une telle caractérisation est NP-complet. Il est donc peu probable qu'il existe un algorithme en temps polynomial.

Deuxièmement, nous montrons que dans le cas spécifique d'un assemblage local d'événements d'épissage alternatif (AS), en utilisant une stratégie basée sur la caractérisation de l'arc compressible, nous pouvons implicitement éviter de tels sous-graphes. Plus précisément, il est possible de trouver les structures (c. Bien qu'il y ait eu de grands efforts dans la littérature pour résoudre les répétitions, il n'y a eu presque aucune exploration sur la façon de les éviter. Cela s'explique par le fait que la plupart des efforts d'assemblage se concentrent sur l'assemblage complet du génome et du transcriptome, dans lequel éviter les répétitions n'est pas une option, et les performances d'un assembleur peuvent être réduites à la manière dont il résout les répétitions. Cependant, dans notre cas, l'évitement des répétitions peut être une technique efficace. En effet, ce fait a été confirmé par nos expériences, où en utilisant des données RNA-seq simulées humaines, nous montrons que le nouvel algorithme améliore considérablement la sensibilité de KisSplice, tout en améliorant également sa précision. Nous avons ensuite comparé notre algorithme à deux des meilleurs assembleurs de transcriptome, à savoir Trinity [2] et Oases [3], dans la tâche spécifique d'appeler des événements AS, et nous montrons que notre algorithme est plus sensible que les deux outils, tout en étant plus précis. De plus, nos résultats montrent que l'avantage d'utiliser le nouvel algorithme proposé dans ce travail est plus évident lorsque les données d'entrée contiennent une teneur élevée en pré-ARNm ou que les événements AS d'intérêt proviennent de gènes hautement exprimés. De plus, nous donnons une indication de l'utilité de notre méthode sur des données réelles.

Troisièmement, nous montrons que la méthode décrite peut également être appliquée dans le contexte de l'assemblage de transcriptome pleine longueur. Nous introduisons une mesure basée sur le modèle proposé pour identifier les transcrits à faible confiance, qui sont ceux qui traversent des régions complexes dans le graphe de de Bruijn. Au sein de ces parties complexes du graphe généré par les répétitions, tout assembleur devra choisir le(s) chemin(s) « bon » parmi les nombreux présents. Ce choix n'est pas simple et peut conduire à des solutions incorrectes (par exemple des transcrits chimériques ou tronqués). Il est donc important de pouvoir identifier les transcrits provenant de régions aussi complexes afin de savoir que la solution présentée n'est pas la seule, et de plus peut ne pas être la bonne. Nous avons comparé notre mesure à deux méthodes de pointe pour l'évaluation de novo du transcriptome, à savoir Rsem-Eval [4] et TransRate [5], pour la tâche spécifique d'identifier les transcrits chimériques dans les ensembles de données réels et simulés. Nous montrons que notre mesure fournit de bons résultats malgré le fait qu'elle utilise uniquement la topologie du graphe, et non la couverture, ni les informations de lecture. Les résultats obtenus suggèrent ainsi que l'exploration de la topologie du sous-graphe autour d'un transcrit, information actuellement ignorée par les méthodes d'évaluation du transcriptome, peut être utile pour en déduire certaines propriétés du transcrit, telles que le niveau de confiance, la qualité, la dureté d'assemblage, etc. Par conséquent, notre mesure peut améliorer les méthodes de pointe pour l'évaluation de novo du transcriptome, car elle est capable de capturer les erreurs d'assemblage manquées par ces outils.

Préliminaires

Soit (Sigma) un alphabet de taille fixe (sigma) . Ici, nous supposons toujours (Sigma =) . Étant donné une séquence (chaîne) (s in Sigma ^*) , soit |s| indiquer sa longueur, s[je] les jel'élément de s, et s[je, j] la sous-chaîne (s[i] s[i+1] ldots s[j]) pour tout (1 le i<j le |s|) .

UNE k-mer est une séquence (s in Sigma ^k) . Étant donné un entier k et un ensemble S de séquences chacune de longueur (n ge k) , on définit envergure(S, k) comme l'ensemble de tous les k-mers qui apparaissent comme une sous-chaîne dans S.

Définition 1

Étant donné un ensemble de séquences (lectures) (R subseteq Sigma ^*) et un entier k, on définit le graphe de Bruijn orienté (G_k(R)=(V,A)) où (V=span(R,k)) et ((u,v) in A) si et seulement si (u[2,k]=v[1,k-1]) .

Etant donné un graphe orienté (G = (V,A)) et un sommet (v in V) , on note son hors-quartier (resp. dans le quartier) par (N^+(v)=< u in V mid (v,u) in A >) (resp. (N^-(v)=< u in V mid (u,v) in A >) ), et son hors-degré (resp. en-diplôme par (d^+(v)=|N^+(v)|) ( (d^-(v)=|N^-(v)|) ). Un simple) chemin (pi = s leadsto t) dans g est une suite de sommets distincts (s = v_0, ldots , v_l = t) telle que, pour chaque (0 le i < l) , ((v_i, v_)) est un arc de g. Si le graphe est pondéré, c'est-à-dire qu'il existe une fonction (w : A ightarrow Q_) associant un poids à chaque arc du graphe, alors le longueur d'un chemin (pi) est la somme des poids des arcs traversés, et est noté (|pi |) .

Un arc ((u,v) in A) est appelé compressible si (d^+(u)=1) et (d^-(v)=1) . L'intuition derrière cette définition vient du fait que chaque chemin passant par vous devrait également passer par v. Il devrait donc être possible de « comprimer » ou de contracter cet arc sans perdre aucune information. Notons que le graphe de Bruijn compressé [2, 3] couramment utilisé par les assembleurs transcriptomiques est obtenu à partir d'un graphe de Bruijn en remplaçant, pour chaque arc compressible (vous, v), les sommets vous, v par un nouveau sommet X, où (N^-(x) = N^-(u)) , (N^+(x) = N^+(v)) et le label est la concaténation des k-mer de vous et le k-mer de v sans la partie chevauchante (voir Fig. 1).

Exemple d'arc compressible dans un graphe de de Bruijn. une L'arc (CTG, ATG) est le seul arc compressible dans le graphe de Bruijn donné ( (k=3) ). b Le graphe de Bruijn compressé correspondant


Conclusion

En utilisant une approche protéomique et transcriptomique intégrée, nous avons révélé des différences dans le protéome hépatique des RMN à longue durée de vie par rapport aux GP à durée de vie plus courte. Nous avons confirmé une manière préférentielle d'utiliser les acides gras pour alimenter la respiration dans les RMN, reflétant une composition distincte de leurs mitochondries. De plus, nous avons identifié une signature progressive du vieillissement qui se manifeste dans le foie des RMN au niveau moléculaire. Il est intéressant de noter que des groupes de protéines fonctionnellement apparentés, y compris des enzymes des voies de détoxification, ont été affectés de la même manière par le vieillissement dans les échantillons de RMN et de foie humain. Ceci souligne un lien direct entre les processus de vieillissement de ces deux espèces. Il reste à démontrer si les voies affectées par le vieillissement influencent l'état de santé des anciennes RMN et limitent leur durée de vie, comme nous l'avons montré chez le nématode C. elegans.


RÉSULTATS

Précision de la prédiction des gènes dans les transcrits de référence

Nous avons utilisé GeneMarkS-T, Prodigal, TransDecoder et ESTscan pour prédire les gènes codant pour des protéines dans des transcrits « complets » et « partiels » de A. thaliana, D. melanogaster, M. musculus et S. pombe (voir la section « Matériaux et méthodes »). Le nombre de gènes prédit dans un ensemble de transcrits dépend de la longueur minimale de gène sélectionnée (mgl). nous avons changé mgl comme paramètre seuil de 90 à 480 pb (par pas de 30 pb). Pour chaque ensemble de prédictions, nous avons calculé Sn et Sp sur la base de l'annotation de transcription et tracé la dépendance de Sn sur 1 − Sp (Figures 2 et 3). Dans ces graphiques, qui ressemblent aux courbes des caractéristiques de fonctionnement du récepteur (ROC), les points en haut à droite ont été obtenus pour mgl égal à 90 pb. Nous ne montrons pas de graphiques pour ESTscan car nous n'avons pas été en mesure d'atteindre des performances suffisamment élevées (c'est-à-dire que pour la souris, nous avions Sn = 0,53 et Sp = 0,54). Nous pensons que l'auto-formation améliorerait les performances d'ESTscan. En l'absence d'une telle option, nous avons été obligés de sélectionner l'un des modèles prédéfinis disponibles, par ex. le modèle humain pour l'analyse des transcrits de souris.

Graphiques de la sensibilité de prédiction génique (Sn) en fonction de la spécificité de prédiction génique (1 − Sp) pour TransDecoder, Prodigal et GeneMarkS-T déterminés sur des ensembles de tests de transcrits de référence « complets » de A. thaliana, D. melanogaster, M. musculus et S. pombe. Nous avons appliqué les trois outils à la fois en mode aveugle et en mode informé (S). Pour construire les courbes, nous avons généré des ensembles de gènes prédits avec une longueur minimale contrôlée par le mgl seuil (voir texte). Comme le mgl les valeurs ont augmenté de 90 à 480 pb (avec un pas de 30 pb) les valeurs de Sn ont diminué.

Graphiques de la sensibilité de prédiction génique (Sn) en fonction de la spécificité de prédiction génique (1 − Sp) pour TransDecoder, Prodigal et GeneMarkS-T déterminés sur des ensembles de tests de transcrits de référence « complets » de A. thaliana, D. melanogaster, M. musculus et S. pombe. Nous avons appliqué les trois outils à la fois en mode aveugle et en mode informé (S). Pour construire les courbes, nous avons généré des ensembles de gènes prédits avec une longueur minimale contrôlée par le mgl seuil (voir texte). Comme le mgl les valeurs ont augmenté de 90 à 480 pb (avec un pas de 30 pb) les valeurs de Sn ont diminué.

Identique à la figure 2 pour les tests sur des transcriptions de référence « partielles » simulées de A. thaliana, D. melanogaster, M. musculus et S. pombe. Les transcriptions « partielles » ont été faites en coupant les séquences à la fois sur les extrémités 5′ et 3′ des transcriptions « complètes » (voir le texte pour le rationnel de cette méthode). Les trois outils ont été utilisés à la fois en mode aveugle et en mode informé (S).

Identique à la figure 2 pour les tests sur des transcriptions de référence « partielles » simulées de A. thaliana, D. melanogaster, M. musculus et S. pombe. Les transcriptions « partielles » ont été faites en coupant les séquences à la fois sur les extrémités 5′ et 3′ des transcriptions « complètes » (voir le texte pour le rationnel de cette méthode). Les trois outils ont été utilisés à la fois en mode aveugle et en mode informé (S).

Pour les transcrits « complets », les versions brin aveugle et spécifique au brin de GeneMarkS-T ont démontré des performances significativement meilleures que les autres outils (Figure 2). Dans des expériences avec des transcrits «partiels» (Figure 3), Prodigal et TransDecoder se sont rapprochés en termes de performances de GeneMarkS-T. Le meilleur (Sn + Sp)/2 que nous ayons vu pour GeneMarkS-T, Prodigal et TransDecoder lorsque le mgl les valeurs étaient respectivement de 150, 210 et 270 pb. L'ajout d'informations sur le brin d'ARN et donc l'utilisation des versions (S) des trois outils de recherche de gènes, a augmenté les valeurs de Sp (figures 2 et 3).

Variation significative de la teneur en G + C dans M. musculus et D. melanogaster Les transcrits (de 0,31 à 0,76 chez la souris et de 0,27 à 0,63 chez la mouche) ont été immédiatement identifiés par GeneMarkS-T qui a regroupé les transcrits dans trois bacs de contenu G + C avec des bordures automatiquement définies (tableau S1). L'auto-formation a été effectuée séparément pour les transcriptions dans chacun des trois groupes. Dans l'étape de prédiction, les paramètres d'algorithme utilisés pour un transcrit donné ont été choisis par rapport au contenu du transcrit G + C. Cette approche a produit de meilleures valeurs de Sn qu'en l'absence de clustering (tableau S1).

Nous avons étudié comment la précision des prédictions dépend du volume de transcriptions en formation. Pour ces expériences, nous avons échantillonné au hasard plusieurs ensembles de transcriptions avec le même volume. Si le volume était supérieur à 600 kb, GeneMarkS-T et Prodigal atteignaient un plateau avec des performances stables et une valeur (Sn + Sp)/2 proche de 96% pour GeneMarkS-T et 94% pour Prodigal (Figure 4). La précision de TransDecoder avait un modèle de changement similaire avec le plateau à 91% atteint au volume de 1 Mb. Une diminution à 100 ko a produit des performances inférieures mais toujours décentes : 90 % pour GeneMarkS-T et Prodigal, et 80 % pour TransDecoder. Le volume minimum de séquence requis pour Prodigal était de 20 kb tandis que la limite GeneMarkS-T était encore plus faible. En dessous de 50 kb, GeneMarkS-T passe automatiquement à l'utilisation de modèles heuristiques de régions codant pour des protéines dont les paramètres pourraient être déterminés pour un fragment de séquence aussi court que 400 pb (15).

Dépendance de (Sn + Sp)/2 des trois outils de prédiction de gènes sur la taille de l'ensemble d'apprentissage de D. melanogaster transcriptions (l'axe X montre la taille totale de l'ensemble, l'échelle logarithmique). Sets of transcripts of the same size were sampled randomly 50 times from the whole set of reference transcripts. Les mgl value that achieved best overall (Sn + Sp)/2 was tool specific (150 bp for GeneMarkS-T, 210 bp for Prodigal and 270 bp for TransDecoder).

Dependence of (Sn + Sp)/2 of the three gene prediction tools on the size of training set of D. melanogaster transcripts (X axis shows the total set size, log scale). Sets of transcripts of the same size were sampled randomly 50 times from the whole set of reference transcripts. Les mgl value that achieved best overall (Sn + Sp)/2 was tool specific (150 bp for GeneMarkS-T, 210 bp for Prodigal and 270 bp for TransDecoder).

In some transcripts GeneMarkS-T predicted several coding regions (with mgl 300 bp). We observed such outcomes in 2.5% of A. thaliana transcripts, 9.4% of S. pombe, 6.0% of D. melanogaster and 20.4% of M. musculus. In the supposed absence of operons such outcomes are possible for three reasons. First, additional predictions could have no connection to carrying genetic code, i.e. pure false positives. Second, a transcript could come from a locus where splicing mechanism generates alternative isoforms. For instance, protein-coding exons related to one isoform could appear outside the protein coding region related to another isoform (e.g. Figure 5A). Third, a transcript could overlap adjacent genes located in the complementary strand. Particulièrement, S. pombe, a species not known for ubiquitous alternative splicing, has short intergenic regions and long UTRs that may overlap adjacent genes (e.g. Figure 5B). Not surprisingly, for S. pombe we observed a significant gain of accuracy after switching to strand-specific versions of the three gene finders (Figures 2 and 3).

Diagrams of two typical events when more than one coding region is predicted in a transcript. We show pre-spliced transcripts: genomic sequences are shown as grey bars exons defined by annotation are shown as wider bars (green colour—UTR, dark green—CDS) predicted protein-coding exons are shown as red bars. (UNE) Two transcripts are originated from the same location of D. melanogaster genome (NM_001275246.1 and NM_206418.3). The FP prediction (the downstream gene in complementary strand) is a part of the coding region of alternative isoform of CapaR gene. (B) The 5′ UTR of S. pombe transcript NM_001020436.2 overlaps with another transcript NM_001020437.2 originated from complementary strand. This transcript topology leads to two predictions in transcript NM_001020436.2: one in the direct strand (FP) as well as one in the complementary strand (TP). The figures were made with the NCBI RefSeq sequence viewer.

Diagrams of two typical events when more than one coding region is predicted in a transcript. We show pre-spliced transcripts: genomic sequences are shown as grey bars exons defined by annotation are shown as wider bars (green colour—UTR, dark green—CDS) predicted protein-coding exons are shown as red bars. (UNE) Two transcripts are originated from the same location of D. melanogaster genome (NM_001275246.1 and NM_206418.3). The FP prediction (the downstream gene in complementary strand) is a part of the coding region of alternative isoform of CapaR gene. (B) The 5′ UTR of S. pombe transcript NM_001020436.2 overlaps with another transcript NM_001020437.2 originated from complementary strand. This transcript topology leads to two predictions in transcript NM_001020436.2: one in the direct strand (FP) as well as one in the complementary strand (TP). The figures were made with the NCBI RefSeq sequence viewer.

If multiple predictions were generated in a transcript GeneMarkS-T selected the one with the maximum log-odd score. This approach produced 93% success rate in selecting the ‘true’ coding region for A. thaliana, 74% for D. melanogaster, 98% for M. musculus and 62% for S. pombe.

Prediction of translation initiation site

To assess the accuracy of TIS prediction by GeneMarkS-T, Prodigal and TransDecoder we used 1392 reference transcripts of M. musculus (with annotated coding regions longer than 300 bp). The TIS annotation in these transcripts was validated by Ribo-seq experiments (see ‘Materials and Methods’ section). GeneMarkS-T was run in three modes: (i) with default settings (ii) with search for the Kozak motif switched off and iii/ with mandatory prediction of complete CDS.

GeneMarkS-T with default settings correctly predicted 68.5% starts in genes where the reading frame was correctly predicted (and, therefore, the 3′ end of the gene). This was higher accuracy in comparison with the two other tools (Table 2). All three tools revealed a tendency to extend the 5′ end of the coding region beyond the 5′ end of the transcript. Notably, TransDecoder adopts the ‘longest-ORF’ rule and selects the 5′-most AUG (with respect to the in-frame stop codon) as the translation initiation site. In comparison, GeneMarkS-T had the largest fraction of TIS predictions located downstream from the 5′-most AUGs. Use of the Kozak motif was responsible for improving Sn of GeneMarkS-T by about 10% (Table 2). Prohibiting predictions of incomplete coding regions would boost the TIS identification accuracy of GeneMarkS-T to 95.0%, however, use of this option is limited to transcripts that are known to be 5′ end complete.

Numbers of protein-coding regions predicted correctly (TP) and incorrectly (FP) by GeneMarkS-T, Prodigal and TransDecoder in D. melanogaster ‘concordant’ transcripts (selected as described in text)

Transcripts built by . No. of transcripts . GeneMarkS-T . Prodigal . TransDecoder .
. . TP . FP . TP . FP . TP . FP .
Cufflinks 7222 7162607098 232 7046 432
Auguste 9444 9423219383 246 9332 480
Exonerate 6971 6953186940 190 6915 454
Velvet 7344 71461987096 312 7030 429
Oases 13 869 13 76910013 659 347 13 598 582
Transcripts built by . No. of transcripts . GeneMarkS-T . Prodigal . TransDecoder .
. . TP . FP . TP . FP . TP . FP .
Cufflinks 7222 7162607098 232 7046 432
Auguste 9444 9423219383 246 9332 480
Exonerate 6971 6953186940 190 6915 454
Velvet 7344 71461987096 312 7030 429
Oases 13 869 13 76910013 659 347 13 598 582

Predictions shorter than the tool-specific mgl (150 bp for GeneMarkS-T, 210 bp for Prodigal and 270 bp for TransDecoder) were filtered out. Bold font highlights best results in a particular row (the largest TP and the smallest FP).

Transcripts built by . No. of transcripts . GeneMarkS-T . Prodigal . TransDecoder .
. . TP . FP . TP . FP . TP . FP .
Cufflinks 7222 7162607098 232 7046 432
Auguste 9444 9423219383 246 9332 480
Exonerate 6971 6953186940 190 6915 454
Velvet 7344 71461987096 312 7030 429
Oases 13 869 13 76910013 659 347 13 598 582
Transcripts built by . No. of transcripts . GeneMarkS-T . Prodigal . TransDecoder .
. . TP . FP . TP . FP . TP . FP .
Cufflinks 7222 7162607098 232 7046 432
Auguste 9444 9423219383 246 9332 480
Exonerate 6971 6953186940 190 6915 454
Velvet 7344 71461987096 312 7030 429
Oases 13 869 13 76910013 659 347 13 598 582

Predictions shorter than the tool-specific mgl (150 bp for GeneMarkS-T, 210 bp for Prodigal and 270 bp for TransDecoder) were filtered out. Bold font highlights best results in a particular row (the largest TP and the smallest FP).

Several ribosome profiling studies ( 12, 23–24) raised concerns about the frequent presence of alternative TIS's located both upstream and downstream of annotated TIS's confirmed by Ribo-seq experiments. However, a recent publication ( 25) indicated that reports of alternative TIS in many cases are likely to be artefacts therefore, the confidence in the Ribo-seq experimental validation of annotated TIS's remains high.

Gene prediction with heuristic models (case for meta-transcriptomics)

To model gene prediction in a metatranscriptome we used the same set of mouse transcripts G + C content of individual transcripts in this set ranged from 27 to 63%. To run GeneMarkS-T on a given transcript we used parameters derived as functions of a single variable, the transcript G + C content. We did not continue the training, assuming that the given transcript is the only sequence from an unknown genome. This assumption is relevant for a typical metatranscriptome. The method of inference of these functions was described earlier for short metagenomics sequences ( 7, 15). We used the functions that reflect dependence of oligonucleotide composition of protein coding regions on G + C content of the sequence the functions were derived for a set of complete prokaryotic genomes ( 15). The results are surprisingly good (last row in Table 2) with correct prediction of 82.4% of genes (1147 out of 1193) also 54.9% of starts were correctly predicted in comparison with 68.6% correct starts predicted with full training of the model.

Model training and gene predictions for transcripts reconstructed from RNA-Seq

A comprehensive assessment of the accuracy of transcript reconstruction from RNA-Seq reads was conducted in the RGASP competition ( 3). We used in this study transcripts reconstructed in ( 3) by Cufflinks, Augustus, Exonerate, Velvet and Oases ( 18–22). It was shown that assembled transcripts frequently contain errors and only a subset of all transcripts could be fully recovered ( 3). Observed average lengths of assembled transcripts were shorter than that of reference transcripts, particularly the average lengths of the de novo assemblies made by Oases and Velvet (Supplementary Figure S1A). Would the errors present in transcript assemblies affect self-training of GeneMarkS-T? To address this question we trained GeneMarkS-T on five sets of D. melanogaster transcripts assembled by the five tools mentioned above. The trained models were used in GeneMarkS-T to predict genes in reference transcripts of D. melanogaster. We observed almost no difference between any of the five graphs of Sn versus 1 − Sp for gene prediction with models trained on D. melanogaster assembled transcripts and the graph depicting Sn versus 1 − Sp for gene prediction with the D. melanogaster model trained on reference transcripts (Figure 6). Thus, GeneMarkS-T training was shown to be robust with respect to use of assembled transcripts instead of ‘ideal’ reference transcripts.

Plots of gene prediction accuracy in D. melanogaster reference transcripts built for GeneMarkS-T trained on sets of different types. The models were trained either on the set of D. melanogaster reference transcripts or on the sets of transcripts assembled by the five transcript assembly tools. Predictions made in reference transcripts were compared with annotation.

Plots of gene prediction accuracy in D. melanogaster reference transcripts built for GeneMarkS-T trained on sets of different types. The models were trained either on the set of D. melanogaster reference transcripts or on the sets of transcripts assembled by the five transcript assembly tools. Predictions made in reference transcripts were compared with annotation.

To assess performance of gene prediction methods in assembled transcripts we used the same five sets of assembled D. melanogaster transcriptions. First, we mapped the assembled transcripts to the corresponding reference transcripts ( 3) to detect and evaluate the differences. We used the results to divide the set of assembled transcripts into three groups: ‘concordant’, ‘conflicting’ and ‘not-aligned’ (see ‘Materials and Methods’ section and Supplementary Figure S2). Many assembled D. melanogaster transcripts fell into ‘conflicting’ category (from 17 to 47%, depending on the tool, see Supplementary Figure S3, ‘A’ bars) Cufflinks, Exonerate and Oases produced larger numbers of ‘conflicting’ transcripts than Augustus and Velvet. Multiple protein-coding regions were predicted more frequently in the ‘conflicting’ transcripts than in the ‘concordant’ transcripts (Supplementary Figure S4). Note, that for GeneMarkS-T events of prediction of multiple coding regions were registered prior to selecting ‘reported’ predictions with highest log-odd score. We have illustrated the distribution of events (multiple, single, none predictions) for GeneMarkS-T (Supplementary Figure S4). The distributions of the same events for the two other gene prediction tools show similar patterns (Table S2). Thus, all the tools predict single coding regions in ‘concordant’ assemblies with higher frequencies than in ‘conflicting’ ones.

To make unambiguous comparison of accuracy of gene prediction in ‘concordant’ transcripts we had to select the sets where gene finders make single gene predictions. As such surrogate sets we chose sets of ‘concordant’ assemblies where GeneMarkS-T predicted single protein-coding regions. Annotation of protein coding regions in these assembled transcripts was accomplished by transfer of the reference transcript annotation. In all the five test sets, GeneMarkS-T generated the largest number of TPs and the fewest number of FPs (Table 3).

In the sets of assembled transcripts where GeneMarkS-T predicted multiple coding regions we have observed high fractions of ‘conflicting’ transcripts (e.g. 90%, for the set of Cufflinks assembled transcripts). Thus, predicting multiple coding regions was an indicator of a higher chance for the transcript to be in the ‘conflicting’ category and to carry some discrepancies in the transcript assembly. Still, this observation should be taken with a caveat that multiple coding regions could appear in the ‘concordant’ transcript encoding alternative isoforms (as illustrated in Figure 5).

Very short coding regions are rare and are rarely predicted. Therefore, if an assembled transcript (complete or incomplete) is short it is likely that no gene will be predicted. Indeed, we observed that the gene finding tools did not predict genes in many transcripts assembled by the de novo methods Velvet and Oases (Supplementary Figure S3). Notably, many of these transcripts were too short (Supplementary Figure S1A).


RÉSULTATS

Accuracy of gene prediction in reference transcripts

We used GeneMarkS-T, Prodigal, TransDecoder and ESTscan to predict protein-coding genes in ‘complete’ as well as ‘partial’ transcripts of A. thaliana, D. melanogaster, M. musculus et S. pombe (see ‘Materials and Methods’ section). The number of genes predicted in a set of transcripts depends on the selected minimum gene length (mgl). We have changed mgl as a threshold parameter from 90 to 480 bp (with 30 bp steps). For each set of predictions we computed Sn and Sp based on the transcript annotation and plotted the dependence of Sn on 1 − Sp (Figures 2 and 3). In these plots, which look similar to receiver operating characteristic (ROC) curves, the top right points were obtained for mgl equal to 90 bp. We do not show plots for ESTscan as we were not able to achieve high enough performance (i.e. for mouse we had Sn = 0.53 and Sp = 0.54). We believe that self-training would improve ESTscan performance. In the absence of such an option we were forced to select one of the available pre-defined models, e.g. the human model for analysis of mouse transcripts.

Plots of gene prediction sensitivity (Sn) as functions of gene prediction specificity (1 − Sp) for TransDecoder, Prodigal and GeneMarkS-T determined on test sets of ‘complete’ reference transcripts of A. thaliana, D. melanogaster, M. musculus et S. pombe. We applied the three tools in both strand blind and strand informed (S) modes. To build the curves we generated sets of predicted genes with minimal length controlled by the mgl threshold (see text). Comme le mgl values increased from 90 to 480 bp (with 30 bp step) the Sn values decreased.

Plots of gene prediction sensitivity (Sn) as functions of gene prediction specificity (1 − Sp) for TransDecoder, Prodigal and GeneMarkS-T determined on test sets of ‘complete’ reference transcripts of A. thaliana, D. melanogaster, M. musculus et S. pombe. We applied the three tools in both strand blind and strand informed (S) modes. To build the curves we generated sets of predicted genes with minimal length controlled by the mgl threshold (see text). Comme le mgl values increased from 90 to 480 bp (with 30 bp step) the Sn values decreased.

Same as in Figure 2 for the tests on simulated ‘partial’ reference transcripts of A. thaliana, D. melanogaster, M. musculus et S. pombe. The ‘partial’ transcripts were made by trimming sequences on both 5′ and 3′ end of the ‘complete’ transcripts (see text for rational of this method). The three tools were used in both strand blind and strand informed (S) modes.

Same as in Figure 2 for the tests on simulated ‘partial’ reference transcripts of A. thaliana, D. melanogaster, M. musculus et S. pombe. The ‘partial’ transcripts were made by trimming sequences on both 5′ and 3′ end of the ‘complete’ transcripts (see text for rational of this method). The three tools were used in both strand blind and strand informed (S) modes.

For ‘complete’ transcripts, both strand-blind and strand-specific versions of GeneMarkS-T demonstrated significantly better performance than the other tools (Figure 2). In experiments with ‘partial’ transcripts (Figure 3) Prodigal and TransDecoder came closer in performance to GeneMarkS-T. The best (Sn + Sp)/2 we saw for GeneMarkS-T, Prodigal and TransDecoder when the mgl values were 150, 210 and 270 bp, respectively. Adding information on RNA strand and thus use of the (S) versions of the three gene finding tools, increased the Sp values (Figures 2 and 3).

Significant variation in G + C content in M. musculus et D. melanogaster transcripts (from 0.31 to 0.76 in mouse and from 0.27 to 0.63 in fly) was immediately identified by GeneMarkS-T which grouped the transcripts into three G + C content bins with automatically defined borders (Table S1). Self-training was done separately for transcripts in each of the three clusters. In the prediction step, algorithm parameters used for a given transcript were chosen with respect to the transcript G + C content. This approach produced better Sn values than in the absence of clustering (Table S1).

We studied how prediction accuracy depends on the volume of transcripts in training. For these experiments we sampled randomly several sets of transcripts with the same volume. If the volume was larger than 600 kb, GeneMarkS-T and Prodigal reached a plateau with steady performance and (Sn + Sp)/2 value close to 96% for GeneMarkS-T and 94% for Prodigal (Figure 4). Accuracy of TransDecoder had a similar pattern of change with the plateau at 91% reached at the volume of 1 Mb. A decrease to 100 kb produced lower but still decent performance: 90% for GeneMarkS-T and Prodigal, and 80% for TransDecoder. The minimum volume of sequence required for Prodigal was 20 kb while the GeneMarkS-T limit was even lower. Below 50 kb GeneMarkS-T automatically switches to use of heuristic models of protein-coding regions whose parameters could be determined for a sequence fragment as short as 400 bp ( 15).

Dependence of (Sn + Sp)/2 of the three gene prediction tools on the size of training set of D. melanogaster transcripts (X axis shows the total set size, log scale). Sets of transcripts of the same size were sampled randomly 50 times from the whole set of reference transcripts. Les mgl value that achieved best overall (Sn + Sp)/2 was tool specific (150 bp for GeneMarkS-T, 210 bp for Prodigal and 270 bp for TransDecoder).

Dependence of (Sn + Sp)/2 of the three gene prediction tools on the size of training set of D. melanogaster transcripts (X axis shows the total set size, log scale). Sets of transcripts of the same size were sampled randomly 50 times from the whole set of reference transcripts. Les mgl value that achieved best overall (Sn + Sp)/2 was tool specific (150 bp for GeneMarkS-T, 210 bp for Prodigal and 270 bp for TransDecoder).

In some transcripts GeneMarkS-T predicted several coding regions (with mgl 300 bp). We observed such outcomes in 2.5% of A. thaliana transcripts, 9.4% of S. pombe, 6.0% of D. melanogaster and 20.4% of M. musculus. In the supposed absence of operons such outcomes are possible for three reasons. First, additional predictions could have no connection to carrying genetic code, i.e. pure false positives. Second, a transcript could come from a locus where splicing mechanism generates alternative isoforms. For instance, protein-coding exons related to one isoform could appear outside the protein coding region related to another isoform (e.g. Figure 5A). Third, a transcript could overlap adjacent genes located in the complementary strand. Particulièrement, S. pombe, a species not known for ubiquitous alternative splicing, has short intergenic regions and long UTRs that may overlap adjacent genes (e.g. Figure 5B). Not surprisingly, for S. pombe we observed a significant gain of accuracy after switching to strand-specific versions of the three gene finders (Figures 2 and 3).

Diagrams of two typical events when more than one coding region is predicted in a transcript. We show pre-spliced transcripts: genomic sequences are shown as grey bars exons defined by annotation are shown as wider bars (green colour—UTR, dark green—CDS) predicted protein-coding exons are shown as red bars. (UNE) Two transcripts are originated from the same location of D. melanogaster genome (NM_001275246.1 and NM_206418.3). The FP prediction (the downstream gene in complementary strand) is a part of the coding region of alternative isoform of CapaR gene. (B) The 5′ UTR of S. pombe transcript NM_001020436.2 overlaps with another transcript NM_001020437.2 originated from complementary strand. This transcript topology leads to two predictions in transcript NM_001020436.2: one in the direct strand (FP) as well as one in the complementary strand (TP). The figures were made with the NCBI RefSeq sequence viewer.

Diagrams of two typical events when more than one coding region is predicted in a transcript. We show pre-spliced transcripts: genomic sequences are shown as grey bars exons defined by annotation are shown as wider bars (green colour—UTR, dark green—CDS) predicted protein-coding exons are shown as red bars. (UNE) Two transcripts are originated from the same location of D. melanogaster genome (NM_001275246.1 and NM_206418.3). The FP prediction (the downstream gene in complementary strand) is a part of the coding region of alternative isoform of CapaR gene. (B) The 5′ UTR of S. pombe transcript NM_001020436.2 overlaps with another transcript NM_001020437.2 originated from complementary strand. This transcript topology leads to two predictions in transcript NM_001020436.2: one in the direct strand (FP) as well as one in the complementary strand (TP). The figures were made with the NCBI RefSeq sequence viewer.

If multiple predictions were generated in a transcript GeneMarkS-T selected the one with the maximum log-odd score. This approach produced 93% success rate in selecting the ‘true’ coding region for A. thaliana, 74% for D. melanogaster, 98% for M. musculus and 62% for S. pombe.

Prediction of translation initiation site

To assess the accuracy of TIS prediction by GeneMarkS-T, Prodigal and TransDecoder we used 1392 reference transcripts of M. musculus (with annotated coding regions longer than 300 bp). The TIS annotation in these transcripts was validated by Ribo-seq experiments (see ‘Materials and Methods’ section). GeneMarkS-T was run in three modes: (i) with default settings (ii) with search for the Kozak motif switched off and iii/ with mandatory prediction of complete CDS.

GeneMarkS-T with default settings correctly predicted 68.5% starts in genes where the reading frame was correctly predicted (and, therefore, the 3′ end of the gene). This was higher accuracy in comparison with the two other tools (Table 2). All three tools revealed a tendency to extend the 5′ end of the coding region beyond the 5′ end of the transcript. Notably, TransDecoder adopts the ‘longest-ORF’ rule and selects the 5′-most AUG (with respect to the in-frame stop codon) as the translation initiation site. In comparison, GeneMarkS-T had the largest fraction of TIS predictions located downstream from the 5′-most AUGs. Use of the Kozak motif was responsible for improving Sn of GeneMarkS-T by about 10% (Table 2). Prohibiting predictions of incomplete coding regions would boost the TIS identification accuracy of GeneMarkS-T to 95.0%, however, use of this option is limited to transcripts that are known to be 5′ end complete.

Numbers of protein-coding regions predicted correctly (TP) and incorrectly (FP) by GeneMarkS-T, Prodigal and TransDecoder in D. melanogaster ‘concordant’ transcripts (selected as described in text)

Transcripts built by . No. of transcripts . GeneMarkS-T . Prodigal . TransDecoder .
. . TP . FP . TP . FP . TP . FP .
Cufflinks 7222 7162607098 232 7046 432
Auguste 9444 9423219383 246 9332 480
Exonerate 6971 6953186940 190 6915 454
Velvet 7344 71461987096 312 7030 429
Oases 13 869 13 76910013 659 347 13 598 582
Transcripts built by . No. of transcripts . GeneMarkS-T . Prodigal . TransDecoder .
. . TP . FP . TP . FP . TP . FP .
Cufflinks 7222 7162607098 232 7046 432
Auguste 9444 9423219383 246 9332 480
Exonerate 6971 6953186940 190 6915 454
Velvet 7344 71461987096 312 7030 429
Oases 13 869 13 76910013 659 347 13 598 582

Predictions shorter than the tool-specific mgl (150 bp for GeneMarkS-T, 210 bp for Prodigal and 270 bp for TransDecoder) were filtered out. Bold font highlights best results in a particular row (the largest TP and the smallest FP).

Transcripts built by . No. of transcripts . GeneMarkS-T . Prodigal . TransDecoder .
. . TP . FP . TP . FP . TP . FP .
Cufflinks 7222 7162607098 232 7046 432
Auguste 9444 9423219383 246 9332 480
Exonerate 6971 6953186940 190 6915 454
Velvet 7344 71461987096 312 7030 429
Oases 13 869 13 76910013 659 347 13 598 582
Transcripts built by . No. of transcripts . GeneMarkS-T . Prodigal . TransDecoder .
. . TP . FP . TP . FP . TP . FP .
Cufflinks 7222 7162607098 232 7046 432
Auguste 9444 9423219383 246 9332 480
Exonerate 6971 6953186940 190 6915 454
Velvet 7344 71461987096 312 7030 429
Oases 13 869 13 76910013 659 347 13 598 582

Predictions shorter than the tool-specific mgl (150 bp for GeneMarkS-T, 210 bp for Prodigal and 270 bp for TransDecoder) were filtered out. Bold font highlights best results in a particular row (the largest TP and the smallest FP).

Several ribosome profiling studies ( 12, 23–24) raised concerns about the frequent presence of alternative TIS's located both upstream and downstream of annotated TIS's confirmed by Ribo-seq experiments. However, a recent publication ( 25) indicated that reports of alternative TIS in many cases are likely to be artefacts therefore, the confidence in the Ribo-seq experimental validation of annotated TIS's remains high.

Gene prediction with heuristic models (case for meta-transcriptomics)

To model gene prediction in a metatranscriptome we used the same set of mouse transcripts G + C content of individual transcripts in this set ranged from 27 to 63%. To run GeneMarkS-T on a given transcript we used parameters derived as functions of a single variable, the transcript G + C content. We did not continue the training, assuming that the given transcript is the only sequence from an unknown genome. This assumption is relevant for a typical metatranscriptome. The method of inference of these functions was described earlier for short metagenomics sequences ( 7, 15). We used the functions that reflect dependence of oligonucleotide composition of protein coding regions on G + C content of the sequence the functions were derived for a set of complete prokaryotic genomes ( 15). The results are surprisingly good (last row in Table 2) with correct prediction of 82.4% of genes (1147 out of 1193) also 54.9% of starts were correctly predicted in comparison with 68.6% correct starts predicted with full training of the model.

Model training and gene predictions for transcripts reconstructed from RNA-Seq

A comprehensive assessment of the accuracy of transcript reconstruction from RNA-Seq reads was conducted in the RGASP competition ( 3). We used in this study transcripts reconstructed in ( 3) by Cufflinks, Augustus, Exonerate, Velvet and Oases ( 18–22). It was shown that assembled transcripts frequently contain errors and only a subset of all transcripts could be fully recovered ( 3). Observed average lengths of assembled transcripts were shorter than that of reference transcripts, particularly the average lengths of the de novo assemblies made by Oases and Velvet (Supplementary Figure S1A). Would the errors present in transcript assemblies affect self-training of GeneMarkS-T? To address this question we trained GeneMarkS-T on five sets of D. melanogaster transcripts assembled by the five tools mentioned above. The trained models were used in GeneMarkS-T to predict genes in reference transcripts of D. melanogaster. We observed almost no difference between any of the five graphs of Sn versus 1 − Sp for gene prediction with models trained on D. melanogaster assembled transcripts and the graph depicting Sn versus 1 − Sp for gene prediction with the D. melanogaster model trained on reference transcripts (Figure 6). Thus, GeneMarkS-T training was shown to be robust with respect to use of assembled transcripts instead of ‘ideal’ reference transcripts.

Plots of gene prediction accuracy in D. melanogaster reference transcripts built for GeneMarkS-T trained on sets of different types. The models were trained either on the set of D. melanogaster reference transcripts or on the sets of transcripts assembled by the five transcript assembly tools. Predictions made in reference transcripts were compared with annotation.

Plots of gene prediction accuracy in D. melanogaster reference transcripts built for GeneMarkS-T trained on sets of different types. The models were trained either on the set of D. melanogaster reference transcripts or on the sets of transcripts assembled by the five transcript assembly tools. Predictions made in reference transcripts were compared with annotation.

To assess performance of gene prediction methods in assembled transcripts we used the same five sets of assembled D. melanogaster transcriptions. First, we mapped the assembled transcripts to the corresponding reference transcripts ( 3) to detect and evaluate the differences. We used the results to divide the set of assembled transcripts into three groups: ‘concordant’, ‘conflicting’ and ‘not-aligned’ (see ‘Materials and Methods’ section and Supplementary Figure S2). Many assembled D. melanogaster transcripts fell into ‘conflicting’ category (from 17 to 47%, depending on the tool, see Supplementary Figure S3, ‘A’ bars) Cufflinks, Exonerate and Oases produced larger numbers of ‘conflicting’ transcripts than Augustus and Velvet. Multiple protein-coding regions were predicted more frequently in the ‘conflicting’ transcripts than in the ‘concordant’ transcripts (Supplementary Figure S4). Note, that for GeneMarkS-T events of prediction of multiple coding regions were registered prior to selecting ‘reported’ predictions with highest log-odd score. We have illustrated the distribution of events (multiple, single, none predictions) for GeneMarkS-T (Supplementary Figure S4). The distributions of the same events for the two other gene prediction tools show similar patterns (Table S2). Thus, all the tools predict single coding regions in ‘concordant’ assemblies with higher frequencies than in ‘conflicting’ ones.

To make unambiguous comparison of accuracy of gene prediction in ‘concordant’ transcripts we had to select the sets where gene finders make single gene predictions. As such surrogate sets we chose sets of ‘concordant’ assemblies where GeneMarkS-T predicted single protein-coding regions. Annotation of protein coding regions in these assembled transcripts was accomplished by transfer of the reference transcript annotation. In all the five test sets, GeneMarkS-T generated the largest number of TPs and the fewest number of FPs (Table 3).

In the sets of assembled transcripts where GeneMarkS-T predicted multiple coding regions we have observed high fractions of ‘conflicting’ transcripts (e.g. 90%, for the set of Cufflinks assembled transcripts). Thus, predicting multiple coding regions was an indicator of a higher chance for the transcript to be in the ‘conflicting’ category and to carry some discrepancies in the transcript assembly. Still, this observation should be taken with a caveat that multiple coding regions could appear in the ‘concordant’ transcript encoding alternative isoforms (as illustrated in Figure 5).

Very short coding regions are rare and are rarely predicted. Therefore, if an assembled transcript (complete or incomplete) is short it is likely that no gene will be predicted. Indeed, we observed that the gene finding tools did not predict genes in many transcripts assembled by the de novo methods Velvet and Oases (Supplementary Figure S3). Notably, many of these transcripts were too short (Supplementary Figure S1A).


Fichier supplémentaire 1 :

Includes 12 supporting figures and four supporting tables. A description of each is given within the file.

Additional file 2:

Performance of four transcriptome assemblers on the Edgren dataset. A table of which true positive breakpoint sequences were assembled by Trinity, Oases, TransABySS and SOAPdenovo-Trans on the Edgren dataset. Oases assembled the highest number of true positive breakpoints with 31.

Additional file 3:

Fusion genes in the BT-474, SK-BR-3, KPL-4 and MCF-7 cell lines. A list of the true positive fusion genes used in the validation of JAFFA on the Edgren and ENCODE dataset, along with a list of the probable true positives, and the fusion calls from JAFFA, FusionCatcher, SOAPfuse, defuse and TopHat-Fusion.

Additional file 4:

Fusion genes in the glioma dataset. A list of the true positive fusion genes, probable true positives and results from JAFFA, SOAPfuse, defuse and TopHat-Fusion for the gliomas dataset.

Additional file 5:

JAFFA commands. This script provides commands to reproduce the results from JAFFA and other tools shown in the manuscript.


Voir la vidéo: What is transcriptome? Introduction part 1: Learn from scratch for beginners. (Janvier 2022).