analyse-rna-seq

Logo

View the Project on GitHub omics-school/analyse-rna-seq

Contenu :

L’objectif de cette partie est d’analyser les données RNA-seq de O. tauri sur votre machine Unix locale.

Voici une vue d’ensemble des étapes pour analyser ces données de séquençage haut débit :

3.1 Préparer l’environnement de travail

Sous Windows, ouvrez un terminal Ubuntu.

Déplacez-vous dans le répertoire /mnt/c/Users/omics/rnaseq_tauri :

$ cd /mnt/c/Users/omics/rnaseq_tauri

Activez l’environnement conda rnaseq-env :

$ conda activate rnaseq-env

Remarque : contrôlez que le nom de l’environnement conda apparait bien à gauche de l’invite de commande sous la forme : (rnaseq-env).

Vous êtes maintenant prêt à analyser des données RNA-seq 🤠

3.2 Analyse manuelle

Pour cette première analyse, choisissez un seul échantillon contenant des reads, c’est-à-dire un fichier parmi :

reads/SRR2960338.fastq.gz
reads/SRR2960341.fastq.gz
reads/SRR2960343.fastq.gz

3.2.1 Contrôler la qualité des reads

Créez le répertoire reads_qc qui va contenir les fichiers produits par le contrôle qualité des fichiers fastq.gz :

$ mkdir -p reads_qc

Lancez FastQC avec la commande (attention, modifiez cette commande) :

$ fastqc reads/ECHANTILLON.fastq.gz --outdir reads_qc

où reads/ECHANTILLON.fastq.gz est le fichier contenant l’échantillon que vous avez choisi. Pensez à l’adapter avant de lancer votre commande !

FastQC va produire deux fichiers (un fichier avec l’extension .html et un autre avec l’extension .zip) dans le répertoire reads_qc. Si par exemple, vous avez analysé le fichier reads/SRR2960338.fastq.gz, vous obtiendrez les fichiers reads_qc/SRR2960338_fastqc.html et reads_qc/SRR2960338_fastqc.zip.

Utilisez l’explorateur de fichiers de Windows pour vous déplacer dans le bon répertoire (voir exemple). Ouvrez ensuite le fichier .html ainsi créé avec Firefox (en double-cliquant sur le fichier). Analysez le rapport créé par FastQC.

3.2.2 Indexer le génome de référence

L’indexation du génome de référence est une étape indispensable pour accélérer l’alignement des reads sur le génome. Elle consiste à créer un annuaire du génome de référence.

Toujours depuis votre shell Ubuntu et dans le répertoire /mnt/c/Users/omics/rnaseq_tauri, créez le répertoire index :

$ mkdir -p index

Lancez l’indexation du génome de référence.

$ bowtie2-build genome/GCF_000214015.3_version_140606.fna index/O_tauri

Les index sont stockés dans des fichiers dont le nom débute par O_tauri dans le répertoire index.

Calculez la taille total des fichiers index avec la commande :

$ du -ch index/O_tauri*

Comparez la taille totale des index à la taille du fichier contenant le genome (genome/GCF_000214015.3_version_140606.fna).

L’indexation du génome n’est à faire qu’une seule fois pour chaque logiciel d’alignement.

3.2.3 Aligner les reads sur le génome de référence

Créez le répertoire map qui va contenir les reads alignés sur le génome de référence :

$ mkdir -p map

Lancez l’alignement :

$ bowtie2 -p 2 -x index/O_tauri -U reads/ECHANTILLON.fastq.gz -S map/bowtie.sam

Les options utilisées sont :

Toutes les options de Bowtie2 sont détaillées dans la documentation.

Cette étape peut prendre 6 ou 7 minutes. Bowtie n’affiche rien à l’écran lorsqu’il fonctionne. Soyez patient.

À la fin de l’alignement, Bowtie2 renvoie des informations qui ressemblent à :

5470272 reads; of these:
  5470272 (100.00%) were unpaired; of these:
    497809 (9.10%) aligned 0 times
    4603112 (84.15%) aligned exactly 1 time
    369351 (6.75%) aligned >1 times
90.90% overall alignment rate

On obtient ainsi :

Il faut être prudent si le nombre de reads non alignés est trop important (par exemple supérieur à 20%).

3.2.4 Convertir en binaire, trier et indexer les reads alignés

Vous allez maintenant utiliser SAMtools pour :

  1. Convertir le fichier .sam (fichier texte) créé par Bowtie2 en fichier .bam, qui est un fichier binaire, et qui donc prend moins de place sur le disque. 
     $ samtools view -@ 2 -b map/bowtie.sam -o map/bowtie.bam

    Cette étape va prendre plusieurs minutes. Comme votre machine dispose de 4 coeurs, nous allons en utiliser 2 (-@ 2) pour accélérer le calcul.
    Avec la commande du et les bonnes options, comparez la taille deux fichiers map/bowtie.sam et map/bowtie.bam. Quel est le ratio de compression entre les deux formats de fichiers ?

  2. Trier les reads alignés suivant l’ordre dans lequel ils apparaissent dans le génome. 
     $ samtools sort -@ 2 map/bowtie.bam -o map/bowtie.sorted.bam

    Cette étape va prendre ici encore quelques minutes.

  3. Indexer le fichier .bam. Cette étape est indispensable pour visualiser l’alignement avec IGV. 
     $ samtools index map/bowtie.sorted.bam

