This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Taking genotypes out of a vcf file

Hi, I would like to take a vcf file and output a tab-delimited file with a line per individual-site that includes the individual, site, GT, alternate and reference alleles, DP, and each PL. It doesn't matter what order sites end up being in.

Any suggestions on how to do this most efficiently and/or any links to code that does something like this?

Example input:

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  samp1 samp2
chr1   100  .       C       T       3106.72 SnpCluster      .       GT:AD:DP:GQ:PL  0/0:1,0:1:3:0,3,42      0/0:3,0:3:9:0,9,132
chr1   120  .       C       G       3106.72 SnpCluster      .       GT:AD:DP:GQ:PL 0/1:3,1:4:30:30,0,123   1/1:0,1:1:3:45,3,0

Example output:

samp1    chr1    100    0/0   C   T   1    0   3   42
samp2    chr1    100   0/1    C   T   4   30   0   123
samp1    chr1    120    0/0   C   G   3    0    9    132
samp2    chr1    120    1/1   C   G   1    45   3    0
vcf
• 4,485 views • link
updated by Ram Ram • written by Jautis Jautis

If you're into python, pyvcf makes this sort of task quite easy.

• 0 views • link
written by ethan.kaufman ethan.kaufman

2 answers

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

file format.js

while(iter.hasNext())
    {
    var ctx = iter.next();
    for(var I = 0;i< ctx.getNSamples();++i)
        {
        var g = ctx.getGenotype(i);
        out.print(g.getSampleName());
        out.print("\t");
        out.print(ctx.contig);
        out.print("\t");
        out.print(ctx.start);
        out.print("\t");
        if(g.isCalled())
            {
            out.print(ctx.getAlleleIndex(g.getAllele(0))+"/"+ctx.getAlleleIndex(g.getAllele(1)));
            out.print("\t");
            out.print(g.getAllele(0).getDisplayString());
            out.print("\t");
            out.print(g.getAllele(1).getDisplayString());
            }
        else
            {
            out.print(".\t.\t.");
            }
        out.print("\t");
        out.print(g.getDP());
        out.print("\t");
        out.print(g.getGQ());
        out.print("\t");
        out.print(g.getAnyAttribute("NR"));
        out.println();
        }
    }

Invoke:

curl -sL "https://raw.githubusercontent.com/KBoehme/VTCTesting/master/src/test/resources/test_data/UnionTests/Test4/input2.vcf" | java -jar dist-1.133/bioalcidae.jar -F VCF -f format.js

Result:

NA00001    chr1    1451489    0/0    T    T    -1    35    18
NA00002    chr1    1451489    1/0    C    T    -1    100    24
NA00001    chr1    1451490    0/0    G    G    -1    35    18
NA00002    chr1    1451490    0/1    G    A    -1    100    24
NA00001    chr1    12939904    0/0    A    A    -1    45    53
NA00002    chr1    12939904    0/1    A    C    -1    100    79
NA00001    chr1    16902943    0/0    T    T    -1    26    93
NA00002    chr1    16902943    0/1    T    C    -1    84    98
NA00001    chr1    16902982    0/0    G    G    -1    52    108
NA00002    chr1    16902982    0/1    G    T    -1    75    115
(...)
• 0 views • link
updated by Ram Ram • written by Pierre Lindenbaum Pierre Lindenbaum

Thank you very much! That seems perfect. However, I'm having difficulties getting a working copy of bioalcidae.jar. When I download the .jar file, I get an error message saying it is invalid or corrupt. Is there anything else I need to do?

• 0 views • link
updated by Ram Ram • written by Jautis Jautis

(?!) there is no jar fileto download, you need to install and compile the program ;

git clone "https://github.com/lindenb/jvarkit.git" && cd jvarkit &&  make bioalcidae && java -jar dist-1.133/bioalcidae.jar -h
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Oh, thank you. I feel like an idiot now

• 0 views • link
written by Jautis Jautis

I'm partial to the vawk tool. It's simple, but gets the job done quickly.

• 0 views • link
written by Matt Shirley Matt Shirley

Log in to answer this question.

More posts like this