Variant calling with reads mapping

How does variant calling with reads mapping work

You need to map the reads to a reference.
You need to call your variants (detection of what is different to the reference).
You filter the variants according to some quality criteria.

We will use Listeria monocytogenes from the EURL proficiency test 2021 in this tutorial.

Our objectives with this tutorial are to:

Show you what data a variant caller that use mapped reads has to deal with before producing results you can trust.
How the choice of a reference influence what/which variants are detected
Provide a basic understanding of what and how a variant caller filters noise/errors in the mapped data in order to provide reliable variant calls that you can trust (or not).
Provide some ideas of how you can adapt variant calling to your purposes.

The data used for this tutorial is deposited in the SAGA directory : /cluster/projects/nn9305k/tutorial/20211108_trusting_snps

Information about the isolates used in this tutorial:

Isolates have been sequenced with NextSeq.
Assemblies have been done using bifrost

We evaluated the assembly quality with quast, using the L. monocytogenes reference genome: EGD-e (CA_000196035.1_ASM19603v1_genomic). The report indicate that the strain Lab9_st166 share a high genome fraction with the reference EGD-e, whereas the strain Lab9_st162 share around 60% of the genome with EGD-e. We will use those two isolates during the tutorial.

What is a variant: read mapping and what does the variant caller has to deal with?

A variant represents a difference to another isolate/reference, for instance a point mutation or SNP, an insertion/deletion (which are the subject for today) or eg. structural variants.

We have mapped those strains with BWA MEM (eg. see tutorial). BWA MEM , a mapping method that performs local alignment. We have used this program against:

the reference EGD-e
another reference that is more closely related to the strain (see note below: How to rapidly find reference for you isolates when you have no idea what to choose ?)
their own assembly. We used the genome sequence .fna provided by prokka annotation for mapping isolates against their own assembly to insure that locus tags are consistent with annotation and ease visualization with IGV Integrative Genomic Viewer (see tutorial).

We used the defaults parameters for the mapping. Note however this might be adjusted to allow more strict matches of the reads to the sequence. We have used both raw reads and trimmed reads for mapping. So we can have an idea of the amount of noise (incl. false variants) that a variant caller must deal with.

You will find the files for the tutorial in SAGA:

mapping for strain: Lab9_st162 /cluster/projects/nn9305k/tutorial/20211108_trusting_snps/Lab9_st162
mapping for strain: Lab9_st162 /cluster/projects/nn9305k/tutorial/20211108_trusting_snps/Lab9_st166
the annotation files for each genome used for reference mapping are provided in the sub-folder prokka_gff

PS: The BAM/SAM format of sequence/reads mapping is described here

Summary: how to load data into IGV

We will now use the program IGV to look at the mapping results.

Loading your assembly: In IGV meny select: Genomes > Load genome from file and choose the assembly/reference to load (.fna)
Add the mapped reads (raw, trimmed) and the annotation file as track In IGV meny select: File > Load from file and choose the mapped read file (.bam) and the .gff file corresponding to the assembly/reference

We can look for genes to manouvre in the genome: in the search box (left of go) type eg. murJ. The allows you to quickly find the gene: Lipid II flippase MurJ, (Locus murJ_1 and murJ_2)

In here you can find information on the color coding for alignments in IGV and how to change the different color coding views.

The color coding used in the viewer is described in IGV manual

In red: insert size that is larger than expected (possible evidence of a deletion).
In blue: inferred insert size that is smaller than expected (possible evidence of an insertion).
Mismatched bases in aliments are assigned a transparency value (proportional to the phred quality score).
Transparency (white transparent box) is used for mapped reads with low mapping quality (eg. mapped at several loci). For more information look here.

Practical: Mapping approach and data noise/errors

We will in this section have a look whether we can see the following things:

Sequencing errors / low quality bases
Adapters (not trimmed, incompletely trimmed)

Visualizing with IGV: Comparing one isolate where both trimmed reads and raw reads have been mapped to its own assembly

