Avertissement

Les réponses aux questions de ce TP peuvent être visibles en appuyant sur le bouton Solution. À n’utiliser que si vous bloquez bien sûr !

1 Rappels sur les ressources de la plateforme Migale

L’objectif de cette partie est de présenter comment lancer des calculs sur l’infrastructure de migale (serveur front).

Créer cette arborescence dans votre espace work

TRAINING/
├── CLUSTER/
└── ANALYSES/

Solution

mkdir -p ~/work/TRAINING/CLUSTER
mkdir -p ~/work/TRAINING/ANALYSES
cd ~/work/TRAINING

Se déplacer dans le répertoire CLUSTER et faire un lien symbolique des fichiers présents dans le répertoire /save_projet/metagenomics_training/cluster/

Solution

cd ~/work/TRAINING/CLUSTER
ln -s /save_projet/metagenomics_training/cluster/* .

Trouver l’environnement conda de blast

Solution

conda info --envs |grep blast
blast-2.12.0             /usr/local/genome/Anaconda3/envs/blast-2.12.0
blast-2.13.0             /usr/local/genome/Anaconda3/envs/blast-2.13.0
blast-2.9.0              /usr/local/genome/Anaconda3/envs/blast-2.9.0
blast-legacy-2.2.26      /usr/local/genome/Anaconda3/envs/blast-legacy-2.2.26
rmblast-2.10.0           /usr/local/genome/Anaconda3/envs/rmblast-2.10.0

Aligner avec blastn query.fasta sur subject.fasta (les indexs blast sont déjà générés), en déportant le calcul sur le cluster en utilisant 4 threads

Solution

qsub -cwd -V -q web.q -N myblast -pe thread 4 -b y "conda activate blast-2.13.0 && blastn -db subject.fasta -query query.fasta -out out.blast -num_threads 4 && conda deactivate"

Comment vérifier que le job s’est bien terminé ?

Solution

# Vérifier le fichier de sortie out.blast
# Vérifier qu'il n'y a pas de message d'erreur dans les fichier .o et .e

Aligner tous les fichiers query*.fasta. Il faut pour cela itérer sur les fichiers d’entrée pour lancer les soumissions

Solution

# Solution 1

for i in query*.fasta ; do qsub -cwd -V -q web.q -N myblast -pe thread 4 -b y "conda activate blast-2.12.0 && blastn -db subject.fasta -query $i -out out-$i.blast -num_threads 4 && conda deactivate" ; done

# Solution 2

for i in query*.fasta ; do echo "conda activate blast-2.12.0 && blastn -db subject.fasta -query $i -out $i.blast -num_threads 4 && conda deactivate" >> myblast.sh ; done
qarray -cwd -V -q web.q -N myblast -pe thread 4 myblast.sh

1.1 Préparation de l’espace de travail et organisation

Note

Pour une question de reproductibilité et de tracabilité, les commandes que nous enverrons sur le cluster seront stockées dans des fichiers .sh, nommés en fonction de l’outil à lancer.

Se déplacer dans le répertoire ANALYSES puis faire un lien symbolique des fichiers présents dans le répertoire /save_projet/metagenomics_training/simulated_reads/

Solution

cd ~/work/TRAINING/ANALYSES
ln -s /save_projet/metagenomics_training/simulated_reads/* .

Créer un répertoire LOGS dans lequel seront stockés les fichiers spécifiques à la soumission des jobs sur le cluster de calcul

Solution

mkdir ~/work/TRAINING/ANALYSES/LOGS

2 Contrôle Qualité des reads

Avertissement

Cette étape sert avant tout à vérifier que les fichiers récupérés et sur lesquels on s’apprête à travailler sont conformes à ce qu’on a demandé au prestataire de séquençage

Choses à vérifier :

on a autant de fichiers que d’échantillons envoyés
on a un nombre de reads suffisants par échantillon
on a des tailles de reads conformes à l’attendu

2.1 Nombre de reads

Combien de reads sont présents dans les fichiers s1_R1.fastq.gz et s1_R2.fastq.gz ?

Étape 1: Utiliser la commande zcat pour lire un fichier compressé et wc -l pour compter le nombre de lignes

Solution

zcat s1_R1.fastq.gz | wc -l

Étape 2: Un read correspondant à 4 lignes dans un fichier FASTQ, il faut diviser le nombre de lignes par 4

Solution

zcat s1_R1.fastq.gz | echo $((`wc -l`/4))

Étape 3: Automatiser le calcul pour tous les fichiers au format FASTQ présents dans le répertoire d’intérêt, et afficher le nom du fichier

Solution

for i in *.fastq.gz ; do echo $i $(zcat $i | echo $((`wc -l`/4))) ; done
# s1_R1.fastq.gz 1000002
# s1_R2.fastq.gz 1000002
# ...

2.2 Contenu des fichiers

Vérifier que les fichiers FASTQ sont conformes et obtenir des statistiques sur leur contenu avec l’outil seqkit

[1]

Étape 1: Trouver l’environnement conda de seqkit

Solution

conda info --envs |grep seqkit
seqkit-2.0.0             /usr/local/genome/Anaconda3/envs/seqkit-2.0.0

Étape 2: Explorer les options de seqkit et identifier les paramètres intéréssants

Solution

conda activate seqkit-2.0.0
seqkit stat -h

#  -a, --all                  all statistics, including quartiles of seq length, sum_gap, N50

Étape 3: Lancer le module stat de l’outil seqkit sur tous les fichiers FASTQ

Solution

echo "conda activate seqkit-2.0.0 && seqkit stat --all *.fastq.gz > seqkit.txt && conda deactivate" > seqkit.sh
qsub -cwd -V -q web.q -N seqkit -o LOGS -e LOGS seqkit.sh

Étape 4: Lire le fichier de sortie

Solution

cat seqkit.txt

file            format  type   num_seqs      sum_len  min_len  avg_len  max_len  Q1   Q2  Q3  sum_gap  N50  Q20(%)  Q30(%)
s1_R1.fastq.gz  FASTQ   DNA   1,000,002  126,000,252      126      126      126  63  126  63        0  126   96.11    92.3
s1_R2.fastq.gz  FASTQ   DNA   1,000,002  126,000,252      126      126      126  63  126  63        0  126   94.77   90.19
s2_R1.fastq.gz  FASTQ   DNA   1,000,002  126,000,252      126      126      126  63  126  63        0  126   96.12   92.31
s2_R2.fastq.gz  FASTQ   DNA   1,000,002  126,000,252      126      126      126  63  126  63        0  126   94.76   90.15
s3_R1.fastq.gz  FASTQ   DNA   1,000,000  126,000,000      126      126      126  63  126  63        0  126   96.11    92.3
s3_R2.fastq.gz  FASTQ   DNA   1,000,000  126,000,000      126      126      126  63  126  63        0  126   94.74   90.13
s4_R1.fastq.gz  FASTQ   DNA     500,030   63,003,780      126      126      126  63  126  63        0  126   96.12   92.31
s4_R2.fastq.gz  FASTQ   DNA     500,030   63,003,780      126      126      126  63  126  63        0  126   94.76   90.17
s5_R1.fastq.gz  FASTQ   DNA   1,000,002  126,000,252      126      126      126  63  126  63        0  126   96.11    92.3
s5_R2.fastq.gz  FASTQ   DNA   1,000,002  126,000,252      126      126      126  63  126  63        0  126   94.77   90.19

2.3 FastQC

FastQC [2] est utilisé pour vérifier la qualité des données issues d’un séquençage NGS. La qualité de base des reads, le contenu en adaptateurs, le nombre de séquences uniques… sont des métriques qu’il faut regarder avec attention.

Lancer FastQC sur chaque échantillon

Étape 1: Trouver l’environnement conda de fastqc

Solution

conda info --envs |grep fastqc
fastqc-0.11.8            /usr/local/genome/Anaconda3/envs/fastqc-0.11.8
fastqc-0.11.9            /usr/local/genome/Anaconda3/envs/fastqc-0.11.9

Étape 2: Faire une boucle pour lancer FastQC sur tous les fichiers FASTQ présents dans le répertoire, en utilisant 4 threads

Solution

mkdir FASTQC
for i in *.fastq.gz ; do echo "conda activate fastqc-0.11.9 && fastqc --outdir FASTQC/ --threads 4 $i && conda deactivate" >> fastqc.sh ; done

Étape 3: Lancer les traitements en parallèle sur le cluster grâce à la commande qarray

Solution

qarray -cwd -V -q web.q -N fastqc -pe thread 4 -o LOGS/ -e LOGS/ fastqc.sh

2.4 MultiQC

MultiQC [3] est un outil qui permet d’agglomérer les résultats de FastQC en un seul rapport. Il est particulièrement utile quand on travaille sur beaucoup d’échantillons.

Lancer MultiQC et écrire les résultats dans le répertoire MULTIQC

Solution

conda info --envs |grep multiqc
multiqc-1.11             /usr/local/genome/Anaconda3/envs/multiqc-1.11
multiqc-1.12             /usr/local/genome/Anaconda3/envs/multiqc-1.12

mkdir MULTIQC
echo "conda activate multiqc-1.12 && multiqc --outdir MULTIQC/ FASTQC/ && conda deactivate" > multiqc.sh
qsub -cwd -V -q web.q -N multiqc -o LOGS -e LOGS multiqc.sh

Visualiser le rapport

Solution

firefox --no-remote MULTIQC/multiqc_report.html &

3 Assignation taxonomique

Avant même de nettoyer les reads, il peut s’avérer utile d’effectuer une classification taxonomique. Cela permet entre entre de détecter une contamination (hôte, environnement, humain…). Il faut donc dans l’idéal utiliser la banque la plus complète possible. Nous allons utiliser ici en l’occurence RefSeq et l’outil Kaiju [4].

Les banques spécifiques à kaiju sont mises à disposition dans le répertoire suivant : /db/outils/kaiju/

Obtenir un Krona représentant la diversité de nos échantiilons

Étape 1: Trouver l’environnement conda de kaiju

Solution

conda info --envs |grep kaiju
kaiju-1.8.0              /usr/local/genome/Anaconda3/envs/kaiju-1.8.0
kaiju-1.8.1              /usr/local/genome/Anaconda3/envs/kaiju-1.8.1
kaiju-1.9.0              /usr/local/genome/Anaconda3/envs/kaiju-1.9.0
kaiju-1.9.2              /usr/local/genome/Anaconda3/envs/kaiju-1.9.2

Étape 2: Assigner les reads bruts avec kaiju contre la banque RefSeq. Écrire les résultats dans le répertoire KAIJU

Les banques dédiées à kaiju sont situées dans /db/outils/. La plus récente aujourd’hui est /db/outils/kaiju-2021-03/refseq/.

Solution

mkdir KAIJU
for i in *_R1.fastq.gz ; do id=$(echo $(basename $i | cut -d '_' -f 1)) ; echo "conda activate kaiju-1.9.2 && kaiju -t /db/outils/kaiju-2021-03/refseq/nodes.dmp -f /db/outils/kaiju-2021-03/refseq/kaiju_db_refseq.fmi -i $i -j ${id}_R2.fastq.gz -o KAIJU/${id} -z 16 && conda deactivate" >> kaiju.sh ; done
qsub -cwd -V -q web.q -N kaiju -o LOGS -e LOGS -pe thread 16 kaiju.sh

Étape 3: Corriger les éventuels taxons qui auraient changé entre la version de la banque Kaiju celle utilisée par krona, en utilisant le script /save_projet/metagenomics_training/scripts/update_taxid_column_from_merged_ncbi_file.py.

Solution

cp /save_projet/metagenomics_training/scripts/update_taxid_column_from_merged_ncbi_file.py .
for i in KAIJU/* ; do echo "python3 update_taxid_column_from_merged_ncbi_file.py --dump /db/ncbi_taxon/current/flat/merged.dmp --file $i --taxid-col 3 > $i.cor" >> update_taxo.sh ; done
qarray -cwd -V -q web.q -o LOGS -e LOGS -N updatetaxonomy update_taxo.sh

Étape 4: Générer le fichier Krona avec la suite Krona-tools [5]

Solution

conda info --envs |grep krona
krona                    /usr/local/genome/Anaconda3/envs/krona
krona-2.8                /usr/local/genome/Anaconda3/envs/krona-2.8

for i in KAIJU/*.cor ; do echo "conda activate krona-2.8 && ktImportTaxonomy -tax /db/outils/KRONA-taxo/current/taxonomy/ -t 3 -q 2 $i -o $i.html  && conda deactivate" >> krona.sh ; done
qarray -cwd -V -q web.q -N krona -o LOGS -e LOGS krona.sh

# ou
for i in KAIJU/*.cor ; do echo "conda activate krona-2.8 && ktImportTaxonomy -tax /db/outils/krona-2.8/taxonomy -t 3 -q 2 KAIJU/*.cor -o KAIJU/all.html  && conda deactivate" >> krona_single.sh ; done
qsub -cwd -V -q web.q -N krona -o LOGS -e LOGS krona_single.sh

Étape 5: Visualiser le fichier HTML généré

Solution

firefox --no-remote KAIJU/kaiju.cor.html &
firefox --no-remote KAIJU/all.html &

4 Nettoyage

4.1 Filtrer et trimmer les reads

FASTP [6] est un outils de préprocessing qui permet d’éliminer les séquences de mauvaise qualité, les bases ambigües, les séquences trop courtes, ainsi de les adaptateurs restants.

Explorez les paramètres de FASTP et lancez FASTP sur les données brutes

Étape 1: Trouver l’environnement conda de fastp

Solution

conda info --envs |grep fastp
fastp-0.23.1             /usr/local/genome/Anaconda3/envs/fastp-0.23.1

Étape 2: Identifier les paramètres obligatoires et les valeurs par défaut

Solution

conda activate fastp-0.23.1
fastp -h
#  -i, --in1                            read1 input file name (string [=])
#  -o, --out1                           read1 output file name (string [=])
#  -I, --in2                            read2 input file name (string [=])
#  -O, --out2                           read2 output file name (string [=])
#
#  -a, --adapter_sequence               if not specified, the adapter will be auto-detected.
#  -M, --cut_mean_quality               the mean quality requirement option. default: 20 (Q20)
#  -n, --n_base_limit                   if one read's number of N base is >n_base_limit, then this read/pair is discarded. Default is 5 (int [=5])
#  -l, --length_required                reads shorter than length_required will be discarded, default is 15. (int [=15])
#  ...

Étape 3: Soumettre le job et écrire les sorties dans un répertoire FASTP

Solution

mkdir FASTP
for i in *_R1.fastq.gz ; do id=$(echo $i | cut -d '_' -f 1) ; echo "conda activate fastp-0.23.1 && fastp --in1 $i --in2 ${id}_R2.fastq.gz --out1 FASTP/${id}_R1.fastq.gz --out2 FASTP/${id}_R2.fastq.gz --verbose --length_required 50 --html FASTP/${id}_fastp.html --json FASTP/${id}_fastp.json --thread 4 && conda deactivate" >> fastp.sh ; done
qsub -cwd -V -q web.q -N fastp -pe thread 4 -o LOGS/ -e LOGS/ fastp.sh

Vérifier l’effet du nettoyage

Solution

# avec seqkit
# avec multiqc
# en comptant le nombre de lignes dans les fichiers...

4.2 Séparer les reads d’ARN ribosomique

L’étape suivante consiste à filtrer et écarter les séquences de rRNA à partir de références, avec l’outil SortMeRNA [7]. Il effecute simplement un mapping des reads sur des bases de données disponibles dans /db/outils/sortmerna/

Utiliser SortMeRNA pour enlever les rRNA

Étape 1: Trouver l’environnement conda de la dernière version de SortMeRNA

Solution

conda info --envs |grep sortmerna
sortmerna-4.3.4          /usr/local/genome/Anaconda3/envs/sortmerna-4.3.4

Étape 2: Explorer les paramètres de SortMeRNA

Solution

conda activate sortmerna-4.3.4
sortmerna -h

#    --paired_in       BOOL        Optional  If one of the paired-end reads is Aligned,
#                                            put both reads into Aligned FASTA/Q file
#                                            Must be used with 'fastx'

Étape 3: Identifier les banques de données de rRNA

Solution

ll /db/outils/sortmerna/

total 74576
drwxr-xr-x 2 root root     8192 sept.  8  2020 idx
-rwxr-xr-x 1 root root  2280449 sept.  8  2020 rfam-5.8s-database-id98.fasta
-rwxr-xr-x 1 root root  8525326 sept.  8  2020 rfam-5s-database-id98.fasta
-rwxr-xr-x 1 root root  3893959 sept.  8  2020 silva-arc-16s-id95.fasta
-rwxr-xr-x 1 root root   752022 sept.  8  2020 silva-arc-23s-id98.fasta
-rwxr-xr-x 1 root root 19437013 sept.  8  2020 silva-bac-16s-id90.fasta
-rwxr-xr-x 1 root root 12911743 sept.  8  2020 silva-bac-23s-id98.fasta
-rwxr-xr-x 1 root root 13259584 sept.  8  2020 silva-euk-18s-id95.fasta
-rwxr-xr-x 1 root root 14945070 sept.  8  2020 silva-euk-28s-id98.fasta

Étape 4: Soumettre le job en utilisant toutes les banques disponibles

Solution

for i in FASTP/*_R1.fastq.gz ; do id=$(echo $(basename $i | cut -d '_' -f 1)) ; echo "conda activate sortmerna-4.3.4 && sortmerna --ref /db/outils/sortmerna/rfam-5.8s-database-id98.fasta --ref /db/outils/sortmerna/rfam-5s-database-id98.fasta --ref /db/outils/sortmerna/silva-arc-16s-id95.fasta --ref /db/outils/sortmerna/silva-arc-23s-id98.fasta --ref /db/outils/sortmerna/silva-bac-16s-id90.fasta --ref /db/outils/sortmerna/silva-bac-23s-id98.fasta --ref /db/outils/sortmerna/silva-euk-18s-id95.fasta --ref /db/outils/sortmerna/silva-euk-28s-id98.fasta --reads $i --reads FASTP/${id}_R2.fastq.gz --threads 4 --workdir SORTMERNA/${id} --fastx --out2 --aligned SORTMERNA/${id}.rRNA --other SORTMERNA/${id}.other -v --paired_in --num_alignments 1 --no-best --idx-dir /db/outils/sortmerna/idx/ && conda deactivate" >> sortmerna.sh ; done
qarray -cwd -V -q web.q -N sortmerna -pe thread 4 -o LOGS/ -e LOGS/ -R y sortmerna.sh

Étape 5: Trouver l’information dans les fichiers de sortie

Solution

tail  SORTMERNA/rRNA.log -n 20

 Results:
    Total reads passing E-value threshold = 22499 (0.54)
    Total reads failing E-value threshold = 4147445 (99.46)
    Minimum read length = 50
    Maximum read length = 151
    Mean read length    = 150

 Coverage by database:
    /db/outils/sortmerna/rfam-5.8s-database-id98.fasta      0.02
    /db/outils/sortmerna/rfam-5s-database-id98.fasta        0.02
    /db/outils/sortmerna/silva-arc-16s-id95.fasta       0.13
    /db/outils/sortmerna/silva-arc-23s-id98.fasta       0.17
    /db/outils/sortmerna/silva-bac-16s-id90.fasta       0.06
    /db/outils/sortmerna/silva-bac-23s-id98.fasta       0.13
    /db/outils/sortmerna/silva-euk-18s-id95.fasta       0.00
    /db/outils/sortmerna/silva-euk-28s-id98.fasta       0.00

 Mon Sep 14 16:12:50 2020

4.3 Enlever des contaminants par mapping

Dans certains projets, il peut être intéressant de filtrer les séquences d’un contaminant (hôte, contaminants connus…) Dans nos échantillons, il n’y a pas réellement d’hôte, mais voici par exemple la procédure pour écarter les séquences d’un des génomes viraux présents dans le jeu de données. Il s’agit de la séquence de Sorex araneus polyomavirus 1 isolate qui est présente dans le fichier /save_projet/metagenomics_training/contaminants/conta.fasta.

Supprimer les reads mappant sur la séquence présente dans le fichier conta.fasta

Étape 1: Copier le fichier /save_projet/metagenomics_training/contaminants/conta.fasta dans un nouveau répertoire DECONTA_MAPPING

Solution

mkdir DECONTA_MAPPING
cp /save_projet/metagenomics_training/contaminants/conta.fasta DECONTA_MAPPING/conta.fasta

Étape 2: Indexer le fichier FASTA avec bwa [8]

Solution

conda info --envs |grep bwa
bwa-0.7.12               /usr/local/genome/Anaconda3/envs/bwa-0.7.12
bwa-0.7.17               /usr/local/genome/Anaconda3/envs/bwa-0.7.17
bwa-mem2-2.2.1           /usr/local/genome/Anaconda3/envs/bwa-mem2-2.2.1

echo "conda activate bwa-0.7.17 && bwa index DECONTA_MAPPING/conta.fasta && conda deactivate" > bwa_index.sh
qsub -cwd -V -q web.q -N bwa_index -o LOG -e LOG bwa_index.sh

Étape 3: Mapper tous les reads sur cette banque

Solution

conda info --envs |grep samtools
msamtools-1.0.0          /usr/local/genome/Anaconda3/envs/msamtools-1.0.0
samtools-1.12            /usr/local/genome/Anaconda3/envs/samtools-1.12
samtools-1.14            /usr/local/genome/Anaconda3/envs/samtools-1.14
samtools-1.9             /usr/local/genome/Anaconda3/envs/samtools-1.9

for i in SORTMERNA/*.other_fwd.fq.gz ; do id=$(echo $(basename $i | cut -d '.' -f 1)) ; echo "conda activate samtools-1.14 && conda activate bwa-0.7.17 --stack && bwa mem -t 4 DECONTA_MAPPING/conta.fasta $i SORTMERNA/${id}.other_rev.fq.gz | samtools view -hbS - | samtools sort -@ 4 - > DECONTA_MAPPING/${id}.bam && samtools index DECONTA_MAPPING/${id}.bam && conda deactivate && conda deactivate" >> mapping_conta.sh ; done
qarray -cwd -V -q web.q -N bwa -o LOGS -e LOGS -pe thread 9 -R y mapping_conta.sh

Étape 4: Analyse de l’alignement avec samtools flagstat

Solution

for i in DECONTA_MAPPING/*.bam ; do conda activate samtools-1.14 && samtools flagstat $i > $i.flagstat && conda deactivate ; done

Étape 5: Extraction des reads ne mappant pas avec bam2fastq

Solution

conda info --envs |grep bam2fastq
bam2fastq-1.1.0          /usr/local/genome/Anaconda3/envs/bam2fastq-1.1.0

for i in DECONTA_MAPPING/*.bam ; do id=$(basename $i | cut -d '.' -f 1) ; echo "conda activate bam2fastq-1.1.0 && bam2fastq --no-aligned --unaligned --output DECONTA_MAPPING/${id}#.fastq $i && conda deactivate" >> bam2fastq.sh ; done
qarray -cwd -V -q web.q -N bam2fastq -o LOGS -e LOGS bam2fastq.sh

for i in DECONTA_MAPPING/*.fastq ; do gzip $i ; done

conda activate seqkit-2.0.0 && seqkit stat DECONTA_MAPPING/*.fastq.gz > DECONTA_MAPPING/seqkit.txt && conda deactivate

5 Assembler les reads

Maintenant que les reads sont nettoyés, nous pouvons passer à l’assemblage. Cette étape consiste en la construction de longs contigs à partir des reads.

5.1 Assemblage

Pour ce jeu de données nous avons choisi d’effectuer un assemblage poolé de tous nos échantillons. Dans certains cas, il peut être intéressant d’assembler les échantillons séparément. Nous vous proposons aujourd’hui d’utiliser MEGAHIT [9].

Assembler les reads avec megahit

Solution

echo "conda activate megahit-1.2.9 && megahit -1 DECONTA_MAPPING/s1_1.fastq.gz,DECONTA_MAPPING/s2_1.fastq.gz,DECONTA_MAPPING/s3_1.fastq.gz,DECONTA_MAPPING/s4_1.fastq.gz,DECONTA_MAPPING/s5_1.fastq.gz -2 DECONTA_MAPPING/s1_2.fastq.gz,DECONTA_MAPPING/s2_2.fastq.gz,DECONTA_MAPPING/s3_2.fastq.gz,DECONTA_MAPPING/s4_2.fastq.gz,DECONTA_MAPPING/s5_2.fastq.gz --num-cpu-threads 8 --memory 20 --tmp-dir /projet/tmp --out-dir MEGAHIT/" > megahit.sh
qsub -cwd -V -q web.q -N megahit -pe thread 8 -o LOGS/ -e LOGS/ megahit.sh

5.2 Évaluation de la qualité de l’assemblage

QUAST [10] permet d’obtenir et d’évaluer les métriques d’assemblage.

Lancer quast sur l’assemblage en ne prenant en compte que les contigs supérieus à 500 bp

Solution

echo "conda activate quast-5.0.2 && HOME=/projet/tmp/ && metaquast.py MEGAHIT/final.contigs.fa -o QUAST --min-contig 500 && conda deactivate" > quast.sh
qsub -cwd -V -q web.q -N quast -o LOGS/ -e LOGS/ quast.sh

Utiliser les références présentes dans le répertoire /save_projet/metagenomics_training/refs/ pour avoir des informations ciblées

Solution

echo "conda activate quast-5.0.2 && HOME=/projet/tmp/ && metaquast.py MEGAHIT/final.contigs.fa -o QUAST_REFS/ -r /save_projet/metagenomics_training/refs/ --min-contig 500 && conda deactivate" > quast_refs.sh
qsub -cwd -V -q web.q -N quast_refs -o LOGS/ -e LOGS/ quast_refs.sh

5.3 Mapping des reads sur les contigs

Maintenant, on aimerait replacer les reads sur les contigs pour obtenir les informations de couverture et de profondeur.

Obtenir les infos de couverture et profondeur à partir des alignements

Étape 1: Créer un répertoire MAPPING_ON_ASSEMBLY

Nous allons avoir besoin de bwa [8] et de samtools [11].

Solution

mkdir MAPPING_ON_ASSEMBLY

Étape 2: Indexer les contigs avec bwa [8]

Solution

echo "conda activate bwa-0.7.17 && bwa index MEGAHIT/final.contigs.fa && conda deactivate" > bwa_index_contigs.sh
qsub -cwd -V -q web.q -N bwa_index_contigs -o LOGS -e LOGS bwa_index_contigs.sh

Étape 3: Mapper les reads

Solution

for i in DECONTA_MAPPING/*_1.fastq.gz ; do id=$(echo $(basename $i | cut -d '_' -f 1)) ; echo "conda activate samtools-1.14 && conda activate bwa-0.7.17 --stack && bwa mem -t 4 MEGAHIT/final.contigs.fa $i DECONTA_MAPPING/${id}_2.fastq.gz | samtools view -hbS - | samtools sort -@ 4 - > MAPPING_ON_ASSEMBLY/${id}.bam && samtools index MAPPING_ON_ASSEMBLY/${id}.bam && conda deactivate" >> bwa_contigs.sh ; done
qarray -cwd -V -q web.q -N bwa_contigs -o LOGS -e LOGS -pe thread 8 -R y bwa_contigs.sh

Combien de reads sont mappés sur les contigs ? Utiliser l’outil samtools flagstat

Solution

for i in MAPPING_ON_ASSEMBLY/*.bam ; do conda activate samtools-1.14 && samtools flagstat $i > $i.flagstat && conda deactivate ; done

Utiliser samtools pour évaluer la profondeur et la couverture

Solution

for i in MAPPING_ON_ASSEMBLY/*.bam ; do conda activate samtools-1.14 && samtools coverage $i > $i.cov && conda deactivate ; done

6 Binning

L’assemblage de shorts reads issus de métagénomes shotguns permet rarement de reconstruire des génomes complets. Cependant le regroupement de contigs par binning permet de regrouper les séquences présentant des propriétés proches et est un bon compromis aux génomes complets. En effet, bien que fragmentés, ces ébauches de génomes sont souvent issues d’organismes proches. L’approche de binning qu’utilise MetaBAT2 [12] est basé sur la fréquence des tétranucléotides ainsi que les abondances des contigs en fonction des différents échantillons.

Lancer metabat2 après avoir exploré les paramètres

Solution

mkdir BINNING
echo "conda activate metabat2-2.15 && jgi_summarize_bam_contig_depths --outputDepth BINNING/depth.txt --pairedContigs BINNING/paired.txt --minContigLength 500 --minContigDepth 2 MAPPING_ON_ASSEMBLY/*.bam && metabat2 -i MEGAHIT/final.contigs.fa -a BINNING/depth.txt -o BINNING/bin -m 1500 -v -t 4 --unbinned && conda deactivate" > metabat2.sh
qsub -cwd -V -q web.q -N metabat -pe thread 8 -o LOGS/ -e LOGS/ metabat2.sh

6.1 Évaluation des bins

Pour l’évaluation des bins, nous utiliserons les deux métriques completeness et contamination estimés par CheckM [13]. Cet outil utilise des sets de gènes marqueurs colocalisés qui sont ubiquitaires et présent en une seule copie au sein d’un taxon phylogénétique donné. Parmi la suite d’outils fourni par CheckM nous allons utiliser le workflow checkm2 predict qui est recommandé pour évaluer l’exhaustivité et la contamination des bins de génomes.

Exécuter le workflow lineage-specific en spécifiant le dossier des bins

Solution

mkdir CHECKM/
echo "conda activate checkm2-0.1.3 && checkm2 predict --threads 8 --extension .fa --input BINNING/ --output-directory CHECKM/ && conda deactivate" > checkm.sh
qsub -cwd -V -q web.q -N checkm -pe thread 8 -o LOGS/ -e LOGS/ checkm.sh

Les résultats seront disponibles dans la table CHECKM/quality_report.tsv.

7 Prédiction de gènes

A partir de contigs au format FASTA, prokka [14] permet d’annoter rapidement les séquences bactériennes, archéennes et virales (ainsi que les métagénomes avec l’option --metagenome). Il est basé sur une suite d’outils : prodigal [15] (détéction des gènes), RNAmmer [16] (rRNA), aragorn [17] (tRNA), signalP (Signal leader peptides) [18], infernal (Non-coding RNA) [19].

Lancer prokka

Solution

echo "conda activate prokka-1.14.6 && prokka --outdir PROKKA --addgenes --kingdom Bacteria --metagenome --cpus 8 MEGAHIT/final.contigs.fa  && conda deactivate" > prokka.sh
qsub -cwd -V -q web.q -N prokka -o LOGS -e LOGS -pe thread 8 prokka.sh

Récupérer les comptages par gène

Créer un répertoire MAPPING_ON_GENES

Solution

mkdir MAPPING_ON_GENES

Récupérer uniquement les features gene dans le fichier d’annotation gff fourni par prokka

Solution

awk '$3 == "gene"' PROKKA/PROKKA_05162023.gff > MAPPING_ON_GENES/genes.gff

Utiliser le mapping sur les contigs pour en extraire avec htseq-count [20] les lecture s’alignant sur les gènes

Solution

echo "conda activate htseq-2.0.1 && htseq-count --with-header -f bam -t gene -i locus_tag MAPPING_ON_ASSEMBLY/*.bam MAPPING_ON_GENES/genes.gff > MAPPING_ON_GENES/genes.htseqcount && conda deactivate" > map_genes.sh

qsub -cwd -V -q web.q -N gene_count -o LOGS -e LOGS map_genes.sh

8 Annotation fonctionnelle sur les gènes

L’objectif de cette étape est d’annoter fonctionellement les gènes codant pour des protéines, à différente niveau de finesse et avec des ontologies ou des hiérarchies fonctionelles permettant l’aggrégation des annotations sur différents niveaux.

8.1 Annotation sur des cluster d’orthologues pré-calculés (egg-nog)

eggnog-mapper [21] est un outil d’annotation fonctionnelle. Il utilise les groupes de proteines orthologues et les phylogénies précalculées de la base eggNOG comme référnce pour transférer les annotations fonctionelles d’orthologues soigneusement identifiées.

Lancer eggnog-mapper

Solution

mkdir EGGNOG
echo "conda activate eggnog-mapper-2.1.5 &&  emapper.py -i PROKKA/PROKKA_05162023.faa --output eggnog -m diamond -d bact --cpu 8 --data_dir /db/outils/eggnog-mapper-2.1.5/ --output_dir EGGNOG && conda deactivate" > eggnog.sh
qsub -cwd -V -q web.q -N eggnog -o LOGS -e LOGS -pe thread 8 eggnog.sh

8.2 Annotation par similarité

Diamond [22] est un outil similaire à blast mais avec des heuristiques accélerant très fortement les comparaisons des gros jeux de données contre des banques.

Lancer diamond

Solution

mkdir DIAMOND
echo "conda activate diamond-2.0.14 && diamond blastp --db /db/outils/diamond-0.9.26/nr.dmnd --out DIAMOND/diamond.out --outfmt 6 --threads 8 --query PROKKA/PROKKA_05162023.faa && conda deactivate" > diamond.sh
qsub -cwd -V -q web.q -N diamond -o LOGS -e LOGS -pe thread 8 diamond.sh

Les références

1. Shen W, Le S, Li Y, Hu F. SeqKit: a cross-platform and ultrafast toolkit for FASTA/Q file manipulation. PloS one. 2016;11:e0163962.

2. Andrews S. FastQC A Quality Control tool for High Throughput Sequence Data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/. 2010. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.

3. Ewels P, Magnusson M, Lundin S, Käller M. MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics. 2016;32:3047‑8.

4. Menzel P, Ng KL, Krogh A. Fast and sensitive taxonomic classification for metagenomics with Kaiju. Nature communications. 2016;7:11257.

5. Ondov BD, Bergman NH, Phillippy AM. Interactive metagenomic visualization in a Web browser. BMC bioinformatics. 2011;12:385.

6. Zhou Y, Chen Y, Chen S, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34:i884‑90. doi:10.1093/bioinformatics/bty560.

7. Kopylova E, Noé L, Touzet H. SortMeRNA: fast and accurate filtering of ribosomal RNAs in metatranscriptomic data. Bioinformatics. 2012;28:3211‑7.

8. Li H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv preprint arXiv:13033997. 2013.

9. Li D, Liu C-M, Luo R, Sadakane K, Lam T-W. MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph. Bioinformatics. 2015;31:1674‑6.

10. Mikheenko A, Prjibelski A, Saveliev V, Antipov D, Gurevich A. Versatile genome assembly evaluation with QUAST-LG. Bioinformatics. 2018;34:i142‑50. doi:10.1093/bioinformatics/bty266.

11. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078‑9.

12. Kang DD, Froula J, Egan R, Wang Z. MetaBAT, an efficient tool for accurately reconstructing single genomes from complex microbial communities. PeerJ. 2015;3:e1165.

13. Parks DH, Imelfort M, Skennerton CT, Hugenholtz P, Tyson GW. CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Research. 2015;25:1043‑55. doi:10.1101/gr.186072.114.

14. Seemann T. Prokka: rapid prokaryotic genome annotation. Bioinformatics. 2014;30:2068‑9. doi:10.1093/bioinformatics/btu153.

15. Hyatt D, Chen G-L, LoCascio PF, Land ML, Larimer FW, Hauser LJ. Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC bioinformatics. 2010;11:1‑11.

16. Lagesen K, Hallin P, Rødland EA, Stærfeldt H-H, Rognes T, Ussery DW. RNAmmer: consistent and rapid annotation of ribosomal RNA genes. Nucleic acids research. 2007;35:3100‑8.

17. Laslett D, Canback B. ARAGORN, a program to detect tRNA genes and tmRNA genes in nucleotide sequences. Nucleic acids research. 2004;32:11‑6.

18. Petersen TN, Brunak S, Von Heijne G, Nielsen H. SignalP 4.0: discriminating signal peptides from transmembrane regions. Nature methods. 2011;8:785‑6.

19. Nawrocki EP, Kolbe DL, Eddy SR. Infernal 1.0: inference of RNA alignments. Bioinformatics. 2009;25:1335‑7.

20. Anders S, Pyl PT, Huber W. HTSeq—a Python framework to work with high-throughput sequencing data. bioinformatics. 2015;31:166‑9.

21. Cantalapiedra CP, Hernández-Plaza A, Letunic I, Bork P, Huerta-Cepas J. eggNOG-mapper v2: functional annotation, orthology assignments, and domain prediction at the metagenomic scale. Molecular biology and evolution. 2021;38:5825‑9.

22. Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nature Methods. 2014;12:59‑60. doi:10.1038/nmeth.3176.

Réutilisation

CC BY-SA 4.0

A work by Migale Bioinformatics Facility
Université Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France
Université Paris-Saclay, INRAE, BioinfOmics, MIGALE bioinformatics facility, 78350, Jouy-en-Josas, France