In [39]:

    

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import datetime
import os
from small_script.myFunctions import *



%matplotlib inline
%load_ext autoreload
%autoreload 2


    

    




The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

In [44]:

    

plt.rcParams['figure.figsize'] = [16.18033, 10]    #golden ratio
plt.rcParams['figure.facecolor'] = 'w'
plt.rcParams['figure.dpi'] = 100


    

In [3]:

    

data = pd.read_csv("/Users/weilu/Research/database/membrane_training_set/proteins-2019-05-01.csv")
data.pdbid = data.pdbid.apply(lambda x: x[2:-1])


    

In [4]:

    

data.head()


    

    
Out[4]:







  

    

      

      
id
      
ordering
      
family_name_cache
      
species_name_cache
      
membrane_name_cache
      
name
      
description
      
comments
      
pdbid
      
resolution
      
...
      
species_id
      
family_id
      
superfamily_id
      
classtype_id
      
type_id
      
secondary_representations_count
      
structure_subunits_count
      
citations_count
      
created_at
      
updated_at
    
  
  

    

      
0
      
1
      
2024.0
      
OmpA family
      
Escherichia coli
      
Gram-neg. outer
      
Outer membrane protein A (OMPA), disordered loops
      
NaN
      
OmpA is required for the action of colicins K ...
      
1qjp
      
1.65
      
...
      
9
      
34
      
26
      
2
      
1
      
3
      
1
      
2
      
2018-08-13 03:49:46 UTC
      
2018-09-21 18:14:03 UTC
    
    

      
1
      
2
      
2028.0
      
Enterobacterial Ail/Lom protein
      
Escherichia coli
      
Gram-neg. outer
      
Outer membrane protein X (OMPX)
      
NaN
      
OmpX from Escherichia coli promotes adhesion t...
      
1qj8
      
1.90
      
...
      
9
      
355
      
26
      
2
      
1
      
7
      
1
      
1
      
2018-08-13 03:49:46 UTC
      
2018-09-21 18:14:03 UTC
    
    

      
2
      
3
      
2033.0
      
Opacity porins
      
Neisseria meningitidis
      
Gram-neg. outer
      
Outer membrane protein NspA
      
NaN
      
Pathogenic Neisseria spp. possess a repertoire...
      
1p4t
      
2.55
      
...
      
24
      
337
      
235
      
2
      
1
      
0
      
1
      
0
      
2018-08-13 03:49:46 UTC
      
2018-09-21 18:14:03 UTC
    
    

      
3
      
4
      
1740.0
      
Influenza virus matrix protein 2
      
Influenza virus
      
Viral
      
M2 proton channel of Influenza A, closed state...
      
NaN
      
NaN
      
3lbw
      
1.65
      
...
      
51
      
263
      
185
      
11
      
1
      
3
      
4
      
0
      
2018-08-13 03:49:46 UTC
      
2018-10-02 17:42:36 UTC
    
    

      
4
      
5
      
2045.0
      
OM protease omptin, OMPT
      
Yersinia pestis
      
Gram-neg. outer
      
