This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Extract unphased genotypes from vcf file and convert to 012 matrix

I have a vcf file in which I want to extract the genotypes and convert them into a 0,1, 2 matrix. Is there a way to extract only the unphased genotypes using Plink or VCF tools or maybe with grep?

Thank you in advance

vcf snp
• 6,172 views • link
updated by Ram Ram • written by Nev • Nev

not clear. what is the result of "I've already extracted all of the genotypes,"

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Sorry, I rephrased. Thanks for pointing that out.

• 0 views • link
written by Nev • Nev

2 answers

Using my tool bioalcidae https://github.com/lindenb/jvarkit/wiki/BioAlcidae

while(iter.hasNext())
    {
    var ctx = iter.next();
    out.print(ctx.contig);
    out.print("\t");
    out.print(ctx.start);
    for(var I = 0;i< ctx.getNSamples();++i)
        {
        var g = ctx.getGenotype(i);
        out.print("\t");
        if(g.isCalled() && !g.isPhased())
            {
            if(g.isHomVar())
                {
                out.print("2");
                }
            else if(g.isHomRef())
                {
                out.print("0");
                }
            else if(g.isHet() && !g.isHetNonRef())
                {
                out.print("1");
                }
            else
                {
                out.print("9");
                }
            }
        else
            {
            out.print("9");
            }
        }
    out.println();
    }

Usage:

java -jar dist-1.133/bioalcidae.jar  -F VCF -f filter.js your.vcf
• 0 views • link
updated by Ram Ram • written by Pierre Lindenbaum Pierre Lindenbaum

Ok, great. Thank you Pierre. Is there citation info for your tool? I will give this a try now!

• 0 views • link
written by Nev • Nev

Is there any simple way to modify the script to print sample names at the top of the output?

• 0 views • link
written by tantrev • tantrev

there is no sample name at the top of the output of this script.

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

It's cool, I figured out how to print the names with the example VCF script on the bioalcidae homepage. :) On a different note, though, I think there's a typo in the script above - I believe the "g.isHomRef()" clause should print a 0, not a 1, and vice-versa with the "g.isHet() && !g.isHetNonRef()" statement.

• 0 views • link
written by tantrev • tantrev

I believe the "g.isHomRef()" clause should print a 0, not a 1, and vice-versa with the "g.isHet() && !g.isHetNonRef()"

ah yes, you're right

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

:) Is there any way to make the script above allele-specific? I'd love to see the exact same output, just with an individual comparison to each potential homozygous A,T,C, & G allele - instead of just the reference allele.

• 0 views • link
written by tantrev • tantrev

how can I get the filter.js file?

• 0 views • link
written by lokman20102003 • lokman20102003

Log in to answer this question.

More posts like this