Réaliser la partie bioinformatique d'une analyse RNAseq

Récupération des données RNAseq :

Les fichiers fastq ont été préparés (téléchargés depuis l'ENA et réduits) et sont déposés sur les ressources du cluster de l'IFB, dans l'espace projet de la formation.

Comme nous venons de copier l'ensemble du dossier dans votre espace projet, les données ont été copiées et sont dans le dossier atelier_rnaseq/01-Bioinfo/data/. Vérifiez que vous êtes bien dans votre répertoire copié et sinon, déplacez-vous y.

Remplacer "<votre_espace_projet>" avant de d'exécuter la commande suivante

# verifiez que vous etes dans votre copie:
pwd
# sinon, allez-y en remplacant <votre_espace_projet> 
# (il ne doit plus y avoir de "<" ni ">")
cd /shared/projects/<votre_espace_projet>/atelier_rnaseq/01-Bioinfo

Listez les données RNAseq, il y a 3 fichiers WT et 3 KO au format fastq:

ls data/

GCF_000001735.4_TAIR10.1_genomic.gtf  KO3_R1.fastq.gz  WT2_R2.fastq.gz
KO1_R1.fastq.gz			      KO3_R2.fastq.gz  WT3_R1.fastq.gz
KO1_R2.fastq.gz			      WT1_R1.fastq.gz  WT3_R2.fastq.gz
KO2_R1.fastq.gz			      WT1_R2.fastq.gz
KO2_R2.fastq.gz			      WT2_R1.fastq.gz

Créez un dossier, results, pour placer les résultats

mkdir results

Les différentes étapes d'une analyse RNAseq sont :

• analyse qualité des séquences
• nettoyage (si besoin)
• mapping des reads sur le génome
• comptage des reads sur chaque gènes

La première étape d'analyse de la qualité des séquences se lance avec le logiciel FastQC. Nous appliquons le nettoyage des reads (même si pas besoin, juste pour montrer) avec le logiciel Fastp. Pour le mapping, nous utilisons l'outil STAR, version 2.7.9a. Pour changer de format, sélectionner des alignements, ou indexer le fichier bam de mapping pour le visualisé, nous utilisons samtools (version 1.15.1) qui permet de travailler avec les fichiers résultants du mapping. Pour le comptage, nous utilisons featureCounts pour le comptage (featureCounts appartient à un ensemble d'outils nommé subread, nous utilisons la version 1.6.1)

Sur le cluster de calcul de l'IFB, il faut activer les outils dont nous avons besoin. Ainsi, nous avons 5 outils à activer avec la commande "module" :

## Charger les outils dont tu as besoin pour ton analyse
module load fastqc/0.11.9 #Controle qualite des sequences
module load fastp/0.23.1 #Nettoyage des sequences (si besoin)
module load star/2.7.9a #Mapping des reads sur le genome
module load samtools/1.15.1 #Selection des alignements
module load subread/1.6.1 #Comptage avec featureCounts

La documentation pour fastQC est là: https://www.bioinformatics.babraham.ac.uk/projects/fastqc/.

On va lancer fastQC sur un échantillon (KO1) en utilisant le répertoire results:

fastqc -h #pour voir toutes les options
fastqc -o results data/KO1_R1.fastq.gz
fastqc -o results data/KO1_R2.fastq.gz

On vérifie les fichiers résultats et en particulier le fichier .html qui va vous donner un recap des résultats. Si besoin (ici pas vraiment besoin), on passe au nettoyage:

La documentation pour fastp est ici: https://github.com/OpenGene/fastp ou:

fastp -h #pour voir toutes les options

Pour lancer sur le même échantillon:

fastp -i data/KO1_R1.fastq.gz -I data/KO1_R2.fastq.gz \
      -o results/KO1_R1.clean.fastq.gz -O results/KO1_R2.clean.fastq.gz \
      --detect_adapter_for_pe 
#Et d'autres paramètres si besoin en fonction du résultat du fastqc

On revérifie les séquences nettoyées avec fastqc

fastqc -h #pour voir toutes les options
fastqc -o results results/KO1_R1.clean.fastq.gz
fastqc -o results results/KO1_R2.clean.fastq.gz

On peut maintenant passer au mapping des reads sur le génome. La documentation de l'outil STAR est ici https://physiology.med.cornell.edu/faculty/skrabanek/lab/angsd/lecture_notes/STARmanual.pdf

L'usage de base de STAR est constituée de 2 étapes :

• indexation du génome
• mapping des reads

L'indexation est un procédé de compression et de ré-organisation des séquences qui permet aux logiciels de gagner beaucoup de temps. Chaque outil a le sien (blast aussi).

Le génome d'Arabidopsis thaliana a déjà été indexé et est accessible dans l'espace des banques de séquence : /shared/bank/arabidopsis_thaliana/TAIR10.1/star-2.7.9a/

Si votre génome n'est pas disponible, l'option "--runMode genomeGenerate" de STAR permet de réaliser l'indexation du génome pour STAR.

Les options de bases pour le lancement de STAR sont décrites ici : https://physiology.med.cornell.edu/faculty/skrabanek/lab/angsd/lecture_notes/STARmanual.pdf#subsection.3.1

• NumberOfThreads : --runThreadN entier
• accès au génome indexé : --genomeDir /path/to/genomeDir
• fichier de reads : --readFilesIn /path/to/read1 [/path/to/read2 ]
• pour des fichiers compressés (.gz) : --readFilesCommand zcat

Nous devons donc les mentionner et adapter leur valeur dans la ligne de lancement.

STAR engendre beaucoup de fichiers en sortie, et l'option STAR --outFileNamePrefix /path/to/output/dir/prefix permet de directmenent les ranger dans un sous-dossier. Avec "prefix" il est possible de spécifier quelques lettres qui correpondent à l'échantillon. Par ex. pour l'échantillon WT1, on peut écrire results/WT1_

STAR a été entraîné par défaut pour le génome humain. Pour Arabidopsis thaliana, nous pouvons (devons) modifier les valeurs par défaut de certaines options :

• maximum intron size --alignIntronMax 1000 (la distribution de la longueur moyenne des introns des gènes d'Arabidopsis thaliana est bimodale, à 80 et à 120 bp), mieux que ~600 Kb par défaut
• maximum gap between two mates --alignMatesGapMax 10000 (defaut ~600 Kb)

Voilà, nous avons tout réuni pour écrire (puis lancer) la commande de mapping. Continuons avec l'échantillon "KO1" :

STAR --genomeDir /shared/bank/arabidopsis_thaliana/TAIR10.1/star-2.7.9a/ --runThreadN 4 \
     --readFilesCommand zcat --outFileNamePrefix results/KO1_ \
     --readFilesIn results/KO1_R1.clean.fastq.gz results/KO1_R2.clean.fastq.gz \
     --outSAMtype BAM SortedByCoordinate --alignIntronMax 1000 --alignMatesGapMax 10000 \
     --sjdbGTFfile /shared/bank/arabidopsis_thaliana/TAIR10.1/gtf/GCF_000001735.4_TAIR10.1_genomic.gtf

Une fois la commande exécutée, il y a plusieurs fichiers dans le dossier "results". Le fichier de sortie d'alignement est celui qui se termine par .bam.

Ce fichier bam n'est pas indexé pour être correctement lu par IGV. Pour cela on utilise samtools index:

samtools index results/KO1_Aligned.sortedByCoord.out.bam

Et enfin, on peut passer au comptage du nombre de reads alignés sur chaque gene avec featureCounts

featureCounts -T 4 -t exon -g gene_id -s 1 \
              -a /shared/bank/arabidopsis_thaliana/TAIR10.1/gtf/GCF_000001735.4_TAIR10.1_genomic.gtf \
              -o results/KO1_feature.out results/KO1_Aligned.sortedByCoord.out.bam 2> results/KO1_counting.logs

Si on récapitule, voici toutes les commandes à lancer pour un échantillon:

#Controle qualite
fastqc -o results data/KO1_R1.fastq.gz
fastqc -o results data/KO1_R2.fastq.gz
#Nettoyage
fastp -i data/KO1_R1.fastq.gz -I data/KO1_R2.fastq.gz \
      -o results/KO1_R1.clean.fastq.gz -O results/KO1_R2.clean.fastq.gz \
      --detect_adapter_for_pe #Et d'autres paramètres si besoin en fonction du résultat du fastqc
#Controle qualite sur fichiers nettoyes
fastqc -o results results/KO1_R1.clean.fastq.gz
fastqc -o results results/KO1_R2.clean.fastq.gz
#Mapping
STAR --genomeDir /shared/bank/arabidopsis_thaliana/TAIR10.1/star-2.7.9a/ --runThreadN 4 \
     --readFilesCommand zcat --outFileNamePrefix results/KO1_ \
     --readFilesIn results/KO1_R1.clean.fastq.gz results/KO1_R2.clean.fastq.gz \
     --outSAMtype BAM SortedByCoordinate --alignIntronMax 1000 --alignMatesGapMax 10000 \
     --sjdbGTFfile /shared/bank/arabidopsis_thaliana/TAIR10.1/gtf/GCF_000001735.4_TAIR10.1_genomic.gtf
#indexation
samtools index results/KO1_Aligned.sortedByCoord.out.bam
#Comptage des reads
featureCounts -T 4 -t exon -g gene_id -s 1 \
              -a /shared/bank/arabidopsis_thaliana/TAIR10.1/gtf/GCF_000001735.4_TAIR10.1_genomic.gtf \
              -o results/KO1_feature.out results/KO1_Aligned.sortedByCoord.out.bam 2> results/KO1_counting.logs

Allons-y pour 2 KO et 2 WT :D, ça donne envie?? (Ne lancez pas ça, on vous montrera plus tard comment faire beaucoup plus simple :))

#KO1
fastqc -o results data/KO1_R1.fastq.gz
fastqc -o results data/KO1_R2.fastq.gz
fastp -i data/KO1_R1.fastq.gz -I data/KO1_R2.fastq.gz \
      -o results/KO1_R1.clean.fastq.gz -O results/KO1_R2.clean.fastq.gz \
      --detect_adapter_for_pe 
fastqc -o results results/KO1_R1.clean.fastq.gz
fastqc -o results results/KO1_R2.clean.fastq.gz
STAR --genomeDir /shared/bank/arabidopsis_thaliana/TAIR10.1/star-2.7.9a/ --runThreadN 4 \
     --readFilesCommand zcat --outFileNamePrefix results/KO1_ \
     --readFilesIn results/KO1_R1.clean.fastq.gz results/KO1_R2.clean.fastq.gz \
     --outSAMtype BAM SortedByCoordinate --alignIntronMax 1000 --alignMatesGapMax 10000 \
     --sjdbGTFfile /shared/bank/arabidopsis_thaliana/TAIR10.1/gtf/GCF_000001735.4_TAIR10.1_genomic.gtf
samtools index results/KO1_Aligned.sortedByCoord.out.bam
featureCounts -T 4 -t exon -g gene_id -s 1 \
              -a /shared/bank/arabidopsis_thaliana/TAIR10.1/gtf/GCF_000001735.4_TAIR10.1_genomic.gtf \
              -o results/KO1_feature.out results/KO1_Aligned.sortedByCoord.out.bam 2> results/KO1_counting.logs

#KO2
fastqc -o results data/KO2_R1.fastq.gz
fastqc -o results data/KO2_R2.fastq.gz
fastp -i data/KO2_R1.fastq.gz -I data/KO2_R2.fastq.gz \
      -o results/KO2_R1.clean.fastq.gz -O results/KO2_R2.clean.fastq.gz \
      --detect_adapter_for_pe 
fastqc -o results results/KO2_R1.clean.fastq.gz
fastqc -o results results/KO2_R2.clean.fastq.gz
STAR --genomeDir /shared/bank/arabidopsis_thaliana/TAIR10.1/star-2.7.9a/ --runThreadN 4 \
     --readFilesCommand zcat --outFileNamePrefix results/KO2_ \
     --readFilesIn results/KO2_R1.clean.fastq.gz results/KO2_R2.clean.fastq.gz \
     --outSAMtype BAM SortedByCoordinate --alignIntronMax 1000 --alignMatesGapMax 10000 \
     --sjdbGTFfile /shared/bank/arabidopsis_thaliana/TAIR10.1/gtf/GCF_000001735.4_TAIR10.1_genomic.gtf
samtools index results/KO2_Aligned.sortedByCoord.out.bam
featureCounts -T 4 -t exon -g gene_id -s 1 \
              -a /shared/bank/arabidopsis_thaliana/TAIR10.1/gtf/GCF_000001735.4_TAIR10.1_genomic.gtf \
              -o results/KO2_feature.out results/KO2_Aligned.sortedByCoord.out.bam 2> results/KO2_counting.logs


#WT1
fastqc -o results data/WT1_R1.fastq.gz
fastqc -o results data/WT1_R2.fastq.gz
fastp -i data/WT1_R1.fastq.gz -I data/WT1_R2.fastq.gz \
      -o results/WT1_R1.clean.fastq.gz -O results/WT1_R2.clean.fastq.gz \
      --detect_adapter_for_pe 
fastqc -o results results/WT1_R1.clean.fastq.gz
fastqc -o results results/WT1_R2.clean.fastq.gz
STAR --genomeDir /shared/bank/arabidopsis_thaliana/TAIR10.1/star-2.7.9a/ --runThreadN 4 \
     --readFilesCommand zcat --outFileNamePrefix results/WT1_ \
     --readFilesIn results/WT1_R1.clean.fastq.gz results/WT1_R2.clean.fastq.gz \
     --outSAMtype BAM SortedByCoordinate --alignIntronMax 1000 --alignMatesGapMax 10000 \
     --sjdbGTFfile /shared/bank/arabidopsis_thaliana/TAIR10.1/gtf/GCF_000001735.4_TAIR10.1_genomic.gtf
samtools index results/WT1_Aligned.sortedByCoord.out.bam
featureCounts -T 4 -t exon -g gene_id -s 1 \
              -a /shared/bank/arabidopsis_thaliana/TAIR10.1/gtf/GCF_000001735.4_TAIR10.1_genomic.gtf \
              -o results/WT1_feature.out results/WT1_Aligned.sortedByCoord.out.bam 2> results/WT1_counting.logs

#KO2
fastqc -o results data/WT2_R1.fastq.gz
fastqc -o results data/WT2_R2.fastq.gz
fastp -i data/WT2_R1.fastq.gz -I data/WT2_R2.fastq.gz \
      -o results/WT2_R1.clean.fastq.gz -O results/WT2_R2.clean.fastq.gz \
      --detect_adapter_for_pe 
fastqc -o results results/WT2_R1.clean.fastq.gz
fastqc -o results results/WT2_R2.clean.fastq.gz
STAR --genomeDir /shared/bank/arabidopsis_thaliana/TAIR10.1/star-2.7.9a/ --runThreadN 4 \
     --readFilesCommand zcat --outFileNamePrefix results/WT2_ \
     --readFilesIn results/WT2_R1.clean.fastq.gz results/WT2_R2.clean.fastq.gz \
     --outSAMtype BAM SortedByCoordinate --alignIntronMax 1000 --alignMatesGapMax 10000 \
     --sjdbGTFfile /shared/bank/arabidopsis_thaliana/TAIR10.1/gtf/GCF_000001735.4_TAIR10.1_genomic.gtf
samtools index results/WT2_Aligned.sortedByCoord.out.bam
featureCounts -T 4 -t exon -g gene_id -s 1 \
              -a /shared/bank/arabidopsis_thaliana/TAIR10.1/gtf/GCF_000001735.4_TAIR10.1_genomic.gtf \
              -o results/WT2_feature.out results/WT2_Aligned.sortedByCoord.out.bam 2> results/WT2_counting.logs