Ce TD a pour objectif de vous initier à l’analyse de données de séquençage haut débit de transcriptomes. Cette stratégie de séquençage est appelée RNA-seq. L’interface Galaxy (lien1 ou lien2) que vous avez découvert au TD précédent sera utilisée pour:

Présentation des données

Commencez par créer un dossier de travail sur votre bureau. Téléchargez ensuite le dossier compressé à l’adresse suivante. Pour décompresser le dossier, ouvrez un terminal (Ctrl + Alt +T) et tapez les commandes suivantes: tar -zxf TD3.tar.gz.

Ce dossier contient plusieurs fichiers préparés à partir d’une étude publiée en 2014 dans PLoS One (Himes et al 2014) et disponible en Open Access. Elle s’intitule RNA-Seq transcriptome profiling identifies CRISPLD2 as a glucocorticoid responsive gene that modulates cytokine function in airway smooth muscle cells. Les données de séquençage haut débit provenant de cette étude ont été déposées dans le SRA sous l’accession SRP033351 et comprend 16 échantillons.

Afin de réduire les temps d’analyse, vous ne traiterez que 8 échantillons. Vous pouvez obtenir des informations dans le fichier sample_table.csv présent dans votre dossier. Les fichiers .csv correspondent à des tableaux dont les colonnes sont séparés par des virgules et peuvent être ouverts dans Excel ou LibreOffice Calc.

Quelles sont les caractéristiques importantes pour ces échantillons?

Traitement dans Galaxy: assemblages, alignements et comptages

La première étape consistera à charger les données dans Galaxy. Il vous faut charger l’ensemble des fichiers .fastq, la séquence chromosomique .fa et le fichier d’annotation .gtf. Ce dernier contient les informations concernant les loci présents sur la séquence chromosomique, c’est-à-dire l’identité des gènes trouvés dans cette région, leur annotation et leurs coordonnées. Pour les fichiers .fastq, pensez à préciser le type, c’est à dire fastqsanger.

Comme vu au précédent TD, vous pouvez vérifier la qualité de séquençage des différents échantillons en utilisant l’outil FastQC.

Vous allez également charger deux fichiers .fastq que vous allez préparer à partir de ceux contenus dans votre dossier. Ces deux fichiers seront la concaténation des fichiers forward d’une part et reverse d’autre part. Ils seront utilisés pour réaliser l’assemblage de novo des transcrits. Pour cela, utilisez le terminal, placez vous dans votre dossier de travail (cd /home/USER/Desktop/dossierdetravail) et tapez les commandes suivantes:cat *_1.fastq >all_1.fastq pour obtenir le premier fichier puis cat *_2.fastq >all_2.fastq. Chargez ensuite ces deux nouveaux fichiers dans Galaxy.

Premier cas: absence de séquence génomique de référence

Assemblage de novo

Pour beaucoup d’organismes, il n’existe pas encore de séquence génomique disponible. Toutefois, une approche RNA-seq permet d’obtenir assez facilement le contenu en transcrit d’un échantillon donné en plus d’une information quantitative (expression des gènes). Dans ce cas, les lectures issues du séquençage doivent être assemblées en séquences plus longues (contigs) dans le but de reconstruire les transcrits dans leur intégralité. Pour maximiser la reconstruction, il peut être utile de combiner les lectures de l’ensemble des échantillons séquencés dans une expérimentation (c’est l’étape que vous venez de réaliser). Ces données combinés sont alors traitées par un assembleur de novo. Le programme que vous utiliserez ici s’appelle Trinity.

Démarrer l’assemblage avec les paramètres par défaut, en précisant les bons fichiers.

Cette étape génère un fichier .fasta contenant les séquences réassemblées. A ce stade, aucune information n’est disponible sur séquence d’où la nécessité de lancer un processus d’annotation, c’est-à-dire attribuer une fonction par des recherches d’homologies. Un moyen simple d’analyser les séquences reconstruites est d’executer un alignement local par BLAST contre les bases de données humaines. Cette étape ne sera pas réalisée ici pour cause de temps. Cependant, si des transcrits s’avèrent être différentiellement exprimés, vous pourrez toujours les analyser individuellement.

