Pratique #3 (2h):Assemblage et mapping, comparaison basée sur les SNPs - DrAnneJ/DrAnneJ.github.io GitHub Wiki

1. Assembler avec megahit

megahit est un assembleur basé sur les graphs de Bruijn. Il est très rapide c'est pour cela que nous allons l'utiliser dans ce TP. Il faut quelques minutes pour réaliser un assemblage alors qu'il faut environ 1h avec spades sur un ordinateur de bureau.

Utiliser les reads nettoyés (_trim_R1.fq.gz et _trim_R2.fq.gz)

avec default settings (k list: 21,41,61,81,99)

Allez dans le dossier des fastq trimés

cd '/home/anne/Desktop/fastq_sub'

Copiez seulement la ligne correspondant à votre souche /!\ puis "enter"

megahit -1 BJ-CH1_trim_R1.fq.gz -2 BJ-CH1_trim_R2.fq.gz --out-dir BJ-CH1_megahit --min-contig-len 500

megahit -1 KpS13_trim_R1.fq.gz -2 KpS13_trim_R2.fq.gz --out-dir KpS13_megahit --min-contig-len 500

megahit -1 SA17_trim_R1.fq.gz -2 SA17_trim_R2.fq.gz --out-dir SA17_megahit --min-contig-len 500

Rq: assemblage avec spades 3.9.0: ne pas le faire car prends du temps...les résultats se trouvent dans le dossier Spades_assemblies

Le résultat final de l'assemblage est nommé "final.contigs.fa" et doit faire la taille attendue pour un génome de Kp.

2. Evaluer la qualité des assemblages

a. Quast sur l'assemblage megahit

QUAST: QUality ASsessment Tool for Genome Assemblies Version: 5.0.2

Usage:

python /home/anne/miniconda2/bin/quast [options] <files_with_contigs>

ALLEZ DANS LE FOLDER AVEC LES ASSEMBLIES avec cd

cd '/home/anne/Desktop/Megahit_assemblies'

Copiez seulement la ligne correspondant à votre souche /!\

python /home/anne/miniconda2/bin/quast -o quast_results_BJ-CH1_megahit -r MAGE_plasm_MGH78578.fa BJ-CH1_megahit_contigs.fa

python /home/anne/miniconda2/bin/quast -o quast_results_KpS13_megahit -r MAGE_plasm_MGH78578.fa KpS13_megahit_contigs.fa

python /home/anne/miniconda2/bin/quast -o uast_results_SA17_megahit -r MAGE_plasm_MGH78578.fa SA17_megahit_contigs.fa

Allez voir le fichier icarus.html

copiez collez le fichier report.tsv qui contient les infos importantes sur votre assemblage megahit

b. Quast sur d'autres assemblages

quast permet de comparer des assemblages différents pour choisir le meilleur

Pour nos 3 souches, il existe des assemblages réalisés avec l'outil CLC qui ont été publiés, nous les avons téléchargés afin de les comparer à notre assemblage megahit

Total DNA was extracted by the phenol-chloroform method. Libraries were constructed using Nextera technology and sequenced on an Illumina’s HiSeq-2000 using a 2 x 100 nucleotides (nt) paired-end strategy. All reads were processed to remove low quality or artefactual nucleotides, using sequentially sickle (github.com/najoshi/sickle), AlienTrimmer (1), and fqDuplicate (ftp.pasteur.fr/pub/gensoft/projects/fqtools). Read pairs were assembled using clc_novo_assemble (www.clcbio.com/products/clc-genomics-workbench) with a minimum contig size of 500 nt.

Nous avons aussi réalisé un assemblage avec l'outil spades

Donc nous avons 3 assemblages réalisés avec 3 outils différents Les résultats de quast sont dans le tableau /home/anne/Desktop/RESULTS/All_quast_results.ods

c. Comparer les résultats

Regarder le nombre de contigs et la longueur totale Quels sont les meilleurs assemblages ?

Rq: on a fait les assemblages Megahit et Spades sur une partie des reads ('subsampling' pour gagner du temps et on a trimmé différemment on ne peut donc pas obtenir les exactement les mêmes assemblages)

3. Mapping avec snippy

a. L'identification des SNPs avec snippy ("SNPs Calling")

Pour identifier les SNPs avec snippy il faut les fichiers suivants "Input Files":

