7. Principal Component Analysis#

  • example with PyTraj

 1#!/usr/bin/env python
 2import os
 3from glob import glob
 4import pytraj as pt
 5import numpy as np
 6import matplotlib.pyplot as plt
 7
 8analysis="combine-all"
 9prm="../input/step3_pbcsetup.parm7"
10crd="combine-nc/all.nc"
11# Atom mask selection
12ambermask='@CA'
13
14os.makedirs('img', exist_ok=True)
15os.makedirs('raw_data', exist_ok=True)
16
17"""
18Pytraj template for MD trajectory analysis
19"""
20
21traj = pt.iterload([crd], top=prm)
22# RMSF Analysis (_trajectory is(are) trajectory file(s), _mask is the Amber format selection mask)
23data = pt.pca(traj, mask=ambermask, n_vecs=10)
24
25# Projection Data
26_data = data[0]
27pc1_data = data[1][0]                                       # Eiganvalues for first mode (percent)
28pc2_data = data[1][0]                                       # Eiganvalues for second mode (percent)
29
30# Percent Variance
31x_label = (pc1_data[0] / np.sum(pc1_data[:])) * 100
32y_label = (pc2_data[1] / np.sum(pc2_data[:])) * 100
33
34flip1 = 1                                                   # Flip the sign of the eigenvector if necessary (1 or -1)
35flip2 = 1                                                   # Flip the sign of the eigenvector if necessary (1 or -1)
36x_data = _data[0] * flip1
37y_data = _data[1] * flip2
38
39plt.scatter(x_data, y_data, marker='o', c=range(traj.n_frames), alpha=0.5)
40
41plt.xlabel(f"PC1 ({ str(np.round(x_label, 1)) } %)")        # Label x-axis
42plt.ylabel(f"PC2 ({ str(np.round(y_label, 1)) } %)")        # Label y-axis
43
44axis_lim = 10                                               # Set axis limit
45plt.xlim(-axis_lim, axis_lim)                               # Set x-axis limit
46plt.ylim(-axis_lim, axis_lim)                               # Set y-axis limit
47cbar = plt.colorbar()                                       # Show colorbar
48cbar.set_label('Frame Number')                              # Label colorbar
49plt.grid(linestyle='--', alpha=0.2)
50plt.savefig(f"img/pca-{ analysis }.png")                        # Save plot
Show Source