The following functions are from last time:
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.
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
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)
Find the closest "k" markers to the specified "geneId"
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.
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
Set up to plot.
%matplotlib inline
import matplotlib.pyplot as plot
import numpy
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')