3.2.5 Compter les reads alignés sur les gènes de O. tauri

Le comptage des reads alignés sur les gènes se fait avec HTSeq.

Toujours depuis votre shell Ubuntu et dans le répertoire /mnt/c/Users/omics/rnaseq_tauri, créez le répertoire count :

$ mkdir -p count

Puis lancez la commande (en une seule ligne) pour compter les reads alignés :

$ htseq-count --stranded=no --type="gene" --idattr="ID" --order=name --format=bam map/bowtie.sorted.bam genome/GCF_000214015.3_version_140606.gff > count/count.txt

HTSeq renvoie le nombre d’annotations trouvées dans le fichier .gff (7659) puis affiche une progression de l’analyse. Les options du programme htseq-count sont décrites dans la documentation.

Déterminez le nombre de reads alignés sur le gène ostta18g01980 avec la commande :

$ grep ostta18g01980 count/count.txt

Vous pouvez également ouvrir le fichier count/count.txt avec la commande less puis chercher le gène ostta18g01980 en tapant /ostta18g01980 puis la touche Entrée (et enfin la touche Q pour quitter).

3.2.6 Visualiser les reads alignés avec IGV

Cette rubrique est facultative. Passez à la section suivante si vous n’avez pas installé IGV et visionné la vidéo correspondante.

Pour visualiser l’alignement des reads sur le génome de référence avec IGV, vous avez besoin des fichiers suivants :

Lancez IGV et visualisez l’alignement des reads sur le génome de référence. Si vous avez oublié comme faire, visionnez la vidéo sur ce sujet qui vous a été proposée précédemment.

Visualisez particulièrement le gène ostta18g01980.

3.3 Automatiser l’analyse : niveau 1

Tout cela est très bien mais les fichiers que vous avez générés (map/bowtie.bam, map/bowtie.sorted.bam, cout/count.txt…) portent des noms qui ne sont pas très informatifs sur l’échantillon dont ils proviennent.

Par ailleurs, entrer toutes ces commandes à la main, les unes après les autres, est pénible et source d’erreurs. Et il y a fort à parier que vous aurez complètement oublié ces commandes dans 1 semaine, voire dans 1 heure. 🤯 C’est parfaitement normal, il n’y a absolument aucun intérêt à se souvenir de toutes ces commandes.

Pour répondre à ces deux problèmes, de gestion de données et d’automatisation, nous allons introduire les notions Bash de variables et de scripts.

Mais d’abord, faites un peu de ménage en supprimant les fichiers créés précédemment :

