Skip to content

Commit 346e966

Browse files
committed
src/dssp.cpp: a few changes/fixes from Claude, improving β-strand detection
1 parent 80c00e5 commit 346e966

File tree

1 file changed

+87
-12
lines changed

1 file changed

+87
-12
lines changed

src/dssp.cpp

Lines changed: 87 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -265,18 +265,18 @@ BridgeType DsspCalculator::calculate_bridge_type(size_t i, size_t j) const {
265265
if (i >= res_infos.size() || j >= res_infos.size())
266266
return BridgeType::None;
267267

268-
// antiparallel
268+
// antiparallel - either condition can create a bridge
269269
bool anti1 = has_hbond_between(res_infos[i], res_infos[j]) && has_hbond_between(res_infos[j], res_infos[i]);
270270
bool anti2 = (i > 0 && j + 1 < res_infos.size() && j > 0 && i + 1 < res_infos.size()) &&
271271
has_hbond_between(res_infos[i - 1], res_infos[j + 1]) && has_hbond_between(res_infos[j - 1], res_infos[i + 1]);
272-
if (anti1 && anti2)
272+
if (anti1 || anti2)
273273
return BridgeType::AntiParallel;
274-
// parallel
274+
// parallel - either condition can create a bridge
275275
bool para1 = (j + 1 < res_infos.size() && i + 1 < res_infos.size()) &&
276276
has_hbond_between(res_infos[i], res_infos[j + 1]) && has_hbond_between(res_infos[j], res_infos[i + 1]);
277277
bool para2 = (i > 0 && j > 0) &&
278278
has_hbond_between(res_infos[i - 1], res_infos[j]) && has_hbond_between(res_infos[j - 1], res_infos[i]);
279-
if (para1 && para2)
279+
if (para1 || para2)
280280
return BridgeType::Parallel;
281281

282282
return BridgeType::None;
@@ -293,29 +293,104 @@ void DsspCalculator::find_bridges_and_strands() {
293293
}
294294
}
295295

