-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathbp2cM.py
45 lines (41 loc) · 1.13 KB
/
bp2cM.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
'''
Author: Dimitris Arnellos
Purpose: Converts base pair lengths to centimorgans
Usage: python bp2cM.py ibdResultsFile plink.map > out.file
'''
import sys, re
firstIndiv = []
secondIndiv = []
segStart = []
segEnd = []
with open(sys.argv[1], 'r') as f:
for line in f:
line = re.match('(\S+)\s+\S+\s+(\S+)\s+\S+\s+\S+\s+(\S+)\s+(\S+)', line)
firstIndiv.append(line.group(1))
secondIndiv.append(line.group(2))
segStart.append(int(line.group(3)))
segEnd.append(int(line.group(4)))
cM = []
bp = []
with open(sys.argv[2], 'r') as f:
for line in f:
line = re.match('\S+\s+\S+\s+(\S+)\s+(\S+)', line)
cM.append(float(line.group(1)))
bp.append(int(line.group(2)))
for i in range(0, len(segEnd)):
cmStart = ""
cmEnd = ""
for j in range(0, len(cM)):
if (j == 0 and segStart[i] < bp[j]):
cmStart = cM[j]
elif (cmStart == "" and segStart[i] < bp[j]):
cmStart = (cM[j-1] + cM[j])/2
elif (cmStart == "" and segStart[i] == bp[j]):
cmStart = cM[j]
if (segEnd[i] < bp[j]):
cmEnd = (cM[j+1] + cM[j])/2
break
elif (segEnd[i] == bp[j]):
cmEnd = cM[j]
break
print(firstIndiv[i], secondIndiv[i], cmEnd - cmStart, sep="\t")