Informations

Défis de la construction de bibliothèques RNA-Seq : les biais de la fragmentation de l'ARN vs la fragmentation de l'ADNc


J'ai récemment regardé une présentation sur RNA-seq qui couvrait certains des choix que l'on peut faire en cours de route, et je n'ai pas bien compris l'un des choix en particulier. Vers le début du processus, vous pouvez choisir une méthode de fragmentation (par exemple, ARN : nébulisation, hydrolyse ; ADNc : sonication, traitement à la Dnase I). Chaque méthode de fragmentation est biaisée à sa manière, et les objectifs de l'étude peuvent influencer le choix de la méthode.

Je n'ai pas bien compris les raisons des biais caractéristiques (comment se passent-ils ?), et quand choisir une méthode par rapport à l'autre (quelles méthodes fonctionnent le mieux pour quels objectifs ?). En d'autres termes, quelles sont les explications mécanistes des biais caractéristiques qui résultent de chaque méthode, et quelle est l'importance du traitement en aval et des défis associés ?

Éditer:

Le cœur de cette question est de comprendre pleinement le choix particulier de quand et comment fragmenter. La raison pour laquelle je n'ai pas encore accepté de réponse est que les réponses fournies jusqu'à présent abordent la partie mécanisme de la question, mais pas les raisons ultimes de choisir l'une par rapport à l'autre.


Bonne question, beaucoup de choses sont encore à découvrir. Voici ce que l'on sait à ce jour :

Les méthodes de fragmentation basées sur des enzymes de restriction ne sont pas aléatoires.

La transcription inverse réalisée avec des oligomères poly dT, qui se lient aux queues poly-A 3', est fortement biaisée vers l'extrémité 3' des transcrits.

La transcription inverse avec des hexamères aléatoires entraîne une sous-représentation des extrémités 3'. Ceci est dû au nombre réduit de positions d'amorçage auxquelles l'enzyme transcriptase inverse peut démarrer la synthèse d'ADNc.


l'ARN messager a une queue poli A à son extrémité 5'. ainsi lorsque le poly-dT hybridera l'ARNm dans la transcription inverse, l'ADNc portera ce poly-dT à son extrémité 3'. Ainsi, comme la transcription inverse commencera toujours en 3' et de l'ARNm, il est plus probable que cette région soit meilleure rev. transcrit. Plus gros l'ARNm est plus gros sera cette préoccupation.

À propos du biais lié à la fragmentation de l'ARN, je ne sais pas, je peux supposer que les deux extrémités sont plus perdues en raison de la fragmentation elle-même

regardez cette vidéo des technologies de la vie http://www.youtube.com/watch?v=0MJIbrS4fbQ


Comparaison de l'ARN-Seq par capture poly (A), appauvrissement de l'ARN ribosomique et puce à ADN pour le profilage de l'expression

Le séquençage de l'ARN (RNA-Seq) est souvent utilisé pour le profilage du transcriptome ainsi que pour l'identification de nouveaux transcrits et d'événements d'épissage alternatifs. En règle générale, les bibliothèques d'ARN-Seq sont préparées à partir d'ARN total en utilisant l'enrichissement poly(A) de l'ARNm (ARNm-Seq) pour éliminer l'ARN ribosomique (ARNr), cependant, cette méthode ne parvient pas à capturer les transcrits non poly(A) ou partiellement dégradés ARNm. Par conséquent, un protocole d'ARNm-Seq ne sera pas compatible pour une utilisation avec des ARN provenant d'échantillons fixés au formol et inclus en paraffine (FFPE).

Résultats

Pour répondre au désir d'effectuer RNA-Seq sur des matériaux FFPE, nous avons évalué deux protocoles de préparation de bibliothèque différents qui pourraient être compatibles pour une utilisation avec de petits fragments d'ARN. Nous avons obtenu des ARN appariés Fresh Frozen (FF) et FFPE de plusieurs tumeurs et les avons soumis à différentes méthodes de profilage d'expression génique. Nous avons testé 11 échantillons de tumeurs mammaires humaines en utilisant : (a) des ARN FF par microarray, ARNm-Seq, Ribo-Zero-Seq et DSN-Seq (Duplex-Specific Nuclease) et (b) des ARN FFPE par Ribo-Zero-Seq et DSN -Séq. Nous avons également réalisé ces différents protocoles RNA-Seq en utilisant 10 tumeurs TCGA comme ensemble de validation.

Les données d'échantillons d'ARN appariés ont montré une concordance élevée dans la quantification des transcrits dans tous les protocoles et entre les ARN FF et FFPE. Dans les deux FF et FFPE, Ribo-Zero-Seq a éliminé l'ARNr avec une efficacité comparable à celle de l'ARNm-Seq, et il a fourni une couverture équivalente ou moins biaisée sur les extrémités 3' du gène. Par rapport à l'ARNm-Seq où 69 % des bases étaient mappées sur le transcriptome, DSN-Seq et Ribo-Zero-Seq contenaient significativement moins de lectures mappées sur le transcriptome (20 à 30 %) dans ces protocoles RNA-Seq, beaucoup sinon la plupart lectures mappées aux régions introniques. Environ 14 millions de lectures en mRNA-Seq et 45 à 65 millions de lectures en Ribo-Zero-Seq ou DSN-Seq ont été nécessaires pour atteindre les mêmes niveaux de détection de gènes qu'une puce à ADN Agilent standard.

Conclusion

Nos résultats démontrent que par rapport à l'ARNm-Seq et aux puces à ADN, Ribo-Zero-Seq offre une efficacité d'élimination d'ARNr équivalente, une uniformité de couverture, des lectures cartographiques basées sur le génome et une quantification de haute qualité constante des transcrits. De plus, Ribo-Zero-Seq et DSN-Seq ont une quantification cohérente des transcrits à l'aide d'ARN FFPE, ce qui suggère que l'ARN-Seq peut être utilisé avec des ARN dérivés de FFPE pour le profilage de l'expression génique.


Introduction

Les progrès récents dans le séquençage de l'ARN (RNA-Seq) ont fourni un moyen de caractérisation et de quantification rapides des transcriptomes. RNA-Seq implique le séquençage direct d'ADN complémentaires (ADNc) à l'aide de technologies de séquençage de nouvelle génération (NGS) à haut débit, suivi de la cartographie des lectures de séquençage sur le génome ou les ensembles de gènes de référence pour l'analyse de l'expression génique et la détection du polymorphisme. Par rapport à d'autres technologies telles que les puces à ADN basées sur l'hybridation et les méthodes basées sur le séquençage de Sanger, RNA-Seq fournit une compréhension plus complète de la complexité du transcriptome et la capacité de détecter une gamme dynamique de niveaux d'expression (Marioni et al., 2008 Wang et al. , 2009 Mader et al., 2011), permettant l'identification de nouveaux transcrits, petits ARN, SNP, produits d'épissage alternatif, transcrits sens et antisens, transcrits de fusion, et peut identifier les sites d'initiation de la transcription (Ozsolak et Milos, 2011).