choose one strain. Upload the assembly, and the mapped reads (both raw reads and trimmed reads) and upload the annotation file.
find the gene murJ annotation murJ_2. You can zoom in an out

Here is an idea of what you should be able to look at:

Raw VS trimmed	Colors Interpretation

The reads with many vertical lines are reads that are poorly mapped to the genome or have bases with low quality scores (where many inconsistencies with the consensus sequence are found).
Look at the pileup, and zoom in and out.
What are the likely errors, what are the adapters …

NB: With Lab9_st162 mapped to EGD-e with trimmed reads, you can see reads with low mapping (eg. mapped at several places) by looking around the gene ssb.

When mapping directly from raw-reads: the variant caller will have to deal with noise, to distinguish between real SNPs in pileup vs sequencing errors and remaining adapters. Soft-clipping at the extremity of the reads, is excluding those extremities as to be considered part of the alignment.

Practical: Influence of reference choice for variant calling using a mapping approach

Variant detection on a single isolate using a reference:

closely related reference
distantly related reference

Visualisation with IGV: Choose one isolate mapped with trimmed reads to two different references sequences Here we use trimmed reads (so we remove some low quality and noise data). This will allow us to look at the effect of reference choice for variant detection PS: if you do not find murJlook at eg. MurA, MurA1 genes and scroll around.

In IGV the symbol I indicates insertions and black - deletions.

Lab9_st162 to EGD-e	Lab9_st162 to 73_clip (close)

The amount of variant discovered depend on how closely related your reference is to the strain(s) used for discovering variants. Questions:

Lab9_st166 is quiet closely related to EGD-e. Did choosing an other closely related strain as reference reduced the number of potential snps?
Lab9_st162 is quite distantly related to EGD-e. How do you think the reference choice will influence your analyses for SNP detection?

The variant caller must discriminate between real variant and “noise” (sequencing error and absence of trimming).

The choice of the reference used in variant calling will influence which variants are detected, as a variant is a difference in base relative to another sequence used as reference. For detecting variants between multiple samples, a common reference has to be used.

Having a common reference that is as closely related as possible to the common set of isolates is important when you want to do phylogenies based on SNPs. Every site that produce indel will be removed from the phylogenetic reconstruction (unless a special model that account for indels is used).

If a reference it too distantly related, it is possible that many of your reads wont map to the reference (as they are too different, they might not be considered for local alignment). This might reduce your ability to detect variants in the loci that not very conserved in the genome.

Mapping quality scores / mapping quality: MAPQ

There is little or no documentation in mapping software on how the different software are computing mapping quality scores.

BWA-MEM has a score for alignment that depends on the parameters used to perform the local alignment and penalties associated to mismatches detection, penalties associated to gap opening and gap extension. Mapping quality scores (MAPQ) in BWA MEM are phred scaled and in the range [0-60], where 60 corresponds to uniquely mapped reads (high quality).

Note that the MAPQ scores are not directly comparable in different mapping software see eg. Olson et al. 2015. Note also that BWA and BWA MEM are two different algorithms with different mapping scores.

Here are some blogs that relate to this problem:

You can also have a look at the article Li et al. 2008 and their supplementary information

Note: How to rapidly find reference for you isolates when you have no idea what to choose ?

You can eg. use the software poppunk on a diverse set of reference isolates. Then if all your isolates belong to the same cluster, you can choose a reference that belongs this same cluster. If not all your isolates belong to the same cluster, you can look at the phylogenetic tree and choose eg. a reference that belong to the same monophyletic group (most recent common ancestor/node shared with all your isolates).

Another solution, is eg. to run a type of typing eg. 7 gene MLST typing and choose a reference that belong to the same ST. This can be done eg. with MLST software

For this tutorial, we used poppunk (version 2.4.0) to search for an appropriate references, as it is fast to do for a large dataset, much faster than a core phylogeny. We used the sequence-dataset composed of the isolates that were used to define the Listeria monocytogenes cgMLST scheme described in (Moura et al. 2017). The strain diversity in this dataset is expected to be representative of L.monocytogenes diversity, therefore allowing us to find alternative reference for SNP calling with our strains.

