def get_mds(conn, metadata, args):
    idx_to_sample = {}
    res = conn.execute(sql.text("select sample_id, name from samples"))
    for row in res:
        idx_to_sample[int(row['sample_id']) - 1] = row['name']

    query = "SELECT DISTINCT v.variant_id, v.gt_types\
    FROM variants v\
    WHERE v.type = 'snp'"
    res = conn.execute(query)

    # keep a list of numeric genotype values
    # for each sample
    genotypes = collections.defaultdict(list)
    import zlib
    unpack = Z.unpack_genotype_blob

    for row in res:

        try:
            gt_types = unpack(row['gt_types'])
        except zlib.error:
            unpack = Z.snappy_unpack_blob
            gt_types = unpack(row['gt_types'])

        # at this point, gt_types is a numpy array
        # idx:  0 1 2 3 4 5 6 .. #samples
        # type [0 1 2 1 2 0 0 ..         ]
        for idx, gt_type in enumerate(gt_types):
            sample = idx_to_sample[idx]
            genotypes[sample].append(gt_type)

    mds = collections.defaultdict(float)
    # convert the genotype list for each sample
    # to a numpy array for performance.
    # masks stores an array of T/F indicating which genotypes are
    # known (True, [0,1,2]) and unknown (False [-1]).
    masks = {}
    for s in genotypes:
        sample = str(s)
        x = np.array(genotypes[sample])
        genotypes[sample] = x
        masks[sample] = \
            np.ma.masked_where(genotypes[sample] != UNKNOWN,
                               genotypes[sample]).mask

    # compute the euclidean distance for each s1/s2 combination
    # using numpy's vectorized sum() and square() operations.
    # we use the mask arrays to identify the indices of known genotypes
    # for each sample.  by doing a bitwise AND of the mask arrays for the
    # two samples, we have a mask array of variants where __both__ samples
    # were called.
    for sample1 in genotypes:
        for sample2 in genotypes:
            pair = (sample1, sample2)
            # which variants have known genotypes for both samples?
            both_mask = masks[str(sample1)] & masks[str(sample2)]
            genotype1 = genotypes[sample1]
            genotype2 = genotypes[sample2]

            # distance between s1 and s2:
            eucl_dist = float(np.sum(np.square((genotype1 - genotype2)[both_mask]))) \
                / \
                float(np.sum(both_mask))

            mds[pair] = eucl_dist

    # report the pairwise MDS for each sample pair.
    print "sample1\tsample2\tdistance"
    for pair in mds:
        print "\t".join([str(pair[0]), str(pair[1]), str(round(mds[pair], 4))])
Compute the pairwise genetic distance between each sample.
dataset/ETHPy150Open arq5x/gemini/gemini/gemini_stats.py/get_mds