Recipe View

This recipe follows the variant calling process in the Biostar Handbook.

2 results • updated 2.5 years ago by Istvan Albert

This recipe follows the variant calling process in the Biostar Handbook.

Ensure that you have the bio package installed.

pip install bio --upgrade

The recipe obtains the reference genome AF086833 (the Ebola Mayinga strain of 1976) and will align it the sequencing data obtained from the 2014 outbreak deposited as SRA id SRR1972739

The recipe produces:

a genome for the reference genome
a GFF file for the reference genome
a BAM alignment file
a VCF file for the sequencing run.

For more information see:

Biostar Handbook handbook
Applied Bioinformatics course

Copy recipe

You need write access to the project to edit.

# Stop on errors. Show all commands.
set -eux 

# Install bio if you don't have it.
#pip install bio --upgrade

# Substitute reference genome accession number.
ACC=AF086833

# Substitute the SRR run number.
SRR=SRR1972739

# The alignment file name.
BAM=${SRR}.bam

# All genotype likelihoods.
ALL=${SRR}.all.vcf

# The variants selected from all genotypes.
VCF=${SRR}.variants.vcf

# The reference genome in FASTA format.
REF=refs/${ACC}.fa

# The reference genome in GFF format.
GFF=refs/${ACC}.gff

#
# Stage 1: Obtain the data
#

# Directory for various reference information.
mkdir -p refs

# Obtain the reference in FASTA format.
bio fetch -format fasta $ACC > $REF

# Obtain the reference in GFF format.
bio fetch -format gff $ACC > $GFF

# Index reference for the bwa aligner.
bwa index $REF

# Index the reference genome so that it can be loaded into IGV.
samtools faidx $REF

# The directory that stores the sequencing reads
mkdir -p reads

# Get sequencing data from a SRR number.
# Use fastq-dump-orig if you get an error
fastq-dump -X 10000 --split-files -O reads $SRR

# The paired end read names.
R1=reads/${SRR}_1.fastq
R2=reads/${SRR}_2.fastq

#
# Stage 2: aligning and calling variants.
#

# Align and generate a sorted BAM file.
bwa mem $REF $R1 $R2 | samtools sort > $BAM

# Index the BAM file.
samtools index $BAM

# Compute all genotypes with bcftools.
bcftools mpileup -O v -f $REF $BAM > $ALL

# From all genotypes call only the variants.
cat $ALL | bcftools call --ploidy 1 -mv -O v -o $VCF 

# We could run both steps in one stage by pipe the two into one call like so.
#bcftools mpileup -O v -f $REF $BAM | bcftools call --ploidy 1 -mv -O v -o $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).

Insert Image

From the web

From your computer

Cancel

Back

Variant calling with bcftools

Insert Image