question on SAMtools consensus calling

[orionzhou]

hey guys, I'm really confused how samtools calls the consensus sequence:

I've generated the sorted and indexed *.bam file. First I generated the pileup of a 100bp region:

samtools view -u 07_bam/HM020.bam chr1:6,269,001-6,269,100 | samtools pileup -f genome.fa -

I'm pasting the interested bases here:

chr1    6269058 A       11      ,,.,.,,,,,.     FEHHH5BHH%H
chr1    6269059 A       11      ,,.,.,,,,,.     F:FHH8@HH%H
chr1    6269060 T       11      ,,.,.,,,,,.     F9HBH48IC%H
chr1    [COLOR="Red"]6269061[/COLOR] T       13      ccCcCcccccC^]C^]c       [COLOR="RoyalBlue"]FEIHH<=HH%HD@[/COLOR]
chr1    [COLOR="Red"]6269062[/COLOR] C       13      ttTtTtttttTTt   [COLOR="RoyalBlue"]EFHDH;4GF@HD9[/COLOR]
chr1    6269063 G       15      ,,.,.,,,,,..,^>,^>,     H:GBH@:HF>H.>?E
chr1    6269064 T       15      ,,.,.,,,,,..,,, H;H=I:<IDAHD4;?
chr1    6269065 C       15      ,,.,.,,,,,..,,, HEHFH/<HGBHB7?C
chr1    6269066 C       15      ,,.,.,,,,,..,,, HADFH>@HH8HF>?C

Obviously positions 6269061 and 6269062 are two SNPs that should be called, and here's the command I entered for consensus calling:

samtools view -u 07_bam/HM020.bam chr1:6,269,001-6,269,100 | samtools pileup -f genome.fa [COLOR="Red"]-c[/COLOR] -

which produces:

