<p>Examining the mpileup output for a single read on its own is a good way to figure out why bases are not appearing as expected.</p>

<p>So if you suspect that the <code>cigar=18S58M</code> read is the missing one (you are correct), prepare a SAM file containing just that read:</p>

<pre><code>@SQ SN:chr1 LN:248956422
NS500355:NS500355:HHYN5BGXB:1:13209:10692:11045 83  chr1    569888  53  18S58M  =   569888  -58 CGTGTGCTCTTCCGATCTCTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAG    EE&lt;A/EE//AE/AA///AAE///EEEEE//6///E//EE/6EE/EEEEEEEA/A/A/E/E//AEAEEEEE6A6AAA    NM:i:0  MD:Z:58 MC:Z:59M17S AS:i:58 XS:i:53RG:Z:../H2O  XA:Z:chrM,-9340,18S58M,1;
</code></pre>

<p>and run it through <code>mpileup</code>:</p>

<pre><code>$ samtools mpileup -Q0 -B oneread.sam
</code></pre>

<p>This will show you that this read's base quality dips below 20 at the position you are interested in. Thus, as finswimmer suspected, your <code>-Q20</code> is removing it at this position.</p>

CC BY-SA 4.0
2019-08-22T09:16:48.480Z
2019-08-22T09:16:48.480Z
John Marshall
<p>Why do you believe it's the read with <code>cigar=18S58M</code>? The mapping qualities given in sam files are different to those in the mpileup output. I guess samtools is doing some recalculation.</p>

<p>You are using the parameters <code>-q20</code> <code>-Q20</code>. So you skip all reads with mapping quality &lt;20. This will result in 7 reads got to mpileup. Furthermore you are skipping bases with base quality &lt;20. Have a look at the base qualities for the read starting at position <code>569904</code>. If I count it correct you end up with a quality of <code>/</code> at position <code>569930</code> which means <code>14</code>. So I guess this read is dropped.</p>

<p><strike>BTW: <code>samtools mpileup</code> is deprecated</strike> (only when creating bcf/vcf) <strike> and should be replaced by <code>bcftools mpileup</code>.</strike> </p>

CC BY-SA 4.0
2019-08-22T13:19:39.947Z
Thank you, I missed the '/' quality for a 'E'. Now I understand.
aerijman
CC BY-SA 4.0
2019-08-23T09:28:19.370Z
The `mpileup` subcommand has been split between the two commands. Textual “mpileup output” stays in samtools (and always will), but VCF/BCF output has moved to bcftools. This groups the latter alongside the other VCF-related facilities and in particular means that both halves of the calling pipeline `bcftools mpileup … | bcftools call …` are in the same program, which simplifies updates and improvements to this pair of commands that work together.
John Marshall
CC BY-SA 4.0
2019-08-22T04:20:27.143Z
2019-08-22T11:02:40.700Z
2019-08-22T11:02:40.700Z
finswimmer
<p>I run the command:</p>

<pre><code>samtools mpileup -O -s -q20-B -Q20 -f hg19.fa -r chr1:569929-569931 myFile.bam
</code></pre>

<p>and get: </p>

<pre><code>chr1    569929  G   7   ...,,., EEEEEEE NTTVTTU 53,48,42,60,30,29,27
chr1    569930  C   6   ...,.,  EEAEEE  NTTTTU  54,49,43,31,30,28
chr1    569931  G   6   ...,.,  EEEEEE  NTTTTU  55,50,44,32,31,29
</code></pre>

<p>I am interested in <strong>position 569930</strong>.</p>

<p>I inspect the bam with:</p>

<pre><code>samtools view myFile.bam chr1:569929-569931
</code></pre>

<p>and I see:</p>