A tabular delimited file containing the ID and the paths for poppunk is created, without headers. We removed potential outlier genome by pruning (it just remove the sequences that are outlier from the analyses).

NB: We have checked that the model fit well (it is not always the case, but for Listeria monocytogenes it works well).

We created a phylogenetic tree to visually find which references sequences are more similar to our own genomes.

# These are the command used for creating the tree to find a reference using poppunk
cd <popunk_analysis_directory>
# create a tab delimted file with ID and path, no headers
# poppunk 2.4.0
conda activate poppunk
# creating the database and distance - used the mourra isolates - removing eventual poor isolates
poppunk --create-db --output listeria_db --r-files r_files.txt --qc-filter prune
# fitting the model
poppunk --ref-db listeria_db --fit-model dbscan
# we get the visualisation
poppunk_visualise --ref-db  listeria_db --microreact --tree nj --output visualisation_out

Then we can chose a reference that is more closely related to our isolates than our initial choice

Lab9_st162	Lab9_st166
Closer reference for 73_CLIP_2008_01408	Closer reference for 9_FSL_R2_561

How to find variants: example using snippy 

Snippy uses a mapping to a reference to find variants (SNPs/indels). The use of a common single reference for variant detection on multiple isolates provides a coordinate system (the position of each SNPs for each isolate is given according to the position in the reference sequence) that allows to combine the variants found in each isolate together in a single variant alignment.

Typed variants can be represented in a variant call format .vcf file. This file report the differences of nucleotides/indels of a/several isolates relative to the reference sequence. It also provides the statistics calculated by the variant caller eg. the likelihood of a variant base to be different to the reference. A set of measures are provided (depending on variant caller), eg. base confidence at mapping, filters used, probabilities calculated in the pileup. An overview about the format can be found here or in wikipedia for a fast overview.

Snippy removes duplicated reads prior to variant calling. (Duplicated reads would lead to an inflated weight when calling variant based on how frequently they are found in a pileup). This is something that you should keep in mind if your read set contains a lot of duplicated reads.

Snippy can be run using an assembly (.fasta) or annotated genome (.gbk: genebank format) as reference. Here we will call variant using a reference provided in the genebank format, as this will ease visualization of the variants files and allow to navigate from loci to loci while using IGV.

Practical Create a directory for one isolate, and create a symlink to the files .

