Skip to content

Commit a38b284

Browse files
isaiahahisaiahah
authored
Merge pull request BioAnalyticResource#6 from isaiahah/hotspots
Create hotspot utils class
2 parents 70086bc + 8fdfd79 commit a38b284

1 file changed

Lines changed: 181 additions & 183 deletions

File tree

‎api/utils/hotspot_utils.py‎

Lines changed: 181 additions & 183 deletions
Original file line numberDiff line numberDiff line change
@@ -1,188 +1,186 @@
1-
# Files with probabilities of SNP significance based on shuffles
2-
# P_SNP_DIR = "/Users/isaiahhazelwood/Documents/UofT/Y3/BCB330/data/hotspots/"
3-
P_SNP_DIR = "/usr/src/hotspots"
4-
ARA_STRUCT_P_SNP = P_SNP_DIR + "/ara_struct_10_0.1_10000.tsv"
5-
POP_STRUCT_P_SNP = P_SNP_DIR + "/pop_struct_10_0.1_10000.tsv"
6-
ARA_SEQ_P_SNP = P_SNP_DIR + "/ara_seq_3_0.1_10000.tsv"
7-
POP_SEQ_P_SNP = P_SNP_DIR + "/pop_seq_3_0.1_10000.tsv"
8-
9-
# File with analyzed homologue pairs
10-
# HOMOLOGUE_DIR = "/Users/isaiahhazelwood/Documents/UofT/Y3/BCB330/data/homologue-info-pop3.0/"
11-
HOMOLOGUE_DIR = "/usr/src/pairs"
12-
ARA_POP_HOMOLOGUE = HOMOLOGUE_DIR + "/ara-pop3.0-all-valid.tsv"
13-
14-
15-
def verify_ara_pop_homologue(ara_id, pop_id):
16-
"""If both Arabidopsis and Poplar IDs are given, verifies they are a
17-
homologous pair. If only only one is given and the other is None
18-
find the matching homologous ID. Returns the IDs and their sequences
19-
20-
:param ara_id: str of TAIR10 gene ID, or None
21-
:param pop_id: str of Pop v3.0 gene ID, or None
22-
:returns: Tuple of IDs and sequences, or None if invalid IDs
23-
:rtype: Tuple[str, str, str, str] or None
24-
25-
"""
26-
if ara_id is None and pop_id is None: # Both invalid inputs
27-
return None
28-
with open(ARA_POP_HOMOLOGUE, "r") as f_ara_pop_homologue:
29-
for line in f_ara_pop_homologue:
30-
# Columns: araid, popid, araseq, popseq, rmsd
31-
cols = line.split("\t")
32-
if cols[0][4:-4].upper() == ara_id and cols[1][4:-4].upper() == pop_id: # Both ID match
33-
return (cols[0][4:-4].upper(), cols[1][4:-4].upper(), cols[2], cols[3])
34-
if cols[0][4:-4].upper() == ara_id and pop_id is None: # Ara ID match, fill Pop
35-
return (cols[0][4:-4].upper(), cols[1][4:-4].upper(), cols[2], cols[3])
36-
if cols[1][4:-4].upper() == pop_id and ara_id is None: # Pop ID match, fill Ara
37-
return (cols[0][4:-4].upper(), cols[1][4:-4].upper(), cols[2], cols[3])
38-
return None # No match
39-
40-
41-
def load_p_snp_data(id, spe, shuffle="struct"):
42-
"""Load the probability of SNP significance at each residue of the given
43-
protein from the cache file.
44-
45-
:param id: str of TAIR10 or Pop v3.0 gene ID
46-
:param spe: Either "ara" or "pop" based on id species
47-
:param shuffle: Either "struct" or "seq" for significance method.
48-
Defaults to "struct"
49-
:returns: List of significance scores at residue positions,
50-
or None if invalid ID.
51-
:rtype: List[float] or None
52-
53-
"""
54-
# Select appropriate file
55-
if spe == "ara":
56-
if shuffle == "struct":
57-
p_snps_file = ARA_STRUCT_P_SNP
58-
elif shuffle == "seq":
59-
p_snps_file = ARA_SEQ_P_SNP
60-
else:
1+
class hotspotUtils:
2+
# Files with probabilities of SNP significance based on shuffles
3+
P_SNP_DIR = "/usr/src/hotspots"
4+
ARA_STRUCT_P_SNP = P_SNP_DIR + "/ara_struct_10_0.1_10000.tsv"
5+
POP_STRUCT_P_SNP = P_SNP_DIR + "/pop_struct_10_0.1_10000.tsv"
6+
ARA_SEQ_P_SNP = P_SNP_DIR + "/ara_seq_3_0.1_10000.tsv"
7+
POP_SEQ_P_SNP = P_SNP_DIR + "/pop_seq_3_0.1_10000.tsv"
8+
9+
# File with analyzed homologue pairs
10+
HOMOLOGUE_DIR = "/usr/src/pairs"
11+
ARA_POP_HOMOLOGUE = HOMOLOGUE_DIR + "/ara-pop3.0-all-valid.tsv"
12+
13+
@staticmethod
14+
def verify_ara_pop_homologue(ara_id, pop_id):
15+
"""If both Arabidopsis and Poplar IDs are given, verifies they are a
16+
homologous pair. If only only one is given and the other is None
17+
find the matching homologous ID. Returns the IDs and their sequences
18+
19+
:param ara_id: str of TAIR10 gene ID, or None
20+
:param pop_id: str of Pop v3.0 gene ID, or None
21+
:returns: Tuple of IDs and sequences, or None if invalid IDs
22+
:rtype: Tuple[str, str, str, str] or None
23+
24+
"""
25+
if ara_id is None and pop_id is None: # Both invalid inputs
6126
return None
62-
elif spe == "pop":
63-
if shuffle == "struct":
64-
p_snps_file = POP_STRUCT_P_SNP
65-
elif shuffle == "seq":
66-
p_snps_file = POP_SEQ_P_SNP
27+
with open(hotspotUtils.ARA_POP_HOMOLOGUE, "r") as f_ara_pop_homologue:
28+
for line in f_ara_pop_homologue:
29+
# Columns: araid, popid, araseq, popseq, rmsd
30+
cols = line.split("\t")
31+
if cols[0][4:-4].upper() == ara_id and cols[1][4:-4].upper() == pop_id: # Both ID match
32+
return (cols[0][4:-4].upper(), cols[1][4:-4].upper(), cols[2], cols[3])
33+
if cols[0][4:-4].upper() == ara_id and pop_id is None: # Ara ID match, fill Pop
34+
return (cols[0][4:-4].upper(), cols[1][4:-4].upper(), cols[2], cols[3])
35+
if cols[1][4:-4].upper() == pop_id and ara_id is None: # Pop ID match, fill Ara
36+
return (cols[0][4:-4].upper(), cols[1][4:-4].upper(), cols[2], cols[3])
37+
return None # No match
38+
39+
@staticmethod
40+
def load_p_snp_data(id, spe, shuffle="struct"):
41+
"""Load the probability of SNP significance at each residue of the given
42+
protein from the cache file.
43+
44+
:param id: str of TAIR10 or Pop v3.0 gene ID
45+
:param spe: Either "ara" or "pop" based on id species
46+
:param shuffle: Either "struct" or "seq" for significance method.
47+
Defaults to "struct"
48+
:returns: List of significance scores at residue positions,
49+
or None if invalid ID.
50+
:rtype: List[float] or None
51+
52+
"""
53+
# Select appropriate file
54+
if spe == "ara":
55+
if shuffle == "struct":
56+
p_snps_file = hotspotUtils.ARA_STRUCT_P_SNP
57+
elif shuffle == "seq":
58+
p_snps_file = hotspotUtils.ARA_SEQ_P_SNP
59+
else:
60+
return None
61+
elif spe == "pop":
62+
if shuffle == "struct":
63+
p_snps_file = hotspotUtils.POP_STRUCT_P_SNP
64+
elif shuffle == "seq":
65+
p_snps_file = hotspotUtils.POP_SEQ_P_SNP
66+
else:
67+
return None
6768
else:
6869
return None
69-
else:
70+
71+
# Load data from file
72+
with open(p_snps_file, "r") as f_p_snps:
73+
for line in f_p_snps:
74+
if line.upper().startswith(id):
75+
return [float(p) for p in line.split("\t")[1].split(",")]
7076
return None
7177