<pre><code>NS500355:NS500355:HHYN5BGXB:2:21311:2303:9255   99  chr1    569869  0   76M =   569900  107 CTCCATAACGCTCCTCATACTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGA    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE&lt;EEEEEEEAEEEEEEE    NM:i:2  MD:Z:5C3C66 MC:Z:76M    AS:i:66 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9321,76M,1;
NS500355:NS500355:HHYN5BGXB:4:22603:25857:2799  99  chr1    569877  45  76M =   569896  95  CGCTCCTCATACTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCAC    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE    NM:i:1  MD:Z:1C74   MC:Z:76M    AS:i:74 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9329,76M,1;
NS500355:NS500355:HHYN5BGXB:1:12105:11868:12011 99  chr1    569882  51  76M =   569900  94  CTCATACTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACC    /AAAAEEEEEAEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEA    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9334,76M,1;
NS500355:NS500355:HHYN5BGXB:2:13301:2681:6114   99  chr1    569888  9   76M =   569952  140 CTAGGCATACTAACCAACACACTAACAATATACCAATGATGGAGCGATGTAACACGAGAAAGCACATACCAACGCC    AAAA/E/EE//AAEEEEE/EAAEAA//EEAEEE/E/EEEE/////A/EE/EE/EE/E//EE/EEEE/E/AEA####    NM:i:4  MD:Z:6C19C15C29G3   MC:Z:76M    AS:i:57 XS:i:52 RG:Z:../H2O XA:Z:chrM,+9340,76M,5;
NS500355:NS500355:HHYN5BGXB:2:23110:10125:9334  99  chr1    569888  51  76M =   569904  92  CTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCC    AAA/AEE6EEEE6E///E/A/AAE/E/AEEEEEEE/EE/EEEAEEEEAEE//EEEEEEAEEEAEEAE/AEE/AEE&lt;    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9340,76M,1;
NS500355:NS500355:HHYN5BGXB:1:13209:10692:11045 83  chr1    569888  53  18S58M  =   569888  -58 CGTGTGCTCTTCCGATCTCTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAG    EE&lt;A/EE//AE/AA///AAE///EEEEE//6///E//EE/6EE/EEEEEEEA/A/A/E/E//AEAEEEEE6A6AAA    NM:i:0  MD:Z:58 MC:Z:59M17S AS:i:58 XS:i:53 RG:Z:../H2O XA:Z:chrM,-9340,18S58M,1;
NS500355:NS500355:HHYN5BGXB:1:21304:26480:15673 83  chr1    569890  11  76M =   569796  -170    AGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCAC    AEEEEEEEAAEEEEEEEEEEEE/EEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAAAAA    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,-9342,76M,1;
NS500355:NS500355:HHYN5BGXB:4:11602:21237:3028  83  chr1    569900  51  76M =   569900  -76 ACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCAC    EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAAAAA    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,-9352,76M,1;
NS500355:NS500355:HHYN5BGXB:3:13504:4369:8812   99  chr1    569901  51  76M =   569901  76  CCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACC    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE6EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9353,76M,1;
NS500355:NS500355:HHYN5BGXB:2:11203:19774:18704 99  chr1    569904  11  76M =   569952  124 ACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTGT    6AAAAEEEEEEEEEEEEEAEEEEEEEEEEEEEEEE6EEEEAEEEAEEEEEEEEEEEEEEEEEEAE/E6EEAEEEEA    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9356,76M,1;
NS500355:NS500355:HHYN5BGXB:2:11309:13954:17567 99  chr1    569905  11  76M =   569972  143 CACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTGTC    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEE    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9357,76M,1;
NS500355:NS500355:HHYN5BGXB:1:21202:13252:14049 99  chr1    569905  11  76M =   569927  98  CACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTGTC    AAAAAEE6EEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEE6EEEEEEEEEEEEEEEEE/EEEE&lt;EEEEEEEEEEE    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9357,76M,1;
NS500355:NS500355:HHYN5BGXB:1:22311:21964:12072 83  chr1    569908  52  5S71M   =   569908  -71 GATCTACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTG    &lt;AEEEA6EEE6EEEEEEEEEEEEEE/EEEEEEEEEEEEEEEEE/EEEEEEEEEEEEEEEEEEAAEEEAEEEAA6AA    NM:i:0  MD:Z:71 MC:Z:71M5S  AS:i:71 XS:i:66 RG:Z:../H2O XA:Z:chrM,-9360,5S71M,1;
NS500355:NS500355:HHYN5BGXB:4:12602:7858:19418  99  chr1    569909  15  76M =   569917  84  CTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTGTCCAAA    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEE/AAEEEEEEEEEEEEEEE    NM:i:1  MD:Z:74G1   MC:Z:76M    AS:i:74 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9361,76M,1;
NS500355:NS500355:HHYN5BGXB:2:12205:19456:7598  83  chr1    569927  0   76M =   569907  -96 TGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTGTCCAAAAAGGCCTTCGATACGGGA    EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAA/EEEEEEEEEAAAAA    NM:i:1  MD:Z:56G19  MC:Z:76M    AS:i:71 XS:i:76 RG:Z:../H2O XA:Z:chrM,-9379,76M,0;
</code></pre>