Les plates-formes de séquençage de nouvelle génération utilisées pour RNA-Seq sont disponibles dans le commerce auprès d'Illumina, Roche, ABI, Helicos BioSciences, etc., et les entreprises améliorent continuellement leurs plates-formes pour augmenter les vitesses, la précision et la profondeur de séquençage à moindre coût. La réduction des coûts et les performances de séquençage élevées permettent des projets tels que les 100 génomes humains 1 à 10 millions de dollars et le Arabidopsis Projet 1001 génomes (Weigel et Mott, 2009). Même si la capacité de séquençage continue d'augmenter, les protocoles de préparation de bibliothèques d'échantillons, laborieux, longs et coûteux, restent une étape limitante. La préparation d'une bibliothèque de séquençage implique la production d'une collection aléatoire de fragments d'ADN modifiés par adaptateur prêts pour la séquence, avec une gamme spécifique de tailles de fragments. Bien que plusieurs procédures pour améliorer la préparation de la bibliothèque Illumina RNA-Seq aient été publiées (Quail et al., 2008 Nagalakshmi et al., 2010 Wilhelm et al., 2010), ces protocoles comportent encore plusieurs étapes laborieuses, notamment la précipitation à l'éthanol, les purifications sur colonne , et l'extraction sur gel pour le fractionnement granulométrique. En plus de prendre du temps, ces étapes comportent un risque élevé de contamination croisée et de mélange d'échantillons inhérent aux protocoles impliquant une manipulation d'échantillons individuels étendue. Récemment, Illumina a introduit une méthode à haut débit (kit de préparation d'échantillons d'ARN TruSeq) remplaçant ces étapes de purification par une méthodologie de nettoyage par réaction de billes magnétiques d'immobilisation réversible en phase solide (SPRI) (Hawkins et al., 1994 Lennon et al., 2010). En utilisant cette méthode, un seul technicien peut créer 96 bibliothèques à partir d'ARN total en 3 jours. Cependant, la quantité de multiplexage est limitée à 24 par le nombre de codes-barres disponibles. Des améliorations similaires peuvent également être observées dans les protocoles de Zhong et al. (2011) et Wang et al. (2011a).

Nous présentons ici plusieurs améliorations à la préparation d'échantillons Illumina pour le protocole RNA-Seq (Illumina Inc., San Diego, USA, Cat. # RS-100-0801) que nous avons apportées pour générer de l'ARN à haut débit et rentable -Bibliothèques Seq d'une manière plus robuste et reproductible, par rapport aux autres protocoles actuels. Nous avons intégré une méthode d'extraction directe d'ARNm à l'aide de billes oligo dT Dynabeads (Invitrogen, Carlsbad, CA, USA) ou de billes oligo dT Sera-Mag (Thermo Scientific, Indianapolis, IN, USA), qui conviennent à l'extraction d'ARN de divers végétaux et animaux. tissus. Un défi pour les protocoles de mise à l'échelle au format 96 puits est l'étape de fragmentation de l'ARN. Plus précisément, il est difficile de contrôler le degré de fragmentation chimique dans l'ARN en raison du temps d'incubation court, conduisant à une reproductibilité réduite, en particulier dans les formats 96 puits. Pour surmonter ce problème, nous avons utilisé la fragmentation enzymatique de l'ADNc. Nous avons également utilisé la méthodologie de nettoyage de réaction de billes magnétiques SPRI pour permettre la manipulation d'échantillons dans un format de 96 puits, similaire au protocole TruSeq et à celui de Zhong et al. (2011). De plus, pour réduire la durée du protocole et le nombre d'étapes de manipulation, nous avons appliqué un protocole de “on Beads” (Fisher et al., 2011) pour plusieurs réactions enzymatiques, notamment la réparation des extrémités, la queue A et la ligature de l'adaptateur. Ces changements réduisent le potentiel d'erreur humaine introduit pendant le processus de préparation de l'échantillon. Enfin, nous avons développé 96 adaptateurs à code-barres uniques pour offrir plus de flexibilité dans le multiplexage. Avec ces modifications et quelques autres petits ajustements, nous avons considérablement augmenté l'efficacité et la reproductibilité, et réduit le coût de préparation de la bibliothèque (de ߣ�×) par rapport aux autres méthodes actuellement disponibles. Notre méthode de préparation de bibliothèques d'ARN-seq (HTR) à haut débit permet à un seul chercheur de créer de manière reproductible 96 bibliothèques d'ARN-Seq, à partir de tissus, en moins de 3 jours. L'analyse des résultats de séquençage de nos bibliothèques a démontré que notre protocole fournit des données dont la qualité correspond ou dépasse celle de la méthode Illumina standard (IL) par composition de séquence, contamination par l'ARN ribosomique et détection de l'expression génique.


Résultats et discussion

Préparation et séquençage de la bibliothèque IVT-seq

Pour générer des bibliothèques IVT-seq (pour plus de détails, veuillez consulter le matériaux et méthodes section), nous avons produit des stocks de glycérol individuels abritant chacun un seul plasmide humain entièrement séquencé de la Mammalian Gene Collection (MGC) [17]. Ensuite, nous avons extrait l'ADN plasmidique et l'avons plaqué à 50 ng par puits dans des plaques de 384 puits. Nous avons mélangé le contenu de trois plaques de 384 puits contenant un total de 1 062 clones d'ADNc (Fichier supplémentaire 1), transformé ce mélange en bactéries et étalé les bactéries sous forme de colonies uniques. Après une incubation d'une nuit, nous avons gratté ces plaques, amplifié les bactéries pendant quelques heures en culture liquide et purifié les plasmides des bactéries en tant que pool (figure 1A). Ensuite, nous avons linéarisé les plasmides et utilisé la polymérase SP6 pour conduire in vitro transcription des séquences d'ADNc clonées (figure 1B). Après un traitement à la DNase I pour éliminer la matrice d'ADN et la purification de l'ARN, nous nous sommes retrouvés avec un pool de 1 062 ARN humains différents dérivés de plasmides entièrement séquencés.

Construction de librairies IVT-seq. (UNE) Préparation d'un pool de 1 062 plasmides d'ADNc humain. Les contenus de trois plaques de 384 puits contenant des plasmides MGC ont été regroupés. Le pool a été amplifié par transformation en Escherichia coli, et les clones résultants ont été purifiés et regroupés. (B) Génération de transcriptions IVT. Un pool de plasmides MGC a été linéarisé et utilisé comme matrice pour un in vitro réaction de transcription. Les enzymes et les nucléotides non incorporés ont été purifiés, laissant un pool de transcrits polyA. (C) Création de librairies IVT-seq. Les quantités listées d'ARN IVT ont été mélangées avec de l'ARN total de foie de souris pour créer six pools avec des quantités finales d'ARN de 1 ug. L'ARN ribosomique a été épuisé dans ces pools à l'aide du kit Ribo-Zero Gold. L'ARN IVT et l'ARN de souris sont maintenant présents dans des pools aux rapports indiqués, suite à l'épuisement de l'ARNr de l'ARN total de souris. Ces pools ont été utilisés pour générer des bibliothèques d'ARN-seq à l'aide du kit/protocole TruSeq d'Illumina. L'ensemble de ce processus a été effectué en double. Les bibliothèques répliquées ont été regroupées séparément et séquencées dans des voies HiSeq 2000 distinctes (deux voies au total). IVT, in vitro transcrit MGC, Mammalian Gene Collection.

