The following functions are from last time:

In [1]: 
def GetAlleles(filename):
    """ Read the alleles in filename """ 
    allele = dict()
    with open(filename, 'rb') as tsvfile:
        for row in tsvfile:
            line = row[:row.find('#')].strip(' \r\n')   # strip comments and line endings
            if (len(line) == 0):                        # skip empty lines
                continue
            field = line.split('\t')
            if (len(field) != 4):                       # print error for lines missing fields
                print "ERROR: %s" % str(field)
                continue
            allele[field[0]] = field[3]
    return allele

def GetCoordinates(filename):
    """ Read the chromosome and positions of alleles in filename """ 
    chrPos = dict()
    with open(filename, 'rb') as tsvfile:
        for row in tsvfile:
            line = row[:row.find('#')].strip(' \r\n')   # strip comments and line endings
            if (len(line) == 0):                        # skip empty lines
                continue
            field = line.split('\t')
            if (len(field) != 4):                       # print error for lines missing fields
                print "ERROR: %s" % str(field)
                continue
            chrPos[field[0]] = (field[1], int(field[2]))
    return chrPos

This is the phasing routine from last time.

In [5]: 
phaseRule = {
    ('--', '--', 'AA'): 'AA', 
    ('--', '--', 'AC'): 'ac', 
    ('--', '--', 'AG'): 'ag', 
    ('--', '--', 'AT'): 'at', 
    ('--', '--', 'CC'): 'CC', 
    ('--', '--', 'CG'): 'cg', 
    ('--', '--', 'CT'): 'ct', 
    ('--', '--', 'DD'): 'DD',
    ('--', '--', 'GG'): 'GG',
    ('--', '--', 'GT'): 'gt', 
    ('--', '--', 'II'): 'II',
    ('--', '--', 'TT'): 'TT', 
    ('--', 'AA', '--'): '-A',
    ('--', 'AA', 'AA'): 'AA', 
    ('--', 'AA', 'AC'): 'CA', 
    ('--', 'AA', 'AG'): 'GA',
    ('--', 'AA', 'AT'): 'TA', 
    ('--', 'AC', 'AA'): 'AA', 
    ('--', 'AC', 'AC'): 'ac',
    ('--', 'AC', 'CC'): 'CC', 
    ('--', 'AG', 'AA'): 'AA', 
    ('--', 'AG', 'AG'): 'ag',
    ('--', 'AG', 'GG'): 'GG',
    ('--', 'AT', 'AA'): 'AA',
    ('--', 'AT', 'AT'): 'at',
    ('--', 'AT', 'TT'): 'TT', 
    ('--', 'CC', '--'): '-C',
    ('--', 'CC', 'AC'): 'AC', 
    ('--', 'CC', 'CC'): 'CC', 
    ('--', 'CC', 'CG'): 'GC', 
    ('--', 'CC', 'CT'): 'TC',
    ('--', 'CG', 'CC'): 'CC',
    ('--', 'CG', 'CG'): 'cg',
    ('--', 'CG', 'GG'): 'GG',
    ('--', 'CT', 'CC'): 'CC', 
    ('--', 'CT', 'CT'): 'ct',
    ('--', 'CT', 'TT'): 'TT',
    ('--', 'DD', 'DD'): 'DD',
    ('--', 'GG', '--'): '-G', 
    ('--', 'GG', 'AG'): 'AG',
    ('--', 'GG', 'CG'): 'CG',
    ('--', 'GG', 'GG'): 'GG', 
    ('--', 'GG', 'GT'): 'TG', 
    ('--', 'GT', 'GG'): 'GG',
    ('--', 'GT', 'GT'): 'gt', 
    ('--', 'GT', 'TT'): 'TT', 
    ('--', 'II', '--'): '-I',
    ('--', 'II', 'II'): 'II', 
    ('--', 'TT', '--'): '-T', 
    ('--', 'TT', 'AT'): 'AT',
    ('--', 'TT', 'CT'): 'CT',
    ('--', 'TT', 'GT'): 'GT', 
    ('--', 'TT', 'TT'): 'TT', 
    ('AA', '--', '--'): 'A-',
    ('AA', '--', 'AA'): 'AA', 
    ('AA', '--', 'AC'): 'AC', 
    ('AA', '--', 'AG'): 'AG',
    ('AA', '--', 'AT'): 'AT', 
    ('AA', 'AA', '--'): 'AA', 
    ('AA', 'AA', 'AA'): 'AA',
    ('AA', 'AC', '--'): 'A-',
    ('AA', 'AC', 'AA'): 'AA',
    ('AA', 'AC', 'AC'): 'AC',
    ('AA', 'AG', '--'): 'A-',
    ('AA', 'AG', 'AA'): 'AA',
    ('AA', 'AG', 'AG'): 'AG',
    ('AA', 'AT', '--'): 'A-',
    ('AA', 'AT', 'AA'): 'AA',
    ('AA', 'AT', 'AT'): 'AT',
    ('AA', 'CC', '--'): 'AC',
    ('AA', 'CC', 'AC'): 'AC',
    ('AA', 'GG', '--'): 'AG',
    ('AA', 'GG', 'AG'): 'AG',
    ('AA', 'TT', '--'): 'AT',
    ('AA', 'TT', 'AT'): 'AT',
    ('AC', '--', 'AA'): 'AA',
    ('AC', '--', 'AC'): 'ac',
    ('AC', '--', 'CC'): 'CC',
    ('AC', 'AA', '--'): '-A',
    ('AC', 'AA', 'AA'): 'AA',
    ('AC', 'AA', 'AC'): 'CA',
    ('AC', 'AC', 'AA'): 'AA',
    ('AC', 'AC', 'AC'): 'ac',
    ('AC', 'AC', 'CC'): 'CC',
    ('AC', 'CC', '--'): '-C',
    ('AC', 'CC', 'AC'): 'AC',
    ('AC', 'CC', 'CC'): 'CC',
    ('AG', '--', 'AA'): 'AA',
    ('AG', '--', 'AG'): 'ag',
    ('AG', '--', 'GG'): 'GG',
    ('AG', 'AA', '--'): '-A',
    ('AG', 'AA', 'AA'): 'AA',
    ('AG', 'AA', 'AG'): 'GA',
    ('AG', 'AG', 'AA'): 'AA',
    ('AG', 'AG', 'AG'): 'ag',
    ('AG', 'AG', 'GG'): 'GG',
    ('AG', 'GG', '--'): '-G',
    ('AG', 'GG', 'AG'): 'AG',
    ('AG', 'GG', 'GG'): 'GG',
    ('AT', '--', 'AA'): 'AA',
    ('AT', '--', 'AT'): 'at',
    ('AT', '--', 'TT'): 'TT',
    ('AT', 'AA', '--'): '-A',
    ('AT', 'AA', 'AA'): 'AA',
    ('AT', 'AA', 'AT'): 'TA',
    ('AT', 'AT', 'AA'): 'AA',
    ('AT', 'AT', 'AT'): 'at',
    ('AT', 'AT', 'TT'): 'TT',
    ('AT', 'TT', '--'): '-T',
    ('AT', 'TT', 'AT'): 'AT',
    ('AT', 'TT', 'TT'): 'TT',
    ('CC', '--', '--'): 'C-',
    ('CC', '--', 'AC'): 'CA',
    ('CC', '--', 'CC'): 'CC',
    ('CC', '--', 'CG'): 'CG',
    ('CC', '--', 'CT'): 'CT',
    ('CC', 'AA', '--'): 'CA',
    ('CC', 'AA', 'AC'): 'CA',
    ('CC', 'AC', '--'): 'C-',
    ('CC', 'AC', 'AC'): 'CA',
    ('CC', 'AC', 'CC'): 'CC',
    ('CC', 'CC', '--'): 'CC',
    ('CC', 'CC', 'CC'): 'CC',
    ('CC', 'CG', '--'): 'C-',
    ('CC', 'CG', 'CC'): 'CC',
    ('CC', 'CG', 'CG'): 'CG',
    ('CC', 'CT', '--'): 'C-',
    ('CC', 'CT', 'CC'): 'CC',
    ('CC', 'CT', 'CT'): 'CT',
    ('CC', 'GG', '--'): 'CG',
    ('CC', 'GG', 'CG'): 'CG',
    ('CC', 'TT', '--'): 'CT',
    ('CC', 'TT', 'CT'): 'CT',
    ('CG', '--', 'CC'): 'CC',
    ('CG', '--', 'CG'): 'cg',
    ('CG', '--', 'GG'): 'GG',
    ('CG', 'CC', 'CC'): 'CC',
    ('CG', 'CC', 'CG'): 'GC',
    ('CG', 'CG', 'CC'): 'CC',
    ('CG', 'CG', 'CG'): 'cg',
    ('CG', 'CG', 'GG'): 'GG',
    ('CG', 'GG', '--'): '-G',
    ('CG', 'GG', 'CG'): 'CG',
    ('CG', 'GG', 'GG'): 'GG',
    ('CT', '--', 'CC'): 'CC',
    ('CT', '--', 'CT'): 'ct',
    ('CT', '--', 'TT'): 'TT',
    ('CT', 'CC', '--'): '-C',
    ('CT', 'CC', 'CC'): 'CC',
    ('CT', 'CC', 'CT'): 'TC',
    ('CT', 'CT', 'CC'): 'CC',
    ('CT', 'CT', 'CT'): 'ct',
    ('CT', 'CT', 'TT'): 'TT',
    ('CT', 'TT', '--'): '-T',
    ('CT', 'TT', 'CT'): 'CT',
    ('CT', 'TT', 'TT'): 'TT',
    ('DD', '--', 'DD'): 'DD',
    ('DD', 'DD', '--'): 'DD',
    ('DD', 'DD', 'DD'): 'DD',
    ('DD', 'DI', 'DI'): 'DI',
    ('DD', 'DI', 'DD'): 'DD',
    ('DD', 'II', 'DI'): 'DI',
    ('DI', 'DD', 'DD'): 'DD',
    ('DI', 'DD', 'DI'): 'DI',
    ('DI', 'DI', 'DD'): 'DD',
    ('DI', 'DI', 'DI'): 'di',
    ('DI', 'DI', 'II'): 'II',
    ('DI', 'II', 'II'): 'II',
    ('DI', 'II', 'DI'): 'DI',
    ('GG', '--', '--'): 'G-',
    ('GG', '--', 'AG'): 'GA',
    ('GG', '--', 'CG'): 'GC',
    ('GG', '--', 'GG'): 'GG',
    ('GG', '--', 'GT'): 'GT',
    ('GG', 'AA', '--'): 'GA',
    ('GG', 'AA', 'AG'): 'GA',
    ('GG', 'AG', '--'): 'G-',
    ('GG', 'AG', 'AG'): 'GA',
    ('GG', 'AG', 'GG'): 'GG',
    ('GG', 'CC', 'CG'): 'GC',
    ('GG', 'CG', 'CG'): 'GC',
    ('GG', 'CG', 'GG'): 'GG',
    ('GG', 'GG', '--'): 'GG',
    ('GG', 'GG', 'GG'): 'GG',
    ('GG', 'GT', '--'): 'G-',
    ('GG', 'GT', 'GG'): 'GG',
    ('GG', 'GT', 'GT'): 'GT',
    ('GG', 'TT', '--'): 'GT',
    ('GG', 'TT', 'GT'): 'GT',
    ('GT', '--', 'GG'): 'GG',
    ('GT', '--', 'GT'): 'gt',
    ('GT', '--', 'TT'): 'TT',
    ('GT', 'GG', '--'): '-G',
    ('GT', 'GG', 'GG'): 'GG',
    ('GT', 'GG', 'GT'): 'TG',
    ('GT', 'GT', 'GG'): 'GG',
    ('GT', 'GT', 'GT'): 'gt',
    ('GT', 'GT', 'TT'): 'TT',
    ('GT', 'TT', '--'): '-T',
    ('GT', 'TT', 'GT'): 'GT',
    ('GT', 'TT', 'TT'): 'TT',
    ('II', '--', 'II'): 'II',
    ('II', 'DD', 'DI'): 'ID',
    ('II', 'DI', 'DI'): 'ID',
    ('II', 'DI', 'II'): 'II',
    ('II', 'II', '--'): 'II',
    ('II', 'II', 'II'): 'II',
    ('TT', '--', '--'): 'T-',
    ('TT', '--', 'AT'): 'TA',
    ('TT', '--', 'CT'): 'TC',
    ('TT', '--', 'GT'): 'TG',
    ('TT', '--', 'TT'): 'TT',
    ('TT', 'AA', '--'): 'TA',
    ('TT', 'AA', 'AT'): 'TA',
    ('TT', 'AT', '--'): 'T-',
    ('TT', 'AT', 'AT'): 'TA',
    ('TT', 'AT', 'TT'): 'TT',
    ('TT', 'CC', '--'): 'TC',
    ('TT', 'CC', 'CT'): 'TC',
    ('TT', 'CT', '--'): 'T-',
    ('TT', 'CT', 'CT'): 'TC',
    ('TT', 'CT', 'TT'): 'TT',
    ('TT', 'GG', '--'): 'TG',
    ('TT', 'GG', 'GT'): 'TG',
    ('TT', 'GT', '--'): 'T-',
    ('TT', 'GT', 'GT'): 'TG',
    ('TT', 'GT', 'TT'): 'TT',
    ('TT', 'TT', '--'): 'TT',
    ('TT', 'TT', 'TT'): 'TT',
    ('--', 'A', '--'): '-A',
    ('--', 'A', 'AA'): 'AA',
    ('--', 'A', 'AC'): 'CA',
    ('--', 'A', 'AG'): 'GA',
    ('--', 'A', 'AT'): 'TA',
    ('--', 'C', '--'): '-C',
    ('--', 'C', 'AC'): 'AC',
    ('--', 'C', 'CC'): 'CC',
    ('--', 'C', 'CG'): 'GC',
    ('--', 'C', 'CT'): 'TC',
    ('--', 'D', 'DD'): 'DD',
    ('--', 'G', '--'): '-G',
    ('--', 'G', 'AG'): 'AG',
    ('--', 'G', 'CG'): 'CG',
    ('--', 'G', 'GG'): 'GG',
    ('--', 'G', 'GT'): 'TG',
    ('--', 'T', '--'): '-T',
    ('--', 'T', 'AT'): 'AT',
    ('--', 'T', 'CT'): 'CT',
    ('--', 'T', 'GT'): 'GT',
    ('--', 'T', 'TT'): 'TT',
    ('AA', 'A', 'AA'): 'AA',
    ('AA', 'C', 'AC'): 'AC',
    ('AA', 'G', 'AG'): 'AG',
    ('AA', 'T', 'AT'): 'AT',
    ('AC', 'A', 'AA'): 'AA',
    ('AC', 'A', 'AC'): 'CA',
    ('AC', 'C', 'AC'): 'AC',
    ('AC', 'C', 'CC'): 'CC',
    ('AG', 'A', 'AA'): 'AA',
    ('AG', 'A', 'AG'): 'GA',
    ('AG', 'G', 'AG'): 'AG',
    ('AG', 'G', 'GG'): 'GG',
    ('AT', 'A', 'AA'): 'AA',
    ('AT', 'A', 'AT'): 'TA',
    ('AT', 'T', 'AT'): 'AT',
    ('AT', 'T', 'TT'): 'TT',
    ('CC', 'A', 'AC'): 'CA',
    ('CC', 'C', 'CC'): 'CC',
    ('CC', 'G', 'CG'): 'CG',
    ('CC', 'T', 'CT'): 'CT',
    ('CG', 'C', 'CC'): 'CC',
    ('CG', 'C', 'CG'): 'GC',
    ('CG', 'G', 'CG'): 'CG',
    ('CG', 'G', 'GG'): 'GG',
    ('CT', 'C', 'CC'): 'CC',
    ('CT', 'C', 'CT'): 'TC',
    ('CT', 'T', 'CT'): 'CT',
    ('CT', 'T', 'TT'): 'TT',
    ('DD', 'D', 'DD'): 'DD',
    ('GG', 'A', 'AG'): 'GA',
    ('GG', 'C', 'CG'): 'GC',
    ('GG', 'G', 'GG'): 'GG',
    ('GG', 'T', 'GT'): 'GT',
    ('GT', 'G', 'GG'): 'GG',
    ('GT', 'G', 'GT'): 'TG',
    ('GT', 'T', 'GT'): 'GT',
    ('GT', 'T', 'TT'): 'TT',
    ('II', 'I', 'II'): 'II',
    ('TT', 'A', 'AT'): 'TA',
    ('TT', 'C', 'CT'): 'TC',
    ('TT', 'G', 'GT'): 'TG',
    ('TT', 'T', 'TT'): 'TT',
    ('--', 'MT', 'A'): 'A',
    ('--', 'MT', 'C'): 'C',
    ('--', 'MT', 'D'): 'D',
    ('--', 'MT', 'G'): 'G',
    ('--', 'MT', 'I'): 'I',
    ('--', 'MT', 'T'): 'T',
    ('A', 'MT', '--'): 'A',
    ('C', 'MT', '--'): 'C',
    ('D', 'MT', '--'): 'D',
    ('G', 'MT', '--'): 'G',
    ('I', 'MT', '--'): 'I',
    ('T', 'MT', '--'): 'T',
    ('A', 'MT', 'A'): 'A',
    ('C', 'MT', 'C'): 'C',
    ('D', 'MT', 'D'): 'D',
    ('G', 'MT', 'G'): 'G',
    ('I', 'MT', 'I'): 'I',
    ('T', 'MT', 'T'): 'T',
    ('--', 'X', 'A'): 'A',
    ('--', 'X', 'C'): 'C',
    ('--', 'X', 'D'): 'D',
    ('--', 'X', 'G'): 'G',
    ('--', 'X', 'I'): 'I',
    ('--', 'X', 'T'): 'T',
    ('AA', 'X', 'A'): 'A',
    ('AC', 'X', 'A'): 'A',
    ('AC', 'X', 'C'): 'C',
    ('AG', 'X', 'A'): 'A',
    ('AG', 'X', 'G'): 'G',
    ('AT', 'X', 'A'): 'A',
    ('AT', 'X', 'T'): 'T',
    ('CC', 'X', 'C'): 'C',
    ('CG', 'X', 'C'): 'C',
    ('CG', 'X', 'G'): 'G',
    ('CT', 'X', 'C'): 'C',
    ('CT', 'X', 'T'): 'T',
    ('DD', 'X', 'D'): 'D',
    ('GG', 'X', 'G'): 'G',
    ('GT', 'X', 'G'): 'G',
    ('GT', 'X', 'T'): 'T',
    ('II', 'X', 'I'): 'I',
    ('TT', 'X', 'T'): 'T',
    ('Y', '--', 'A'): 'A',
    ('Y', '--', 'C'): 'C',
    ('Y', '--', 'D'): 'D',
    ('Y', '--', 'G'): 'G',
    ('Y', '--', 'I'): 'I',
    ('Y', '--', 'T'): 'T',
    ('Y', 'A', '--'): 'A',
    ('Y', 'C', '--'): 'C',
    ('Y', 'D', '--'): 'D',
    ('Y', 'G', '--'): 'G',
    ('Y', 'I', '--'): 'I',
    ('Y', 'T', '--'): 'T',
    ('Y', 'A', 'A'): 'A',
    ('Y', 'C', 'C'): 'C',
    ('Y', 'D', 'D'): 'D',
    ('Y', 'G', 'G'): 'G',
    ('Y', 'I', 'I'): 'I',
    ('Y', 'T', 'T'): 'T'
}