<p>Based on mapq and positions (respect to 569930) I see that the read with <code>cigar=18S58M</code> is skipped by mpileup.  </p>

<p><strong>WHY?</strong> The soft clip is way far from position chr1:569930 </p>

<samtools><mpileup>
samtools mpileup skipping read
# samtools mpileup skipping read

(Asked by: username_0 on 2019-08-22)

I run the command:
samtools mpileup -O -s -q20-B -Q20 -f hg19.fa -r chr1:569929-569931 myFile.bam

and get: 
chr1    569929  G   7   ...,,., EEEEEEE NTTVTTU 53,48,42,60,30,29,27
chr1    569930  C   6   ...,.,  EEAEEE  NTTTTU  54,49,43,31,30,28
chr1    569931  G   6   ...,.,  EEEEEE  NTTTTU  55,50,44,32,31,29

I am interested in position 569930.
I inspect the bam with:
samtools view myFile.bam chr1:569929-569931

and I see:
NS500355:NS500355:HHYN5BGXB:2:21311:2303:9255   99  chr1    569869  0   76M =   569900  107 CTCCATAACGCTCCTCATACTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGA    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
Based on mapq and positions (respect to 569930) I see that the read with cigar=18S58M is skipped by mpileup.  
WHY? The soft clip is way far from position chr1:569930

**Tags:** samtools, mpileup

---

## 2 answers

### Answer by username_1 on 2019-08-22. Score: 3 | Accepted answer

Examining the mpileup output for a single read on its own is a good way to figure out why bases are not appearing as expected.
So if you suspect that the cigar=18S58M read is the missing one (you are correct), prepare a SAM file containing just that read:
@SQ SN:chr1 LN:248956422
NS500355:NS500355:HHYN5BGXB:1:13209:10692:11045 83  chr1    569888  53  18S58M  =   569888  -58 CGTGTGCTCTTCCGATCTCTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAG    EE
and run it through mpileup:
$ samtools mpileup -Q0 -B oneread.sam

This will show you that this read's base quality dips below 20 at the position you are interested in. Thus, as finswimmer suspected, your -Q20 is removing it at this position.


---

### Answer by username_2 on 2019-08-22. Score: 3

Why do you believe it's the read with `cigar=18S58M`? The mapping qualities given in sam files are different to those in the mpileup output. I guess samtools is doing some recalculation.

You are using the parameters `-q20` `-Q20`. So you skip all reads with mapping quality <20. This will result in 7 reads got to mpileup. Furthermore you are skipping bases with base quality <20. Have a look at the base qualities for the read starting at position `569904`. If I count it correct you end up with a quality of `/` at position `569930` which means `14`. So I guess this read is dropped.

BTW: `samtools mpileup` is deprecated (only when creating bcf/vcf)  and should be replaced by `bcftools mpileup`.

> **Comments:**
> **username_0**: Thank you, I missed the '/' quality for a 'E'. Now I understand.  
> **username_1**: The `mpileup` subcommand has been split between the two commands. Textual “mpileup output” stays in samtools (and always will), but VCF/BCF output has moved to bcftools. This groups the latter alongside the other VCF-related facilities and in particular means that both halves of the calling pipeline `bcftools mpileup … | bcftools call …` are in the same program, which simplifies updates and improvements to this pair of commands that work together.  

---

