This is a test version of Biostars. For the public version, visit https://www.biostars.org.
count reads with samtools view

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!

rna-seq next-gen sequence alignment
• 3,925 views • link
updated by Ram Ram • written by schelarina • schelarina

2 answers

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
• 0 views • link
written by tiago211287 tiago211287
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
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Log in to answer this question.

More posts like this