samtools mpileup skipping read

Question

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<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<    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<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<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    <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;

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

score 3 · Answer 1 · edited Aug 22 '19 at 11:02

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.~~

edited Aug 22 '19 at 11:02

terdon

10,071
5
22
48

answered Aug 22 '19 at 04:20

finswimmer

1,342
7
15

Thank you, I missed the '/' quality for a 'E'. Now I understand. – aerijman Aug 22 '19 at 13:19
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 Aug 23 '19 at 09:28

score 3 · Accepted Answer · answered Aug 22 '19 at 09:16

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<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;

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.

samtools mpileup skipping read

2 Answers2