# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# <img src="https://www.python.org/static/opengraph-icon-200x200.png" align="left" width="70" height="70"> <h1>Python Basics</h1>

# [Python](https://www.python.org/) is an interpreted language originally designed for general-purpose programming, but it also became a full-fledged tool for scentific computing thanks to its highly modularized and extensible design, and a growing scentific user/developer community. The first part of this tutorial is to show or remind you some basic syntax of Python.

print('hello world!')

# The real power of Python is with the packages. There are many built-in ones, just like the one we will see below, and also many open-source packages developed by the Python community. There are several ways to import a package (module). For starters, you can import the module by its name:

import numpy as np
np.pi

# Another common way to import is using a wildcard (namely `*`):

from numpy import *
pi

# ### Data types

# Just like other languages, Python has several basic data types: `int`, `float`, `bool`, `str`, to define a variable of a certain type, one can simply write out the initial value:

i = 1
f = 3.1415
b = True or False
s = 'hello world!'

# To create a list:

a = [0, 1, 2, 3, 4]
a

# Notably the array in Python starts from 0. A list can be also generated by `range` function, minimally the only input is the **length** of the list. So to generate a list like above, we can write:

a = range(5)
a

# Note that in Python 3K, this will generate an iterator (a [Range](https://docs.python.org/3/library/stdtypes.html#typesseq-range) class to be precise). Following code will convert it to a `list` instance:

a = list(range(5))
a

# *Optionally*, the most robust and universal way to convert a iterable object to a list is through [list comprehension](https://docs.python.org/3/tutorial/datastructures.html#list-comprehensions):

a = [_ for _ in range(5)]
a

# ### Control flow

# The list comprehension mentioned above is essentially an one liner of a `for` loop. The full version of the `for` loop looks like the following:

for i in range(5):
    print(i)

# Which iterates through numbers from 0 to 4. `if` block is used for making choices. Let's first initialize two integers:

a = 3; b = -3

# Then write out a control block as follows:

if a > b:
    print('a is bigger than b')
elif a < b:
    print('a is smaller than b')
else:
    print('a equals to b')

# ### NumPy

# [NumPy](http://www.numpy.org/) is a Python package for numerical computing and linear algebra. It is also a required installation for ProDy.

from numpy import *

# The most fundamental data type in numpy is `ndarray` (n-dimensional array). It can be used to represent vectors or matrices. There are various ways to initialize an vector in numpy:

v = arange(10, dtype=float)
v

v = zeros(5)
v

v.shape

# Or a matrix:

V = zeros([5, 5])
V

V.shape

V = ones([5, 5])
V

I = eye(5)
I

# ### Matplotlib

# [Matplotlib](https://matplotlib.org/) is a package used by ProDy to visualize the data. It is compatible with Python objects as well as numpy objects and provides [MATLAB](https://www.mathworks.com/products/matlab.html)-like syntax for plotting data points, matrices, 3D coordinates, etc.

from matplotlib.pyplot import *

x = arange(0, 2*pi, 0.1)
y = sin(x)
y

plot(x, y)
xlabel('x')
ylabel('y')
title('sin')

y2 = cos(x)
plot(x, y, x, y2)
xlabel('x')
ylabel('y')
legend(['sin', 'cos'])

plot(y, y2)

# <img src="http://prody.csb.pitt.edu/_static/ubiquitin.png" align="left" width="55" height="55"> <h1>Prody Basics</h1>

# This tutorial aims to teach basic data structures and functions in Prody. First we need to import required packages:

from prody import *
from numpy import *
from matplotlib.pyplot import *
# %matplotlib inline
confProDy(auto_show=False)
confProDy(auto_secondary=True)

# These import commands will load numpy, matplotlib, and ProDy into the memory. `confProDy` is used to modify the default behaviors of ProDy. Here we turned off `auto_show` so that the plots can be made in the same figure, and we turn on `auto_secondary` to parse the secondary structure information whenever we load a PDB into ProDy. See [here](http://prody.csb.pitt.edu/manual/reference/prody.html?highlight=confprody#prody.confProDy) for a complete list of behaviors that can be changed by this function. This function only needs to be called once, and the setting will be remembered by ProDy.

# ### Load PDB files and visualization

# ProDy comes with many functions that can be used to fetch data from [Protein Data Bank](https://www.rcsb.org/). Let's first parse a structure:

p38 = parsePDB('1p38')
p38

# `parsePDB` will download the PDB file and load it into the memory. Let's inspect the variable `p38`:

p38

# To visualize the structure:

showProtein(p38);
legend();

