This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How "bases mapped (cigar)" value in samtools stats is calculated?

Hi:

I am trying to use samtools stats to estimate "bases mapped (cigar)" value from some target regions. But I can't figure out how the number is calcualted out. Below is a small test:

This is a sam file subset "test_subset.sam", contains only 15 reads that all covers posiotion chr22:31656041. All 15 reads are 34M in cigar.:

@HD     VN:1.3
@SQ     SN:chr1 LN:249250621
@SQ     SN:chr2 LN:243199373
@SQ     SN:chr3 LN:198022430
@SQ     SN:chr4 LN:191154276
@SQ     SN:chr5 LN:180915260
@SQ     SN:chr6 LN:171115067
@SQ     SN:chr7 LN:159138663
@SQ     SN:chr8 LN:146364022
@SQ     SN:chr9 LN:141213431
@SQ     SN:chr10        LN:135534747
@SQ     SN:chr11        LN:135006516
@SQ     SN:chr12        LN:133851895
@SQ     SN:chr13        LN:115169878
@SQ     SN:chr14        LN:107349540
@SQ     SN:chr15        LN:102531392
@SQ     SN:chr16        LN:90354753
@SQ     SN:chr17        LN:81195210
@SQ     SN:chr18        LN:78077248
@SQ     SN:chr19        LN:59128983
@SQ     SN:chr20        LN:63025520
@SQ     SN:chr21        LN:48129895
@SQ     SN:chr22        LN:51304566
@SQ     SN:chrX LN:155270560
@SQ     SN:chrM LN:16571
USI-EAS28:3:8:1631:168#0        16      chr22   31656010        37      34M     *       0       0       CGTGTCCGTGGAGAGTGCCTGCTCCAACTACGCC      \ZZNON\QTUYEGa``R^^``\^`a_`a`X^\`_     XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:34
USI-EAS28:1:82:1303:1490#0      16      chr22   31656011        37      34M     *       0       0       GTGTCGGTGGAGAGTGCCTGCTCCAACTACGCCA      BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB     XT:A:U  NM:i:1  X0:i:1  X1:i:0  XM:i:1  XO:i:0  XG:i:0  MD:Z:5C28
USI-EAS28:2:36:650:781#0        16      chr22   31656012        37      34M     *       0       0       TGTCCGTGGAGATTGCCTGCTCCAACTACGCCAC      BBBS]NUOWRKLPSDYUXXXZVVVM\ZN___ZT\     XT:A:U  NM:i:1  X0:i:1  X1:i:0  XM:i:1  XO:i:0  XG:i:0  MD:Z:12G21
USI-EAS28:3:71:1674:707#0       16      chr22   31656012        37      34M     *       0       0       TGTCCGTGGAGAGTGCCTGCTCCAACTACGCCAC      _]^]WY[a```a``\[UaaX`_aa`[`^aa`_a`     XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:34
HWI-EAS276:4:52:1678:656#0      16      chr22   31656012        37      34M     *       0       0       TGTCCGTGGAGAGTGCCTGCTCCAACTACGCCAC      Z_```_U_a\a\aZa`\P__O_`]S]JV`a``Y`     XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:34
HWI-EAS276:4:55:1651:892#0      16      chr22   31656012        37      34M     *       0       0       TGTCCGTGGAGAGTGCCTGCTCCAACTACGCCAC      [a\`__V\aTa_aaa`YL_aa^``WaR\aaaaZ_     XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:34
HWI-EAS276:4:57:1535:1224#0     16      chr22   31656014        37      34M     *       0       0       TCCGTGGAGAGTGCCTGCTCCAACTACGCCACCG      `aa_X_aaa`_Z_[\O]][[aa]a`\aaa_\`aa     XT:A:U  NM:i:1  X0:i:1  X1:i:0  XM:i:1  XO:i:0  XG:i:0  MD:Z:33A0
HWI-EAS276:4:6:1768:1347#0      0       chr22   31656017        37      34M     *       0       0       GTGGAGAGTGCCTGCTCCAACTACGCCACCACTG      a\\R``_^R\^_FVaNa_P\`\Wa^__M_\\a]^     XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:34
USI-EAS28:2:3:1230:1901#0       0       chr22   31656019        25      34M     *       0       0       GGAGAGTGCCTGCTCCAACNACGCCACCACTGAG      ]U]\``U]``a]aa`a[[[BBBBBBBBBBBBBBB     XT:A:U  NM:i:2  X0:i:1  X1:i:0  XM:i:2  XO:i:0  XG:i:0  MD:Z:19T12T1
USI-EAS28:3:93:186:495#0        16      chr22   31656020        37      34M     *       0       0       GAGAGTGCCTGCTCCAACTACGCCACCACTGTGC      BBBBBBBBBBBBBBB[Z_Y[_X_OH]__ZPMIV`     XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:34
USI-EAS28:3:55:1400:1957#0      0       chr22   31656021        37      34M     *       0       0       AGAGTGCCTGCTCCAACTACGCCACCACTGTGCA      _]_`]a`bW\a`b`Q\aaaa^^[^[^^U`]^ZSP     XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:34
HWI-EAS276:4:20:11:866#0        0       chr22   31656023        37      34M     *       0       0       AGTGCCTGCTCCAACTACGCCACCACTGTGCAAG      ]I]W^a_T^`aabb]bb^YW_aa`````\W[]]O     XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:34
USI-EAS28:3:56:764:965#0        16      chr22   31656029        37      34M     *       0       0       TGCTCCAACTACGCCACCACTGTGCAAGTGAAAG      ][a_a_a\^W`a`baZbbaTa]`abb^SaTbbb_     XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:34
USI-EAS28:2:7:707:1027#0        0       chr22   31656030        37      34M     *       0       0       GCTCCAACTACGCCACAACTGTGCAAGTGAAAGA      ZVWY`bb]R`^T]_aVD[```^_N`aS[[\]^ZB     XT:A:U  NM:i:1  X0:i:1  X1:i:0  XM:i:1  XO:i:0  XG:i:0  MD:Z:16C17
USI-EAS28:2:38:109:1946#0       16      chr22   31656031        37      34M     *       0       0       CTCCAACTACGCCACCACTGTGCAAGTGAAAGAG      \X`_Z``^UaL`aUaa^baUaaaZ__aI_ababa     XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:34

