filtering bam file

Question

[schelarina]

Hello,

I want to filter a bam file based on the number of reads mapped to each position in the genome.

For example I want to keep only the positions that are covered by >=40 stacked reads (with different length) but starting at the same position

Thanks

[schelarina]

I tried

samtools view -h -b file.bam | awk '{a[$4]++}END{for (i in a)if (a[I]>=40)print $0;}' > file.bam

Answer 1

[michael.ante]

Hi,

you can use the bedTools' genomeCoverageBed bamToBed and mergeBed to get all start-positions' coverage; filter all positions with coverage >= 40, and build a bed file. This bed-file can be used as input for samtools view (-L bedfile) intersectBed with the first bed-file and extract the read-names.

Cheers,

Michael