Recipe View

How to generate variant calls from sequencing data

This recipe demonstrates variant call generation from sequencing deposited at SRA

1 result • updated 5.1 years ago by Istvan Albert

This recipe demonstrates variant call generation from sequencing deposited at SRA

For more information see the

Bioinformatics Data Analysis online course.
Biostar Handbook for system setup and other information.

Copy recipe

You need write access to the project to edit.

# Stop on errors. Show the commands.
set -uex

# Reference genome accession number.
ACC=AF086833

# SRR run number.
SRR=SRR1553500

# The directory to store the reference in.
mkdir -p refs

# The name of the file that stores the reference genome.
REF=refs/$ACC.fa

# The resulting alignment file name.
BAM=$SRR.bam

# Get the reference sequence.
efetch -db=nuccore -format=fasta -id=$ACC | seqret -filter -sid $ACC > $REF

# Index reference for the aligner.
bwa index $REF >> log.txt 2>> log.txt

# Index the reference genome for IGV
samtools faidx $REF

# Get data from an Ebola sequencing run.
fastq-dump -X 100000 --split-files $SRR >> log.txt

# Shortcut to read names.
R1=${SRR}_1.fastq
R2=${SRR}_2.fastq

# Tag the alignments. GATK will only work when reads are tagged.
TAG="@RG\tID:$SRR\tSM:$SRR\tLB:$SRR"

# Align and generate a BAM file.
bwa mem -R $TAG $REF $R1 $R2 2>> log.txt | samtools sort > $BAM 2>> log.txt

# Index the BAM file.
samtools index $BAM

# Determine the genotype likelihoods for each base.
bcftools mpileup -Ovu -f $REF $BAM > genotypes.vcf

# Call variants with bcftools.
bcftools call -vm -Ov genotypes.vcf > variants1.vcf 2>> log.txt

# Call variants with freebayes.
freebayes --ploidy 1 -f $REF $BAM > variants2.vcf

# Call variants with GATK. See the book for GATK installation.
# These commands are turned off in this example.

# It needs yet another so called Dictionary index.
# picard CreateSequenceDictionary REFERENCE=$REF  OUTPUT=refs/$ACC.dict

# Generate the variants.
#gatk HaplotypeCaller -R $REF -I $BAM -O variants3.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

Edit the content of each interface element.

[acc]
label = "Genome accession number"
display = "TEXTBOX"
value = "AF086833"
regex = "^\\w{1,20}$"
help = "Must be an NCBI accession number"

[srr]
label = "The SRR run number"
display = "TEXTBOX"
value = "SRR1553500"
regex = "^\\w{1,20}$"
help = "Must be an SRA run id number"

[N]
label = "The number of reads"
display = "INTEGER"
value = 100000
range = [ 1, 1000000,]
help = "How many simulated read pairs to generate."

[settings]
id = 87
recipe_uid = "recipe-variant-example"
uid = "recipe-variant-example"
name = "How to generate variant calls from sequencing data"
template = "How_to_generate_variant_calls_from_sequencing_data_bio-data-analysis_87.sh"
image = "How_to_generate_variant_calls_from_sequencing_data_bio-data-analysis_87.png"
project_uid = "2a988f7f"
url = "http://localhost8000"
help = "This recipe demonstrates variant call generation from sequencing deposited at SRA\n\nFor more information see the\n\n* [Bioinformatics Data Analysis][link]  online course.\n* [Biostar Handbook][book] for system setup and other information.\n\n[link]: https://www.biostarhandbook.com/edu/course/4/\n[book]: https://www.biostarhandbook.com"

You need write access to the original recipe to edit.

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