Contact MapContact maps and interface contacts
誰碰到誰，以及碰了多久。介面分析最耐用的一把尺。
01一句話只有一句
數出介面上有多少對重原子靠得夠近，以及它們維持了多久。
02Why should I care?我為什麼要讀這頁
因為在你手上通過可靠性閘門的那幾個量裡，只有它同時做到「CV 個位數」與「拆得開到股別」，所以它是這個 project 建議的 Gate 3 打分軸。
你在 5F9R 完整三元複合體上跑了 4 條獨立 replica，得到的跨 replica CV 是：總 protein–DNA 接觸 4.1%、target-DNA 接觸 9.4%。兩個都落在 reproducible tier。對照組是同一份資料裡的 interface iRMSD 23.9%、DNA RMSF 59.1%、groove 接觸計數 145.5%。接觸數不是「看起來比較好」，是差了一個量級。同一份資料裡跟它同級的還有 BSA（2.1%）與鹽橋（9.5%）。BSA 在原理上就拆不開（它是一個面積），鹽橋在你現有的分析裡則是只報了總數、沒有分股。差別在這裡，不在 CV。
而這件事對你來說應該一點都不陌生。你在 HEA 表面上不會報「第 472 號原子的位置偏移了 0.03 Å」，你報的是 coordination number、generalized coordination number、d-band center —— 全部都是聚合量。理由完全一樣：個別原子的瞬時位置是熱運動，聚合起來的配位環境才是描述符。接觸數就是把配位數的概念搬到蛋白–核酸介面上。而它繼承的警告也一樣：GCN 好用是因為它被聚合了，沒有人會拿單一原子的單一快照 CN 去宣稱那是活性位。
把 protein–DNA engagement 當聚合量報，永遠不要逐殘基報。 底下所有內容都是在解釋這句話為什麼成立、以及怎麼在 paper 裡把它寫得站得住。
03What這是什麼
定義。 給一個截止距離 rc，若蛋白殘基 i 與核酸殘基 j 之間最短的重原子對距離小於 rc，就說這一對在該幀有接觸。整個介面在該幀的狀態可以寫成一個布林矩陣 Cij(t) ∈ {0, 1}。
從這個矩陣可以往三個方向取：
| 取法 | 定義 | 性質 |
|---|---|---|
| 瞬時聚合數 | Nc(t) = Σij Cij(t) | 一條時間序列，可以對 replica 取平均，就是打分用的量 |
| 佔據率（occupancy） | Oij = ⟨Cij(t)⟩t | 0 到 1，就是「contact map」那張圖本身 |
| 逐殘基接觸度 | ki = Σj Oij | 每個殘基的介面參與度，本頁警告的主角 |
原子對還是殘基對，這是兩種不同的量，不要混講。
- 殘基對布林（上面的定義）：一對殘基不管靠得多密，都只算 1。優點是不被大殘基灌水，缺點是門檻附近的行為最劇烈 —— 一對殘基從 4.01 Å 走到 3.99 Å，整個貢獻從 0 跳到 1。
- 原子對計數：把所有在 rc 內的重原子對全部數進去。優點是它幾乎是連續的（一次只多一兩對），因此對 cutoff 的敏感度低很多；缺點是 Arg 這種大側鏈天生貢獻多，逐殘基比較時會有系統性偏差。
打分建議用原子級計數，因為它更平滑；但兩種都要在 methods 裡定義清楚，而且 CV 要分別報。
你自己那份 4-replica 分析數的既不是殘基對布林、也不是嚴格的原子對數，而是「與該核酸鏈任一重原子距離在截止內的蛋白重原子個數」，截止取 0.45 nm。總 protein–DNA 接觸的均值落在 400 上下，而且是 target 加 non-target 兩股，sgRNA 不含在裡面。CV 是綁在定義上的，換一個定義就要重新量一次，不能沿用。
Q（fraction of native contacts） 是另一個常見變形：先從參考結構定義一組「native contacts」，再看軌跡中保留了幾成。它的進階版用一個平滑的 sigmoid 開關函數取代硬截止，讓 Q 變成連續可微的量 —— 對於「接觸在門檻附近跳動」這個問題，這是最正統的解法。
04Why為什麼重要
為什麼聚合起來就穩了，逐殘基卻是雜訊？ 這不是玄學，是一個可以講清楚的統計事實。
總接觸數是幾百項的和。和的變異數不只是各項變異數相加，還要加上兩兩之間的共變異：
在一個緊密的蛋白–DNA 介面上，這些共變異項大量是負的。一個 Lys 側鏈離開 A 磷酸，通常代表它靠上了隔壁的 B 磷酸；DNA 骨架在溝槽裡滑一格，換掉的是接觸的身分，不是接觸的總數。負的共變異把和的變異壓下去，於是總數穩定（CV 4.1%），而構成它的每一項本身劇烈跳動（逐殘基 top-8 只共享 3/8）。
這兩件事同時成立，而且不矛盾。 能不能把這句話講清楚，就是這頁分析在 paper 裡站不站得住的關鍵。
物理上它代表什麼？ 接觸數是交互作用的粗粒化影子。protein–DNA 介面的主導項是靜電（磷酸骨架對 Lys/Arg）加上形狀互補，而接觸數量的是「有多少機會發生這些交互作用」，不對能量下任何承諾。它比 BSA 多了一個關鍵優勢：可分解。BSA 的 CV 更低（2.1%），但它是一個不透明的數字，你沒辦法問它「是 target 股還是 non-target 股」。接觸數可以按鏈拆成 target strand、non-target strand、sgRNA 三塊，也可以另外按區域（PAM-proximal / seed / PAM-distal）再切一次，而 grip 突變體的整個故事就活在這些分解裡。
為什麼這對 Cas9 特別重要。 eSpCas9(1.1)（三個突變 K848A／K1003A／R1060A）的設計邏輯就是純粹的 contact engineering：把 HNH、RuvC、PI 三個 domain 之間那條正電溝槽裡的殘基中性化，削弱對 non-target 股（R-loop 裡被 sgRNA 擠開、沒有配對的那一股）的握持。你自己量到的是 —— 三個 grip 位點裡只有 K848 一直貼著 non-target 股（4 條 replica 平均 0.7 nm），K1003 與 R1060 平均在一到兩奈米外（1.6 nm 與 1.5 nm）。但要誠實地加一句：這三個都是逐位點距離，跨 4 條 replica 的 CV 落在 38 到 46%，三個全是 noise tier —— R1060 在其中一條 replica 貼到 0.5 nm，在另一條卻是 2.0 nm。所以能講的只有「K848 一直貼著、另外兩個平均遠一些」這個量級判斷，不能倒過來說「隔那麼遠所以碰不到」，也不能拿小數點後兩位去排序。這正好是本站方法學主張的自我示範：逐位點量不能拿來下二分結論。即便如此，這個層次的敘事也只有接觸／距離分析講得出來，MM-PBSA 給你一個總數，講不出這件事。
05How怎麼做
定義層（先決定，再算，不要邊算邊改）：
- 先逐字寫下你數的是什麼。 你現有的 4-replica 分析用的是 0.45 nm 重原子截止（
CONTACT_CUT = 0.45），數的是靠近該核酸鏈的蛋白重原子個數。要延續那個 4.1% 就照用；要改成殘基對最短距離 ≤ 4.0 Å 的布林版本也可以，但那是另一個量，CV 要重新量。兩種都算一份最省事。 - 按鏈分解：target strand / non-target strand / sgRNA。這三項互斥，你的「總 protein–DNA 接觸」是前兩項相加。PAM duplex 不是第四項 —— 它是這兩條 DNA 鏈上的一段區域，跟按鏈的分解不在同一個軸上，混進同一張堆疊圖會重複計數。要看它就另開一個按區域的軸（PAM-proximal / seed / PAM-distal），並在圖說寫明是兩種不同的切法。
- 排除所有無實驗座標的建模補出區段（5F9R 的 non-target 股 res 1 到 11 就是這種：該鏈 30 nt 只有 19 個有密度），或至少單獨標色。
工具：
- mdtraj：
md.compute_contacts(traj, contacts=pairs, scheme='closest-heavy')直接回傳每個殘基對的最短重原子距離，再自己套門檻。核酸不要用scheme='ca'，核酸沒有 Cα，那個選項在這裡沒有意義。 - MDAnalysis：
MDAnalysis.analysis.contacts.Contacts（支援 hard cutoff 與 soft switching function），或直接用MDAnalysis.lib.distances.capped_distance自己組 —— 對 535k 原子的系統，先把 selection 縮到介面附近再算，否則你在浪費時間算兩個離很遠的原子。 - GROMACS：
gmx mindist搭配 index group 可以拿到殘基對最短距離；gmx select的動態選擇可以直接數「距離某組多近的原子數」。
統計層（這一段決定 reviewer 買不買單）：
- 每條 replica 對 production 時間平均，得到一個數字；4 條算 mean ± SD 與 CV。打分軸報的是這個，不是單條軌跡的曲線。
- cutoff 敏感度掃描：3.5 / 4.0 / 4.5 / 5.0 Å 各算一次（4.5 Å 是你現在的基準線），畫「結論 vs cutoff」。如果變體 A 在四個 cutoff 下都大於變體 B，你的結論就防彈了。這是全站投報率最高的一個檢查，成本幾乎是零。
- 佔據率門檻要明說：「occupancy > 0.5 視為穩定接觸」這種話一定要寫出來，而且要說明換成 0.3 或 0.7 結論會不會變。
- 逐殘基的東西只能當定性地圖，而且必須附跨 replica 一致性數字：top-N 重疊、pairwise Jaccard、Spearman ρ。
成本： 對 4 條 × 87 ns 的軌跡，全部分析加起來是分鐘等級。相對於產生軌跡的 H100 時數，接觸分析在帳面上等於免費。在「便宜」是硬約束的 workflow 裡，這件事本身就是它拿五星的理由之一。
06When什麼時候用
- 當 Gate 3 的打分軸。 這是它最主要的用途，也是這頁存在的理由。用總 protein–DNA 接觸數加上 target strand 分項，附 CV 與 cutoff 掃描。
- 比較 WT 與 grip 變體的握持幅度。 但要記得：eSpCas9(1.1) 在 on-target 底物上是 null，本來就該 null。要看到效應，底物必須帶 mismatch。用 on-target 底物比出「沒差別」不是發現，是設計錯誤。
- 把介面拆給讀者看。 按鏈的 target / non-target / sgRNA 分項，再配一張按區域（PAM-proximal / seed / PAM-distal）的，是 Cas9 論文裡最有敘事力的組合。兩張圖不要合成一張。
- 跟 BSA 交叉驗證。 兩者都是聚合幾何量，方向應該一致。若一個升一個降，先懷疑分析而不是先寫故事。
- 檢查起始結構品質。 把 MD 平均 contact map 跟晶體 contact map 算 Jaccard，是一個很便宜的「我的系統有沒有偏離已知結構」的檢查。
07When NOT什麼時候別用
-
不要逐殘基報，一次都不要。 為什麼會壞：你的 4 條 replica 的 top-8 接觸殘基只共享 3 個（[450, 695, 765]），pairwise Jaccard 平均 0.48。排名的變異來自「哪一個 Lys 側鏈在這一秒剛好轉向誰」，那是熱運動不是生物學。而且這份逐殘基統計是只對 target 股算的 —— 也就是整個介面上最穩的那一股，連在最有利的條件下都只剩 3/8。怎麼看出它壞了：算跨 replica 的 top-N 重疊、pairwise Jaccard 與 Spearman ρ。你的數字是 Jaccard 平均 0.48、ρ = +0.77 —— 排序趨勢還在（所以聚合量可信），但誰進 top-8 已經是擲骰子（所以逐殘基不可信）。這兩個數字必須一起看，只報 ρ 會讓人誤以為逐殘基也很穩。
它還壞得很好看，這才是真正的陷阱。 你的 top-8 名單裡，695 出現在全部四條、497 出現在其中三條 —— 這兩個正好是 SpCas9-HF1 四個突變位點（N497A／R661A／Q695A／Q926A）裡的兩個。逐殘基地圖確實落在文獻已知的工程位點上，所以它看起來是對的，所以你會更想拿它排序。它「命中已知位點」只證明你的介面選對了，不證明名單的成員身分可重現。
-
不要在起始結構沒過 QC 的情況下比較接觸。 為什麼會壞：元兇是起始結構，不是 metric。 同一套 metric、同一套 protocol，從 Chai-1 預測結構出發時 top-8 共享 0/8、ρ ≈ −0.03（等於隨機）；換成 5F9R 晶體結構出發就變成 3/8、ρ = +0.77。原因很清楚：Chai-1 在 HNH（pLDDT 62.9）與 target DNA（64 到 65）都是低信心區，那裡的介面幾何本來就是猜的，MD 只是把不同的猜測往不同方向鬆弛。怎麼看出它壞了：先看起始結構在介面區的 pLDDT／B-factor，再看跨 replica 一致性有沒有隨起始結構改變而改變。若會，問題在輸入端。
-
不要用 n = 2 判定可重現性。 為什麼會壞：小 n 會系統性高估一致性 —— 可比對的配對少、變異被低估。你自己踩過這個坑：2 條 replica 時逐殘基 top-8 看起來 6/8，非常漂亮；跑到 4 條掉到 3/8。怎麼看出它壞了：把 4 條裡所有 C(4,2) = 6 種兩兩組合各算一次一致性，看散布有多寬。你自己那 6 個 pairwise Jaccard 是 0.23 到 0.78 —— 如果你當初只跑到那個 0.78 的組合，你會很有把握地宣告可重現。散布跨越了你的判斷門檻，n = 2 的結論就不存在。
-
不要只用一個 cutoff 就下結論。 為什麼會壞：介面上有大量原子對的距離就落在 4 到 5 Å 這個帶裡，硬截止會讓它們在 0 與 1 之間反覆切換；兩個變體的排序完全可能在 4.0 Å 是 A > B、在 4.5 Å 翻成 B > A。怎麼看出它壞了：掃 3.5 到 5.0 Å，畫排序 vs cutoff。若翻轉，改用平滑開關函數，或退回 BSA 這種不需要選 cutoff 的量。
-
不要把 non-target strand 的接觸當可靠訊號。 為什麼會壞：non-target 股在 R-loop 裡是被擠開的單股，沒有配對約束，構形空間大得多，取樣需求也高得多。這不只是模擬的毛病 —— 5F9R 裡 non-target 股 res 1 到 11 根本沒有電子密度，pacesa2022 的 18 nt checkpoint 結構裡 non-target 股同樣是無序的。實驗看不清楚的那一段，你的 MD 不會替你變清楚。你的數字說明了一切 —— non-target 接觸 CV 17.7%（marginal），而 groove 接觸計數 CV 145.5%（純雜訊）。怎麼看出它壞了：分開報 target 與 non-target 的 CV，不要合併就完事。 合併之後的 4.1% 有一部分是被 target 股撐起來的漂亮數字，把它當成「整個介面都很穩」是自欺。
-
不要把接觸數當結合自由能。 為什麼會壞：接觸數沒有能量權重（一對 3.0 Å 的鹽橋和一對 3.9 Å 的凡德瓦接觸都算 1）、沒有去溶劑化代價、沒有熵。對這個 project 更致命的是方向性問題：高保真變體對 on/off-target 的結合親和力跟 WT 差不多，差別在 HNH 構形檢查點的動力學。接觸數再穩，也量不到動力學門檻。怎麼看出它壞了：拿 SpCas9-HF1 當測試案例，然後盯著看你的分數到底在追什麼。它的四個突變 N497A／R661A／Q695A／Q926A 就是把四根對 target 股磷酸骨架的直接接觸拿掉的，所以接觸數一定會降 —— 降的是你自己動手刪掉的那幾根。真正該對得上的是它的表型：on-target 活性幾乎沒掉、off-target 掉到 GUIDE-seq 測不到。那個 on／off 的鑑別力並沒有寫在接觸數的差值裡。若你的排序只是把變體按「刪了幾個接觸」排，它量的是擾動大小，不是專一性。
ρ = +0.77 這個數字看起來很棒，很容易讓人以為逐殘基分析過關了。但同一份資料的 top-8 只共享 3 個。排序趨勢可重現，成員身分不可重現。 在 paper 裡只報前者、不報後者，是這個領域最常見的一種善意的不誠實。兩個一起報，你的可信度反而更高。
08Project Lens跟我的 project 多相關
為什麼是五星（為什麼不是更低）： 因為在目前手上所有的觀測量裡，只有它同時滿足三個條件。
- 通過可靠性閘門。 總 protein–DNA 接觸 CV 4.1%、target-DNA 9.4%，都在 reproducible tier。與它同級的只有 BSA（2.1%）與鹽橋（9.5%）。
- 可分解成機制敘事。 BSA 的 CV 更漂亮，但它是一個不可拆的數字。接觸數可以拆成 target / non-target / sgRNA，而 grip 突變體的整個故事就住在這個分解裡 —— 例如你量到 eSpCas9(1.1) 三個 grip 位點只有 K848 一直貼著 non-target 股（平均 0.7 nm），K1003 與 R1060 平均在一到兩奈米外（1.6 與 1.5 nm；逐位點距離自身的 CV 是 38 到 46%，逐 replica 從 0.5 到 2.3 nm 都出現過，所以只能當量級判斷用）。這種話 BSA 講不出來，MM-PBSA 也講不出來。
- 成本近乎為零。 在「便宜」是硬約束的 project 裡，一個軌跡跑完就順手拿到、不用額外 GPU 時數的打分軸，價值被嚴重低估。
五星的意思是「在現有工具裡的優先序」，不是「它能回答終極問題」。 誠實地說，它有一個結構性天花板：接觸數是平衡結構的描述，而這個 project 已經確立的核心判斷是專一性由 HNH 構形檢查點的動力學決定，不是由平衡結合強度決定。接觸數量不到那件事。補這一段要靠 HNH 距離、enhanced sampling 或 kinetics 導向的設計，但那些目前都還沒有 4.1% 這種可靠性。在你有更好的東西之前，接觸數是唯一一個既通過可靠性閘門、又拆得出機制故事的打分軸。
具體怎麼用：
- Gate 3 打分：總 protein–DNA 原子級接觸數 + target strand 分項，跨 4 replica 平均，附 SD、CV 與 3.5–5.0 Å 的 cutoff 敏感度掃描。
- 正文圖：按鏈的分項長條圖（target / non-target / sgRNA），誤差棒是跨 replica SD；PAM 區另外用「按區域」的圖講，不要跟按鏈的堆在一起。
- SI：逐殘基 contact map 熱圖，圖說必須同時寫出 top-8 共享 3/8、pairwise Jaccard 平均 0.48（散布 0.23–0.78）、ρ = +0.77，並註明僅供定性參考。
- 絕對不要做的事：拿逐殘基接觸度去排變體。
09Reviewer Thinkingreviewer 會問什麼
你的 contact cutoff 是怎麼選的？換成 4.5 Å 結論還一樣嗎？點開看參考答案
這一問幾乎必問，而且是最好回答的一問 —— 只要你事先做了掃描。答案的形狀是：「主文用 4.5 Å 重原子截止（與我們可靠性分析所用的定義相同），並在 SI 給出 3.5 / 4.0 / 4.5 / 5.0 Å 四組結果；變體排序在四個 cutoff 下一致。」如果你沒做，就會被逼著承認那個數字是慣例而不是論證，整張圖的可信度就降一級。順帶一提，這跟你在 Bader 分析裡要說明分割方式、或在鄰居列表裡要說明截止半徑是同一種責任。
接觸數跟結合自由能是什麼關係？你憑什麼用它排序？點開看參考答案
不要硬拗成能量。正確的回答是把它定位成幾何 engagement 的可重現代理，而不是 ΔG 的估計。可以講的是：它與 BSA 高度一致、CV 只有 4.1%、而且可分解到股別；不能講的是它等價於結合強度。如果 reviewer 追問「那你為什麼不做 MM-PBSA」，答案是 MM-PBSA 在這個系統上又貴又不穩，而且它排的是平衡 ΔG —— 而高保真變體的專一性差別根本不在平衡 ΔG 上。
你的逐殘基 contact map 在獨立軌跡之間可重現嗎？點開看參考答案
這是致命一問，而你有一個很強的答案：不可重現，我們知道，所以我們沒有用它下結論。 給出 top-8 共享 3/8、Jaccard 0.48、ρ = +0.77，說明排序趨勢可重現但成員身分不可重現，因此所有量化結論都建立在聚合量上，逐殘基圖僅作為定性地圖放在 SI。這種答法會把一個弱點變成方法學嚴謹的證據。
你用的是原子對還是殘基對？大殘基會不會主導你的數字？點開看參考答案
必須明講，而且最好兩種都報。原子對計數比較平滑、對 cutoff 不敏感，適合當打分軸，但 Arg、Lys、Trp 這類大側鏈天生貢獻多；殘基對布林不受側鏈大小影響，但在門檻附近跳動劇烈。若你要比較 WT 與把 Arg 換成 Ala 的變體，這個區別會直接影響結論方向 —— 側鏈變小本來就會少掉一堆原子對接觸，那不見得是 engagement 真的變弱。
為什麼 non-target strand 的接觸比 target strand 不穩？這是物理還是取樣不足？點開看參考答案
兩者都有，而且要誠實區分。物理上，non-target 股在 R-loop 裡被擠開成單股、沒有配對約束，構形空間本來就大得多；取樣上，這意味著它需要更長的軌跡或更多 replica 才能收斂。你的數字（target 9.4% vs non-target 17.7%，groove 接觸計數 145.5%）支持前者為主，但不能排除後者。可防守的說法是把兩者分開報，並明說 non-target 的量目前只達 marginal tier。
10Common Mistakes最常犯的錯
- 報一張漂亮的逐殘基 contact map，卻沒有任何跨 replica 一致性數字。 後果：整張圖在統計上是裝飾品，reviewer 一問就垮。正解：圖說裡直接寫 top-N 重疊、Jaccard 與 Spearman ρ。
- 只給一個 cutoff。 後果：你的結論可能只是那個 cutoff 的產物，而你自己不知道。正解：掃 3.5 到 5.0 Å，把「結論 vs cutoff」放 SI。
- 把原子對計數與殘基對布林混著講。 後果：數量級差好幾倍，讀者無法判斷你在講什麼，也無法重現。正解：methods 裡逐字定義，兩種都報時用不同符號。
- 佔據率門檻沒寫。 後果：「穩定接觸有 23 個」這句話沒有資訊量。正解：寫出門檻，並說明換門檻結論是否翻轉。
- 用預測結構起始卻不做介面 QC。 後果：你以為在比較變體，其實在比較結構預測的猜測。正解：先看介面區的 pLDDT，低信心區（你的 HNH 62.9、target DNA 64 到 65）的接觸結論一律降級處理。
- 把 target 與 non-target 合併成一個「protein–DNA 接觸」就完事。 後果：漂亮的 4.1% 掩蓋了 non-target 的 17.7%，而 grip 變體的效應恰恰在 non-target 上。正解：永遠分股報。
11Further Reading讀哪幾篇
12Summary三句話
Contact map 把介面拆成「誰碰到誰、碰了多久」，聚合起來是這個 project 目前最可靠也最便宜的介面觀測量。在你的 5F9R 四條 replica 上，總 protein–DNA 接觸 CV 只有 4.1%、target-DNA 9.4%，但同一份資料的逐殘基 top-8 只共享 3 個、pairwise Jaccard 0.48。所以規則簡單到可以貼在螢幕上：接觸數當聚合量報，逐殘基只當地圖看。