Mercurial > repos > chemteam > mdanalysis_ramachandran_protein
comparison pca_cosine.py @ 0:0f270722aca6 draft
"planemo upload for repository https://github.com/galaxycomputationalchemistry/galaxy-tools-compchem/ commit 1b23e024af45cc0999d9142d07de6897d4189ec2"
author | chemteam |
---|---|
date | Mon, 24 Aug 2020 16:27:56 -0400 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:0f270722aca6 |
---|---|
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('--itraj', help='input traj') | |
16 parser.add_argument('--istr', help='input str') | |
17 parser.add_argument('--itrajext', help='input traj ext') | |
18 parser.add_argument('--istrext', help='input str ext') | |
19 parser.add_argument('--icomponents', help='number of principle components') | |
20 parser.add_argument('--iindex', help='index of the PC') | |
21 parser.add_argument('--output', help='output') | |
22 parser.add_argument('--cosout', help='cosine output') | |
23 return parser.parse_args() | |
24 | |
25 | |
26 args = parse_command_line(sys.argv) | |
27 | |
28 u = mda.Universe(args.istr, args.itraj, | |
29 topology_format=args.istrext, format=args.itrajext) | |
30 | |
31 components = int(args.icomponents) | |
32 pca_index = int(args.iindex) | |
33 | |
34 PSF_pca = pca.PCA(u, select='backbone') | |
35 PSF_pca.run() | |
36 n_pcs = np.where(PSF_pca.cumulated_variance > 0.95)[0][0] | |
37 atomgroup = u.select_atoms('backbone') | |
38 | |
39 pca_space = PSF_pca.transform(atomgroup, n_components=components) | |
40 cosine = mda.analysis.pca.cosine_content(pca_space, pca_index) | |
41 | |
42 PCA = list(pca_space) | |
43 | |
44 with open(args.output, 'w') as f: | |
45 writer = csv.writer(f, delimiter='\t') | |
46 writer.writerows(PCA) | |
47 | |
48 with open(args.cosout, 'w') as f1: | |
49 f1.write(str(cosine)) |