<jupyter_start><jupyter_text>Now use MD Analysis to calculate no. of frames a center PEG residues and terminal PEG residue is with 4 Angstroms of BSA (0.9 nm restrained system) Import MDAnalysis<jupyter_code>from prot_polymer_analysis import get_protresd_list, aa_frmcount, grptwocnt_aa, gtwo_trjcnt 
from prot_polymer_analysis import frac_cont, bavg_frac_cnt, prot_poly_cntmovie, AA_list_org
# Import MDAnalysis
import MDAnalysis as mda
import MDAnalysis.analysis.distances as maa_dist<jupyter_output><empty_output><jupyter_text>First table will be total fractional contacts and oligomer occupancy values for each Rg value   Distance-based analysis Find residues that have at least one atom within a cutoff $d = 4.0$ Angstrom near water molecules in BSA/water simulationCalculate the number of surface bsa residues from a 1 ns BSA/water simulation<jupyter_code>#Units of Angstroms 
dmax = 4.0 
def middle_of_band(band_start, band_stop, plot_min=0, plot_max=60):
    half_way = (band_stop - band_start) / 2
    mid_band = band_start + half_way
    plot_fraction = (mid_band - plot_min) / (plot_max - plot_min)

    return plot_fraction<jupyter_output><empty_output><jupyter_text>0.9 nm PEG unrestrained Rg 100 ns trajectory C2 Load the rg = 1.2 nm (3 PEG N = 20 oligomer/BSA system)<jupyter_code># Set up the MD Simulation, Make sure you do gmx trjconv -s topol.tpr -f confout.gro  -o new_conf.pdb 
# -dump 0 -n bsaplga_nk.ndx to generate 
# a new pdb file that contains unique chain identifiers 
uC2_n20PEG_NR = mda.Universe("bsapeg_n20waterT2/C2_pegNoRes/0.9nmpeg_prodC2/bpeg_0.9nmNR.pdb"
                             , "bsapeg_n20waterT2/C2_pegNoRes/0.9nmpeg_prodC2/C2NR_0.9nmbpeg.xtc")
uC2_n20PEG_NR<jupyter_output><empty_output><jupyter_text>Check that we are on the first frame<jupyter_code>uC2_n20PEG_NR.trajectory.frame
pn20C2_lenNR = len(uC2_n20PEG_NR.trajectory)
pn20C2_lenNR
# Select one polymer chain, heavy atoms only 
#all_pn20 = u_pn20.select_atoms("(resname sPLG PLG tPLG and segid B) and not type H")

#Select all the PLGA residues, heavy atoms only
allC2_n20PEG_NR = uC2_n20PEG_NR.select_atoms("resname sPEG PEG tPEG and not type H")
#list(all_n20PEG_NR)
# Select BSA residues, heavy atoms only 
protC2_09nmNR = uC2_n20PEG_NR.select_atoms("protein and not type H")
protC2_09nmNR<jupyter_output><empty_output><jupyter_text>Calculate AA frame counts for PLGA residues, 1.2 nm RG restraint, 100ns trajectory<jupyter_code>#dmax = 4.0, protein group(4653 atoms), plga atom group (543 atoms), took 381.6 s (6 min 36s on 4 cores)
start = 0
end = pn20C2_lenNR - 1
s_time = timeit.default_timer()
h2diC2_09nmNR = aa_frmcount(protC2_09nmNR, allC2_n20PEG_NR, dmax, uC2_n20PEG_NR, start, end)
timeit.default_timer() - s_time
len(h2diC2_09nmNR.keys())
pr_res_PEG = list(protC2_09nmNR.residues)
ss_res_PEG = [str(row) for row in pr_res_PEG]
rkgC2_NR = {key:h2diC2_09nmNR[key][1] for key, value in h2diC2_09nmNR.items()}
plgC2_09nmNR_occ = pd.DataFrame(data=ss_res_PEG, columns=["BSA_des_res"])
plgC2_09nmNR_occ['mda_occ_0.9nmNRC2'] = plgC2_09nmNR_occ['BSA_des_res'].map(rkgC2_NR)
plgC2_09nmNR_occ['mda_occ_0.9nmNRC2'] = plgC2_09nmNR_occ['mda_occ_0.9nmNRC2'].replace('nan', np.nan).fillna(0)
plgC2_09nmNR_occ['mda_occ_0.9nmNRC2'] = plgC2_09nmNR_occ['mda_occ_0.9nmNRC2'].round(2)
plgC2_09nmNR_occ
plgC2_09nmNR_occ['mda_occ_0.9nmNRC2'][plgC2_09nmNR_occ['mda_occ_0.9nmNRC2'] != 0]
bsaC2_NR = np.array(list(protC2_09nmNR.resids)) # shape is 4652
mC2_occ_09NR = np.array(list(plgC2_09nmNR_occ['mda_occ_0.9nmNRC2'])) # shape is 583
mC2_occ = np.zeros(shape=(4653))
at_ind = np.where(bsaC2_NR[:-1] != bsaC2_NR[1:])[0]
at_in_nw = np.sort(np.append([0,4653],at_ind))
nw_v = 0
for i in range(583):
    b = at_in_nw[i+1] +1
    mC2_occ[nw_v:b] = mC2_occ_09NR[i]
    nw_v = at_in_nw[i+1] + 1 
np.nonzero(mC2_occ)<jupyter_output><empty_output><jupyter_text>Visualize Occupanct on protein<jupyter_code>protC2_09nmNR.occupancies = mC2_occ
protC2_09nmNR.occupancies
with mda.Writer("prot09nm_NoResC2.pdb") as pdb:
    pdb.write(protC2_09nmNR)
# Frame count and occupancy for each residue
#h2di<jupyter_output><empty_output><jupyter_text>Residue Importance: 0.9 nm unrestrained C2<jupyter_code># Need to fix function, the residue number are not counting the other 2 PLGA oligomers cuz of same resid number
trjmapC2_09nmPEGNR = prot_poly_cntmovie(protC2_09nmNR, allC2_n20PEG_NR, dmax, uC2_n20PEG_NR, 0, 10000)
#trj_ppmap_12nm_chC = prot_poly_cntmovie(prot, all_pn20_C, dmax, u_pn20, 0, 10000)
np.save('0.9nmPEG_NoResC2.npy', trjmapC2_09nmPEGNR)    # .npy extension is added if not given
trjmapC2_09nmPEGNR = np.load('0.9nmPEG_NoRes.npy', allow_pickle=True)
trjmapC2_09nmPEGNR[0].shape
np.sum(trjmapC2_09nmPEGNR[1000][0])
kj = np.zeros(shape=(10000, 583))
kj[:,582].shape
ppC2_09nmNR_ct = np.zeros(shape=(10000, 583))
for i in range(10000):
    for j in range(583):
        ppC2_09nmNR_ct[i][j] = np.sum(trjmapC2_09nmPEGNR[i][j])
