Informations

Comptes de gènes normalisés à edgeR/DESeq2 ?


Les deux packages d'expression différentielle très populaires supposent un nombre brut de matrice de gènes. Cela est logique car le modèle statistique modélise la profondeur de la bibliothèque. Mais que se passe-t-il si je normalise mes comptes avec, par exemple, la norme ERCC enrichie en ARN ? Cela affecterait-il la puissance statistique? Faux positif? Pourquoi exactement n'est-il pas recommandé de donner une matrice de comptage normalisée ?


Le conseil est de ne pas utiliser du tout les spikes ERCC en raison des variations introduites par le pipetage aux volumes qu'ils recommandent.

Le fil explique également comment utiliser DESeq et EdgeR avec la normalisation des pics, le processus étant considérablement plus facile avec DESeq, où vous pouvez utiliser les calcSizeFactors sur une matrice de nombre de lectures de pics uniquement. Avec edgeR, vous devrez passer des valeurs à l'aide du paramètre lib.sizes dans les fonctions appropriées.

Si vous souhaitez utiliser des comptes fractionnaires via limma-voom, cela devrait fonctionner, je pense; J'ai eu de bons résultats en utilisant voom sur les comptes RSEM, qui sont fractionnaires.


Analyse des données RNA-Seq dans R - Enquêtez sur les gènes exprimés de manière différentielle dans vos données !

Dans ce didacticiel, un binôme négatif a été utilisé pour effectuer une analyse différentielle de l'expression génique dans R à l'aide des packages DESeq2, pheatmap et tidyverse. Le flux de travail pour les données RNA-Seq est :

  • Obtention des fichiers de séquençage FASTQ à partir de l'installation de séquençage
  • Évaluer la qualité des lectures de séquençage
  • Effectuer l'alignement du génome pour identifier l'origine des lectures
  • Générez la matrice de comptage des lectures alignées, c'est-à-dire le nombre de lectures alignées sur les exons de chaque gène.

Résumé

Fond

Plusieurs packages R existent pour la détection de gènes différentiellement exprimés à partir de données RNA-Seq. Le processus d'analyse comprend trois étapes principales, à savoir la normalisation, l'estimation de la dispersion et le test d'expression différentielle. Les étapes de contrôle de la qualité tout au long de ce processus sont recommandées mais pas obligatoires, et le fait de ne pas vérifier les caractéristiques de l'ensemble de données peut conduire à des résultats erronés. De plus, les méthodes de normalisation et les modèles statistiques ne sont pas échangeables entre les packages sans des transformations adéquates dont les utilisateurs ne sont souvent pas conscients. Ainsi, des pipelines d'analyse dédiés sont nécessaires pour inclure des étapes de contrôle qualité systématiques et éviter les erreurs dues à une mauvaise utilisation des méthodes proposées.

Résultats

SARTools est un pipeline R pour l'analyse différentielle des données de comptage d'ARN-Seq. Il peut gérer des conceptions impliquant deux ou plusieurs conditions d'un seul facteur biologique avec ou sans facteur de blocage (comme un effet de lot ou un appariement d'échantillons). Il est basé sur DESeq2 et edgeR et est composé d'un package R et de deux modèles de script R (respectivement pour DESeq2 et edgeR). En réglant un petit nombre de paramètres et en exécutant l'un des scripts R, les utilisateurs ont accès aux résultats complets de l'analyse, y compris des listes de gènes différentiellement exprimés et un rapport HTML qui (i) affiche des tracés de diagnostic pour le contrôle qualité et la vérification des hypothèses de modèle et (ii) assure le suivi de l'ensemble du processus d'analyse, des valeurs des paramètres et des versions des packages R utilisés.

Conclusion

SARTools fournit des contrôles de qualité systématiques de l'ensemble de données ainsi que des tracés de diagnostic qui aident à ajuster les paramètres du modèle. Il donne accès aux principaux paramètres de DESeq2 et edgeR et empêche les utilisateurs non formés d'abuser de certaines fonctionnalités des deux packages. En gardant une trace de tous les paramètres du processus d'analyse, il répond aux exigences de la recherche reproductible.

Citation: Varet H, Brillet-Guéguen L, Coppée J-Y, Dillies M-A (2016) SARTools: A DESeq2- and EdgeR-Based R Pipeline for Comprehensive Differential Analysis of RNA-Seq Data. PLoS ONE 11(6) : e0157022. https://doi.org/10.1371/journal.pone.0157022

Éditeur: Ken Mills, Queen's University Belfast, ROYAUME-UNI

A reçu: 6 avril 2016 Accepté: 23 mai 2016 Publié : 9 juin 2016

Droits d'auteur: © 2016 Varet et al. Il s'agit d'un article en libre accès distribué selon les termes de la licence d'attribution Creative Commons, qui permet une utilisation, une distribution et une reproduction sans restriction sur tout support, à condition que l'auteur et la source d'origine soient crédités.

Disponibilité des données: Toutes les données pertinentes sont disponibles sur https://github.com/PF2-pasteur-fr/SARToolsPaperData.

Le financement: Ce travail a été soutenu par l'infrastructure France Génomique Nationale, financée dans le cadre du programme « Investissements d'Avenir » géré par l'Agence Nationale pour la Recherche (contrat ANR-10-INBS-09).

Intérêts concurrents : Les auteurs ont déclaré qu'ils n'existaient pas de conflit d'intérêts.


Comptage vs FPKM dans RNA-seq

La plupart du temps, la raison pour laquelle les gens effectuent l'ARN-seq est de quantifier les niveaux d'expression des gènes. En théorie, l'ARN-seq est une donnée au niveau du rapport, et vous devriez être légitimement en mesure de comparer le gène A dans l'échantillon 1 contre l'échantillon 2 ainsi que le gène A contre le gène B dans l'échantillon 1.

Il existe deux manières principales de mesurer l'expression d'un gène, ou d'un transcrit, ou autre, dans les données RNA-seq :

  1. compte sont simplement le nombre de lectures chevauchant une caractéristique donnée telle qu'un gène.
  2. FPKM ou Fragments par kilobase d'exons par million de lectures sont beaucoup plus compliqués. Fragment signifie fragment d'ADN, de sorte que les deux lectures qui constituent une lecture appariée comptent pour une. Par kilobase d'exon signifie que le nombre de fragments est ensuite normalisé en divisant par la longueur totale de tous les exons du gène (ou du transcrit). Ce peu de magie permet de comparer le gène A au gène B même s'ils sont de longueurs différentes. Par million de lectures signifie que cette valeur est ensuite normalisée par rapport à la taille de la bibliothèque. Ce peu de magie permet de comparer le gène A de l'échantillon 1 à l'échantillon 2 même si la bibliothèque d'ARN-seq de l'échantillon 1 contient 60 millions de paires de lectures et la bibliothèque de l'échantillon 2 n'a que 30 millions de paires de lectures.