# If you would like to display the 3D structure using other packages or your own code, you can get the 3D coordinates via:

p38.getCoords()

p38.getCoords().shape

showContactMap(p38.ca);

# `AtomGroup` is essentially a collection of protein atoms. Each atom can be indexed/queried/found by the following way:

p38[10]

# This will give you the **11-th** atom from `p38`, noting that Python index starts from **0**. We can also examine the spatial location of this atom by querying the coordinate:

p38[10].getCoords()

showProtein(p38);
ax3d = gca()
x, y, z = p38[10].getCoords()
ax3d.plot([x], [y], [z], 'r*', markersize=30);

# We could select a chain, e.g. chain A, of the protein by indexing using its identifier, like the following:

p38['A']

p38['A'].getSequence()

# In many cases, it is more convenient to examine the structure with **[residue numbers](https://pdb101.rcsb.org/learn/guide-to-understanding-pdb-data/primary-sequences-and-the-pdb-format)**, and `AtomGroup` supports indexing with a chain ID and a residue number:

p38[10].getResnum()

p38['A', 5]

# This will give you the residue with the residue number of 10, which is an arginine in `p38`. Please note the difference between this line and the previous one. 

p38['A', 5].getNames()

# Note that some ProDy objects may not support indexing using a chain identifier or a residue number. In such cases, we can first obtain a hierarchical view of the object:

hv = p38.getHierView()

# And then use `HierView` to index using a chain identifier and residue number as it will always be supported:

hv['A', 5]

# ### Retrive data from AtomGroup

# Many properties of the protein can be acquired by functions named like "getxxx". For instance, we can obtain the B-factors by:

betas = p38.getBetas()
betas.shape

# In this way, we can obtain the B-factor for every single atom. However, in some cases, we only need to know the B-factors of alpha-carbons, 

p38.ca

betas = p38.ca.getBetas()
betas.shape

plot(betas);
ylabel('B-factor');
xlabel('Residue index');

# If we would like to use residue numbers in the PDB, instead of the indices, as the x-axis of the plot, it would be much more convenient to use the ProDy plotting function, `showAtomicLines`.

showAtomicLines(betas, atoms=p38.ca);
ylabel('B-factor');
xlabel('Residue number');

# We can also obtain the secondary structure information as an array:

p38.ca.getSecstrs()

# To make it easier to read, we can convert the array into a string using the Python's built-in function, `join` :

''.join(p38.ca.getSecstrs())

# `C` is for coil, `H` for alpha helix, `I` for pi helix, `G` for 3-10 helix, and `E` for beta strand (sheet). To get a complete list of "get" function, you can type `p38.get<TAB>`.

p38

# Speaking of which, in [measure](http://prody.csb.pitt.edu/manual/reference/measure/index.html?highlight=measure#module-prody.measure) module, you can find various functions for calculations for structural properties. For example, you can calculate the phi angle of 11th residue:

calcPhi(p38['A', 10])

round(calcPhi(p38['A', 10]), 3)

# A [dihedral angle](https://en.wikipedia.org/wiki/Dihedral_angle) is the angle between two intersecting planes. In chemistry it is the angle between planes through two sets of three atoms, having two atoms in common. In proteins, there are two most interested dihedral angles, namely Phi and Psi, and they are illustrated as follows.
# <img src="https://pitt.box.com/shared/static/taq12ir6r4s7rskfyycm9t2zjva0eawi.png" width="400" height="400">
# <img src="https://pitt.box.com/shared/static/q23s5crz0svu8yu2v0br7gab2fa7y3o0.png" width="400" height="400">

# Note that the residue at N-terminus or C-terminus does not have a Phi or Psi angle, respectively. If we calculate the Phi and Psi angle for every non-terminal residue, we can obtain a [Ramachandran plot](https://en.wikipedia.org/wiki/Ramachandran_plot) for a protein. An example of Ramachandran plot for human [PCNA](https://en.wikipedia.org/wiki/Proliferating_cell_nuclear_antigen) is shown as follows:
# <img src="https://upload.wikimedia.org/wikipedia/commons/4/43/1axc_PCNA_ProCheck_Rama.jpg" align="bottom" width="300" height="300">
#
# Three favored regions (in red)--**upper left: beta sheet; center left: alpha helix; center right: left-handed helix**. Each blue data point corresponds to the two dihedrals of a residue. We will reproduce this plot for ubiquitin (we will only reproduce the points).