ppC2_09nmtotNR = np.zeros(shape=(583))
for i in range(583):
    ppC2_09nmtotNR[i] = np.sum(ppC2_09nmNR_ct[:,i])
#pp_12nmtot
np.sum(ppC2_09nmtotNR)
y_pos = np.arange(583) + 1
wid = np.zeros(shape=583)
wid += 1.5
#wid
fig = plt.figure(figsize=(12,12))
fig.canvas.layout.width = '800px'
fig.canvas.layout.height = '700px'
plt.bar(y_pos, ppC2_09nmtotNR/np.sum(ppC2_09nmtotNR), align='center', color='#1D77CF',width=wid, alpha=0.4, label='0.9 nm PEG')
plt.title("BSA in water with PEG Unrestrained, 100 ns", fontsize=18)
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)
plt.xlim([0,600])
plt.ylim([0,0.2])
plt.legend(fontsize=14)
plt.ylabel(r'Normalized Total No. of PEG contacts C2', fontsize=15)
plt.xlabel(r'BSA Residue ID', fontsize=15)<jupyter_output><empty_output><jupyter_text>Total number of residues that are within 4 angstroms of a PEG oligomer residue within a 100 ns trajectory block C2<jupyter_code>a_a = ["GLY","ALA","VAL","LEU","ILE","MET","PHE","TRP","PRO","SER","THR","CYS","TYR","ASN","GLN","ASP"
               ,"GLU","LYS","ARG","HIS"]
# This code chunk gets the BSA residues and their corresponding number in a pandas dataframe 
red_bsa = []
bh = np.arange(0,584)
for i in range(583):
    b_str = str(list(protC2_09nmNR.residues[i:i+1]))
    if str(bh[i+1]) in b_str: 
        red_bsa.append(str(b_str[10:13])+" "+str(bh[i+1]))
prC2_res = list(protC2_09nmNR.residues)
ssC2_res = [str(row) for row in prC2_res]
rkg = {key:h2diC2_09nmNR[key][0] for key, value in h2diC2_09nmNR.items()}
plgC2_09nmaaNR = pd.DataFrame(data=ssC2_res, columns=["BSA_des_res"])
plgC2_09nmaaNR['mda_plga_frm_0.9nmNRC2'] = plgC2_09nmaaNR['BSA_des_res'].map(rkg)
plgC2_09nmaaNR['BSA_des_res'] = red_bsa
plgC2_09nmaaNR['mda_plga_frm_0.9nmNRC2'] = plgC2_09nmaaNR['mda_plga_frm_0.9nmNRC2'].replace('nan', np.nan).fillna(0)
plgC2_09nmaaNR.tail()
# Read in data from the oputput of wrapper.sh, where the frame count is given for each BSA residue that was within 
# 4 angstroms of PLGA trimer 
wat_data = pd.read_csv('occ_BSA1ns.txt', sep=" ", header=None, usecols=None ,index_col=None)
wat_data.columns = ["BSA_res_no","No. of frames (VMD)"]
wat_data = wat_data.drop("BSA_res_no", axis=1)

#pr_res = list(prot_09nm.residues)
#ss_res = [str(row) for row in pr_res]

wat_data['BSA_des_res'] = ssC2_res
wat_data = wat_data[['BSA_des_res',"No. of frames (VMD)"]]
#wat_data.head()

# load MDAnalysis values from MDA_BSA1ns.txt file(129003 atoms SOL group was used to calc. frame counts for txt.
# file)
h2ob_dict = json.load(open("MDA_BSA1ns.txt"))
wat_data['Mda_frames'] = wat_data['BSA_des_res'].map(h2ob_dict)


# From MD Analysis
#Get the count of bsa residues that have 1001 or 1002 frames ( I ran a 1 ns NPT simulation of 1 BSA in water )

#aa_count = pd.DataFrame(data=a_a)
c_list = []

for i in range(len(a_a)):
    count = 0
    for index, row in wat_data.iterrows():
        if a_a[i] in row["BSA_des_res"]:
            if row['Mda_frames'] == 1001: 
                count += 1
                #c_list.append(str(str(a_a[i])+"  "+str(row['No. of frames']))) 
            elif row['Mda_frames'] == 1000:
                count += 1
                #c_list.append(str(str(a_a[i])+"  "+str(row['No. of frames'])))
    c_list.append(str(str(a_a[i])+"  "+str(count)))
    
#c_list

# From VMD
#Get the count of bsa residues that have 1001 or 1002 frames ( I ran a 1 ns NPT simulation of 1 BSA in water )

#aa_count = pd.DataFrame(data=a_a)
vmd_list = []

for i in range(len(a_a)):
    count = 0
    for index, row in wat_data.iterrows():
        if a_a[i] in row["BSA_des_res"]:
            if row["No. of frames (VMD)"] == 1001: 
                count += 1
                #c_list.append(str(str(a_a[i])+"  "+str(row['No. of frames']))) 
            elif row["No. of frames (VMD)"] == 1002:
                count += 1
                #c_list.append(str(str(a_a[i])+"  "+str(row['No. of frames'])))
    vmd_list.append(str(str(a_a[i])+"  "+str(count)))

# Main difference is that Alanine 583 is counted for all 1001 frames. It seems VMD is unable to calc dist for that res
#vmd_list

#hydrophobic_res = ['ALA', 'ILE', 'LEU', 'VAL', 'GLY', 'PRO','PHE', 'TRP','MET']
#polar_res = ['ASN', 'CYS', 'GLN', 'SER', 'THR','TYR']
#neg_res = ['ASP', 'GLU']
#pos_res = ['ARG', 'HIS', 'LYS']
# aromatic_res = ['PHE', 'TRP', 'TYR', 'HIS']
#all_res = [pos_res, neg_res, polar_res, hydrophobic_res]

# Put the AA count in a pandas dataframe 
dg , ji = AA_list_org(c_list)
aa_count = pd.DataFrame(data=dg, index=None, columns=['Amino_acids'])
new_lf = pd.Series(data=ji, index=None)
vmg, vmdj = AA_list_org(vmd_list)
n2lf = pd.Series(data=vmdj, index=None)
aa_count['No_of_surf_res (MDAnalysis)'] = new_lf
aa_count['No_of_surf_res (VMD)'] = n2lf
aplC2_list = []

# Some residues don't have any contact with the 3 N = 20 PLGA oligomers within 100 ns,
# Put residues that do have contact with BSA in a separate list
for index, r_pl in plgC2_09nmaaNR.iterrows():
    if r_pl['mda_plga_frm_0.9nmNRC2'] != 0:
        aplC2_list.append(r_pl['BSA_des_res'])
        
