Applying genotype calls to multi-way alignment incidence matrix

    :param alnfile: alignment incidence file (h5),
    :param gtypefile: genotype calls by GBRS (tsv),
    :param grpfile: gene ID to isoform ID mapping info (tsv)
    :return: genotyped version of alignment incidence file (h5)

def stencil(**kwargs):
    """
    Applying genotype calls to multi-way alignment incidence matrix

    :param alnfile: alignment incidence file (h5),
    :param gtypefile: genotype calls by GBRS (tsv),
    :param grpfile: gene ID to isoform ID mapping info (tsv)
    :return: genotyped version of alignment incidence file (h5)
    """
    alnfile = kwargs.get('alnfile')
    gtypefile = kwargs.get('gtypefile')
    grpfile = kwargs.get('grpfile')
    if grpfile is None:
        grpfile2chk = os.path.join(DATA_DIR, 'ref.gene2transcripts.tsv')
        if os.path.exists(grpfile2chk):
            grpfile = grpfile2chk
        else:
            print >> sys.stderr, '[gbrs::stencil] A group file is *not* given. Genotype will be stenciled as is.'

    # Load alignment incidence matrix ('alnfile' is assumed to be in multiway transcriptome)
    alnmat = emase.AlignmentPropertyMatrix(h5file=alnfile, grpfile=grpfile)

    # Load genotype calls
    hid = dict(zip(alnmat.hname, np.arange(alnmat.num_haplotypes)))
    gid = dict(zip(alnmat.gname, np.arange(len(alnmat.gname))))
    gtmask = np.zeros((alnmat.num_haplotypes, alnmat.num_loci))
    gtcall_g = dict.fromkeys(alnmat.gname)
    with open(gtypefile) as fh:
        if grpfile is not None:
            gtcall_t = dict.fromkeys(alnmat.lname)
            for curline in dropwhile(is_comment, fh):
                item = curline.rstrip().split("\t")
                g, gt = item[:2]
                gtcall_g[g] = gt
                hid2set = np.array([hid[c] for c in gt])
                tid2set = np.array(alnmat.groups[gid[g]])
                gtmask[np.meshgrid(hid2set, tid2set)] = 1.0
                for t in tid2set:
                    gtcall_t[alnmat.lname[t]] = gt
        else:
            for curline in dropwhile(is_comment, fh):
                item = curline.rstrip().split("\t")
                g, gt = item[:2]
                gtcall_g[g] = gt
                hid2set = np.array([hid[c] for c in gt])
                gtmask[np.meshgrid(hid2set, gid[g])] = 1.0

    alnmat.multiply(gtmask, axis=2)
    for h in xrange(alnmat.num_haplotypes):
        alnmat.data[h].eliminate_zeros()
    outfile = kwargs.get('outfile')
    if outfile is None:
        outfile = 'gbrs.stenciled.' + os.path.basename(alnfile)
    alnmat.save(h5file=outfile)