chr1    6269058 A       A       57      0       50      11      ,,.,.,,,,,.     F@HHH5BHH%H
chr1    6269059 A       A       57      0       50      11      ,,.,.,,,,,.     F:FHH8@HH%H
chr1    6269060 T       N       0       0       0       11      ,,.,.,,,,,.     &&&&&),&&%&
chr1    [COLOR="Red"]6269061[/COLOR] T       N       0       0       0       13      ccCcCcccccC^]C^]c       [COLOR="RoyalBlue"]!!!!!$(!!%!!![/COLOR]
chr1    [COLOR="Red"]6269062[/COLOR] C       C       15      0       15      13      ttTtTtttttTTt   [COLOR="RoyalBlue"]!!!!!$(!!0!!![/COLOR]
chr1    6269063 G       G       72      0       48      15      ,,.,.,,,,,..,^>,^>,     H:GBH@:HF>H.<?E
chr1    6269064 T       T       72      0       48      15      ,,.,.,,,,,..,,, H;H=I:<IDAHD4;?
chr1    6269065 C       C       72      0       48      15      ,,.,.,,,,,..,,, HEHFH/<HGBHB7?C
chr1    6269066 C       C       72      0       48      15      ,,.,.,,,,,..,,, HADFH>@HH8HF>?C

It is the same BAM file I'm using, why are the read qualities at these positions differ so much?

PS: the upstream and downstream SNPs all get called properly, also, the mapping qualities of concerning reads are like 29-60, with no Repeat alignments

[lh3]

All reads are mapped with mapping quality 0. You cannot call a SNP from that. As to 062, the mapping quality is also too low to justify it as a SNP. Also note that the two sites may be due to wrong alignments. Samtools largely does the right thing.

[orionzhou]

Thanks for your reply, however, the mapping qualities of the underlying reads are not 0 (which I manually checked, between 29-60, never below 20).

Also, taking a look at the (RMS) mapping qualities of surrounding bases reveals that a sudden change occurs at these 3 locations:

chr1	6269001	T	T	57	0	43	10	......,,,,	7BH@H:HGHH
chr1	6269002	T	T	57	0	43	10	.$.....,,,,	7@HBH8HGFH
chr1	6269003	G	G	54	0	41	9	.$....,,,,	BH?H9HGFH
chr1	6269004	A	A	51	0	42	8	....,,,,	H9D.HDEH
chr1	6269005	C	C	51	0	42	8	....,,,,	EAHAHEHD
chr1	6269006	T	T	51	0	42	8	....,,,,	IAHADAEH
chr1	6269007	G	G	51	0	42	8	....,,,,	H?H>DEFH
chr1	6269008	T	T	51	0	42	8	....,,,,	H;H>HBGH
chr1	6269009	G	G	37	0	42	8	...T,,,,	HAH.HDHH
chr1	6269010	T	T	51	0	42	8	....,,,,	H;H@HHEH
chr1	6269011	A	A	51	0	42	8	....,,,,	E;H;HHHH
chr1	6269012	A	A	54	0	44	9	....,,,,^],	A@H;IHHFE
chr1	6269013	C	C	54	0	44	9	....,,,,,	GAH=HGHF9
chr1	6269014	A	A	54	0	44	9	....,,,,,	FBHAHHHGE
chr1	6269015	A	A	57	0	46	10	....,,,,,^].	H<E@HHHGBE
chr1	6269016	T	T	57	0	46	10	....,,,,,.	HAHEHHHEEI
chr1	6269017	T	T	57	0	46	10	....,,,,,.	H?H?HGEF4G
chr1	6269018	G	G	57	0	46	10	....,,,,,.	HAH>HGDH5H
chr1	6269019	T	T	57	0	46	10	.$...,,,,,.	E<H;HHHF4H
chr1	6269020	A	A	54	0	44	9	...,,,,,.	AE1HHHH4H
chr1	6269021	A	A	51	0	42	9	.$.$.,,,,,.	:BFHHHH,H
chr1	6269022	T	T	45	0	40	7	.,,,,,.	7GDBF,H
chr1	6269023	G	G	51	0	46	8	.,,,,,.^],	2HHCF5GD
chr1	6269024	A	A	51	0	46	8	.,,,,,.,	@FHBB5HF
chr1	6269025	C	C	51	0	46	8	.,,,,,.,	1HH=F5HA
chr1	6269026	T	T	48	0	44	8	.,,,,,.,	;EH9F,HA
chr1	6269027	T	T	51	0	46	9	.,,,,,.,^].	>HHDF*HD>
chr1	6269028	T	T	54	0	48	9	.,,,,,.,.	AEGHG4HF>
chr1	6269029	C	C	54	0	48	9	.$,,,,,.,.	6HGGH4HGH
chr1	6269030	G	G	51	0	46	8	,,,,,.,.	FHHH5H=H
chr1	6269031	C	C	51	0	46	8	,,,,,.,.	HHHH4HFH
chr1	6269032	T	T	51	0	46	8	,,,,,.,.	IHFF/HDH
chr1	6269033	T	T	48	0	44	8	,,,,,.,.	HHHH+H8H
chr1	6269034	C	C	51	0	46	8	,$,,,,.,.	ECEE:HCH
chr1	6269035	A	A	48	0	48	7	,,,,.,.	<EEAHCH
chr1	6269036	A	A	48	0	48	7	,,,,.,.	BEE=HCG
chr1	6269037	A	A	48	0	48	7	,,,,.,.	CEE>GCH
chr1	6269038	A	A	48	0	48	7	,,,,.,.	GHHEGHH
chr1	6269039	A	A	48	0	48	7	,,,,.,.	HGH9HHH
chr1	6269040	A	A	48	0	48	7	,,,,.,.	HHHAHHF
chr1	6269041	G	G	48	0	48	7	,,,,.,.	HCEAH@H
chr1	6269042	A	A	48	0	48	7	,,,,.,.	GBEEHEE
chr1	6269043	A	A	48	0	48	8	,,,,.,.^],	GACEHEG%
chr1	6269044	G	G	48	0	48	8	,$,,,.,.,	EECBH:G%
chr1	6269045	A	A	45	0	51	7	,,,.,.,	CCEHEH%
chr1	6269046	A	A	45	0	51	7	,,,.,.,	EEAHEH%
chr1	6269047	A	A	45	0	51	7	,,,.,.,	CE:HEH%
chr1	6269048	A	A	48	0	53	7	,,,.,.,	CE@HEG?
chr1	6269049	A	A	51	0	50	8	,,,.,.,^>,	HH>HDG57
chr1	6269050	A	A	54	0	49	9	,,,.,.,,^>,	FH>HHGD78
chr1	6269051	A	A	54	0	49	9	,,,.,.,,,	HH@HHGA88
chr1	6269052	A	A	57	0	47	10	,,,.,.,,,^>,	GH9HHG999;
chr1	6269053	A	A	57	0	47	10	,,,.,.,,,,	DHEGHGA;;;
chr1	6269054	A	A	57	0	47	11	,,,.,.,,,,^0,	@HFHHG8>>>%
chr1	6269055	T	T	57	0	47	11	,,,.,.,,,,,	FFBH7G/9EA%
chr1	6269056	C	C	60	0	48	12	,$,,.,.,,,,,^].	DFEFHH@BHG%E
chr1	6269057	A	A	57	0	50	11	,,.,.,,,,,.	FCEHH8@HH%H
chr1	6269058	A	A	57	0	50	11	,,.,.,,,,,.	F@HHH5BHH%H
chr1	6269059	A	A	57	0	50	11	,,.,.,,,,,.	F:FHH8@HH%H
chr1	[COLOR="Red"]6269060[/COLOR]	T	N	0	0	[COLOR="Red"]0[/COLOR]	11	,,.,.,,,,,.	&&&&&),&&%&
chr1	[COLOR="Red"]6269061[/COLOR]	T	N	0	0	[COLOR="Red"]0[/COLOR]	13	ccCcCcccccC^]C^]c	!!!!!$(!!%!!!
chr1	[COLOR="Red"]6269062[/COLOR]	C	C	15	0	[COLOR="Red"]15[/COLOR]	13	ttTtTtttttTTt	!!!!!$(!!0!!!
chr1	6269063	G	G	72	0	48	15	,,.,.,,,,,..,^>,^>,	H:GBH@:HF>H.<?E
chr1	6269064	T	T	72	0	48	15	,,.,.,,,,,..,,,	H;H=I:<IDAHD4;?
chr1	6269065	C	C	72	0	48	15	,,.,.,,,,,..,,,	HEHFH/<HGBHB7?C
chr1	6269066	C	C	72	0	48	15	,,.,.,,,,,..,,,	HADFH>@HH8HF>?C
chr1	6269067	T	T	75	0	49	17	,,.,.,,,,,..,,,^].^]c	E1H>H54HH=H@/?:@+
chr1	6269068	T	T	75	0	49	17	,,.,.,,,,,..,,,.,	HEHHHA<HEFHD9?H?,
chr1	6269069	G	G	78	0	50	17	,,.,.,,,,,..,,,.,	HEHGH:?HHDHA?@H8B
chr1	6269070	T	T	66	0	50	18	,,.,.,,,,,..,,,A,^].	HA@GH@?HH>GB;=H/FE
chr1	6269071	A	A	84	0	51	19	,,.,.,,,,,..,,,.,.^],	HCFHH:<HE?F867G;BH7
chr1	6269072	C	C	84	0	51	19	,,.,.,,,,,..,,,.,.,	HEEGE<BHH8C1@=H@FHD
chr1	6269073	A	A	84	0	51	19	,,.,.,,,,,..,,,.,.,	H:HHA:@HG<CB?>HD@H?
chr1	6269074	T	T	84	0	51	19	,,.,.,,,,,..,,,.,.,	EE=GH?@HFHHF?@>DCHA
chr1	6269075	G	G	84	0	51	19	,,.,.,,,,,..,,,.,.,	HEHHGDAHEFH4?<H;@H7
chr1	6269076	T	T	84	0	51	19	,$,.,.,,,,,..,,,.,.,	EEHHF:@HHBGH?@H@FGC
chr1	6269077	G	G	81	0	51	18	,.,.,,,,,..,,,.,.,	EAHH/AGFHH=?@HCBH?
chr1	6269078	G	G	81	0	51	18	,.,.,,,,,..,,,.,.,	E=HG?<HHHGH@BHDDGB
chr1	6269079	A	A	65	0	51	18	,.,.,,,,g..,,,.,.,	BHHH>AHH0HH@AH?DHB
chr1	6269080	A	A	81	0	51	18	,.,.,,,,,..,,,.,.,	EBHH:@FHDHF?BHCBED
chr1	6269081	G	G	81	0	51	18	,.,.,,,,,..,,,.,.,	E==HA@HBAH8@?HEDHB
chr1	6269082	C	C	84	0	52	19	,.,.,,,,,..,,,.,.,^],	B=HG<BGGACBB>HGCHHA
chr1	6269083	T	T	87	0	52	20	,.,.,,,,,..,,,.,.,,^].	ABHEB?HH?BFA@HC8HHDE
chr1	6269084	G	G	87	0	52	20	,.,.,,,,,..,,,.,.,,.	B8GGB@HH?HGA@HC8FHEH
chr1	6269085	A	A	87	0	52	20	,.,.,,,,,..,,,.,.,,.	AGGDG=HI;B=?=FGAHHAH
chr1	6269086	C	C	87	0	52	20	,.,.,,,,,..,,,.,.,,.	AFGGF?HH?AB=<H;@GHDH
chr1	6269087	T	T	84	0	52	20	,.,.,,,,,..,,,.,.,,.	@DHH8:HH8?G+@HGAGHCH
chr1	6269088	T	T	87	0	52	20	,.,.,,,,,..,,,.,.,,.	4CFHC<HHAAA9BHF?HHGH
chr1	6269089	G	G	87	0	52	20	,.,.,,,,,..,,,.,.,,.	9F=@HBHH=?8BBHB>HHEI
chr1	6269090	A	A	87	0	52	20	,.,.,,,,,..,,,.,.,,.	=HEGHBHH><E>BHG=HHGH
chr1	6269091	A	A	87	0	52	20	,.,.,,,,,..,,,.,.,,.	0FHHFBHH:E?BBH7=HCHH
chr1	6269092	A	A	84	0	52	20	,.,.,,,,,..,,,.,.,,.	;GGHA?HH=H?>?IE+HCDH
chr1	6269093	C	C	87	0	52	20	,.,.,,,,,..,,,.,.,,.	BDFBGCHH=HFBAHE>FEHE
chr1	6269094	T	T	84	0	52	20	,.,.,,,,,..,,,.,.,,.	*HDF;>HA@HF?AEG/HDCH
chr1	6269095	T	T	87	0	52	20	,.,.,,,,,..,,,.,.,,.	ADHHH=HHFFG@AHB>HCGH
chr1	6269096	C	C	87	0	52	20	,.,.,,,,,..,,,.,.,,.	;HHBHBHHHHFABHG>HHHH
chr1	6269097	A	A	87	0	52	20	,.,.,,,,,..,,,.,.,,.	EGH?HBHHHGE@AH<FHHGH
chr1	6269098	T	T	90	0	53	21	,.,.,,,,,..,,,.,.,,.^].	EFGFI?HHEDDBBHEFHHHHE
chr1	6269099	G	G	90	0	53	21	,.,.,,,,,..,,,.,.,,..	EFHEHBHFFGC@@HAEHGHHH
chr1	6269100	T	T	93	0	53	22	,.,.,,,,,..,,,.,.,,..^].	EDHFH?HH@BFABH4GHHEHHE