# This chunk of code gets an AA count from the above list, in order 
# to get a total number of residues that contact BSA
cpl_lC2 = []

for index, r_a in aa_count.iterrows():
    count = 0
    for i in range(len(aplC2_list)):
        if r_a['Amino_acids'] in aplC2_list[i]:
            count += 1
    cpl_lC2.append(count)      
        
aa_count['pegC2_0.9nmNR_100ns'] = cpl_lC2
aa_count
# This gives the total number of residues that are within 4 angstroms of a PLGA oligomer residue
# within a 100 ns trajectory block
aa_count['pegC2_0.9nmNR_100ns'].sum()
# This gives the total number of residues that are within 4 angstroms of a water molecule
# within a 1 ns trajectory block
aa_count['No_of_surf_res (MDAnalysis)'].sum()
# This gives the total fraction of contacts within the 1.2 nm Rg 100 ns trajectory
aa_count['pegC2_0.9nmNR_100ns'].sum()/aa_count['No_of_surf_res (MDAnalysis)'].sum()<jupyter_output><empty_output><jupyter_text>Calculate mean occupancy and the standard deviation for 1.2 nm trajectoryNumpy mean and std function was used to calculate mean occupancy and std dev using occ values from aa_frmcount output<jupyter_code># Mean occupancy and std deviation 
ll_moC2 = [value[1] for key, value in h2diC2_09nmNR.items()]
print("Mean Occpancy (0.9 nm Rg): "+str(np.mean(ll_moC2)), "Occ. std. dev.: "+str(np.std(ll_moC2)))<jupyter_output>Mean Occpancy (0.9 nm Rg): 0.11570754716981133 Occ. std. dev.: 0.18748667736221486<jupyter_text>Calc. fractional contacts for each AA group type<jupyter_code>cd_12nm = frac_cont(h2diC2_09nmNR)
cd_12nm
cd = frac_cont(h2diC2_09nmNR)
kklh = []
for key, value in cd.items():
    kklh.append(value[1])
# Must substract aromatic residues, since they are already counted
sum(kklh) - cd['Aromatic'][1]
no_surf = aa_count['No_of_surf_res (MDAnalysis)'].sum()
no_surf
fcntC2_rg09nmNR, prgrpC2_09nmNR, aaC2_matx09nmNR = bavg_frac_cnt(5, protC2_09nmNR, allC2_n20PEG_NR, dmax, 
                                                                 uC2_n20PEG_NR, no_surf, 0, 10000)
fcntC2_rg09nmNR
fc09nmC2_meanNR = np.array([np.mean(fcntC2_rg09nmNR['Negative']), np.mean(fcntC2_rg09nmNR['Positive'])
                        ,np.mean(fcntC2_rg09nmNR['Polar']),np.mean(fcntC2_rg09nmNR['Hydrophobic'])
                        , np.mean(fcntC2_rg09nmNR['Aromatic'])])
fc09nmC2_meanNR
fc09nmC2_stdNR = np.array([np.std(fcntC2_rg09nmNR['Negative']), np.std(fcntC2_rg09nmNR['Positive'])
                       ,np.std(fcntC2_rg09nmNR['Polar']),np.std(fcntC2_rg09nmNR['Hydrophobic'])
                       , np.std(fcntC2_rg09nmNR['Aromatic'])])
fc09nmC2_stdNR
x_pos = np.arange(5)
aa_types = ["Negative", "Positive", "Polar", "Hydrophobic", "Aromatic"]
fig = plt.figure(figsize=(7,7))
fig.canvas.layout.width = '500px'
fig.canvas.layout.height = '400px'
plt.bar(x_pos, fc09nmC2_meanNR, yerr=fc09nmC2_stdNR, ecolor='black',capsize=5, color='c')
plt.title(r'Fractional Contacts 0.9 nm Rg restrained Config 2', fontsize=15)
plt.xticks(x_pos, labels=aa_types, fontsize=12)
plt.ylabel(r'Fractional Contacts', fontsize=15)<jupyter_output><empty_output><jupyter_text>Total fraction of contacts: averages and std dev calc from 5 20ns blocks<jupyter_code># Average of total fraction of contacts
np.mean(fcntC2_rg09nmNR['total_frac'])
# Std Deviation of total fraction of contacts
np.std(fcntC2_rg09nmNR['total_frac'])<jupyter_output><empty_output><jupyter_text>Avg no. PEG residues per BSA AA residue group C2<jupyter_code>prgrpC2_09nmNR
mean_09nmC2 = np.zeros(shape=5)
std_09nmC2 = np.zeros(shape=5)
count = 0
for key, value in prgrpC2_09nmNR.items():
    mpl_09nmC2 = []
    var_09nmC2 = []
    for i in prgrpC2_09nmNR[str(key)].flat:
        mpl_09nmC2.append(i[0])
        var_09nmC2.append((i[1])**2)
    
    # calc frac cont averages
    mean_09nmC2[count] = np.mean(mpl_09nmC2)
    
    # calc frac cont std dev: https://stats.stackexchange.com/questions/25848/how-to-sum-a-standard-deviation 
    std_09nmC2[count] = np.std(mpl_09nmC2)
    # std_12nm[count] = np.sqrt(np.sum(var_12nm)/5)
    
    count += 1

mean_09nmC2
std_09nmC2
x_pos = np.arange(5)
aa_types = ["Negative", "Positive", "Polar", "Hydrophobic", "Aromatic"]
fig = plt.figure(figsize=(7,7))
fig.canvas.layout.width = '500px'
fig.canvas.layout.height = '400px'
plt.bar(x_pos, mean_09nmC2, yerr=std_09nmC2, ecolor='black',capsize=5)
plt.title(r'No. of PEG residues 0.9 nm Rg restrained', fontsize=15)
plt.xticks(x_pos, labels=aa_types, fontsize=12)
plt.ylabel(r'No. of PEG residues', fontsize=15)<jupyter_output><empty_output><jupyter_text>Protein/polymer contact map movie<jupyter_code>fig = plt.figure(figsize=(10,10))

# Set the axis and the plot titles pp

plt.title("BSA/PEG contact map 0.9 nm Unrestrained", fontsize=22, loc='left')
plt.xlabel("PLGA Residue No.", fontsize=22)
plt.ylabel("BSA Residue No.", fontsize=20)

 # Set the axis range 
plt.ylim(583, 0)
plt.xlim(0, 60)

# Plot bands for each chain 
BANDS = (
    (0, 20, "purple", "B"),
    (20, 40, "blue", "C"),
    (40, 60, "green", "D"),
)
    
