Shotgun metagenomics
Migale bioinformatics facility
Valentin Loux - Cédric Midoux - Olivier Rué
2022-06-21
library(data.table)
library(DT)
library(tidyverse)
options(DT.options = list(pageLength = 10,
scrollX = TRUE,
language = list(search = 'Filter:'),
dom = 'lftipB',
buttons = c('copy', 'csv', 'excel')))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/mkdir -p ~/work/TRAINING/CLUSTER
mkdir -p ~/work/TRAINING/ANALYSES
cd ~/work/TRAINING- Se déplacer dans le répertoire
CLUSTERet faire un lien symbolique des fichiers présents dans le répertoire/save_projet/metagenomics_training/cluster/
cd ~/work/TRAINING/CLUSTER
ln -s /save_projet/metagenomics_training/cluster/* .- Trouver l’environnement conda de blast
conda info --envs |grep blast
blast-2.12.0 /usr/local/genome/Anaconda3/envs/blast-2.12.0
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
qsub -cwd -V -q web.q -N myblast -pe thread 4 -b y "conda activate blast-2.12.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é ?
# 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 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.sh2 Analyse d’un jeu de données métagénomiques simulé
2.1 Présentation du jeu de données
Nous avons créé un jeu de données pour ce TP avec InsilicoSeq [1] avec les caractéristiques suivantes :
- 5 échantillons
- 1 M ou 2 M de reads par échantillon
- 50 Bacteria, 10 Archaea, 4 Viruses
- Modèle d’erreur Illumina Hiseq 2x125 bp
- Abondance uniforme ou log-normal
conda activate insilicoseq-1.5.4
# après un premier tirage aléatoire permettant de récupérer 50 génomes bactériens, 4 génomes viraux et 10 génomes d'archées
iss generate --seed 13062022 -k bacteria viruses archaea -U 50 4 10 --model hiseq --output s1 -n 2000000 --abundance_file hiseq_ncbi_abundance.txt
iss generate --seed 14062022 --genomes s1_genomes.fasta --model hiseq --output s1 -n 2000000 --abundance_file s1_abundance.txt -z
iss generate --seed 14062022 --genomes s2_genomes.fasta --model hiseq --output s2 -n 2000000 --abundance_file s2_abundance.txt -z
iss generate --seed 14062022 --genomes s1_genomes.fasta --model hiseq --abundance uniform --output s3 -n 2000000 -z
iss generate --seed 14062022 --genomes s2_genomes.fasta --model hiseq --abundance uniform --output s4 -n 1000000 -z
iss generate --seed 14062022 --genomes s1_genomes.fasta --model hiseq --output s5 -n 2000000 --abundance_file s1_abundance.txt -zVoici la liste des génomes récupérés dont les reads ont été extraits :
2.2 Préparation de l’espace de travail et organisation
.sh, nommés en fonction de l’outil à lancer.
- Se déplacer dans le répertoire
ANALYSESpuis faire un lien symbolique des fichiers présents dans le répertoire/save_projet/metagenomics_training/simulated_reads/
cd ~/work/TRAINING/ANALYSES
ln -s /save_projet/metagenomics_training/simulated_reads/* .- Créer un répertoire
LOGdans lequel seront stockés les fichiers spécifiques à la soumission des jobs sur le cluster de calcul
mkdir ~/work/TRAINING/ANALYSES/LOG3 Analyse préliminaire des reads bruts
3.1 Inspection des fichiers
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
3.1.1 Nombre de reads
- Combien de reads sont présents dans les fichiers
s1_R1.fastq.gzets1_R2.fastq.gz?
Étape 1: Utiliser la commande zcat pour lire un
fichier compressé et wc -l pour compter le nombre de
lignes
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
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
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
# ...3.1.2 Contenu des fichiers
- Vérifier que les fichiers FASTQ sont conformes et obtenir des statistiques sur leur contenu avec l’outil seqkit [2]
Étape 1: Trouver l’environnement conda de seqkit
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
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
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 LOG -e LOG seqkit.shÉtape 4: Lire le fichier de sortie
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
3.2 Contrôle Qualité des reads
3.2.1 FastQC
FastQC [3] 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
conda info --envs |grep fastqc
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
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
qarray -cwd -V -q web.q -N fastqc -pe thread 4 -o LOG/ -e LOG/ fastqc.sh3.2.2 MultiQC
MultiQC [4] 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
conda info --envs |grep multiqc
multiqc-1.11 /usr/local/genome/Anaconda3/envs/multiqc-1.11
mkdir MULTIQC
echo "conda activate multiqc-1.11 && multiqc --outdir MULTIQC/ FASTQC/ && conda deactivate" > multiqc.sh
qsub -cwd -V -q web.q -N multiqc -o LOG -e LOG multiqc.sh- Visualiser le rapport
firefox --no-remote MULTIQC/multiqc_report.html &3.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 [5].
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
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É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/.
mkdir KAIJU
for i in *_R1.fastq.gz ; do id=$(echo $(basename $i | cut -d '_' -f 1)) ; echo "conda activate kaiju-1.8.1 && 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 LOG -e LOG -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/simulated_reads/update_taxid_column_from_merged_ncbi_file.py.
cp /save_projet/metagenomics_training/simulated_reads/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 LOG -e LOG -N updatetaxonomy update_taxo.shÉtape 4: Générer le fichier Krona avec la suite Krona-tools [6]
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-2.8/taxonomy -t 3 -q 2 $i -o $i.html && conda deactivate" >> krona.sh ; done
qarray -cwd -V -q web.q -N krona -o LOG -e LOG 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 LOG -e LOG krona_single.shÉtape 5: Visualiser le fichier HTML généré
firefox --no-remote KAIJU/kaiju.cor.html &
firefox --no-remote KAIJU/all.html &4 Nettoyage
4.1 Filtrer et trimmer les reads
FASTP [7] 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 adapters restants.
- Explorez les paramètres de FASTP et lancez FASTP sur les données brutes
Étape 1: Trouver l’environnement conda de fastp
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
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
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 LOG/ -e LOG/ fastp.sh*.json peuvent être interprétés par
MULTIQC.
- Vérifier l’effet du nettoyage
# 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 [8]. 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
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
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
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
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 LOG/ -e LOG/ -R y sortmerna.shÉtape 5: Trouver l’information dans les fichiers de sortie
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 20204.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
mkdir DECONTA_MAPPING
cp /save_projet/metagenomics_training/contaminants/conta.fasta DECONTA_MAPPING/conta.fastaÉtape 2: Indexer le fichier FASTA avec bwa [9]
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
conda info --envs |grep samtools
msamtools-1.0.0 /usr/local/genome/Anaconda3/envs/msamtools-1.0.0
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 LOG -e LOG -pe thread 9 -R y mapping_conta.shÉtape 4: Analyse de l’alignement avec
samtools flagstat
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
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 LOG -e LOG 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 deactivate5 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 [10].
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 LOG/ -e LOG/ megahit.sh5.2 Évaluation de la qualité de l’assemblage
QUAST [11] 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
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 LOG/ -e LOG/ quast.sh- Utiliser les références présentes dans le répertoire
/save_projet/metagenomics_training/refs/pour avoir des informations ciblées
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 LOG/ -e LOG/ quast_refs.sh5.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.
Étape 1: Créer un répertoire
MAPPING_ON_ASSEMBLY
Nous allons avoir besoin de bwa [9] et de samtools [12].
mkdir MAPPING_ON_ASSEMBLYÉtape 2: Indexer les contigs avec bwa [9]
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 LOG -e LOG bwa_index_contigs.shÉtape 3: Mapper les reads
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 LOG -e LOG -pe thread 8 -R y bwa_contigs.sh- Combien de reads sont mappés sur les contigs ? Utiliser
l’outil
samtools flagstat
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
for i in MAPPING_ON_ASSEMBLY/*.bam ; do conda activate samtools-1.14 && samtools coverage $i > $i.cov && conda deactivate ; done6 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’utiliseMetaBAT2 [13] 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
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 LOG/ -e LOG/ metabat2.sh6.1 Évaluation des bins
Pour l’évaluation des bins, nous utiliserons les deux métriques
completeness et contamination estimés par
CheckM [14].
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
lineage-specific qui est recommandé pour évaluer l’exhaustivité
et la contamination des bins de génomes.
- Exécuter le workflow
lineage-specificen spécifiant le dossier des bins
mkdir CHECKM/
echo "conda activate checkm-genome-1.1.3 && checkm lineage_wf --extension fa --threads 8 --file CHECKM/report.tsv --tab_table BINNING/ CHECKM/ && conda deactivate" > checkm.sh
qsub -cwd -V -q web.q -N checkm -pe thread 8 -o LOG/ -e LOG/ checkm.sh- Explorer plus en détail les résultats avec
checkm qa
conda activate checkm-genome-1.1.3
checkm qa CHECKM/lineage.ms CHECKM/ -o2
checkm qa CHECKM/lineage.ms CHECKM/ -o6
conda deactivate7 Prédiction de gènes
A partir de contigs au format FASTA,
prokka [15]
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 [16] (détéction des gènes),
RNAmmer (rRNA) [17], aragorn
(tRNA) [18],
signalP (Signal leader peptides) [19], infernal
(Non-coding RNA) [20].
- Lancer prokka
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 LOG -e LOG -pe thread 8 prokka.sh7.1 Récupérer les comptages par gène
- Créer un répertoire dédié
mkdir MAPPING_ON_GENES- Récupérer uniquement les features
genedans le fichier d’annotation gff fourni par prokka
awk '$3 == "gene"' PROKKA/PROKKA_06152022.gff > MAPPING_ON_GENES/genes.gff- Utiliser le mapping sur les contigs pour en extraire avec htseq-count [21] les lecture s’alignant sur les gènes
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 LOG -e LOG map_genes.sh8 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 [22] 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.
mkdir EGGNOG
echo "conda activate eggnog-mapper-2.1.5 && emapper.py -i PROKKA/PROKKA_06152022.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 LOG -e LOG -pe thread 8 eggnog.sh8.2 Annotation par similarité
Diamond [23] 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.
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_06152022.faa && conda deactivate" > diamond.sh
qsub -cwd -V -q web.q -N diamond -o LOG -e LOG -pe thread 8 diamond.sh9 References
A work by Migale Bioinformatics Facility
https://migale.inrae.fr
Our two affiliations to cite us:
Université Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France
Université Paris-Saclay, INRAE, BioinfOmics, MIGALE bioinformatics facility, 78350, Jouy-en-Josas, France
