what does "hetero SNP/INDEL" look like in BWA-SW+Samtools

shuang · 10-11-2011, 07:20 AM

I'm using BWA-SW+Samtools to analyze SNP/INDEL of Sanger sequencing data. My samples may contain heterozygotes.. but not sure yet.

Can anyone give me an example of what "hetero SNP/INDEL" looks like in an alignment output from BWA-SW and a pileup output from Samtools?

Heisman · 10-11-2011, 08:04 AM

Not sure if alignment will tell you that, but if you use SAMtools pileup the output in column 4 will tell you if it's a heterozygous or homozygous call (following this code:http://biocorp.ca/IUB.php). If you use mpileup than if AC1=1 it's heterozygous and if AC1=2 it's homozygous (for a single individual), or you can look at the GT filed in the last column (0/1 or 0|1 is heterozygous and 1/1 or 1|1 is homozygous).

shuang · 10-11-2011, 09:26 AM

Thanks for the great information.

However, how to interpretate when AC1=2, but GT=0/1? Below is an example from my output.

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 5128_RF1-1_RN1.bam
chr_8 50762920 . T C 4.77 . DP=1;AF1=1;AC1=2;DP4=0,0,0,1;MQ=60;FQ=-30 GT:PL:GQ 0/1:33,3,0:3

Heisman · 10-11-2011, 09:30 AM

That's a weird situation because you have 1 read, so I'm pretty sure the normal rules don't apply. Do you ever see that with 2 or more reads?

shuang · 10-11-2011, 09:38 AM

I only pile upped 1 Sanger read at a time. Below is the actual full output for this sequence. A similar situation happens on the last INDEL, too. I do notice such confusion usually happens on a SNP with low QUAL.

How do I retrieve correct SNP/INDEL from such outputs? Or any other solution?

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 5128.bam
chr_8 50762920 . T C 4.77 . DP=1;AF1=1;AC1=2;DP4=0,0,0,1;MQ=60;FQ=-30 GT:PL:GQ 0/1:33,3,0:3
chr_8 50763047 . T C 26 . DP=1;AF1=1;AC1=2;DP4=0,0,0,1;MQ=60;FQ=-30 GT:PL:GQ 1/1:56,3,0:6
chr_8 50763143 . T A 26 . DP=1;AF1=1;AC1=2;DP4=0,0,0,1;MQ=60;FQ=-30 GT:PL:GQ 1/1:56,3,0:6
chr_8 50763248 . G A 26 . DP=1;AF1=1;AC1=2;DP4=0,0,0,1;MQ=60;FQ=-30 GT:PL:GQ 1/1:56,3,0:6
chr_8 50763265 . G T 26 . DP=1;AF1=1;AC1=2;DP4=0,0,0,1;MQ=60;FQ=-30 GT:PL:GQ 1/1:56,3,0:6
chr_8 50763559 . c cG 3.81 . INDEL;DP=1;AF1=1;AC1=2;DP4=0,0,0,1;MQ=60;FQ=-37.5 GT:PL:GQ 0/1:39,3,0:4

Heisman · 10-11-2011, 09:44 AM

Oh, I see. I haven't looked at Sanger sequencing in this way before so I'm not sure if relying on the quality scores is a good way to determine heterozygous from homozygous.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Samtools "is recognized as '*'" "truncated file" error	axiom7	Bioinformatics	3	11-26-2014 03:53 AM
suitable tools for hetero SNP/INDEL discovery?	shuang	Bioinformatics	0	10-04-2011 09:18 AM
The position file formats ".clocs" and "_pos.txt"? Ist there any difference?	elgor	Illumina/Solexa	0	06-27-2011 08:55 AM
"Systems biology and administration" & "Genome generation: no engineering allowed"	seb567	Bioinformatics	0	05-25-2010 01:19 PM
SEQanswers second "publication": "How to map billions of short reads onto genomes"	ECO	Literature Watch	0	06-30-2009 12:49 AM

10-11-2011, 07:20 AM	#1
shuang Senior Member Location: IL Join Date: Jul 2011 Posts: 100	what does "hetero SNP/INDEL" look like in BWA-SW+Samtools I'm using BWA-SW+Samtools to analyze SNP/INDEL of Sanger sequencing data. My samples may contain heterozygotes.. but not sure yet. Can anyone give me an example of what "hetero SNP/INDEL" looks like in an alignment output from BWA-SW and a pileup output from Samtools?

10-11-2011, 08:04 AM	#2
Heisman Senior Member Location: St. Louis Join Date: Dec 2010 Posts: 535	Not sure if alignment will tell you that, but if you use SAMtools pileup the output in column 4 will tell you if it's a heterozygous or homozygous call (following this code:http://biocorp.ca/IUB.php). If you use mpileup than if AC1=1 it's heterozygous and if AC1=2 it's homozygous (for a single individual), or you can look at the GT filed in the last column (0/1 or 0\|1 is heterozygous and 1/1 or 1\|1 is homozygous).

10-11-2011, 09:26 AM	#3
shuang Senior Member Location: IL Join Date: Jul 2011 Posts: 100	Thanks for the great information. However, how to interpretate when AC1=2, but GT=0/1? Below is an example from my output. #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 5128_RF1-1_RN1.bam chr_8 50762920 . T C 4.77 . DP=1;AF1=1;AC1=2;DP4=0,0,0,1;MQ=60;FQ=-30 GT:PL:GQ 0/1:33,3,0:3

10-11-2011, 09:30 AM	#4
Heisman Senior Member Location: St. Louis Join Date: Dec 2010 Posts: 535	That's a weird situation because you have 1 read, so I'm pretty sure the normal rules don't apply. Do you ever see that with 2 or more reads?

10-11-2011, 09:38 AM	#5
shuang Senior Member Location: IL Join Date: Jul 2011 Posts: 100	I only pile upped 1 Sanger read at a time. Below is the actual full output for this sequence. A similar situation happens on the last INDEL, too. I do notice such confusion usually happens on a SNP with low QUAL. How do I retrieve correct SNP/INDEL from such outputs? Or any other solution? #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 5128.bam chr_8 50762920 . T C 4.77 . DP=1;AF1=1;AC1=2;DP4=0,0,0,1;MQ=60;FQ=-30 GT:PL:GQ 0/1:33,3,0:3 chr_8 50763047 . T C 26 . DP=1;AF1=1;AC1=2;DP4=0,0,0,1;MQ=60;FQ=-30 GT:PL:GQ 1/1:56,3,0:6 chr_8 50763143 . T A 26 . DP=1;AF1=1;AC1=2;DP4=0,0,0,1;MQ=60;FQ=-30 GT:PL:GQ 1/1:56,3,0:6 chr_8 50763248 . G A 26 . DP=1;AF1=1;AC1=2;DP4=0,0,0,1;MQ=60;FQ=-30 GT:PL:GQ 1/1:56,3,0:6 chr_8 50763265 . G T 26 . DP=1;AF1=1;AC1=2;DP4=0,0,0,1;MQ=60;FQ=-30 GT:PL:GQ 1/1:56,3,0:6 chr_8 50763559 . c cG 3.81 . INDEL;DP=1;AF1=1;AC1=2;DP4=0,0,0,1;MQ=60;FQ=-37.5 GT:PL:GQ 0/1:39,3,0:4

10-11-2011, 09:44 AM	#6
Heisman Senior Member Location: St. Louis Join Date: Dec 2010 Posts: 535	Oh, I see. I haven't looked at Sanger sequencing in this way before so I'm not sure if relying on the quality scores is a good way to determine heterozygous from homozygous.