text_y = 0.98 # Close to the top
for start, stop, color, band in BANDS:
    plt.axvspan(start, stop,color=color, alpha=0.15)
    text_x = middle_of_band(start,stop)
    plt.text(
        text_x,
        text_y,
        "PLGA chain " + band,
        color=color,
        fontsize=18,
        transform=fig.gca().transAxes,
        horizontalalignment='center',
        verticalalignment='center',
        style='italic',
    )
    
plt.text(0.93, 1, "Time [ns]:", fontsize=20, transform=fig.gca().transAxes, horizontalalignment='right', verticalalignment='bottom')

# Set tick label size
fig.gca().tick_params(axis='both', which='major', labelsize=20)

ims = []
for i in range(10000):
    data = trj_ppmap_12nm[i]
    im = plt.imshow(data, aspect='auto', cmap='Greys')
    t_sim = plt.text(1, 1, str(i/100), fontsize=20, transform=fig.gca().transAxes, horizontalalignment='right', verticalalignment='bottom')
    ims.append([im, t_sim])
    
ani = animation.ArtistAnimation(fig, ims, blit=True, repeat=False)
ani.save('1.2nm_res.mp4',writer='ffmpeg', fps=50, bitrate=100000)
#plt.tight_layout()
#plt.show()<jupyter_output><empty_output><jupyter_text>1.1 nm PEG unrestrained Rg 100 ns trajectory C2 Load the rg = 1.5 nm (3 PLGA N = 20 oligomer/BSA system)<jupyter_code># Set up the MD Simulation
u11nmNR_n20PEGC2 = mda.Universe("bsapeg_n20waterT2/C2_pegNoRes/1.1nmpeg_prodC2/bpeg_1.1nmNR.pdb"
                                , "bsapeg_n20waterT2/C2_pegNoRes/1.1nmpeg_prodC2/C2NR_1.1nmbpeg.xtc")
u11nmNR_n20PEGC2
pn20C2_11nmNR = len(u11nmNR_n20PEGC2.trajectory)
pn20C2_11nmNR
#Select all the PLGA residues, heavy atoms only 
n20PEG_11nmNRC2 = u11nmNR_n20PEGC2.select_atoms("resname sPEG PEG tPEG and not type H")
n20PEG_11nmNRC2
# Select BSA residues, heavy atoms only 
protNR_11nmPEGC2 = u11nmNR_n20PEGC2.select_atoms("protein and not type H")
protNR_11nmPEGC2<jupyter_output><empty_output><jupyter_text>Contact Analysis<jupyter_code>#dmax = 4.0, protein group(4653 atoms), plga atom group (543 atoms), took 381.6 s (6 min 36s on 4 cores)
start = 0
end = pn20C2_11nmNR - 1
s_time = timeit.default_timer()
h2diC2_11nmNR = aa_frmcount(protNR_11nmPEGC2, n20PEG_11nmNRC2, dmax, u11nmNR_n20PEGC2, start, end)
timeit.default_timer() - s_time
#h2di_11nm
len(h2diC2_11nmNR.keys())
pr_res_PEG = list(protNR_11nmPEGC2.residues)
ss_res_PEG = [str(row) for row in pr_res_PEG]
rkgC2_11NR = {key:h2diC2_11nmNR[key][1] for key, value in h2diC2_11nmNR.items()}
plgC2_09nmNR_occ['mda_occ_1.1nmNRC2'] = plgC2_09nmNR_occ['BSA_des_res'].map(rkgC2_11NR)
plgC2_09nmNR_occ['mda_occ_1.1nmNRC2'] = plgC2_09nmNR_occ['mda_occ_1.1nmNRC2'].replace('nan', np.nan).fillna(0)
plgC2_09nmNR_occ['mda_occ_1.1nmNRC2'] = plgC2_09nmNR_occ['mda_occ_1.1nmNRC2'].round(2)
plgC2_09nmNR_occ
pr_res11nmC2 = list(protNR_11nmPEGC2.residues)
ss_res11nmC2 = [str(row) for row in pr_res11nmC2]
rkgC2_11nmNR = {key:h2diC2_11nmNR[key][0] for key, value in h2diC2_11nmNR.items()}
plgC2_1_1nmaaNR = pd.DataFrame(data=ss_res11nmC2, columns=["BSA_des_res"])
plgC2_1_1nmaaNR['mda_plga_frm_1.1nmNRC2'] = plgC2_1_1nmaaNR['BSA_des_res'].map(rkgC2_11nmNR)
plgC2_1_1nmaaNR['BSA_des_res'] = red_bsa
plgC2_1_1nmaaNR['mda_plga_frm_1.1nmNRC2'] = plgC2_1_1nmaaNR['mda_plga_frm_1.1nmNRC2'].replace('nan', np.nan).fillna(0)
plgC2_1_1nmaaNR.head()
bsa11C2_NR = np.array(list(protNR_11nmPEGC2.resids)) # shape is 4652
mC2_occ_11nmNR = np.array(list(plgC2_09nmNR_occ['mda_occ_1.1nmNRC2'])) # shape is 583
mC2_occ_11New = np.zeros(shape=(4653))
at_indC2 = np.where(bsa11C2_NR[:-1] != bsa11C2_NR[1:])[0]
atC2_in_nw = np.sort(np.append([0,4653],at_indC2))
nw_v = 0
for i in range(583):
    b = atC2_in_nw[i+1] +1
    mC2_occ_11New[nw_v:b] = mC2_occ_11nmNR[i]
    nw_v = atC2_in_nw[i+1] + 1 
mC2_occ_11New[0:33]
protNR_11nmPEGC2.occupancies[3089:3099]
np.nonzero(mC2_occ_11New)<jupyter_output><empty_output><jupyter_text>Visualize Occupanct on protein<jupyter_code>protNR_11nmPEGC2.occupancies = mC2_occ_11New
protNR_11nmPEGC2.occupancies[0:33]
with mda.Writer("prot11nm_NoResC2.pdb") as pdb:
    pdb.write(protNR_11nmPEGC2)<jupyter_output><empty_output><jupyter_text>Residue Importance: 1.1 nm unrestrained<jupyter_code># Need to fix function, the residue number are not counting the other 2 PLGA oligomers cuz of same resid number
trjmapC2_11nmPEGNR = prot_poly_cntmovie(protNR_11nmPEGC2, n20PEG_11nmNRC2, dmax, u11nmNR_n20PEGC2, 0, 10000)
#trj_ppmap_12nm_chC = prot_poly_cntmovie(prot, all_pn20_C, dmax, u_pn20, 0, 10000)
trjmapC2_11nmPEGNR = np.load('1.1nmNR_PEGres.npy', allow_pickle=True)
np.save('1.1nmNR_PEGresC2.npy', trjmapC2_11nmPEGNR)    # .npy extension is added if not given
trjmapC2_11nmPEGNR[0].shape
np.sum(trjmapC2_11nmPEGNR[1000][0])
kj = np.zeros(shape=(10000, 583))
kj[:,582].shape
ppC2_11nmNR_ct = np.zeros(shape=(10000, 583))
for i in range(10000):
    for j in range(583):
        ppC2_11nmNR_ct[i][j] = np.sum(trjmapC2_11nmPEGNR[i][j])
