This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Extract flanking regions from bam files based on coverage threshold

Hi all,

I have bam files for 10 individuals mapped onto a reference genome. I'd like to extract SNPs and their flanking regions (50 bp upstream and 50 bp downstream) based on a coverage threshold. F.e. minimum coverage of 20 for all 10 bam files combined. I'd also prefer to have not more than 1 polymorhpic site in the flanking regions.

Any suggestions on how to approach this? I assume the first step would be to merge the bam files?

Thanks!

snp flanking region coverage bam mapping
• 2,114 views • link
written by laura.bertola • laura.bertola

One way to do this would be to use samtools mpileup supplying a list of SNP regions as an option, then filter the output based on your minimum coverage requirement. Related: A: Mpileup With A Region File

• 0 views • link
written by mark.ziemann mark.ziemann

Thanks for the suggestion! If I understand correctly, I'd use samtools mpileup to generate a list of SNPs first, then convert this to a list of regions and repeat the mpileup step? This makes sense, although I'd thought there may be a more straight forward way of achieving this.

• 0 views • link
written by laura.bertola • laura.bertola

Most SNP callers like GATK, VarScan, BCFtools, will provide a coverage for each SNP that you can filter on. After filtering for SNPs that satisfy your threshold, you can use bedtools slop to expand them by 50 bases and supply these regions to samtools mpileup.

• 0 views • link
written by mark.ziemann mark.ziemann

Thanks! I've been looking into this and I think it is working!

• 0 views • link
written by laura.bertola • laura.bertola

0 answers

No answers yet.

Log in to answer this question.

More posts like this