In biomics-pasteur-fr/RNADiff: R Tools for RNA-Seq Projects at Institut Pasteur (Biomics Platform)

Mode d'emploi du package RNADiff pour l'analyse différentielle : comparaison de 2 conditions avec DESeq2

Version de RNADiff : r packageVersion("RNADiff")

Auteurs : Marie-Agnès Dillies et Hugo Varet (Plate-forme Transcriptome et Epigénome - Institut Pasteur)

1 Introduction

Ce document est le mode d'emploi du package RNADiff développé à la plate-forme pour l'analyse différentielle de données RNA-seq en routine. Seuls les projets qui consistent à comparer deux conditions biologiques (répliquées ou non) sont concernés par ce documents. Bien que le terme ne soit pas approprié, nous appellerons pipeline cette suite de commandes R. Ce pipeline permet de réaliser l'analyse statistique des données et de produire un rapport et des slides de présentation au format PDF.

L'analyse différentielle des données et la génération des documents de résultats se font en deux temps :

• l'analyse statistique proprement dite,
• la génération du rapport et des slides rendus aux utilisateurs : ceux-ci contiennent la description des étapes et méthodes d'analyse ainsi que les figures générées par l'analyse statistique.

Dans un premier temps, on met au point l'analyse statistique à l'aide d'un script R. Cette première étape permet de vérifier la pertinence des paramètres par défaut, et de les ajuster si nécessaire. Une fois ces paramètres fixés et l'analyse réalisée, les paramètres sont automatiquement sauvegardés par R et récupérés par knitr pour la génération du rapport et des slides.

Dans ce document, la section 2 décrit l'ensemble des logiciels et fichiers nécessaires à la réalisation complète de l'analyse. La section suivante détaille la procédure de l'analyse statistique et la section 4 décrit la production du rapport et des slides).

2 De quoi a-t-on besoin pour réaliser l'analyse ?

Les logiciels

Le processus complet nécessite que soient installés sur votre ordinateur les logiciels suivants :