Dans un tableur, récupérer l’ensemble des informations obtenues par BLAST pour chaque transcrit.

Alignement à l’assemblage de novo

L’estimation de l’abondance de chaque transcrit dans chacun des échantillons se fait après mapping des lectures sur le transcriptome reconstruit. L’expression est obtenue en fonction de la couverture du transcrit par les lectures. Plus des lectures sont obtenues pour un transcrit dans une condition donnée, plus celui-ci est exprimé dans cette condition.

Commencez par aligner les lectures de chaque échantillon au transcriptome assemblé de novo en utilisant l’outil Bowtie2 comme vu au TD précédent. Cependant, plutôt que de lancer les 8 alignements manuellement, vous pouvez créer, dans le panneau de droite après avoir coché l’option Operations on multiple datasets, un dataset collection. Sélectionnez l’ensemble des fichiers SRR.fq (sans les fastq concaténés ayant à l’assemblage), et exécuter Build dataset list. Dans l’outile Bowtie2, sélectionnez l’assemblage de novo au format .fasta comme séquence de réféfence. Vous obtiendrez au final 8 alignements au format .bam.

Vous pouvez visualiser pour un échantillon par exemple, le mapping des lectures sur le transcriptome dans IGV. Pour cela charger le transcriptome assemblé de novo au format .fasta avec la fonction “load genome from file” et l’alignement avec “load from file”.

Traitez ensuite la collection d’alignements avec l’outil de comptage eXpress pour estimer les niveaux de chaque transcrit. Ce dernier outil génère 2 fichiers, l’un contenant les niveaux d’expression, mesurés en FPKM (Fragments per Kilobase per Million of reads), et l’autre les données de distribution et d’estimation d’erreur utilisés utilisé pour l’estimation de l’abondance. Si vous utilisez l’instance PRABI, vous ne pourrez utiliser que RSEM qui s’utilise en 3 étapes: (i) préparation de la référence, (ii) génération d’une carte gène-isoforme, (iii) quantification.

Dans votre dossier de travail, créez un répertoire appelé galaxy_express_denovo (mkdir galaxy_express_denovo) et enregistrer les tables contenant les niveaux d’expression dedans. Renommez-les avec le nom de l’échantillon (SRRXXXX).

Second cas: séquence de référence disponible

Lorsqu’une séquence génomique est disponible l’analyse est plus simple car il suffit d’aligner les lectures issues du séquençage sur cette séquence. De plus, des données d’annotation y sont généralement associées, sous forme fichier .gtf. Ce type de fichier contient des informations structurelles sur les gènes (exons, introns,…) ainsi que leur identité. Ces fichiers sont nécessaire pour délimiter des entités sur lesquels les lectures seront cartographiées, permettant ensuite d’estimer leur abondance. Vous allez utiliser le pipeline Tophat-Cufflinks pour compter l’expression des transcrits déjà connus à partir des données RNA-seq.

Analyser chacun des échantillons en utilisant Tophat sur le dataset collection que vous avez préparé précédemment avec les paramètres par défaut. La séquence de référence à utiliser est le fichier .fa appelé Hs.GRCh37_chr1.fa.

Le programme Tophat utilise Bowtie2 (vu au 1er TD) pour aligner les lectures sur la séquence de référence mais est capable d’optimiser l’alignement par rapport à la présence des introns. Il retourne plusieurs fichiers, dont un alignement intitulé accepted hits. Un autre fichier contient les données de jonction, c’est-à-dire la présence des introns dans la séquence de référence. Vous pouvez également consulter le fichier de statistiques qui rapporte notamment le % de lectures correctement alignées.

Vous pouvez visualiser pour un échantillon par exemple, le mapping des lectures sur le génome dans IGV. Pour cela charger l’extrait génomique au format .fasta avec la fonction “load genome from file” et l’alignement que vous avez enregistré avec son index avec “load from file”. Pensez à charger également le fichier d’annotation GFF pour voir la position des gènes sur la séquence de référence.

Pour quantifier les niveaux d’expression sur chacun des exons décrits dans le fichier .gtf mais aussi obtenir la quantification sur les gènes dans leur intégralité, vous allez utiliser le programme Cufflinks. Exécuter le programme sur la collection d’alignements délivrée par Tophat en précisant le fichier Homo_sapiens.GRCh37.75_subset_ok.gtf comme fichier d’annotation (dans l’option, Use reference annotation). Cufflinks retourne plusieurs fichiers, en particulier les niveaux d’expression mesurés par transcrit et par gène.