[lh3]

I don't know why unless i see the reads. But anyway i would not trust SNPs on reads containing 5 mismatches! They are likely to be false alignments.

EDIT: I see the reason. Your alignment can be:

aaattcgtcctt
aaatCTgtcctt

or

aaattc-gtcctt
aaat-ctgtcctt

Would you prefer two mismatches or two indels? We may prefer the former, but in the low-complexity region you have shown, the chance of having two indels historically is not slim. The latest samtools only calls a SNP when it is unlikely to be explained by any combinations of indels.

[orionzhou]

Thanks for your reply

Really appreciate your comments, and sorry for the confusion - this is not human resequencing data, the resequenced plant species is quite distant to the reference plant - that's why we see so many mismatches (if not sequencing errors).

Click me for the BAM file of this 100bp region, after running:

samtools view -b 07_bam/HM020.bam chr1:6,269,001-6,269,100 > test.bam

It's not possible to attach the reference sequence, but surprisingly, running "samtools pileup" without a reference sequence seems to work just fine:

 samtools pileup -c test.bam

chr1    6269056 N       C       33      235     48      12      c$ccCcCccccc^]C DFEFHH@BHG%H
chr1    6269057 N       A       30      230     50      11      aaAaAaaaaaA     FEEHH8@HH%H
chr1    6269058 N       A       30      232     50      11      aaAaAaaaaaA     FEHHH5BHH%H
chr1    6269059 N       A       30      225     50      11      aaAaAaaaaaA     F:FHH8@HH%H
chr1    [COLOR="Red"]6269060[/COLOR] N       [COLOR="Red"]T[/COLOR]       30      217     50      11      ttTtTtttttT     F9HBH48IC%H
chr1    [COLOR="Red"]6269061[/COLOR] N       [COLOR="Red"]C[/COLOR]       36      255     52      13      ccCcCcccccC^]C^]c       FEIHH<=HH%HD@
chr1    [COLOR="Red"]6269062[/COLOR] N       [COLOR="Red"]T[/COLOR]       39      252     50      13      ttTtTtttttTTt   EFHDH;4GF@HD9
chr1    6269063 N       G       45      248     48      15      ggGgGgggggGGg^>g^>g     H:GBH@:HF>H.>?E
chr1    6269064 N       T       45      250     48      15      ttTtTtttttTTttt H;H=I:<IDAHD4;?
chr1    6269065 N       C       45      255     48      15      ccCcCcccccCCccc HEHFH/<HGBHB7?C

