-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy path20_motif_finder_exercise.py
executable file
·101 lines (85 loc) · 3.05 KB
/
20_motif_finder_exercise.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
#!/usr/bin/python3
#
# This file is part of Progesterone pipeline.
#
# Progesterone pipeline is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# Progesterone pipeline is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with Progesterone pipeline. If not, see <https://www.gnu.org/licenses/>.
#
from utils.utils import *
def read_pfm(jaspar_motifs_file, tf_name):
motif = None
with open(jaspar_motifs_file) as handle:
for m in motifs.parse(handle,"jaspar"):
if m.name==tf_name:
motif = m
break
return motif
#########################################
def read_or_download_sequence(chipseq_regions_dir, assembly, chromosome, tf, start, end):
seqfile = "{}/{}_{}_{}_{}_{}.txt".format(chipseq_regions_dir,tf, assembly, chromosome, start, end)
if (os.path.exists(seqfile) and os.path.getsize(seqfile)>0):
outf = open(seqfile, "r")
seq = outf.read()
outf.close()
else:
seq = ucsc_fragment_sequence(assembly, chromosome, start, end)
outf = open(seqfile, "w")
outf.write(seq.replace("\n",""))
outf.close()
return seq
#########################################
# noinspection PyUnreachableCode
def main():
tf = "ESR1"
if False:
assembly = "hg19"
chromosome = "4"
start = 174447651- 10000
end = 174447651
else:
assembly = "rn4"
# Hand2 at chr16:36324327-36326037
chromosome = "16"
start = 36326037
end = 36326037 + 50000
jaspar_motifs_file = "/storage/databases/jaspar/JASPAR2018_CORE_vertebrates_non-redundant_pfms_jaspar.txt"
for f in [ jaspar_motifs_file]:
if not os.path.exists(f):
print(f,"not found")
exit()
motif = read_pfm(jaspar_motifs_file, tf)
# add something so that the counts are not 0
pwm = motif.counts.normalize(pseudocounts=1)
pssm = pwm.log_odds()
# seq from UCSC
print("intial range:", start, end)
seq = ucsc_fragment_sequence(assembly,chromosome,start,end)
bpseq = Seq(seq,unambiguous_dna)
for position, score in pssm.search(bpseq, threshold=5.0):
if position>0: # motif is on the direct strand
offset = position
print("range:", start, end)
print("offset %d: score = %5.1f" % (offset, score))
matched_seq = bpseq[position:position+motif.length]
else: # motif is on the negative strand
offset = len(bpseq)+position
print("offset %d: score = %5.1f (on the compl strand)" % (offset, score))
matched_seq = bpseq[position:position+motif.length].reverse_complement()
print(motif.consensus," <--- consensus")
print(matched_seq.upper()," <--- motif found")
print(bpseq[position:position+motif.length]," <--- direct strand")
print()
#########################################
########################################
if __name__ == '__main__':
main()