Dans votre dossier de travail, créez un répertoire appelé galaxy_cufflinks (mkdir galaxy_cufflinks) et enregistrer les tables contenant les niveaux d’expression des transcrits dedans. Renommez-les avec le nom de l’échantillon (SRRXXXXX).

Analyse dans R: visualisation des données

Préparation des tables de données

Vous avez généré une table d’expression pour chaque échantillon sous Galaxy. Avec l’aide de R, vous allez créer deux tables combinant ces données d’expression. La première contiendra les valeurs de TPM mesurés sur le transcriptome assemblé de novo et la seconde les données obtenues à partir de la séquence de référence. A ce stade, votre dossier sur le bureau contient deux autres répertoires contenant ces différentes tables.

Ouvrez RStudio , une interface graphique qui vous aidera à utiliser le langage R. Il s’agit bel et bien d’un langage de programmation, orienté objet et destiné initialement pour des études statistiques. Sa facilité de manipulation lui ont permis d’être utilisé plus largement et le langage s’est enrichi progressivement par la publication de paquets (packages) qui élargissent ses champs d’utilisation. En particulier l’extension Bioconductor a été développée plus spécifiquement pour l’analyse in silico de données biologiques (bioinformatique). Bioconductor est une base de donnée regroupant de nombreux paquets pour des applications bioinformatiques, dont l’analyse de données NGS, mais aussi des feuilles de routes et des exemples.

Créez un nouveau script dans lequel vous rentrerez les différentes commandes indiquées ci-après. Rappel: une seule commande par ligne, tapez run ou Ctrl+entrée pour exécuter la ligne. 

# Commencez par charger le fichier de description 'sample_table.csv'
ech <- read.csv("sample_table.csv", row.names = 1)

# lister les fichiers dans le répertoire
fichiers.denovo <- list.files("galaxy_express_denovo/", pattern = ".txt")
# lire récursivement chacun des fichiers et récupérer les colonnes
# intéressantes
rn.denovo <- rownames(read.table(paste("galaxy_express_denovo/", fichiers.denovo[1], 
    sep = ""), header = T, row.names = 2))
table.denovo <- do.call("cbind", lapply(fichiers.denovo, function(x) read.table(paste("galaxy_express_denovo/", 
    x, sep = ""), header = T, row.names = 2)[, "fpkm"]))
rownames(table.denovo) <- rn.denovo
colnames(table.denovo) <- fichiers.denovo

Executer une commande similaire pour créer une table d’expression des transcrits obtenue après mapping sur la séquence de référence. 

# lister les fichiers dans le répertoire
fichiers.ref <- list.files("galaxy_cufflinks/")
# lire récursivement chacun des fichiers et récupérer les colonnes
# intéressantes
table.ref <- do.call("cbind", lapply(fichiers.ref, function(x) read.table(paste("galaxy_cufflinks/", 
    x, sep = ""), header = T, row.names = 1)$FPKM))
rownames(table.ref) <- rownames(read.table(paste("galaxy_cufflinks/", fichiers.ref[1], 
    sep = ""), header = T, row.names = 1))
colnames(table.ref) <- fichiers.ref

Visualisation

La visualisation des données est très importante car elle va notamment permettre de s’assurer de la validité des approches. L’intérêt ici est de voir comment se comportent les gènes et les différents échantillons. Vous commencerez par observation les niveaux de corrélation entre variables (échantillons). Ceci devrait permettre de s’assurer que les niveaux d’expression reflètent bien les différences de traitement dans le design de l’expérimentation. 

# Chargez le paquet pour la visualisation
library(gplots)
heatmap.2(cor(log2(table.denovo + 1)), trace = "none", margins = c(15, 15))
heatmap.2(cor(log2(table.ref + 1)), trace = "none", margins = c(15, 15))

