|
0
|
1 #!/usr/bin/env python
|
|
|
2 #Guruprasad Ananda
|
|
|
3 """
|
|
|
4 Least Common Ancestor tool.
|
|
|
5 """
|
|
|
6 import sys, string, re, commands, tempfile, random
|
|
|
7
|
|
|
8 def stop_err(msg):
|
|
|
9 sys.stderr.write(msg)
|
|
|
10 sys.exit()
|
|
|
11
|
|
|
12 def main():
|
|
|
13 try:
|
|
|
14 inputfile = sys.argv[1]
|
|
|
15 outfile = sys.argv[2]
|
|
|
16 rank_bound = int( sys.argv[3] )
|
|
|
17 """
|
|
|
18 Mapping of ranks:
|
|
|
19 root :2,
|
|
|
20 superkingdom:3,
|
|
|
21 kingdom :4,
|
|
|
22 subkingdom :5,
|
|
|
23 superphylum :6,
|
|
|
24 phylum :7,
|
|
|
25 subphylum :8,
|
|
|
26 superclass :9,
|
|
|
27 class :10,
|
|
|
28 subclass :11,
|
|
|
29 superorder :12,
|
|
|
30 order :13,
|
|
|
31 suborder :14,
|
|
|
32 superfamily :15,
|
|
|
33 family :16,
|
|
|
34 subfamily :17,
|
|
|
35 tribe :18,
|
|
|
36 subtribe :19,
|
|
|
37 genus :20,
|
|
|
38 subgenus :21,
|
|
|
39 species :22,
|
|
|
40 subspecies :23,
|
|
|
41 """
|
|
|
42 except:
|
|
|
43 stop_err("Syntax error: Use correct syntax: program infile outfile")
|
|
|
44
|
|
|
45 fin = open(sys.argv[1],'r')
|
|
|
46 for j, line in enumerate( fin ):
|
|
|
47 elems = line.strip().split('\t')
|
|
|
48 if len(elems) < 24:
|
|
|
49 stop_err("The format of the input dataset is incorrect. Taxonomy datatype should contain at least 24 columns.")
|
|
|
50 if j > 30:
|
|
|
51 break
|
|
|
52 cols = range(1,len(elems))
|
|
|
53 fin.close()
|
|
|
54
|
|
|
55 group_col = 0
|
|
|
56 tmpfile = tempfile.NamedTemporaryFile()
|
|
|
57
|
|
|
58 try:
|
|
|
59 """
|
|
|
60 The -k option for the Posix sort command is as follows:
|
|
|
61 -k, --key=POS1[,POS2]
|
|
|
62 start a key at POS1, end it at POS2 (origin 1)
|
|
|
63 In other words, column positions start at 1 rather than 0, so
|
|
|
64 we need to add 1 to group_col.
|
|
|
65 if POS2 is not specified, the newer versions of sort will consider the entire line for sorting. To prevent this, we set POS2=POS1.
|
|
|
66 """
|
|
|
67 command_line = "sort -f -k " + str(group_col+1) +"," + str(group_col+1) + " -o " + tmpfile.name + " " + inputfile
|
|
|
68 except Exception, exc:
|
|
|
69 stop_err( 'Initialization error -> %s' %str(exc) )
|
|
|
70
|
|
|
71 error_code, stdout = commands.getstatusoutput(command_line)
|
|
|
72
|
|
|
73 if error_code != 0:
|
|
|
74 stop_err( "Sorting input dataset resulted in error: %s: %s" %( error_code, stdout ))
|
|
|
75
|
|
|
76 prev_item = ""
|
|
|
77 prev_vals = []
|
|
|
78 remaining_vals = []
|
|
|
79 skipped_lines = 0
|
|
|
80 fout = open(outfile, "w")
|
|
|
81 block_valid = False
|
|
|
82
|
|
|
83
|
|
|
84 for ii, line in enumerate( file( tmpfile.name )):
|
|
|
85 if line and not line.startswith( '#' ) and len(line.split('\t')) >= 24: #Taxonomy datatype should have at least 24 columns
|
|
|
86 line = line.rstrip( '\r\n' )
|
|
|
87 try:
|
|
|
88 fields = line.split("\t")
|
|
|
89 item = fields[group_col]
|
|
|
90 if prev_item != "":
|
|
|
91 # At this level, we're grouping on values (item and prev_item) in group_col
|
|
|
92 if item == prev_item:
|
|
|
93 # Keep iterating and storing values until a new value is encountered.
|
|
|
94 if block_valid:
|
|
|
95 for i, col in enumerate(cols):
|
|
|
96 if col >= 3:
|
|
|
97 prev_vals[i].append(fields[col].strip())
|
|
|
98 if len(set(prev_vals[i])) > 1:
|
|
|
99 block_valid = False
|
|
|
100 break
|
|
|
101
|
|
|
102 else:
|
|
|
103 """
|
|
|
104 When a new value is encountered, write the previous value and the
|
|
|
105 corresponding aggregate values into the output file. This works
|
|
|
106 due to the sort on group_col we've applied to the data above.
|
|
|
107 """
|
|
|
108 out_list = ['']*24
|
|
|
109 out_list[0] = str(prev_item)
|
|
|
110 out_list[1] = str(prev_vals[0][0])
|
|
|
111 out_list[2] = str(prev_vals[1][0])
|
|
|
112
|
|
|
113 for k, col in enumerate(cols):
|
|
|
114 if col >= 3 and col < 24:
|
|
|
115 if len(set(prev_vals[k])) == 1:
|
|
|
116 out_list[col] = prev_vals[k][0]
|
|
|
117 else:
|
|
|
118 break
|
|
|
119 while k < 23:
|
|
|
120 out_list[k+1] = 'n'
|
|
|
121 k += 1
|
|
|
122
|
|
|
123 j = 0
|
|
|
124 while True:
|
|
|
125 try:
|
|
|
126 out_list.append(str(prev_vals[23+j][0]))
|
|
|
127 j += 1
|
|
|
128 except:
|
|
|
129 break
|
|
|
130
|
|
|
131 if rank_bound == 0:
|
|
|
132 print >>fout, '\t'.join(out_list).strip()
|
|
|
133 else:
|
|
|
134 if ''.join(out_list[rank_bound:24]) != 'n'*( 24 - rank_bound ):
|
|
|
135 print >>fout, '\t'.join(out_list).strip()
|
|
|
136
|
|
|
137 block_valid = True
|
|
|
138 prev_item = item
|
|
|
139 prev_vals = []
|
|
|
140 for col in cols:
|
|
|
141 val_list = []
|
|
|
142 val_list.append(fields[col].strip())
|
|
|
143 prev_vals.append(val_list)
|
|
|
144
|
|
|
145 else:
|
|
|
146 # This only occurs once, right at the start of the iteration.
|
|
|
147 block_valid = True
|
|
|
148 prev_item = item #groupby item
|
|
|
149 for col in cols: #everyting else
|
|
|
150 val_list = []
|
|
|
151 val_list.append(fields[col].strip())
|
|
|
152 prev_vals.append(val_list)
|
|
|
153
|
|
|
154 except:
|
|
|
155 skipped_lines += 1
|
|
|
156 else:
|
|
|
157 skipped_lines += 1
|
|
|
158
|
|
|
159 # Handle the last grouped value
|
|
|
160 out_list = ['']*24
|
|
|
161 out_list[0] = str(prev_item)
|
|
|
162 out_list[1] = str(prev_vals[0][0])
|
|
|
163 out_list[2] = str(prev_vals[1][0])
|
|
|
164
|
|
|
165 for k, col in enumerate(cols):
|
|
|
166 if col >= 3 and col < 24:
|
|
|
167 if len(set(prev_vals[k])) == 1:
|
|
|
168 out_list[col] = prev_vals[k][0]
|
|
|
169 else:
|
|
|
170 break
|
|
|
171 while k < 23:
|
|
|
172 out_list[k+1] = 'n'
|
|
|
173 k += 1
|
|
|
174
|
|
|
175 j = 0
|
|
|
176 while True:
|
|
|
177 try:
|
|
|
178 out_list.append(str(prev_vals[23+j][0]))
|
|
|
179 j += 1
|
|
|
180 except:
|
|
|
181 break
|
|
|
182
|
|
|
183 if rank_bound == 0:
|
|
|
184 print >>fout, '\t'.join(out_list).strip()
|
|
|
185 else:
|
|
|
186 if ''.join(out_list[rank_bound:24]) != 'n'*( 24 - rank_bound ):
|
|
|
187 print >>fout, '\t'.join(out_list).strip()
|
|
|
188
|
|
|
189 if skipped_lines > 0:
|
|
|
190 print "Skipped %d invalid lines." % ( skipped_lines )
|
|
|
191
|
|
|
192 if __name__ == "__main__":
|
|
|
193 main() |
