Mercurial > repos > chemteam > mdanalysis_distance

comparison pca_cosine.py @ 3:489b25966bb9 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/galaxycomputationalchemistry/galaxy-tools-compchem/ commit 3ff06e3182c3a1546ea0a3b29e0d4383e12169e1
author chemteam
date Wed, 03 Apr 2019 15:46:56 -0400
parents
children 312f912de69d
comparison
equal deleted inserted replaced
2:46892d756cec 3:489b25966bb9
1 #!/usr/bin/env python
2
3 import argparse
4 import csv
5 import sys
6
7 import MDAnalysis as mda
8 import MDAnalysis.analysis.pca as pca
9
10 import numpy as np
11
12
13 def parse_command_line(argv):
14 parser = argparse.ArgumentParser()
15 parser.add_argument('--idcd', help='input dcd')
16 parser.add_argument('--ipdb', help='input pdb')
17 parser.add_argument('--icomponents', help='number of principle components')
18 parser.add_argument('--iindex', help='index of the PC')
19 parser.add_argument('--output', help='output')
20 parser.add_argument('--cosout', help='cosine output')
21 return parser.parse_args()
22
23
24 args = parse_command_line(sys.argv)
25
26 u = mda.Universe(args.ipdb, args.idcd, topology_format="PDB", format="DCD")
27
28 components = int(args.icomponents)
29 pca_index = int(args.iindex)
30
31 PSF_pca = pca.PCA(u, select='backbone')
32 PSF_pca.run()
33 n_pcs = np.where(PSF_pca.cumulated_variance > 0.95)[0][0]
34 atomgroup = u.select_atoms('backbone')
35
36 pca_space = PSF_pca.transform(atomgroup, n_components=components)
37 cosine = mda.analysis.pca.cosine_content(pca_space, pca_index)
38
39 PCA = list(pca_space)
40
41 with open(args.output, 'w') as f:
42 writer = csv.writer(f, delimiter='\t')
43 writer.writerows(PCA)
44
45 with open(args.cosout, 'w') as f1:
46 f1.write(str(cosine))