1515
1616
1717def verify_ara_pop_homologue (ara_id , pop_id ):
18- """If both Arabidopsis and Poplar IDs are given, verifies they are a
18+ """If both Arabidopsis and Poplar IDs are given, verifies they are a
1919 homologous pair. If only only one is given and the other is None
2020 find the matching homologous ID. Returns the IDs and their sequences
2121
@@ -25,36 +25,30 @@ def verify_ara_pop_homologue(ara_id, pop_id):
2525 :rtype: Tuple[str, str, str, str] or None
2626
2727 """
28- if ara_id is None and pop_id is None : # Both invalid inputs
28+ if ara_id is None and pop_id is None : # Both invalid inputs
2929 return None
30- with open (ARA_POP_HOMOLOGUE , 'r' ) as f_ara_pop_homologue :
30+ with open (ARA_POP_HOMOLOGUE , "r" ) as f_ara_pop_homologue :
3131 for line in f_ara_pop_homologue :
3232 # Columns: araid, popid, araseq, popseq, rmsd
33- cols = line .split ('\t ' )
34- if (cols [0 ][4 :- 4 ].upper () == ara_id and
35- cols [1 ][4 :- 4 ].upper () == pop_id ): # Both ID match
36- return (cols [0 ][4 :- 4 ].upper (), cols [1 ][4 :- 4 ].upper (),
37- cols [2 ], cols [3 ])
38- if (cols [0 ][4 :- 4 ].upper () == ara_id and
39- pop_id is None ): # Ara ID match, fill Pop
40- return (cols [0 ][4 :- 4 ].upper (), cols [1 ][4 :- 4 ].upper (),
41- cols [2 ], cols [3 ])
42- if (cols [1 ][4 :- 4 ].upper () == pop_id and
43- ara_id is None ): # Pop ID match, fill Ara
44- return (cols [0 ][4 :- 4 ].upper (), cols [1 ][4 :- 4 ].upper (),
45- cols [2 ], cols [3 ])
46- return None # No match
47-
48-
49- def load_p_snp_data (id , spe , shuffle = "struct" ):
50- """Load the probability of SNP significance at each residue of the given
33+ cols = line .split ("\t " )
34+ if cols [0 ][4 :- 4 ].upper () == ara_id and cols [1 ][4 :- 4 ].upper () == pop_id : # Both ID match
35+ return (cols [0 ][4 :- 4 ].upper (), cols [1 ][4 :- 4 ].upper (), cols [2 ], cols [3 ])
36+ if cols [0 ][4 :- 4 ].upper () == ara_id and pop_id is None : # Ara ID match, fill Pop
37+ return (cols [0 ][4 :- 4 ].upper (), cols [1 ][4 :- 4 ].upper (), cols [2 ], cols [3 ])
38+ if cols [1 ][4 :- 4 ].upper () == pop_id and ara_id is None : # Pop ID match, fill Ara
39+ return (cols [0 ][4 :- 4 ].upper (), cols [1 ][4 :- 4 ].upper (), cols [2 ], cols [3 ])
40+ return None # No match
41+
42+
43+ def load_p_snp_data (id , spe , shuffle = "struct" ):
44+ """Load the probability of SNP significance at each residue of the given
5145 protein from the cache file.
5246
5347 :param id: str of TAIR10 or Pop v3.0 gene ID
5448 :param spe: Either "ara" or "pop" based on id species
55- :param shuffle: Either "struct" or "seq" for significance method.
49+ :param shuffle: Either "struct" or "seq" for significance method.
5650 Defaults to "struct"
57- :returns: List of significance scores at residue positions,
51+ :returns: List of significance scores at residue positions,
5852 or None if invalid ID.
5953 :rtype: List[float] or None
6054
@@ -74,14 +68,14 @@ def load_p_snp_data(id, spe, shuffle = "struct"):
7468 p_snps_file = POP_SEQ_P_SNP
7569 else :
7670 return None
77- else :
71+ else :
7872 return None
7973
8074 # Load data from file
81- with open (p_snps_file , 'r' ) as f_p_snps :
75+ with open (p_snps_file , "r" ) as f_p_snps :
8276 for line in f_p_snps :
8377 if line .upper ().startswith (id ):
84- return [float (p ) for p in line .split (' \t ' )[1 ].split (',' )]
78+ return [float (p ) for p in line .split (" \t " )[1 ].split ("," )]
8579 return None
8680
8781
@@ -98,34 +92,34 @@ def mark_significant(null_probs, p):
9892
9993
10094def match_residues (aln ):
101- """For each index in the first protein which aligns with a residue in the
95+ """For each index in the first protein which aligns with a residue in the
10296 second protein (not a gap), provide the aligned index in the second protein.
10397
10498 :param aln: Tuple of Arabidopsis and Poplar protein sequences
10599 :returns: Dict from index in first protein to index in second protein
106100 :rtype: Dict[int, int]
107101 """
108102 matchings = {}
109- curr_prot1 = 1 # Current index in first protein
110- curr_prot2 = 1 # Current index in second protein
103+ curr_prot1 = 1 # Current index in first protein
104+ curr_prot2 = 1 # Current index in second protein
111105 # Iterate over all positions in the proteins
112- for i in range (len (aln [0 ])):
106+ for i in range (len (aln [0 ])):
113107 # If both not gaps, match the indices
114- if ( aln [0 ][i ] != '-' and aln [1 ][i ] != '-' ) :
115- matchings [curr_prot1 ] = curr_prot2
108+ if aln [0 ][i ] != "-" and aln [1 ][i ] != "-" :
109+ matchings [curr_prot1 ] = curr_prot2
116110 # If the position in the first protein is not a gap, increment index
117- if aln [0 ][i ] != '-' :
111+ if aln [0 ][i ] != "-" :
118112 curr_prot1 += 1
119113 # If the position in the second protein is not a gap, increment index
120- if aln [1 ][i ] != '-' :
114+ if aln [1 ][i ] != "-" :
121115 curr_prot2 += 1
122116 return matchings
123117
124118
125119def significant_in_both (sig1 , sig2 , aln_matching ):
126- """Mark a residue as significant in both if it aligns with a residue
120+ """Mark a residue as significant in both if it aligns with a residue
127121 in the other protein and those residues are both marked significant.
128-
122+
129123 :param sig1: Boolean list of significant residues in protein 1.
130124 :param sig2: Boolean list of significant residues in protein 2.
131125 :param aln_matching: Dictionary of aligned indices from protein 1 to 2.
@@ -136,7 +130,7 @@ def significant_in_both(sig1, sig2, aln_matching):
136130 both_sig1 = [False ] * len (sig1 )
137131 both_sig2 = [False ] * len (sig2 )
138132 # For each aligned index in protein 1
139- for i in aln_matching :
133+ for i in aln_matching :
140134 # If that reside and the counterpart in protein 2 are both significant,
141135 # mark those residues as significant in both.
142136 # Dictionary stores 1-indexed positions, list is 0-indexed
@@ -158,11 +152,11 @@ def get_sig_index(sig):
158152
159153# Find hotspot clusters
160154def cluster_components (hotspots , neighbours ):
161- """Determine clusters of hotspots residues. Clusters are connected
155+ """Determine clusters of hotspots residues. Clusters are connected
162156 components in the neighbourhood graph, found using DFS.
163157
164158 BFS algorithm: Initialize a frontier of positions to explore.
165- Add the starting residue to the frontier. While the frontier is not empty,
159+ Add the starting residue to the frontier. While the frontier is not empty,
166160 remove the last item from the frontier and add its significant non-explored
167161 residues to the frontier.
168162
@@ -173,25 +167,24 @@ def cluster_components(hotspots, neighbours):
173167 :rtype: List[List[int]]
174168 """
175169 clusters = []
176- # Track explored residues. Set non-hotspot residues to explored, so the
170+ # Track explored residues. Set non-hotspot residues to explored, so the
177171 # algorithm will not visit them.
178172
179173 explored = [not hotspot for hotspot in hotspots ]
180174 residue_frontier = []
181175 # explored and residue_fronter is 0-indexed
182176 # neighbors and clusters is 1-indexed.
183- for i in range (len (hotspots )): # O-indexed i
177+ for i in range (len (hotspots )): # O-indexed i
184178 if not explored [i ]:
185- residue_frontier .append (i ) # Expand: add to frontier, set explored
179+ residue_frontier .append (i ) # Expand: add to frontier, set explored
186180 explored [i ] = True
187181 curr_cluster = []
188182 while len (residue_frontier ) > 0 :
189183 curr_res = residue_frontier .pop ()
190- curr_cluster .append (curr_res + 1 ) # Add to cluster
191- for neighbour in neighbours [i + 1 ]: # Push neighbours
184+ curr_cluster .append (curr_res + 1 ) # Add to cluster
185+ for neighbour in neighbours [i + 1 ]: # Push neighbours
192186 if not explored [neighbour - 1 ]:
193187 residue_frontier .append (neighbour - 1 )
194188 explored [neighbour - 1 ] = True
195- clusters .append (curr_cluster ) # Save component on DFS finish
189+ clusters .append (curr_cluster ) # Save component on DFS finish
196190 return clusters
197-
0 commit comments