Counting number of duplicated reads for each gene or exon

Question

[arta]

Hey,

I am analyzing RNA-seq data and I am interested in duplicated reads. I know I can count number of duplicated reads overall using picard markDuplicates.

samtools view -f 1024 dedup_reads.bam | wc -l

But i am interesting in the distribution of these duplicated. I have both sam and bam files.

Here is simple what i would to have

               # total reads              #duplicated reads          
   Gene1         30                                 10
   Gene2        100                                20
   Gene3         20                                 0

I googled but i couldn't find any tools, is there any tools, softwares or packages? Or should I implement myself.

Answer 1

[Pierre Lindenbaum]

awk '{printf("%s:%d-%d\n",$1,int($2)+1,$3);}' input.bed  | while read F; 
do echo -n "$F " && samtools view -F 4 -c  input.bam $F | tr -d '\n' && echo -n " " &&  samtools view -f 1024 -F 4 -c  input.bam $F
done

updated : oopss added $F after each "samtools view"

[arta]

Great, thanks Pierre !!!