how to count un-mapped reads at specific position using samtools?

Hello all.

I want to get count of unmapped reads at specific position. After searching the google, I am able to do that using below command with -f option.

samtools view -f 4 [bam]

But this is not what I am looking for. It's because I wanted to find the count of unmapped reads at the specific positions , not all over the sequences, therefore I added other options and I finally got the unmapped reads at the specific positions as shown in command below.

samtools view -f 4 -L /DATA2/sclee1/URC_WES/ml/allel_info/ms/ms_AYA06_mpileup_modified.bed 01U_T_Filtered_Sorted_Makrdup_Readgroup_Realigned_Recal.bam

However, I have to do additional works to count reads at each site because command above just outputs a list of reads information (not total count). Therefore, it is a little nuisance and I didn't know how to do it in more sophisticated ways.

Could you give me any ideas or solutions for doing this?

sequencing genome sequence

• 2,850 views • link

updated 3.5 years ago by Ram

• written 11.2 years ago by mangfu100

1 answer

Unmapped alignments may or may not have a position associated. Even when they do have such a position, it's not necessarily meaningful (for paired-end reads, it's often the position of a read's mate if it's mapped).

Anyway, for counting purely with samtools there's the -c option. This will not, however, give you per-region counts. For that you'd need to use something like htseq-count or featureCounts. Of course, these will typically exclude unmapped reads for the reasons I mentioned in the previous paragraph.

In short, you might want to say what your actual biological goal is here.

• 0 views • link

updated 3.5 years ago by Ram

• written 11.2 years ago by Devon Ryan

Thank you for response.

Before answering, I am using the paired end reads.

My goal is that I would like to calculate the ratio between uniquely mapped reads and unmapped reads at specific positions. For example, at chr1:1000, there are total # of uniquely mapped reads are 50 while unmapped reads are 5, then the ratio should be 5/50 = 0.1

Anyway, I got to know that by using samtools -f options, I am able to get the unmapped reads...but I don't know how to check the reads are uniquely mapped or not.

Could you give me any advice for me?

• 0 views • link

updated 3.5 years ago by Ram

• written 11.2 years ago by mangfu100

There are multiple definitions of "uniquely mapped", you'll have to figure out a definition that makes sense for you.

Again, please state the biological goal. The ratio of mapped to unmapped alignments in a region is not a biological goal, it's a questionable method.

• 0 views • link

written 11.2 years ago by Devon Ryan