This is a test version of Biostars. For the public version, visit https://www.biostars.org.
filtering bam file

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

rna-seq next-gen alignment
• 2,799 views • link
updated by Ram Ram • written by schelarina • 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
• 0 views • link
updated by Ram Ram • written by schelarina • schelarina

1 answer

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

• 0 views • link
written by michael.ante michael.ante

Log in to answer this question.

More posts like this