ppC2_11nmtotNR = np.zeros(shape=(583))
for i in range(583):
    ppC2_11nmtotNR[i] = np.sum(ppC2_11nmNR_ct[:,i])
#pp_12nmtot
np.sum(ppC2_11nmtotNR)
y_pos = np.arange(583) + 1
wid = np.zeros(shape=583)
wid += 1.5
#wid
fig = plt.figure(figsize=(12,12))
fig.canvas.layout.width = '800px'
fig.canvas.layout.height = '700px'
plt.bar(y_pos, ppC2_11nmtotNR/np.sum(ppC2_11nmtotNR), align='center',width=wid, color='#562A8B', alpha=0.3, label='1.1 nm PEG')
#plt.bar(y_pos, pp_09nmtotNR, align='center',width=wid, alpha=0.5, color='#1D77CF',label='0.9 nm PEG')
plt.title("BSA in water with PEG Unrestrained, 100 ns, Config 2", fontsize=18)
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)
plt.xlim([0,600])
plt.ylim([0,0.2])
plt.legend(fontsize=14)
plt.ylabel(r'Normalized Total No. of PEG contacts', fontsize=15)
plt.xlabel(r'BSA Residue ID', fontsize=15)<jupyter_output><empty_output><jupyter_text>Total number of residues that are within 4 angstroms of a PEG oligomer residue within a 100 ns trajectory block<jupyter_code>aplC2_11nmNR = []

# Some residues don't have any contact with the 3 N = 20 PLGA oligomers within 100 ns,
# Put residues that do have contact with BSA in a separate list
for index, r_pl in plgC2_1_1nmaaNR.iterrows():
    if r_pl['mda_plga_frm_1.1nmNRC2'] != 0:
        aplC2_11nmNR.append(r_pl['BSA_des_res'])
        
# This chunk of code gets an AA count from the above list, in order 
# to get a total number of residues that contact BSA
cplC2_11nmNR = []

for index, r_a in aa_count.iterrows():
    count = 0
    for i in range(len(aplC2_11nmNR)):
        if r_a['Amino_acids'] in aplC2_11nmNR[i]:
            count += 1
    cplC2_11nmNR.append(count)      
        
aa_count['pegC2_1.1nmNR_100ns'] = cplC2_11nmNR
#aa_count.drop('No_of_surf_res (VMD)', axis=1, inplace=True)
aa_count
# This gives the total number of residues that are within 4 angstroms of a PLGA oligomer residue
# within a 100 ns trajectory block
aa_count['pegC2_1.1nmNR_100ns'].sum()
# This gives the total number of residues that are within 4 angstroms of a water molecule
# within a 1 ns trajectory block
aa_count['No_of_surf_res (MDAnalysis)'].sum()
# This gives the total fraction of contacts within the 1.2 nm Rg 100 ns trajectory
aa_count['pegC2_1.1nmNR_100ns'].sum()/aa_count['No_of_surf_res (MDAnalysis)'].sum()
# Mean occupancy and std deviation 
ll_mo11 = [value[1] for key, value in h2diC2_11nmNR.items()]
print("Mean Occpancy (1.1 nm Rg): "+str(np.mean(ll_mo11)), "Occ. std. dev.: "+str(np.std(ll_mo11)))
cd_11nmNRC2 = frac_cont(h2diC2_11nmNR)
cd_11nmNRC2<jupyter_output><empty_output><jupyter_text>Calc. fractional contacts for each AA group type<jupyter_code>fcntC2_rg11nmNR, prgrpC2_11nmNR, aamatxC2_11nmNR = bavg_frac_cnt(5, protNR_11nmPEGC2, n20PEG_11nmNRC2,
                                                             dmax ,u11nmNR_n20PEGC2, no_surf, 0, 10000)
fcntC2_rg11nmNR
fc11nmC2_meanNR = np.array([np.mean(fcntC2_rg11nmNR['Negative']), np.mean(fcntC2_rg11nmNR['Positive'])
                        ,np.mean(fcntC2_rg11nmNR['Polar']),np.mean(fcntC2_rg11nmNR['Hydrophobic'])
                        , np.mean(fcntC2_rg11nmNR['Aromatic'])])
fc11nmC2_meanNR
fc11nmC2_stdNR = np.array([np.std(fcntC2_rg11nmNR['Negative']), np.std(fcntC2_rg11nmNR['Positive'])
                       ,np.std(fcntC2_rg11nmNR['Polar']),np.std(fcntC2_rg11nmNR['Hydrophobic'])
                       , np.std(fcntC2_rg11nmNR['Aromatic'])])
fc11nmC2_stdNR
x_pos = np.arange(5)
width = 0.35
aa_types = ["Negative", "Positive", "Polar", "Hydrophobic", "Aromatic"]
fig = plt.figure(figsize=(7,7))
fig.canvas.layout.width = '500px'
fig.canvas.layout.height = '400px'
plt.bar(x_pos, fc09nmC2_meanNR, width, yerr=fc09nmC2_stdNR, ecolor='black',capsize=5, color='royalblue')
plt.bar(x_pos+width, fc11nmC2_meanNR, width, yerr=fc11nmC2_stdNR, ecolor='black',capsize=5, color='c')
plt.title(r'Fractional Contacts Rg restrained, Config 2', fontsize=15)
plt.xticks(x_pos+width/2, labels=aa_types, fontsize=12)
plt.legend(['Rg = 0.9 nm', 'Rg = 1.1 nm'], frameon=False)
plt.ylabel(r'Fractional Contacts', fontsize=15)<jupyter_output><empty_output><jupyter_text>Total fraction of contacts: averages and std dev calc from 5 20 ns blocks<jupyter_code>np.mean(fcntC2_rg11nmNR['total_frac'])
np.std(fcntC2_rg11nmNR['total_frac'])<jupyter_output><empty_output><jupyter_text>Avg no. PEG residues per BSA AA residue group C2<jupyter_code>prgrpC2_11nmNR
mean_11nmC2 = np.zeros(shape=5)
std_11nmC2 = np.zeros(shape=5)
count = 0
for key, value in prgrpC2_11nmNR.items():
    mpl_11nm = []
    var_11nm = []
    for i in prgrpC2_11nmNR[str(key)].flat:
        mpl_11nm.append(i[0])
        var_11nm.append((i[1])**2)
    
    # calc frac cont averages
    mean_11nmC2[count] = np.mean(mpl_11nm)
    
    # calc frac cont std dev: https://stats.stackexchange.com/questions/25848/how-to-sum-a-standard-deviation 
    std_11nmC2[count] = np.std(mpl_11nm)
    #std_15nm[count] = np.sqrt(np.sum(var_15nm)/5)
    
    count += 1

