Question: Extract unphased genotypes from vcf file and convert to 012 matrix
1
gravatar for Nev
5.2 years ago by
Nev20
United States
Nev20 wrote:

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.

 

snp vcf • 3.4k views
ADD COMMENTlink modified 5.2 years ago by Pierre Lindenbaum131k • written 5.2 years ago by Nev20

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

ADD REPLYlink written 5.2 years ago by Pierre Lindenbaum131k

Sorry, I rephrased. Thanks for pointing that out.

ADD REPLYlink written 5.2 years ago by Nev20
3
gravatar for Pierre Lindenbaum
5.2 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum131k wrote:

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

 

ADD COMMENTlink modified 3.5 years ago • written 5.2 years ago by Pierre Lindenbaum131k

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

ADD REPLYlink written 5.2 years ago by Nev20

https://github.com/lindenb/jvarkit/wiki/Citing

ADD REPLYlink written 5.2 years ago by Pierre Lindenbaum131k

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

ADD REPLYlink written 3.5 years ago by tantrev30

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

ADD REPLYlink written 3.5 years ago by Pierre Lindenbaum131k

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.

ADD REPLYlink written 3.5 years ago by tantrev30

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

ADD REPLYlink written 3.5 years ago by Pierre Lindenbaum131k

:) 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.

ADD REPLYlink modified 3.5 years ago • written 3.5 years ago by tantrev30
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1160 users visited in the last hour