Then the bed file "test.bed" is like below, only one record:

chr22   31656041    31656041

Then I run:

samtools stats --target-regions ./test.bed ./test_subset.sam

Below is the output of the above command:

# This file was produced by samtools stats (1.22.1+htslib-1.22.1) and can be plotted using plot-bamstats
# This file contains statistics for all reads.
# The command line was:  stats --target-regions ./test.bed ./test_subset.sam
# CHK, Checksum [2]Read Names   [3]Sequences    [4]Qualities
# CHK, CRC32 of reads which passed filtering followed by addition (32bit overflow)
CHK 26fde71f    57af31e8    868ff8b1
# Summary Numbers. Use `grep ^SN | cut -f 2-` to extract this part.
SN  raw total sequences:    15  # excluding supplementary and secondary reads
SN  filtered sequences: 0
SN  sequences:  15
SN  is sorted:  1   # sorted by coordinate
SN  1st fragments:  15
SN  last fragments: 0
SN  reads mapped:   15
SN  reads mapped and paired:    0   # paired-end technology bit set + both mates mapped
SN  reads unmapped: 0
SN  reads properly paired:  0   # proper-pair bit set
SN  reads paired:   0   # paired-end technology bit set
SN  reads duplicated:   0   # PCR or optical duplicate bit set
SN  reads MQ0:  0   # mapped and MQ=0
SN  reads QC failed:    0
SN  non-primary alignments: 0
SN  supplementary alignments:   0
SN  total length:   510 # ignores clipping
SN  total first fragment length:    510 # ignores clipping
SN  total last fragment length: 0   # ignores clipping
SN  bases mapped:   510 # ignores clipping
SN  bases mapped (cigar):   168 # more accurate
SN  bases trimmed:  0
SN  bases duplicated:   0
SN  mismatches: 6   # from NM fields
SN  error rate: 3.571429e-02    # mismatches / bases mapped (cigar)
SN  average length: 34
SN  average first fragment length:  34
SN  average last fragment length:   0
SN  maximum length: 34
SN  maximum first fragment length:  34
SN  maximum last fragment length:   0
SN  average quality:    55.6
SN  insert size average:    0.0
SN  insert size standard deviation: 0.0
SN  inward oriented pairs:  0
SN  outward oriented pairs: 0
SN  pairs with other orientation:   0
SN  pairs on different chromosomes: 0
SN  percentage of properly paired reads (%):    0.0
SN  bases inside the target:    1
SN  percentage of target genome with coverage > 0 (%):  100.00

The "bases mapped (cigar)" value is 168? How this is calculated?

Below is the IGV plot, it only contains these 15 reads:

enter image description here

samtools

0 answers

No answers yet.

Log in to answer this question.