mean_11nmC2
std_11nmC2
x_pos = np.arange(5)
width = 0.35
aa_types = ["Negative", "Positive", "Polar", "Hydrophobic", "Aromatic"]
fig = plt.figure(figsize=(7,7))
fig.canvas.layout.width = '500px'
fig.canvas.layout.height = '400px'
plt.bar(x_pos, mean_09nmC2, width, yerr=std_09nmC2, ecolor='black',capsize=5, color='royalblue')
plt.bar(x_pos+width, mean_11nmC2, width, yerr=std_11nmC2, ecolor='black',capsize=5, color='c')
plt.title(r'No. of PEG residues Rg restrained, Config 2', fontsize=15)
plt.xticks(x_pos+width/2, labels=aa_types, fontsize=12)
plt.legend(['Rg = 0.9 nm', 'Rg = 1.1 nm'], frameon=False)
plt.ylabel(r'No. of PLGA residues', fontsize=15)<jupyter_output><empty_output><jupyter_text>Protein/polymer contact map movie<jupyter_code>fig = plt.figure(figsize=(10,10))

# Set the axis and the plot titles pp

plt.title("BSA/PLGA contact map 1.5 nm res.", fontsize=22, loc='left')
plt.xlabel("PLGA Residue No.", fontsize=22)
plt.ylabel("BSA Residue No.", fontsize=20)

 # Set the axis range 
plt.ylim(583, 0)
plt.xlim(0, 60)

# Plot bands for each chain 
BANDS = (
    (0, 20, "purple", "B"),
    (20, 40, "blue", "C"),
    (40, 60, "green", "D"),
)
    
text_y = 0.98 # Close to the top
for start, stop, color, band in BANDS:
    plt.axvspan(start, stop,color=color, alpha=0.15)
    text_x = middle_of_band(start,stop)
    plt.text(
        text_x,
        text_y,
        "PLGA chain " + band,
        color=color,
        fontsize=18,
        transform=fig.gca().transAxes,
        horizontalalignment='center',
        verticalalignment='center',
        style='italic',
    )
    
plt.text(0.94, 1, "Time [ns]:", fontsize=20, transform=fig.gca().transAxes, horizontalalignment='right', verticalalignment='bottom')

# Set tick label size
fig.gca().tick_params(axis='both', which='major', labelsize=20)

ims = []
for i in range(10000):
    data = trj_ppmap_15nm[i]
    im = plt.imshow(data, aspect='auto', cmap='Greys')
    t_sim = plt.text(1.03, 1, str(i/100), fontsize=20, transform=fig.gca().transAxes, horizontalalignment='right', verticalalignment='bottom')
    ims.append([im, t_sim])
    
ani = animation.ArtistAnimation(fig, ims, blit=True, repeat=False)
ani.save('1.5nm_res.mp4', writer='ffmpeg', fps=50, bitrate=100000)
#plt.tight_layout()
#plt.show()<jupyter_output><empty_output><jupyter_text>1.4 nm unrestrained Rg PEG 100 ns trajectory C2  Contact Analyis<jupyter_code># Set up the MD Simulation
u14nmNR_n20PEGC2 = mda.Universe("bsapeg_n20waterT2/C2_pegNoRes/1.4nmpeg_prodC2/bpeg_1.4nmNR.pdb"
                              , "bsapeg_n20waterT2/C2_pegNoRes/1.4nmpeg_prodC2/C2NR_1.4nmbpeg.xtc")
u14nmNR_n20PEGC2
pn20C2_len14nmNR = len(u14nmNR_n20PEGC2.trajectory)
pn20C2_len14nmNR
#Select all the PLGA residues, heavy atoms only 
pn20C2_14nmNR = u14nmNR_n20PEGC2.select_atoms("resname sPEG PEG tPEG and not type H")
pn20C2_14nmNR
# Select BSA residues, heavy atoms only 
protC2_14nmNR = u14nmNR_n20PEGC2.select_atoms("protein and not type H")
protC2_14nmNR
#dmax = 4.0, protein group(4653 atoms), plga atom group (543 atoms), took 381.6 s (6 min 36s on 4 cores)
start = 0
end = pn20C2_len14nmNR - 1
s_time = timeit.default_timer()
h2diC2_14nmNR = aa_frmcount(protC2_14nmNR, pn20C2_14nmNR, dmax, u14nmNR_n20PEGC2, start, end)
timeit.default_timer() - s_time
#h2di_14nm
len(h2diC2_14nmNR.keys())
prC2_res14nmNR = list(protC2_14nmNR.residues)
ssC2_res14nmNR = [str(row) for row in prC2_res14nmNR]
rkgC2_14nmNR = {key:h2diC2_14nmNR[key][0] for key, value in h2diC2_14nmNR.items()}
plgC2_14nmaaNR = pd.DataFrame(data=ssC2_res14nmNR, columns=["BSA_des_res"])
plgC2_14nmaaNR['mda_plga_frm_1.4nmNRC2'] = plgC2_14nmaaNR['BSA_des_res'].map(rkgC2_14nmNR)
plgC2_14nmaaNR['BSA_des_res'] = red_bsa
plgC2_14nmaaNR['mda_plga_frm_1.4nmNRC2'] = plgC2_14nmaaNR['mda_plga_frm_1.4nmNRC2'].replace('nan', np.nan).fillna(0)
plgC2_14nmaaNR.head()
# Extract mean occupancy values
pr_res_14ur = list(protC2_14nmNR.residues)
ss_res_14ur = [str(row) for row in pr_res_14ur]
rkgC2_14NR = {key:h2diC2_14nmNR[key][1] for key, value in h2diC2_14nmNR.items()}
plgC2_09nmNR_occ['mda_occ_1.4nmNRC2'] = plgC2_09nmNR_occ['BSA_des_res'].map(rkgC2_14NR)
plgC2_09nmNR_occ['mda_occ_1.4nmNRC2'] = plgC2_09nmNR_occ['mda_occ_1.4nmNRC2'].replace('nan', np.nan).fillna(0)
plgC2_09nmNR_occ['mda_occ_1.4nmNRC2'] = plgC2_09nmNR_occ['mda_occ_1.4nmNRC2'].round(2)
plgC2_09nmNR_occ
bsa14C2_NR = np.array(list(protC2_14nmNR.resids)) # shape is 4652
mC2_occ_14nmNR = np.array(list(plgC2_09nmNR_occ['mda_occ_1.4nmNRC2'])) # shape is 583
mC2_occ_14New = np.zeros(shape=(4653))
at_ind = np.where(bsaC2_NR[:-1] != bsa14C2_NR[1:])[0]
at_in_nw = np.sort(np.append([0,4653],at_ind))
nw_v = 0
for i in range(583):
    b = at_in_nw[i+1] +1
    mC2_occ_14New[nw_v:b] = mC2_occ_14nmNR[i]
    nw_v = at_in_nw[i+1] + 1 