72-
# Load data from file
73-
with open(p_snps_file, "r") as f_p_snps:
74-
for line in f_p_snps:
75-
if line.upper().startswith(id):
76-
return [float(p) for p in line.split("\t")[1].split(",")]
77-
return None
78-
79-
80-
def mark_significant(null_probs, p):
81-
"""Mark residues with p-significant SNPs.
82-
83-
:param null_probs: List of significance scores at residue positions.
84-
:param p: p-value for significance
85-
:returns: Boolean list of residues significant at given p-value.
86-
:rtype: List[bool]
87-
88-
"""
89-
return [(prob >= p) for prob in null_probs]
90-
91-
92-
def match_residues(aln):
93-
"""For each index in the first protein which aligns with a residue in the
94-
second protein (not a gap), provide the aligned index in the second protein.
95-
96-
:param aln: Tuple of Arabidopsis and Poplar protein sequences
97-
:returns: Dict from index in first protein to index in second protein
98-
:rtype: Dict[int, int]
99-
"""
100-
matchings = {}
101-
curr_prot1 = 1 # Current index in first protein
102-
curr_prot2 = 1 # Current index in second protein
103-
# Iterate over all positions in the proteins
104-
for i in range(len(aln[0])):
105-
# If both not gaps, match the indices
106-
if aln[0][i] != "-" and aln[1][i] != "-":
107-
matchings[curr_prot1] = curr_prot2
108-
# If the position in the first protein is not a gap, increment index
109-
if aln[0][i] != "-":
110-
curr_prot1 += 1
111-
# If the position in the second protein is not a gap, increment index
112-
if aln[1][i] != "-":
113-
curr_prot2 += 1
114-
return matchings
115-
116-
117-
def significant_in_both(sig1, sig2, aln_matching):
118-
"""Mark a residue as significant in both if it aligns with a residue
119-
in the other protein and those residues are both marked significant.
120-
121-
:param sig1: Boolean list of significant residues in protein 1.
122-
:param sig2: Boolean list of significant residues in protein 2.
123-
:param aln_matching: Dictionary of aligned indices from protein 1 to 2.
124-
:returns: Two Boolean lists of significant residues in both proteins
125-
:rtype: List[bool], List[bool]
126-
"""
127-
# Create significance array for overlap, initialize to false
128-
both_sig1 = [False] * len(sig1)
129-
both_sig2 = [False] * len(sig2)
130-
# For each aligned index in protein 1
131-
for i in aln_matching:
132-
# If that reside and the counterpart in protein 2 are both significant,
133-
# mark those residues as significant in both.
134-
# Dictionary stores 1-indexed positions, list is 0-indexed
135-
if sig1[i - 1] and sig2[aln_matching[i] - 1]:
136-
both_sig1[i - 1] = True
137-
both_sig2[aln_matching[i] - 1] = True
138-
return (both_sig1, both_sig2)
139-
140-
141-
def get_sig_index(sig):
142-
"""Return the 1-based indexes marked as significant (true in input list).
143-
144-
:param sig: Boolean list of significant residues.
145-
:returns: List of 1-indexed positions marked significant.
146-
:rtype: List[int]
147-
"""
148-
return [(i + 1) for i in range(len(sig)) if sig[i]]
149-
150-
151-
# Find hotspot clusters
152-
def cluster_components(hotspots, neighbours):
153-
"""Determine clusters of hotspots residues. Clusters are connected
154-
components in the neighbourhood graph, found using DFS.
155-
156-
BFS algorithm: Initialize a frontier of positions to explore.
157-
Add the starting residue to the frontier. While the frontier is not empty,
158-
remove the last item from the frontier and add its significant non-explored
159-
residues to the frontier.
160-
161-
:param hotspots: List of hotspot indices
162-
:param neighbours: Adjacency list of neighbourhood graph, mapping each index
163-
to a tuple of neighbouring indices
164-
:returns: List of clusters, each a list of indices
165-
:rtype: List[List[int]]
166-
"""
167-
clusters = []
168-
# Track explored residues. Set non-hotspot residues to explored, so the
169-
# algorithm will not visit them.
170-
171-
explored = [not hotspot for hotspot in hotspots]
172-
residue_frontier = []
173-
# explored and residue_fronter is 0-indexed
174-
# neighbors and clusters is 1-indexed.
175-
for i in range(len(hotspots)): # O-indexed i
176-
if not explored[i]:
177-
residue_frontier.append(i) # Expand: add to frontier, set explored
178-
explored[i] = True
179-
curr_cluster = []
180-
while len(residue_frontier) > 0:
181-
curr_res = residue_frontier.pop()
182-
curr_cluster.append(curr_res + 1) # Add to cluster
183-
for neighbour in neighbours[i + 1]: # Push neighbours
184-
if not explored[neighbour - 1]:
185-
residue_frontier.append(neighbour - 1)
186-
explored[neighbour - 1] = True
187-
clusters.append(curr_cluster) # Save component on DFS finish
188-
return clusters
78+
@staticmethod
79+
def mark_significant(null_probs, p):
80+
"""Mark residues with p-significant SNPs.
81+
82+
:param null_probs: List of significance scores at residue positions.
83+
:param p: p-value for significance
84+
:returns: Boolean list of residues significant at given p-value.
85+
:rtype: List[bool]
86+
87+
"""
88+
return [(prob >= p) for prob in null_probs]
89+
90+
@staticmethod
91+
def match_residues(aln):
92+
"""For each index in the first protein which aligns with a residue in the
93+
second protein (not a gap), provide the aligned index in the second protein.
94+
95+
:param aln: Tuple of Arabidopsis and Poplar protein sequences
96+
:returns: Dict from index in first protein to index in second protein
97+
:rtype: Dict[int, int]
98+
"""
99+
matchings = {}
100+
curr_prot1 = 1 # Current index in first protein
101+
curr_prot2 = 1 # Current index in second protein
102+
# Iterate over all positions in the proteins
103+
for i in range(len(aln[0])):
104+
# If both not gaps, match the indices
105+
if aln[0][i] != "-" and aln[1][i] != "-":
106+
matchings[curr_prot1] = curr_prot2
107+
# If the position in the first protein is not a gap, increment index
108+
if aln[0][i] != "-":
109+
curr_prot1 += 1
110+
# If the position in the second protein is not a gap, increment index
111+
if aln[1][i] != "-":
112+
curr_prot2 += 1
113+
return matchings
114+
115+
@staticmethod
116+
def significant_in_both(sig1, sig2, aln_matching):
117+
"""Mark a residue as significant in both if it aligns with a residue
118+
in the other protein and those residues are both marked significant.
119+
120+
:param sig1: Boolean list of significant residues in protein 1.
121+
:param sig2: Boolean list of significant residues in protein 2.
122+
:param aln_matching: Dictionary of aligned indices from protein 1 to 2.
123+
:returns: Two Boolean lists of significant residues in both proteins
124+
:rtype: List[bool], List[bool]
125+
"""
126+
# Create significance array for overlap, initialize to false
127+
both_sig1 = [False] * len(sig1)
128+
both_sig2 = [False] * len(sig2)
129+
# For each aligned index in protein 1
130+
for i in aln_matching:
131+
# If that reside and the counterpart in protein 2 are both significant,
132+
# mark those residues as significant in both.
133+
# Dictionary stores 1-indexed positions, list is 0-indexed
134+
if sig1[i - 1] and sig2[aln_matching[i] - 1]:
135+
both_sig1[i - 1] = True
136+
both_sig2[aln_matching[i] - 1] = True
137+
return (both_sig1, both_sig2)
138+
139+
@staticmethod
140+
def get_sig_index(sig):
141+
"""Return the 1-based indexes marked as significant (true in input list).
142+
143+
:param sig: Boolean list of significant residues.
144+
:returns: List of 1-indexed positions marked significant.
145+
:rtype: List[int]
146+
"""
147+
return [(i + 1) for i in range(len(sig)) if sig[i]]
148+
149+
@staticmethod
150+
def cluster_components(hotspots, neighbours):
151+
"""Determine clusters of hotspots residues. Clusters are connected
152+
components in the neighbourhood graph, found using DFS.
153+
154+
BFS algorithm: Initialize a frontier of positions to explore.
155+
Add the starting residue to the frontier. While the frontier is not empty,
156+
remove the last item from the frontier and add its significant non-explored
157+
residues to the frontier.
158+
159+
:param hotspots: List of hotspot indices
160+
:param neighbours: Adjacency list of neighbourhood graph, mapping each index
161+
to a tuple of neighbouring indices
162+
:returns: List of clusters, each a list of indices
163+
:rtype: List[List[int]]
164+
"""
165+
clusters = []
166+
# Track explored residues. Set non-hotspot residues to explored, so the
167+
# algorithm will not visit them.
168+
169+
explored = [not hotspot for hotspot in hotspots]
170+
residue_frontier = []
171+
# explored and residue_fronter is 0-indexed
172+
# neighbors and clusters is 1-indexed.
173+
for i in range(len(hotspots)): # O-indexed i
174+
if not explored[i]:
175+
residue_frontier.append(i) # Expand: add to frontier, set explored
176+
explored[i] = True
177+
curr_cluster = []
178+
while len(residue_frontier) > 0:
179+
curr_res = residue_frontier.pop()
180+
curr_cluster.append(curr_res + 1) # Add to cluster
181+
for neighbour in neighbours[i + 1]: # Push neighbours
182+
if not explored[neighbour - 1]:
183+
residue_frontier.append(neighbour - 1)
184+
explored[neighbour - 1] = True
185+
clusters.append(curr_cluster) # Save component on DFS finish
186+
return clusters

0 commit comments

Comments
 (0)