Rappels sur les ressources de la plateforme Migale

Lancer un job

Créer cette arborescence votre espace work

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

WORK_DIR="TRAINING"
mkdir -p ~/work/$WORK_DIR/CLUSTER
mkdir -p ~/work/$WORK_DIR/ANALYSES
cd ~/work/$WORK_DIR

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

cd ~/work/$WORK_DIR/CLUSTER
ln -s /db/tutorials/cluster/* .

Trouver l’environnement conda de blast

conda info --envs |grep blast
blast                    /usr/local/genome/Anaconda2-5.1.0/envs/blast
blast-2.9.0              /usr/local/genome/Anaconda2-5.1.0/envs/blast-2.9.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 formation.q -N myblast -pe thread 4 -b y "conda activate blast-2.9.0 && blastn -db subject.fasta -query query.fasta -out out.blast -num_threads 4 && conda deactivate"

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 formation.q -N myblast -pe thread 4 -b y "conda activate blast-2.9.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.9.0 && blastn -db subject.fasta -query $i -out $i.blast -num_threads 4 && conda deactivate" >> myblast.sh ; done
qarray -cwd -V -q formation.q -N myblast -pe thread 4 myblast.sh

Analyse d’un jeu de données métagénomiques

Jeu de donnée utilisé dans la publi :

Shotgun metagenome data of a defined mock community using Oxford Nanopore, PacBio and Illumina Technologies

Accession number : SRX4901583 Dans la publi, pour les données Illumina, les caractéristiques du séquençage sont : Illumina HiSeq 2500 2x150 bp : 426,735,646 reads bruts.

REDUCED_DATA/
├── R1.fastq.gz
├── R2.fastq.gz
REFS/
├── GCF_002563805.1_ASM256380v1_genomic.fna.gz
├── GCF_002563885.1_ASM256388v1_genomic.fna.gz
├── GCF_002813125.1_ASM281312v1_genomic.fna.gz
├── GCF_002846515.1_ASM284651v1_genomic.fna.gz
├── GCF_004379815.1_ASM437981v1_genomic.fna.gz
├── GCF_900086985.1_IMG-taxon_2623620618_annotated_assembly_genomic.fna.gz
├── GCF_900090235.1_IMG-taxon_2623620557_annotated_assembly_genomic.fna.gz
├── GCF_900090295.1_IMG-taxon_2623620609_annotated_assembly_genomic.fna.gz
├── GCF_900091445.1_IMG-taxon_2623620567_annotated_assembly_genomic.fna.gz
├── GCF_900215475.1_IMG-taxon_2623620617_annotated_assembly_genomic.fna.gz
├── GCF_900215505.1_IMG-taxon_2615840601_annotated_assembly_genomic.fna.gz
└── Psychrobacter_sp_LV10R520-6.fasta.gz

Récupération des fichiers

Le fichier au format SRA a été téléchargé puis les fichiers FASTQ ont été générés grâce à l’outil sra-toolkit “SRA Tools” (2020). Voici les commandes qui ont permis la récupération des fichiers :

cd /projet/metagenomics_training/FULL
conda activate sra-tools-2.10.1
prefetch SRR8073716 --max-size 200G
fastq-dump SRR8073716.sra --split-3 --skip-technical
conda deactivate
qsub -pe thread 8 -R y -b y "pigz -p 8 *.fastq"

Les génomes de référence correspondant aux génomes utilisés dans la Mock ont été téléchargés grâce à l’outil ncbi-genome-download.

conda activate ncbi-genome-download
ncbi-genome-download --taxids "1798204,1798206,1798205,1945889,1415571,1415567,1415574,47857,47858,356852,1761789,1666906" bacteria --formats all -o refs
# Psychrobacter not available at NCBI, downloaded from JGI
conda deactivate

Voici la liste des génomes utilisés :

IMG Taxon ID	Organism	Phylum	Class
2615840527	Muricauda sp. ES.050	Bacteroidetes	Flavobacteria
2615840533	Thioclava sp. ES.032	Proteobacteria	Alphaproteobacteria
2615840601	Cohaesibacter sp. ES.047	Proteobacteria	Alphaproteobacteria
2615840646	Propionibacteriaceae bacterium ES.041	Actinobacteria	Actinobacteria
2615840697	Marinobacter sp. LV10R510-8	Proteobacteria	Gammaproteobacteria
2616644829	Marinobacter sp. LV10MA510-1	Proteobacteria	Gammaproteobacteria
2617270709	Psychrobacter sp. LV10R520-6	Proteobacteria	Gammaproteobacteria
2623620557	Micromonospora echinaurantiaca DSM 43904	Actinobacteria	Actinobacteria
2623620567	Micromonospora echinofusca DSM 43913	Actinobacteria	Actinobacteria
2623620609	Micromonospora coxensis DSM 45161	Actinobacteria	Actinobacteria
2623620617	Halomonas sp. HL-4	Proteobacteria	Gammaproteobacteria
2623620618	Halomonas sp. HL-93	Proteobacteria	Gammaproteobacteria

Pour accélérer le temps de calcul lors de cette formation, nous avons extrait un sous-ensemble de reads du jeu de données, correspondant à 1% du nombre total de reads.

conda activate seqtk-1.3
seqtk sample -s100 SRR8073716_1.fastq.gz 0.01 | pigz - > R1.fastq.gz
seqtk sample -s100 SRR8073716_2.fastq.gz 0.01 | pigz - > R2.fastq.gz
conda deactivate

Préparation de l’espace de travail

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

cd ~/work/$WORK_DIR/ANALYSES
ln -s /projet/metagenomics_training/REDUCED_DATA/* .

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

mkdir ~/work/$WORK_DIR/ANALYSES/LOG

Pour une question de traçabilité, les commandes que nous enverrons sur le cluster seront stockées dans des fichiers .sh, nommés en fonction de l’outil à lancer. Ainsi, on pourra retracer facilement ce qui a été lancé.

Contrôle qualité des fichiers

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, à savoir :

autant de fichiers que d’échantillons
un nombre de reads suffisants
des tailles de reads conformes
une cohérence entre les informations contenues dans chaque fichier FASTQ pour un séquençage paired-end

Nombre de reads

Combien de reads sont présents dans les fichiers R1.fastq.gz et R2.fastq.gz ?

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

zcat 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 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
# R1.fastq.gz 2131034
# R2.fastq.gz 2131034

Sur le jeu de données réel

for i in /projet/metagenomics_training/FULL/*.fastq.gz ; do echo $i $(zcat SRR8073716_2.fastq.gz | echo $((`wc -l`/4))) ; done
SRR8073716_1.fastq.gz 213367823
SRR8073716_2.fastq.gz 213367823

Contenu des fichiers

Vérifiez que les fichiers FASTQ sont conformes et obtenez des statistiques sur leur contenu avec l’outil seqkit (Shen et al. (2016)).

Étape 1: Trouver l’environnement conda de seqkit

conda info --envs |grep seqkit

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

echo "conda activate seqkit-0.11.0 && seqkit stat *.fastq.gz > seqkit.txt && conda deactivate" > seqkit.sh
qsub -cwd -V -q formation.q -N seqkit -o LOG -e LOG seqkit.sh
cat seqkit.txt

#file                   format  type     num_seqs         sum_len  min_len  avg_len  max_len
#R1.fastq.gz           FASTQ   DNA   2,131,034  321,786,134      151      151      151
#R2.fastq.gz           FASTQ   DNA   2,131,034  321,786,134      151      151      151

Sur le jeu de données réel

conda activate seqkit-0.11.0
qsub -cwd -V -q formation.q -N seqkit -b y "seqkit stat /projet/metagenomics_training/FULL/*.fastq.gz > seqkit.txt"
file                   format  type     num_seqs         sum_len  min_len  avg_len  max_len
SRR8073716_1.fastq.gz  FASTQ   DNA   213,367,823  32,218,541,273      151      151      151
SRR8073716_2.fastq.gz  FASTQ   DNA   213,367,823  32,218,541,273      151      151      151
conda deactivate

Analyse préliminaire des reads bruts

Contrôle Qualité des reads

FastQC

FastQC Andrews (2010) est utilisé classiquement 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 échantillon par échantillon, en parallèle

Étape 1: Trouver l’environnement conda de fastqc

conda info --envs |grep fastqc

É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.8 && 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 formation.q -N fastqc -pe thread 4 -o LOG/ -e LOG/ fastqc.sh

Sur le jeu de données réel

for i in $FULL/*.fastq.gz ; do echo "conda activate fastqc-0.11.8 && fastqc --outdir $FULL/FASTQC/ --threads 4 $i && conda deactivate" >> fastqc.sh ; done
qarray -cwd -V -q formation.q -N fastqc -pe thread 4 -o LOG/ -e LOG/ fastqc.sh

Les fichiers suivants sont les résultats de FastQC sur le jeu de données complet :

MultiQC

MultiQC Ewels et al. (2016) est un outil qui permet de rassembler les rapports de plusieurs échantillons et outils en un seul rapport. Il est particulièrement utile quand on travaille sur plusieurs échantillons.

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

mkdir MULTIQC
echo "conda activate multiqc-1.8 && multiqc --outdir MULTIQC/ FASTQC/ && conda deactivate" > multiqc.sh
qsub -cwd -V -q formation.q -N multiqc -o LOG -e LOG multiqc.sh

Visualiser le rapport

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

Sur le jeu de données réel

qsub -cwd -V -q formation.q -N multiqc -o LOG -e LOG -b y "conda activate multiqc-1.8 && multiqc --outdir $FULL/MULTIQC/ $FULL/FASTQC/ && conda deactivate"

Le fichier suivant est le résultat de MultiQC sur le jeu de données complet :

rapport web multiQC

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 nr et l’outil Kaiju Menzel, Ng, and Krogh (2016).

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

Étape 1: Trouver l’environnement conda de kaiju

conda info --envs |grep kaiju

É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/kaiju/.

mkdir KAIJU
echo "conda activate kaiju-1.7.3 && kaiju -t /db/outils/kaiju-2020-05-25/refseq/nodes.dmp -f /db/outils/kaiju-2020-05-25/refseq/kaiju_db_refseq.fmi -i R1.fastq.gz -j R2.fastq.gz -o KAIJU/kaiju -z 4 && conda deactivate" > kaiju.sh
qsub -cwd -V -q formation.q -N kaiju -o LOG -e LOG -pe thread 4 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 update_taxid_column_from_merged_ncbi_file.py.

echo "update_taxid_column_from_merged_ncbi_file.py --dump /db/ncbi_taxon/current/flat/merged.dmp --file KAIJU/kaiju --taxid-col 3 > KAIJU/kaiju.cor" > update_taxo.sh
qsub -cwd -V -q formation.q -o LOG -e LOG -N updatetaxonomy update_taxo.sh

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

echo "conda activate krona-2.71 && ktImportTaxonomy -tax /db/outils/krona-2.7 -t 3 -q 2 KAIJU/kaiju.cor -o KAIJU/kaiju.cor.html  && conda deactivate" > krona.sh
qsub -cwd -V -q formation.q -N krona -o LOG -e LOG krona.sh

Visualiser le fichier HTML généré

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

Bonus: Assigner les reads avec Kraken

# Kraken vs bacteria ~ 40 Go mem
conda activate kraken2-2.0.8_beta
kraken2 --db /db/outils/kraken-2.0/bacteria/ --threads 8 --quick --classified-out /projet/tmp/kraken-classified --unclassified-out /projet/tmp/kraken-unclassified --gzip-compressed /projet/metagenomics_training/TP/R1.fastq.gz /projet/metagenomics_training/TP/R2.fastq.gz --output /projet/tmp/kraken --report /projet/tmp/kraken-report
conda deactivate
conda activate krona-2.71
ktImportTaxonomy -q 2 -t 3 /projet/tmp/kraken -o /projet/tmp/kraken.html
conda deactivate

Sur le jeu de données réel, lancé contre nr

conda activate kaiju-1.7.3
kaiju -t /db/outils/kaiju/nr/nodes.dmp -f /db/outils/kaiju/nr/kaiju_db_nr_euk.fmi -i $FULL/SRR8073716_1.fastq.gz -j $FULL/SRR8073716_2.fastq.gz -o $FULL/KAIJU_RAW/kaijuNR -z 20
conda deactivate
update_taxid_column_from_merged_ncbi_file.py --dump /db/ncbi_taxon/current/flat/merged.dmp --file $FULL/KAIJU_RAW/kaijuNR --taxid-col 3 > $FULL/KAIJU_RAW/kaijuNR.cor
conda activate krona-2.71
ktImportTaxonomy -tax /db/outils/krona-2.7 -t 3 -q 2 $FULL/KAIJU_RAW/kaijuNR.cor -o $FULL/KAIJU_RAW/kaijuNR.cor.html
conda deactivate

Voici le Krona sur le jeu de données réel:

Nettoyage

Filtrer et trimmer les reads

Nous avons vu avec FASTQC Andrews (2010) qu’il existait une variabilité de qualité des données.
FASTP Zhou et al. (2018) est un outils de préprocessing tout en un. Il va nous permettre d’éliminer les séquences de mauvaise qualité, les bases ambiguës, 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

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

mkdir FASTP
echo "conda activate fastp-0.20.0 && fastp --in1 R1.fastq.gz --in2 R2.fastq.gz --out1 FASTP/R1.fastq.gz --out2 FASTP/R2.fastq.gz --verbose --length_required 50 --html FASTP/fastp.html --json FASTP/fastp.json --thread 4 && conda deactivate" > fastp.sh
qsub -cwd -V -q formation.q -N fastp -pe thread 4 -o LOG/ -e LOG/ fastp.sh

Les données nettoyées sont placées dans le dossier FASTP/ avec un rapport d’analyse.
Voici le rapport FASTP sur le jeu de données complet :

FASTP

Le rapport JSON fastp.json peut être agrégé avec MULTIQC.

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 Kopylova, Noé, and Touzet (2012).

Filtrez les rRNA parmi les données nettoyées avec FASTP.
Les bases de données de rRNA sont disponibles dans /db/outils/sortmerna/.
En cas d’ambiguité entre les deux reads d’un paire, considérez la paire comme rRNA.

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

conda info --envs |grep sortmerna

Étape 2: Explorer les paramètres de sortmerna-4.3.2

conda activate sortmerna-4.3.2
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/

# drwxr-xr-x 2 root root 8.0K Sep  8 15:53 idx
# -rwxr-xr-x 1 root root 2.2M Sep  8 15:47 rfam-5.8s-database-id98.fasta
# -rwxr-xr-x 1 root root 8.2M Sep  8 15:47 rfam-5s-database-id98.fasta
# -rwxr-xr-x 1 root root 3.8M Sep  8 15:47 silva-arc-16s-id95.fasta
# -rwxr-xr-x 1 root root 735K Sep  8 15:47 silva-arc-23s-id98.fasta
# -rwxr-xr-x 1 root root  19M Sep  8 15:47 silva-bac-16s-id90.fasta
# -rwxr-xr-x 1 root root  13M Sep  8 15:47 silva-bac-23s-id98.fasta
# -rwxr-xr-x 1 root root  13M Sep  8 15:47 silva-euk-18s-id95.fasta
# -rwxr-xr-x 1 root root  15M Sep  8 15:47 silva-euk-28s-id98.fasta

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

mkdir SORTMERNA
echo "conda activate sortmerna-4.3.2 && 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 FASTP/R1.fastq.gz --reads FASTP/R2.fastq.gz --threads 4 --workdir SORTMERNA --fastx --out2 --aligned SORTMERNA/rRNA --other SORTMERNA/other -v --paired_in --num_alignments 1 --idx-dir /db/outils/sortmerna/idx/ && conda deactivate" > sortmerna.sh
qsub -cwd -V -q formation.q -N sortmerna -q maiage.q -pe thread 4 -o LOG/ -e LOG/ -R y sortmerna.sh

Étape 5: Résultats

Vous pouvez voir la répartition des reads mappés contre les différentes banques dans le fichier log

 Results:
    Total reads passing E-value threshold = 22499 (0.54)
    Total reads failing E-value threshold = 4147445 (99.46)
    Minimum read length = 151
    Maximum read length = 50
    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

Étape 6: Renommer les fichiers de sortie de sortmerna

cd SORTMERNA
mv other_fwd.fq.gz other_R1.fastq.gz
mv other_rev.fq.gz other_R2.fastq.gz
cd ../

Résultats sur le jeu de données complet

Results:
    Total reads passing E-value threshold = 2204810 (0.53)
    Total reads failing E-value threshold = 415297800 (99.47)
    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.12
    /db/outils/sortmerna/silva-arc-23s-id98.fasta       0.15
    /db/outils/sortmerna/silva-bac-16s-id90.fasta       0.06
    /db/outils/sortmerna/silva-bac-23s-id98.fasta       0.15
    /db/outils/sortmerna/silva-euk-18s-id95.fasta       0.00
    /db/outils/sortmerna/silva-euk-28s-id98.fasta       0.00

Enlever des contaminants par mapping

Dans certains projets, il peut être intéressant de filtrer les séquences d’un contaminant, comme par exemple celles de l’hôte. Dans ce jeu de données, il n’y a pas d’hôte, mais voici par exemple la procédure pour écarter les séquences de Muricauda sp. (ref : GCF_002813125.1).

Étape 1: Faire un lien symbolique du fichier /projet/metagenomics_training/REFS/GCF_002813125.1_ASM281312v1_genomic.fna dans un nouveau répertoire MAPPING

mkdir MAPPING
cd MAPPING
ln -s /projet/metagenomics_training/REFS/GCF_002813125.1_ASM281312v1_genomic.fna
cd ..

Étape 2: Indexer le génome avec bwa (H. Li 2013)

echo "conda activate bwa-0.7.17 && bwa index MAPPING/GCF_002813125.1_ASM281312v1_genomic.fna && conda deactivate" > bwa_index.sh
qsub -cwd -V -q formation.q -N bwa_index -o LOG -e LOG bwa_index.sh

Étape 3: Mapping

echo "conda activate samtools-1.12 && conda activate bwa-0.7.17 --stack && conda activate bedtools-2.29.0 --stack && bwa mem -t 4 MAPPING/GCF_002813125.1_ASM281312v1_genomic.fna SORTMERNA/other_R1.fastq.gz SORTMERNA/other_R2.fastq.gz | samtools view -hbS -f 12 - | samtools sort -n - -@4 | bedtools bamtofastq -i - -fq MAPPING/no_muricauda_R1.fastq -fq2 MAPPING/no_muricauda_R2.fastq ; gzip MAPPING/no_muricauda_R1.fastq ; gzip MAPPING/no_muricauda_R2.fastq ;" > bwa.sh

qsub -cwd -V -q formation.q -N bwa -o LOG -e LOG -pe thread 10 -R y bwa.sh

Normalisation avant assemblage (digital normalization)

Lorsque la profondeur de séquençage est importante, il peut être intéressant de supprimer les reads sur-représentés pour faciliter l’assemblage et gagner en temps de calcul.

Étape 1: Pour cette étape de normalisation, nous utilisons normalize-by-median.py de KHMER. Cependant, cet outil nécessite des reads enchevêtrés (interleaved). On commence donc par cette étape :

mkdir KHMER
#paste <(zcat MAPPING/no_muricauda_R1.fastq.gz) <(zcat MAPPING/no_muricauda_R2.fastq.gz) | paste - - - - | awk -v OFS="\n" -v FS="\t" '{print($1."/1",$3,$5,$7,$2."/2",$4,$6,$8)}' |gzip - > KHMER/no_muricauda.fastq.gz

paste <(zcat MAPPING/no_muricauda_R1.fastq.gz) <(zcat MAPPING/no_muricauda_R2.fastq.gz) | paste - - - - | awk -v OFS="\n" -v FS="\t" '{print($1,$3,$5,$7,$2,$4,$6,$8)}' |gzip - > KHMER/no_muricauda.fastq.gz

ou (moins rapide)

echo "conda activate khmer-3.0.0a3 && interleave-reads.py --output KHMER/R1R2.fastq.gz --gzip MAPPING/no_muricauda_R1.fastq.gz MAPPING/no_muricauda_R2.fastq.gz && conda deactivate" > interleave.sh
qsub -cwd -V -q formation.q -N interleave -o LOG/ -e LOG/ interleave.sh

Étape 2: On peut normaliser en écartant les k-mers sur-représentés (k = 32bp et cutoff à 20 occurences)

echo "conda activate khmer-3.0.0a3 && normalize-by-median.py KHMER/no_muricauda.fastq.gz -o KHMER/no_muricauda_normalized_R1R2.fastq.gz -M 16G --gzip --report KHMER/khmer.txt && conda deactivate" > normalize.sh
qsub -cwd -V -q formation.q -N normalize -o LOG/ -e LOG/ normalize.sh

Étape 3: Résultats

cat KHMER/khmer.txt 
# total,kept,f_kept
# 100000,99998,1.0
# ...
# 3400742,2770098,0.8146

conda activate seqkit-0.10.1
seqkit stat KHMER/no_muricauda_normalized_R1R2.fastq.gz
# file                                         format  type   num_seqs      sum_len  min_len  avg_len  max_len
# KHMER/no_muricauda_normalized_R1R2.fastq.gz  FASTQ   DNA   2,770,098  417,985,358       50    150.9      151

conda deactivate

Assembler les reads

Assemblage

Maintenant que les reads sont nettoyés, nous pouvons passer à l’étape centrale de l’analyse métagénomique : l’assemblage. Cette étape consiste en la reconstruction de longs contigs, voir de génome complet, à partir de nos reads.

Spades

SPADES Bankevich et al. (2012) est l’un des assembleur les plus populaire. Il est primordial d’utiliser l’option --meta pour les projets de métagénomiques.

Étape 1: Séparer les reads enchevérés pour revenir à une situation classique R1.fastq.gz et R2.fastq.gz

zcat KHMER/no_muricauda_normalized_R1R2.fastq.gz | paste - - - - - - - - | tee | cut -f 1-4 | tr "\t" "\n" | egrep -v '^$' | gzip - > KHMER/no_muricauda_R1.fastq.gz
zcat KHMER/no_muricauda_normalized_R1R2.fastq.gz | paste - - - - - - - - | tee | cut -f 5-8 | tr "\t" "\n" | egrep -v '^$' | gzip - > KHMER/no_muricauda_R2.fastq.gz

Étape 2: Trouver l’environnement conda de spades

conda info --envs |grep spades

Étape 3: Soumettre le job

mkdir SPADES
echo "conda activate spades-3.14.0 && spades.py --threads 8 --memory 10 --tmp-dir /projet/tmp --meta -1 KHMER/no_muricauda_R1.fastq.gz -2 KHMER/no_muricauda_R2.fastq.gz -o SPADES" > spades.sh
qsub -cwd -V -q formation.q -N spades -pe thread 8 -o LOG/ -e LOG/ spades.sh

Megahit

En supplément, voici la commande pour un assemblage avec MEGAHIT D. Li et al. (2015)

echo "conda activate megahit-1.2.9 && megahit -1 KHMER/no_muricauda_R1.fastq.gz -2 KHMER/no_muricauda_R2.fastq.gz --preset meta-large --num-cpu-threads 8 --memory 20 --tmp-dir /projet/tmp --out-dir MEGAHIT/" > megahit.sh
qsub -cwd -V -q formation.q -N megahit -pe thread 8 -o LOG/ -e LOG/ megahit.sh

Évaluation de la qualité de l’assemblage

On utilse QUAST Mikheenko et al. (2018) pour évaluer les métriques et la qualité de l’assemblage.

Sans référence

echo "HOME=/projet/tmp/$USER ; conda activate quast-5.0.2 && metaquast.py SPADES/contigs.fasta -o QUAST --min-contig 1000 --gene-finding --rna-finding --conserved-genes-finding && conda deactivate" > quast.sh

echo "HOME=/projet/tmp/$USER ; conda activate quast-5.0.2 && metaquast.py --threads 8 SPADES/contigs.fasta -o QUAST --min-contig 1000  && conda deactivate" > quast.sh
qsub -cwd -V -q formation.q -pe thread 8 -N quast -o LOG/ -e LOG/ quast.sh

Avec références

mkdir REFS
cd REFS
ln -s /projet/metagenomics_training/REFS/12REFS/* .
cd ../
#echo "HOME="\$WORK_DIR" ; conda activate quast-5.0.2 && metaquast.py SPADES/contigs.fasta -o QUAST_REFS/ -r REFS/ --min-contig 1000 --gene-finding --rna-finding --conserved-genes-finding && conda deactivate" > quast_refs.sh
echo "HOME=/projet/tmp/$USER ; conda activate quast-5.0.2 && metaquast.py --threads 8 SPADES/contigs.fasta -o QUAST_REFS/ -r REFS/ --min-contig 1000  && conda deactivate" > quast_refs.sh
qsub -cwd -V -q formation.q -pe thread 8 -N quast_refs -o LOG/ -e LOG/ quast_refs.sh

Résultats sur le jeu de données complet sans référence et avec références.

Mapping des reads sur les contigs

Étape 1: Créer un répertoire MAPPING_ON_ASSEMBLY et faire un lien symbolique des contigs issus de SPADES

mkdir MAPPING_ON_ASSEMBLY
cd MAPPING_ON_ASSEMBLY
ln -s ../SPADES/contigs.fasta .

Étape 2: Indexer les contigs avec bwa

echo "conda activate bwa-0.7.17 && bwa index contigs.fasta && conda deactivate" > bwa_index_contigs.sh
qsub -cwd -V -q formation.q -N bwa_index_contigs -o ../LOG -e ../LOG bwa_index_contigs.sh

Étape 3: Mapping

cd ../
echo "conda activate samtools && conda activate bwa-0.7.17 --stack && bwa mem -t 4 MAPPING_ON_ASSEMBLY/contigs.fasta MAPPING/no_muricauda_R1.fastq.gz MAPPING/no_muricauda_R2.fastq.gz | samtools view -F 12 -hbS - | samtools sort -@ 4 - > MAPPING_ON_ASSEMBLY/contigs.bam ; samtools index MAPPING_ON_ASSEMBLY/contigs.bam && conda deactivate" > bwa_contigs.sh
qsub -cwd -V -q formation.q -N bwa_contigs -o LOG -e LOG -pe thread 9 -R y bwa_contigs.sh

Étape 4: Combien de reads sont mappés sur les contigs ?

cd MAPPING_ON_ASSEMBLY
conda activate samtools
samtools flagstat contigs.bam
conda deactivate

3404773 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
4031 + 0 supplementary
0 + 0 duplicates
3382088 + 0 mapped (99.33% : N/A)
3400742 + 0 paired in sequencing
1700371 + 0 read1
1700371 + 0 read2
3288712 + 0 properly paired (96.71% : N/A)
3368028 + 0 with itself and mate mapped
10029 + 0 singletons (0.29% : N/A)
74498 + 0 with mate mapped to a different chr
73251 + 0 with mate mapped to a different chr (mapQ>=5)

Étape 5: Utiliser samtools pour évaluer la profondeur et la couverture

conda activate samtools
samtools coverage contigs.bam > contigs.bam.cov
conda deactivate

Binning

MetaBAT

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 s’approcher 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’utiliseMetaBAT Kang et al. (2015) (Metagenome Binning with Abundance and Tetra-nucleotide frequencies) est basé sur la fréquence des tétranucléotides ainsi que les abondances des contigs en fonction des différents échantillons. Bien que MetaBAT soit optimisé pour une utilisation avec des échantillons nombreux et varié, il est possible de l’utiliser sur un unique échantillon ou même uniquement sur la distribution en 4-mer des contigs (c’est à dire, contigs FASTA sans alignement BAM).

cd ../
echo "conda activate metabat2-2.12.1 && metabat2 --numThreads 8 --verbose --unbinned --inFile SPADES/contigs.fasta --outFile METABAT/metabat --abdFile MAPPING_ON_ASSEMBLY/contigs.bam.cov --seed 42 && conda deactivate" > metabat.sh
qsub -cwd -V -q formation.q -N metabat -pe thread 8 -o LOG/ -e LOG/ metabat.sh

cd METABAT
conda activate seqkit-0.11.0 && seqkit stat *.fa  && conda deactivate
cd ../

Sur le jeu de données réel

file                         format  type  num_seqs    sum_len  min_len   avg_len  max_len
METABAT/metabat.10.fa        FASTA   DNA         79  3,972,900    2,833  50,289.9  244,528
METABAT/metabat.1.fa         FASTA   DNA        305  1,043,115    1,680     3,420    8,355
METABAT/metabat.2.fa         FASTA   DNA         50    255,987    2,744   5,119.7    7,662
METABAT/metabat.3.fa         FASTA   DNA        290  1,565,251    2,502   5,397.4   29,816
METABAT/metabat.4.fa         FASTA   DNA        112  1,437,788    2,627  12,837.4   25,098
METABAT/metabat.5.fa         FASTA   DNA         52  2,282,155    2,920  43,887.6  116,560
METABAT/metabat.6.fa         FASTA   DNA         74  2,982,386    3,032  40,302.5  153,121
METABAT/metabat.7.fa         FASTA   DNA         69    656,990    3,491   9,521.6   17,523
METABAT/metabat.8.fa         FASTA   DNA        281  3,994,721    2,590  14,216.1   96,608
METABAT/metabat.9.fa         FASTA   DNA         81  1,424,958    2,805  17,592.1   26,447
METABAT/metabat.unbinned.fa  FASTA   DNA      2,521  5,095,982    1,000   2,021.4   84,317

Évaluation des bins

Pour l’évaluation des bins, nous utiliserons les deux métriques completeness et contamination estimés par CheckM Parks et al. (2015). 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 utiliseront le workflow lineage-specific qui est recommandé pour évaluer l’exhaustivité et la contamination des bins de génomes.

On exécute le workflow lineage-specific en spécifiant le dossier des bins.

mkdir CHECKM/
echo "HOME=/projet/tmp ; conda activate checkm-genome-1.0.18 && checkm lineage_wf --extension fa --threads 8 --file CHECKM/report.tsv --tab_table METABAT/ CHECKM/ && conda deactivate" > checkm.sh
qsub -cwd -V -q formation.q -N checkm -pe thread 8 -o LOG/ -e LOG/ checkm.sh

On peut explorer plus en détail les résultats avec checkm qa

qsub -cwd -V -q formation.q -N checkm -o LOG -e LOG -b y "conda activate checkm-genome-1.0.18 && checkm qa CHECKM/lineage.ms CHECKM/ -o1 && conda deactivate"
# checkm qa CHECKM/lineage.ms CHECKM/ -o2
# checkm qa CHECKM/lineage.ms CHECKM/ -o6

Résultats attendus

Bin Id	Marker lineage	# genomes	# markers	# marker sets	0	1	2	3	4	Completeness	Contamination	Strain heterogeneity
metabat.1	f__Rhodobacteraceae (UID3340)	84	568	330	415	152	1	0	0	26.77	0.30	0.00
metabat.2	root (UID1)	5656	56	24	56	0	0	0	0	0.00	0.00	0.00
metabat.3	k__Bacteria (UID203)	5449	104	58	71	33	0	0	0	46.03	0.00	0.00
metabat.4	k__Bacteria (UID203)	5449	104	58	84	20	0	0	0	31.03	0.00	0.00
metabat.5	k__Bacteria (UID203)	5449	104	58	64	40	0	0	0	56.90	0.00	0.00
metabat.6	o__Pseudomonadales (UID4679)	87	674	364	16	651	6	1	0	97.39	1.85	11.11
metabat.7	o__Rhizobiales (UID3447)	356	416	249	330	86	0	0	0	11.85	0.00	0.00
metabat.8	c__Gammaproteobacteria (UID4444)	263	506	232	18	478	10	0	0	94.07	3.36	60.00
metabat.9	k__Bacteria (UID203)	5449	104	58	90	14	0	0	0	22.41	0.00	0.00
metabat.10	c__Gammaproteobacteria (UID4444)	263	507	232	20	476	11	0	0	99.14	0.97	54.55
metabat.unbinned	k__Bacteria (UID203)	5449	104	58	25	52	22	4	1	66.58	22.23	27.50

Et représenter ces résultats sous forme de plots avec checkm bin_qa_plot

qsub -cwd -V -q formation.q -N checkm -o LOG -e LOG -b y "conda activate checkm-genome-1.0.18 && checkm bin_qa_plot -x fa CHECKM/ METABAT/ CHECKM/ && conda deactivate"

Résultats attendus

Bar plot of bin quality

Si l’on connaît la taxonomie associée à notre bin, on peut le spécifié en utilisant le workflows checkm taxonomy_wf

mkdir CHECKM_Micromonospora/
qsub -cwd -V -q formation.q -N checkm_gamma -pe thread 8 -o LOG/ -e LOG/ -b y "conda activate checkm-genome-1.0.18 && checkm taxonomy_wf genus Micromonospora --extension fa --threads 8 --file CHECKM_Micromonospora/report.tsv --tab_table METABAT/ CHECKM_Micromonospora/ && checkm bin_qa_plot --extension fa CHECKM_Micromonospora/ METABAT/ CHECKM_Micromonospora/ && conda deactivate"

Prédiction de gènes

Prodigal

A partir de contigs au format FASTA, prodigal permet de prédire rapidement la structure des gènes codants pour des protéines (CDS). prodigal permet la prédiction de gènes sur des contigs de bactéries, d’archées, et de virus. prodigal est très rapide et prédit autant les gènes complets (de start à stop) que ceux morcelés. L’option -p meta permet de signaler à prodigal qu’il s’agit d’un assemblage métagénomique.

Lancer prodigal sur le jeu de données :

mkdir PRODIGAL
echo "conda activate prodigal-2.6.3 && prodigal -f gff -i SPADES/contigs.fasta -p meta -a PRODIGAL/genes.faa -d PRODIGAL/genes.fna -o PRODIGAL/genes.gff && conda deactivate" > prodigal.sh
qsub -cwd -V -q formation.q -N prodigal -o LOG -e LOG prodigal.sh

Prokka (Non conseillé)

A partir de contigs au format FASTA, prokka permet d’annoter rapidement les séquences bactériennes, archéennes et virales (ainsi que les métagénomes avec l’option --metagenome). L’utilisation de prokka est basée sur une suite d’outils : <strong class="tool>prodigal (détéction des gènes), RNAmmer (rRNA), aragorn (tRNA), signalP (Signal leader peptides), infernal (Non-coding RNA).

Au delà de l’annotation syntaxique (structure des gènes), réalisés par prodigal à partir de modèles HMM, Prokka assigne une fonction à chaque gène codant pour une protéine. Pour cela il aligne d’abord chaque CDS contre une banque de protéines bacteriennes expertisées issues de Swissprot, IS Finder et NCBI Antimicrobial Resistance Databank (blastp) et Pfam (hmmsearch). Les protéines présentant un “hit” avec une e-value inférieur à 1e-6 sont annotés avec la fonction, le nome de gène et éventuellement l’EC number de ce hit. Les protéines n’ayant pas de similarité suffisante pour inférer une fonction sont annotées comme hypothetical protein.

Ces deux approches permettent d’annoter de façon homogène et très rapide les familles de gènes conservées. Il est également possible de fournir à prokka un proteome (sous la forme d’un multi-fasta, cf. documentation) déjà annoté qu’il utilisera comme banque de référence pour l’annotation fonctionnelle. Cela peut être une stratégie intéressante pour l’annotation des bins. Prokka fournit donc une annotation syntaxique et fonctionnelle bien plus complète que Prodigal seul.

Cependant, l’auteur de Prokka ne recommande pas son utilisation sur des assemblages métagénomiques, car la version actuelle de Prokka ne traite que les CDS complets ( Cf https://github.com/tseemann/prokka/issues/399 , “Prokka is not really designed for metagenomes.”).

echo "conda activate prokka-1.14.5 && prokka --outdir PROKKA --addgenes --kingdom Bacteria --metagenome --cpus 8 SPADES/contigs.fasta  && conda deactivate" > prokka.sh
qsub -cwd -V -q formation.q -N prokka -o LOG -e LOG -pe thread 8 prokka.sh

Récupérer les comptages par gène

Récupérer uniquement les features gene dans le fichier d’annotation gff fourni par prokka
Utiliser le mapping sur les contigs pour en extraire avec htseq-countles lecture s’alignant sur les gènes.
Se référer à la documentation pour choisir la methode de calcul.

mkdir COUNT

echo "conda activate htseq-0.12.4 && htseq-count --format bam --type CDS --idattr ID MAPPING_ON_ASSEMBLY/contigs.bam PRODIGAL/genes.gff  > COUNT/genes_count.tsv && conda deactivate" > map_genes.sh

qsub -cwd -V -q formation.q -N gene_count -o LOG -e LOG -pe thread 8 map_genes.sh

The script outputs a table with counts for each feature, followed by the special counters, which count reads that were not counted for any feature for various reasons. The names of the special counters all start with a double underscore, to facilitate filtering. (Note: The double unscore was absent up to version 0.5.4). The special counters are:

__no_feature: reads (or read pairs) which could not be assigned to any feature (set S as described above was empty).

__ambiguous: reads (or read pairs) which could have been assigned to more than one feature and hence were not counted for any of these, unless the –nonunique all option was used (set S had more than one element).

__too_low_aQual: reads (or read pairs) which were skipped due to the -a option, see below

__not_aligned: reads (or read pairs) in the SAM file without alignment

__alignment_not_unique: reads (or read pairs) with more than one reported alignment.

Annotation fonctionnelle sur les gènes

L’objectif de cette étape est d’annoter fonctionnellement les gènes codant pour des protéines, avec des ontologie ou des hiérarchies fonctionnelles permettant l’agrégation des annotation sur différents niveaux. Cette étape est coûteuse en temps de calcul.

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

eggnog-mapper

mkdir EGGNOG
echo "conda activate eggnog-mapper-2.1.1 &&  emapper.py -i PRODIGAL/genes.faa --output eggnog -m diamond -d bact --cpu 8 --data_dir /db/outils/eggnog-mapper-2.1.1/ --output_dir EGGNOG && conda deactivate" > eggnog.sh
qsub -cwd -V -q formation.q -N eggnog -o LOG -e LOG -pe thread 8 eggnog.sh

Annotation par similarité sur NR (diamond)

Diamond

Diamond est un outil similaire à BLAST mais avec des heuristiques accélérant très fortement les comparaisons des gros jeux de données contre des banques. Nous pouvons annoter fonctionnellement notre jeu de données par similarité contre la banque NCBI NR.

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

References

Andrews, S. 2010. “FastQC a Quality Control Tool for High Throughput Sequence Data.” Http://Www.bioinformatics.babraham.ac.uk/Projects/Fastqc/. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.

Bankevich, Anton, Sergey Nurk, Dmitry Antipov, Alexey A Gurevich, Mikhail Dvorkin, Alexander S Kulikov, Valery M Lesin, et al. 2012. “SPAdes: A New Genome Assembly Algorithm and Its Applications to Single-Cell Sequencing.” Journal of Computational Biology 19 (5): 455–77.

Ewels, Philip, Måns Magnusson, Sverker Lundin, and Max Käller. 2016. “MultiQC: Summarize Analysis Results for Multiple Tools and Samples in a Single Report.” Bioinformatics 32 (19): 3047–48.

Kang, Dongwan D, Jeff Froula, Rob Egan, and Zhong Wang. 2015. “MetaBAT, an Efficient Tool for Accurately Reconstructing Single Genomes from Complex Microbial Communities.” PeerJ 3: e1165.

Kopylova, Evguenia, Laurent Noé, and Hélène Touzet. 2012. “SortMeRNA: Fast and Accurate Filtering of Ribosomal RNAs in Metatranscriptomic Data.” Bioinformatics 28 (24): 3211–17.

Li, Dinghua, Chi-Man Liu, Ruibang Luo, Kunihiko Sadakane, and Tak-Wah Lam. 2015. “MEGAHIT: An Ultra-Fast Single-Node Solution for Large and Complex Metagenomics Assembly via Succinct de Bruijn Graph.” Bioinformatics 31 (10): 1674–76.

Li, Heng. 2013. “Aligning Sequence Reads, Clone Sequences and Assembly Contigs with BWA-MEM.” arXiv Preprint arXiv:1303.3997.

Menzel, Peter, Kim Lee Ng, and Anders Krogh. 2016. “Fast and Sensitive Taxonomic Classification for Metagenomics with Kaiju.” Nature Communications 7: 11257.

Mikheenko, Alla, Andrey Prjibelski, Vladislav Saveliev, Dmitry Antipov, and Alexey Gurevich. 2018. “Versatile genome assembly evaluation with QUAST-LG.” Bioinformatics 34 (13): i142–50. https://doi.org/10.1093/bioinformatics/bty266.

Parks, Donovan H, Michael Imelfort, Connor T Skennerton, Philip Hugenholtz, and Gene W Tyson. 2015. “CheckM: Assessing the Quality of Microbial Genomes Recovered from Isolates, Single Cells, and Metagenomes.” Genome Research 25 (7): 1043–55.

Shen, Wei, Shuai Le, Yan Li, and Fuquan Hu. 2016. “SeqKit: A Cross-Platform and Ultrafast Toolkit for FASTA/q File Manipulation.” PloS One 11 (10): e0163962.

“SRA Tools.” 2020. http://ncbi.github.io/sra-tools/; NCBI.

Zhou, Yanqing, Yaru Chen, Shifu Chen, and Jia Gu. 2018. “Fastp: An Ultra-Fast All-in-One FASTQ Preprocessor.” Bioinformatics 34 (17): i884–90. https://doi.org/10.1093/bioinformatics/bty560.

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

Shotgun metagenomics

Migale facility

Valentin Loux - Cédric Midoux - Olivier Rué

2021-06-25

Rappels sur les ressources de la plateforme Migale

Lancer un job

Analyse d’un jeu de données métagénomiques

Récupération des fichiers

Préparation de l’espace de travail

Contrôle qualité des fichiers

Nombre de reads

Contenu des fichiers

Analyse préliminaire des reads bruts

Contrôle Qualité des reads

FastQC

MultiQC

Assignation taxonomique

Nettoyage

Filtrer et trimmer les reads

Séparer les reads d’ARN ribosomique

Enlever des contaminants par mapping

Normalisation avant assemblage (digital normalization)

Assembler les reads

Assemblage

Spades

Megahit

Évaluation de la qualité de l’assemblage

Mapping des reads sur les contigs

Binning

MetaBAT

Évaluation des bins

Prédiction de gènes

Prodigal

Prokka (Non conseillé)

Récupérer les comptages par gène

Annotation fonctionnelle sur les gènes

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

eggnog-mapper

Annotation par similarité sur NR (diamond)

Diamond

References