Recipe View

Bacterial variant caller

1 result • updated 16 months ago by Istvan Albert

Bacterial variant caller

This recipe contains the code presented by Torsten Seemann in the blog post titled A Unix one-liner to call bacterial variants

The code from the above website has been generalized a bit with additional utility.

This recipe provides code that:

Downloads a reference genome (default value: NZ_CP008918 Pasteurella multocida strain ATCC 43137, complete genome)
Downloads a sequencing run from SRA (default value: SRR4124989)
Generates statistics on the sequencing data
Runs minimap2 to align the sequences and create a SAM output
Runs bcftools mpileup to generate the genotype likelihoods of each base
Runs bcftools call to filter for multiallelic variants only
Runs bcftools norm to normalize each variant to a standard form
Runs bcftools filter to remove variants with low quality (QUAL) or low coverage (DP)

Edit description

# Stop on errors.
set -uexo pipefail

# The reference genome.
ACC={{acc.value}}

# Select sequencing run.
SRR={{srr.value}}

# The directory that stores the reference sequence.
mkdir -p ref

# The file that stores the reference genome.
REF=ref/${ACC}.fa

# Download the reference accession number.
efetch -db nuccore -id $ACC -format fasta > $REF

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

# Download the sequencing data.
fastq-dump --origfmt  --split-files -O reads $SRR

# Optional step. Generate a report on the data size
echo "*** Sequence statistics"
seqkit stat reads/${SRR}_?.fastq
echo "*** "
#
# The one-liner code comes from:
#
# https://thegenomefactory.blogspot.com/2018/10/a-unix-one-liner-to-call-bacterial.html
#

# The number of CPUs to use.
CPUS=4

# The names for the read pairs.
R1=reads/${SRR}_1.fastq
R2=reads/${SRR}_2.fastq

# Unix trick: A single command line that is too long may be broken
# into several lines as long as each ends with a \ symbol.
#
# Call the variants.
minimap2 -a -x sr -t $CPUS $REF $R1 $R2 \
 | samtools sort -l 0 --threads $CPUS \
 | bcftools mpileup -Ou -B --min-MQ 60 -f $REF - \
 | bcftools call -Ou -v -m - \
 | bcftools norm -Ou -f $REF -d all - \
 | bcftools filter -Ov -e 'QUAL<40 || DP<10 || GT!="1/1"' \
 > variants.vcf

# Generate a statistics on variants.
bcftools stats variants.vcf > variants-stats.txt

# Plot the variant statistics.
plot-vcfstats -P -p plots variants-stats.txt

Back

[acc]
label = "Genome Accession Number"
display = "TEXTBOX"
value = "NZ_CP008918"
regex = "^\\w{1,20}$"
help = "An NCBI accession number"

[srr]
label = "SRA Run Number"
display = "TEXTBOX"
value = "SRR4124989"
help = "An SRR Run number"

[settings]
name = "Bacterial Variant Calling"
template = "Bacterial_Variant_Calling_cookbook_6.sh"
image = "Bacterial_Variant_Calling_cookbook_6.png"
id = 6
recipe_uid = "bacterialsnp"
uid = "bacterialsnp"
help = "Bacterial variant caller\n\nThis recipe contains the code presented by Torsten Seemann in the blog post titled [A Unix one-liner to call bacterial variants][code]\n\nThe code from the above website has been generalized a bit with additional utility.\n\nThis recipe provides code that:\n\n1. Downloads a reference genome (default value: [NZ_CP008918][NZ_CP008918] *Pasteurella multocida strain ATCC 43137, complete genome*)\n1. Downloads a sequencing run from SRA (default value: [SRR4124989][SRR4124989])\n1. Generates statistics on the sequencing data\n1. Runs `minimap2` to align the sequences and create a SAM output\n1. Runs `bcftools mpileup` to generate the genotype likelihoods of each base\n1. Runs `bcftools call` to filter for multiallelic variants only\n1. Runs `bcftools norm` to normalize each variant to a standard form\n1. Runs `bcftools filter` to remove variants with low quality (QUAL) or low coverage (DP)\n\n[NZ_CP008918]: https://www.ncbi.nlm.nih.gov/nuccore/NZ_CP008918.1\n[SRR4124989]: https://www.ncbi.nlm.nih.gov/sra/?term=SRR4124989\n\n[code]: https://thegenomefactory.blogspot.com/2018/10/a-unix-one-liner-to-call-bacterial.html"
url = "http://localhost8000"

Back

Interface preview shows the resulting view of the combined interface elements.

Interface element specification. Shows the code and the resulting interface element.

Double click on an element to insert it into the interface.

Integer values

                        [size]
                        label = "Window size"
                        display = "INTEGER"
                        value = 100
                        range = [1, 100]
                        help = "Selects the smoothing window."

Float values

                        [cutoff]
                        label = "P_Value Cutoff"
                        display = "FLOAT"
                        value = 0.05
                        range = [0, 1]
                        help = "Selects the cutoff."

Text box

                    [sra]
                    label = "Run Number"
                    display = "TEXTBOX"
                    value = "SRR519926"
                    regex = 'SRR\d+'
                    help = "Please provide SRR Run number"

Dropdown menu

                    [color]
                    label = "Select color"
                    display = "DROPDOWN"
                    choices = [ ["R","Red"], ["B","Blue"] ]
                    value = "R"
                    help = "Select a color of your choice"

Check box

                    [validate]
                    label = "Cross Validate"
                    display = "CHECKBOX"
                    value = true
                    help = "Apply cross validation on the results."

Upload Field

                    [file]
                    label = "Upload a file "
                    display = "UPLOAD"
                    help = "Upload a file to analyze"

Pick Data

                    [data]
                    label = "Pick data"
                    source = "PROJECT"
                    help = "Pick data from this project."

Radio Buttons

                    [species]
                    label = "Select species"
                    display = "RADIO"
                    choices = [ ["c","Cat"], ["d","Dog"] ]
                    value = "c"
                    help = "Select the species."

Name

Recipe display name

Identifier

Unique identifier for the recipe.

Authorization

Determines who can run a recipe.

Bacterial variant caller

This recipe contains the code presented by Torsten Seemann in the blog post titled [A Unix one-liner to call bacterial variants][code]

The code from the above website has been generalized a bit with additional utility.

This recipe provides code that:

1. Downloads a reference genome (default value: [NZ_CP008918][NZ_CP008918] *Pasteurella multocida strain ATCC 43137, complete genome*)
1. Downloads a sequencing run from SRA (default value: [SRR4124989][SRR4124989])
1. Generates statistics on the sequencing data
1. Runs `minimap2` to align the sequences and create a SAM output
1. Runs `bcftools mpileup` to generate the genotype likelihoods of each base
1. Runs `bcftools call` to filter for multiallelic variants only
1. Runs `bcftools norm` to normalize each variant to a standard form
1. Runs `bcftools filter` to remove variants with low quality (QUAL) or low coverage (DP)

[NZ_CP008918]: https://www.ncbi.nlm.nih.gov/nuccore/NZ_CP008918.1
[SRR4124989]: https://www.ncbi.nlm.nih.gov/sra/?term=SRR4124989

[code]: https://thegenomefactory.blogspot.com/2018/10/a-unix-one-liner-to-call-bacterial.html

A detailed explanation of what the recipe does (markdown OK).

Image :

Optional image for the recipe ( 500px Maximum ).

Rank:

Used to order recipes (optional).

Back