Pour approximer ce qui se passe dans une réaction ARN-seq totale, nous avons soumis cet ARN IVT à l'épuisement de l'ARNr, puis préparé des bibliothèques à l'aide du protocole Illumina TruSeq (figure 1C, IVT uniquement). Pour tenir compte des effets possibles du transporteur, nous avons également mélangé l'ARN IVT avec diverses quantités d'ARN total de souris dérivé du foie. L'ajout de l'ARN de souris a donné à ces échantillons une plus grande diversité (transcriptions d'environ 10 000 gènes contre 1 062) et ressemblait davantage à un échantillon biologique réel. De plus, en ajoutant de l'ARN de fond d'une espèce différente (souris) que l'ARN IVT (humain), nous avons facilité la différenciation entre les transcrits IVT et les séquences de souris lors de l'analyse en aval. Étant donné que l'ARN IVT ne contenait pas de séquences d'ARNr alors que l'ARN de souris en contenait, la quantité d'ARN de souris serait significativement réduite par l'étape d'épuisement de l'ARNr. Pour tenir compte de cela, nous avons mélangé l'IVT et l'ARN de souris de telle sorte que, après l'épuisement de l'ARNr, nous aurions des pools finaux avec des ratios IVT:souris de 1:1, 1:2 et 1:10. Enfin, pour tenir compte des ARN de souris potentiellement mappés sur le génome humain de référence et nos séquences IVT, nous avons préparé un pool composé d'ARN de souris seul. Nous avons regroupé les six bibliothèques résultantes et les avons séquencées à l'aide d'un Illumina HiSeq 2000. Nous avons effectué tout ce processus en double.

Cartographie et couverture des données IVT-seq

Après le séquençage et le démultiplexage, nous avons aligné toutes les données sur le génome humain de référence (hg19) à l'aide du RNA-seq Unified Mapper (RUM) [14]. Pour toutes les analyses, nous n'avons utilisé que les données des lectures mappées de manière unique sur la référence, à l'exclusion de tous les multi-mappers (données contenues dans les fichiers RUM_Unique et RUM_Unique.cov). Sur les 1 062 transcrits IVT originaux, nous avons trouvé 11 alignés sur plusieurs loci génomiques, tandis que 88 alignés sur des loci qui se chevauchent. Pour éviter tout effet de confusion dans nos analyses, nous avons filtré ces transcriptions de toutes les analyses, nous laissant 963 transcriptions IVT non superposées et alignées de manière unique. Nous avons constaté une excellente corrélation dans les niveaux d'expression entre les réplicats (niveau de transcription R 2 entre les réplicats >0.95 Fichier supplémentaire 2 : Figure S1A). Deuxièmement, au moins 90 % des 963 transcrits IVT ont été exprimés avec des fragments par kilobase d'exons par million de lectures cartographiées (FPKM) valeurs ≥5 dans tous les ensembles de données IVT-seq, à l'exception de la souris uniquement (tableau 1). Dans les échantillons IVT uniquement, plus de 80 % des séquences IVT étaient exprimées au-dessus de 100 FPKM (Fichier supplémentaire 2 : Figure S1B). Parce que nous avons préparé les plasmides MGC et les transcrits IVT en tant que pools, il est probable que les transcrits IVT présentant une couverture faible ou nulle étaient initialement présents à de faibles concentrations de plasmides avant les étapes de transformation et d'IVT. En utilisant la technique IVT-seq, nous avons pu détecter spécifiquement la grande majorité des transcrits IVT humains avec une couverture élevée à la fois en l'absence et en présence de l'ARN de souris de fond.

Bien que nous voyions des lectures alignées sur les transcriptions IVT humaines dans les données de souris uniquement, ces transcriptions représentent collectivement environ 2 % des lectures (tableau 1). Ces transcriptions avec une couverture plus élevée sont probablement le résultat de lectures de souris alignées sur des séquences humaines très similaires. Nous avons exclu ces séquences de nos analyses.

Variation intra-transcript de la couverture ARN-seq des transcrits IVT

Considérons d'abord les données IVT uniquement. Étant donné que ces transcrits ont été générés à partir d'une réaction IVT utilisant des séquences d'ADNc, ces données ne sont pas affectées par l'épissage ou une autre régulation post-transcriptionnelle. Ainsi, la plupart des régions des transcriptions devraient être « exprimées » et présentes à des niveaux similaires. Les exceptions seraient les séquences répétitives qui correspondent à plusieurs emplacements du génome et peuvent être mal représentées, et les extrémités des ADNc, qui sont sujettes à un biais de fragmentation. Pour tenir compte de cela, nous avons créé un jeu de données simulé qui modélise le processus de fragmentation et ne s'écarte des données uniformes que par le caractère aléatoire induit par la fragmentation. Nous avons généré deux de ces ensembles de données à l'aide du Benchmarker for Evaluating the Effectiveness of RNA-Seq Software (BEERS) [14]. Le premier ensemble de données contenait tous les transcrits IVT exprimés à peu près au même niveau d'expression (environ 500 FPKM). Pour le second, nous avons utilisé les valeurs FPKM des échantillons IVT uniquement comme graines, créant un ensemble de données simulé avec des niveaux d'expression correspondant étroitement aux données réelles (Fichier supplémentaire 3 : Figure S2). Ces ensembles de données sont respectivement appelés simulés et simulés en quantité appariée (QM). Les données simulées fournissent un résultat idéal, tandis que les données QM nous permettent de contrôler les artefacts résultant du niveau d'expression (par exemple, les transcrits avec une expression inférieure peuvent montrer plus de variabilité). Ensuite, nous avons aligné les deux jeux de données simulés à l'aide de RUM, avec les mêmes paramètres que pour les données biologiques. Ainsi, les deux ensembles de données simulés servent également de contrôles pour tout artefact introduit par l'alignement (par exemple, une faible couverture dans les régions répétées). Pour plus de détails sur la création de données simulées, consultez le matériaux et méthodes section.

En utilisant les données IVT dérivées du transcrit BC015891 comme exemple représentatif, le tracé de couverture théorique idéal à partir des données simulées montre une couverture presque uniforme sur toute la longueur du transcrit, sans aucun des pics et vallées extrêmes caractéristiques des ensembles de données biologiques (Figure 2A) . Cependant, nos données observées ont montré un degré élevé de variabilité, avec des pics et des vallées au sein d'un exon (figure 2B). De plus, ces modèles étaient reproductibles dans nos réplicats (Fichier supplémentaire 4 : Figure S3). Nous avons vu de nombreux autres cas de changements extrêmes de couverture : plus de 50 % des transcriptions IVT ont montré des changements plus du double de la couverture intra-transcription attribuables à la préparation et au séquençage de la bibliothèque (tableau 2 et fichier supplémentaire 5 : figure S4). Par exemple, BC009037 a montré des baisses soudaines à des niveaux d'expression extrêmement faibles dans ses deux exons (figure 2C). Les deux ensembles de données simulés n'ont montré aucun de ces modèles, ce qui indique que cette variabilité de la couverture n'est pas le résultat d'artefacts d'alignement. De plus, l'absence de ce modèle dans les données simulées par QM indique que ces différences de variation de couverture n'étaient pas dues au bruit d'échantillonnage introduit par les transcriptions avec une couverture faible ou élevée. Dans le cas de BC016283, les pics et les creux de couverture ont conduit à des différences de plus de cinq fois dans les niveaux d'expression entre les exons (figure 2D). Encore une fois, ces modèles étaient reproductibles à travers les réplicats (Fichier supplémentaire 4 : Figure S3). La polymérase SP6 ne peut pas tomber puis se rattacher à un point ultérieur du transcrit, laissant une région non transcrite. Par conséquent, étant donné que ces modèles ont montré des creux suivis de pics, ils ne peuvent pas être le résultat d'artefacts de in vitro transcription. De plus, nous avons séquencé les produits IVT directement et avons trouvé que la grande majorité était transcrite avec peu ou pas de biais. Prises ensemble, ces données suggèrent que ces modèles de couverture sont principalement le résultat de biais techniques introduits lors de la construction de la bibliothèque, plutôt que de la biologie. Ces résultats sont cohérents avec une étude précédente qui utilisait l'ARN IVT comme standard dans les expériences d'ARN-seq [16], suggérant que notre méthodologie IVT-seq est adaptée pour identifier la variabilité technique des données de séquençage.

