Recipe View

Mission Impossible: Multi Sample Variant Calling in 60 seconds

Your mission, should you choose to accept it, is to reproduce major findings of a Science paper in …

1 result • updated 2.4 years ago by Istvan Albert

Your mission, should you choose to accept it, is to reproduce major findings of a Science paper in just 60 seconds.

The goal of the recipe is to demonstrate that, when properly documented it is possible to reproduce major scientific findings with relative ease and perhaps under a minute.

This recipe reproduces the main results published in the work:

Genomic surveillance elucidates Ebola virus origin and transmission during the 2014 outbreak published in 2014 in the Science research journal.

On a properly set up system, and depending on internet bandwidth, the recipe completes under a minute [1].

The following steps are performed in the recipe:

Get the Ebola KJ660346 GenBank file from NCBI.
Reformat the GenBank into FASTA and GFF.
Download sequencing run information (SRR numbers) using the project id associated with the publication
Download FASTQ data based on SRR numbers
Align each FASTQ file to the reference with bwa
Create a multisample VCF file with bcftools for all alignment files.
Annotate the variants with snpEff to summarize the types of variations seen

The main results of the analysis can be found in the files:

all_effects_report.html
all_variants_annotated.vcf

All that ... in less than a minute. How is this even possible?

It is not magic, it is standing on the shoulders of giants that created tools like bwa, bcftools, SRA and NCBI!

The recommended way to enjoy this analysis is to have the following video playing while the analysis is running:

https://www.youtube.com/watch?v=yLo9k8Gifb0

We found that the theme song makes even the otherwise bland analysis feel super exciting.

1. Disclaimer: It is the second time when the process will take less than a minute, as the SRA files will be already cached on the system. The first time it will be more sensitive to the current network speed. 😉

Copy recipe

You need write access to the project to edit.

# Stop on any error. Show the commands.
set -uex

# Install the bio package if not present.
#pip install bio --upgrade

# Download the database for ebola.
#snpeff download -v ebola_zaire

# NCBI Bioproject ID.
ID=PRJNA257197

# The number of samples to process.
N=5

# NCBI accession number for the reference genome.
ACC=KJ660346

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

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

# The directory that stores the reference files.
mkdir -p refs

# Get the reference sequence in GenBank format.
bio fetch -format fasta $ACC  > $FASTA

# Get the reference sequence in GenBank format.
bio fetch -format gff $ACC > $GFF

# Index the FASTA file for bwa.
bwa index $FASTA 2>> log.txt

# Index the FASTA file for IGV.
samtools faidx $FASTA >> log.txt

# Obtain the SRR run number for the data deposited into the BioProject
esearch -db sra -query $ID | efetch -format runinfo > runinfo.csv

# Keep the top N SRR numbers for now.
cat runinfo.csv  | cut -f 1 -d , | grep SRR | head -$N  > samples.txt

# Make a directory for the sequencing data.
mkdir -p reads

# Limit each dataset to these many reads.
LIMIT=100000

# Download the sequencing data in parallel.
cat samples.txt | parallel "fastq-dump -X $LIMIT -O reads --split-files {} >> log.txt"

# Generate alignments
mkdir -p bam

# Run the bwa mem processes in parallel.
# See also the book titled: The Art of Bioinformatics Scripting

# Put an echo in front of the command to understand what it does.
cat samples.txt | parallel "bwa mem $FASTA reads/{}_1.fastq reads/{}_2.fastq 2>> log.txt | samtools sort > bam/{}.bam 2>> log.txt"

# Index the BAM files.
cat samples.txt | parallel "samtools index bam/{}.bam 2>> log.txt"

# Call variants on all BAM files in parallel.
bcftools mpileup -Ou -f $FASTA  bam/*.bam 2>> log.txt | bcftools call --ploidy 1 -vm -Ou |  bcftools norm -Ov -f $FASTA -d all - > all_variants.vcf

# Annotate the variants with SNPeff.
cat all_variants.vcf | snpEff ebola_zaire -s all_effects_report.html > all_variants_annotated.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).

Your mission, should you choose to accept it, is to reproduce major findings of a Science paper in just 60 seconds.

The goal of the recipe is to demonstrate that, when properly documented it is possible to reproduce major scientific findings with relative ease and perhaps under a minute.

This recipe reproduces the main results  published in the work:

* [Genomic surveillance elucidates Ebola virus origin and transmission during the 2014 outbreak][pub] published in 2014 in the Science research journal.

On a properly set up system, and depending on internet bandwidth, the recipe completes under a minute \[1\].

The following steps are performed in the recipe:

1. Get the Ebola KJ660346 GenBank file from NCBI.
1. Reformat the GenBank into FASTA and GFF.
1. Download sequencing run information (SRR numbers) using the project id associated with the publication
1. Download FASTQ data based on SRR numbers
1. Align each FASTQ file to the reference with `bwa`
1. Create a multisample VCF file with `bcftools` for all alignment files.
1. Annotate the variants with `snpEff` to summarize the types of variations seen

The main results of the analysis can be found in the files:

1. all_effects_report.html
 1. all_variants_annotated.vcf

All that ... in less than a minute. How is this even possible?

It is not magic, it is standing on the shoulders of giants that created tools like `bwa`, `bcftools`, SRA and NCBI!

The recommended way to enjoy this analysis is to have the following video playing while the analysis is running:

https://www.youtube.com/watch?v=yLo9k8Gifb0

[![mission impossible][1]][video]

We found that the theme song makes even the otherwise bland analysis feel super exciting.

*1\. Disclaimer: It is the second time when the process will take less than a minute, as the SRA files will be already cached on the system. The first time it will be more sensitive to the current network speed.* 😉

[pub]: http://science.sciencemag.org/content/early/2014/08/27/science.1259657.full

[video]: https://www.youtube.com/watch?v=yLo9k8Gifb0
[1]: /media/images/mission-impossible.png

Insert Image

From the web

From your computer

Cancel

Back