-
Notifications
You must be signed in to change notification settings - Fork 42
/
fasta_mask_regions.py
executable file
·110 lines (83 loc) · 2.41 KB
/
fasta_mask_regions.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
#!/usr/bin/env python3
"""Mask regions from fasta file given in bedfile
Usage:
<program> input_fasta input_bed flank_size output_fasta
beddfile format:
sequence_name<TAB>from<TAB>to<NEWLINE>
"""
# Modules
from collections import defaultdict
import gzip
import sys
# Classes
class Fasta(object):
"""Fasta object with name and sequence
"""
def __init__(self, name, sequence):
self.name = name
self.sequence = sequence
def write_to_file(self, handle):
handle.write(">" + self.name + "\n")
handle.write(self.sequence + "\n")
def __repr__(self):
return self.name + " " + self.sequence[:31]
# Functions
def myopen(_file, mode="rt"):
if _file.endswith(".gz"):
return gzip.open(_file, mode=mode)
else:
return open(_file, mode=mode)
def fasta_iterator(input_file):
"""Takes a fasta file input_file and returns a fasta iterator
"""
with myopen(input_file) as f:
sequence = []
name = ""
begun = False
for line in f:
line = line.strip()
if line.startswith(">"):
if begun:
yield Fasta(name, "".join(sequence))
name = line[1:]
sequence = ""
begun = True
else:
sequence += line
if name != "":
yield Fasta(name, "".join(sequence))
def overlap(r1, r2):
"""Return True if ranges overlap, even by only one position
"""
return (r1[0] <= r2[1]) and (r1[1] >= r2[0])
# Parsing user input
try:
input_fasta = sys.argv[1]
input_bed = sys.argv[2]
flank_size = int(sys.argv[3])
output_fasta = sys.argv[4]
except:
print(__doc__)
sys.exit(1)
# Get regions
regions = defaultdict(set)
count = 0
with open(input_bed) as infile:
for line in infile:
count += 1
l = line.strip().split("\t")
chrom = l[0]
start = int(l[1])
stop = int(l[2])
regions[chrom].update(range(start-flank_size, stop+flank_size+1))
# Read, mask, write fasta sequences
sequences = fasta_iterator(input_fasta)
with open(output_fasta, "wt") as outfile:
for s in sequences:
print(s, s.name in regions)
print(f"masking {len(regions[s.name])} positions")
seq = list(s.sequence)
for p in regions[s.name]:
seq[p] = "N"
s.sequence = "".join(seq)
s.write_to_file(outfile)