forked from sreeramkannan/Shannon
-
Notifications
You must be signed in to change notification settings - Fork 0
/
filter_FP_2s.py
208 lines (177 loc) · 7.5 KB
/
filter_FP_2s.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
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
import os, sys
def run_cmd(s1):
print(s1);
os.system(s1) # + '>> temp_out.txt')
def write_filtered_tr(depth_file, in_tr_file, out_tr_file, log_file):
import sys, pdb
tr_hits = {}
THRESH = 0.9
for a in open(depth_file): # as depth_file:
#a=depth_file.readlines()
fields = a.strip().split()
tr_hits[fields[0]] = tr_hits.get(fields[0],0) +1;
tr_len = {}; tr_name = ''
with open(out_tr_file,'w') as tr_file, open(log_file,'w') as lf: #output transcripts
for line in open(in_tr_file): #in_tr_file
fields = line.strip().split()
if fields[0][0] == '>': tr_name = fields[0][1:]; continue
tr_len[tr_name]=len(fields[0]); trh = tr_hits.get(tr_name,0)
lf.write(tr_name + '\t' + str(trh) + '\t' + str(len(fields[0])) + '\n')
if tr_hits.get(tr_name,0) >= len(fields[0]) * THRESH:
tr_file.write('>' + tr_name + '\n')
tr_file.write(fields[0] + '\n')
def convert_vector(l1):
nl = [];
for i in l1:
if nl and nl[-1][1]==(i-1):
nl[-1][1]=i; continue
else:
nl.append([i,i])
return nl
def break_seqs(tr_l,tr_r,l1):
K = 25
insert = 200
L = 100
MIN_SEQ = 200 #Only write sequences longer than MIN_SEQ
#l1 = len(seq)
seqs = []
if len(tr_l)==0 or len(tr_r)==0: return seqs
#Add end buffers
if tr_l[0][0] < K: tr_l[0][0] =1
if tr_l[-1][1] > l1-insert: tr_l[-1][1] = l1
if tr_r[0][0] < insert: tr_r[0][0]=1
if tr_r[-1][1] < K : tr_r[-1][1] = l1
l_max =len(tr_l)-1; r_max = len(tr_r)-1
l_ptr = 0; r_ptr =0
while(1):
print(str(l_ptr) + ',' + str(r_ptr))
if l_ptr > l_max or r_ptr > r_max: break
curr_l = tr_l[l_ptr]; curr_r = tr_r[r_ptr]
print(curr_l)
print(curr_r)
l1 = max(curr_l[0],curr_r[0])
u1 = min(curr_l[1],curr_r[1])
if l1+ MIN_SEQ<=u1: seqs.append([l1,u1])
if curr_l[1] <= u1:
l_ptr +=1
else:
tr_l[l_ptr][0] = l1
if curr_r[1] <= u1:
r_ptr +=1
else:
tr_r[r_ptr][0] = l1
return seqs
def break_seqs_new(tr_l,tr_r,l1):
K = 25
insert = 200
L = 100
MIN_SEQ = 200 #Only write sequences longer than MIN_SEQ
#l1 = len(seq)
seqs = []
if len(tr_l)==0 or len(tr_r)==0: return seqs
#Add end buffers
if tr_l[0][0] < K: tr_l[0][0] =1
if tr_l[-1][1] > l1-insert: tr_l[-1][1] = l1
if tr_r[0][0] < insert: tr_r[0][0]=1
if tr_r[-1][1] < K : tr_r[-1][1] = l1
l_max =len(tr_l)-1; r_max = len(tr_r)-1
l_ptr = 0; r_ptr =0
while(1):
print(str(l_ptr) + ',' + str(r_ptr))
if l_ptr > l_max or r_ptr > r_max: break
curr_l = tr_l[l_ptr]; curr_r = tr_r[r_ptr]
print(curr_l)
print(curr_r)
l1 = max(curr_l[0],curr_r[0])
u1 = min(curr_l[1],curr_r[1])
if l1<=u1:
#There is intersection
new_l1 = min(curr_l[0],curr_r[0])
new_u1 = max(curr_l[1],curr_r[1])
seqs.append([new_l1,new_u1])
l_ptr +=1; r_ptr+=1
else:
#No intersection
if curr_l[0] < l1:
l_ptr +=1
if curr_r[0] < u1:
r_ptr +=1
return seqs
def read_in(depth_file):
tr_dict = {}
curr_name = ''; curr_list = []
for line in open(depth_file): # as left_depth_file:
fields = line.strip().split()
if fields[0] == curr_name:
curr_list.append(int(fields[1]))
else:
nl = convert_vector(curr_list)
if curr_name: tr_dict[curr_name]=nl
curr_name = fields[0]
curr_list = [int(fields[1])]
nl = convert_vector(curr_list)
if curr_name: tr_dict[curr_name]=nl
return tr_dict
def write_new_tr(ld_file, rd_file, in_tr_file, out_tr_file, log_file):
import sys, pdb
MIN_SEQ = 200 #Only write sequences longer than MIN_SEQ
tr_left = read_in(ld_file)
tr_right = read_in(rd_file)
tr_len = {}; tr_name = ''
with open(out_tr_file,'w') as tr_file, open(log_file,'w') as lf: #output transcripts
for line in open(in_tr_file): #in_tr_file
fields = line.strip().split()
if fields[0][0] == '>': tr_name = fields[0][1:]; continue
tr_len[tr_name]=len(fields[0]);
bs = break_seqs_new(tr_left.get(tr_name,[]),tr_right.get(tr_name,[]),tr_len.get(tr_name,[]))
seq = fields[0]
ctr = 0
for (l,u) in bs:
if l + MIN_SEQ <= u:
aug = '' if (l==1 and u==len(seq)) else ('_frag_' + str(ctr))
tr_file.write('>' + tr_name + aug + '\n')
tr_file.write(seq[l-1:u] + '\n')
ctr +=1
lf.write(tr_name + '\t' + str(bs) + '\n')
def filter_FP(rec_fasta, read_1, read_2, out_dir, flags = '-f --ff'):
#take rec_fasta, read_1, read_2 files of type.
#Output in out_fasta. Temp dir = out_dir.
#Output fasta is in out_dir/reconstructed.fasta
#Flags is '-f --ff' for fasta and forward-forward
#Flags is '-q --fr' for fastq and forward-reverse
hisat_dir = '' #include dir/ if specifying directory
#Build hisat file
run_cmd(hisat_dir + 'hisat-build ' + rec_fasta + ' ' + out_dir + '/rec.hisat')
#Align
run_cmd(hisat_dir + 'hisat --no-spliced-alignment --no-discordant '+ flags + ' -x ' + out_dir + '/rec.hisat -1 ' + read_1 + ' -2 ' + read_2 + ' -S ' + out_dir + '/rec.sam' )
#Process SAM / BAM file to get depth information
run_cmd('samtools view -bS -f 0x2 ' + out_dir + '/rec.sam > ' + out_dir + '/rec.bam')
#run_cmd('samtools sort '+ out_dir + '/rec.bam ' + out_dir + '/rec_sort')
#run_cmd('samtools depth ' + out_dir + '/rec_sort.bam > ' + out_dir + '/rec.depth')
run_cmd('samtools view -b -f 64 ' + out_dir + '/rec.bam > ' + out_dir + '/c1.bam')
run_cmd('samtools view -b -f 128 ' + out_dir + '/rec.bam > ' + out_dir + '/c2.bam')
run_cmd('samtools view -b -F 16 ' + out_dir + '/c1.bam > ' + out_dir + '/c1f.bam')
run_cmd('samtools view -b -f 16 ' + out_dir + '/c1.bam > ' + out_dir + '/c1r.bam')
run_cmd('samtools view -b -f 32 ' + out_dir + '/c2.bam > ' + out_dir + '/c2r.bam')
run_cmd('samtools view -b -F 32 ' + out_dir + '/c2.bam > ' + out_dir + '/c2f.bam')
run_cmd('samtools merge ' + out_dir + '/l.bam ' + out_dir + '/c1f.bam ' +out_dir + '/c2r.bam ' )
run_cmd('samtools merge ' + out_dir + '/r.bam ' + out_dir + '/c2f.bam ' +out_dir + '/c1r.bam ' )
run_cmd('samtools sort '+ out_dir + '/l.bam ' + out_dir + '/l_sort')
run_cmd('samtools sort '+ out_dir + '/r.bam ' + out_dir + '/r_sort')
run_cmd('samtools depth ' + out_dir + '/l_sort.bam > ' + out_dir + '/l.depth')
run_cmd('samtools depth ' + out_dir + '/r_sort.bam > ' + out_dir + '/r.depth')
#Use the BAM file along with our tool to get new fasta file
ld_file = out_dir + '/l.depth'; rd_file = out_dir + '/r.depth'
in_tr_file = rec_fasta
out_tr_file = out_dir + '/rec.fasta';
log_file = out_dir + '/rec.log'
write_new_tr(ld_file, rd_file, in_tr_file, out_tr_file, log_file)
run_cmd('cp ' + rec_fasta + ' ' + out_dir +'/reconstructed_org.fasta')
run_cmd('mv ' + out_tr_file + ' ' + out_dir + '/reconstructed.fasta')
if __name__ == '__main__':
rec_fasta = sys.argv[1]
read_1 = sys.argv[2]
read_2 = sys.argv[3]
out_dir = sys.argv[4]
#flags = ' '.join(sys.argv[5:])
filter_FP(rec_fasta, read_1, read_2, out_dir)