Zhanxw at 03:18, 19 June 2012

2012-06-19T03:18:36Z

← Older revision		Revision as of 23:18, 18 June 2012
Line 89:		Line 89:
	fo = open(fout, 'w')		fo = open(fout, 'w')
	for ln in myopen(fin):		for ln in myopen(fin):
	f = ln.strip().split('\t')		f = ln.strip().split()
	fo.write('\t'.join(f[:6]))		l = len(f)
	fo.write('\t')		f = f[:6] + [ f[i2] + ' '+f[i2 +1] for i in xrange(3, l/2 )]
	if len(f[6:]) != len(flag):		fo.write('\t'.join(f[:6]))
	die('Inconsistent length of ped and map files')		fo.write('\t')
	newMarker = [m for i, m in enumerate(f[6:]) if flag[i]]		if len(f[6:]) != len(flag):
	fo.write('\t'.join(newMarker))		die('Inconsistent length of ped and map files')
	fo.write('\n')		newMarker = [m for i, m in enumerate(f[6:]) if flag[i]]
	#print marker[:10]		fo.write('\t'.join(newMarker))
	#die('test')		fo.write('\n')
			#print marker[:10]
			#die('test')
	return True		return True

Zhanxw: Created page with ' #!/usr/bin/python import sys, os def usage(): print("%s -m old.map -p old.ped -d old.dat -o new.prefix" % sys.argv[0] ) print("lift old.map (optio…'

2011-08-10T15:34:27Z

Created page with '<source lang="python"> #!/usr/bin/python import sys, os def usage(): print("%s -m old.map -p old.ped -d old.dat -o new.prefix" % sys.argv[0] ) print("lift old.map (optio…'

New page

<source lang="python">
#!/usr/bin/python
import sys, os

def usage():
print("%s -m old.map -p old.ped -d old.dat -o new.prefix" % sys.argv[0] )
print("lift old.map (optionally .ped and .dat) to new.prefix.map (.ped/.dat if possible)")

def die(msg):
print msg
sys.exit(2)

def myopen(fn):
import gzip
try:
h = gzip.open(fn)
ln = h.read(2) # read arbitrary bytes so check if @param fn is a gzipped file
except:
# cannot read in gzip format
return open(fn)
h.close()
return gzip.open(fn)

def map2bed(fin, fout):
fo = open(fout, 'w')
for ln in myopen(fin):
chrom, rs, mdist, pos = ln.split()
chrom = 'chr' + chrom
pos = int(pos)
fo.write('%s\t%d\t%d\t%s\n' % (chrom, pos-1, pos, rs))
fo.close()
return True

# global var:
LIFTED_SET = set()
UNLIFTED_SET = set()
def liftBed(fin, fout, funlifted):
params = dict()
params['LIFTOVER_BIN'] = '/net/wonderland/home/zhanxw/mycode/liftOver/liftOver'
params['OLD'] = fin
params['CHAIN'] = '/home/zhanxw/mycode/liftOver/hg18ToHg19.over.chain.gz'
params['NEW'] = fout
params['UNLIFTED'] = fout + '.unlifted'
from string import Template
cmd = Template('$LIFTOVER_BIN $OLD $CHAIN $NEW $UNLIFTED')
cmd = cmd.substitute(params)
os.system(cmd)
#record lifted/unliftd rs
for ln in myopen(params['UNLIFTED']):
if len(ln) == 0 or ln[0] == '#':continue
UNLIFTED_SET.add(ln.strip().split()[-1])
for ln in myopen(params['NEW']):
if len(ln) == 0 or ln[0] == '#':continue
LIFTED_SET.add(ln.strip().split()[-1])

return True

def bed2map(fin, fout):
fo = open(fout, 'w')
for ln in myopen(fin):
chrom, pos0, pos1, rs = ln.split()
chrom = chrom.replace('chr', '')
fo.write('%s\t%s\t0.0\t%s\n' % (chrom, rs, pos1))
fo.close()
return True

def liftDat(fin, fout):
fo = open(fout, 'w')
for ln in myopen(fin):
if len(ln) == 0 or ln[0] != 'M':
fo.write(ln)
else:
t, rs = ln.strip().split()
if rs in LIFTED_SET:
fo.write(ln)
fo.close()
return True

def liftPed(fin, fout, fOldMap):
# two ways to do it:
# 1. write unlifted snp list
# use PLINK to do this job using --exclude
# 2. alternatively, we can write our own method
# we will use method 2
marker = [i.strip().split()[1] for i in open(fOldMap)]
flag = map(lambda x: x not in UNLIFTED_SET, marker)
# print marker[:10]
# print flag[:10]
fo = open(fout, 'w')
for ln in myopen(fin):
f = ln.strip().split('\t')
fo.write('\t'.join(f[:6]))
fo.write('\t')
if len(f[6:]) != len(flag):
die('Inconsistent length of ped and map files')
newMarker = [m for i, m in enumerate(f[6:]) if flag[i]]
fo.write('\t'.join(newMarker))
fo.write('\n')
#print marker[:10]
#die('test')
return True

def makesure(result, succ_msg, fail_msg = "ERROR"):
if result:
print 'SUCC: ', succ_msg
else:
print 'FAIL: ', fail_msg
sys.exit(2)

if __name__ == '__main__':
import argparse
parser = argparse.ArgumentParser(description='')
parser.add_argument('-m', dest='mapFile', required = True)
parser.add_argument('-p', dest='pedFile')
parser.add_argument('-d', dest='datFile')
parser.add_argument('-o', dest='prefix', required = True)
args = parser.parse_args()

oldBed = args.mapFile + '.bed'
makesure(map2bed(args.mapFile, oldBed),
'map->bed succ')

newBed = args.prefix + '.bed'
unlifted = args.prefix + '.unlifted'
makesure(liftBed(oldBed, newBed, unlifted),
'liftBed succ')

newMap = args.prefix + '.map'
makesure(bed2map(newBed, newMap),
'bed->map succ')

if args.datFile:
newDat = args.prefix + '.dat'
makesure(liftDat(args.datFile, newDat),
'liftDat succ')

if args.pedFile:
newPed = args.prefix + '.ped'
makesure(liftPed(args.pedFile, newPed, args.mapFile),
'liftPed succ')

</source>

LiftMap.py - Revision history

Zhanxw at 03:18, 19 June 2012

Zhanxw: Created page with ' #!/usr/bin/python import sys, os def usage(): print("%s -m old.map -p old.ped -d old.dat -o new.prefix" % sys.argv[0] ) print("lift old.map (optio…'