Variations intra-transcript de la couverture ARN-seq. (UNE) Couverture RNA-seq simulée pour un transcrit IVT représentatif (BC015891). Le graphique de couverture ARN-seq (noir) est affiché en fonction du modèle de gène (vert), car il est mappé sur le génome de référence. Les blocs correspondent aux exons et les lignes indiquent les introns. Les chevrons dans les lignes introniques indiquent la direction de la transcription. Les nombres sur l'axe des y font référence à la profondeur de lecture de RNA-seq à une position de nucléotide donnée. (B) Le graphique de couverture ARN-seq réel pour BC015891 dans l'échantillon IVT uniquement. Diagrammes de couverture représentatifs pour les transcriptions IVT (C) BC009037 et (RÉ) BC016283 sont affichés selon les mêmes conventions utilisées ci-dessus. Toutes les transcriptions sont affichées dans la direction 5ʹ à 3ʹ.

Variation entre les échantillons de la couverture ARN-seq des transcrits IVT

En plus de cette variabilité au sein des transcrits, nous avons également trouvé de nombreuses régions de transcrits montrant une extrême variabilité de la couverture entre les échantillons (figure 3). Par exemple, le sixième exon de BC003355 variait énormément par rapport au reste du transcrit dans toutes les dilutions IVT:souris. Fait intéressant, le modèle global de variation par rapport au reste du transcrit à travers les dilutions a été maintenu entre les répétitions. Presque aucune lecture dans la carte de l'échantillon de souris uniquement sur ce transcrit, ce qui élimine la possibilité que cette variabilité soit due à un alignement incorrect de l'ARN de souris.

Variations inter-échantillons de la couverture ARN-seq. Graphiques de couverture ARN-seq dans tous les échantillons pour les exons 4 à 11 du transcrit IVT BC003355. Les rectangles noirs identifient l'exon six, qui montre une extrême variabilité de la couverture par rapport au reste du transcrit lorsqu'il est visualisé sur tous les échantillons. Le rapport ARN IVT sur ARN de souris est indiqué à gauche des graphiques de couverture de chaque échantillon. Les parcelles de couverture (rouge pour le premier réplicat bleu pour le deuxième réplicat) sont affichées en fonction du modèle de gène (noir), tel qu'il est mappé sur le génome de référence. Les blocs dans le modèle de gène correspondent aux exons et les lignes indiquent les introns. Les chevrons à l'intérieur des lignes introniques indiquaient la direction de la transcription. Les nombres sur les axes y font référence à la profondeur de lecture de RNA-seq à une position de nucléotide donnée.

Y compris BC003355, nous avons trouvé 86 régions de couverture élevée et imprévisible (hunc) réparties sur 65 transcriptions (Fichier supplémentaire 6). Par conséquent, plus de 6 % des 963 transcrits IVT contenaient des régions présentant des variations sauvages mais reproductibles de la couverture ARN-seq entre les échantillons. Lors de l'identification de ces régions hunc, nous avons utilisé un filtre en deux étapes pour éliminer les régions variables résultant des lectures de souris mappées sur des séquences humaines très similaires. Tout d'abord, nous avons éliminé toutes les régions hunc provenant de transcrits avec FPKM ≥5 dans l'un ou l'autre ensemble de données de souris uniquement. Ensuite, pour tenir compte du désalignement localisé des lectures de souris, nous avons filtré toutes les régions hunc avec une couverture moyenne ≥10 dans l'un ou l'autre ensemble de données de souris uniquement. Nous avons également supprimé ces régions hunc avec une couverture de souris uniquement ≥10 dans les 100 paires de bases (pb) flanquantes de chaque côté. Compte tenu des critères rigoureux que nous avons utilisés pour identifier ces régions de forte densité (voir matériaux et méthodes section pour plus de détails), il est probable qu'il s'agisse d'une sous-estimation. Pour aborder la possibilité que les ARN de souris puissent interagir avec des ARN humains homologues et interférer avec eux dans trans, nous avons testé les séquences entourant ces régions à l'aide de la suite MEME [18], mais nous n'avons trouvé aucun motif de séquence que ces régions ont en commun. De plus, la profondeur de couverture dans ces régions n'a pas suivi une relation linéaire avec l'augmentation de l'ARN de souris, ce qui suggère qu'il ne s'agit pas simplement d'une interaction directe avec l'ARN de fond. Il n'y a pas de cause claire pour ces régions hunc, d'autant plus que nous avons préparé tous les échantillons à partir du même pool d'ARN IVT et que la seule différence entre les échantillons était les rapports relatifs d'ARN IVT à l'ARN de foie de souris. Nous avons également recherché des régions hunc divergentes entre les deux réplicats, mais n'en avons trouvé aucune. Si de telles régions existent, elles pourraient être identifiées et surmontées en créant des bibliothèques en double. Les régions hunc que nous avons identifiées ci-dessus avec des modèles d'expression maintenus entre les réplicats présentent un plus grand défi, car elles n'ont pas pu être identifiées et filtrées en créant des bibliothèques en double. Ceci est particulièrement problématique pour l'utilisation de valeurs d'expression au niveau de l'exon pour identifier des événements d'épissage alternatifs ou une expression différentielle. La variation intra-transcript et inter-échantillons que nous voyons dans nos données IVT-seq suggère que la génération de bibliothèques introduit de forts biais techniques, qui pourraient confondre les tentatives d'étude de la biologie sous-jacente.

Sources de variabilité dans la couverture RNA-seq