whereas running it with a reference makes me unhappy:

samtools pileup -f genome.fa -c test.bam

chr1    6269056 C       C       60      0       48      12      ,$,,.,.,,,,,^]. DFEFHH@BHG%E
chr1    6269057 A       A       57      0       50      11      ,,.,.,,,,,.     FCEHH8@HH%H
chr1    6269058 A       A       57      0       50      11      ,,.,.,,,,,.     F@HHH5BHH%H
chr1    6269059 A       A       57      0       50      11      ,,.,.,,,,,.     F:FHH8@HH%H
chr1    [COLOR="Red"]6269060[/COLOR] T       N       0       0       0       11      ,,.,.,,,,,.     &&&&&),&&%&
chr1    [COLOR="Red"]6269061[/COLOR] T       N       0       0       0       13      ccCcCcccccC^]C^]c       !!!!!$(!!%!!!
chr1    [COLOR="Red"]6269062[/COLOR] C       C       15      0       15      13      ttTtTtttttTTt   !!!!!$(!!0!!!
chr1    6269063 G       G       72      0       48      15      ,,.,.,,,,,..,^>,^>,     H:GBH@:HF>H.<?E
chr1    6269064 T       T       72      0       48      15      ,,.,.,,,,,..,,, H;H=I:<IDAHD4;?
chr1    6269065 C       C       72      0       48      15      ,,.,.,,,,,..,,, HEHFH/<HGBHB7?C

[orionzhou]

Thanks a lot for your help!

ok I see - that makes sense, but would that make it really difficult to call 2 consecutive SNPs - since they could always be explained by 2 indels:

AG
GT

or

AG-
-GT

but anyway, that only makes it more stringent - which is good.

[lh3]

Actually a more proper way is to say there is a variant at the locus but what is actually happening (two SNPs or two indels) is unknown. Unfortunately, no callers give such calls so far as I know.

[orionzhou]

Exactly. Thanks a lot for your time and have a good day!

[bioinfosm]

Why does using a reference generate a different pileup file?

Really appreciate your comments, and sorry for the confusion - this is not human resequencing data, the resequenced plant species is quite distant to the reference plant - that's why we see so many mismatches (if not sequencing errors).

Click me for the BAM file of this 100bp region, after running:

samtools view -b 07_bam/HM020.bam chr1:6,269,001-6,269,100 > test.bam

It's not possible to attach the reference sequence, but surprisingly, running "samtools pileup" without a reference sequence seems to work just fine:

 samtools pileup -c test.bam

chr1    6269056 N       C       33      235     48      12      c$ccCcCccccc^]C DFEFHH@BHG%H
chr1    6269057 N       A       30      230     50      11      aaAaAaaaaaA     FEEHH8@HH%H
chr1    6269058 N       A       30      232     50      11      aaAaAaaaaaA     FEHHH5BHH%H
chr1    6269059 N       A       30      225     50      11      aaAaAaaaaaA     F:FHH8@HH%H
chr1    [COLOR="Red"]6269060[/COLOR] N       [COLOR="Red"]T[/COLOR]       30      217     50      11      ttTtTtttttT     F9HBH48IC%H
chr1    [COLOR="Red"]6269061[/COLOR] N       [COLOR="Red"]C[/COLOR]       36      255     52      13      ccCcCcccccC^]C^]c       FEIHH<=HH%HD@
chr1    [COLOR="Red"]6269062[/COLOR] N       [COLOR="Red"]T[/COLOR]       39      252     50      13      ttTtTtttttTTt   EFHDH;4GF@HD9
chr1    6269063 N       G       45      248     48      15      ggGgGgggggGGg^>g^>g     H:GBH@:HF>H.>?E
chr1    6269064 N       T       45      250     48      15      ttTtTtttttTTttt H;H=I:<IDAHD4;?
chr1    6269065 N       C       45      255     48      15      ccCcCcccccCCccc HEHFH/<HGBHB7?C

