How to find heterozygote by ssahaSNP

[shuang]

I'm trying to use ssahaSNP for SNP/INDEL discovery. However, I only received homozygote cases. How do I find heterozygote via this program? Below are the 4 commands I used based the instruction from http://computing.bio.cam.ac.uk/local/doc/ssaha2.pdf


ssahaSNP allinput.fastq Sorbil.fasta -tags 1 > all_ssahaSNP.out

egrep ALIGN all_ssahaSNP.out | awk '{print $1,$2,$3,$4,$5,$6,$7,$8,$9,$10,$11,$12}' > all_alignment.dat

egrep ssaha:SNP all_ssahaSNP.out > all_SNP.dat

parse_SNP all_SNP.dat all_alignment.dat > all_parseSNP.out

[sdvie]

quote from the ssaha manual:

ssaha2SNP is a polymorphism detection tool. It detects homozygous SNPs
and indels by aligning shotgun reads to the finished genome sequence.
From the best alignment, SNP candidates are screened, taking into account
the quality value of the bases with variation as well as the quality
values in the neighbouring bases, using neighbourhood quality standard
(NQS).

Looks like this executable is not meant for the detection of heterozygous SNPs.

Try to transform your alignment output file to another format, such as bam, to be able to use a different SNP calling program that detects both homo- and heterozygous SNPs.

Cheers!

[shuang]

which program would you recommend to detect both homo and hetero SNP/INDEL?

[sdvie]

Once you have your alignments in bam format, there is a variety of variant callers that output the so-called "variant call format" vcf that contains the information of zygosity that you are looking for:

Try:

- samtools + bcftools
- samtools + varscan
- GATK

cheers!

PS.: Keep reading this excellent forum to find out more!

[shuang]

My purpose is
1. find all SNP/INDEL
2. If a (Sanger) read which covers a SNP/INDEL position but it's not detected as SNP/INDEL, mark the position of that sequence as wild-type.


I tried samtools + bcftools, however, I couldn't accomplish the 2nd purpose. Let's why I used ssahaSNP. ssahaSNP provides leftmost and rightmost position of each alignment on the reference genome. By comparing a SNP position with those 2 positions, I could tell if a read covered that SNP position. Then, if that read was not on the list of SNP, it's wild-type.
However, I wasn't aware that ssahaSNP didn't detect hetero.

Could any of those 3 methods you suggest accomplish the 2nd purpose? Could you explain more about how? Thanks!

- samtools + bcftools
- samtools + varscan
- GATK