Recipe View

A curious case of variant calling errors

Variant miscalling demonstration

1 result • updated 2.4 years ago by Istvan Albert • needs authorization

Variant miscalling demonstration

This recipe is built to demonstrate an interesting case where a two-base deletion in a certain region of the Ebola genome can lead to variant callers systematically miscalling that region.

This recipe requires the bio package, install it with pip install bio --upgrade

The recipe proceeds as follows:

downloads the ebola genome corresponding to accession number AF086833
extracts the first 3000 basepairs of it and call that the REFERENCE.fa
applies mutations and calls the resulting file GENOME.fa
generates simulated data from GENOME.fa and aligns it against REFERENCE.fa
also aligns GENOME.fa against REFERENCE.fa as a "whole genome alignment"
performs snp calling with bcftools and freebayes

The region around tother variations is prone to so called "misalignment" and that often produces errors that substantially affect snp calling accuracy. The mutations that are present in the data are the following:

AGGGTGGACAACAGAAGAACA
|||||||||||||||-|||--
AGGGTGGACAACAGA-GAA--

Here is an example, bcftools generates three different variants in the location, freebayes just one. (the quality of the calls could be filtered further!)

Doesn't always happen, but either snp caller of the two may indicate additional variants at position 1998 (you may need to run the recipe multiple times).

Copy recipe

You need write access to the project to edit.

# Show the commands as are executed.
set -uex

# Install bio with:
# pip install bio --upgrade

# Create a directory to store the reference and all indices.
mkdir -p refs

# The accession number.
ACC=AF086833

# The reference genome we align to.
REF=refs/REFERENCE.fa

# The "real" genome that is different from the reference.
GENOME=GENOME.fa

# The read alignment for the simulated reads.
RESULTS=results.bam

# The semi-global alignment of the genome to the reference.
REALITY=reality.bam

# Error rate for the simulation.
ERR=0.05

# How many reads to generate.
N=2500

# Obtain the data save the first 3kb of it as our reference FASTA.
bio fetch -format fasta $ACC | bio fasta --end 3kb > $REF

# We create the modified genome here. 
# We are replacing a sequence pattern with a different version.
# You are welcome to experiment with your changes.
MATCHED=AGGGTGGACAACAGAAGAACA
REPLACE=AGGGTGGACAACAGAGAA

# Modify the genome.
cat $REF | biosed -filter -target $MATCHED -replace $REPLACE > $GENOME

# Show the replacement
bio align $MATCHED $REPLACE -matrix NUC.4.4

# Create a genome wide alignment of the refence vs the genome.
# This is only possilble because we "know" what the actual genome is.
minimap2 -a $REF $GENOME | samtools sort > $REALITY

# Read names
R1=read1.fq
R2=read2.fq

# Generate simulated reads with no other mutations.
wgsim -N $N -e $ERR -r 0 $GENOME $R1 $R2

# Index the reference genome with bwa index.
bwa index $REF

# READ tagging means adding more information into each alignment.
TAG="@RG\tID:variants\tSM:${ACC}\tLB:simulation"

# Align the reads and generate athe BAM file.
bwa mem -R $TAG $REF $R1 $R2 | samtools sort > $RESULTS

# Create the indexes for the BAM files.
samtools index $RESULTS
samtools index $REALITY

# Variant call with bcftools.
bcftools mpileup -Ov -f $REF $RESULTS | bcftools call -mv -Ov -o bcftools.vcf

# Variant call with freebayes. Filter at 0.1 pval.
freebayes -P 0.1 -f $REF $RESULTS > freebayes.vcf

You need write access to the original recipe to edit.

Click the buttons on the right to create new fields.

Add text field Add float field Add data field Add checkbox Add dropdown Add upload field Add integer field Add radio button

Name

Recipe display name

Identifier

Unique identifier for the recipe.

Image :

Optional image for the recipe ( 500px Maximum ).

Rank:

Used to order recipes (optional).

Variant miscalling demonstration

This recipe is built to demonstrate an interesting case where a two-base deletion in a certain region of the Ebola genome can lead to variant callers systematically miscalling that region.

This recipe requires the [bio package][bio], install it with `pip install bio --upgrade`

The recipe proceeds as follows:

1. downloads the ebola genome corresponding to accession number `AF086833` 
 2. extracts the first 3000 basepairs of it and call that the `REFERENCE.fa`
 3. applies mutations and calls the resulting file `GENOME.fa`
 4. generates simulated data from `GENOME.fa` and aligns it against `REFERENCE.fa` 
 5. also aligns `GENOME.fa` against `REFERENCE.fa` as a "whole genome alignment" 
 7. performs snp calling with `bcftools` and `freebayes`

The region around tother variations is prone to so called "misalignment" and that often produces errors that substantially affect snp calling accuracy.  The mutations that are present in the data are the following:

AGGGTGGACAACAGAAGAACA
    |||||||||||||||-|||--
    AGGGTGGACAACAGA-GAA--

Here is an example, bcftools generates three different variants in the location, freebayes just one. (the quality of the calls could be filtered further!)
![enter image description here][1]

Doesn't always happen, but either snp caller of the two may indicate additional variants at position 1998 (you may need to run the recipe multiple times).

[bio]: https://ialbert.github.io/bio/

[1]: /media/images/grab2.png
  [2]: /media/images/Screenshot_2020-12-01_223817.png

Insert Image

From the web

From your computer

Cancel

Back