PCAPrincipal component analysis of trajectories
把幾萬維的運動壓成兩三個「主要在做什麼」的方向。
01一句話只有一句
把上萬維的原子運動壓成幾個主要的集體方向。
02Why should I care?我為什麼要讀這頁
因為這件事你早就會了，只是換了一個矩陣。
在 DFT 裡你對角化 Hessian（力常數矩陣），拿到 normal mode 與 phonon：特徵值給你頻率，特徵向量告訴你哪些原子一起動。PCA 做的是同一個動作，只是它對角化的不是二次微分，而是軌跡的位移共變異數矩陣。用一句話記住差別：phonon 是解出來的，principal component 是取樣出來的。 這一個字，是這頁後半所有警告的唯一來源。
實務上的價值很直接。你的 5F9R 三元複合體有 535,265 個原子；就算只留蛋白的 Cα，SpCas9 常引用的長度是 1368 個殘基，也還有 3 × 1368 = 4104 個維度。你不可能盯著 4104 條曲線找故事。PCA 讓你問一個能回答的問題：在這 87 ns 裡，這台機器主要在做哪一件事？ 如果答案是「HNH 相對 REC lobe 的開合」，你就拿到一個可以畫自由能面、可以拿來比 WT 與變體、可以當 enhanced sampling collective variable 的量。如果答案是「建模補出來的 DNA 末端在甩尾」，你也學到了很重要的事，只是那件事跟 Cas9 無關。
03What這是什麼
把對齊過的每一幀座標攤平成一個 3N 維向量 x(t)，共變異數矩陣是
i、j 跑遍 3N 個 Cartesian 分量，所以 C 是 3N × 3N 的實對稱矩陣。對角化它：
- 特徵向量 vk 就是第 k 個 principal component：一個 3N 維方向，意思是「所有原子按這個比例一起位移」的一種集體運動。
- 特徵值 λk 是軌跡投影到該方向上的變異數，單位是長度平方，習慣由大到小排。
- 投影 pk(t) = vk · ( x(t) − ⟨x⟩ ) 是一條一維時間序列。把 p1 對 p2 畫成散佈圖或密度圖，就是那張到處都看得到的 PC1-PC2 構形分佈圖。
三個必須先記住的性質：
- 所有特徵值的和，等於所有原子的均方位移總和。 trace(C) = Σk λk = Σi ⟨Δxi²⟩ = Σi RMSFi²。也就是說 PCA 沒有創造任何新資訊，它只是把 RMSF 的總量從「分配到原子上」改成「分配到方向上」。 這是理解 PCA 最快的一句話，也解釋了為什麼 RMSF 的毛病 PCA 大多會繼承。
- 對齊之後會冒出 6 個接近零的特徵值。 3 個平移加 3 個轉動被 fitting 扣掉了 —— 正是你在分子 normal mode 計算裡要丟掉的那 6 個零頻模態，或晶體在 Γ 點趨於零的 3 支 acoustic mode。看到它們表示對齊做對了，看不到就是你根本沒對齊。但這個檢查只有在「疊合用的原子」與「進矩陣的原子」是同一組時才乾淨：如果你用蛋白疊合卻把 DNA 一起放進共變異數矩陣，那 6 個值不會壓到零，這時只能看它們有沒有明顯小一截。
- 前幾個 PC 佔多少百分比的變異，是取樣的性質，不是系統的性質。 軌跡越長、走過的構形越多，總變異就會被重新分配。所以「前兩個 PC 佔 XX%」這個數字不能跨軌跡比較，更不能拿來當收斂的證據 —— 事實上取樣越差，這個百分比通常看起來越漂亮。
跟 Hessian 的量化關係。 在諧近似下，質量加權座標的均方位移是 ⟨qk²⟩ = kBT / ωk²，所以
這條式子有一個前提，很容易被跳過： 上面定義的 C 是純 Cartesian 的，要換到頻率必須先把座標乘上 √m 做成質量加權版本再對角化，否則 λk 跟 ωk 對不起來。多數軌跡分析工具預設不做質量加權（gmx covar 要加 -mwa），所以只要你打算報頻率或熵，就得先確認這一步。只做構形分析、不換算頻率的話，不加權完全可以。
共變異數矩陣的大特徵值 = 低頻 = 軟模，這就是 MD 文獻裡 quasi-harmonic analysis 的核心。但要注意它跟 Hessian 是反過來的： Hessian／動力學矩陣的特徵值是 ω²，大特徵值是高頻硬模；共變異數矩陣的特徵值是擺動幅度的平方，大特徵值才是軟模。兩者在諧近似下互為倒數（H ≈ kBT C⁻¹），所以看到「最大的特徵值」時，第一件事永遠是先確認你手上是哪一個矩陣。這也解釋了為什麼 PC1 通常長得像 domain 開合而不是鍵長振動：跟晶格裡頻率最低的那支聲學支對應到最大幅的原子位移，是完全同一件事。
用取樣出來的共變異數矩陣反推有效簡正模態與構形熵，MD 文獻習慣叫 quasi-harmonic analysis。這跟固態／DFT 語境裡的 quasi-harmonic approximation（QHA）不是同一件事。 QHA 是在好幾個不同體積各做一次諧性聲子計算、讓頻率隨體積變化，用來取熱膨脹與有限溫度自由能；這裡的 quasi-harmonic 只有一次對角化，非諧與溶劑效應是靠「軌跡實際走過的分布」隱含帶進來的。名字撞在一起，物理內容不同，寫方法段時不要混用。
PCA 是「用取樣代替 Hessian」的 normal mode 分析。好處是它天生含非諧、含溶劑、在真實溫度下；壞處是它只知道你走過的地方，而 Hessian 不需要你走到任何地方。
04Why為什麼重要
因為 Cas9 的核心問題發生在 domain 尺度，而 PCA 是少數天生就在這個尺度上工作、又不需要你事先指定「要看哪裡」的分析。
HNH 從檢查點位置翻轉並 docking 到催化構形，是一個涉及數百個殘基協同位移的大幅運動。逐殘基的量（RMSF、逐位點距離）在你自己的可靠性分析裡已經被判進 noise tier：DNA RMSF 跨 4 條 replica 的 CV 是 59.1%，逐位點距離 38 到 46%。而聚合量是可信的：BSA 2.1%、總 protein–DNA 接觸 4.1%、鹽橋 9.5%。PCA 在構造上就是一個聚合量 —— 每一個 PC 都是全體原子的加權組合，投影值 pk(t) 是幾千個原子位移的加權和，個別原子的雜訊在這個加權裡會被平均掉。這是它值得留在 workflow 裡的最強理由。
但這個好處只在「基底已經固定」的前提下成立。 基底本身是從同一個有雜訊的矩陣估出來的：兩個特徵值靠得近的時候，對應的特徵向量會在它們張出的平面裡任意旋轉，換一條 replica 就換一組模態。所以「投影量是聚合的」與「模態方向是穩的」是兩件事，前者成立不保證後者成立 —— 這就是本頁 When NOT 第一條與 Project Lens 只給三星的根源。
第二個理由是它是通往其他所有東西的入口。你之後要做的 free energy landscape、Markov state model、metadynamics、umbrella sampling，第一件事都是「選一組 collective variable」。選錯 CV 的代價極高：你會花掉 GPU 時數去加速一個跟你在乎的事情正交的方向。PCA 是產生 CV 候選最便宜的方式。
但要注意 enhanced sampling 不是只有這一條路，這是一個真正的分岔。 Palermo 等人 2017 年處理 Cas9 構形活化用的是 GaMD，那是不需要事先指定 CV 的方法：它對整體位能面加一個 boost 再重加權回去。他們的做法反而是先用 targeted MD 造出一條路徑、沿路徑取 24 個等距 frame 當起點各自撒種子，累積約 15 μs。所以真正要問的是：你有沒有把握選對 CV。 有把握就走 CV-based（metadynamics、umbrella sampling），PCA 是產生候選的工具；沒把握就走 CV-free（GaMD），代價是 boost 會扭曲原始能量面、重加權後的自由能仍帶系統誤差，而且拿不到速率常數。
第三個理由比較少人講：PCA 是一個誠實的取樣體檢儀。 因為它的答案完全由軌跡走過的地方決定，所以反過來，它的不穩定程度就直接量到了你的取樣有多不足。把 4 條 replica 各自做一次 PCA、比對前 10 個特徵向量張出的子空間有多重疊，你就得到一個關於「87 ns 到底夠不夠」的定量答案 —— 而這件事你目前沒有別的便宜方法可以回答。
05How怎麼做
五個步驟，前兩步決定成敗：
- 對齊。 移除整體平移與轉動，跟 RMSF 完全一樣的要求，而且錯的後果更嚴重：沒對齊的話 PC1 幾乎一定是整個分子在水盒裡翻滾。fit selection 要寫進方法段。
- 選原子。 逐殘基故事用 Cα（蛋白）加 P 或 C1′（核酸）。要不要把核酸放進來是一個真實的決策點：核酸柔軟度遠高於蛋白核心，放進來之後 PC1 常常被 DNA 末端主宰。建議兩種都算，並在報告裡說明你選哪一種、為什麼。
- 組矩陣並對角化。 4104 × 4104 的實對稱矩陣在筆電上是秒級的事，瓶頸永遠在讀軌跡不在線性代數。
- 投影。 把每一幀投到 PC1、PC2，畫散佈圖或 2D 直方圖。
- 看模態長什麼樣。 產生沿著 PC1 的極端構形做成動畫。這一步不可省 —— 一個 PC 只有在你看過它實際在動什麼之後才有物理意義。
工具：
# GROMACS：covar 預設 -fit yes，會自己把每一幀疊到 -s 給的結構，
# 所以這裡不必先跑 trjconv；要質量加權才加 -mwa（預設 no）
gmx covar -f prod.xtc -s prod.tpr -n ca.ndx -o eigenval.xvg -v eigenvec.trr
# 投影到 PC1、PC2
gmx anaeig -v eigenvec.trr -f prod.xtc -s prod.tpr -first 1 -last 2 -2d proj_1_2.xvg
# 沿 PC1 的極端構形，丟進 VMD 看它到底在動什麼
gmx anaeig -v eigenvec.trr -f prod.xtc -s prod.tpr -first 1 -last 1 -extr pc1.pdb -nframes 20
# MDAnalysis
from MDAnalysis.analysis import align, pca
# 注意參考結構：這樣寫是疊到第 0 幀。理論上該疊到平均結構
# （先算平均、再對它疊一次），第 0 幀只是「剛好那一幀」
align.AlignTraj(u, u, select="name CA", in_memory=True).run()
p = pca.PCA(u, select="name CA").run()
p.results.cumulated_variance[:10] # 累積變異
proj = p.transform(u.select_atoms("name CA"), n_components=3)
這裡有一個工具間不一致的坑： gmx covar 預設會替你對齊，MDAnalysis 的 pca.PCA 預設 align=False 不會，跟 rms.RMSF 是同一個坑。換工具的時候務必回頭確認這一項，不然你會拿到兩張完全不同的圖卻以為是物理差異。R 使用者用 Bio3D 的 pca.xyz()。
跨系統與跨 replica 的正確做法（這段最常被做錯）： 每個系統各自對角化出來的基底不一樣，A 的 PC1 跟 B 的 PC1 沒有可比性。要比較就必須建立共同子空間：把要比的軌跡串接起來做一次對角化拿到共同基底，再把每條軌跡分別投影上去。這樣畫出來的 PC1-PC2 圖才是同一個座標系。
收斂診斷（成本近乎為零，但幾乎沒人報）：
- Cosine content。 隨機擴散的 principal component 是純餘弦函數，週期數等於 PC 序號的一半 —— GROMACS 手冊把這件事直接寫在
gmx anaeig -proj的說明裡。所以一個「看起來像從 A 狀態走到 B 狀態」的 PC1 時間序列，可能只是布朗運動。這個量從 0（完全不像餘弦）到 1（完美餘弦），越接近 1 越無法與隨機擴散區分。指令：先gmx anaeig ... -first 1 -last 3 -proj proj.xvg，再gmx analyze -f proj.xvg -cc cosine.xvg。它沒有公認的閾值，所以正確做法是把 PC1 到 PC3 的數值連同軌跡長度一起報，並且對 cosine content 偏高的那幾個 PC 不下物理結論。 - RMSIP（root mean square inner product）。 拿兩組各前 10 個特徵向量，算它們張出的子空間有多重疊，1 是完全相同、0 是完全正交。至少做兩種比較：同一條軌跡的前半對後半，以及 replica 兩兩之間。文獻上流通的粗略基準在 0.7 上下，但這是慣例不是定理，不同作者用的門檻不一樣、取前幾個特徵向量也會改變答案 —— 所以報數字（連同你取了幾個模態）比報「通過」重要得多。
報告時要寫的東西： fit selection、原子選擇、丟掉多少 equilibration、累積變異曲線、cosine content、跨 replica 與前後半的 RMSIP、以及 PC1 的極端構形圖。少了最後三項，reviewer 有充分理由認為你只是在展示一條軌跡的偶然。
06When什麼時候用
- 看完 RMSD 與 RMSF 之後的第三個動作，用來回答「這條軌跡主要在做什麼」。
- 要挑 collective variable 給 free energy landscape、metadynamics、umbrella sampling 或 MSM 時，這是最便宜的候選產生器。
- 要把 domain 級的運動做成動畫給人看時，
-extr出來的極端構形比任何文字都有效。 - 要診斷取樣夠不夠時。 用 cosine content 與 RMSIP 把 PCA 當儀器而不是當結論，這是它在你現階段最誠實、最可防守的用途。
- 要在共同子空間裡比較 WT 與變體的構形分佈時 —— 前提是合併對角化、取樣足夠、而且你已經先報過 RMSIP。
07When NOT什麼時候別用
- 取樣不足時，不要把 PC 當成系統的柔性模態。 這是本頁最重要的一條。PCA 只知道軌跡走過的地方，取樣不足時前幾個 PC 反映的是「這條軌跡剛好走了哪裡」，換一顆亂數種子就換一組模態。量級推估一下你的處境：你有 4 × 87 ns，假設每 10 ps 存一幀就是每條 8700 幀，看起來比 4104 個維度多；但 MD 的幀之間高度相關，如果 domain 級運動的相關時間是 ~1 ns 的量級，那 87 ns 只提供了 10¹ 到 10² 量級的獨立樣本，去估一個 4104 × 4104 的共變異數矩陣。怎麼看出它壞了： 算 PC1 到 PC3 的 cosine content，接近 1 就是隨機擴散；再把軌跡切前後半各做一次 PCA，RMSIP 掉下來就代表模態還沒定下來。
- 不要跨系統直接比 PC。 WT 的 PC1 與變體的 PC1 是兩個不同對角化的產物，把它們畫在同一張圖上、比較投影範圍或變異百分比，在數學上沒有定義。怎麼看出它壞了： 把兩邊 PC1 的極端構形動畫並排播放，如果一邊在動 HNH、另一邊在動 DNA 末端，那張比較圖就是垃圾。正解： 合併軌跡做一次對角化取得共同基底，或把 B 投影到 A 的特徵向量上，並且報 RMSIP。
- PC1 變異大，不代表有功能意義。 變異大只代表軟，而最軟的通常是表面 loop、鏈端與建模補出來的區段。你的系統補出來的 5′ 單股末端 RMSF 有 6 到 8 Å，晶體本體只有 3.9 Å。RMSF 是平方進 trace(C) 的，所以這十幾個殘基對總變異的貢獻遠高於它們的數量佔比；更關鍵的是它們同調地甩，那份變異會集中在少數幾個方向上 —— 這正是主宰 PC1 的配方。怎麼看出它壞了： 看
-extr的動畫。正解： 排除無實驗座標的區段重算，並且兩種結果都報。 - 不要把 PC1-PC2 的密度圖直接當自由能面讀。 −kBT ln P(p1, p2) 只有在取樣收斂時才是自由能。沒收斂的時候，它是「訪問頻率的對數」，而有限長度的軌跡投影上去很容易出現幾個看似分離的盆地 —— 那可能是取樣的坑洞，不是能量的坑洞。怎麼看出它壞了： 四條 replica 各自的密度圖如果落在不同位置，那些「盆地」就是 replica 標籤，不是狀態。細節見自由能景觀。
- 不要把 PC 投影當成排序變體的分數。 它跟 RMSF 有完全相同的問題：沒有跟任何實驗 observable 校準過，跨系統不可比，而且它的誤差棒你根本沒量過。要打分請用 CV 在 10% 以下的聚合量，例如介面埋藏面積或接觸數。
- PCA 是線性方法，不要指望它抓到彎曲的反應路徑。 如果真實的構形變化在高維空間裡是一條曲線，把它投影到兩個直線方向上，路徑的兩端可能被壓在同一個位置。怎麼看出它壞了： PC1-PC2 圖上兩個「不同狀態」的結構其實 RMSD 很大卻投影重疊。要處理這件事得換工具（tICA、擴散映射，或直接做 MSM）。
軌跡越短，走過的構形越少，PC1 佔的變異百分比就越高、投影圖上的「狀態」就越乾淨。視覺說服力與統計可靠度在這裡是反向的。 看到一張四個盆地分得清清楚楚的 PC1-PC2 圖，第一個該問的不是「這四個狀態是什麼」，而是「你跑了幾條軌跡、每條的投影分別落在哪」。
08Project Lens跟我的 project 多相關
為什麼不是更低： 因為它是這個 project 目前少數同時滿足「近乎免費」與「天生聚合」兩個條件的分析。你的可靠性分析已經把逐殘基與逐位點的量判進 noise tier（DNA RMSF CV 59.1%、逐位點距離 38 到 46%、groove 接觸計數 145.5%），而 PCA 在構造上就是把所有原子加權平均成一個方向，屬於聚合量那一族。它也是通往 free energy landscape、MSM 與 CV-based enhanced sampling 的入口 —— 不做 PCA，你連 collective variable 的候選都生不出來（改走 CV-free 的 GaMD 可以繞開這一步，但那是另一組成本與另一組系統誤差）。加上它的成本只有幾秒 CPU，在「便宜」是硬約束的專案裡，這種 metric 沒有理由不留著。
為什麼不是四星或五星： 因為到目前為止，你還沒有量過 PCA 在你自己系統上的跨 replica 可重現性（你那份 4-replica 可靠性分析的清單裡沒有 PCA 這一項，目前尚未量測）。而先驗上不樂觀：PCA 屬於「對取樣極度敏感」那一族，跟已經被判為 noise 的 interface iRMSD（23.9%）是近親。在沒有 RMSIP 數字之前，任何以 PC 為基礎的物理主張都缺一個必要的支撐。它也不可能進 Gate 3 的打分層 —— 跨系統不可比這件事是結構性的，不是取樣多一點就能解決的。
具體怎麼用： 分成兩層。診斷層（現在就做，優先）：對現有 4 條 replica 算合併 PCA，報 cosine content 與兩兩 RMSIP，用它回答「87 ns 到底夠不夠」。這件事本身就是 methodology paper 裡有價值的一段，而且成本是零。敘事層（做，但標明定性）：用 PC1 的極端構形說明 HNH 相對 REC 的開合，配合 HNH 距離這個有明確物理意義的量一起講。打分層：不放。
09Reviewer Thinkingreviewer 會問什麼
你的 principal component 收斂了嗎？cosine content 與 RMSIP 是多少？點開看參考答案
這是這一頁最致命的一問，因為多數 MD 論文的 PCA 完全沒有報這兩個數字。誠實的回答是直接給出來：PC1 到 PC3 的 cosine content、前半對後半的 RMSIP、以及 replica 兩兩之間的 RMSIP（都取前 10 個特徵向量）。如果數字不好看，正確的做法是把 PCA 的角色降級成定性描述並明說，而不是不報。先自己說出限制，永遠比被問出來好。
你的 PC1 是不是被建模補出來的區段或鏈端主宰？點開看參考答案
把 -extr 產生的 PC1 極端構形動畫放進補充材料，這一問就自動回答了。同時要報一組排除無實驗座標區段之後重算的 PCA，看主要模態是否維持。你的 5′ 端 res 1 到 11（RMSF 6 到 8 Å，本體 3.9 Å）是最明顯的候選污染源。
WT 和變體的 PC 你是怎麼比的？兩邊是同一組基底嗎？點開看參考答案
如果兩個系統各自對角化，答案就是「沒有可比性」，這一問會直接廢掉那張比較圖。正確回答是：把兩條軌跡串接後做一次對角化得到共同基底，兩邊都投影到同一組特徵向量上，並且報兩個系統各自子空間之間的 RMSIP 作為「它們是否在同一個構形空間裡活動」的量化證據。
你把 PC1-PC2 的密度圖標成自由能面，你的取樣支持這個標法嗎？點開看參考答案
只有在你能證明取樣收斂時，−kBT ln P 才是自由能。可防守的做法是把四條 replica 用不同顏色疊在同一張投影圖上：如果它們互相重疊並且各自都跨越了那些「盆地」，那個標法還有救；如果每條各佔一角，那張圖只能標成 population density，而且不能推論能障。
PCA 是線性的，你怎麼確定它沒有把兩個不同的構形狀態壓在同一個投影位置？點開看參考答案
誠實的回答是先承認這個風險存在，再給檢查：挑出投影位置相近的幀，算它們之間真實的 pairwise RMSD，如果分布是雙峰的就代表發生了投影重疊。若真的重疊，要換到能處理非線性或動力學相關性的方法（tICA、擴散映射、MSM），而不是硬用 PC 講故事。
10Common Mistakes最常犯的錯
- 沒對齊就做 PCA。 後果：PC1 是整個分子在水盒裡翻滾，所有解釋全錯。正解：先 superpose；檢查特徵：對角化後找不到那 6 個接近零的特徵值。
- 只用一條軌跡的 PCA 下物理結論。 後果：換一顆種子就換一組模態，結論不可重現。正解：至少 4 條 replica，報 RMSIP 與 cosine content。
- 拿兩個各自對角化的系統比 PC1。 後果：比較在數學上沒有定義。正解：合併軌跡做一次對角化取得共同基底。
- 報「前兩個 PC 佔 XX% 變異」當作收斂或重要性的證據。 後果：這個百分比在取樣越差時看起來越好，論證方向是反的。正解：只把它當描述性數字，收斂另外用 RMSIP 與 cosine content 講。
- 把補出來的鏈端留在原子選擇裡。 後果：末端甩尾主宰 trace(C)，PC1 變成一段沒有實驗座標的 DNA。正解：排除後重算，兩種都報。
- 只給 PC1-PC2 散佈圖，不給模態動畫。 後果：讀者無法判斷那個方向的物理意義，reviewer 只能假設最壞情況。正解：
-extr出極端構形，放進補充材料。 - 直接把密度圖標成自由能（kcal/mol）。 後果：把取樣的坑洞講成能量的坑洞。正解：先證明收斂，否則只標 population。
11Further Reading讀哪幾篇
12Summary三句話
PCA 把軌跡的共變異數矩陣對角化，特徵向量是集體運動方向、特徵值是該方向上的變異數，數學上就是把 Hessian 換成取樣的 quasi-harmonic 版 normal mode 分析。它天生是聚合量、成本只有幾秒 CPU，適合把 domain 級的故事講清楚，也是挑 collective variable、通往自由能面與 MSM 最便宜的入口。但它的答案完全由軌跡走過的地方決定，所以在 4 × 87 ns 這個尺度上，它更該被當成診斷取樣的儀器，而不是產生結論的工具。