This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Mpileup With A Region File

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

samtools variant chip-seq mpileup
• 15,041 views • link
updated by Devon Ryan Devon Ryan • written by Dataminer Dataminer

1 answer

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
• 222 views • link
written by Devon Ryan Devon Ryan

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

• 0 views • link
written by Dataminer Dataminer

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.

• 368 views • link
written by Devon Ryan Devon Ryan

Log in to answer this question.

More posts like this