Compare samples in one multisample VCF file

Question

[Alice]

Hello biostars!

I was wondering, how to compare a few samples in my multisample vcf file?

Basically, I would like to get a venn diagram, but vcf-compare is working only on multiple vcf files, bcftools stats and plot-vcfstats do not help either (they do not compare positions across the samples).

I can just split my vcf on three separate files and compare them using the approach described above, but I believe there should be more simple solution.

Answer 1

[Pierre Lindenbaum]

Using Bioalcidae and the following script?

var counts={};
while(iter.hasNext()) {
var ctx  = iter.next();
for(var i=0;i< ctx.getNSamples();++i)
    {
    var gi = ctx.getGenotype(i);
    for(var j=i+1;j< ctx.getNSamples();++j)
        {
        var gj = ctx.getGenotype(j);
        if(  gi.sameGenotype(gj) ) {
            var key = gi.getSampleName()+"-"+gj.getSampleName();
            var c= counts[key];
            if(c==null) c=0;
            ++c;
            counts[key]=c;
            }
        }
    }
}

for(var i in counts) {
    out.println("["+i+"]\t"+counts[i]);
    }

run:

gunzip -c input.vcf.gz| java -jar dist-2.0.1/bioalcidae.jar -F vcf -f script.js | column -t

[S1-S10]  35
[S1-S2]   16
[S1-S3]   20
[S10-S2]  16
[S10-S3]  20
[S2-S3]   19
[S1-S4]   16
[S10-S4]  16
[S3-S4]   19
[S2-S4]   15

[Alice]

Cool, I did not know about this package! And thanks for the script! I will try it asap.

I also need to get a list of snips themselves. As I believe, I can easily catch them using this file after BioAlcidae.

[Pierre Lindenbaum]

It's not a R package. It's a java program.

You should test it before validating my answer :-)

[Alice]

Ok, I see.

No, I mean I will parse this file into R.