RMSFRoot-Mean-Square Fluctuation
每個殘基在模擬裡晃多大,一條沿著序列的柔軟度曲線。
01一句話只有一句
每個殘基在模擬過程中平均晃動多大。
02Why should I care?我為什麼要讀這頁
因為這是你第一個能親眼看到 HNH 檢查點的量,而且它幾乎不用錢。
軌跡跑完之後,RMSF 是不需要任何額外算力就能拿到的結果。你在自己的系統上算出來的東西是:apo 蛋白裡 HNH 的 RMSF 是 1.86 Å,全域最高,其他 domain 只有 0.8 到 1.3 Å;換成完整三元複合體之後 HNH 仍然是最高的 1.61 Å。這條曲線把「Cas9 有一個會翻轉的催化 domain」從課本上的一句話,變成你自己機器上的一個數字。
用你熟悉的語言講:RMSF 就是均方位移 ⟨u²⟩ 的平方根。你在 DFT 裡要拿這個量,得先算 phonon、在諧近似下由 phonon DOS 積出 Debye-Waller factor。MD 直接對一個 300 K、有水、非諧的系綜取樣,省掉諧近似,代價是誤差來源從「ENCUT 與 k-mesh 收斂」換成「軌跡長度與 replica 數」。這個對應關係值得記住,因為它決定了你該對這個數字抱多少信心。
03What這是什麼
對第 i 個原子(通常取 Cα):
⟨·⟩t 是對時間平均,⟨ri⟩ 是該原子的平均位置。單位通常是 Å 或 nm。
三件事決定了這個數字的意義:
- 它對時間平均、對原子逐一保留。 這是它跟 RMSD 的唯一差別:同一條軌跡,RMSD 把原子壓掉留下時間,RMSF 把時間壓掉留下原子。
- 它必須在移除整體平移與轉動之後才有意義。 沒有先做結構疊合(superposition),整個分子在水盒裡翻滾的貢獻會被算進每一個殘基。
- 它與晶體學的 B-factor 是同一個物理量的兩種寫法:Bi = (8π²/3) · RMSFi²(RMSF 用 Å,B 得 Ų)。這給了你一個免費的 sanity check:把 RMSF 換算成 B-factor,跟晶體結構的 B-factor 欄位比對趨勢。
04Why為什麼重要
因為 Cas9 的專一性是構形檢查點問題,而檢查點的前提是那個 domain 得能動。
Dagdas 等人用單分子 FRET 證明 Cas9 在「綁上 DNA」與「切下去」之間有一個構形檢查點:HNH 必須翻轉並 docking 到催化構形,REC lobe 會先感測 R-loop 遠端的配對品質才放行。RMSF 是你在自己的軌跡裡看到這件事的最直接方式。如果 HNH 在你的模擬裡是全域最軟的區域,那你的系統至少在定性上抓到了正確的物理。
更有價值的是它作為跨方法交叉驗證的角色。你的專案裡有兩個完全獨立、而且都有數字的方法同時指向 HNH:Chai-1 預測的低 pLDDT(HNH 只有 62.9,L1 linker 765 到 809 更低到 46)與 MD 的高 RMSF(apo 1.86 Å、完整三元 1.61 Å,皆為全域最高)。兩者的失敗模式完全不同卻收斂到同一個結論,這比任何單一方法的高精度都更值得寫進 paper。你甚至量到了逐殘基 RMSF 與 pLDDT 的反相關(ρ = −0.36):結構預測「不確定」的地方,就是分子真的在動的地方。
ESM-2 的序列保守性常被當成第三個獨立訊號,但目前還沒有「ESM-2 指向 HNH」的量化結果——現有的掃描只涵蓋一個小型變體 panel 與一段不含 HNH 的逐位區間,另有一份全序列飽和掃描尚未依 domain 邊界聚合。要補齊很便宜,也不需要 GPU:把飽和掃描依 domain 邊界聚合成逐 domain 平均分數就好。做完之前,請寫「兩個方法」。
RMSF 在 domain 尺度上是可靠的訊號,在逐殘基尺度上通常是雜訊。這兩件事必須分開講,混在一起就是本頁所有錯誤的來源。
05How怎麼做
演算法只有三步:
- 選 fit selection,把每一幀疊合到參考結構。 這一步決定了答案。用整個蛋白疊合,一個大幅擺動的 domain 會顯示成高 RMSF(含它的剛體運動);只用那個 domain 自己疊合,剛體運動被扣掉,剩下的是它的內部柔軟度。兩者都對,但回答的是不同問題,而且必須在方法段寫清楚。
- 選原子。 逐殘基報告通常用 Cα;要含側鏈就用該殘基所有重原子的平均。核酸沒有 Cα,慣例是用 P 或 C1',選哪個要寫出來。
- 對時間取平均。 先丟掉平衡段(equilibration),只用 production 的部分。
工具(都是幾秒鐘的事):
gmx rmsf -f traj.xtc -s topol.tpr -res -o rmsf.xvg,-res直接給逐殘基,-fit控制是否先疊合。- MDAnalysis 的
MDAnalysis.analysis.rms.RMSF。它不會幫你對齊,官方文件明確要求先 align 到平均結構,否則結果沒有意義。 - mdtraj 有
md.rmsf,但請先確認你那個版本對reference與superpose的處理,不同版本行為不完全一致。
成本: 對 500k 原子、87 ns 的軌跡,只讀 solute 的話是分鐘等級。相對於跑出那條軌跡的 GPU 時數,RMSF 幾乎是免費的。
報告方式(這部分決定 reviewer 買不買單):
- 一定要跑多條獨立 replica,報 mean ± SD 或 CV,不要只給一條曲線。
- 在圖上標出 domain 邊界,並且分開列出 domain 級的平均值。
- 明確標出哪些區段是建模補出來的。
- 寫出 fit selection、原子選擇、丟掉多少 equilibration。
06When什麼時候用
07When NOT什麼時候別用
- 沒有先移除整體平移與轉動時,不要用。 分子在水盒裡的翻滾會被算成每個殘基的「柔軟度」,整條曲線一起被抬高。這個錯誤的特徵是:連核心區域都有異常高的 RMSF,而且數值隨軌跡長度單調上升。
- 只有一條軌跡時,不要報逐殘基結論。 你自己的數據就是最好的反例:5F9R 完整三元的 DNA RMSF 跨 4 條獨立 replica 的 CV 是 59.1%。同一個系統、同樣的設定、只換亂數種子,逐殘基數字可以差到將近一倍。任何從單條軌跡讀出「第 X 號殘基特別柔軟」的敘述,在統計上都沒有支撐。
- 建模補出來的區段與鏈端,不要放進比較。 你的系統裡補上去的 5′ 單股末端(res 1 到 11)RMSF 有 6 到 8 Å,而晶體本體只有 3.9 Å。這些區段沒有實驗座標約束,它們的高 RMSF 反映的是「這裡本來就沒資料」,不是生物學。把它們留在圖裡會汙染平均值,也會讓 reviewer 直接懷疑整張圖。
- 不要把 RMSF 高低當成功能重要性的代理。 催化中心通常是剛性的,柔軟的往往是表面 loop。「這個殘基 RMSF 最高所以最重要」是這個領域最常見的過度解讀之一。RMSF 告訴你哪裡在動,不告訴你哪裡重要。
- 不要跨系統比較絕對值。 不同的軌跡長度、力場、水模型、溫度、fit selection 都會系統性地移動 RMSF。跨 paper 比較絕對數字幾乎一定是錯的;能比的只有同一套設定下的相對排序。
- 不要把 RMSF 當成 ranking 分數。 它沒有跟任何實驗 observable 校準過。要拿來排序變體,你得先證明它和你真正在乎的東西(cleavage rate、specificity ratio)有相關性,而這件事目前沒人做到。
在 2 條 replica 的時候,逐殘基的一致性看起來很漂亮;跑到 4 條就掉下來了。小 n 會系統性地高估一致性。 這不是 RMSF 特有的問題,但 RMSF 因為便宜、看起來直觀,特別容易讓人在 n=1 或 n=2 就下結論。詳見取樣與收斂。
08Project Lens跟我的 project 多相關
為什麼是四星而不是五星: 因為它已經被自己的數據判了一半死刑。逐殘基與 DNA 的 RMSF 在 4 條 replica 之間 CV 高達 59.1%,屬於明確的 noise tier,不可能拿來當 Gate 3 的打分軸。任何以逐殘基 RMSF 排序變體的設計都會被自己的可靠性分析打臉。
為什麼還有四星:
- 它是跨方法交會處的一角。 Chai-1 pLDDT 與 MD RMSF 兩個獨立方法同時指向 HNH,這個交會本身就是可以寫進 paper 的論點,而 RMSF 是其中最物理、最容易對 reviewer 解釋的一個。(ESM-2 那條腿的分析還沒做,見上方補充框。)
- 成本近乎為零。 在「便宜」是硬約束的 workflow 裡,一個軌跡跑完就順手拿到、不用額外 GPU 時數的 metric,就算只能當 QC 也值得保留。
- 它是必要的健康檢查。 每一條 production 軌跡都應該先看 RMSF 再看別的東西,這是防止把建系統的錯誤一路帶到結論的最後一道便宜關卡。
具體怎麼用: 放在 workflow 的 QC 層與敘事層,不放在打分層。報告時只給 domain 級聚合值加跨 replica 誤差棒,逐殘基曲線可以畫但要標明僅供定性參考。要打分,用介面埋藏面積或聚合接觸數那類 CV 在 10% 以下的量。
09Reviewer Thinkingreviewer 會問什麼
你跑了幾條獨立軌跡?逐殘基 RMSF 的統計誤差是多少?點開看參考答案
這是最致命的一問,因為多數 Cas9 MD 論文只跑一到兩條。誠實的回答是給出跨 replica 的 CV,並且明說哪些尺度可信、哪些不可信。你的數據可以直接回答:聚合量 CV 在 2 到 10%,逐殘基與 DNA RMSF 在 24 到 146%,所以結論只建立在前者上。先自己說出限制,比被問出來好。
你的 fit selection 是什麼?有沒有把 domain 的剛體運動誤算成柔軟度?點開看參考答案
如果你用整個蛋白疊合,那 HNH 的高 RMSF 有相當一部分是它整塊在擺動,而不是它內部變柔軟。兩種說法在生物學上意義不同。正確做法是兩種都報:全域疊合看整體位移,domain 內疊合看內部柔軟度,並說明結論建立在哪一個上。
你的 RMSF 高峰是不是都落在鏈端與你自己補上去的區段?點開看參考答案
這一問幾乎可以廢掉一張圖。務必在圖上用不同顏色標出所有無實驗座標的區段,並且在計算 domain 平均時把它們排除。你的 5′ 端 res 1 到 11(6 到 8 Å,本體 3.9 Å)就是典型例子。
RMSF 和實驗 B-factor 對得起來嗎?點開看參考答案
可以用 B = (8π²/3)·RMSF² 換算後比對趨勢。要注意兩件事:晶體 B-factor 含有靜態無序與晶格效應,數值不會一對一;而且晶體堆積會壓制表面 loop 的運動,所以模擬的表面 RMSF 系統性偏高是正常的。比趨勢(Spearman 相關)不比絕對值。
你怎麼確定高 RMSF 不是模擬還沒收斂?點開看參考答案
分段檢查:把 production 切成前半後半分別算 RMSF,看有沒有系統性漂移。如果後半明顯高於前半,那是還在鬆弛,不是平衡態的柔軟度。這個檢查很便宜,但幾乎沒人報。
10Common Mistakes最常犯的錯
- 只跑一條軌跡就報逐殘基 RMSF。 後果:讀者無從判斷哪些峰是真的。正解:至少 3 到 4 條獨立 replica,報 mean ± SD,並公開 CV。
- 忘記移除整體轉動。 後果:整條曲線被抬高,核心區也「柔軟」。正解:先 superpose;檢查特徵是核心 β-sheet 是否異常高。
- 把建模補出的 loop 和鏈端留在圖裡與平均值裡。 後果:少數殘基把 domain 平均值拉走,結論失真。正解:標色排除,並在圖說寫明排除了哪些殘基。
- 比較 WT 與變體時用了不同的 fit selection。 後果:差異全部來自對齊方式,不是突變。正解:固定 fit selection,並在補充材料放兩者的對齊定義。
- 由「RMSF 最高」直接推論「功能最重要」。 後果:邏輯跳躍,reviewer 一眼看穿。正解:只講「這個區域在模擬中最具構形自由度」,功能推論要另外的證據支持。
- 跨 paper 比較 RMSF 絕對值。 後果:力場、長度、溫度都不同,數字不可比。正解:只比同一套設定下的相對排序,或自己重跑對照組。
11Further Reading讀哪幾篇
12Summary三句話
RMSF 是把一條軌跡對時間壓平、保留每個殘基柔軟度的曲線,它與晶體 B-factor 是同一個物理量。它在 domain 尺度上可靠、成本近乎為零,是你在自己機器上看見 HNH 構形檢查點的最便宜方式。但逐殘基與核酸的 RMSF 在獨立軌跡之間可以差到六成,所以它適合當 QC 與敘事,不適合當排序變體的打分軸。