296-
// Find strands from bridge patterns
296+
// Find strands from bridge patterns - more permissive approach
297297
for (size_t i = 0; i < bridges_.size(); ++i) {
298298
for (size_t j = i + 1; j < bridges_.size(); ++j) {
299299
const Bridge& b1 = bridges_[i];
300300
const Bridge& b2 = bridges_[j];
301301
if (b1.type == b2.type && no_chain_breaks_between(b1.partner1, b1.partner2) && no_chain_breaks_between(b2.partner1, b2.partner2)) {
302-
if (b1.partner1 == b2.partner1 + 1 && b1.partner2 == b2.partner2 - 1 && b1.type == BridgeType::AntiParallel) {
303-
ss_info[b1.partner1].ss_type = SecondaryStructure::Strand;
304-
ss_info[b1.partner2].ss_type = SecondaryStructure::Strand;
305-
ss_info[b2.partner1].ss_type = SecondaryStructure::Strand;
306-
ss_info[b2.partner2].ss_type = SecondaryStructure::Strand;
302+
// More permissive conditions for strand formation
303+
// Allow bridges that are nearby (within 4 residues) to form strands
304+
bool forms_strand = false;
305+
306+
if (b1.type == BridgeType::AntiParallel) {
307+
// For antiparallel: check if bridges are reasonably close and in correct orientation
308+
if (abs((int)b1.partner1 - (int)b2.partner1) <= 4 && abs((int)b1.partner2 - (int)b2.partner2) <= 4) {
309+
forms_strand = true;
310+
}
311+
} else if (b1.type == BridgeType::Parallel) {
312+
// For parallel: check if bridges are reasonably close and in same direction
313+
if (abs((int)b1.partner1 - (int)b2.partner1) <= 4 && abs((int)b1.partner2 - (int)b2.partner2) <= 4) {
314+
forms_strand = true;
315+
}
307316
}
308-
if (b1.partner1 == b2.partner1 + 1 && b1.partner2 == b2.partner2 + 1 && b1.type == BridgeType::Parallel) {
317+
318+
if (forms_strand) {
319+
// Mark the bridge positions as strands
309320
ss_info[b1.partner1].ss_type = SecondaryStructure::Strand;
310321
ss_info[b1.partner2].ss_type = SecondaryStructure::Strand;
311322
ss_info[b2.partner1].ss_type = SecondaryStructure::Strand;
312323
ss_info[b2.partner2].ss_type = SecondaryStructure::Strand;
324+
325+
// Extend strands between the bridges
326+
size_t start1 = std::min(b1.partner1, b2.partner1);
327+
size_t end1 = std::max(b1.partner1, b2.partner1);
328+
size_t start2 = std::min(b1.partner2, b2.partner2);
329+
size_t end2 = std::max(b1.partner2, b2.partner2);
330+
331+
// Mark residues between bridges as strands (if no chain breaks)
332+
for (size_t k = start1; k <= end1; ++k) {
333+
if (ss_info[k].ss_type == SecondaryStructure::Loop && no_chain_breaks_between(start1, k)) {
334+
ss_info[k].ss_type = SecondaryStructure::Strand;
335+
}
336+
}
337+
for (size_t k = start2; k <= end2; ++k) {
338+
if (ss_info[k].ss_type == SecondaryStructure::Loop && no_chain_breaks_between(start2, k)) {
339+
ss_info[k].ss_type = SecondaryStructure::Strand;
340+
}
341+
}
313342
}
314343
}
315344
}
316345
}
317346

318-
// Mark isolated bridges
347+
348+
// Simple strand extension: extend strands by 1-2 residues at termini
349+
std::vector<bool> extend_here(ss_info.size(), false);
350+
351+
for (size_t i = 0; i < ss_info.size(); ++i) {
352+
if (ss_info[i].ss_type == SecondaryStructure::Strand) {
353+
// Mark 1-2 residues before strand start for extension
354+
if (i > 0 && ss_info[i-1].ss_type == SecondaryStructure::Loop && no_chain_breaks_between(i-1, i)) {
355+
extend_here[i-1] = true;
356+
if (i > 1 && ss_info[i-2].ss_type == SecondaryStructure::Loop && no_chain_breaks_between(i-2, i-1)) {
357+
extend_here[i-2] = true;
358+
}
359+
}
360+
// Mark 1-2 residues after strand end for extension
361+
if (i + 1 < ss_info.size() && ss_info[i+1].ss_type == SecondaryStructure::Loop && no_chain_breaks_between(i, i+1)) {
362+
extend_here[i+1] = true;
363+
if (i + 2 < ss_info.size() && ss_info[i+2].ss_type == SecondaryStructure::Loop && no_chain_breaks_between(i+1, i+2)) {
364+
extend_here[i+2] = true;
365+
}
366+
}
367+
}
368+
}
369+
370+
// Apply extensions
371+
for (size_t i = 0; i < ss_info.size(); ++i) {
372+
if (extend_here[i]) {
373+
ss_info[i].ss_type = SecondaryStructure::Strand;
374+
}
375+
}
376+
377+
// Convert bridges to strands if they're adjacent to existing strands
378+
for (const auto& bridge : bridges_) {
379+
bool convert_to_strand = false;
380+
381+
// Check if bridge is adjacent to strands
382+
if (bridge.partner1 > 0 && ss_info[bridge.partner1 - 1].ss_type == SecondaryStructure::Strand) convert_to_strand = true;
383+
if (bridge.partner1 + 1 < ss_info.size() && ss_info[bridge.partner1 + 1].ss_type == SecondaryStructure::Strand) convert_to_strand = true;
384+
if (bridge.partner2 > 0 && ss_info[bridge.partner2 - 1].ss_type == SecondaryStructure::Strand) convert_to_strand = true;
385+
if (bridge.partner2 + 1 < ss_info.size() && ss_info[bridge.partner2 + 1].ss_type == SecondaryStructure::Strand) convert_to_strand = true;
386+
387+
if (convert_to_strand) {
388+
ss_info[bridge.partner1].ss_type = SecondaryStructure::Strand;
389+
ss_info[bridge.partner2].ss_type = SecondaryStructure::Strand;
390+
}
391+
}
392+
393+
// Mark remaining isolated bridges
319394
for (const auto& bridge : bridges_) {
320395
if (ss_info[bridge.partner1].ss_type == SecondaryStructure::Loop)
321396
ss_info[bridge.partner1].ss_type = SecondaryStructure::Bridge;

0 commit comments

Comments
 (0)