File size: 7,385 Bytes
e781e66 | 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 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 | import numpy as np
from sklearn.metrics import precision_score, recall_score
from pathlib import Path
def _read_relevant_lines(file_path: Path):
with open(file_path, 'r') as f:
lines = f.readlines()
lines = list(filter(lambda line: not line.lstrip().startswith("#"), lines)) # Ignore comment lines
return lines
def parse_bpseq_file(bpseq_file_path: Path):
lines = _read_relevant_lines(bpseq_file_path)
seq_len = len(lines)
seq = ''
pair_mat = np.zeros((seq_len, seq_len), dtype=np.float32)
for line in lines:
res_idx, res_tkn, pair_idx = line.rstrip().split()
seq += res_tkn
if pair_idx != '0':
pair_mat[int(res_idx) - 1, int(pair_idx) - 1] = 1.0
return seq, pair_mat
def parse_ct_file(ct_file_path: Path):
lines = _read_relevant_lines(ct_file_path)
seq_len = int(lines[0].split()[0])
pair_mat = np.zeros((seq_len, seq_len))
seq = ""
for line in lines[1:]:
i, nuc, *_, j, _ = line.split()
i, j = int(i), int(j)
if j > 0:
pair_mat[i - 1, j - 1] = 1.0
seq += nuc
return seq, pair_mat
def save_to_ct(ct_file_path: Path, sec_struct: np.array, seq: str):
with open(ct_file_path, "w") as f:
f.write(f"{len(seq)}\t{ct_file_path.stem}\n")
for i in range(len(seq)):
f.write(f"{i + 1}\t") # Base number: Index n
f.write(f"{seq[i]}\t") # Base
f.write(f"{i}\t") # Index n - 1
f.write(f"{i + 2}\t") # Index n + 1
f.write(f"{sec_struct[i].argmax() + 1 if sec_struct[i].sum() > 0.0 else 0}\t") # Pair index
f.write(f"{i + 1}\n") # Natural numbering
ST_SEQ_LINE_IDX = 0
ST_DBN_LINE_IDX = 1
def parse_st_file(st_file_path: Path):
lines = _read_relevant_lines(st_file_path)
seq = lines[ST_SEQ_LINE_IDX].rstrip()
db_notation = lines[ST_DBN_LINE_IDX].rstrip()
pair_mat = dot_bracket_to_2d_mat(db_notation)
return seq, pair_mat
def parse_sec_struct_file(sec_struct_file_path: Path):
if sec_struct_file_path.suffix == ".ct":
return parse_ct_file(sec_struct_file_path)
elif sec_struct_file_path.suffix == ".bpseq":
return parse_bpseq_file(sec_struct_file_path)
elif sec_struct_file_path.suffix == ".st":
return parse_st_file(sec_struct_file_path)
else:
raise NotImplementedError("Given secondary structure file type is not supported!")
def dot_bracket_to_2d_mat(db_notation: str):
seq_len = len(db_notation)
pair_mat = np.zeros((seq_len, seq_len))
# Initialize bracket stacks
stacks= {}
stacks["("] = stacks[")"] = []
stacks["["] = stacks["]"] = []
stacks["{"] = stacks["}"] = []
stacks["<"] = stacks[">"] = []
# Iterate through the dot-bracket notation and fill the 2D matrix accordingly
for i in range(seq_len):
current_tkn = db_notation[i]
if current_tkn in ("(", "[", "{", "<"):
stacks[current_tkn].append(i)
elif current_tkn in (")", "]", "}", ">"):
j = stacks[current_tkn].pop()
pair_mat[i, j] = 1.0
elif db_notation[i] == ".":
pass
else:
raise RuntimeError(f"Encountered unexpected symbol in dot-bracket notation string! (index {i}: '{db_notation[i]}')")
# Symmetrize pairing matrix
pair_mat = pair_mat + pair_mat.transpose()
pair_mat = np.minimum(pair_mat, 1.0)
return pair_mat
_SHARP_LOOP_DIST_THRESHOLD = 4
def _generate_sharp_loop_mask(seq_len):
mask = np.eye(seq_len, k=0, dtype=bool)
for i in range(1, _SHARP_LOOP_DIST_THRESHOLD):
mask = mask + np.eye(seq_len, k=i, dtype=bool) + np.eye(seq_len, k=-i, dtype=bool)
return mask
CANONICAL_PAIRS = ['AU', 'UA', 'GC', 'CG', 'GU', 'UG']
def _generate_canonical_pairs_mask(seq: str):
seq = seq.replace('T', 'U')
mask = np.zeros((len(seq), len(seq)), dtype=bool)
for i, nt_i in enumerate(seq):
for j, nt_j in enumerate(seq):
if f'{nt_i}{nt_j}' in CANONICAL_PAIRS:
mask[i, j] = True
return mask
def _clean_sec_struct(sec_struct: np.array, probs: np.array):
clean_sec_struct = np.copy(sec_struct)
tmp_probs = np.copy(probs)
tmp_probs[sec_struct < 1] = 0.0
while np.sum(tmp_probs > 0.0) > 0:
i, j = np.unravel_index(np.argmax(tmp_probs, axis=None), tmp_probs.shape)
tmp_probs[i, :] = tmp_probs[j, :] = 0.0
clean_sec_struct[i, :] = clean_sec_struct[j, :] = 0
tmp_probs[:, i] = tmp_probs[:, j] = 0.0
clean_sec_struct[:, i] = clean_sec_struct[:, j] = 0
clean_sec_struct[i, j] = clean_sec_struct[j, i] = 1
return clean_sec_struct
def prob_mat_to_sec_struct(probs: np.array, seq: str, threshold: float = 0.5, allow_nc_pairs: bool = False, allow_sharp_loops: bool = False):
assert np.all(np.isclose(probs, np.transpose(probs))), "Probability matrix must be symmetric!"
seq_len = probs.shape[-1]
allowed_pairs_mask = np.logical_not(np.eye(seq_len, dtype=bool))
if not allow_sharp_loops:
# Prevent pairings that would cause sharp loops
allowed_pairs_mask = np.logical_and(allowed_pairs_mask, ~_generate_sharp_loop_mask(seq_len))
if not allow_nc_pairs:
# Prevent non-canonical pairings
allowed_pairs_mask = np.logical_and(allowed_pairs_mask, _generate_canonical_pairs_mask(seq))
probs[~allowed_pairs_mask] = 0.0
sec_struct = np.greater(probs, threshold).astype(int)
sec_struct = _clean_sec_struct(sec_struct, probs)
return sec_struct
def _relax_ss(ss_mat: np.array) -> np.array:
# Pad secondary structure (because of cyclical rolling)
ss_mat = np.pad(ss_mat, ((1, 1), (1, 1)), mode='constant')
# Create relaxed pairs matrix
relax_pairs = \
np.roll(ss_mat, shift=1, axis=-1) + np.roll(ss_mat, shift=-1, axis=-1) +\
np.roll(ss_mat, shift=1, axis=-2) + np.roll(ss_mat, shift=-1, axis=-2)
# Add relaxed pairs into original matrix
relaxed_ss = ss_mat + relax_pairs
# Ignore cyclical shift and clip values
relaxed_ss = relaxed_ss[..., 1: -1, 1: -1]
relaxed_ss = np.clip(relaxed_ss, 0, 1)
return relaxed_ss
def ss_recall(target_ss: np.array, pred_ss: np.array, allow_flexible_pairings: bool = True) -> float:
if allow_flexible_pairings:
pred_ss = _relax_ss(pred_ss)
seq_len = target_ss.shape[-1]
upper_tri_idcs = np.triu_indices(seq_len, k=1)
return recall_score(target_ss[upper_tri_idcs], pred_ss[upper_tri_idcs], zero_division=0.0)
def ss_precision(target_ss: np.array, pred_ss: np.array, allow_flexible_pairings: bool = True) -> float:
if allow_flexible_pairings:
target_ss = _relax_ss(target_ss)
seq_len = target_ss.shape[-1]
upper_tri_idcs = np.triu_indices(seq_len, k=1)
return precision_score(target_ss[upper_tri_idcs], pred_ss[upper_tri_idcs], zero_division=0.0)
EPSILON = 1e-5
def ss_f1(target_ss: np.array, pred_ss: np.array, allow_flexible_pairings: bool = True) -> float:
precision = ss_precision(target_ss, pred_ss, allow_flexible_pairings=allow_flexible_pairings)
recall = ss_recall(target_ss, pred_ss, allow_flexible_pairings=allow_flexible_pairings)
# Prevent division with 0.0
if precision + recall < EPSILON:
return 0.0
return (2 * precision * recall) / (precision + recall)
|