🔬 Theoretically and Practically Faster Algorithms for Protein Structure Alignment

This repository contains an implementation of algorithms for the Protein Structure Alignment problem, formulated as the sequential Largest Common Point-set (sLCP) problem under the bottleneck distance.

The implementation accompanies the paper:

Theoretically and Practically Faster Algorithms for Protein Structure Alignment
Masahito Tsukahara and Tetsuo Shibuya

---

📖 Overview

Protein structure alignment is a fundamental task in structural bioinformatics. Given two proteins represented as ordered sequences of Cα atoms, the goal is to find a largest order-preserving common substructure that can be aligned within a given distance threshold.

This work studies the problem through the sequential LCP formulation under the bottleneck distance. It provides both:

• a faster polynomial-time exact algorithm from a theoretical perspective, and
• a practically faster approximation algorithm with the same solution-size guarantee as the previous guaranteed framework.

The approximation algorithm is designed for realistic protein data and is evaluated on structures from the Protein Data Bank (PDB).

---

✨ Key Contributions

• ⚙️ Improved exact algorithm:
  We improve the theoretical running time of a reliable exact algorithm from
  O(n^14) to O(Δ_ε n^13 log n)
  by exploiting the bounded-density property of protein structures.

• 🔎 Faster guaranteed approximation:
  We introduce a filtering technique for the approximation framework of Li and Ng and Tsukahara--Shibuya.

• 🛡️ Preserved theoretical guarantee:
  The improved approximation algorithm preserves the guarantee that, for approximation factor 0 < delta <= 3, it finds a subsequence whose size is at least the optimum under threshold epsilon, while allowing alignment distance (1 + delta) epsilon.

• 🚀 Substantial practical speedup:
  On PDB-derived benchmark instances, the proposed method achieves over a 24-fold speedup at n = 130 compared with the fastest previous guaranteed baseline.

• 🧬 Larger practical input range:
  Previous guaranteed methods became costly around n = 120--130.
  The proposed method processes:

  • n = 120 in 5.83 minutes on average,
  • n = 130 in 9.50 minutes on average,
  • n = 200 in 2.33 hours on average,
  • n = 210 in 3.33 hours on average.

• 🌍 Expanded PDB coverage:
  Since the median PDB entity length is around 210 residues in the analyzed dataset, the proposed method extends the practically accessible range from roughly the shortest 25% of PDB entities to nearly 50%.

---

📂 Repository Structure

.
├── src/             # C++ implementation and Python visualizer
├── include/         # Project and third-party C++ headers
├── inputs/          # Random PDB input generator
├── samples/         # Original and 10-residue sample PDB files
├── results/         # Transformed coordinates produced by the aligner
├── ReadMe.md        # Project overview and usage
├── Makefile         # Build instructions
└── LICENSE          # License information

---

📦 Getting Started

The commands below provide a minimal end-to-end test: install the dependencies, build the executable, align the two small sample structures, and visualize the transformed coordinates. Run all commands from the repository root unless otherwise noted.

Requirements

• Git
• GNU Make
• A C++20 compiler (the supplied Makefile uses g++)
• Eigen 3.4.0
• p-ranav/argparse v3.2
• Python 3 with numpy and plotly (only required for visualization)
• Python 3 with biopython (only required for random PDB input generation)

Eigen and argparse are required to compile the aligner; they are not optional.

1. Clone the repository

git clone https://github.com/masat03110/ProteinStructureAlignmentPlus.git
cd ProteinStructureAlignmentPlus

2. Install the C++ header dependencies

The source code expects the following two paths to exist:

include/eigen-3.4.0/Eigen/Dense
include/argparse/argparse.hpp

One way to create that layout is:

git clone --depth 1 --branch 3.4.0 \
  https://gitlab.com/libeigen/eigen.git include/eigen-3.4.0

git clone --depth 1 --branch v3.2 \
  https://github.com/p-ranav/argparse.git /tmp/argparse-v3.2
cp -R /tmp/argparse-v3.2/include/argparse include/

If these libraries are already installed elsewhere, edit CXXFLAGS in the Makefile so that both headers above are discoverable by the compiler.

3. Build

make

This creates ./aligner. To remove the executable and rebuild it:

make clean
make

4. Run the small example

./aligner samples/1pyv_trimmed.pdb samples/7my8_trimmed.pdb \
  --Algorithm AlignFastPlus --epsilon 0.5 --delta 1.0

The example contains 10 C-alpha atoms in each input and should finish quickly. A successful run prints information similar to:

Number of CA atoms in samples/1pyv_trimmed.pdb: 10
Number of CA atoms in samples/7my8_trimmed.pdb: 10
Epsilon value: 0.500000
Delta value: 1.000000
Using AlignFastPlus algorithm.
Solution size with AlignFastPlus: 6
AlignFastPlus completed.
Results are outputted to outputP.txt and outputQ.txt in the results directory.

Small floating-point differences across compilers and platforms are possible. The command overwrites results/outputP.txt and results/outputQ.txt.

5. Visualize the result

Create an isolated Python environment and install the visualization dependencies:

python3 -m venv .venv
source .venv/bin/activate
python -m pip install --upgrade pip
python -m pip install numpy plotly

Generate index.html using the same epsilon value as the alignment:

python src/plot.py --epsilon 0.5
python -m http.server 8000

Open http://localhost:8000/index.html in a browser. Stop the local server with Ctrl+C.

6. Generate random PDB inputs