np.nonzero(mC2_occ_14New)<jupyter_output><empty_output><jupyter_text>Visualize Occupanct on protein<jupyter_code>protC2_14nmNR.occupancies = mC2_occ_14New
protC2_14nmNR.occupancies
with mda.Writer("prot14nm_NoResC2.pdb") as pdb:
    pdb.write(protC2_14nmNR)<jupyter_output><empty_output><jupyter_text>Residue Importance: 1.4 nm unrestrained<jupyter_code># Need to fix function, the residue number are not counting the other 2 PLGA oligomers cuz of same resid number
trjmapC2_14nmPEGNR = prot_poly_cntmovie(protC2_14nmNR, pn20C2_14nmNR, dmax, u14nmNR_n20PEGC2, 0, 10000)
#trj_ppmap_12nm_chC = prot_poly_cntmovie(prot, all_pn20_C, dmax, u_pn20, 0, 10000)
np.save('1.4nmPEG_NoResC2.npy', trjmapC2_14nmPEGNR)    # .npy extension is added if not given
trjmapC2_14nmPEGNR = np.load('1.4nm_NoRes.npy', allow_pickle=True)
trjmapC2_14nmPEGNR
trjmapC2_14nmPEGNR[0].shape
np.sum(trjmapC2_14nmPEGNR[1000][0])
kj = np.zeros(shape=(10000, 583))
kj[:,582].shape
ppC2_14nmNR_ct = np.zeros(shape=(10000, 583))
for i in range(10000):
    for j in range(583):
        ppC2_14nmNR_ct[i][j] = np.sum(trjmapC2_14nmPEGNR[i][j])
ppC2_14nmtotNR = np.zeros(shape=(583))
for i in range(583):
    ppC2_14nmtotNR[i] = np.sum(ppC2_14nmNR_ct[:,i])
#pp_12nmtot
np.nonzero(ppC2_14nmtotNR)
a_peg = np.sum(ppC2_14nmtotNR)
b_peg = np.sum(ppC2_11nmtotNR)
c_peg = np.sum(ppC2_09nmtotNR)
plt.close('all')
y_pos = np.arange(583) + 1
wid = np.zeros(shape=583)
wid += 1.5
#wid
fig = plt.figure(figsize=(12,12))
fig.canvas.layout.width = '800px'
fig.canvas.layout.height = '700px'
#plt.bar(y_pos, ppC2_09nmtotNR/c_peg, align='center',width=wid, alpha=0.5, color='#1D77CF',label='0.9 nm PEG')
#plt.bar(y_pos+0.25, ppC2_11nmtotNR/b_peg, align='center',width=wid, color='#562A8B', alpha=0.3, label='1.1 nm PEG')
plt.bar(y_pos+0.3, ppC2_14nmtotNR/a_peg, align='center',width=wid, color='#4E4C4D', alpha=0.3, label='1.4 nm PEG')
#plt.bar(y_pos+0.25, pp_11nmtot/b_peg, align='center',width=wid, color='#562A8B', alpha=0.3, label='1.1 nm PEG')
#plt.bar(y_pos, pp_09nmtot/c_peg, align='center',width=wid, alpha=0.5, color='#1D77CF',label='0.9 nm PEG')
plt.title("BSA in water with PEG Unrestrained, 100 ns, Config. 2", fontsize=16)
plt.xticks(fontsize=16)
plt.yticks(fontsize=16)
plt.xlim([0,600])
plt.ylim([0,0.2])
plt.legend(fontsize=16)
plt.ylabel(r'Normalized Total No. of PEG contacts', fontsize=16)
plt.xlabel(r'BSA Residue ID', fontsize=16)<jupyter_output><empty_output><jupyter_text>Total number of residues that are within 4 angstroms of a PEG oligomer residue within a 100 ns trajectory block<jupyter_code>apl_14nmNR = []

# Some residues don't have any contact with the 3 N = 20 PLGA oligomers within 100 ns,
# Put residues that do have contact with BSA in a separate list
for index, r_pl in plgC2_14nmaaNR.iterrows():
    if r_pl['mda_plga_frm_1.4nmNRC2'] != 0:
        apl_14nmNR.append(r_pl['BSA_des_res'])
        
# This chunk of code gets an AA count from the above list, in order 
# to get a total number of residues that contact BSA
cpl_14nmNR = []

for index, r_a in aa_count.iterrows():
    count = 0
    for i in range(len(apl_14nmNR)):
        if r_a['Amino_acids'] in apl_14nmNR[i]:
            count += 1
    cpl_14nmNR.append(count)      
        
aa_count['pegC2_1.4nmNR_100ns'] = cpl_14nmNR
#aa_count.drop('No_of_surf_res (VMD)', axis=1, inplace=True)
aa_count
# This gives the total number of residues that are within 4 angstroms of a PLGA oligomer residue
# within a 100 ns trajectory block
aa_count['pegC2_1.4nmNR_100ns'].sum()
# This gives the total number of residues that are within 4 angstroms of a water molecule
# within a 1 ns trajectory block
aa_count['No_of_surf_res (MDAnalysis)'].sum()
# This gives the total fraction of contacts within the 1.2 nm Rg 100 ns trajectory
aa_count['pegC2_1.4nmNR_100ns'].sum()/aa_count['No_of_surf_res (MDAnalysis)'].sum()
# Mean occupancy and std deviation 
ll_mo14nm = [value[1] for key, value in h2diC2_14nmNR.items()]
print("Mean Occpancy (1.4 nm Rg): "+str(np.mean(ll_mo14nm)), "Occ. std. dev.: "+str(np.std(ll_mo14nm)))
cd_14NRnm = frac_cont(h2diC2_14nmNR)
cd_14NRnm<jupyter_output><empty_output><jupyter_text>Calc. fractional contacts for each AA group type<jupyter_code>fcntC2_rg14nmNR, prgrpC2_14nmNR, aa_matxC2_14nmNR = bavg_frac_cnt(5, protC2_14nmNR, pn20C2_14nmNR, dmax
                                                                  , u14nmNR_n20PEGC2, no_surf, 0, 10000)
fcntC2_rg14nmNR
fc14nmC2_meanNR = np.array([np.mean(fcntC2_rg14nmNR['Negative']), np.mean(fcntC2_rg14nmNR['Positive'])
                        ,np.mean(fcntC2_rg14nmNR['Polar']),np.mean(fcntC2_rg14nmNR['Hydrophobic'])
                        , np.mean(fcntC2_rg14nmNR['Aromatic'])])