Il existe trois sources potentielles de biais technique dans la préparation des bibliothèques : la biologie moléculaire spécifique à l'ARN (fragmentation de l'ARN, transcription inverse), la méthode de sélection de l'ARN (depletion de l'ARNr, sélection polyA) et la biologie moléculaire spécifique au séquençage (ligature de l'adaptateur, enrichissement de la bibliothèque, pont PCR). Pour identifier les biais introduits uniquement par la biologie moléculaire spécifique au séquençage, nous avons créé une bibliothèque d'ADN-seq à partir des mêmes plasmides MGC utilisés comme modèles pour les bibliothèques IVT-seq (Fichier supplémentaire 7 : Figure S5). Ce faisant, nous avons sauté les étapes spécifiques à la biologie moléculaire IVT ou ARN. Nous avons également préparé deux bibliothèques IVT-seq supplémentaires en utilisant la sélection polyA ou aucune sélection, au lieu de l'épuisement de l'ARNr. En comparant nos données de bibliothèque de plasmides et les données IVT-seq à l'aide de diverses méthodes de sélection, nous avons pu identifier quels modèles de couverture étaient le résultat de la biologie moléculaire spécifique à l'ARN, de la méthode de sélection de l'ARN ou d'un aspect commun du protocole de génération de bibliothèque.

Nous avons séquencé la bibliothèque de plasmides à l'aide d'un Illumina MiSeq et aligné les données résultantes sur le génome humain de référence en utilisant la même méthode que les bibliothèques IVT-seq. Dans ces données de plasmide, nous avons vu 924 des séquences de clones d'ADNc avec des valeurs FPKM 5, contre environ 870 dans les deux échantillons IVT uniquement (tableau 1). Cette petite baisse de couverture était probablement due au fait que l'ARN IVT passe par plus d'étapes de regroupement pendant la construction de la bibliothèque que les plasmides. De plus, les plasmides ne sont pas affectés par les efficacités de transcription et de transcription inverse. De plus, les données plasmidiques mappées sur les séquences d'ADNc avec une couverture moyenne normalisée de 42,08, ce qui se situe dans la plage de valeurs de couverture que nous voyons pour les échantillons IVT-seq. Nous avons séquencé les bibliothèques de sélection sans sélection et polyA sur un HiSeq 2500. Ces données montrent également des valeurs de couverture de clone d'ADNc similaires aux autres bibliothèques IVT-seq.

Les données de plasmide représentent «l'entrée» dans la réaction IVT et les données de non-sélection représentent la mesure la plus proche de sa sortie directe. En mesurant le rapport 3'/5' en profondeur de couverture pour chaque transcrit IVT, nous avons pu évaluer la processivité de la polymérase SP6. Dans une réaction parfaitement processive, ce rapport 3'/5' serait de 1, indiquant que la polymérase n'est pas tombée de la matrice d'ADNc et a conduit à la formation de produits tronqués. Les rapports médians 3′/5′ pour le plasmide et aucune donnée de sélection étaient de 1 et de 0,98, respectivement, indiquant que l'arrêt prématuré de la réaction IVT n'était pas un facteur dans nos analyses.

Effet de différentes méthodes de sélection d'ARN sur les modèles de couverture

Notre analyse est illustrée par un examen des graphiques de couverture pour BC003355 dans tous nos différents ensembles de données. Le degré élevé de variabilité que nous avons noté dans le graphique de couverture de ce gène à partir de nos données appauvries en ARNr était absent dans les données de non-sélection et de plasmide (figure 4A). Alors que les données polyA ont également montré moins de pics et de vallées que les données d'ARN-seq total appauvri en ARNr, elles ont été marquées par le biais 3′ bien documenté. Ces données suggèrent que l'étape d'épuisement de l'ARNr est probablement responsable d'une grande partie des biais de couverture observés.

Sources de biais dans la couverture RNA-seq. (UNE) Graphiques de couverture ARN-seq pour le transcrit IVT BC003355 à partir de données simulées (noir), plasmide (bleu), aucune sélection (vert), appauvri en ARNr (rouge) et polyA (orange). Le modèle de gène est affiché en noir, sous tous les graphiques de couverture. Les blocs correspondent aux exons et les lignes indiquent les introns. Les chevrons dans les lignes introniques indiquent la direction de la transcription. (B) Distributions des coefficients de variation entre les données affichées ci-dessus, avec l'ajout des données simulées par QM (gris). Notez que bien que le graphique soit coupé à un coefficient de variation de 1,3, les queues des distributions appauvries en ARNr et polyA s'étendent à 2,13 et 2,7, respectivement. (C) Tailles d'effet pour les différences de distribution des coefficients de variation entre les bibliothèques de séquençage et les données simulées. Les tailles d'effet sont calculées comme les ratios par transcription des coefficients de variation entre une bibliothèque donnée et l'ensemble de données simulé. QM, quantité correspondante.

Pour quantifier la variabilité pour chaque méthode de sélection, nous avons calculé le coefficient de variation au niveau de base unique de la couverture pour tous les transcrits IVT dans chacun de ces ensembles de données (figure 4B). En utilisant un test de somme de rang de Wilcoxon (plasmide n = 924, aucune sélection n = 870, ARNr appauvri n = 869), nous avons constaté que les données appauvries en ARNr présentaient une variabilité significativement plus élevée que les données sans sélection et plasmide (P <2.2e -16 ). De plus, les bibliothèques appauvries en ARNr et polyA étaient plus de 60 % plus variables en moyenne que la bibliothèque de plasmides (figure 4C). Cela suggère qu'une partie importante de la variabilité observée de la couverture entre les transcrits dans les données IVT-seq est le résultat de la biologie moléculaire spécifique à l'ARN, en particulier l'étape de sélection de l'ARN. De plus, après avoir pris en compte le biais introduit par les séquences elles-mêmes (données plasmidiques) et le biais introduit par la réaction IVT (pas de données de sélection), nous avons constaté que 50% des transcrits présentaient une variation de deux fois et 10% avaient une variation de 10 fois dans l'intervalle. expression du transcrit (tableau 2 et fichier supplémentaire 5 : figure S4). Bien qu'il soit bien apprécié que la sélection polyA introduit un biais, nous avons constaté que les données appauvries en ARNr introduisaient autant sinon plus. Aucun des deux ensembles de données simulés n'a montré de transcriptions avec un changement double ou supérieur dans l'expression intra-transcription. Encore une fois, cela suggère que les variations observées au sein de la transcription ne sont pas le résultat d'artefacts d'alignement ou d'échantillonnage dus à une expression faible ou élevée. Une source de biais communément reconnue provient de l'amorçage aléatoire lors de la préparation de la bibliothèque [10]. Lorsque nous avons examiné les différentes bibliothèques, nous avons constaté que des fragments de toutes les données RNA-seq présentaient des fréquences de nucléotides caractéristiques d'un biais d'amorçage aléatoire (Fichier supplémentaire 8 : Figure S6). Comme prévu, les données du plasmide n'ont montré aucun biais de ce type, car elles étaient dérivées directement de l'ADN et ne nécessitaient pas d'étape de génération d'ADNc. Cependant, les différences significatives entre toutes les bibliothèques d'ARN suggèrent que le biais de l'amorçage aléatoire n'est pas le seul facteur. Le plasmide et aucune donnée de sélection ne contiennent toujours une bonne quantité de variabilité lorsqu'ils sont visualisés à côté des données simulées (figure 4A noire). Lorsque nous avons examiné l'ensemble de données dans son intégralité, les données relatives au plasmide et à l'absence de sélection présentaient une variation significativement plus élevée que l'un ou l'autre des ensembles de données simulés (données simulées par le test de somme des rangs de Wilcoxon n = 963, données simulées par QM n = 869, plasmide n = 924, aucune sélection n = 870 P <2.2e -16 ). Ces données suggèrent que la biologie moléculaire spécifique au séquençage commune à toutes les bibliothèques que nous avons préparées (ligature d'adaptateur, amplification de bibliothèque par PCR) est également responsable d'une partie de la variabilité de la couverture observée et du biais de séquençage.

Les biais associés aux caractéristiques des séquences dépendent de la méthode de sélection de l'ARN

Compte tenu de ces différences significatives dans la variabilité de la couverture, nous avons cherché à identifier les caractéristiques de séquence qui pourraient contribuer à ce biais. Nous avons considéré trois caractéristiques de séquence quantifiables : l'entropie de l'hexamère, la teneur en GC et la similarité de séquence avec l'ARNr (voir matériaux et méthodes for a detailed description of these metrics). For each sequencing strategy (plasmid, no selection, rRNA-depleted, polyA), we tested if any of the three sequence characteristics had a significant effect on variability in sequencing coverage, as measured by the coefficient of variation. While we are primarily focused on coverage variability as an indicator of sequencing bias, we also looked at depth of coverage, as measured by FPKM.

For each sequencing strategy, we sorted the transcripts by coverage variability or depth. Next, we selected the 100 most and 100 least extreme transcripts from each list. We compared the values of the sequence characteristics between the 100 most and 100 least extreme transcripts using a Wilcoxon rank-sum test. Significant P-values indicate a significant association of the sequence characteristic with coverage variability and/or depth. The results of our analysis are displayed as box-plots (Figure 5 and Additional file 9: Figure S7).

Effects of sequence characteristics on coverage variability. Distributions of (UNE) hexamer entropy, (B) GC-content, and (C) rRNA sequence similarity for the 100 transcripts with the highest and lowest coefficients of variation for transcript coverage from the plasmid, no selection, rRNA-depleted, and polyA libraries. Asterisks indicate the significance of a Wilcoxon signed-rank test comparing values for the listed sequence characteristics between each pair of groups from the same sample. *P <0.05 **P <0.01 ***P <0.001.

To check for any confounding effects between coverage depth and variability, we tested the least and most expressed transcripts for any correlations with variability in coverage (Additional file 10: Figure S8). The polyA library showed a significant correlation (P <2.2e -16 ) between coverage variability and depth, which indicates sequence features could be affecting coverage through variability (or vice versa). The rRNA-depleted data showed a slight, significant correlation (P = 0.04933). It is possible some feature of RNA selection affects both variability and coverage, given that we saw no significant correlations for the two remaining samples. This indicates that coverage variability and depth are independent for the plasmid and no selection data.

All three sequence characteristics had a significant association with variability and depth-of-coverage in at least one of the sequencing strategies. In particular, lower hexamer entropy, a measure of sequence complexity [19–21], was strongly associated with higher variance in all of the RNA libraries (no selection P = 4.712e -05 rRNA depletion P = 3.956e -11 polyA P = 0.003921 Figure 5A). This suggests that bias associated with hexamer entropy is due partially to RNA-specific procedures in library preparation. Furthermore, an association with lower hexamer entropy indicates there are more repeat sequences in the transcripts with higher variability. This could be indicative of complex RNA secondary structures, as repeated motifs could facilitate hairpin formation. Furthermore, the absence of this association from the plasmid data suggests that this observation was not due to mapping artifacts. The plasmid data contained the same sequences as the RNA-seq data, and would be subject to the same biases introduced by our exclusion of multi-mapped reads.

Higher GC-content was strongly associated with lower coverage variability in the no selection and polyA data (P = 5.627e -13 P = 4.914e -05 Figure 5B), suggesting that the effects of GC-bias on within-transcript variability could arise, in part, due to some RNA-specific aspects of library preparation. Also, it appears that GC-bias was not a significant contributing factor to either depth of coverage or the extreme variability in the rRNA-depleted data. Meanwhile, lower GC-content was associated with higher coverage in the plasmid data (P = 3.776e -05 ), and lower coverage depth in the no selection and polyA libraries (no selection P = 8.531e -05 polyA P = 0.0009675 Additional file 9: Figure S7B). Given that this trend switched directions between the plasmid library and the RNA libraries, this also suggests that some RNA-specific aspect of library preparation is introducing GC-bias distinct from the high GC-bias associated with Illumina sequencing [22].

Interestingly, higher rRNA sequence similarity was associated with higher coverage variability in the rRNA-depleted library (P = 9.006e -05 ) and lower variability in the no selection library (P = 0.0367 Figure 5C). It is unsurprising that similarity to rRNA sequences contributed to variability in the rRNA-depleted data, given that rRNA depletion is based upon pair-binding between probes and rRNA sequences. While it is unclear why this trend was reversed in the no selection library, it is striking given the significant increase in within-transcript variability between the no selection and rRNA-depleted libraries (Figure 4B). Furthermore, we saw a slight but highly significant correlation (Pearson R 2 = 0.308452 P <2.2e -16 ) between a transcript sequence’s similarity to rRNA and the magnitude of the difference in coverage between the no selection and rRNA-depleted libraries (Additional file 11: Figure S9 and Additional file 12). While the majority of the factors contributing to the extreme bias in sequence coverage we saw in the rRNA-depleted data remain unclear, our data suggest this could be partially due to depletion of sequences homologous to rRNA.

Taken together, our data demonstrate the utility and potential of the IVT-seq method to identify sources of technical bias introduced by sequencing platforms and library preparation protocols.


Discussion

RNA-seq has become a standard method for gene expression quantification and in most cases the sequencing library preparation involves amplification steps. Ideally, we would like to count the number of RNA molecules in the sample and thus would want to keep only one read per molecule. A common strategy applied for amplification correction in SNP-calling and ChIP-Seq protocols 23,24 is to simply remove reads based on their 5′ ends, so called read duplicates. Here, we show that this strategy is not suitable for RNA-seq data, because the majority of such SE-duplicates is likely due to sampling. For highly transcribed genes, it is simply unavoidable that multiple reads have the same 5′ end, also if they originated from different RNA-molecules. We find that only

20% (Smart-Seq) of the read duplicates cannot be explained by a simple sampling model with random fragmentation. This fraction decreases even more, if we factor in that the fragmentation of mRNA or cDNA during library preparation is clearly non-random, as evidenced by a strong correlation between the 5′ read positions of the ERCC-spike-ins across samples. Because local sequence content has little or no detectable effect on fragmentation, we cannot predict fragmentation, but we can quantify the observed effect. For example, we find that a fragmentation bias that halves the number of break points can fit the observed proportion of duplicates for TruSeq libraries well. For the Smart-Seq datasets, fragmentation biases would have to be much higher to explain the observed numbers of read duplicates. Furthermore, the fit between model estimates and the observed duplicate fractions is worse than for the TruSeq data and the model estimates for fragmentation bias are also inconsistent between the datasets (38.5 for the UHRR and 8 for the scHCT116).

Since computational methods cannot distinguish between fragmentation and PCR duplicates, the removal of read duplicates could introduce a bias rather than removing it. Using the ERCC-spike-ins, we can indeed show that removing duplicates computationally does not improve a fit to the known concentrations, but rather makes it worse, especially if only single-end reads are available (Fig. 5). This is in line with our observation that most single end duplicates are due to sampling and fragmentation. Hence, removing duplicates is similar to a saturation effect known for microarrays 25,26,27 .

Moreover, the Smart-Seq protocol, which was designed for small starting amounts, involves PCR amplification before the final fragmentation of the sequencing library. Thus in the case of Smart-Seq, computational methods cannot identify PCR duplicates that occur during the pre-amplification step. When we use unique molecular identifiers (UMIs), we find that 66% of the reads are PCR duplicates and only 34% originate from independent mRNA molecules. In contrast, when using paired-end mapping for a comparable Smart-Seq library, we identify 13% as duplicates and 87% as unique. This might in part be due to the fact that in UMI-Seq we sequence mainly 3′ ends of transcripts, thus decreasing the complexity of the library, which in turn increases the potential for PCR duplicates for a given sequencing depth (Fig. 4a, Supplementary Figure S1). However, it is unlikely that library complexity can explain the 53% difference in duplicate occurrence. This difference is more likely to be due to PCR-duplicates that are generated during pre-amplification and thus remain undetectable by computational means.

All in all, computational methods are limited when it comes to removing PCR-duplicates, but how much noise or bias do PCR duplicates introduce? In other words, we want to know how PCR-duplicates impact the power and the false discovery rate for the detection of differentially expressed genes. Both, power and FDR, are determined by the gene-wise mean expression and dispersion. Based on simulated differential expression using the empirically determined mean and dispersion distributions, we find that computational removal of duplicates has either a negligible or a negative impact on FDR and power and we therefore recommend not to remove read duplicates. In contrast, if PCR duplicates are removed using UMIs, both FDR and power improve. Even though the effects in the bulk data analysed here are relatively small: FDR is improved by 4% and the power by 2%, UMIs will become more important when using smaller amounts of starting material as it is the case for single-cell RNA-seq 6,28 .

The major differences in power are between the datasets with the TruSeq and the UMI-seq data achieving a power of around 80%, the UHRR-Smart-Seq 52% and the single cell Smart-Seq data (scHCT116) only 27%. Note that this apparently bad performance of the single cell Smart-Seq data is at least in part due to an unfair comparison. While all the other datasets were produced using commercially available mRNA and thus represent true technical replicates, the single cell data necessarily represent biological replicates and thus are expected to have a larger inherent variance and thus lower power.

However, also the UHRR Smart-Seq bulk data achieves with 52% a much lower power than the other bulk datasets. One possible explanation for the differences in power is the total number of PCR-cycles involved in the library preparation. With every PCR-cycle the power to detect a log 2-fold change of 0.5 appears to drop by 2.4% (Fig. 6c). The only exception is the UMI-seq dataset, that gives a power of 81%, even if duplicates are not removed, which is comparable to the power reached with TruSeq data despite the UMI-seq method having 12 more PCR-cycles. Technically UMI-seq is most similar to the Smart-Seq method. The biggest difference between the two methods is that all UMI-seq libraries are pooled before PCR-amplification, suggesting that the PCR-noise is due to the different PCR-reactions and not due to amplification efficiency per-se.

We conclude that computational removal of duplicates is not recommendable for differential expression analysis and if sufficient starting material is available so that only few PCR-cycles are necessary, the loss in power due to PCR duplicates is negligible. However, if more amplification is needed, power would be improved if all samples are pooled early on and for really low amounts as for single cell data also the gain in power that is achieved by removing PCR-duplicates using UMIs will become important.


The effect of methanol fixation on single-cell RNA sequencing data

Single-cell RNA sequencing (scRNA-seq) has led to remarkable progress in our understanding of tissue heterogeneity in health and disease. Recently, the need for scRNA-seq sample fixation has emerged in many scenarios, such as when samples need long-term transportation, or when experiments need to be temporally synchronized. Methanol fixation is a simple and gentle method that has been routinely applied in scRNA-seq. Yet, concerns remain that fixation may result in biases which may change the RNA-seq outcome.

Researchers from the Hong Kong University of Science and Technology adapted an existing methanol fixation protocol and performed scRNA-seq on both live and methanol fixed cells. Analyses of the results show methanol fixation can faithfully preserve biological related signals, while the discrepancy caused by fixation is subtle and relevant to library construction methods. By grouping transcripts based on their lengths and GC content, the researchers found that transcripts with different features are affected by fixation to different degrees in full-length sequencing data, while the effect is alleviated in Drop-seq result. This deep analysis reveals the effects of methanol fixation on sample RNA integrity and elucidates the potential consequences of using fixation in various scRNA-seq experiment designs.

Basic evaluation of fixation effect on sequencing data

(UNE) Workflow and experimental scheme. (B) Size distributions of cDNA libraries. Traces from single-cell libraries were merged to obtain a general pattern for live (top) and fixed (bottom) samples. Although the intensity of the

1500bp peak is diminished in fixed cells, there is no visible degradation. (C) Correlation matrix showing the transcriptome similarity of cells randomly chosen from live and fixed samples. The upper triangle of the matrix shows the Pearson correlation coefficient and the bottom triangle visualized correlation trend. Correlations are consistently high for both inter- and intra-treatment comparisons of live vs. fixed. There is no obvious bias revealed by measuring correlation between single-cell transcriptomes for all pairwise comparisons. (RÉ) Correlation factors of all single cells were calculated pairwise and clustered by Euclidean distance. Correlations are consistently high for both inter- and intra-treatment comparisons of live vs. fixed (R2 >0.7). The mixed annotation bar indicates the transcriptome similarities do not distinguish cell treatments during sample preparation.


CONCLUSION

Bias in general describes systematic errors that reflect method-related distortion from the truth. Bias might be detected easily when data generated by different methods are subjected to analytical comparison. Because it is systematic in nature, bias is not subject to variation in repeated experiments. Here, based on biochemical and deep sequencing data analyses, we have presented various steps in sRNA-seq protocols likely to cause severe distortions in the relative expression levels of individual sRNAs.

The efficacy of RNA end-modification to permit cDNA construction depends not only on differences in strategies and reagents but also on changes in, among others, buffer compositions and additives that could increase molecular crowding. Therefore, meta-analysis of data generated with non-identical methods, reagents or even buffers is likely to be error-prone ( 36, 39, 63, 78). To achieve comparable results and to evaluate relative changes in the quantification of the resulting data, experiments and analyses should be conducted under strictly identical conditions. It has to be stressed that bias in RNA-specific variation is not only restricted to inter-sample comparison, but also effects the relative ranking of intra-sample RNA expression.

In addition to the more obvious quantitative aspects of bias, there are also qualitative aspects because RNA-seq methods are also employed for RNA discovery. Individual RNAs might ‘drop out’ from detection because they are not amenable to cDNA construction under chosen conditions. Certainly, to increase the likelihood of productive RNA end-modification and increase the complexity of sRNA-seq library content, the parallel application of different strategies and reagents is required ( 33, 39, 63, 79). In particular, increased variability of adapter sequences (i.e. adapter pools) helps to increase the diversity of RNAs accessed ( 63). A recent report emphasized the need for parallel application of different RNA- and DNA-based adapter oligonucleotides ( 39).

The source of expression level distortion and the failure to detect certain RNA species at all, appear not to be limited to reactions of RNA/cDNA-end modification, but also apply to cDNA amplification by PCR ( 30–32). Differences in reaction buffers, RNA G/C-content, secondary structures, RNA length and primers might each lead to bias during PCR. Therefore, the analysis of sRNA-seq data to examine relative changes in gene expression or identification reflects both true RNA abundance and biases related to the methods applied.

The sRNA-seq approach is also attractive to medical research. Deep sequencing might permit the screening of multiple patient samples to monitor the progression of disease or treatment. Therefore, it is urgently necessary to establish rigid international standardizations to minimize distortion and to avoid misleading conclusions during RNA-seq data interpretation. In summary, qualitative and quantitative RNomics remain more challenging than anticipated.


From RNA to sequencing reads

The main goal of RNA-seq in most contexts is the accurate quantification of the original RNAs’ abundances in a sample, whether that refers to ‘bulk’ RNA from a homogenized cell population, or single cells. In practice, this amounts to correctly interpreting the number of sequencing reads that are obtained for each transcript. This problem is non-trivial due to several confounding factors preventing precise quantification, most of which are owed to the complexity of RNA-seq sample preparation.

Several steps are necessary to convert the RNAs in cell lysates into sequencing reads. Common to the vast majority of protocols are selecting which RNA is to be sequenced, the cDNA production steps of reverse transcription (RT) (often referred to as first-strand synthesis) and second-strand synthesis. The reason for selecting RNA to be sequenced is that the vast majority of RNA in cell lysates is ribosomal RNA, which is normally undesired. Removing it allows for more reads to be used towards the detection of less abundant RNA species of interest, such as mRNA. This is achieved by removing rRNA (‘ribodepletion’) or positive selection of RNAs of interest. RNA is replaced with DNA because RNA is problematic to work with it is subject to degradation through RNases and metal ion catalyzed hydrolysis at higher temperatures. It has a propensity to form secondary structures and cannot easily be amplified due to a lack of suitable enzymes and its compromised stability during thermal cycling. Synthesis of the second cDNA strand is necessary to enable adapter ligation for next generation sequencing, unless special adaptations are used [ 8]. Other protocols use the RT step to add adapter sequences directly, for instance by using a RT primer with overhanging adapter sequences. This idea is taken further in scRNA-seq protocols where the RT primer is often oligo-(dT)s (to capture polyadenylated mRNAs) with an overhang including adapter sequences, cell barcodes and unique molecular identifiers (see section UMIs e.g. 10x Chromium [ 9], Drop-seq [ 10] or InDrop [ 11]).

The only RNA-seq strategy that avoids cDNA conversion is direct RNA sequencing, as implemented by the ill-fated Helicos sequencing machine [ 12] or nanopore sequencing [ 13]. The latter is promising as a future technology producing long reads for single molecules it records the base sequence of individual nucleic acid strands as they are electrophoretically pulled through channels in a membrane. The system is plagued with high error rates, though, and most studies have exploratory character and/or use additional second generation (e.g. Illumina) sequencing to bolster sequencing quality [ 14].

RNA-seq libraries are usually fragmented by various means and size-selected in order to produce more sequencing reads at optimal length. This can occur before or after cDNA production. Direct fragmentation of RNA often uses metal-ion catalyzed hydrolysis at high temperatures (e.g. TruSeq) and cDNA fragmentation often uses physical methods (e.g. sonication) or enzymatic methods. ‘Tagmentation’ is a convenient enzymatic way to combine fragmentation and adapter ligation [ 15]. It uses transposase Tn5 to internally cleave double-stranded DNA and ligate oligonucleotides to both resulting ends in the same reaction. The material is usually further amplified by PCR. Often, an extended first PCR cycle is used to synthesize the second-strand. An alternative to PCR is linear amplification by in vitro transcription (IVT), as implemented by the CEL-seq protocol [ 16].

Each of these steps can skew the representation of original transcripts by sequencing reads. It is worth noting that there is a difference between variability and bias. Statistically, the average of a repeatedly sampled value needs to deviate from the true value to make it an actual bias random variation en soi is not enough. Biases in RNA-seq can have very different effects and it is important to understand, classify and quantify these. Two key properties that help categorize biases are their escalader (local – bias is specific to one gene or individual positions, or global – bias occurs across genes in a systematic overall pattern) and their visibility (can be seen on a coverage plot, e.g. Figure 1A), which are explained in more detail below. These properties are not always independent.

(UNE) RNA-seq coverages by sequencing read along an example gene (Ube2s) for two biological replicates. Abrupt changes in exonic read densities (vertical dashed lines) often coincide across samples, suggesting that the local sequence environment is responsible for this type of bias. Data from GEO, accession numbers GSM710183 and GSM710184. (B) RNA-seq coverage along a typical transcript can be subject to bias at different scales the schematic illustration depicts an absence of visible bias (top left), a local bias (bottom left), a global bias (top right) and a combination of the latter two (bottom right). (C) Global bias depends on transcript lengths. Schematic illustration of the length-dependent effects compared to a short reference transcript with no visible bias (top left). Upon considering longer transcripts in the same sample, a global bias can appear (bottom left), which does not necessarily lead to a skewed overall representation of the transcripts (the dashed horizontal line indicates average coverage equal to the reference). However, different lengths often do lead to unequal representation of transcripts due to global bias that might be invisible or visible in terms of coverage (top and bottom right, respectively).

(UNE) RNA-seq coverages by sequencing read along an example gene (Ube2s) for two biological replicates. Abrupt changes in exonic read densities (vertical dashed lines) often coincide across samples, suggesting that the local sequence environment is responsible for this type of bias. Data from GEO, accession numbers GSM710183 and GSM710184. (B) RNA-seq coverage along a typical transcript can be subject to bias at different scales the schematic illustration depicts an absence of visible bias (top left), a local bias (bottom left), a global bias (top right) and a combination of the latter two (bottom right). (C) Global bias depends on transcript lengths. Schematic illustration of the length-dependent effects compared to a short reference transcript with no visible bias (top left). Upon considering longer transcripts in the same sample, a global bias can appear (bottom left), which does not necessarily lead to a skewed overall representation of the transcripts (the dashed horizontal line indicates average coverage equal to the reference). However, different lengths often do lead to unequal representation of transcripts due to global bias that might be invisible or visible in terms of coverage (top and bottom right, respectively).

In the next section, we introduce the two major methods for quantifying the abundance of RNA in a sample. We discuss how the sample preparation process introduces bias for coverage-based approaches, avoids these biases for UMI-based approaches, and how these approaches compare otherwise.


1. Introduction—the Importance of Transcriptional Profiling at Single Cell Resolution

Transcriptome profiling has been a popular tool in molecular biological research for more than a decade. Mostly implemented by microarray technology, it has led to numerous insights and discoveries. These range from cell type specific single factors that were identified in differential expression screens to findings based on large portions of the entire transcriptome, such as disease signatures (e.g., [1]) and the interrelations of epigenetic modifications and gene expression (e.g., [2]).

Following its recent introduction, RNA-sequencing (RNA-seq) [3,4] is rapidly replacing microarrays as the method of choice for the aforementioned endeavors. Besides superior accuracy in the quantification of expression, RNA-seq offers other advantages, such as the possibility to detect novel transcripts, splice variants or allele-specific expression [5]. mRNA-profiling was performed on single cells early [6] and RNA-seq is following suit [7].

The analysis of the transcriptome in individual cells offers a number of advantages compared to cell-averaging experiments. Tissues or other cellular assemblages are heterogeneous even if a single, “traditional” cell type is concerned. This is particularly apparent from the study of immune cell types that are often defined based on the expression of surface markers. Improving experimental technologies reveal continuous expression ranges and fairly unrestricted combinatorial expression of surface markers [8]. This means that boundaries between cell types are blurred and that every individual cell is different. A similar picture emerges with tumors. While tumors have been known to be heterogeneous mixtures of cell types for a long time [9], pioneering studies demonstrate the potential of genome sequencing in individual cells [10]. RNA-seq will thus provide a powerful means to facilitate functional characterization of the differences among the cells in a tumor.

While tumor cell heterogeneity can be attributed to accumulating mutations, even genetically identical cells, under identical conditions, display high variability in their gene and protein expression levels. This is usually referred to as “noise”, defined as variance or standard deviation over mean [11]. A number of studies have probed into the origins and mechanisms of noise and found it to be mostly due to the stochastic effects associated with the low numbers of involved molecules [12].

While standard microarray or RNA-seq experiments yield mean expression levels, distributions in single cells demonstrate that only a negligible portion of cells express mRNAs close to the actual mean levels. Depending on the skewness of the distributions, this could mean that cells that express certain mRNAs at “outlier” levels are functionally important, yet remain undetected by traditional experiments. Knowledge about the shape of the distributions can also be used to understand mechanisms that are involved in transcriptional regulation [13].


  • Edwin Antony, Ph.D.
  • Yuna Ayala, Ph.D.
  • Angel Baldan, Ph.D.
  • Tomasz Heyduk, Ph.D.
  • Michelle Pherson, Ph.D.
  • Nicola Pozzi, Ph.D.
  • Tracey Baird
  • Enrico Di Cera, M.D.
  • Joel Eissenberg, Ph.D.
  • Zachary Montague
  • Angela Spencer

Abdul Waheed, Ph.D., Emeritus Research Professor of Biochemistry, gifted $1 million to the Department of&hellip

Congratulations to Kaush Amunugama, Ph.D. Student in Dave Ford's lab, on receiving an ASBMB 2021 Graduate/Postdoctoral&hellip

Congratulations to David Ford, Ph.D., Professor of Biochemistry and Director of the Center for Cardiovascular&hellip

Congratulations to Sergey Korolev, Ph.D., whose application was the first to be funded through the&hellip

A recently published commentary by David Ford, Professor of Biochemistry, was highlighted in ASBMBToday. The&hellip


Voir la vidéo: ECL8202, Cours Étapes de développement dun modèle hiérarchique bayésien (Janvier 2022).