invalidTrio = {}

def phase(mom, pop, kid, ChrPos):
    """ Use trio alleles to phase kid. Returns a new kid allele dictionary """
    modified = 0
    unchanged = 0
    invalid = 0
    for gId, call in kid.iteritems():
        try:
            chromo = ChrPos[gId][0]
            momA = mom[gId] if (chromo != 'Y') else 'Y'
            popA = pop[gId] if (chromo != 'MT') else 'MT'
        except KeyError:
            # print gId
            continue
        if ((chromo == 'X') and (len(call) == 1)):
            popA = 'X'
        try: 
            newCall = phaseRule[momA, popA, call]
        except KeyError:
            if (call != '--'):
                invalidTrio[momA, popA, call] = invalidTrio.get((momA, popA, call), 0) + 1
                invalid += 1
            newCall = '--'
        if (newCall != call):
            kid[gId] = newCall
            modified += 1
        else:
            unchanged += 1
    print modified, unchanged, invalid

Let's read in all genomes and phase all valid trios

In [6]: 
alleles = dict()
for sample in "ABCDEFGHI":
    alleles[sample] = GetAlleles("genome%s.txt" % sample)
coordinates = GetCoordinates("genomeA.txt")

# Note the order matters here since the child's genotypes are changed by phase
phase(alleles['F'], alleles['G'], alleles['H'], coordinates)
phase(alleles['F'], alleles['G'], alleles['I'], coordinates)
phase(alleles['D'], alleles['E'], alleles['G'], coordinates)