fc14nmC2_meanNR
fc14nmC2_stdNR = np.array([np.std(fcntC2_rg14nmNR['Negative']), np.std(fcntC2_rg14nmNR['Positive'])
                       ,np.std(fcntC2_rg14nmNR['Polar']),np.std(fcntC2_rg14nmNR['Hydrophobic'])
                       , np.std(fcntC2_rg14nmNR['Aromatic'])])
fc14nmC2_stdNR
plt.close("all")
x_pos = np.arange(5)
width = 0.28
aa_types = ["Negative", "Positive", "Polar", "Hydrophobic", "Aromatic"]
fig = plt.figure(figsize=(7,7))
fig.canvas.layout.width = '800px'
fig.canvas.layout.height = '700px'
plt.bar(x_pos, fc09nmC2_meanNR, width, yerr=fc09nmC2_stdNR, ecolor='black',capsize=5, color='royalblue')
plt.bar(x_pos+width, fc11nmC2_meanNR, width, yerr=fc11nmC2_stdNR, ecolor='black',capsize=5, color='c')
plt.bar(x_pos+(2*width), fc14nmC2_meanNR, width, yerr=fc14nmC2_stdNR, ecolor='black',capsize=5, color='lightslategray')
plt.title(r'Fractional Contacts Rg PEG/BSA Unrestrained, Config. 2', fontsize=16)
plt.xticks(x_pos+width, labels=aa_types, fontsize=16)
plt.yticks(fontsize=16)
plt.ylim(0,0.4)
plt.legend(['Rg = 0.9 nm', 'Rg = 1.1 nm', 'Rg = 1.4 nm'], frameon=False, fontsize=16)
plt.ylabel(r'Fractional Contacts', fontsize=16)<jupyter_output><empty_output><jupyter_text>Total fraction of contacts: averages and std dev calc from 5 20 ns blocks<jupyter_code>np.mean(fcntC2_rg14nmNR['total_frac'])
np.std(fcntC2_rg14nmNR['total_frac'])
prgrpC2_14nmNR<jupyter_output><empty_output><jupyter_text>Avg no. PLGA residues per BSA AA residue group<jupyter_code>prgrpC2_14nmNR
mean_14nmC2 = np.zeros(shape=5)
std_14nmC2 = np.zeros(shape=5)
count = 0
for key, value in prgrpC2_14nmNR.items():
    mpl_14nm = []
    var_14nm = []
    for i in prgrpC2_14nmNR[str(key)].flat:
        mpl_14nm.append(i[0])
        var_14nm.append((i[1])**2)
    
    # calc frac cont averages
    mean_14nmC2[count] = np.mean(mpl_14nm)
    
    # calc frac cont std dev: https://stats.stackexchange.com/questions/25848/how-to-sum-a-standard-deviation 
    std_14nmC2[count] = np.std(mpl_14nm)
    #std_2nm[count] = np.sqrt(np.sum(var_2nm)/5)
    
    count += 1

mean_14nmC2
std_14nmC2
x_pos = np.arange(5)
width = 0.28
aa_types = ["Negative", "Positive", "Polar", "Hydrophobic", "Aromatic"]
fig = plt.figure(figsize=(7,7))
fig.canvas.layout.width = '800px'
fig.canvas.layout.height = '700px'
plt.bar(x_pos, mean_09nmC2, width, yerr=std_09nmC2, ecolor='black',capsize=5, color='royalblue')
plt.bar(x_pos+width, mean_11nmC2, width, yerr=std_11nmC2, ecolor='black',capsize=5, color='c')
plt.bar(x_pos+(2*width), mean_14nmC2, width, yerr=std_14nmC2, ecolor='black',capsize=5, color='lightslategray')
plt.title(r'No. of PEG residues Rg Unrestrained, Config 2', fontsize=16)
plt.xticks(x_pos+width, labels=aa_types, fontsize=16)
plt.yticks(fontsize=16)
plt.ylim(0,45)
plt.legend(['Rg = 0.9 nm', 'Rg = 1.1 nm', 'Rg = 1.4 nm'], frameon=False, fontsize=16)
plt.ylabel(r'No. of PEG residues', fontsize=16)<jupyter_output><empty_output><jupyter_text>Protein/polymer contact map movie<jupyter_code>fig = plt.figure(figsize=(10,10))

# Set the axis and the plot titles pp

plt.title("BSA/PLGA contact map 2 nm restrained", fontsize=22, loc='left')
plt.xlabel("PLGA Residue No.", fontsize=22)
plt.ylabel("BSA Residue No.", fontsize=20)

 # Set the axis range 
plt.ylim(583, 0)
plt.xlim(0, 60)

# Plot bands for each chain 
BANDS = (
    (0, 20, "purple", "B"),
    (20, 40, "blue", "C"),
    (40, 60, "green", "D"),
)
    
text_y = 0.98 # Close to the top
for start, stop, color, band in BANDS:
    plt.axvspan(start, stop,color=color, alpha=0.15)
    text_x = middle_of_band(start,stop)
    plt.text(
        text_x,
        text_y,
        "PLGA chain " + band,
        color=color,
        fontsize=18,
        transform=fig.gca().transAxes,
        horizontalalignment='center',
        verticalalignment='center',
        style='italic',
    )
    
plt.text(0.94, 1, "Time [ns]:", fontsize=20, transform=fig.gca().transAxes, horizontalalignment='right', verticalalignment='bottom')

# Set tick label size
fig.gca().tick_params(axis='both', which='major', labelsize=20)

ims = []
for i in range(10000):
    data = trj_ppmap_2nm[i]
    im = plt.imshow(data, aspect='auto', cmap='Greys')
    t_sim = plt.text(1.03, 1, str(i/100), fontsize=20, transform=fig.gca().transAxes, horizontalalignment='right', verticalalignment='bottom')
    ims.append([im, t_sim])
    
ani = animation.ArtistAnimation(fig, ims, blit=True, repeat=False)
ani.save('2nm_res.mp4', writer='ffmpeg', fps=50, bitrate=100000)
#plt.tight_layout()
#plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13649477100613056, 'fraction_numerical': 0.07202101900983979, 'mean_word_length': 5.170216147488874, 'pattern_counts': {'":': 0, '<': 120, '<?xml version=': 0, '>': 120, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '30', 'dir': 'jupyter-structured-clean-dedup', 'id': '589357', 'license': 'MIT', 'n_tokens_mistral': 18277, 'n_tokens_neox': 16545, 'n_words': 3150, 'path': 'Prot_poly_analysis/C2_bpeg_Unrestrained_11_22_2019.ipynb', 'repo_name': 'UWPRG/Nance_Enzyme_Encap_MD', 'size': '38822'}