Plasminogen activator PLA (coagulase/fibrinoly...
      
NaN
      
NaN
      
2x55
      
1.85
      
...
      
299
      
36
      
27
      
2
      
1
      
2
      
1
      
0
      
2018-08-13 03:49:46 UTC
      
2018-09-21 18:14:03 UTC
    
  

5 rows × 31 columns

In [8]:

    

d = data.query("classtype_id == 1")
dd = d.groupby("superfamily_id").apply(pd.DataFrame.sample, 5, replace=True).reset_index(drop=True)


    

In [10]:

    

dd.shape


    

    
Out[10]:





(600, 31)

In [13]:

    

ddd = dd.drop_duplicates().reset_index(drop=True)


    

In [15]:

    

ddd.to_csv("/Users/weilu/Research/database/membrane_training_set/chosen_more_data.csv")


    

In [17]:

    

a = pd.read_csv("/Users/weilu/Research/database/membrane_training_set/chosen_more_data.csv", index_col=0)


    

In [19]:

    

a.shape


    

    
Out[19]:





(344, 31)

In [22]:

    

# print(a.pdbid.tolist())
pdb_list = a.pdbid.tolist()


    

In [34]:

    




    

In [27]:

    

pre = "/Users/weilu/Research/server/may_2019/four_body_helix_more_data/dompdb/"
filtered_list = []
for pdb in pdb_list:
    location = pre + f"{pdb}.pdb"
    a = get_inside_or_not_table(location)
    # ratio = sum(a)/len(a)
    ratio = sum(a)/(len(a)+1e-6)
    if ratio < 0.4:
        print("not good", pdb, ratio)
    else:
        filtered_list.append(pdb)
    # print(pdb, ratio)


    

    




not good 4l6v 0.32476319306527307
not good 6hqb 0.3249001327231689
not good 6btm 0.08144796343236216
not good 1kqf 0.0
not good 6f0k 0.091787439170109
not good 4ysx 0.0
not good 6adq 0.13461538418392505
not good 3j9u 0.0
not good 4b2q 0.0
not good 3j9t 0.0
not good 4u2q 0.08460471555533326
not good 5goa 0.051709027154043526
not good 6nr2 0.2149877147236023
not good 3udc 0.15355805185933313
not good 3t9n 0.1691729316948386
not good 5y4o 0.18834080633031028
not good 4y7k 0.21199999915200002
not good 4y7j 0.33195020609149295
not good 3jad 0.24629080045611038
not good 6his 0.20671834571907405
not good 5sxv 0.26298701213315906
not good 6dg7 0.21227621429085367
not good 5njy 0.27009646215403066
not good 2gif 0.2625968989703518
not good 5xap 0.35792349677879304
not good 6nzi 0.2660731945933994
not good 6dmo 0.26476578384443406
not good 6baj 0.2607421872453689
not good 6eti 0.2614840984779433
not good 2onk 0.0
not good 5tv4 0.22280701715297013
not good 4ksd 0.22504230099404202
/Users/weilu/Research/server/may_2019/four_body_helix_more_data/dompdb/3bo0.pdb (' ', 1, ' ')
not good 3bo0 0.0
/Users/weilu/Research/server/may_2019/four_body_helix_more_data/dompdb/5gae.pdb (' ', 1, ' ')
not good 5gae 0.0
not good 3j09 0.22406638973158174
not good 2c9m 0.2434607643425948
not good 4nab 0.28433734905501523
not good 4bbj 0.25602409599996373
not good 2yn9 0.23464249724608005
not good 6cxh 0.10994764369123654
not good 4phz 0.0
not good 3chx 0.11049723726381978
not good 1yew 0.10471204161070147
not good 2w8d 0.0
not good 1p49 0.0912408757459108
not good 5fgn 0.22947761151216864
not good 2mi2 0.31730768925665687
/Users/weilu/Research/server/may_2019/four_body_helix_more_data/dompdb/5xls.pdb (' ', 409, ' ')
not good 5xls 0.0
not good 2iub 0.10876132897655187
not good 2bbj 0.11976047868335186
not good 3jcf 0.11747850969203866
not good 3viv 0.0045662100248118264
not good 3bpp 0.0
not good 6haw 0.0
not good 6g94 0.0709459457062637
not good 4rfs 0.0
not good 5x3x 0.0
not good 1kn9 0.004273504255241435
not good 1b12 0.012552301202710036
not good 1t7d 0.0
not good 2leg 0.0
not good 2zup 0.010638297815753735
not good 2yvx 0.2737556554892406
not good 2zy9 0.29039812578361096
not good 4nyk 0.12801932336227217
not good 4fz1 0.1108374381506467
not good 5svl 0.1651651646691737
not good 3h9v 0.15805471076579117
not good 3j1z 0.38851351220096786
not good 6g9o 0.11699164329109799
not good 6el1 0.0626780624994927
not good 6ek7 0.024725274657348147
not good 2wcd 0.18947368354570637
not good 3blc 0.018315018247930335
not good 6gyb 0.018018017855693534
not good 3jqo 0.1507537680866645
not good 2bhv 0.010582010526021108
not good 6ezn 0.05707762544046204
not good 5fn3 0.03458646611340381
not good 4uis 0.03571428565354713
not good 5fn4 0.03458646611340381
not good 2ybb 0.0
not good 6nbq 0.0
not good 5xth 0.0
not good 4fyg 0.00807537011026195
not good 4jza 0.0
not good 5iji 0.18061673929243727
not good 5jeq 0.18061673929243727
not good 5jgp 0.18141592840081447
not good 5jef 0.20091324109172037
not good 3kds 0.0
not good 1lv7 0.003984063729147156
not good 3b8n 0.0
not good 3pjv 0.007812499938964845
not good 5ejz 0.3241758237305277
not good 4p00 0.3104395600131325
not good 4hg6 0.29852744270612125
not good 4p02 0.32005494461530914
not good 6qpc 0.37781629050638427
not good 6dz7 0.35017421541781496
not good 6qpb 0.29001367949382534
not good 6n26 0.25705329073024047
not good 6n25 0.25147928919680684
not good 6n28 0.25147928919680684
not good 4rdq 0.2295081960942399
not good 5nf8 0.0
not good 5syt 0.3963963955036117
not good 4mnd 0.3137254894271434
not good 5d92 0.34795321535686197
not good 5d91 0.36417910339051013
not good 5eke 0.17361111050829475
not good 5sy1 0.0
not good 5sv1 0.28251120949546543
not good 5sv0 0.29017857013313136
not good 5zfu 0.30222222087901235
not good 5zfv 0.29910714152184315
not good 5uph 0.0
not good 5xhq 0.3934426221445848
not good 5mm1 0.26011560618463697
not good 5mm0 0.29295774565364013
not good 5v7v 0.02120717777943364
not good 6bml 0.3972602726121224
not good 6bmm 0.3655172401189061
not good 6b85 0.19069767353163872
not good 6f2d 0.03940886680094154
not good 6c5w 0.1473684205355494
not good 6dt0 0.18220338905846023
not good 6d7w 0.14909090854876034
not good 6dnf 0.18666666583703703
not good 6bug 0.0
not good 6ds5 0.0437499997265625
not good 6mlu 0.07843137203639626
not good 6nt8 0.2721088426118747
not good 6nt5 0.2945205469365735
not good 6nt6 0.2877697831375188

In [28]:

    

len(filtered_list)


    

    
Out[28]:





209

In [29]:

    

pdb_list = filtered_list


    

In [35]:

    

pre = "/Users/weilu/Research/server/may_2019/four_body_helix_more_data/"
for pdb in pdb_list:
    toLocation = pre + f"dompdb/{pdb}.pdb"
    location = pre + f"cleaned_pdbs/{pdb}.pdb"
    extractTransmembrane(toLocation, location)


    

In [37]:

    

with open("/Users/weilu/Research/server/may_2019/four_body_helix_more_data/optimization/protein_list", "w") as out:
    for pdb in pdb_list:
        # print(pdb)
        out.write(pdb+"\n")


    

In [41]:

    

pre = "/Users/weilu/Research/server/may_2019/four_body_helix_more_data/"
for pdb in pdb_list:
    location = f"{pre}database/dompdb/{pdb}.pdb"
    toLocation = f"{pre}database/S20_seq/{pdb}.seq"
    seq,resseqs = getSeqFromPDB(location, considerGap=False)
    with open(toLocation, "w") as out:
        out.write(seq+'\n')


    

In [42]:

    

pre = '/Users/weilu/Research/server/may_2019/four_body_helix/optimization/'
location = pre + "gammas/protein_list_phi_six_letter_four_body_helix_docking3.5_6.5_5.0_10_True_gamma"
gamma_k1000 = np.loadtxt(location)

pre = '/Users/weilu/Research/server/may_2019/four_body_helix_repeat_n_decoys2000/optimization/'
location = pre + "gammas/protein_list_phi_six_letter_four_body_helix_docking3.5_6.5_5.0_10_True_gamma"
gamma_k2000 = np.loadtxt(location)

pre = '/Users/weilu/Research/server/may_2019/four_body_helix_repeat_n_decoys4000/optimization/'
location = pre + "gammas/protein_list_phi_six_letter_four_body_helix_docking3.5_6.5_5.0_10_True_gamma"
gamma_k4000 = np.loadtxt(location)

pre = '/Users/weilu/Research/server/may_2019/four_body_helix_repeat_n_decoys8000/optimization/'
location = pre + "gammas/protein_list_phi_six_letter_four_body_helix_docking3.5_6.5_5.0_10_True_gamma"
gamma_k8000 = np.loadtxt(location)

pre = '/Users/weilu/Research/server/may_2019/four_body_helix_repeat_n_decoys8000/optimization/'
location = pre + "gammas/protein_list_phi_six_letter_four_body_helix_docking3.5_6.5_5.0_10_True_gamma_filtered"
gamma_k8000_f = np.loadtxt(location)


    

In [47]:

    

pre = '/Users/weilu/Research/server/may_2019/four_body_helix_more_data/optimization/'
location = pre + "gammas/protein_list_phi_six_letter_four_body_helix_docking3.5_6.5_5.0_6_True_gamma"
gamma_more = np.loadtxt(location)


    

In [49]:

    

plt.plot(gamma_k1000, label="k1000")
plt.plot(gamma_k2000, label="k2000")
plt.plot(gamma_k4000, label="k4000")
plt.plot(gamma_more, label="gamma_more")
# plt.plot(gamma_k8000, label="k8000")
# plt.plot(gamma_k8000_centered, label="k8000_center")
plt.legend()


    

    
Out[49]:





<matplotlib.legend.Legend at 0x1a270207f0>






    

In [48]:

    

plt.plot(gamma_k1000, label="k1000")
plt.plot(gamma_k2000, label="k2000")
plt.plot(gamma_k4000, label="k4000")
plt.plot(gamma_more, label="gamma_more")
# plt.plot(gamma_k8000, label="k8000")
# plt.plot(gamma_k8000_centered, label="k8000_center")
plt.legend()
plt.ylim(-5,5)


    

    
Out[48]:





(-5, 5)






    

In [50]:

    

import glob
a = glob.glob("/Users/weilu/Research/server/may_2019/four_body_helix_more_data/optimization/data/*.dat")
data_all = []
for d in a:
    name = d.split("/")[-1].split("-")[0]
    tmp = pd.read_csv(d, sep=" ", names=["res_pair1_chain", "res_pair1_index_1", "res_pair1_index_2", 
                               "res_pair2_chain", "res_pair2_index_1", "res_pair2_index_2", 
                               "res_type",  "six_letter_string", "total_phi",
                              "d00", "d11", "d01", "d10"])
    data_all.append(tmp.assign(protein=name))
data = pd.concat(data_all)
data = data.reset_index(drop=True)
data["six_letter_string"] = data["six_letter_string"].astype(str).str.pad(width=6, fillchar="0")


    

In [52]:

    

data.head()


    

    
Out[52]:







  

    

      

      
res_pair1_chain
      
res_pair1_index_1
      
res_pair1_index_2
      
res_pair2_chain
      
res_pair2_index_1
      
res_pair2_index_2
      
res_type
      
six_letter_string
      
total_phi
      
d00
      
d11
      
d01
      
d10
      
protein
    
  
  

    

      
0
      
A
      
15
      
19
      
A
      
77
      
81
      
GAVF
      
110101
      
0.691777
      
8.218560
      
9.569242
      
6.418541
      
4.126931
      
2c3e
    
    

      
1
      
A
      
17
      
21
      
A
      
283
      
287
      
ASGL
      
120100
      
0.687627
      
10.557247
      
5.287187
      
6.421049
      
4.385104
      
2c3e
    
    

      
2
      
A
      
17
      
21
      
A
      
286
      
290
      
ASVY
      
020101
      
0.995874
      
8.894876
      
6.960624
      
5.734986
      
4.060895
      
2c3e
    
    

      
3
      
A
      
18
      
22
      
A
      
76
      
80
      
AKNY
      
010021
      
0.987386
      
8.097395
      
8.315829
      
4.241660
      
3.940909
      
2c3e
    
    

      
4
      
A
      
19
      
23
      
A
      
73
      
77
      
ATNV
      
020110
      
0.995154
      
9.297452
      
6.720541
      
4.126931
      
5.918079
      
2c3e
    
  

In [53]:

    

data["six_letter_string"].value_counts()


    

    
Out[53]:





020200    710
010300    505
030100    430
110200    427
120100    422
010201    413
020101    334
110101    295
210100    244
020110    242
010210    240
130000    210
220000    204
000400    200
120001    195
110110    176
000301    161
010111    157
010102    145
100300    144
100201    142
200200    135
011200    131
120010    128
040000    125
030001    124
021100    109
210001    108
000310    101
200101     98
         ... 
011002      8
100012      7
100030      7
101011      6
000022      6
202000      5
012100      5
000031      4
001012      4
101002      4
022000      3
012010      3
201010      2
012001      2
002101      2
001003      2
301000      2
002011      2
000004      2
000013      2
102100      2
002002      1
102010      1
002110      1
102001      1
002020      1
112000      1
000040      1
002200      1
001021      1
Name: six_letter_string, Length: 118, dtype: int64

In [ ]:

    

def getSeqFromPDB(location, considerGap=True):
    x = PDBParser().get_structure("x", location)
    seq = ""
    resseqs = []
    preResId = 0
    for res in x.get_residues():
        resId = res.get_id()[1]
        if considerGap and resId != preResId + 1:
            seq += " "
            resseqs.append(-1)
        seq += three_to_one(res.get_resname())
        resseqs.append(res.get_id()[1])
        preResId = resId
    return seq,resseqs