(En fait, comme cet article le montrera, il y a plus de différences entre les deux méthodes que celles-ci - j'y reviendrai dans la conclusion.)

À mon avis, la normalisation par la longueur exonique et la taille de la bibliothèque semble être une évidence, j'utilise donc des FPKM et je n'avais jamais compris pourquoi quelqu'un utiliserait des comptes. Mais si vous voulez vraiment défendre votre analyse, vous devez pouvoir répondre à n'importe quelle question avec “Oui, j'ai essayé et voici ce que j'ai trouvé,” et donc je voulais répéter mon analyse en utilisant des comptes. Pendant ce temps, un collègue qui s'intéresse aux décomptes m'a dit que les FPKM appliquaient trop de normalisation, faisant disparaître une partie de la différence entre un échantillon et un autre. Pourquoi serait-ce le cas ? J'ai décidé tant que j'allais répéter mes analyses à l'aide de décomptes, autant faire une comparaison côte à côte avec les FPKM afin de vraiment comprendre en quoi le comportement diffère.

Pour comparer les deux, je me suis tourné vers mon jeu de données RNA-seq : Human BodyMap 2.0. Pour les besoins de cet exercice, je n'examinerai que les transcriptions connues.

comment calculer les FPKM

comment calculer les comptes

Vous pouvez calculer le nombre à l'aide de bedtools multicov, mais vous avez besoin d'un fichier d'annotation de transcription au format BED pour indiquer à bedtools où chercher – contrairement aux boutons de manchette avec le paramètre -N 1, multicov ne va pas découvrir de nouvelles transcriptions pour vous. Pour rendre les comptes directement comparables aux FPKM que j'ai calculés plus tôt, je voulais utiliser ce même fichier d'annotation de transcription et le convertir du format GTF au format BED.

Dès le départ, les choses se compliquent. J'ai remarqué que le fichier d'annotation de transcription d'origine comporte une ligne pour chaque combinaison d'une transcription avec un exon ou une séquence codante ou un codon de démarrage ou d'arrêt. Considérons PRNP, qui n'a que deux exons (l'exon 1 est l'UTR 5 & 8242 et l'exon 2 est la séquence codante et 3 & 8242UTR) et vraiment un seul transcrit - il n'y a pas de variations d'épissage majeures dont je suis au courant. Il a 18 lignes dans ce fichier.

C'est parce que 4 versions distinctes de PRNP ont en quelque sorte été intégrées à la base de données de l'Ensemble en tant que transcriptions distinctes.

Dans tous les cas, si c'est si mauvais pour PRNP, vous pouvez imaginer combien de lignes sont présentes pour les gènes qui ont légitimement beaucoup de variantes d'épissage :

Cela nous pose un problème. Maintenant, si nous voulait compte pour chaque exon possible, nous pourrions simplement utiliser l'outil gtf2bed dans bedops qui convertira ce fichier GTF d'origine en un fichier BED, ligne par ligne :

5 minutes). Mais il est plus probable que notre unité d'analyse soit les transcrits ou les symboles génétiques. Si nous devions faire des comptes par exon, puis regrouper par transcrit ou symbole de gène et prendre la somme des comptes d'exons, nous comptions quatre fois chaque exon dans PRNP et chaque exon dans TTN beaucoup plus de fois que cela ! Ce dont nous avons besoin, c'est de convertir le fichier GTF en une ligne par, disons, symbole de gène, si le symbole de gène est notre unité d'analyse.

Il s'avère que les ea-utils d'Erik Aronesty contiennent un script Perl pour faire exactement cela. Il s'appelle gtf2bed tout comme l'outil bedops ci-dessus, donc pour plus de clarté, je l'ai renommé gtf2bed_2.pl . Pour le télécharger et l'exécuter :

1 minute. (Avertissement : si vous utilisez cet article comme pipeline, notez que l'utilisation du fichier BED résultant sans modification peut donner des résultats très absurdes pour les gènes répertoriés dans plusieurs loci – voir la discussion sur SNORD60 plus loin dans cet article).

Si vous ouvrez le fichier lit résultant, vous verrez que les trois premières colonnes sont simplement le chromosome, le site de début de transcription (le plus ancien) et le site de fin de transcription (le plus récent) pour le gène – ou en d'autres termes, l'union de tous les éléments transcrits. sites dans ce gène sur tous les transcrits possibles.

Maintenant, gtf2bed_2.pl observe un format BED12 très correct et donc Est-ce que préservez les informations sur la structure de l'exon sous la forme des colonnes blockSizes et blockStarts. Mais multicov ne lit que les trois premières colonnes. Par conséquent, lorsque vous comptez sur ce fichier que nous venons de créer, vous comptez à la fois les introns et les exons. Pour autant que je sache en interrogeant autour de moi, c'est ainsi que tout le monde fait son RNA-seq count.

Comparez cela avec les FPKM, où les boutons de manchette ne compteront que les lectures exoniques et se normaliseront par une longueur totale d'exon d'un gène (ou transcription) si vous comptez (au moins selon ce pipeline / à moins que vous ne fassiez d'autres trucs plus sophistiqués) vous incluez des lectures introniques. Ainsi, les comptes, contrairement aux FPKM, seront affectés par la quantité de contamination pré-ARNm (et donc la couverture intronique) que vous avez dans vos bibliothèques.

Cela dit, j'ai ensuite lancé multicov, comme ceci :

Ce qui a pris environ 50 heures de temps CPU.

Soit dit en passant, générer la liste des fichiers BAM pour cette commande est ennuyeux cette fois je l'ai fait avec echo -n :

Le fichier résultant contiendra le fichier de lit original à 12 colonnes créé par gtf2bed_2.pl plus, dans ce cas, 16 colonnes supplémentaires pour chacun des 16 BAM que j'ai appelés compte sur l'utilisation de multicov.

quelques covariables

Étant donné que les FPKM ne sont, en théorie, que des nombres normalisés par la taille de la bibliothèque et la taille de la transcription, j'ai pensé que je devrais également avoir ces deux valeurs à portée de main pour cette analyse. J'ai calculé la taille de la bibliothèque en nombre de lectures dans chaque BAM avec samtools view -c :

Ce qui prend un temps étonnamment long (

30 min/BAM), d'où la nécessité de soumettre chacun comme un travail.

L'autre covariable que je voulais était la longueur de chaque gène. Mais lequel longueur, demandez-vous? En utilisant le fichier BED que je viens de créer, il est facile d'obtenir la longueur du premier site de début de transcription au dernier site de fin de transcription possible :

Si vous voulez la longueur exonique, c'est légèrement plus délicat. De toute évidence, Cufflinks connaît ces informations sous une forme quelconque puisqu'elles sont utilisées pour la normalisation. Vous pouvez le retirer comme ceci :

Mais genes.fpkm_tracking n'a pas cela pour les gènes, probablement parce que choisir une longueur comme "la longueur pour un gène avec plusieurs transcrits est maladroit. Avec des scripts plus sophistiqués et une fusion de bedtools, vous pouvez obtenir la longueur de l'union de tous les exons possibles dans un gène, une sorte d'analogue au fichier gene.lengths.txt que nous venons de créer, qui est la longueur de l'union de toutes les transcriptions possibles. Mais je n'en aurai pas nécessairement besoin pour aujourd'hui.

Pour les besoins de l'argument, j'ai également calculé les longueurs moyennes directes pour chaque symbole de gène, aussi grossier soit-il. J'ai d'abord saisi les symboles et la longueur des gènes en bash :

et puis juste du SQL enveloppé dans R :

Avec les dénombrements, les FPKM et les covariables en main, j'ai cherché à comprendre comment et pourquoi ces mesures différaient les unes des autres.

Tout d'abord, les trucs de configuration ennuyeux:

Question la plus fondamentale : les décomptes et les FPKM sont-ils corrélés ? Je l'espère bien sûr ! Nous pouvons poser cette question de plusieurs manières. Tout d'abord, posons cette question à toutes les combinaisons de symboles de gènes et de tissus.

Cela est étrange. Dans l'espace linéaire (corrélation de Pearson’s), les comptes et les FPKM sont significativement mais à peine corrélés, avec rho = .006. Dans l'espace des rangs (corrélation de Spearman), ils sont assez fortement corrélés, rho = 0,81. Que pourraient bien ces données voir Comme?

C'est tellement extrême : de ce point de vue, il semble y avoir fondamentalement deux types de gènes : ceux avec quelques comptes mais

0 FPKM, et ceux avec quelques FPKM mais

0 compte. Incroyable que nous ayons vu une quelconque corrélation.

Cela est même vrai si nous prenons la valeur moyenne pour chaque gène à travers les multiples tissus considérés ici :

Les deux valeurs aberrantes les plus extrêmes étaient IGHJ6 et SNORD60, alors je les ai recherchées individuellement.

IGHJ6 ne mesure que 61 pb, à chr14 : 106 329 408-106 329 468, il n’est donc pas étonnant qu’il puisse avoir des comptes faibles mais des FPKM élevés. SNORD60, en revanche, est aussi un gène court, un snoRNA de seulement 83 pb à chr16:2 205 024-2 205 106. Alors, qu'est-ce que l'accord SNORD60′?

J'ai d'abord regardé les données brutes :

13 à 21 millions de lectures mais zéro FPKM dans de nombreux tissus. Il n'a pas fallu longtemps pour trouver la source du problème : dans le fichier BED que j'ai utilisé pour créer les comptes, SNORD60 fait 204 Mo de long :

Ce qui s'avère être parce que dans le fichier GTF d'origine, il est répertorié avec trois exons dans des loci génomiques complètement différents.

Ainsi, lorsque j'ai exécuté gtf2bed_2.pl pour convertir ce GTF en fichier BED, il a simplement choisi la base de départ la plus basse et la base de fin la plus élevée comme points de terminaison d'une transcription.

Il s'est avéré étonnamment difficile de trouver un moyen de filtrer de tels cas. L'histogramme des longueurs de gènes dans mon fichier BED est tout aussi extrême que les tracés précédents :

À la recherche d'une coupure pour filtrer les gènes dont la longueur est évidemment une erreur, j'ai recherché sur Google le « gène humain le plus long » et j'ai trouvé le DMD, qui mesure près de 2,3 Mo. L'histogramme des gènes ≤ 2,3Mb est légèrement meilleur que le premier histogramme :

C'est plus proche de la distribution exponentielle à laquelle je m'attendrais, bien que je soupçonne qu'il y a encore des gènes faussement longs dans cette distribution aussi.

Si ce sous-ensemble, de gènes < 2,3 Mo, est plus rationnel et a au moins éliminé certaines des erreurs les plus scandaleuses, j'aurais espéré qu'il serait possible d'expliquer une grande partie de la variabilité des nombres par rapport aux FPKM au sein de ce sous-ensemble :

Mais non, modèle linéaire des FPKM

counts donne un R^2 de seulement 0,008. L'inclusion de la longueur des gènes dans le modèle n'a pas aidé :

Et diviser explicitement les comptes par la longueur des gènes n'a que peu aidé, nous amenant à un R^2 de 0,016 :

Cet ensemble de données comprend 52 686 symboles de gènes Ensembl, alors je me suis demandé si les données seraient peut-être mieux comptées si nous ne considérions que les 23 705 gènes hg19 RefSeq. Cela n'a aidé qu'un peu, nous amenant à un R^2 de .026 :

Et quand je suis revenu à toutes les combinaisons gène-tissu avec cet ensemble de données plus limité, j'ai finalement obtenu un rho de 0,26 pour une corrélation de Pearson et de 0,83 pour Spearman.

Ce n'est toujours pas une corrélation aussi étroite que je l'avais espéré, étant donné que ces deux mesures sont censées mesurer globalement la même chose "l'expression des gènes" dans le même ensemble de données. À titre de comparaison, lorsque j'exécute mon pipeline QC d'expression génique standard sur des données RNA-seq pour différents échantillons mais que j'appelle en utilisant le même pipeline, je trouve souvent une corrélation de Pearson entre les échantillons de 0,85 ou mieux. Alors qu'ici, pour les mêmes données appelées avec deux pipelines différents, j'obtiens un Pearson’s de seulement 0,26. C'est peut-être un autre rappel malheureux de l'irreproductibilité des découvertes d'expression génique. Les technologies utilisées (y compris les différents pipelines bioinformatiques) introduisent plus de variabilité que ce qui est présent dans les échantillons sous-jacents eux-mêmes.

J'ai pensé qu'une explication possible pourrait être la différence entre la longueur exonique et la longueur totale du gène. Ici, les comptes sont évalués sur la longueur totale des gènes, et je les ai ensuite divisés par la longueur totale des gènes, tandis que les FPKM sont évalués sur les exons et normalisés par la longueur exonique. Au sein de cet ensemble de gènes relativement bien comportés ≤ 2,3Mb et dans RefSeq, la corrélation entre la longueur totale et la longueur exonique n'est encore que de 0,19 dans l'espace linéaire et de 0,49 dans l'espace des rangs :

Ce qui suggère qu'au moins une partie du problème ici est simplement que les comptes, qui incluent les exons et les introns, mesurent quelque chose de très différent des FPKM, qui n'incluent que les exons.

Il semble donc que ces deux mesures ne font que mesurer quelque chose de différent et obtiennent des réponses différentes (comme en témoigne la faible corrélation entre elles). Cela suggère qu'au plus l'une des deux méthodes "comptes et FPKM" convient pour comparer le gène A au gène B. Au moins à un niveau de ratio, c'est-à-dire. On peut soutenir que, étant donné que la corrélation de Spearman est plus forte, les deux pourraient convenir pour les analyses de niveau ordinal.

Cela ne fait que comparer le gène A au gène B. Mais souvent, la réponse que nous recherchons dans nos analyses est de trouver des gènes dont le niveau d'expression est en corrélation avec une variable d'intérêt, par exemple un génotype, un traitement médicamenteux ou un point dans le temps. De tels résultats ne seront reproductibles entre les comptages et les FPKM que dans la mesure où les comptages et les FPKM pour chaque gène individuel sont corrélés entre les échantillons. Dans ce cas, nos “échantillons” sont les 16 tissus différents de Human BodyMap 2.0. Pour évaluer dans quelle mesure le niveau de chaque gène est reproductible dans différents tissus, j'ai fait un « tracé du volcan », d'abord, des corrélations de Pearson :

Les résultats sont bien meilleurs que ce à quoi je m'attendais :

La corrélation de Pearson % de gènes
positif (p < .05) 83%
aucun (p > .05) 6%
négatif (p < .05) 0.01%
N / A* 11%

*Les valeurs NA résultent des rangées où tous les tissus avaient 0 comptage ou tous avaient 0 FPKM, donc le test de corrélation a échoué.

Étonnamment, lorsque j'ai réexécuté cela avec Spearman’s, les résultats étaient pratiquement identiques (tous les chiffres du tableau ci-dessus étaient à une fraction de pour cent).

Ainsi, pour la plupart des gènes, la différence entre les différents niveaux d'expression d'échantillons de ce gène est au moins nominalement reproductible entre les deux métriques considérées ici : les nombres et les FPKM. Cependant, j'hésite à attribuer trop d'importance à cette conclusion car ce que j'utilise ici comme exemple de jeu de données est l'expression à travers différents tissus, par opposition à différents personnes. Les différences d'expression génique entre les tissus sont assez importantes et assez fondamentales pour la biologie, et je m'attendrais à ce que les différences entre les individus soient beaucoup plus subtiles. Si les mêmes différences interindividuelles apparaissent dans les décomptes que dans les FPKM, je ne peux pas le dire dans cet exemple.

conclusions

Le nom “FPKM” – fragments par kilobase d'exon par million de lectures – implique que FPKM est une mesure de l'expression génique normalisée par la longueur exonique et la taille de la bibliothèque, contrairement aux comptes bruts. Cependant, au cours de cet exemple, j'ai réalisé qu'il existe plusieurs autres différences entre les comptes et les FPKM :

  • Lorsqu'une lecture chevauche plusieurs définitions d'exons ou plusieurs définitions de transcription, Cufflinks décide à quelle(s) transcription(s) attribuer la lecture lorsqu'il calcule les FPKM. Le calcul des comptes, du moins dans le pipeline simple que j'ai présenté ici, n'est pas aussi sophistiqué.
  • En conséquence, les comptes ne sont normalement évalués que par le symbole du gène. S'ils étaient évalués par transcription, de nombreuses lectures seraient comptées deux fois (ou même comptées des dizaines de fois) car de nombreux gènes ont une multiplicité de transcriptions. En comparaison, il existe relativement peu de loci génomiques où deux gènes distincts se chevauchent.
  • Les FPKM ne comptent que les alignements exoniques, les comptes (au moins ce pipeline) incluent les introns. La longueur totale d'un gène (y compris les introns) n'est que modestement corrélée à sa longueur exonique (rho = 0,19), ce qui fait une grande différence.
  • Les pipelines générant des comptages ne sont généralement pas capables de découvrir des transcriptions. Au lieu de cela, vous devez leur fournir une liste de loci génomiques avec des gènes connus (avec les FPKM, c'est facultatif). Il est important de faire attention à ce que la fusion des transcrits en une ligne par gène ne crée pas de résultats absurdes comme nous l'avons vu pour SNORD60 ci-dessus.

Toutes ces différences semblent contribuer à expliquer pourquoi les FPKM et les décomptes que j'ai appelés ici "sur exactement le même ensemble de données" ont si peu de corrélation les uns avec les autres (R^2 < .01 même après avoir supprimé les valeurs aberrantes de la longueur des gènes ). Malgré cela, les FPKM et les comptes pour n'importe quel gène peut être un peu plus reproductible, bien que cette analyse ait considéré différents tissus (qui ont d'énormes différences dans l'expression des gènes) et non des individus différents (qui ont des différences subtiles dans l'expression des gènes).

Étant donné que les décomptes et les FPKM semblent mesurer des choses assez différentes, il est à débattre de la mesure la plus valide. Je vais me lancer et argumenter un peu pour les FPKM. Les bibliothèques d'ARNm-seq sont enrichies pour les ARNm, généralement par sélection polyA, éliminant ainsi, espérons-le, la plupart des couvertures introniques. Étant donné que vous utilisez une méthode de laboratoire spécifiquement pour obtenir uniquement des ARNm, votre pipeline doit correspondre à cela et ne compter que les exons. De toute évidence, les FPKM représentent également une méthode plus sophistiquée impliquant l'attribution de lectures à des transcriptions particulières et la normalisation de la longueur exonique et de la taille de la bibliothèque, toutes bonnes choses. Je n'ai entendu personne nier cela, l'argument que j'ai entendu pour les décomptes était qu'il s'agissait d'une mesure différente qui peut avoir plus de variabilité et plus de puissance pour certaines choses. Mais rien de ce que j'ai vu ici ne m'a convaincu que cette variabilité supplémentaire reflète quelque chose de significatif que vous voudriez analyser.

Cela dit, ma motivation initiale pour ce post est que vous voulez toujours faire l'analyse dans les deux sens afin de pouvoir répondre à toutes les questions.

À propos d'Eric Vallabh Minikel

Eric Vallabh Minikel est dans une quête permanente pour prévenir la maladie à prions. Il est un scientifique basé au Broad Institute du MIT et à Harvard.


Matériaux et méthodes

Dans cette section, nous décrivons les méthodes de normalisation utilisées dans notre étude ainsi que les critères spécifiques utilisés dans notre comparaison. Nous discutons également de l'étude TCGA et de l'étude de simulation utilisée pour évaluer les méthodes. Nous avons examiné les effets de différentes méthodes de normalisation sur l'analyse d'expression différentielle à l'aide de trois workflows d'analyse (Fig. 1). Ici, les flux de travail 1 et 2 comparent l'effet des méthodes de normalisation de la taille de la bibliothèque sur l'analyse de l'expression différentielle, tandis que le flux de travail 3 compare différentes méthodes pour estimer les artefacts latents suivis par la normalisation de l'échantillon pour ces facteurs inconnus, et considère également l'impact de ne pas tenir compte de la perte de degrés de liberté dus à la normalisation pour l'analyse d'expression différentielle. Notez que les workflows 1 à 3 sont utilisés par l'étude cervicale TCGA, et pour l'étude de simulation, seul le workflow 3 est pris en compte.

Utilisation de l'ensemble de données CESC pour comparer l'effet des méthodes de normalisation de la taille de la bibliothèque sur l'analyse de l'expression différentielle à l'aide du workflow 1 (la matrice de conception contient le principal facteur d'intérêt) et du workflow 2 (la matrice de conception contient l'ID de lot ainsi que le principal facteur d'intérêt). Le workflow 3 compare les différentes méthodes pour estimer les artefacts latents suivis d'une normalisation de l'échantillon pour ces facteurs inconnus, et considère l'impact de la non-prise en compte de la perte de degrés de liberté due à la normalisation pour l'analyse de l'expression différentielle en utilisant à la fois l'ensemble de données CESC et les données simulées. Pour les données simulées et CESC, nous avons envisagé deux approches pour déterminer les gènes DE (c. comprend la principale variable d'intérêt. Notez que dans le flux de travail 3, la méthode de normalisation UQ n'a pas été prise en compte pour les données simulées, et il n'y a pas non plus d'artefact technique connu (par exemple, l'ID de lot).

Normalisation de la longueur des gènes

L'impact de la longueur des gènes sur l'estimation de l'abondance des gènes est un biais technique non observé avec les données de puces à ADN, mais observé dans l'achèvement des études RNA-Seq. En particulier, les gènes les plus gros auront inévitablement un nombre de lectures plus élevé que les gènes plus petits en raison de la différence de longueur ou de taille de leurs gènes [15]. Une méthode souvent utilisée pour corriger ce biais est l'utilisation de RPKM/FPKM (lectures/fragments par kilo-base par million de lectures mappées) [15,29,30]. Une autre approche pour ajuster la longueur du gène est la méthode TPM (transcriptions par million), qui prend en compte à la fois la longueur du gène et les corrections de la longueur de lecture du séquençage. 24]. ERPKM est une amélioration de RPKM qui remplace la longueur du gène par une longueur de lecture effective (c. Ces approches redimensionnent le nombre de gènes pour corriger les différences de longueur des gènes, comme illustré dans S1A Fig.

Toutes ces méthodes reposent sur des approches de normalisation basées sur les comptes totaux ou effectifs, et ont tendance à être peu performantes lorsque les échantillons ont des distributions de transcrits hétérogènes [12,31]. La mise à l'échelle des comptes par longueur de gène peut donner des estimations biaisées de l'expression différentielle et l'association clairement positive entre la longueur des gènes et les comptes n'est pas entièrement supprimée en appliquant la normalisation de la longueur des gènes [10,12,30,31]. Cependant, les valeurs TPM et RPKM/FPKM sont appropriées à utiliser si l'objectif est de comparer les niveaux d'expression entre les gènes (c'est-à-dire la comparaison entre les gènes), l'analyse de l'expression différentielle nécessite que les niveaux d'expression soient comparés entre les échantillons .

Normalisation de la taille de la bibliothèque

Une source de variation entre les échantillons est la différence de taille de bibliothèque, où la taille de bibliothèque est le nombre total de lectures générées pour un échantillon donné. La différence de taille de bibliothèque peut être due à de nombreux facteurs, y compris des différences dans le multiplexage des échantillons (allocation des échantillons aux voies dans une cellule d'écoulement) ou des différences globales dans les niveaux d'expression génique (S1B Fig). L'objectif de la normalisation de la taille des bibliothèques est de rendre les tailles des bibliothèques comparables en mettant à l'échelle le nombre de lectures brutes dans chaque échantillon par un seul facteur spécifique à l'échantillon reflétant sa taille de bibliothèque. Il existe trois méthodes couramment utilisées : le quartile supérieur (UQ), la moyenne tronquée des valeurs M (TMM) et l'expression logarithmique relative (RLE).

  • Quartile supérieur (QU): Selon cette méthode de normalisation, après avoir supprimé les gènes ayant des comptes de lecture nuls pour tous les échantillons, les comptes de gènes restants sont divisés par le quartile supérieur des comptes différents de zéro dans le calcul des facteurs de normalisation associés à leur échantillon et multipliés par le quartile supérieur moyen sur tous les échantillons de l'ensemble de données [12]. Cette méthode de normalisation est implémentée dans le EDASeq et bordR Paquets de bioconducteurs [32,33].
  • Moyenne tronquée des valeurs M (TMM): Cette méthode de normalisation est basée sur l'hypothèse que la plupart des gènes ne sont pas différentiellement exprimés (DE). Pour chaque échantillon, le facteur TMM est calculé tandis qu'un échantillon est considéré comme échantillon de référence et les autres comme échantillons de test. Pour chaque échantillon d'essai, après exclusion des gènes les plus exprimés et des gènes avec les plus grands rapports logarithmiques, le TMM est calculé comme la moyenne pondérée des rapports logarithmiques entre ce test et la référence. En raison de l'hypothèse de faible DE, le TMM devrait être proche de 1. Si ce n'est pas le cas, sa valeur fournit une estimation du facteur de correction qui doit être appliqué aux tailles des bibliothèques [21]. Cette méthode de normalisation est implémentée dans le bordR Package bioconducteur comme méthode de normalisation par défaut [33].
  • Expression du journal relatif (RLE): Similaire à la TMM, cette méthode de normalisation est basée sur l'hypothèse que la plupart des gènes ne sont pas DE. Pour un échantillon donné, le facteur d'échelle RLE est calculé comme la médiane du rapport, pour chaque gène, de ses comptes de lecture sur sa moyenne géométrique sur tous les échantillons. En supposant que la plupart des gènes ne sont pas des DE, la médiane du rapport pour un échantillon donné est utilisée comme facteur de correction de tous les comptes de lecture pour répondre à cette hypothèse [34]. Cette méthode de normalisation est incluse dans le DESeq et DESeq2 Paquets de bioconducteurs [34,35].

Normalisation à travers l'échantillon

Comme les méthodes de normalisation de la taille de la bibliothèque corrigent principalement la profondeur de séquençage et ne s'ajustent pas aux autres variations techniques, des méthodes de normalisation des échantillons ont été proposées pour corriger d'autres artefacts techniques afin d'améliorer la qualité des données et la capacité à détecter les gènes biologiquement pertinents. Cependant, de telles sources techniques de variation deviennent problématiques à traiter lorsqu'elles sont corrélées ou confondues avec un facteur biologique principal d'intérêt. Par conséquent, une bonne conception expérimentale est essentielle lors de la réalisation d'études HTS. Ainsi, en complétant la normalisation pour les artefacts techniques, nous supposons souvent que le principal facteur d'intérêt est indépendant de tous les artefacts [16,17,25]. De plus, de nombreuses sources possibles de variations techniques ne sont pas enregistrées ou inconnues du chercheur. Par conséquent, en plus de la normalisation des artefacts techniques connus, l'évaluation et l'ajustement des variables potentielles inconnues ou latentes sont également justifiés [18].

Artefact technique connu

Contrairement aux méthodes de modélisation plus complexes, l'approche impliquant l'ajustement direct pour les artefacts techniques connus dans le modèle statistique approprié (par exemple, les modèles de régression linéaire). À titre d'exemple, un modèle avec les principaux facteurs biologiques d'intérêt peut être inclus dans la matrice de conception avec les artefacts techniques connus, ce qui peut améliorer la précision des comparaisons d'intérêt [16,36]. De plus, certaines méthodes de normalisation proposées ne sont pas robustes aux valeurs aberrantes dans les petits échantillons, où une méthode empirique de Bayes flexible, appelée Combat, a été proposé pour fournir des ajustements d'effet de lot plus robustes avec de petites tailles d'échantillon [26,27]. Cependant, le Combat L'approche met moins l'accent sur le modèle biologique et réduit principalement la variation globale même sans spécifier de modèle biologique. De plus, le Combat L'approche utilise une méthode empirique de Bayes pour éviter de sur-corriger les artefacts techniques connus, ce qui est essentiel pour une utilisation avec de petits lots (ou taille d'échantillon).

Artefact technique inconnu

Récemment, des méthodes de normalisation ont été développées pour évaluer et supprimer les variations techniques inconnues en estimant les facteurs latents pour capturer ces sources de variation. Plusieurs de ces méthodes reposent sur la décomposition en valeur singulière (SVD) ou sur d'autres techniques d'analyse factorielle pour déterminer la variation indésirable directement à partir des données. Un problème avec l'application de ces approches est la difficulté de distinguer la variation technique indésirable des facteurs biologiques d'intérêt. Par conséquent, la méthode de suppression des variations indésirables (RUV) a été proposée pour ajuster les variations techniques inconnues en effectuant une analyse factorielle sur les gènes témoins négatifs (c'est-à-dire les gènes connus pour ne pas être liés au principal facteur d'intérêt) [17,37]. Therefore, variations in the expression levels of these genes can be assumed to be unwanted variations. Housekeeping genes [38] or spike-in controls [39] are the examples of negative controls. However, RUV method does not need the negative control genes or samples [17]. Other commonly used methods to address this problem in identifying the unknown technical variations are the surrogate variable analysis (SVA) [18] and principal component analysis (PCA) [40]. It should be noted that in the case of RUV, SVA and PCA methods, it is possible that some of the estimated latent factors are not technical artifacts but rather represent true biology presented in the data. Thus, it is important to adjust for any known biological factors of interest and known technical artifacts prior to estimation of latent factors. The correct usage of these methods in estimating the latent technical artifacts has the potential to increase statistical power in downstream differential expression analysis, while note that increasing the number of estimated batch effects also can reduce power due to the additional bias of degrees of freedom [27,41].

After estimating the latent factors using RUV, SVA and PCA approaches, an appropriate statistical approach (e.g., linear model or ComBat) is used to obtain the normalized data.

  • Remove Unwanted Variation (RUV): Under this approach, the factors of unknown technical variations are estimated and removed by performing the factor analysis on suitable sets of negative control genes or samples by keeping the primary factor of interest. Therefore, RUV [17] method is divided into three sub-methods: RUVg, RUVs, and RUVr. The RUVg and RUVs are used when negative control genes and negative control samples (i.e., samples whose read counts are not influenced by the primary factor of interest) exist. However, RUVr (i.e., residual RUV) does not require the existence of negative control genes or samples. SVD is then computed on the residual matrix to estimate the factors of unknown technical variations. The number of factors of unwanted variation, k, should be guided by considerations that include samples sizes, extent of technical effects captured by the first k factors, and extent of differential expression [17,25].
  • Surrogate Variable Analysis (SVA): In this approach, the unknown technical variations or “surrogate variables” (SV’s) are estimated by applying SVD on the computed residual matrix and selecting significant eigenvectors [18,25,42]. The first step in SVA is to determine the number of SVs using one of the two methods, “BE” or “Leek” as noted in [18,26,43]. The “BE” method is based on a permutation procedure originally proposed by Buja and Eyuboglu [43], while the “Leek” method provides an interface to the asymptotic approach proposed by Leek [42], where under the specific assumptions, the right singular vectors are asymptotically consistent for latent artifacts as the number of features grows large. Once the number of SVs is calculated, then using the two-step algorithm following Leek and Storey [18] to estimate unknown technical artifacts.
  • Principal Component Analysis (PCA): This approach is completed by applying SVD to the scaled residual matrix to estimate the factors of unknown technical variations [44]. One can determine the number of PCs to include in the model by multiple methods, including: PCs that explain a given percent of the variation PCs that are associated with the biological factors of interest (i.e., confounders) top PCs regardless of association with the primary factors of interest application of the Tracy-Widom test for determining eigenvalues significantly different from zero (noting that the assumption of independence is not valid) [45–47] or determine the PCs to include based on a permutation testing approach similar to that implemented in the SVA method.

Issues of loss of degrees of freedom

For practical purposes it is more convenient to perform downstream analyses on the batch adjusted or normalized data without further consideration of technical artifacts effects. However, adjustment for technical artifacts reduces the effective degrees of freedom in the dataset and thus changes the null distribution of the test statistics. Not accounting for this change in the degrees of freedom due to normalization for batch effects may lead to increase the false positive rates, especially when the primary factors of interest are not equally represented in all batches or batch effects act as a confounder [41]. It should be noted that whatever our normalization approach is, one is in essence reducing the degrees of freedom in the data which in turn should lower the statistical power of the test. For example, let’s assume there are two studies (Study A and Study B) with the same sample size. Let’s also assume that Study A implemented good experimental design and did not need to normalize for any known or unknown technical artifacts, while Study B did not implement adequate experimental design and thus needed extensive across sample normalization. In this situation, Study A will have more power to detect a true biological effect compared to Study B. However, the loss of degrees of freedom associated with Study B due to normalization is often overlooked in the implementation of analysis approach involving across sample normalization followed by association analysis of gene expression levels with the biological factor of interest. In the following sections, we assess the extent of the impact of the loss of degrees of freedom on the type I error rate on association testing via a simulation study.

Comparison of methods

TCGA cervical study.

To compare the normalization methods and their impact on the differential expression analysis, publically available data from the TCGA cervical study (CESC) was used [28]. Level 3 RNA-Seq data (summarized gene expression levels) and clinical patient data were downloaded via Genomic Data Commons (GDC) (https://gdc.cancer.gov/) (July 2017). The large-scale study the size of the TCGA unavoidably generated technical artifacts. These factors (e.g., tissue source site, plate ID, sequence center) were tracked for each sample with this information contained within the original TCGA ID. These known factors can also be downloaded from MBatch, a web-based analysis tool for normalization of TCGA data developed by MD Anderson. (https://bioinformatics.mdanderson.org/tcgabatcheffects).

For the CESC study, gene expression data was measured on 60,433 genes and 178 cervical tissue samples (144 squamous cell carcinomas, 31 adenocarcinomas and 3 adenosquamous cancers). The integrative clustering analysis reported in the main CESC TCGA paper used mRNA, DNA methylation, miRNA and copy number variation data identified two squamous-carcinoma-enriched groups and one adenocarcinomas-enriched group [28]. The two squamous-carcinoma groups differ largely based on gene expression levels where one squamous cluster had high expression of keratin gene family members (keratin-high) and the other squamous cluster had low expression of keratin genes (keratin-low). Hence, for comparison of the normalization methods and impact on the differential expression analysis results, we restricted our analysis to the squamous cell carcinomas and set out to determine DE genes between the keratin-high (N = 47) and keratin-low (N = 86) tumors groups. After filtering non-expressed or low-expressed genes based on counts per million (CPM), 20,884 genes with CPM values above 1 in at least 3 libraries remain.

Simulation study.

An extensive simulation study was completed to compare the performance of SVA (“BE” and “Leek”) and PCA (based on different percent of variation) methods to identify the number of latent factors (i.e., SVs), and determine SVs, where the residual RUV method was also included. Then, the performance of different across latent factor identification methods were followed by normalization compared using Euclidean distance. The simulation study also investigated the impact of not accounting for the loss of degrees of freedom due to normalization on testing (i.e., impact on the type I error rate). The empirical type I error rate was computed for each “null” gene in which the proportion of the simulated datasets (out of 1,000 simulations) with differential gene expression p-value less than 0.05. Then, for the set of simulated null genes, the average type I error rate was computed by averaging the individual “null” gene type I error rates. In the simulation of the data, we considered two main scenarios: (I) only the batch effect(s) simulated (no primary biological factor) and (II) batch effects plus the effect of a biological variable of interest, where batch and biological effect were uncorrelated. The technical artifact was simulated to represent different mechanisms: discrete number of batches or runs of the samples (e.g., two groups) or a trend effect due to time of run with a continuous effect. The primary biological effect simulated was a binary factor, such as a treatment group and a control group. Note that in the simulation study, the genes with p-values less than 0.05 are considered to be DE.

Nous avons g = 1,…,g gènes, m = 1,…,N samples, k = 1,2 biological groups, and je = 1,…,L batches. Laisser Xgnkl be the count for gene g in biological group k, sample m, and batch je, with a Negative Binomial distribution: Xgnkl

NegBin(??gnkl,??g). The parameters ??gkl et ??g are the mean and dispersion, respectively. Under each scenario (I or II), we changed the sample size (N = 20, 50, 100, 200), the percentage of genes were affected by the batch effect(s) (5%, 10%, 15%), and the percentage of DE genes (0%, 3%, 5%, 10%, 15%). For each scenario, 1,000 datasets were generated, where for each dataset we simulated expression levels for G = 1,000 genes. The baseline parameters and (no batch or biological effects, “null” hypothesis) were estimated using the maximum likelihood estimation (MLE) for the keratin-high samples from the CESC data. The different “non-null” simulated datasets were generated as follows.

Under the scenario I, we considered the binary batch (i.e., sequencing in different labs), continuous batch (e.g., time of day), and both binary and continuous batches. While under the scenario II, we took into account the effect of primary biological factor along with the batch effects (both binary and continuous). Then, the non-null genes (i.e., DE genes) were affected by the batch effect(s) or/and primary biological factor were generated using a mean shift and dispersion , where ??j = (??1j,??2j,…,??Gj) et ??j = (??1j,??2j,…,??Nj) represent the j e effect of batch or primary biological variable, and J represents the total number of variables as batch and primary biological in the study (see S1 Table).

For the scenario I, to generate only the binary batch variable (i.e., je = 1,2),??1 was generated from the Bernouli distribution and ??1 was generated from the Normal distribution (?? = 0,?? = 2). While for the continuous batch variable (i.e., je = 1,…,N),??1 was generated from the standard Uniform distribution, U(0,1), and ??1 were generated from the Normal distribution (?? = 0,?? = 6). If we consider both binary and continuous variables (J = 2), then ??j et ??j for j = 1,2 were generated as explained before. Lastly, for the scenario II, to take into account the effect of biological factor along with the batch effects (both binary and continuous), then ??3 et ??3 were generated from the Bernoulli distribution and the Normal distribution (?? = 0,?? = 2), respectively. The reason of changing the values of ?? is that to provide the moderate effects for batch and primary biological variable under each scenario. The code to generate the simulated data is available at S1 File.


Usage

a DESeqDataSet, or matrix of counts

logical, whether to blind the transformation to the experimental design. blind=TRUE should be used for comparing samples in an manner unbiased by prior information on samples, for example to perform sample QA (quality assurance). blind=FALSE should be used for transforming data for downstream analysis, where the full use of the design information should be made. blind=FALSE will skip re-estimation of the dispersion trend, if this has already been calculated. If many of genes have large differences in counts due to the experimental design, it is important to set blind=FALSE for downstream analysis.

by default, this is not provided and calculated automatically. if provided, this should be a vector as long as the number of rows of object, which is log2 of the mean normalized counts from a previous dataset. this will enforce the intercept for the GLM, allowing for a "frozen" rlog transformation based on a previous dataset. You will also need to provide mcols(object)$dispFit .

a single value, the variance of the prior on the sample betas, which if missing is estimated from the data

in case dispersions have not yet been estimated for object , this parameter is passed on to estimateDispersions (options described there).


Differential gene expression analysis

Differential expression analysis means taking the normalised read count data and performing statistical analysis to discover quantitative changes in expression levels between experimental groups. For example, we use statistical testing to decide whether, for a given gene, an observed difference in read counts is significant, that is, whether it is greater than what would be expected just due to natural random variation.

Methods for differential expression analysis

There are different methods for differential expression analysis such as edgeR and DESeq based on negative binomial (NB) distributions or baySeq and EBSeq which are Bayesian approaches based on a negative binomial model. It is important to consider the experimental design when choosing an analysis method. While some of the differential expression tools can only perform pair-wise comparison, others such as edgeR, limma-voom, DESeq and maSigPro can perform multiple comparisons.

In Figure 11, below, we outline the RNA-seq processing pipeline used to generate data for Expression Atlas.

Figure 11 RNA-seq processing pipeline used to generate gene expression data in Expression Atlas.

In this pipeline raw reads (FASTQ files) undergo quality assessment and filtering. The quality-filtered reads are aligned to the reference genome via HISAT2. The mapped reads are summarised and aggregated over genes via HTSeq. For baseline expression, the FPKMs are calculated from the raw counts by iRAP. These are averaged for each set of technical replicates, and then quantile normalised within each set of biological replicates using limma.

Finally, they are averaged for all biological replicates (if any). For differential expression, genes expressed differentially between the test and the reference groups of each pairwise contrast are identified using DESeq2.


In projects that involve samples from different biological conditions, statistical analyses can be used to identify quantitative changes in gene expression between the different conditions. We perform this analysis using the DESeq2 framework (Love et al., 2014). The main output is a table that contains the average expression, fold-change, and associated statistics such as the P and corrected P values for each gene.


FPKM vs raw read count for differential expression testing

I'm a plant pathology student and we used RNA sequencing to examine differential expression of genes related to a specific pathway in response to a stressor in two plant varieties (a tolerant and susceptible line). I was recently criticized while presenting my RNA-seq results for using FPKM values rather than raw read counts when examining differential expression of a given gene across conditions and with the expression of other genes.

We used Illumina sequencing and mapped using Tophat followed by Cufflinks, then Cuffdiff. This professor was under the impression that FPKM values are too normalized, and that you loose accuracy if the gene is incorrectly annotated and/or the gene is very short or very long. My question, is this a criticism that you have heard before? Is it necessarily wrong to present data as FPKM values? I'm just curious what the opinion is on this subject.

Ok --- so it seems this question comes up beaucoup (I'm going to sit down and write a blog post about this at some point). There are a couple of things worth pointing out from your question.

First, most packages ne pas support the use of TPM or FPKM for differential expression testing. This means that e.g. it's completely wrong to feed them to programs expecting counts (e.g. DESeq2 or EdgeR). One reason for this is that these measures are normalized. What I mean by this is that, for example, if you sum the TPM of all genes/transcripts in a sample, the sum will always be 1,000,000 this is a direct result of the way TPM is calculated. FPKM will behave in a similar manner (though, as many have pointed out, TPM should always be preferred to FPKM which has a more arbitrary and less stable normalization term).

Some tools do make use of FPKM for differential expression testing (the Tuxedo tools), but they also maintain extra information about samples (e.g. the total number of reads) that allow for the proper comparison of these normalized measures across samples. On the other hand, most popular stand-alone DE tools (e.g. DESeq2, EdgeR, etc.) expect counts as input. This is because they perform their own, internal, normalization to help account for effect size (1000 reads coming from a transcript means something different when there are 10 million reads in my sample vs. 50 million). My personal recommendation would be to use salmon (disclosure: I'm the author of this tool P) to process your sample, and then feed the "NumReads" column (appropriately rounded) to a DE tool like DESeq2. That being said, the pipeline that it sounds like you followed (i.e. TopHat => Cufflinks => CuffDiff) is ne pas unreasonable --- these tools were meant to be used together. While some recent surveys have shown that they can be outperformed by other methods, I don't believe there is anything systematically wrong on that front.

On a finer point, regarding what was said by the professor that

FPKM values are too normalized, and that you loose accuracy if the gene is incorrectly annotated and/or the gene is very short or very long.


Anders, S., and Huber, W. (2010). Differential expression analysis for sequence count data. Génome Biol. 11:R106. doi: 10.1186/gb-2010-11-10-r106

Bacher, R., Chu, L., Leng, N., Gasch, A. P., Thomson, J. A., Stewart, R. M., et al. (2017). SCnorm: robust normalization of single-cell RNA-seq data. Nat. Méthodes 14:584. doi: 10.1038/nmeth.4263

Bullard, J. H., Purdom, E., Hansen, K. D., and Dudoit, S. (2010). Evaluation of statistical methods for normalization and differential expression in mRNA-Seq experiments. BMC Bioinformatique 11:94. doi: 10.1186/1471-2105-11-94

Cole, M. B., Risso, D., Wagner, A., DeTomaso, D., Ngai, J., Purdom, E., et al. (2018). Performance assessment and selection of normalization procedures for single-cell RNA-seq. bioRxiv [Preprint]. doi: 10.1101/235382

Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., et al. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatique 29, 15�. doi: 10.1093/bioinformatics/bts635

Gao, S. (2018). Data analysis in single-cell transcriptome sequencing. Méthodes Mol. Biol. 1754:18.

Gao, S., Ou, J., and Xiao, K. R. (2014). language and Bioconductor in Bioinformatics Applications(Chinese Edition). Tianjin: Tianjin Science and Technology Translation Publishing Ltd.

Gao, S., Tian, X., Chang, H., Sun, Y., Wu, Z., Cheng, Z., et al. (2017). Two novel lncRNAs discovered in human mitochondrial DNA using PacBio full-length transcriptome data. Mitochondrion 38, 41�. doi: 10.1016/j.mito.2017.08.002

Gao, S., Zhang, N., Duan, G. Y., Yang, Z., Ruan, J. S., and Zhang, T. (2009). Prediction of function changes associated with single-point protein mutations using support vector machines (SVMs). Hum. Mutat. 30, 1161�. doi: 10.1002/humu.21039

Glusman, G., Caballero, J., Robinson, M., Kutlu, B., and Hood, L. (2013). Optimal scaling of digital transcriptomes. PLoS Un 8:e77885. doi: 10.1371/journal.pone.0077885

Jiang, L., Schlesinger, F., Davis, C. A., Zhang, Y., Li, R., Salit, M., et al. (2011). Synthetic spike-in standards for RNA-seq experiments. Génome Res. 21:1543. doi: 10.1101/gr.121095.111

Li, P., Piao, Y., Shon, H. S., and Ryu, K. H. (2015). Comparing the normalization methods for the differential analysis of Illumina high-throughput RNA-Seq data. BMC Bioinformatique 16:347. doi: 10.1186/s12859-015-0778-7

Lovén, J., Orlando, D. A., Sigova, A. A., Lin, C. Y., Rahl, P. B., Burge, C. B., et al. (2012). Revisiting global gene expression analysis. Cellule 151, 476�. doi: 10.1016/j.cell.2012.10.012

Lun, A. T., Karsten, B., and Marioni, J. C. (2016). Pooling across cells to normalize single-cell RNA sequencing data with many zero counts. Génome Biol. 17:75. doi: 10.1186/s13059-016-0947-7

Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatique 26, 139�. doi: 10.1093/bioinformatics/btp616

Wu, Z., Liu, W., Jin, X., Yu, D., Wang, H., Glusman, G., et al. (2018). NormExpression: an R package to normalize gene expression data using evaluated methods. bioRxiv [Preprint]. doi: 10.1101/251140

Zhang, M., Zhan, F., Sun, H., Gong, X., Fei, Z., et al. (2014). �stq_clean: an optimized pipeline to clean the Illumina sequencing data with quality control,” in Proceedings of the IEEE International Conference on Bioinformatics and Biomedicine (BIBM), (Piscataway, NJ: IEEE).

Keywords : gene expression, normalization, evaluation, R package, scRNA-seq

Citation: Wu Z, Liu W, Jin X, Ji H, Wang H, Glusman G, Robinson M, Liu L, Ruan J and Gao S (2019) NormExpression: An R Package to Normalize Gene Expression Data Using Evaluated Methods. Devant. Genet. 10:400. doi: 10.3389/fgene.2019.00400

Received: 24 December 2018 Accepted: 12 April 2019
Published: 30 April 2019.

Tuo Zhang, Cornell University, United States

Yudong Cai, Shanghai University, China
Naibin Duan, Shandong Academy of Agricultural Sciences, China

Copyright © 2019 Wu, Liu, Jin, Ji, Wang, Glusman, Robinson, Liu, Ruan and Gao. Il s'agit d'un article en libre accès distribué sous les termes de la Creative Commons Attribution License (CC BY). L'utilisation, la distribution ou la reproduction dans d'autres forums est autorisée, à condition que le ou les auteurs originaux et le ou les titulaires des droits d'auteur soient crédités et que la publication originale dans cette revue soit citée, conformément aux pratiques académiques acceptées. Aucune utilisation, distribution ou reproduction non conforme à ces conditions n'est autorisée.