a reference genome in FASTA or GENBANK format (can be in multiple contigs)
sequence read file(s) in FASTQ or FASTA format (can be .gz compressed) format
a folder to put the results in

Snippy va générer les fichiers suivants "Output Files":

Extension	Description
.tab	A simple tab-separated summary of all the variants
.vcf	The final annotated variants in VCF format
.bam	The alignments in BAM format. Includes unmapped, multimapping reads. Excludes duplicates.
.bam.bai	Index for the .bam file

Les différents types de variants sont:

Type	Name	Example
snp	Single Nucleotide Polymorphism	A => T
mnp	Multiple Nuclotide Polymorphism	GC => AT
ins	Insertion	ATT => AGTT
del	Deletion	ACGG => ACG
complex	Combination of snp/mnp	ATTC => GTTA

b. Mise en pratique de snippy

Si vous avez réussi à obtenir vos propres fastq trimés allez dans votre dossier:

cd '/home/anne/Desktop/fastq_sub'

Sinon allez dans ce dossier:

cd '/home/anne/Desktop/RESULTS/FASTQ_trim'

i. On va mapper avec snippy nos reads trimmés contre un génome de ref annoté (fichier .gbk)

Nous avons 1 génome de ref:

MAGE_plasm_MGH78578.gbk

reads trimmés

BJ-CH1_trim_R1.fq.gz

BJ-CH1_trim_R2.fq.gz

KpS13_trim_R1.fq.gz

KpS13_trim_R2.fq.gz

SA17_trim_R1.fq.gz

SA17_trim_R2.fq.gz

on peut aussi chercher les snps sur des contigs, dans ce cas on ne peut pas évaluer leur qualité -> à votre avis pourquoi ?

ii. contre la référence `MAGE_plasm_MGH78578.gbk`

Copiez seulement la ligne correspondant à votre souche /!\

snippy --outdir snippy_SA17vsMGH --ref MAGE_plasm_MGH78578.gbk --pe1 SA17_trim_R1.fq.gz --pe2 SA17_trim_R2.fq.gz

snippy --outdir snippy_KpS13vsMGH --ref MAGE_plasm_MGH78578.gbk --pe1 KpS13_trim_R1.fq.gz --pe2 KpS13_trim_R2.fq.gz

snippy --outdir snippy_BJ-CH1vsMGH --ref MAGE_plasm_MGH78578.gbk --pe1 BJ-CH1_trim_R1.fq.gz --pe2 BJ-CH1_trim_R2.fq.gz

The snps.txt file in the SNIPPY output contains summary information -> open it !

Open the snps.tab file [with gedit then ctrlA ctrlC open LibrOffice Calc and ctrlV] -> understand the output !