172563 787018 432
171188 788392 483
168617 790963 303

Find the closest "k" markers to the specified "geneId"

In [8]: 
def kNearest(k, geneId, coordinates):
    ''' returns tuples of (dist, id) for the specified geneId's k-nearest neighbors '''
    neighbor = [(2**31, 'foo') for i in xrange(k)]
    chromo, position = coordinates[geneId]
    for gid, chrpos in coordinates.iteritems():
        c, p = chrpos
        if (c != chromo):
            continue
        if (gid == geneId):
            continue
        delta = abs(p - position)
        if (delta < neighbor[-1][0]):
            neighbor.pop()
            neighbor.append((delta, gid))
            neighbor.sort()
    return neighbor

Let's find a het/het/het combination to infer. Recall that in class we tested with the 500th het/het/het, and I asked the class to see f they could resolve the phase of the 310th and 10006th. Once you figure those out, find one of your own to describe in class.

In [43]: 
N = 0
for gId, call in alleles['G'].iteritems():
    if (call.islower()):
        N += 1
        if (N == 10006):  # 310, 500, 10006
            target = gId
            
print "%d unresolved hets" % N
print target, alleles['G'][target], coordinates[target]
knn = kNearest(10, target, coordinates)
knn.append((0,target))
print knn

53446 unresolved hets
rs1459283 ct ('4', 21708295)
[(390, 'rs4235291'), (695, 'rs1459282'), (1643, 'rs7669417'), (3755, 'rs1503976'), (6739, 'rs6850480'), (6994, 'rs2167246'), (11001, 'rs9884923'), (14570, 'rs9291429'), (18566, 'rs11721958'), (18865, 'rs1354390'), (0, 'rs1459283')]

Set up to plot.

In [44]: 
%matplotlib inline
import matplotlib.pyplot as plot
import numpy
In [45]: 
cmap = { 'A':'r', 'C':'y', 'G':'b', 'T':'c', 'D':'g', 'I':'m', '-':'k' }
plot.figure(figsize=(10,6))
chromo, position = coordinates[target]
mindiff = 0
for i, sample in enumerate("ABCDEFGHI"):
    for delta, gId in knn:
        call = alleles[sample][gId]
        diff = coordinates[gId][1]-position
        offset = 0.2 if ((sample in 'GHI') and (not call.islower())) else 0.1
        glyph = 's' if (sample in 'GHI') else 'o'
        plot.plot([diff], [9-i+offset], cmap[call[0].upper()]+glyph)
        plot.plot([diff], [9-i-offset], cmap[call[-1].upper()]+glyph)
        mindiff = min(diff, mindiff)
plot.title("%s: Chr %s, %d" % (target, chromo, position))
for i, sample in enumerate("ABCDEFGHI"):
    plot.text(mindiff-600, 8.9-i, sample, color='k')
In []: