-
Notifications
You must be signed in to change notification settings - Fork 0
/
unitig_filter.py
executable file
·181 lines (143 loc) · 6.64 KB
/
unitig_filter.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
#!/usr/bin/env python
import numpy as np
import sys
from Bio import SeqIO
paf_file = sys.argv[1]
unitigs = sys.argv[2]
output_stats = sys.argv[3]
output_file = sys.argv[4]
unitigs_corrected = open(output_file, 'w')
stats_file = open(output_stats, 'a')
# open up sequence map for illumina
illumina_dict = SeqIO.index_db(unitigs + ".idx", unitigs, "fasta")
def find_sequence_illumina(seqid, start, end):
return illumina_dict[seqid].seq[start:end + 1]
def write_sequence(f, seq):
for i in range(0, len(seq), 60):
f.write(str(seq[i:i + 60]) + "\n")
############ read input paf file #######################################################################################
mapping_freq = {}
nano_lengths = {}
with open(paf_file, 'r') as input_file:
current_id = ''
current_length = 0
current_nano_ids = []
current_profile = []
block_counter = 0
for line in input_file:
split_line = line.rstrip().split('\t')
illu_id = split_line[0]
illu_length = int(split_line[1])
illuminaStart = int(split_line[2])
illuminaEnd = int(split_line[3]) - 1
nano_id = split_line[5]
nano_length = split_line[6]
if nano_id not in nano_lengths:
nano_lengths[nano_id] = int(nano_length)
if current_id == illu_id: # still in the same block
if nano_id not in current_nano_ids:
block_counter += 1
current_nano_ids.append(nano_id)
for pos in range(illuminaStart, illuminaEnd + 1):
current_profile[pos] += 1
else: # new block
if current_id != '': # close previous block, avoid beginning
mapping_freq[current_id] = [block_counter, current_length, max(current_profile)]
# start new block
current_profile = [0] * illu_length
for pos in range(illuminaStart, illuminaEnd + 1):
current_profile[pos] += 1
block_counter = 1
current_id = illu_id
current_length = illu_length
del current_nano_ids[:]
current_nano_ids.append(nano_id)
mapping_freq[current_id] = [block_counter, current_length, max(current_profile)] # last block
#print("finished reading input file")
############# define filter cut-off ####################################################################################
def filter_cutoff(mapping_dict):
blocks_maxcov = [x[2] for x in list(mapping_dict.values())]
# define cut-off
q1 = np.percentile(blocks_maxcov, 25)
q3 = np.percentile(blocks_maxcov, 75)
iqr = q3 - q1
upper_whisker = q3 + 1.5 * iqr
stats_file.write(">>> unitig filter \n")
stats_file.write("upper_outlier: " + str(upper_whisker)+'\n')
stats_file.write("Q3: " + str(q3)+'\n')
return upper_whisker, blocks_maxcov, q3
cutoff, blocks_maxcov, q3 = filter_cutoff(mapping_freq)
########### evaluate and cut outliers #################################################################################
def cut_peaks(unitig_id, coverage_profile, file):
subrange = (0, 0)
below_cutoff = False
unitig_fragments = []
for i, cov in enumerate(coverage_profile):
if cov <= q3: # set q3 as cutoff
if not below_cutoff: # start new subrange below cutoff
subrange = (i, i)
else: # prolonging subrange
subrange = (subrange[0], i)
below_cutoff = True
else:
if below_cutoff and subrange[1] - subrange[0] + 1 >= 500:
file.write('>' + unitig_id + '_' + str(len(unitig_fragments)) + ' ' + str(subrange[1] - subrange[0] + 1)
+ ' ' + str(subrange[0]) + ' ' + str(subrange[1]) + '\n')
write_sequence(file, find_sequence_illumina(unitig_id, subrange[0], subrange[1]))
unitig_fragments.append(subrange)
subrange = (0, 0)
below_cutoff = False
if below_cutoff and subrange[1] - subrange[0] + 1 >= 500:
file.write('>' + unitig_id + '_' + str(len(unitig_fragments)) + ' ' + str(subrange[1] - subrange[0] + 1)
+ ' ' + str(subrange[0]) + ' ' + str(subrange[1]) + '\n')
write_sequence(file, find_sequence_illumina(unitig_id, subrange[0], subrange[1]))
unitig_fragments.append(subrange)
return unitig_fragments
########### detect outlier and write new unitig file ###########################################################
with open(paf_file, 'r') as input_file:
current_id = ''
current_profile = []
all_counter = 0
outlier_counter = 0
rescued_outlier = 0
total_rescued_length=0
for line in input_file:
split_line = line.rstrip().split('\t')
illu_id = split_line[0]
illu_length = int(split_line[1])
illuminaStart = int(split_line[2])
illuminaEnd = int(split_line[3]) - 1
if current_id == illu_id: # still in the same block
if mapping_freq[illu_id][2] > cutoff: # outlier unitig block
for pos in range(illuminaStart, illuminaEnd + 1):
current_profile[pos] += 1
else: # normal unitig block
continue
else: # new block
all_counter += 1
if current_id != '' and mapping_freq[current_id][2] > cutoff: # close previous block if outlier
fragments = cut_peaks(current_id, current_profile, unitigs_corrected)
if len(fragments) > 0:
rescued_outlier += 1
for entry in fragments:
total_rescued_length+=(entry[1] - entry[0] + 1)
outlier_counter += 1
if mapping_freq[illu_id][2] > cutoff: # new outlier block
current_profile = [0] * illu_length
for pos in range(illuminaStart, illuminaEnd + 1):
current_profile[pos] += 1
current_id = illu_id
else: # new normal unitig
unitigs_corrected.write('>' + illumina_dict[illu_id].description + '\n')
write_sequence(unitigs_corrected, illumina_dict[illu_id].seq)
current_id = illu_id
if current_id != '' and mapping_freq[current_id][2] > cutoff: # close last block if outlier
fragments = cut_peaks(current_id, current_profile, unitigs_corrected)
if len(fragments) > 0:
rescued_outlier += 1
outlier_counter += 1
stats_file.write("#all unitigs: " + str(all_counter)+'\n')
stats_file.write("#outliers: " + str(outlier_counter)+'\n')
stats_file.write("#rescued outliers: " + str(rescued_outlier)+'\n')
unitigs_corrected.close()
stats_file.close()