feat(alignment): add RMSD calculation and FASTA export

Author: AmirMohammad Cheraghali

src/components/SequenceAlignmentModal.tsx

@@ -1,5 +1,5 @@
 import React, { useMemo, useState } from 'react';
-import { X, GitCommitVertical, AlertTriangle, FileText, BarChart2, Hash, Percent } from 'lucide-react';
+import { X, GitCommitVertical, AlertTriangle, FileText, BarChart2, Hash, Percent, Download, Activity, Ruler } from 'lucide-react';
 import type { ChainInfo, SuperposedStructure } from '../types';
 import clsx from 'clsx';
 
@@ -20,30 +20,26 @@ interface AlignedResult {
             identity: number; // %
             similarity: number; // %
             gaps: number; // count
-            length: number; // total alignment length
-            score: number; // raw score
+            length: number; // length
+            score: number;
+            rmsd?: number; // Angstrom
         };
         alignment: {
             seq1: string; // Primary
             seq2: string; // Target
-            matchStr: string; // Match string ( | : . )
+            matchStr: string;
         };
     }[];
 }
 
-// Scientific Residue Coloring (Clustal-inspired simplified)
+// Scientific Residue Coloring
 const RESIDUE_COLORS: Record<string, string> = {
-    // Hydrophobic (Blue)
     'A': 'text-blue-400', 'V': 'text-blue-400', 'L': 'text-blue-400', 'I': 'text-blue-400',
     'M': 'text-blue-400', 'F': 'text-blue-400', 'W': 'text-blue-400', 'P': 'text-blue-400',
-    // Polar (Green)
     'G': 'text-green-400', 'S': 'text-green-400', 'T': 'text-green-400', 'C': 'text-yellow-400',
     'N': 'text-green-400', 'Q': 'text-green-400', 'Y': 'text-green-400',
-    // Positive (Red)
     'K': 'text-red-400', 'R': 'text-red-400', 'H': 'text-red-400',
-    // Negative (Magenta)
     'D': 'text-fuchsia-400', 'E': 'text-fuchsia-400',
-    // Gap
     '-': 'text-neutral-700'
 };
 