(on peut modifier l'extension .tab en .tsv pour que le fichier suivre directement dans LibrOffice Calc)

4. Utiliser snippy pour faire une phylogénie basée sur les SNPs communs

In order to compare the identified SNPs with each other we need to know if a certain position exists in all isolates = "core sites". A "core site" is a genomic position that is present in all the samples. A core site can have the same nucleotide in every sample ("monomorphic") or some samples can be different ("polymorphic" or "variant").

If you call SNPs for multiple isolates from the same reference, you can produce an alignment of "core SNPs" which can be used to build a high-resolution phylogeny (ignoring possible recombination).

SNIPPY will concatenate the core SNPs, i.e. ignoring sites that are identical in all isolates and in the reference. For the phylogeny we ignore the complications of "ins", "del" variant types, and just use variant sites, these are the "core SNP genome".

Visualiser l'arbre basé sur les core SNPs

Open Dendroscope (icon on the desktop), go to File>Open and choose core_vsMGH.aln_phyml_tree.txt

Open Dendroscope choose core_vsNTUHK.aln_phyml_tree.txt

Ce fichier texte est en fait au format NEWICK. https://en.wikipedia.org/wiki/Newick_format

Il y a de nombreuses options pour modifier l'aspect de l'arbre.

Rq1: ces souches sont très éloignées les unes des autres on peut juste dire que NTUHK, KpS13 et SA17 sont plus proches entre elles qu'elles ne le sont de BJ-CH1 et MGH. En pratique si vous recherchez une épidémie les souches seront très proches (quelques SNPs et non des milliers comme ici). La phylogénie moléculaire est une discipline à part entière...

Nous avons récemment investigué une vraie épidémie à Salmonella Para-Typhi. Ouvrez le fichier tree_paratyphi_A.newick

Choisissez "Draw as rectangular phylogram" au lieu de "cladogram"

The main difference between cladogram and phylogenetic tree is that cladogram is an evolutionary tree with branches with equal distance, showing the relationship between a group of clades whereas phylogram is an evolutionary tree showing an estimate of phylogeny where the distance of the each branch is proportional to the amount of inferred evolutionary change.

Un des meilleurs outils de visualisation d'arbre est iTOL

https://itol.embl.de/upload.cgi

Copiez collez le contenu du fichier Newick et jouez avec les options de iTOL

5. Utiliser TABLET pour visualiser un mapping

Snippy a généré 2 fichiers de mapping que nous allons explorer (leur nom doit être identique et ils doivent être dans le même folder): snps.bam et snps.bam.bai

Ouvrir TABLET en tapant "tablet" puis "enter" dans le terminal

1- fasta assembly file

the .FASTA file for your reference genome sequence, it can be loaded by clicking on "Open Assembly" from the Ribbon Bar

The Primary assembly file refers to the main file containing your alignment data (.bam)

/home/anne/Desktop/Tablet_files/BJ1vsMGH_snps.bam

The Reference/consensus file refers to reference data (.fasta)

/home/anne/Desktop/Tablet_files/MGHref.fasta

2- GFF3 file

if you have the gene annotation file, it can be loaded

Select Import Features from the Ribbon Bar to import a GFF3 formatted file containing additional feature information that you want to associated with an assembly.

/home/anne/Desktop/Tablet_files/MAGE_plasm_MGH78578.gff

Once a feature file has been loaded, the Features Browser tab will become active whenever a contig is selected that contains features.

3- Faire des réglages

Faire un click droit au niveau de la fenêtre CIGAR -> select visible tracks -> décocher CIGAR et cocher CDS Read packing -> padded (paired end)

Dézoomez pour voir le mapping des reads sur la référence

Colour schemes -> choisir "Direction"

This view allows you to understand how paired-reads work !!

4- Chercher des régions non couvertes

Chercher des régions non couvertes (= régions variables entre les génomes)
Quelles sont les fonctions de ces régions variables ?

To obtain the sequence you can right click and select ‘Copy reference subsequence to clipboard’ then blastn on:

https://blast.ncbi.nlm.nih.gov/Blast.cgi

ex: KPN_01428, KPN_03404 (on peut aller directement au gène d'intérêt en allant sur le 2ème onglet) phage, pili

Aller sur MAGE

http://www.genoscope.cns.fr/agc/microscope/mage/viewer.php?O_id=897

et explorer les gènes KPN_01428, KPN_03404 en tapant le numero du gène dans VIEW puis enter

Exercice: trouvez un autre gène présent dans la référence qui est absent dans votre mapping et donner sa fonction probable

5- Identify SNPs manually

Passer à Colour shemes "nucleotides" Jouer sur le Zoom et sur Variants pour mieux les visualisezr

Colour schemes -> choisir "Variants"

Cherchez des vrais et des faux variants

Can you find any interesting SNPs? Which genes (if any) are they in? How reliable do they look? (Hint – look at the number of reads mapping, their orientation - which strand they are on and how good the base- calls are).

Exercice: trouvez un SNP de qualité et reporter sa position et le gène dans lequel il se trouve

Identify errors: Note that sequencing errors in Illumina data are quite common. We rely on depth of sequencing to average out these errors. This is why people often mention that a minimum median coverage of 20-30x across the genome is required for accurate SNP-calling with Illumina data.

Exercice: trouvez un faux SNP (erreur de séquençage) et reporter sa position et le gène dans lequel il se trouve

Recap: SNP/Indel identification

Genuine SNPs/Indels should be present on both read 1 and read 2 (see pairs)
Genuine SNPs/Indels should be present on both strands (arrows in two different directions)
Genuine SNPs/Indels should be supported by a good (i.e. 20-30x)

depth of coverage (at least 90% concordant)

Very important mutations (i.e. ones relied upon in a paper) should be confirmed via PCR/Sanger sequencing.