Skip to content

Allow MNPs to be reported as SNPs #19

Open
jamescasbon opened this Issue Oct 17, 2011 · 4 comments

2 participants

@jamescasbon

I have a good quality GG/TC MNP:

$ freebayes -f /data/reference/ucsc/hg19/hg19.fa -r chr22:42526560..42526563 -F 0.3  reads.bam 2> /dev/null| cut -f1-6 | grep -v "##"
chr22   42526561    .   GG  TC  49314.7

But I don't want MNPs, I want two SNPs. So I'll use the --no-mnps flag, right? Doesn't work:

$ freebayes -f /data/reference/ucsc/hg19/hg19.fa -r chr22:42526560..42526563 -F 0.3 --no-mnps reads.bam 2> /dev/null| cut -f1-6 | grep -v "##"
chr22   42526561    .   GGTGGTGGGGCATCCTCAGG    TCTGTAGGGGCATCCTCAGG,TCTCGGTAGGGGAGCCTCAGG,TCTGTAGGGAGCCTCAGG,TCTGGCAGGGGAGCCTCAGG,TCCGGTAGGGGAGCCTCAGG,TCTGTACGGGGAGCCTCAGG,TCAGGTAGGGGAGCCTCAGG,TCTGGGTAGGGGAGCCTCAGG,TCTAGTAGGGGAGCCTCAGG,TCTGGTTAGGGGAGCCTCAGG,TCCGGTAGGGGAGCCTTAGC,TCTGGGAGGGGAGCCTCAGG,TCTGTTAGGGGAGCCTCAGG,TCTGTAGGGAGCTCAGG,TCTTGTAGGGGAGCCTCAGG,TCTCGGTAGGGGAGCCGTCAGG,TCTCGGTACGGGGAGCCTCAGG,TTCGGTACGGGGAGCCGCTCAGG  14761.8

OK, so how about adding in --no-complex? That doesn't even give the best SNPs:

$ freebayes -f /data/reference/ucsc/hg19/hg19.fa -r chr22:42526560..42526563 -F 0.3 --no-mnps --no-complex reads.bam 2> /dev/null| cut -f1-6 | grep -v "##"
chr22   42526561    .   G   C   62.4447
chr22   42526562    .   G   C   0.0781914

Weird, so what about trying to force it to produce the best alleles?

$ freebayes -f /data/reference/ucsc/hg19/hg19.fa -r chr22:42526560..42526563 -F 0.3 --use-best-n-alleles 1 --no-mnps --no-complex reads.bam 2> /dev/null| cut -f1-6 | grep -v "##"
chr22   42526561    .   G   C   62.4447
chr22   42526562    .   G   C   0.0781914

So how I can get freebayes to report the GG/TC as two individual SNPs?

@ekg
Owner
ekg commented Oct 17, 2011
@jamescasbon

I see. I want to compare this to known snps and each site is annotated as a single SNP:

http://genome.ucsc.edu/cgi-bin/hgc?hgsid=218425869&o=42526560&t=42526561&g=snp132&i=rs28695233
http://genome.ucsc.edu/cgi-bin/hgc?hgsid=218425869&o=42526561&t=42526562&g=snp132&i=rs75276289

So I was hoping --no-mnps would make these single calls so that I could process them downstream without splitting the calls. You should close this ticket - but maybe it's a case for a contrib script?

@ekg
Owner
ekg commented Oct 18, 2011
@ekg
Owner
ekg commented Dec 14, 2011

Please see: https://github.com/ekg/vcflib/blob/master/vcfallelicprimitives.cpp

This is a first attempt to allow the flattening of small haplotypes back into their primitive allelic components. Thus an MNP will yield a set of SNPs. Genotypes are not currently handled, but this extension could be developed.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Something went wrong with that request. Please try again.