count reads with samtools view

Question

[schelarina]

Hello,

I want to count reads starting at a specific position e.g position 10000 in chromosome 1, output 50 reads start at position 10000

Here is the command that I run and gave me the correct output that I need:

samtools view -F 16 file.bam | grep chromosome_1 | awk '{print $4}' | grep -w 10000 | wc -l

I have a list of positions in a file.txt and for every position I would like to count how many reads start at that position and save the output in new file in this way: position 10000 counts 50

How can I do this in automatic way without entering the command for each position?

Thanks!

Answer 1

[tiago211287]

you can use in your txt, assuming there is one number per line:

   samtools view -F 16 file.bam | grep chromosome_1 | awk '{print $4}' > positions.txt
   sort positions.txt | uniq -c > output.txt

Then you can use grep again to look for some specific position.

grep 1000 output.txt

Answer 2

[Pierre Lindenbaum]

put your positions in a bed file, use the BAM index and use awk to count the lines.

sed -e 's/\t/:/' -e 's/\t/-/' input.bed |\
read F
do
   samtools view -F 4 - "${F}" |\
    awk -v "p=$F" 'BEGIN{n=0;split(p,a,/\-/);} {if($4==a[2]) n++;} END{ printf("%s\t%d\n",p,n);}'
done