chain = p38['A']
Phi = []; Psi = []; c = []
for res in chain.iterResidues():
    try:
        phi = calcPhi(res)
        psi = calcPsi(res)
    except:
        continue
    else:
        Phi.append(phi)
        Psi.append(psi)
        if res.getResname() == 'GLY':
            c.append('black')
        else:
            secstr = res.getSecstrs()[0]
            if secstr == 'H':
                c.append('red')
            elif secstr == 'G':
                c.append('darkred')
            elif secstr == 'E':
                c.append('blue')
            else:
                c.append('grey')

# In the above code, we use an exception handler to exclude the terminal residues from the calculation.

scatter(Phi, Psi, c=c, s=10);
xlabel('Phi (degree)');
ylabel('Psi (degree)');

# ### Selection

# In theory you could retrieve any set of atoms by indexing the `AtomGroup`, but it would be cumbersome to do so. To make it more convienient, ProDy provides VMD-like syntax for selecting atoms. Here lists a few common selection strings, and for a more complete tutorial on selection, please see [here](http://prody.csb.pitt.edu/tutorials/prody_tutorial/selection.html).

ca = p38.select('calpha')
ca

p38.ca

bb = p38.select('backbone')
p38.bb

# We could also perform some simple selections right when the structure is being parsed. For example, we can specify that we would like to obtain only alpha-carbons of chain A of the p38 as follows:

chainA_ca = parsePDB('1p38', chain='A', subset='ca')

# We could find the chain A using selection (as an alternative to the indexing method shown above):

chA = p38.select('calpha and chain A')
chA

# Selection also works for finding a single residue or multiple residues:

res = p38.ca.select('chain A and resnum 10')
res.getResnums()

res = p38.ca.select('chain A and resnum 10 11 12')
res.getResnums()

head = p38.ca.select('resnum < 50')
head.numAtoms()

# We can also select a range of residues by:

fragment = p38.ca.select('resnum 50 to 100')

# If we have data associated to the full length of the protein, we can slice the data using the `sliceAtomicData`:

subbetas = sliceAtomicData(betas, atoms=p38.ca, select=fragment)

# We can visualize the data of this range using `showAtomicLines`:

showAtomicLines(subbetas, atoms=fragment);
xlabel('Residue number');
ylabel('B factor');

# Or highlight the subset in the plot of the whole protein:

showAtomicLines(betas, atoms=p38.ca, overlay=True);
showAtomicLines(subbetas, atoms=fragment, overlay=True);
xlabel('Residue number');
ylabel('B factor');

# Selection also allows us to extract particular amino acid types:

args = p38.ca.select('resname ARG')
args

# Again, combined with `sliceAtomicData` and `showAtomicLines`, we can highlight these residues in the plot of the whole protein:

argbetas = sliceAtomicData(betas, atoms=p38.ca, select=args)
showAtomicLines(betas, atoms=p38.ca, overlay=True);
showAtomicLines(argbetas, atoms=args, linespec='r*', overlay=True);
xlabel('Residue number');
ylabel('B factor');

# ### Compare and align structures

# You can also compare different structures using some of the methods in proteins module. Let’s parse another p38 MAP kinase structure. 

bound = parsePDB('1zz2')

# You can find similar chains in structure 1p38 and 1zz2 using `matchChains` function

results = matchChains(p38, bound)
results[0]

# In Python, a tuple (or any indexable objects) can be unpacked by:

apo_chA, bnd_chA, seqid, overlap = results[0]
apo_chA

bnd_chA

# The first two terms are the mapping of the proteins to each other. Then the third term is the sequence identity:

seqid

# and the forth term is the sequence coverage:

overlap

# If we calculate RMSD right now, we will obtain the value for the unsuperposed proteins:

calcRMSD(bnd_chA, apo_chA)

# After superposition, the RMSD will be much improved,

bnd_chA, transformation = superpose(bnd_chA, apo_chA)
calcRMSD(bnd_chA, apo_chA)

confProDy(auto_show=False)
showProtein(bnd_chA);
showProtein(apo_chA);
legend();

# To visualize the superposition of the full proteins, we need to apply transform matrix to the entire structure:

showProtein(p38);
showProtein(bound);
legend();

# # Advanced Visualization

# Using `matplotlib`, we only obtained a very simple linear representation of proteins. ProDy also support a more sophisticated way of visualizing proteins in 3D via [py3Dmol](http://3dmol.csb.pitt.edu/):

import py3Dmol
showProtein(p38)

# The limitation is that `py3Dmol` only works in an iPython notebook. You can always write out the protein to a PDB file and visualize it in an external program:

writePDB('bound_aligned.pdb', bnd_chA)

StarcoderdataJupyterScriptsDedupFiltered