Question: help with calling invariant sites in parallel: scaffolds not present in all bam files are skipped
gravatar for truebeliever24
11 weeks ago by
truebeliever2420 wrote:


I'm using bcftools to call invariant sites in parallel. It appears that only scaffolds that are present across all bam files are included, while any scaffold missing in an individual in the dataset is omitted. What can I do to ensure all scaffolds are included in the invariant site calling? Thanks for your help-I am using the following tutorial:

Here's my code:

split chromosomes from bam files

function bam_chromosomes { samtools idxstats $1 | cut -f 1 | grep -v '*' }

genome=/path_to_reference export genome

apply operations to each scaffold individually in parallel

function parallel_call { bcftools mpileup --fasta-ref ${genome} --regions $2 --output-type u --bam-list $1 | \ bcftools call --multiallelic-caller --variants-only --output-type v - > ${1/.bam/}.$2.vcf }

export -f parallel_call

chrom_set=bam_chromosomes test.bam

parallel --verbose -j 25 parallel_call $1 ::: ${chrom_set}

scaffold bcftools bam invariant • 175 views
ADD COMMENTlink written 11 weeks ago by truebeliever2420

while parallel is a great tool, you'd better use a workflow manager like nextflow or snakemake.

ADD REPLYlink written 11 weeks ago by Pierre Lindenbaum134k

see also : Generating shell scripts

ADD REPLYlink written 11 weeks ago by Pierre Lindenbaum134k
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 2100 users visited in the last hour