Visualize Variant Density using Python

Question

[kristina.mahan]

I found these instructions on visualizing variant density using python: http://alimanfoo.github.io/2016/06/10/scikit-allel-tour.html

Step 1 is to open HDF5 callset:

callset_fn = 'data/2016-06-10/ag1000g.phase1.ar3.h5'
callset = h5py.File(callset_fn, mode='r')
callset

Step 2 is to pick a chromosome to work with:

chrom = '3L'

In my file contigs in the #CHROM column are labelled: unitig_2993_pilon

When I pick a contig to work with:

chrom = 'unitig_2993_pilon'

and then proceed:

variants = allel.VariantChunkedTable(callset[chrom]['variants'], names=['POS', 'REF', 'ALT', 'DP', 'MQ', 'QD', 'num_alleles'], index='POS')
variants

I get this error: KeyError: "Unable to open object (object 'unitig_2993_pilon' doesn't exist)"

What am I doing wrong? Do I have to pick a chromosome to work with?

[WouterDeCoster]

Could you take a look at callset.keys()?

[kristina.mahan]

callset.keys()

[u'calldata', u'samples', u'variants']

[WouterDeCoster]

So there is no unitig_2993_pilon in there. Maybe it's under callset['variants'] but that's just a guess since I don't know the format. Either you should search for how the dictionary is formatted, or just explore which keys are present.