@@ -70,17 +66,21 @@ const CONSERVATIVE_GROUPS = [
 const getMatchChar = (a: string, b: string): string => {
     if (a === '-' || b === '-') return ' ';
     if (a === b) return '|';
-    // Check conservative
     for (const group of CONSERVATIVE_GROUPS) {
         if (group.includes(a) && group.includes(b)) return ':';
     }
     return '.';
 };
 
-const alignSequences = (seq1: string, seq2: string) => {
+const alignSequences = (chain1: ChainInfo, chain2: ChainInfo) => {
+    const seq1 = chain1.sequence;
+    const seq2 = chain2.sequence;
+    const coords1 = chain1.coords;
+    const coords2 = chain2.coords;
+
     const match = 10;
     const mismatch = -2;
-    const gap = -5; // Penalty
+    const gap = -5;
 
     const n = seq1.length;
     const m = seq2.length;
@@ -108,18 +108,36 @@ const alignSequences = (seq1: string, seq2: string) => {
     let similarityCount = 0;
     let gapCount = 0;
 
+    // RMSD Calc
+    let sumSqDist = 0;
+    let atomPairs = 0;
+
     while (i > 0 || j > 0) {
         if (i > 0 && j > 0 && scoreMatrix[i][j] === scoreMatrix[i - 1][j - 1] + (seq1[i - 1] === seq2[j - 1] ? match : mismatch)) {
             const c1 = seq1[i - 1];
             const c2 = seq2[j - 1];
             align1 = c1 + align1;
             align2 = c2 + align2;
+
+            // Stats
             if (c1 === c2) {
                 identityCount++;
                 similarityCount++;
             } else if (getMatchChar(c1, c2) === ':') {
                 similarityCount++;
             }
+
+            // RMSD: Calculate dist between coords1[i-1] and coords2[j-1]
+            if (coords1 && coords2 && coords1[i - 1] && coords2[j - 1]) {
+                const p1 = coords1[i - 1];
+                const p2 = coords2[j - 1];
+                const dx = p1.x - p2.x;
+                const dy = p1.y - p2.y;
+                const dz = p1.z - p2.z;
+                sumSqDist += dx * dx + dy * dy + dz * dz;
+                atomPairs++;
+            }
+
             i--; j--;
         } else if (i > 0 && scoreMatrix[i][j] === scoreMatrix[i - 1][j] + gap) {
             align1 = seq1[i - 1] + align1;
@@ -137,6 +155,9 @@ const alignSequences = (seq1: string, seq2: string) => {
     const length = align1.length;
     const matchStr = align1.split('').map((c, k) => getMatchChar(c, align2[k])).join('');
 
+    // RMSD Final
+    const rmsd = atomPairs > 0 ? Math.sqrt(sumSqDist / atomPairs) : undefined;
+
     return {
         seq1: align1,
         seq2: align2,
@@ -146,7 +167,8 @@ const alignSequences = (seq1: string, seq2: string) => {
             similarity: (similarityCount / length) * 100,
             gaps: gapCount,
             length,
-            score: scoreMatrix[n][m]
+            score: scoreMatrix[n][m],
+            rmsd
         }
     };
 };
@@ -168,12 +190,10 @@ export const SequenceAlignmentModal: React.FC<SequenceAlignmentModalProps> = ({
             const chainMatches: AlignedResult['chainMatches'] = [];
 
             primaryStructure.forEach(pChain => {
-                // Try to find same chain name, else index?
-                // For exact matches we prefer name.
                 const targetChain = ov.chains?.find(c => c.name === pChain.name) || ov.chains?.[0];
 
                 if (targetChain) {
-                    const result = alignSequences(pChain.sequence, targetChain.sequence);
+                    const result = alignSequences(pChain, targetChain);
                     chainMatches.push({
                         primaryChain: pChain.name,
                         targetChain: targetChain.name,
@@ -202,6 +222,21 @@ export const SequenceAlignmentModal: React.FC<SequenceAlignmentModalProps> = ({
         if (!selectedChain && availableChains.length > 0) setSelectedChain(availableChains[0]);
     }, [availableChains]);
 
+    const handleExport = (result: AlignedResult, match: AlignedResult['chainMatches'][0]) => {
+        const text = `>Primary | Chain ${match.primaryChain}
+${match.alignment.seq1}
+>Overlay: ${result.overlayName} | Chain ${match.targetChain}
+${match.alignment.seq2}
+`;
+        const blob = new Blob([text], { type: 'text/plain' });
+        const url = URL.createObjectURL(blob);
+        const a = document.createElement('a');
+        a.href = url;
+        a.download = `alignment_primary_${match.primaryChain}_vs_${result.overlayName}_${match.targetChain}.fasta`;
+        a.click();
+        URL.revokeObjectURL(url);
+    };
+
     if (!isOpen) return null;
 
     return (
@@ -250,7 +285,7 @@ export const SequenceAlignmentModal: React.FC<SequenceAlignmentModalProps> = ({
                         return (
                             <div key={result.overlayId} className="bg-black/40 border border-white/10 rounded-xl overflow-hidden">
                                 {/* Stats Dashboard */}
-                                <div className="grid grid-cols-2 md:grid-cols-4 divide-x divide-white/10 border-b border-white/10 bg-white/5">
+                                <div className="grid grid-cols-2 md:grid-cols-5 divide-x divide-white/10 border-b border-white/10 bg-white/5">
                                     <StatBox
                                         label="Identity"
                                         value={`${match.stats.identity.toFixed(1)}%`}
@@ -263,6 +298,13 @@ export const SequenceAlignmentModal: React.FC<SequenceAlignmentModalProps> = ({
                                         icon={<Hash size={14} />}
                                         color="text-blue-400"
                                     />
+                                    <StatBox
+                                        label="RMSD"
+                                        value={match.stats.rmsd ? `${match.stats.rmsd.toFixed(2)} Å` : "N/A"}
+                                        icon={<Activity size={14} />}
+                                        color="text-cyan-400"
+                                        subtext={match.stats.rmsd ? "Cα Atoms" : "No Coordinates"}
+                                    />
                                     <StatBox
                                         label="Gaps"
                                         value={match.stats.gaps.toString()}
@@ -280,9 +322,18 @@ export const SequenceAlignmentModal: React.FC<SequenceAlignmentModalProps> = ({
 
                                 <div className="p-6">
                                     <div className="flex items-baseline justify-between mb-4">
-                                        <h3 className="text-lg font-bold text-neutral-200">
-                                            {result.overlayName} <span className="text-neutral-500 text-sm font-normal">(Chain {match.targetChain})</span>
-                                        </h3>
+                                        <div className="flex items-center gap-3">
+                                            <h3 className="text-lg font-bold text-neutral-200">
+                                                {result.overlayName} <span className="text-neutral-500 text-sm font-normal">(Chain {match.targetChain})</span>
+                                            </h3>
+                                            <button
+                                                onClick={() => handleExport(result, match)}
+                                                className="flex items-center gap-1.5 px-3 py-1 bg-white/5 hover:bg-cyan-500/20 text-xs font-bold text-cyan-400 rounded-md transition-colors border border-white/10 hover:border-cyan-500/50"
+                                            >
+                                                <Download size={12} />
+                                                Export FASTA
+                                            </button>
+                                        </div>
                                         <div className="flex gap-4">
                                             {RESIDUE_GROUPS.map(g => (
                                                 <div key={g.name} className="flex items-center gap-2" title={g.desc}>
@@ -301,14 +352,12 @@ export const SequenceAlignmentModal: React.FC<SequenceAlignmentModalProps> = ({
                                             <div className="w-24 shrink-0" />
                                             <div className="flex relative h-4 w-full">
                                                 {Array.from({ length: Math.ceil(match.stats.length / 10) }).map((_, i) => (
-                                                    <span key={i} className="absolute text-[10px] text-neutral-500 border-l border-neutral-700 pl-1 h-3" style={{ left: `${i * 10 * 12}px` }}> {/* Approx 12px per char? No, char width varies. Monospace char width is fixed but usually 1ch ~8-10px depending on font size. We need flex basis. Better to use simple indices in a flex row */}
+                                                    <span key={i} className="absolute text-[10px] text-neutral-500 border-l border-neutral-700 pl-1 h-3" style={{ left: `${i * 10 * 12}px` }}>
                                                         {i * 10 + 1}
                                                     </span>
                                                 ))}
                                             </div>
                                         </div>
-                                        {/* Revisit Ruler: Absolute positioning on text is hard without fixed width font metrics. 
-                                            Let's just put markers every 10 chars in the flex flow. */}
 
                                         <div className="space-y-1">
                                             {/* Primary Sequence */}
@@ -368,20 +417,9 @@ const SequenceRow = ({ label, sequence, compareSeq }: { label: string, sequence:
         <span className="w-24 shrink-0 text-xs font-bold text-neutral-500 uppercase tracking-wider select-none pl-2 group-hover:text-neutral-300 transition-colors">
             {label}
         </span>
-        <div className="flex font-mono text-sm tracking-widest"> {/* tracking-widest ~ 0.1em. need fixed width */}
+        <div className="flex font-mono text-sm tracking-widest">
             {sequence.split('').map((char, i) => {
                 let colorClass = RESIDUE_COLORS[char] || 'text-neutral-300';
-                if (compareSeq) {
-                    // If it's the target seq, maybe dim non-mismatches?
-                    // Or just use standard colors. Standard colors look "scientific".
-                    // If difference:
-                    const c1 = compareSeq[i];
-                    if (c1 !== '-' && char !== '-' && c1 !== char) {
-                        // Mismatch implies importance? Or just typical color?
-                        // Keep typical color but maybe background highlight? 
-                        // Let's stick to text color for clean look.
-                    }
-                }
                 return (
                     <span key={i} className={`w-[1ch] inline-block text-center ${colorClass}`}>
                         {char}

src/types.ts

@@ -46,6 +46,7 @@ export interface ChainInfo {
     type?: 'protein' | 'nucleic' | 'unknown'; // Added: To distinguish polymer type
     atoms?: AtomInfo[]; // Added: For small chemicals, list atoms directly
     bFactors?: number[]; // Added: Per-residue B-factor for coloring sync
+    coords?: { x: number; y: number; z: number }[]; // Added: C-alpha/P coordinates for RMSD
 }
 
 export interface StructureInfo {

src/utils/pdbUtils.ts

@@ -142,6 +142,7 @@ export const extractChainsFromComponent = (component: any): { chains: ChainInfo[
 
         const resMap: number[] = [];
         const bFactors: number[] = [];
+        const coords: { x: number, y: number, z: number }[] = [];
 
         try {
             c.eachResidue((r: any) => {
@@ -158,16 +159,30 @@ export const extractChainsFromComponent = (component: any): { chains: ChainInfo[
                     resMap.push((maxSeq > -Infinity ? maxSeq : 0) + 1);
                 }
 
-                // B-Factor Extraction
+                // B-Factor & Coordinate Extraction
                 let bSum = 0;
                 let bCount = 0;
+                let caAtom: any = null;
+                let firstAtom: any = null;
+
                 r.eachAtom((a: any) => {
                     bSum += a.bfactor;
                     bCount++;
+
+                    if (a.atomname === 'CA' || a.atomname === 'P') caAtom = a;
+                    if (!firstAtom) firstAtom = a;
                 });
                 const avgB = bCount > 0 ? bSum / bCount : 0;
                 bFactors.push(avgB);
 
+                // Store coords (CA/P or fallback to first atom)
+                const targetAtom = caAtom || firstAtom;
+                if (targetAtom) {
+                    coords.push({ x: targetAtom.x, y: targetAtom.y, z: targetAtom.z });
+                } else {
+                    coords.push({ x: 0, y: 0, z: 0 }); // Should theoretically not happen
+                }
+
                 // Determine Type
                 if (r.isNucleic()) nucleicCount++;
                 else if (r.isProtein()) proteinCount++;
@@ -224,7 +239,8 @@ export const extractChainsFromComponent = (component: any): { chains: ChainInfo[
             residueMap: resMap,
             type: chainType,
             atoms: atomList.length > 0 ? atomList : undefined,
-            bFactors: bFactors
+            bFactors: bFactors,
+            coords: coords // Added field
         });
     });