Recipe View

How to generate variant calls for multiple samples

This recipe demonstrates variant call generation for several samples at once.

1 result • updated 5.1 years ago by Istvan Albert

This recipe demonstrates variant call generation for several samples at once.

The resulting file will be a VCF that stores cumulative information on all samples.

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

# NCBI Bioproject ID.
ID=PRJNA257197

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

# The number of samples to process.
N=5

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

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

# Obtain the reference genome.
efetch -db nuccore -id $ACC -format fasta > $REF

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

# Index the FASTA file for IGV.
samtools faidx $REF >> 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  > srr.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 srr.txt | parallel "fastq-dump -X $LIMIT -O reads --split-files {} >> log.txt"

# Generate alignments
mkdir -p bam

# Run the bwa mem processes in parallel for all samples.
cat srr.txt | parallel "bwa mem $REF reads/{}_1.fastq reads/{}_2.fastq 2>> log.txt | samtools sort > bam/{}.bam 2>> log.txt"

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

# Call variants on all BAM files at the same time.
bcftools mpileup -Ou -f $REF  bam/*.bam 2>> log.txt | bcftools call --ploidy 1 -vm -Ou |  bcftools norm -Ov -f $REF -d all - > variants.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"

[proj]
label = "The NCBI project run number"
display = "TEXTBOX"
value = "PRJNA257197"
regex = "^\\w{1,20}$"
help = "The bioproject that stores the data."

[N]
label = "How many samples to process"
display = "INTEGER"
value = 5
range = [ 1, 100,]
help = "How many SRR runs to process in the bioproject."

[settings]
id = 98
recipe_uid = "recipe-multisample"
uid = "recipe-multisample"
name = "How to generate variant calls for multiple samples"
template = "How_to_generate_variant_calls_for_multiple_samples_bio-data-analysis_98.sh"
image = "How_to_generate_variant_calls_for_multiple_samples_bio-data-analysis_98.png"
project_uid = "2a988f7f"
url = "http://localhost8000"
help = "This recipe demonstrates variant call generation for several samples at once.\n\nThe resulting file will be a VCF that stores cumulative information on all samples.\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