# Mapping statistics from bam file using bbtools and sambamba

Below are the statistics for RNA-seq mapped and unmapped paired-end reads to rice genome using reformat.sh from bbtools on bam files. It gives 77% mapped and 5% unmapped, what about the remaining 18% reads? How I can get information for all reads?

Command;
reformat.sh in=Leaf_T1_F_R10_S1_L001.bam out=Leaf_T1_F_R10_S1_L001.mapped.bam, mappedonly

Unspecified format for output Leaf_T1_F_R10_S1_L001.mapped.bam,; defaulting to fastq.
Found sambamba.
Input is being processed as unpaired
Output:                 52020538 reads (77.10%) 7208830835 bases (77.41%)

Time:                         171.212 seconds.
Bases Processed:       9312m 54.39m bases/sec

Command;
reformat.sh in=Leaf_T1_F_R10_S1_L001.bam out=Leaf_T1_F_R10_S1_L001.unmapped.bam, unmappedonly

Unspecified format for output Leaf_T1_F_R10_S1_L001.unmapped.bam,; defaulting to fastq.
Found sambamba.
Input is being processed as unpaired
Output:                 3435326 reads (5.09%) 465456267 bases (5.00%)

Time:                         142.978 seconds.
Bases Processed:       9312m 65.13m bases/sec


With sambamba; Can I say 94.91% reads are mapped and 5.09% are unmapped?

67471075 + 0 in total (QC-passed reads + QC-failed reads)
12015211 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
64035749 + 0 mapped (94.91%:N/A)
55455864 + 0 paired in sequencing
51006046 + 0 properly paired (91.98%:N/A)
51419222 + 0 with itself and mate mapped
601316 + 0 singletons (1.08%:N/A)
263816 + 0 with mate mapped to a different chr
237111 + 0 with mate mapped to a different chr (mapQ>=5)


Try samtools flagstat or sambamba flagstat to get similar type of information about your bam file. It may be easier to interpret. Both of these tools can be installed using conda if needed.

EDIT:

Yes, the results of sambamba flagstat you showed in your edited post indicate exactly what you said: 94.91% mapped and 100% - 94.91% = 5.09% unmapped:

64035749 + 0 mapped (94.91%:N/A)

• Hi Timur, Thanks, I have added the results from sambamba as well, I think it gives the accurate result. Could you please respond to query related to sambamba result? Sep 27 '20 at 17:58
• @bioinfonext Edited the answer. Sep 27 '20 at 18:50

More of a comment than an answer... In my opinion the terminology is a bit confusing since the term reads actually refers to a more generic records.

Look at this test file:

samtools view test.sam | cut -f 1-4
HSQ9103:404:C6F0VANXX:1:2208:4363:50381 0   chr7    5529094
HSQ9103:404:C6F0VANXX:1:2208:4363:50381 0   chr7    5529194
HSQ9103:404:C6F0VANXX:1:2208:4363:50381 0   chr7    5529194
HSQ9103:404:C6F0VANXX:1:2208:4363:50381 0   chr7    5529194
HSQ9103:403:C6F0HANXX:8:2215:4792:20065 4   chr7    5529871


I would say that here there are 5 records (of which 4 are alignments) but only 2 reads.

samtools flagstat reports 80% mapping rate but maybe I would say this is 50%.

samtools flagstat test.sam