-
Notifications
You must be signed in to change notification settings - Fork 9
/
12av_removeMAFoverlaps_chr2.py
161 lines (135 loc) · 5.29 KB
/
12av_removeMAFoverlaps_chr2.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
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
#-----------------------------------------------------------+
# |
# 12av_removeMAFoverlaps.py - script to take longest alignment from the overlapping alignments |
# |
#-----------------------------------------------------------+
# |
# AUTHOR: Vikas Gupta |
# CONTACT: [email protected] |
# STARTED: 09/06/2013 |
# UPDATED: 09/06/2013 |
# |
# DESCRIPTION: |
# Short script to convert and copy the wheat BACs |
# Run this in the parent dir that the HEX* dirs exist |
# |
# LICENSE: |
# GNU General Public License, Version 3 |
# http://www.gnu.org/licenses/gpl.html |
# |
#-----------------------------------------------------------+
# Example:
# python ~/script/python/100b_fasta2flat.py -i 02_Stegodyphous_cdna.refined.fa.orf.tr_longest_frame
### import modules
import os,sys,getopt, re, classMAF
### global variables
global ifile, bed
global chro
#chro = ['chr1', 'chr2' , 'chr3', 'chr4', 'chr5', 'chr6']
chro = ['chr2']
### make a logfile
import datetime
now = datetime.datetime.now()
o = open(str(now.strftime("%Y-%m-%d_%H%M."))+'logfile','w')
### write logfile
def logfile(infile):
o.write("Program used: \t\t%s" % "100b_fasta2flat.py"+'\n')
o.write("Program was run at: \t%s" % str(now.strftime("%Y-%m-%d_%H%M"))+'\n')
o.write("Infile used: \t\t%s" % infile+'\n')
def help():
print '''
python 100b_fasta2flat.py -i <ifile>
'''
sys.exit(2)
### main argument to
def options(argv):
global ifile, bed
ifile = ''
bed = False
try:
opts, args = getopt.getopt(argv,"hi:b",["ifile=","bed="])
except getopt.GetoptError:
help()
for opt, arg in opts:
if opt == '-h':
help()
elif opt in ("-i", "--ifile"):
ifile = arg
elif opt in ("-b", "--bed"):
bed = True
logfile(ifile)
def parseMAF():
o = open(ifile+'.nonRed','w')
if bed == True:
b = open(ifile+'.bed','w')
objs = ''
new_black = True
ref_align = True
obj_count=1
obj = classMAF.MAF(obj_count)
objs_list = [obj]
for chrom in chro:
for line in open(ifile, 'r'):
line = line.strip()
if line == '':
obj_count += 1
ref_align = True
obj = classMAF.MAF(obj_count)
objs_list.append(obj)
else:
if not line.startswith('s '):
obj.addData(line)
else:
if ref_align == True:
obj.addRefAlign(line)
ref_align = False
else:
obj.addtargetAlign(line)
chr_obj_pos = {}
for ob in objs_list[:-1]:
if ob.RefChro() == chrom:
chr_obj_pos[int(ob.RefStart())] = ob
chr_obj = []
for key in sorted(chr_obj_pos):
chr_obj.append(chr_obj_pos[key])
hash_coords={}
hash_coords_obj={}
obj_list = {}
obj_count = 0
for obj in chr_obj:
obj_count += 1
if obj.RefChro() == chrom:
for i in range(int(obj.RefStart()), int(obj.RefEnd())+1, 1):
if i in hash_coords:
if hash_coords_obj[i].RefSize() > obj.RefSize():
if obj in obj_list:
del obj_list[obj]
break
else:
hash_coords[i] = obj.RefSize()
if hash_coords_obj[i] in obj_list:
del obj_list[hash_coords_obj[i]]
if obj not in obj_list:
obj_list[obj]=''
hash_coords_obj[i] = obj
else:
hash_coords[i] = obj.RefSize()
if obj not in obj_list:
obj_list[obj] = ''
hash_coords_obj[i] = obj
print len(obj_list)
print chrom, '{:9,.0f}'.format(i)
print obj.TargetChro(), obj.TargetStart()
for ob in obj_list:
o.write(ob.Data()+'\n')
if bed == True:
b.write(ob.RefChro()+'\t'+str(ob.RefStart())+'\t'+str(ob.RefEnd())+'\t'+ob.TargetChro()+'_'+str(ob.TargetStart())+'_'+str(ob.TargetEnd())+'\t'+str(ob.identity())+'\t'+str(ob.RefStrand())+'\n')
o.close()
if bed == True:
b.close()
if __name__ == "__main__":
options(sys.argv[1:])
### parseMAF
parseMAF()
### close the logfile
o.close()