Question: Is it possible to use scikit-allel to generate per sample genotype counts per gene (from gff)?

I've loaded a gff file:

geneset = allel.FeatureTable.from_gff3('~/geneset.gtf', attributes = ["gene_id","gene_name"])

Converted to pandas table following tutorial:

def geneset_to_pandas(geneset):
"""Life is a bit easier when a geneset is a pandas DataFrame."""
items = []
for n in geneset.dtype.names:
v = geneset[n]
# convert bytes columns to unicode (which pandas then converts to object)
if v.dtype.kind == 'S':
v = v.astype('U')
items.append((n, v))
return pandas.DataFrame.from_dict(dict(items))

geneset = geneset_to_pandas(geneset)

Loaded vcf file and converted to genotype array:

callset = allel.read_vcf(called.vcf.gz)
gt = allel.GenotypeArray(callset['calldata/GT'])

Started iterating through geneset to create variables for use with query (of vcf file), but am having trouble aggregating the genotype counts by the positions in the geneset:

for index,row in geneset_Homo_sapiens_GRCh37_75.iterrows():
end = row['end']
start = row['start']
chrom = row['seqid']
gene_name = row['gene_name']
for i,(x,y) in enumerate(zip(callset['variants/CHROM'],callset['variants/POS'])):
if (var_chrom == chrom and var_pos >= start and var_pos <= end):
print var_chrom, var_pos, gene_name, gt[i]

Answer: Here's a [gist with a worked example using human data][1]. A few extracts...

To load data from a GFF3 file into a pandas DataFrame you can now use the scikit-allel gff3_to_dataframe() function, e.g.:

geneset = allel.gff3_to_dataframe('GRCh38_latest_genomic.gff.gz', attributes=['Name'])
# select only gene records
genes = geneset[geneset['type'] == 'gene']
# select only genes on Chromosome 22
genes_chr22 = genes[genes.seqid == 'NC_000022.11']

To extract genotypes for a given sample and a given gene, it's necessary to obtain the variant start and stop indices corresponding to the first and last variants within the gene. If the data are for a single chromosome, this can be done using a SortedIndex on the variant positions.

If you have previously parsed the VCF out into Zarr format grouped by chromosome as described [here](http://alimanfoo.github.io/2017/06/14/read-vcf.html), assuming you are working with a single chromosome (e.g., Chromosome 22), this will set up a SortedIndex for all variants in the chromosome:

zarr_path = 'ALL.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.zarr'
import zarr
callset = zarr.open(zarr_path, mode='r')
chrom = '22'
pos = allel.SortedIndex(callset[chrom]['variants/POS'])

This will load genotypes for a given sample over all variants in the chromosome:

# pick an arbitrary sample to work with
sample_idx = 42 # N.B., this is the 43rd sample, zero-based indexing
# load genotypes for the sample
gv = allel.GenotypeVector(callset[chrom]['calldata/GT'][:, sample_idx])

This will then iterate over genes and compute genotype counts for each gene:

# setup some arrays to hold per-gene genotype counts for our sample of interest
import numpy as np
n_genes = len(genes_chr22)
n_hom_ref = np.zeros(n_genes, dtype=int)
n_het = np.zeros(n_genes, dtype=int)
n_hom_alt = np.zeros(n_genes, dtype=int)
n_variants = np.zeros(n_genes, dtype=int)

# iterate over genes
for i, (_, gene) in enumerate(genes_chr22.iterrows()):
try:
# locate data for this gene - this maps genomic coordinates onto variant start and stop indices
loc_gene = pos.locate_range(gene.start, gene.end)
except KeyError:
# no data for the gene, leave counts as zero
pass
else:
# extract genotypes for the gene
gv_gene = gv[loc_gene]
# compute genotype counts
n_hom_ref[i] = gv_gene.count_hom_ref()
n_het[i] = gv_gene.count_het()
n_hom_alt[i] = gv_gene.count_hom_alt()
# also store number of variants in the gene
n_variants[i] = loc_gene.stop - loc_gene.start

[1]: http://nbviewer.jupyter.org/gist/alimanfoo/6e6b7854a735907e2837f306adf8680b