Question: In [this thread][1] I saw a line that would run for multiple threads on a computed

parallel 'bcftools mpileup -Ou -f ref.fa -r {} in.bam | bcftools call -m -v -Oz -o {}.vcf.gz' ::: {1..22} X Y

Using the -r flag "-r, --regions chr|chr:pos|chr:from-to|chr:from-[,…]" it is possible to specify a region to look for. Then from there, since my reference genome is not separated in chromosomes (only scaffold level => [see the list here][2])

I get that the command is to divide the task between the 24 chromosomes ({1..22} X Y), but how would it be done if I don't have that for my genome? Also, I'm not sure to understand what the `{}.vcf.gz` or simply the `{}` in the output file. Should I put it in the ``? Finally, the --jobs flag in `parallel` will determine the number of scores, but in this case, since it's not flagged, how do I know how many cores it's using (if I have access to 16 cores max)?

This is how I think it should look like, but how do I specify the division of the task (is it better to list the scaffolds (where there are 1000s of them)?). Do I need to put a `{}` in the output file?

parallel 'bcftools mpileup -C 50 -E -f --regions {} -b > ' ::: {SCAFFOLDS??}

[1]: https://github.com/samtools/bcftools/issues/949
[2]: https://genome.ucsc.edu/cgi-bin/hgTracks?chromInfoPage=&hgsid=794599813_cKcep24UuGc7AxHiB7JudsBFFrra

Answer: The trick here is to use the random access (so the BAM index) to pull out only the reads you want from the BAM file in order to save time. This code snippet will parallelize the mpileup per scaffold given you have a BED file with the chr-start-end. **I did only do a minimal testing** and the script looks a bit bulky but should actually work.

The idea is the following: The `samtools view` will use the region information to only extract the overlapping reads and send them to `mpileup`. Therefore `mpileup` only processes exactly these reads, so in your case the contigs one-by-one, but parallelized via GNU parallel. Eventually, the intermediate BCFs are concatenated.

YOUR_BAM="your.bam"
CONTIGS="your_file_with_the_contigs.bed"
NUMBER_OF_JOBS = 8
## Index bam file:
samtools index ${YOUR_BAM}

## A function for bcftools where samtools view extracts only the indicated regions:
function MPILEUP {

BAM=$1
REGION=$2
samtools view -bu -r ${REGION} $BAM \
| bcftools mpileup (options...) -o output_fromBAM_${BAM%.bam}_region_${REGION}.bcf -O b -

}; export -f MPILEUP

cat $CONTIGS | awk '{print $1":"$2"-"$3}' | parallel -j ${NUMBER_OF_JOBS} "MPILEUP ${YOUR_BAM} {}"

## Concatenate everything:
bcftools concat -O b output_fromBAM_${YOUR_BAM%.bam}_region_*.bcf > concat_fromBAM_${YOUR_BAM%.bam}.bcf

## remove temporary files:
rm output_fromBAM_${YOUR_BAM%.bam}_region_*.bcf

Here `NUMBER_OF_JOBS` would be 8 so 8 parallel cores working on the task. I would check with `top` if your harddrive can handle this in terms of I/O bottleneck. This works well on a HPC, not so sure on a workstation with HDD.
Also, since `mpileup` will read BAM from `stdin`the sample name in the 10th column of the VCF/BCF will be `-` so you would need to check if mpielup can handle an external sample name to feed in via the command line, or you have to `awk` it in manually.