Mpileup With A Region File

Question

[Dataminer]

Hi!

I have two sorted .bam files and I have a set of regions (chr str stp) in a bed format (3 columns).

I would like to check for variants only in these regions rather through out the genome.

What I know, is that, following commands can be used

samtools mpileup -uD -r 2L:100,000-150,000 -f hg19.fa file1.sorted.bam file2.sorted.bam | bcftools view -bvcg - > RAL_samtools.raw.bcf
bcftools view RAL_samtools.raw.bcf | vcfutils.pl varFilter -D100 > RAL_samtools.vcf

but how can I put a region file and also filter the results based on regions that have more than 5X coverage on an average.

Any help is welcome.

Thank you

Answer 1

[Devon Ryan]

Just use the -l regions.bed option in samtools mpileup:

samtools mpileup -uDl regions.bed -f hg19.fa file1.sorted.bam file2.sorted.bam | bcftools view -bvcg - > RAL_samtools.raw.bcf
bcftools view RAL_samtools.raw.bcf | vcfutils.pl varFilter -D100 > RAL_samtools.vcf

[Dataminer]

Thank you, the command is running. One more thing, is it possible to add more cores to it?

[Devon Ryan]

The base samtools mpileup command is single threaded, so no. Having said that, you could just split the regions.bed file into multiple chunks and then run multiple instances of the command at once (perhaps merging the VCF files at the end). I'd be surprised if no one's yet written a little wrapper script to do that.