-
Notifications
You must be signed in to change notification settings - Fork 0
/
3dna.py
executable file
·130 lines (93 loc) · 2.9 KB
/
3dna.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
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
#!/usr/bin/python
import sys
import math
# display usage if arguments are missing
if len(sys.argv) < 3:
print 'usage: 3dna.py <input file> <resolution>'
sys.exit(0)
# open specified genome file
filename=sys.argv[1]
# if specified, only process a subset of the file (default to 100)
chromosome_sample_count = 0
chromosome_sample_limit = 100
if len(sys.argv) > 2:
chromosome_sample_limit = int(sys.argv[2])
# globals
current_chromosome = ""
current_chromosome_lines = 0
chromosomes = {}
points = []
vectors = []
radius = 200
# open file for analysis
gfile = open(filename)
for line in gfile:
# ignore comments
if line[0] != '#':
# split line
row = line.split('\t')
# count chromosomes
if row[1] != current_chromosome:
# if this is the end of the interesting chromosomes, exit analysis
if not row[1].isdigit():
chromosomes[current_chromosome] = current_chromosome_lines
break
# skip output if this is the first chromosome
if current_chromosome != "":
chromosomes[current_chromosome] = current_chromosome_lines
current_chromosome_lines = 0
current_chromosome = row[1]
# count lines
if current_chromosome_lines > chromosome_sample_limit:
current_chromosome_lines = current_chromosome_lines + 1
else:
current_chromosome_lines = chromosome_sample_limit
# close & reopen file for actual processing
gfile.close()
gfile = open(filename)
current_chromosome = ""
angle = 0
for line in gfile:
# ignore comments
if line[0] != '#':
# split line
row = line.split('\t')
# increment chromosome
if row[1] != current_chromosome:
# output last chromosome
if current_chromosome != "":
print "hull(){"
print "cylinder(r=50,h=20);"
print "translate(v=[0,0,%d]){" % (int(current_chromosome) * 10)
print "linear_extrude(height=%d,center=true,convexity=10,twist=-100){" % 20 #(int(current_chromosome) * 5)
print "polygon(points=%s,paths=[%s]);" % (points, vectors)
print "}}}"
# reset the sample count
chromosome_sample_count = 0
# if the next row isn't a chromosome we handle, eject
if not row[1].isdigit():
break
else:
# create a new layer
current_chromosome = row[1]
points = []
vector = 0
vectors = []
# get angle inc. from chromosome dictionary
angle_inc = (3.14 * 2) / chromosomes[current_chromosome]
# if we've hit the sample limit, move along
if chromosome_sample_count < chromosome_sample_limit:
# turn marker into integer
point_val = int(ord(row[3][0])) + int(ord(row[3][1]))
# convert to x,y
point_val_x = round((radius - point_val) * math.cos(angle),1)
point_val_y = round((radius - point_val) * math.sin(angle),1)
point = [point_val_x,point_val_y]
points.append(point)
vectors.append(vector)
# increment angle
angle = angle + angle_inc
# increment vector
vector = vector + 1
# increment sample count
chromosome_sample_count = chromosome_sample_count + 1