• le logiciel R, avec les packages DESeq2, edgeR, genefilter, xtable (disponibles en ligne) et RNADiff développé à la plate-forme (cf GitHub pour l'installation et le choix des versions),
• LaTeX, afin de générer le rapport et la présentation,
• Acrobat Reader ou un visionneur pour la visualisation des fichiers PDF.

Le script R d'analyse

Le seul script à modifier par l'utilisateur se nomme script_DESeq2_1factor.R. Il contient le code nécessaire à l'analyse statistique, la production des figures, la sauvegarde des résultats et la génération du rapport et des slides au format LaTeX. Les fichiers dits "systèmes" sont inclus dans le package RNADiff :

• les fonctions R développées pour réaliser chaque étape de l'analyse statistique,
• le script knitr report1factor.Rnw permettant de générer le rapport de résultats,
• le script knitr slides1factor.Rnw permettant de générer les slides de présentation des résultats,
• le fichier NGS.bib contenant les références bibliographiques nécessaires à la rédaction du rapport.

Les fichiers de données

La réalisation de l'analyse suppose le calcul préalable des comptages bruts par gène (ou par transcrit). RNADiff autorise l'importation de fichiers de comptages générés par HTSeq-count ou par featureCounts. Si le comptage a été réalisé par HTSeq-count, les fichiers de comptage sont au format suivant :

• un fichier de comptages par échantillon sans en-tête (noms de colonnes),
• au format texte tabulé,
• la première colonne contient les identifiants uniques des gènes ou des transcrits,
• la seconde colonne contient les comptages bruts associés à ces gènes ou transcrits (nombres entiers positifs).

Si le comptage a été réalisé par featureCounts, les fichiers de comptage sont au format suivant :

• un fichier de comptages par échantillon avec en-tête (noms de colonnes),
• au format texte tabulé,
• la première colonne contient les identifiants uniques des gènes ou des transcrits,
• la septième colonne contient les comptages bruts associés à ces gènes ou transcrits (nombres entiers positifs).

Ces fichiers de comptages seront rassemblés dans un dossier appelé raw par défaut. Un autre fichier texte (appelé target.txt par défaut) est utilisé pour décrire le plan de l'expérience, c'est-à-dire le nom de la condition biologique associée à chaque échantillon. Ce fichier au format texte tabulé contient au moins trois colonnes (avec en-têtes) et doit se terminer par un retour à la ligne après la dernière ligne (dans le cas contraire, un message d'avertissement apparaîtra à l'exécution du script R) :

• première colonne : identifiant (unique) de l'échantillon. C'est ce label qui apparaîtra sur toutes les figures pour désigner l'échantillon,
• deuxième colonne : nom du fichier de comptages bruts associé à l'échantillon,
• troisième colonne : condition biologique de l'échantillon.

Voici un exemple du contenu de ce type de fichier :

| label | files | group | |:------|:-----------------------------|:------| | s1c1 | count_file_sample1_cond1.txt | cond1 | | s2c1 | count_file_sample2_cond1.txt | cond1 | | s1c2 | count_file_sample1_cond2.txt | cond2 | | s2c2 | count_file_sample2_cond2.txt | cond2 |

Les fichiers optionnels

geneLengthFile

L'utilisateur peut souhaiter normaliser ses données par la longueur des gènes, afin de comparer l'expression de différents gènes au sein d'un échantillon donné. En aucun cas cette normalisation ne sera utilisée pour l'analyse différentielle. Elle servira seulement à produire une table de données normalisées par gène et par condition. Cette procédure est décrite plus en détail dans le rapport d'analyse.

L'information sur la longueur des gènes sera fournie sous forme d'un fichier texte tabulé à deux colonnes, l'une contenant l'identifiant du gène (identique à l'identifiant utilisé dans les tableaux de comptages bruts), l'autre la longueur du gène en nucléotides. Le nom de ce fichier sera spécifié lors du paramétrage du script d'analyse.

Informations sur les gènes

Plutôt que d'avoir seulement les identifiants des gènes dans les tableaux de résultats, il est possible d'ajouter des informations supplémentaires contenues un fichier externe. Celui-ci devra impérativement être au format texte tabulé et la première colonne devra contenir l'identifiant unique de chaque gène tel qu'écrit dans les fichiers de comptages bruts. Les colonnes suivantes pourront contenir, par exemple, le nom du gène, le numéro d'accession, etc.

Architecture conseillée

Une architecture conseillée pour l'ensemble de ces fichiers et dossiers pourrait être celle décrite dans la figure ci-dessous. Elle permet de séparer les fichiers donnés en entrée de l'analyse et les fichiers de résultats spécifiques d'une version. Pour chaque projet, la séparation entre fichier bruts et fichiers de résultats permet également de clarifier l'organisation et d'accéder plus facilement aux informations recherchées.

Une fois que l'ensemble des fichiers nécessaires a été rassemblé et organisé, l'analyse peut commencer.

3 Etape 1 : l'analyse statistique

Créez tout d'abord un dossier dédié à l'analyse que vous allez réaliser. Supposons que le projet à analyser soit le projet SLX080. Pour simplifier, appelons ce dossier SLX080. L'analyse statistique proprement dite est réalisée à l'aide du fichier script_DESeq2_1factor.r. Pour la première version d'analyse, commencez par copier ce fichier dans le répertoire version1 du projet, puis renommez-le, par exemple avec le nom du projet suivi du nom de la version de l'analyse (même si vous prévoyez plusieurs versions). Dans la suite, nous supposerons que ce script s'appelle maintenant script-version1.R et se trouve dans le répertoire version1 du dossier SLX080.

L'analyse se déroule en deux étapes : l'initialisation des paramètres puis l'exécution du script. Double-cliquez sur le fichier script-version1.R (ou faites un clic droit sur le fichier, ouvrir avec..., R). Le script s'ouvre dans un éditeur de texte, et une fenêtre de commande R apparaît. Placez-vous dans l'éditeur de texte pour modifier le script.

Paramétrage de l'analyse

Les paramètres à configurer se trouvent au début du script :

• workspace : chemin d'accès à l'espace de travail de la session R (répertoire version1 du dossier SLX080),
• projectName : nom du projet qui apparaîtra ensuite dans les titres du rapport et des slides ainsi que sur l'ensemble des tableaux et figures générés. Par souci de cohérence on peut tout simplement utiliser "SLX080",
• analysisVersion : nom de version de l'analyse, par exemple "version1" (ne peut pas contenir de ".", à renseigner même si une seule version est prévue),
• author : nom de l'auteur du rapport statistique,
• researcher : nom du chercheur en charge du projet,
• chief : nom du chef du laboratoire du chercheur,
• rawDir : chemin d'accès au dossier contenant les fichiers de comptages bruts ("raw" par défaut),
• varInt : nom de la variable d'intérêt dans le fichier target.txt ("group" par défaut),
• condRef : condition biologique de référence : sera au dénominateur pour le calcul des fold-changes et des log-ratios. Le nom indiqué ici doit être présent dans la colonne varInt du fichier target.txt,
• batch : nom de l'effet batch sur lequel ajuster l'analyse statistique (NULL par défaut),
• targetFile : chemin d'accès au fichier de description du plan d'expérience décrit en section 2.3 ("target.txt" par défaut),
• infoFile : chemin d'accès au fichier au format texte tabulé contenant les informations sur les gènes, la première colonne contenant les identifiants,
• outfile : TRUE pour exporter les figures au format PDF ou FALSE pour les afficher dans la console R (faire apparaître les figures à l'écran peut être intéressant lorsque l'on exécute le script pas à pas),
• colors : liste des deux couleurs à utiliser pour les figures,
• cooksCutoff : NULL pour laisser DESeq2 utiliser son propre seuil, Inf pour ne pas détecter les gènes outliers,
• independentFiltering : TRUE pour exécuter l'independent filtering : filtre les gènes pour lesquels la moyenne des comptages est inférieure à un certain seuil, ces gènes étant peu susceptibles d'être différentiellement exprimés,
• allComp : TRUE pour tester toutes les comparaisons ou FALSE pour tester seulement les comparaisons à la modalité de référence,
• alpha : seuil de significativité (appliqué sur la p-valeur ajustée) pour la sélection des gènes ou transcrits différentiellement exprimés (0.05 par défaut),
• adjMethod : méthode d'ajustement des p-valeurs brutes pour le contrôle du taux de faux positifs ("BH" ou "BY"),
• type.trans : méthode de transformation à utiliser en vue de la classification et de l'ACP ("VST" pour une variance stabilizing transformation ou "rlog" pour regularized log transformation),
• locfunc : fonction à utiliser pour l'estimation des size factors ("median" par défaut, ou "shorth" du package genefilter),
• geneLengthFile : chemin d'accès au fichier contenant les longueurs des gènes ou des transcrits, en nucléotides (par défaut, aucun fichier de ce type n'est fourni et la valeur de ce paramètre est NULL),
• interestingFeatures : liste de gènes pour lesquels il est intéressant de tracer l'évolution à travers les conditions biologiques,
• featuresToRemove : liste des noms des séquences à supprimer (lignes inutiles des fichiers de comptages ou ARN ribosomique toujours présentes dans le jeu de données par exemple). La présence de ces séquences n'est pas connue à l'avance, elle peut être détectée pendant l'analyse par la présence d'une séquence majoritaire captant plus de 50% de la totalité des reads dans la plupart des échantillons. On peut spécifier ces séquences en listant toutes les séquences d'ARN ribosomiques connues (c("ARNr1", "ARNr2", "ARNr3")) ou bien en indiquant une chaine de caractères commune à tous les noms de ces séquences ("ARNr"). Dans ce dernier cas, l'ensemble des séquences dont le nom contient cette chaine de caractères sera exclu de l'analyse.
• fitType : type de modélisation utilisée pour la relation moyenne-dispersion ("parametric" par défaut, ou "local").

Ces paramètres seront automatiquement repris par les scripts de production du rapport et des slides, sans intervention de votre part. Pour permettre la génération du rapport et des slides, le script doit avoir été exécuté au moins une fois avec le paramètre outfile=TRUE (valeur par défaut). Lorsque tous les paramètres ont été ajustés, n'oubliez pas de sauvegarder le script pour enregistrer les modifications.

Exécution du script

Pour exécuter le script, deux solutions sont possibles :

• se placer dans la console R et lancer la commande source("script-version1.R"),
• ou bien copier-coller le contenu du fichier script-version1.R dans la console R.

Les résultats suivants apparaissent également au fur et à mesure de l'exécution :

• le contenu du fichier target.txt rappelant les fichiers de comptages lus et la condition biologique associée à chaque échantillon,
• les premières lignes de la table de comptages générée à partir de l'ensemble des fichiers de données brutes (cette table doit contenir une colonne par échantillon et les noms des lignes doivent être les identifiants des gènes),
• le nombre de gènes exclus de l'analyse : ce sont les gènes ou transcrits n'ayant aucun read dans tous les échantillons,
• les coefficients SERE calculés entre toutes les paires d'échantillons,
• les coefficients de normalisation par échantillon (size factors),
• le nombre de gènes différentiels obtenus.

Quelques messages d'avertissement peuvent également apparaître à la fin de l'analyse, en particulier si certains comptages nuls subsistent dans le tableau de données.

Fichiers générés lors de l'exécution du script

Pendant l'exécution, R génère les fichiers suivants dans le dossier figures si le paramètre outfile vaut TRUE :

• *-barplotNull.pdf : nombre de comptages nuls par échantillon,
• *-barplotTC.pdf : nombre total de comptages par échantillon,
• *-cluster.pdf : classification des échantillons sur comptages transformées (VST ou rlog),
• *-clusterSERE.pdf : classification des échantillons à partir des coefficients SERE,
• *-coefficientsOfVariation.pdf : coefficient de variation par condition,
• *-densplot.pdf : estimation de la distribution des comptages pour chaque échantillon,
• *-diagEigenValuesPCA.pdf : graphique de diagnostic des valeurs propres de l'ACP,
• *-diagSizeFactors.pdf : diagnostic de l'estimation des paramètres de normalisation,
• *-diagSizeFactorsTC.pdf : size factors vs nombre total de reads,
• *-diagLogNormalityDisp.pdf : diagnostic de log-linéarité des dispersions,
• *-dispersionEstimates.pdf : relation moyenne-dispersion,
• *-heatmap.pdf : heatmap construit sur les gènes differentiels,
• *-histoRawp.pdf : histogrammes des p-valeurs brutes,
• *-majSeq.pdf : séquence/gène majoritaire pour chaque échantillon,
• *-nbDiffSeuil.pdf : figure du nombre de gènes différentiels selon le seuil de FDR (usage interne uniquement),
• *-MAplot.pdf : MA-plot,
• *-pairwiseScatter.pdf : relation entre tous les échantillons,
• *-PCA.pdf : premier plan de l'ACP,
• *-raw-boxplot.pdf : boîtes à moustaches sur les comptages brutes,
• *-norm-boxplot.pdf : boîtes à moustaches sur les données normalisées,
• *-plotEvolution.pdf : évolution des comptages pour certains gènes intéressants définis dans les paramètres,
• *-vennDiagram.pdf : diagramme de Venn des comparaisons testées (maximum 5),
• *-volcanoPlot.pdf : volcano-plot.

Les tables de résultats (avec fold-changes, p-valeurs, etc.) au format texte tabulé sont créées dans le dossier tables :

• *.TestVsRef.complete.xls : tous les gènes présents dans les fichiers de comptages,
• *.TestVsRef.down.xls : gènes sous-exprimés dans Test par rapport à Ref,
• *.TestVsRef.up.xls : gènes sur-exprimés dans Test par rapport à Ref.

Si les longueurs des gènes ont été fournies, alors le fichier SLX080-version1.LengthNorm.xls est également créé. Celui-ci contient les données normalisées par la longueur des gènes. Enfin, un fichier au format .RData contenant les paramètres d'analyse et les résultats est sauvegardé. Avec ces fichiers maintenant disponibles vous pouvez générer le rapport et les slides

4 Etape 2 : le rapport et les slides

Le rapport est généré automatiquement à partir du script R et s'appelle report-SLX080-version1.pdf. En revanche, les slides doivent être personnalisées via le fichier slides-SLX080-version1.tex (commentaires de certaines figures ou conclusion de l'analyse). Une fois les slides modifiées, il suffit de compiler avec LaTeX. Pour cela :

• ouvrez le fichier slides-SLX080-version1.tex avec votre éditeur LaTeX préféré (TeXworks, TeXShop ou TeXnicCenter par exemple),
• compilez le fichier TeX,
• avec TeXworks ou TeXShop, compilez au moins une fois avec BibTeX pour initialiser les références,
• recompilez encore une ou deux fois, jusqu'à ce que toutes les références soient indexées correctement,
• les slides sont prêtes.

biomics-pasteur-fr/RNADiff documentation built on Aug. 27, 2020, 12:44 a.m.