cd <directory>
# example simlink
ln -s /cluster/projects/nn9305k/tutorial/20211108_trusting_snps/snippy/input/Lab9_st162/* .
ln -s /cluster/projects/nn9305k/tutorial/20211108_trusting_snps/snippy/GCF_000196035.1.gbk .

# request ressources
screen
srun --account=nn9305k --mem=32G --cpus-per-task 4 --time=1:00:00 --pty bash -i
# activate snippy
conda activate snippy
# Run snippy (the --report option is heavy - please someone do without?)
snippy --outdir out_Lab9_st162 --ref GCF_000196035.1.gbk --R1 Lab9_st162_R1.fastq.gz --R2 Lab9_st162_R2.fastq.gz --cpus 4 --report

Look at the filtered variant file. The report is generated with samtools view.

less snps.filt.svc

Headers of interest
1. vcf File format version
2. command line to generate the variant file
3. descriptions of the abbreviations in variant calls (read depth, likelihood of , type of variant: eg indel, SNP)
Variant call
1. contig (chrom)
2. position(relative to reference)
3. variant allele
4. and other fields described in the header.

Bonus: detecting variants with snippy multi

Snippy run can be run on multiple samples using the same reference for variant detection.

You can eg. use Haukon’s script here. Or the example bellow that I use to create snippy multi-input files. In this example, reads pairs for each isolate are organized in their own folder. The folder name is named after the isolate ID. The input folder snippy_input_dir contains all the directories of all isolates we want to run with snippy.

# example of creating the snippy multi-input file
snippy_input_dir="/cluster/projects/nn9305k/active/...input/"
cd $snippy_input_dir

# creating an array of ID based on a list of folders
id=($(ls -d Lab9_st16*))
# checking that the this is correct
echo ${id[@]}

function create_snippy () {
	for ((i=0 ; i < ${#id[@]} ; ++i ))
	do
	R1=${snippy_input_dir}/${id[i]}/*_R1.fastq.gz
	R2=${snippy_input_dir}/${id[i]}/*_R2.fastq.gz
	printf "%s\t%s\t%s\n" ${id[i]} ${R1} ${R2}
	done
}
create_snippy > snippy_input.tab

Running snippy multi:

srun --account=nn9305k --mem-per-cpu=32G --time=1:00:00 --pty bash -i
snippy-multi snippy_input.tab --ref <path>/GCF_000196035.1.gbk > created_snippy_run_script.sh
# veryfing that the script has been created
less created_snippy_run_script.sh
# launching the script
bash ./created_snippy_run_script.sh
# wait ...
# cleaning the alignment for further processing (eg. gubbins)
snippy-clean_full_aln core.full.aln > clean.full.aln

Snippy produces a multiple core alignment. The common reference serve as common coordinate system and allow to reconstruct an alignment from the common part of the genomes (core). All the variant detected at all positions of the reference are merged into one file for all isolates under study.

How to trust your variants?: example using snippy

Principles of variant calling in snippy

Snippy uses several filters are used to remove errors/noise variants. Freebayes is the variant caller used by snippy.

reads are uniquely mapped with BWA-MEM (corresponding to a mapping score of 60).
Duplicated reads tagged and then removed from variant calling.
the number of bases that are allowed to be soft-clipped at the extremities of the reads before discarding alignment can be adjusted.
The variant caller Freebayes is used to determine the most likely haplotype in short-read alignments; by extension which base(s) is/are more likely at each position in the reference. The proportion/frequency of variants that differ to the reference can be adjusted to strengten discriminating between noise and potential true variant.
The set of variants is then filtered, to ensure removal of variants detected due to noise/errors. The filters rely on :
- QUAL = the probability that there is a polymorphism at the loci $\1-P(homozygous|Data)$. Snippy use a default of 100(%)
- DP = the coverage/depth (default to 10 in snippy)
- A default base-quality of minimum phred score=13 is used for a variant calling (error rate of ~5%) . This differs from the default settings of Freebayes (phred score=20, which correspond to an error rate of 1%. See eg. here.

NB: depending on your objective, you might consider to adjust the minimum coverage/depth, the minimum fraction/proportion and variant probability. Eg. if you (1) have a very deep sequencing (2) have contaminants or eg. (3) want to call “real” rare intra-isolate variants and not solely the majority variants (4) OR you want to solely obtain the majority variants.

Practical

We can now look at the variants typed by snippy.

We ran snippy using the option --report to create a variant report. It used samtools tview to create this report.

Have a fast look at the report

less snps.report.txt

Lets look at the variant calling files in IGV and compare it with previous mappings. Guide: IGV vcf viewing

load the reference genome ref.fa (in the the sub-folder reference from snippy output)
load the files:
- mapped reads snps.bam
- the reference annotation .gff (in the the sub-folder reference from snippy output)
- the unfiltered variants snps.raw.vcf
- the filtered variants snps.filt.vcf
- the final variant file snps.vcf

Look at the genes/loci: eg for the isolates, and also navigate through the genome view. Try to answer questions below.

Lab_st162: eg. look at: MurA (we looked at previously), MurZ, lmo2553 and eg. regulatory_3308 , lmo 2821 …
Lab9_st166: eg. look at: no variants MurA, around lmo 2566, lmo 2652, lmor14 (eg.poor mapping quality), lmo2396 (many messy mapped reads) …
scroll through the genome to get a general impression

Questions:

Do the variant appear randomly distributed ?
Where are variants filtered due to low coverage?
How often to it appears that variants are filtered due a low fraction of reads that differ to the reference?
How often do the 3 types of .vcf files appear not to be concordant ?
What should/could we have done in a better way/otherwise ?

eg. MurZ coverage	eg. lmo2553 Filtering effect