forked from reedacartwright/novo-muta
-
Notifications
You must be signed in to change notification settings - Fork 1
/
driver.py
54 lines (47 loc) · 1.64 KB
/
driver.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
46
47
48
49
50
51
52
53
54
#!/usr/bin/env python
"""
Driver file for TrioModel. This accepts an input text file where each parameter
is deliminated by a tab and each model object is on a new line. Leave an
optional field empty to use the default value. Parameters include:
#A #C # G #T child read counts (each count is deliminated by a tab)
#A #C # G #T mother read counts
#A #C # G #T father read counts
population mutation rate
germline mutation rate
somatic mutation rate
sequencing error rate
dirichlet multinomial dispersion (optional)
dirichlet multinomial bias (optional)
The output is formated with the read counts followed by the probability.
"""
import sys
from family.trio_model import TrioModel
from family import utilities as ut
# run python driver.py <parameters.txt>
handle = open(sys.argv[1])
for line in handle:
values = line.strip('\n').split('\t')
child_read_arr = values[:4]
mom_read_arr = values[4:8]
dad_read_arr = values[8:12]
child_read = [int(count) for count in child_read_arr]
mom_read = [int(count) for count in mom_read_arr]
dad_read = [int(count) for count in dad_read_arr]
reads = [child_read, mom_read, dad_read]
rates_arr = values[12:16]
rates = [float(rate) for rate in rates_arr]
disp = float(values[16]) if values[16] else 1000 # default dispersion value
bias = float(values[17]) if values[17] else None
trio_model = TrioModel(
reads=reads,
pop_muta_rate=rates[0],
germ_muta_rate=rates[1],
soma_muta_rate=rates[2],
seq_err_rate=rates[3],
dm_disp=disp,
dm_bias=bias
)
proba = trio_model.trio()
print(reads, end='\t')
print(proba)
handle.close()