Use inputs/sample_pdb_segments.py to select random protein entities from the RCSB Protein Data Bank and extract contiguous C-alpha segments. The first argument specifies the segment length, and the second specifies the number of files to generate. For example, 10 20 requests 20 files containing 10 consecutive C-alpha atoms each:

python inputs/sample_pdb_segments.py 10 20

Generated files can be passed directly to the aligner:

./aligner inputs/first_generated_file.pdb inputs/second_generated_file.pdb \
  --Algorithm AlignFastPlus --epsilon 2.0 --delta 3.0

Replace the two example filenames with names printed by the generator. See inputs/README.md for dependency installation, options, sampling details, and limitations.

For the complete command-line help, run:

./aligner -h
python src/plot.py -h
python inputs/sample_pdb_segments.py -h

---

⚙️ Command-Line Options

./aligner PATH1 PATH2 [--epsilon VALUE] [--delta VALUE] [--Algorithm NAME]

Argument	Meaning	Default
PATH1, PATH2	Input PDB files	required
--epsilon	Distance threshold in angstroms; use a value greater than 0	0.5
--delta	Approximation parameter; use 0 < delta <= 3	1.0
--Algorithm	Algorithm implementation listed below; the option name is case-sensitive	AlignFastLIS

delta is used by the approximation algorithms and does not affect BasicAlign. The program currently does not reject every out-of-range numeric value, so callers should observe the ranges above.

Available algorithms

Name	Role	Practical note
BasicAlign	Simplified basic/exact baseline	Intended only for very small inputs because of its high running time
LiNgAlign	Simplified Li--Ng guaranteed-approximation baseline	Useful for comparison; can be substantially slower than the proposed methods
AlignFastLIS	Faster method using the LIS-based subroutine	Default implementation
AlignFastLCS	Faster method using the LCS-based subroutine	Alternative implementation for comparison
AlignFastPlus	Proposed filtered version of the faster method	Recommended for the supplied quick example

Running time grows rapidly with input length. Start with the 10-residue *_trimmed.pdb files before trying the full sample structures.

---

🧬 Input Format and Limitations

The executable accepts PDB files and extracts the coordinates of C-alpha (CA) atoms from ATOM records. The order of those records defines the point sequence used by the algorithm.

The current parser is intentionally simple:

• mmCIF input is not supported;
• HETATM records are ignored;
• multi-model, multi-chain, alternate-location, and insertion-code behavior is not configurable; and
• malformed or empty input files may not produce a descriptive error.

For predictable results, preprocess an input into a single ordered chain/model containing standard ATOM records. The files in samples/ show the expected format. See samples/README.md for their sources and citations.

Random benchmark inputs can be downloaded and generated with inputs/sample_pdb_segments.py. See inputs/README.md for installation, sampling, and reproducibility details.

---

📄 Output

The program reports the computed solution size on standard output and writes:

results/outputP.txt
results/outputQ.txt

Each line of these files is an x y z coordinate in angstroms. The files contain the transformed coordinates of all C-alpha atoms from the two input structures, in their original sequence order. They are used by src/plot.py to display the superposition.

The files are not lists of matched residue indices: their number of lines can therefore be larger than the reported solution size. The current command- line program does not export the selected correspondence itself.

src/plot.py reads these files from results/ by default and writes index.html to the directory from which it is invoked. A different results directory can be selected with --results_path.

---

🔁 Reproducibility Scope

The commands above reproduce the small end-to-end demonstration included in this repository. The paper-scale timing experiment used random consecutive C-alpha segments, 1000 trials per setting, and substantially larger inputs. The present repository does not include the complete benchmark dataset, sampling scripts, random seeds, or batch timing scripts needed to regenerate the aggregate timing table automatically. The table below is therefore a summary of the experiments reported in the accompanying paper, rather than an automated benchmark produced by make.

📊 Experimental Summary

The algorithms were evaluated on random consecutive Cα segments sampled from PDB proteins.

Basic settings:

• distance threshold: epsilon = 2.0 Å
• approximation factor: delta = 3
• equal input lengths: n = m
• 1000 trials for each setting

Metric	Previous Guaranteed Baseline	Proposed Method
Mean time at n = 120	1.54 hours	5.83 minutes
Mean time at n = 130	3.80 hours	9.50 minutes
Speedup at n = 130	Baseline	over 24x
Largest tested practical size	around 130	around 210
Mean time at n = 200	not evaluated	2.33 hours
Mean time at n = 210	not evaluated	3.33 hours
Theoretical guarantee	guaranteed	preserved

---

🎯 Applications

• Pairwise protein structure alignment
• Common 3D substructure detection
• Comparative structural bioinformatics
• Evaluation of protein model similarity under distance constraints

---

📝 Citation

If you use this repository, please cite the accompanying paper:

@inproceedings{tsukahara_shibuya_2026_protein_alignment,
  title     = {Theoretically and Practically Faster Algorithms for Protein Structure Alignment},
  author    = {Tsukahara, Masahito and Shibuya, Tetsuo},
  year      = {2026}
}

---

🤝 Contributing

Contributions, issues, and feature requests are highly welcome!

Feel free to open a pull request or submit an issue for discussion.

---

📫 Contact

For questions, discussions, or collaborations:

Masahito Tsukahara

Email: tsuka03@hgc.jp

---

⭐ Acknowledgments

Special thanks to the open-source community and the Protein Data Bank (PDB) for providing the invaluable datasets that made this research possible.