La visualisation que vous obtenez s’appelle heatmap. Elle permet d’afficher l’intensité d’une mesure dans une condition croisée. Dans le cas présent, vous visualisez une matrice de corrélation, c’est-à-dire les valeurs de corrélation (covariance) pour toutes les combinaisons de variable (échantillons) pris deux à deux. Plus les variables sont corrélées positivement, plus le coefficient de variation tend vers 1.

Observez les deux matrices et confrontez par rapport à la nature des échantillons (traités, non traités).

Il est également possible de visualiser, toujours à l’aide de heatmaps, les niveaux de chacun des transcrits. 

heatmap.2(table.denovo, trace = "none", margins = c(15, 15))
heatmap.2(table.ref, trace = "none", margins = c(15, 15))

Les graphiques ainsi obtenus ne sont pas très exploitables. En effet, les FPKM mesurés sont compris entre 0 et + l’infini. Dans le cadre du TD, les échantillons de séquençage ont été tronqués (ils ne contiennent pas l’ensemble des données de séquençage). Ceci impact grandement les niveaux d’expression, car vous voyez que beaucoup de transcrits sont à 0 dans la table obtenue avec la référence. Au contraire, à partir du transcriptome de novo, chaque transcrit dispose d’un niveau d’expression.

Statistiques

Le design expérimental utilisé vise à comparer les réponses de lignées de cellules musculaires à des traitements anti-inflammatoires, ici le dexaméthasone. Vous allez donc analyser si ce traitement influe sur les niveaux d’expression des transcrits mesurés par RNA-seq. Cette analyse statistique restera assez simple ici dans le cadre du TP et se bornera à construire un modèle linéaire qui va faire ressortir les effets du design expérimental sur l’expression des gènes. Ce type d’approche ne prend tout son sens que lorsque des réplicats biologiques sont disponibles. 

# Chargement du module reshape2
library(reshape2)

# Préparation de la table de novo
table.denovo.m <- melt(cbind.data.frame(t(table.denovo), Traitement = ech[rownames(t(table.denovo)), 
    "dex"], Cell = ech[rownames(t(table.denovo)), "cell"]))
# Résultats du modèle linéaire
summary(lm(value ~ Traitement * variable, data = table.denovo.m))
# Visualisation par boxplot-exemples de transcrits variant significativement
boxplot(table.denovo["comp4_c0_seq1", ] ~ ech[colnames(table.denovo), "dex"])
boxplot(table.denovo["comp45_c0_seq1", ] ~ ech[colnames(table.denovo), "dex"])
boxplot(table.denovo["comp6_c0_seq1", ] ~ ech[colnames(table.denovo), "dex"])

# Préparation de la table référence
table.ref.m <- melt(cbind.data.frame(t(table.ref), Traitement = ech[rownames(t(table.ref)), 
    "dex"], Cell = ech[rownames(t(table.ref)), "cell"]))
summary(lm(value ~ Traitement * variable, data = table.ref.m))
# Visualisation par boxplot-exemples de transcrits variant significativement
boxplot(table.ref["ENST00000376957", ] ~ ech[colnames(table.ref), "dex"])
boxplot(table.ref["ENST00000376819", ] ~ ech[colnames(table.ref), "dex"])
boxplot(table.ref["ENST00000240185", ] ~ ech[colnames(table.ref), "dex"])

L’analyse des p-values (significativité du test de comparaison de moyennes) indique que l’expression de quelques gènes varient en fonction du traitement appliqué. Vous pouvez récupérer les séquences des transcrits assemblés de novo et les analyser par BLAST sur ENSEMBL contre le génome humain. Vous devriez trouver des similarités en termes d’expression mesurée et de séquence, par exemple entre ENST00000376957 et comp4_c0_seq1.

Bilan

Les données utilisées au cours du TD sont insuffisantes pour approfondir l’analyse et découvrir de nouveaux mécanismes impliqués dans les réponses des lignées cellulaires au traitement. Toutefois, vous noterez que des niveaux d’expression ont pu être comparés entre deux approches différentes dans le traitement des données RNA-seq. Ceci vous a permis de pratiquer l’analyse de données RNA-seq comme elle peut se faire à l’heure actuelle. Les analyses de ce type sont notamment mises en place dans des études de pharmacogénomique pour mieux comprendre les différences de réponse aux traitements entre patients.

Ce document a été généré avec RStudio dans le langage Markdown