$ rm -f reads_qc/*fastqc* index/*bt2 map/bowtie* count/count*

💣 Attention à l’utilisation de la commande rm qui supprime définitivement les fichiers. 💣

3.3.1 Variables

Une variable va simplement contenir de l’information qui sera utilisable autant de fois que nécessaire.

Création de variables :

$ toto=33
$ t="salut"

Attention : Il faut coller le nom de la variable et son contenu au symbole =.

Affichage de variables :

$ echo $toto
33
$ echo "$t Pierre"
salut Pierre

La commande echo affiche une chaîne de caractère, une variable, ou les deux.

Pour utiliser une variable (et accéder à son contenu), il faut précéder son nom du caractère $. Attention, ce symbole n’est pas à confondre avec celui qui désigne l’invite de commande de votre shell Linux.

Enfin, une bonne pratique consiste à utiliser une variable avec le symbole $ et son nom entre accolades :

$ echo ${toto}
33
$ echo "${t} Pierre"
salut Pierre

3.3.2 Script

Un script est un fichier texte qui contient des instructions Bash. Par convention, il porte l’extension .sh. L’objectif premier d’un script Bash est d’automatisr l’exécution de plusieurs commandes Bash, la plupart du temps pour manipuler ou analyser des fichiers.

Dans un script Bash, tout ce qui suit le symbole # est considéré comme un commentaire et n’est donc pas traité par Bash.

3.3.3 Analyse RNA-seq

Testez le script script1.sh sur un seul de vos échantillons. Pour cela :

Vérifiez que le déroulement du script se passe bien. Vous avez le temps de prendre un café (~ 15 ‘), voir plusieurs ☕ 🍪 ☕ 🍪.

Évaluez approximativement le temps nécessaire au script 1 pour s’exécuter. ⏱️ À partir de cette valeur, extrapoler le temps nécessaire qu’il faudrait pour analyser les 3 échantillons.

Utilisez enfin la commande tree pour contempler votre travail (ici avec l’échantillon SRR2960338) :

$ tree
.
├── count
│   └── count-SRR2960338.txt
├── genome
│   ├── GCF_000214015.3_version_140606.fna
│   ├── GCF_000214015.3_version_140606.gff
│   └── md5sum.txt
├── index
│   ├── O_tauri.1.bt2
│   ├── O_tauri.2.bt2
│   ├── O_tauri.3.bt2
│   ├── O_tauri.4.bt2
│   ├── O_tauri.rev.1.bt2
│   └── O_tauri.rev.2.bt2
├── map
│   ├── bowtie-SRR2960338.sorted.bam
│   └── bowtie-SRR2960338.sorted.bam.bai
├── reads
│   ├── SRR2960338.fastq.gz
│   ├── SRR2960341.fastq.gz
│   └── SRR2960343.fastq.gz
├── reads_qc
│   ├── SRR2960338_fastqc.html
│   └── SRR2960338_fastqc.zip
└── script1.sh

6 directories, 18 files

Affichez ensuite la taille occupée par chaque sous-répertoire ainsi que la taille totale avec la commande :

$ du -csh *
140K    count
14M     genome
26M     index
315M    map
1,3G    reads
1,1M    reads_qc
4,0K    script1.sh
1,7G    total

3.4 Automatiser l’analyse : niveau 2

Le script précédent est pratique mais il ne conserve pas les informations liées à l’alignement générées par Bowtie2 (nombre de reads non-alignés, alignés une fois…).

Le script 2 répond à ce problème. Téléchargez-le avec la commande :

$ wget https://raw.githubusercontent.com/omics-school/analyse-rna-seq/master/script2.sh

Ouvrez ce script avec nano. Vous remarquerez que la solution proposée pour conserver les informations liées à l’alignement est un peu particulière (2> map/bowtie-${sample}.out). Nous allons en discuter, mais dans un premier temps essayer de comprendre l’explication donnée ici.

3.5 Automatiser l’analyse : niveau 3 (ninja)

Le script précédent est intéressant mais il ne prend en compte qu’un seul échantillon à la fois. Quel ennui !

On aimerait avoir un seul script qui traiterait tous les échantillons qu’on souhaite analyser. Cela est possible avec une boucle. Une boucle permet de répéter un ensemble d’instructions.

Voici un exemple en Bash :

$ for prenom in gaelle bertrand pierre
> do
> echo "Salut ${prenom} !"
> done
Salut gaelle !
Salut bertrand !
Salut pierre !

En sacrifiant un peu de lisibilité, la même commande peut s’écrire sur une seule ligne :

$ for prenom in gaelle bertrand pierre; do echo "Salut ${prenom} !"; done
Salut gaelle !
Salut bertrand !
Salut pierre !

Notez l’utilisation du symbole ; pour séparer les différents éléments de la boucle.

Une leçon de Software Carpentry aborde la notion de boucle. Prenez quelques minutes pour parcourir cette leçon et comprendre de quoi il s’agit.

Le script 3 utilise une boucle pour automatiser l’analyse de plusieurs échantillons. Téléchargez-le avec la commande :

$ wget https://raw.githubusercontent.com/omics-school/analyse-rna-seq/master/script3.sh

Ouvrez ce script avec nano. Observez la structure du script et essayez de comprendre son fonctionnement.

La ligne set -euo pipefail tout au début du script va arrêter celui-ci :

C’est une mesure de sécurité importante pour votre script. Si vous le souhaitez, vous pouvez lire l’article de Aaron Maxwell à ce sujet : Use the Unofficial Bash Strict Mode (Unless You Looove Debugging)

Si vous pensez en avoir le temps, lancez le script 3. Comme ce script va automatiser toute l’analyse, il va fonctionner environ 45 minutes et vous aurez peut-être besoin de fermer votre terminal. Pour ne pas arrêter brutalement l’analyse, lancez le script de cette manière :

$ nohup bash script3.sh &

Mais pour autant, n’éteignez pas votre ordinateur !

Le message

nohup: ignoring input and appending output to 'nohup.out'

vous rappelle que les messages qui apparaissaient habituellement à l’écran seront redirigés dans le fichier nohup.out.

3.6 Comparer les versions des logiciels utilisés dans Galaxy (si vous avez du temps)

Connectez-vous maintenant à votre compte sur Galaxy. Essayez de retrouver les versions des logiciels que vous avez utilisés (FastQC, Bowtie2, SAMtools, HTSeq).

Pour ce faire, dans votre History, cliquez sur le nom d’un résultat d’analyse, puis cliquez sur le petit i entouré (ℹ️) et lisez les informations de la section Job Dependencies.

Comparez les versions des logiciels disponibles dans Galaxy et de ceux que vous avez utilisés sur votre machine.