whereas running it with a reference makes me unhappy:

samtools pileup -f genome.fa -c test.bam

chr1    6269056 C       C       60      0       48      12      ,$,,.,.,,,,,^]. DFEFHH@BHG%E
chr1    6269057 A       A       57      0       50      11      ,,.,.,,,,,.     FCEHH8@HH%H
chr1    6269058 A       A       57      0       50      11      ,,.,.,,,,,.     F@HHH5BHH%H
chr1    6269059 A       A       57      0       50      11      ,,.,.,,,,,.     F:FHH8@HH%H
chr1    [COLOR="Red"]6269060[/COLOR] T       N       0       0       0       11      ,,.,.,,,,,.     &&&&&),&&%&
chr1    [COLOR="Red"]6269061[/COLOR] T       N       0       0       0       13      ccCcCcccccC^]C^]c       !!!!!$(!!%!!!
chr1    [COLOR="Red"]6269062[/COLOR] C       C       15      0       15      13      ttTtTtttttTTt   !!!!!$(!!0!!!
chr1    6269063 G       G       72      0       48      15      ,,.,.,,,,,..,^>,^>,     H:GBH@:HF>H.<?E
chr1    6269064 T       T       72      0       48      15      ,,.,.,,,,,..,,, H;H=I:<IDAHD4;?
chr1    6269065 C       C       72      0       48      15      ,,.,.,,,,,..,,, HEHFH/<HGBHB7?C

[orionzhou]

From my understanding, when calling variants using a reference, samtools will only make confident calls: if there's another explanation for the variants (in this case, 2 indels rather than 2 consecutive SNPs), the program will not call anything but leave unknown states for these positions.