From 314a12386bfc770573e4ce709ac321788d0b7120 Mon Sep 17 00:00:00 2001 From: Thomas Forest Date: Wed, 18 Sep 2019 18:50:30 +0200 Subject: [PATCH] cleaner output --- src/__pycache__/atom.cpython-37.pyc | Bin 11068 -> 0 bytes src/__pycache__/geom.cpython-37.pyc | Bin 2211 -> 0 bytes src/__pycache__/pdb.cpython-37.pyc | Bin 3248 -> 0 bytes src/__pycache__/structure.cpython-37.pyc | Bin 1436 -> 0 bytes src/chem.py | 507 +++++++++++++++++++++++ src/{atom.py => chem.py~} | 33 +- src/dssp.py | 166 ++++---- src/geom.py~ | 66 --- src/pdb.py | 6 +- 9 files changed, 615 insertions(+), 163 deletions(-) delete mode 100644 src/__pycache__/atom.cpython-37.pyc delete mode 100644 src/__pycache__/geom.cpython-37.pyc delete mode 100644 src/__pycache__/pdb.cpython-37.pyc delete mode 100644 src/__pycache__/structure.cpython-37.pyc create mode 100644 src/chem.py rename src/{atom.py => chem.py~} (95%) delete mode 100644 src/geom.py~ diff --git a/src/__pycache__/atom.cpython-37.pyc b/src/__pycache__/atom.cpython-37.pyc deleted file mode 100644 index cf2b4a735c9fcaa901a814d851ca0bbe6d4e9616..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 11068 zcmds7TWlQHd7d*fd*wwGMN_0D%UjuT!_<{5+DQtRk*&K_)3hQfmbJ08UF{i>OD=cG zGeenTciY4?TDL#~wGX`*hARa{fabwZc_@k^E&9~7K+vb6m?sB)(rb~2LMY&V-+yL@ zyGtrH`q*8}Kj%JY{`>!*IXyB`bnv@%_5GWdE;`PCu+sf`D7=m<{18ApS~Z+9?`}g~ zb1pdA)BXpJ_P5=#i&{o!QOly{>734g;FNt`&_&!cdPooBp4B6I6!)Cor^j&5>-~Bh z_kup4CvY$7r}QN5L;9ee!hKjz>!)!a(TDV5+()&#>dZXz2)l6>W>l=cQe(0Eb5VF5 zS9lm8aLP(MWmhZhVyb61i!EiaL^k$MMXl8uXYj{+2G*Y^%1h8Z@Ku_?l@g2T`Tg+Az`a(A17|*0HsF zaZ%U9sKQ#v*9fx!j`GzJwY85gdTgh^mS-P+ymPbi&hsd=@x~RW?EP6{(@0h^AUmZ0cMJEjvK7hbzG1_v$`6@ zc@u;cTdFNp>rHkS*%k)tmFDeuwB8H@6V+QyytWPu-^}C)xc@jSs10ZGrogsQf_(&I z1p5ib2@U|nLo&GDtVTg7h_+yeefdhI-mFKJ3Txp6fa7M>&cCSm-x?8YcE>J&(3r^W zS1NI_Qdw>3jRwoZmCE`?wPAbAFjg~@1SD?{;KK~1j8{lB<=uQ%T571IKho)ST%iVd zGG|i&A#7IovWf8J65-1y!dFOyub2qmP$GQ8iSUi+gP`A3oVyZ)b-fYPxHtZ|Jzi&U zh3^2=$q3Hod|RD$Zl3LUWSV8C?JcW0LNXHCASPrV=Rt<$hh(xWN@SHB?s?qv!gPi8 z;p@;Yc(_^sw}?oDmYSqtoB>$-dM%jo zg+IzgGG16|Kw8YRT+?U6AVEl}XwpQ0YFOp*9~HPYmSE;smk&}sLiCV#1VwjH zp!T8Wp>iZ$NE6b7v>+Ww1LA*iqT>?(RD4u=RCxLseMCR2OZr#zQT?2LUVllypnp|= z8DoiR_W>mKA*2xqn^`D(r{+%qT)6;nmDQ`4%iftASIgdo`EL=R?ZW(xvUl+k+AhtX zD|_?j&@s;;Z=SnW_AXyU+vO{7mA$Lyu9UrN7YVLlm1{Rpdi&ITtWL%1N~{)Qbv9OK zVs$Q7-;C9TSY3?OrC4FW8?kybR_9}NIaXI=buCucV|61|Z$n#WK#O2AA%&2!oS{s= zY;Pw9$?y#D_dKqU;^U|>{D-|GYSi7D>Jng(3sOQ*&Qf0eAu3uySwpEmbXIa^T)S6c zy?}DRM_rt;PT;hBq4 z7I{8`+tErEpE2KC$;CRTLRvdE#3Iq1mwi5B$msSvK{m6Bv?Q6B- znWe~wjl!sDXF0P?qKb9jnk$x)--RG5RU3^`v=o%04bu!ucT8)wR4s*TL9KqLUMm4j zRCCADDb{d#*!3g6!y+#HKKuLgGg$I}KS_<_j=5X=yO>VLr^qIs^P#f9k5EVv z?Y80mtYpk5Z8rk{Z!4{)9H8-&$bH~New2x_KpxEavKaH_{I;?!9M{{r?RLj%b+E5z z_X^wYwo4+n;c~6`=1+0=7eN}_GrU6)+B+3jZb6)b*gbRV(Rpm$yaW&@IRDPw*JlcX za&wY>bZy`}28$*DVn^-pA)K{rd3J_a3g9D ziTm)(Wr$7z9=4~8(=Z_IEGRqB2D9XKd~z1N_1vP{$**6gmLvy78EDHP%=h)Mz>#TZ z8~#H(PrrT6>FGz09hLOMxS<~^0;&VTjMaYDokzbzkK&>j?^_UK^oXiR)Tc;wsySnS zt!k=_ilA&ynBTyt*snuCV}F@2b33}X7Qh*6Ms@Qh*WeIynP7q7Er722wYuONY^L%R zwHs$oo4Q^Mpe`tKhN2zkZcAGjXKf)KgI;Jsb0bu;Es6A>tH6;ol+k$DhH6RofKeQx z?oAf0Vsv;1fMFK@MRgdkViR2zIxaX}q<6C`y7H<<5_^b-Yin9(^Cnc^>Zcl|Rd-+RDw{Bym2rT;U zVRusX+fW&oVvvTO%qti)D5T2SV$d_q6njm<9uzdxjm2KtCu-5vJ$+E~jD-AUcg4jT zK9UF?F0?L!0v~&DeBiQZA9&8#kU4Sd;!GG-kv#BVfXr3w=sfHe-4DaiD;}~uAlQj> z5UWaaV>R(Ug#86W<>5xFR&69qTOO&k5KY&TYMgCUHA3n^G#1I`w(Ctvh7D6c2+Vn*&d`?ubJr`n~}&1KFj+OM}m%F1?p(WQ!p1f zR@@GEtzPzOt#Bqc$b?i?<=mOKD%5S=!mPQ(sRIfmC(c`~^%L}^{zM*#9Hx+qpd-hy zuNS3Jt(D@FOqzNVdpdd_rH-EX59T9U1p3r;1jVI;Q^ z*@qYIn-6f$tSIw=&Vqf}rpym;&!O)JI^W4O5$AS%U1V=toV}SLYas;j<{WXu%+RAB z13Xs0R-ga;^UoivU$+(mexJF3j+xAk^B0eRk82fOUkbFTHg;&Vi|6A)ZK;kS_2`}u z0ZF*}VsnrbBP2yS@L%vW$r?K1;87Qkg4ZWqXk=Vjk8#h|Q@w;q=ceQ#pQPmZTT}*g zh-;>_+wq#G>B4pOZeyRbssn`eK|+X8 zMP~(lpJ9=d3^0B!$v{GY1Rw*;lLMwOnYl1CY7TM>R91$`8=G%q4#Bz!v?%7c2&g|z znc%kxXvqy<)GQP1h^&+(L-#Z~xUhu6FxNx)PvrOKNAf50`|@MBKV>+kMx61-+$5)q zIXUWy->)Ex;`_s!|E~xM3F>Swz z#_a3@UOM}B&D(3dZIkC)5qR=G@NE06o;~zoVxYZmJol5#mv{JIpUddLzJDJ+``*(v zOI)`XKYR4=IsZxXzVN5CJjrPLcM zNJ}?fzzUu(kpd4ZCeqRrl%?;AvP&%#9SB~cmVSl;NsnZlxA_lk4=Tr|-u{W{IV9d} zwlr%t!F+HU7LgI%JkN|KLeOsCUE9Z8zwM#*fknF8&LBij+a4+<3h3)dd#ukgFh#<{ zH$O>cQ@zN#R%8(LGaZ<=X!m0>LUP8F{K1 z@Dd{=7)I3K;ehk9$@?WyQa?!<7d7z|(<0cRv~Skfuui~;ZqL^n-}QR4@ACbIF`@f@ zh~GUEq26cSQ#pDS<9W>U0j0uu2nubh&W^EHn^Qa=yRUAJcKrKl6B$?3oSXmMVX6!5P z8L8o%&8~%UEX!I@$S*<8xcV<>3cyLfpp@B>(~NGhJDH!h_08p+oUvmMN{akeg^RMbW>`C7cxE^taTyiyj^eKtd!WoB0mnfIoKE!PPsu1YZQW$6kv4pb-q7OFk~n zQGRQ!)wPpNGAWnMI$K8=%7jZn5QWF)B;F5CqR<_B{{OZOvw_vkU4njp|FcAawA{C;!bWH_ReY~h742Gv!b6ESdXbm{0R@chDJF>3gwu=jcBWn zP=HNToA=w)v>!t6zAvRMgwdcGqcT9~q;1zs(1Gu5p?;QRMSm;iK-MAUb~KY0NeL3I zA?1)7oF8)4q>m-^Z5NcYJ^@)NLbvB(-wOYKWyK~2eiQgN6@pU)q=~5Qt$ZSwvkS-q z(%OhZjx%ef-i%!voJCwJISV`#VJP;*p{JQ2xLhptb-N8u+z@KMrtn{llJDl+29 z)??MM{R2WrtCNJF?ds3_PLsZDg^zd+bo8RXSw_$s3zW*+=1qcBo zq~d2za;lv*lbd&tD#{||#!>I!SP5nLX|8#fC1ea>q2T-F4mxwNN=y}@Z&Uh^;T-7m z(Kjc3v>xfa^*QkKBymA|vFHIK$a@OiDTRjK|08Bn;HDG;-ge*B^=3Mh9P%B-MblY=izoq4q1R;mCuix55jzVw7R0#8Bypj+` zzD+RSVU;QuqWE2w%LF72K0UCmvn1CmFUL7sc%@sgbTRK^Y@Ct4a?uW9mcLxEW#(Xi z3B8g*3B5igy-HZHpI+lQKH=v&JZ#}L<(dyMrtf1WiZWG<*bh2yJRD)9P6rh(7C#or znxkY?IqiabEvH>cnDJPZAi$E6r@4b0%9odyGk?I+DT01Ejt_{&s`S0XKk^ksuEA=I z4$L26=H6S~LeCSodW;d#6s(VIJ!g)|sn5=kXMJS zQcNGKlaJM$`9l<+h=SMAxhD$7EehmtCy~0J5(A=pfdR(Jhe0yjUGqIYVxQ#T;{wj_ zHfop`$QN>wa@Y;!!h*)h(E-Sdy2<&{8QdNbdg;Nz(lwOP~YrV6vw|GC883&Bk z7s8*xAso5$*Q_{zIOW10zy*FYYddMGO6+Lfyf+_j-h1=jZ;zYJ2En&k|MBIYE+Kzu zVev6ApP|dAAe?YICcSE;J;oWgej+`K+uXrubC=gJI^5%Rj4p5R6^u2$%9|KHzQ)%v z);WDn+8dY9NKQh1gZ^+cvCto*%SRvu<@63IaOWMeMP>{#GcqV-3$kToP{?+{3ahY( z*32m!=(kSfEQz`b zG?3jVSvX0eRLO2GvR9GnhuK7^t`uQ+7-dN(pO))qQK&M}KNVTh&oilgktVQo!?Xq! z9_B)Akp8qy4!Uu;14fYxTDi;Frgj{TY3{b>;dI6SLT9wl0h~K1^hH7|YqQc$0yQda z`K?fp=vir}nMkhO5;g2ASLUYbHEp5PZWH;g zP_U?=hzV0jvU57KFj9=D8->a_r=a#ZZINvAQ_pdz2gqlk-O?3)n! zNd)}j9sdhpl0S&EXH&>&=h>6vhL7)Vz=K>x+*g^O=AA%@VA?h-R~psm4XuMneHTTk zf8fs((i{;{E+PqxjME{O!YoYzh)|wqF@Uqfoi2T*H)jV^|8o(C{?RxLVjgwO$f&~i%?JD)!{X$YS+cw^{8wbhCTd&h?xiA<;LtC^Y zQXY=Fwjdu0yG0NESo9F|Qb5kFe@}0`_9UQ(o{Jv(dn3vAx+zL&9%tU0d32sy<18_n{N<2ewVO%gQ)oMIT2qE6^Za-%J|4 zP8h3ny(^QWQ2C9dC)Ii?8|z*1`C9ihD;;WXtT4+Vn=BB=DCr|5}V137P8cC;{^kUI3KHK?p_wg2SvwVsTYlQvc#;1GRJCoHx58Epd8Ak>$ zKcSFhCIV42WgcJ3oKC0$5vU-uC|)`)V?u^xMcU7?0cDKSH$WJ-c-biPDmMN`GxNMZ ze=`i@Wmm~$V`T+RN8z9F^^gsiKc@VUogs72%vb#6&xnh$$zJjO-e7FHrw1uZlE6#@Q|F z6_?&1cT^(IATR7?w!)4|En~tjzs4`OutS*7j`+|~GXq0wqTIEY*ttR9p`E;?|I4q? zPkO%me|k{l6YG!fAy@zFtK%;YOnlEc=WjBx)I!|kKU765pKG6ovja|6X@FR%>lF8_ zjlG>GUVeLPZ)5Mt=o93`s1B0x?U614O3YD-2<6J79KA{89+9sT zStIf;NX^ddKqVbN5?MiaH263&O#N}M(F~%P#6{wz;jtfoo>iha4W){b81N&)tkg&n zDWW(~VS4yo{rv66{m<%QC|s4efe>!gXW>;=OjRIN z`Z7|@ESD-7bu()N6Eic4RhBzwCxOcHK^G`2vRo94@I=xEN;-7pH$gmmcdRkXbtNsD zffN3f6E#}{gQRlxiYn}MeQk(m1CC>Wft-gVvxu4T=r&*!2r%>lib?tigz-Fn7JgM; z;fvfstKe7U^SaIRHw@s3S>?<4y~QnqdU*4Sv24y8j?tgHi7Q(B6|B%%MNPvB<%`0} zO3A0sPP5)$+^X*ZN*`};J^t9;{>A3*&WGUCxYSuw0HC#+=EC zJ2(>=%F;8LDQD~8y&@lIAnxzs0&pA=7#3DJapfpBm2_emJaaAFA7GUaEQUr3KwAXl6;>E%QPa{21H_dl5H~lxGl94q z=FS0Z?&5&^;yF8RKG)VC**N58>d3d#PSWn14~`pNA&%6zYWjI9M>m#B*dniRm5BXU zsY+xh*1Adhf{I|@{h|9}|9$s1du^QU_esqvO`3d zdcWp|Z4}}EQupGPzWtiis6)v%q8LlcHk1qDj{F9Z`$TF)uBuJ0Vpv;}5ApJ<2T>fQ zO}`!V;zl!y4>PmT5jc61KzMmvlfs7ba!nrviJ|u1GvU5B6RA)Ff_zlJ3Oy?gei?eD zF>d$-E0S<8GXV%$j(mnv@=<=9jM4S4GyK{>LIU|cB4l|je_ixhaiSM^{-5eBmKIdRm=68){@C*$C1Ox;K6F|j$T8d zt0JO}rlyPNZFB}~h+ESo$K3&KXk5}Q!bH!yanj%_nv-;bXwa0ni8Dzji2MP_{1RbT|-o1e)?obH7Sa S<~$tR^dm}M?U7?ZpYb2k)&)NR diff --git a/src/__pycache__/structure.cpython-37.pyc b/src/__pycache__/structure.cpython-37.pyc deleted file mode 100644 index ee79b270cecd5469277e267eddedf020446bea25..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 1436 zcmbW1L2uJA6vyqvNz-)Q#)NOz|1Z zyOHMs>0lO~vbRiyN?b7|E`wZ_zBBK*A$-} zXWa3p1>C*`?8$LZfXKIAQhdqd}G5n2Bll@qBVVuk@rLJGGNH^P;VE-Pp;(%rA5!i+O9Bi8aN< z_M09omM)R9tjk79#IQaKzS`U5bp2)mr7a9RXc;fCFQ4fI_qXkLScM9 zH_QOP(E+ij-y`8(-x`=WZB1Pu6dJ}2t6dTut98(DNoU*<3*C?a*RokIwBOW2%R+-> z%oqj~|4tWmjHp{oG9L5Q7<=F3VTLqga7G3z|6lqR&0?*`n*;{>)}_Zsrc`ITK-y>+ zH^o0A(HBqGg|9=h0f7jmSPRAX8R>|;zM=r39165>>4$F$n^8QzJ|0}8W`fWHVqXNeQlv^mcsI^@MN59&oO p&h)(M&KM`;Ng6D^k)h8OX}aWP_Yd$2*T)}e!s!XN#csTt{shF0{bc|E diff --git a/src/chem.py b/src/chem.py new file mode 100644 index 0000000..37bc8ee --- /dev/null +++ b/src/chem.py @@ -0,0 +1,507 @@ +import math +#from structure import * +from geom import * + +class Atom: + + def dist_atoms(self, atom2): + return(math.sqrt((self.coord_x-atom2.coord_x)**2 + + (self.coord_y-atom2.coord_y)**2 + + (self.coord_z-atom2.coord_z)**2)) + + + def __init__(self, atom_id, atom_name, res_name, chain_id, + res_seq_nb, insertion_code, coordinates): + self.atom_id = atom_id + self.atom_name = atom_name + self.res_name = res_name + self.chain_id = chain_id + self.res_seq_nb = res_seq_nb + self.insertion_code = insertion_code + self.coord_x = coordinates[0] + self.coord_y = coordinates[1] + self.coord_z = coordinates[2] + self.coords = coordinates + +class Structure: + def __init__(self, res): + self.res = res + +class Turn(Structure): + + def __init__(self, turn_type, res): + self.turn_type = turn_type + Structure.res = res + +class Bridge(Structure): + + def __init__(self, bridge_type, res_num, res_partner, indices): + self.bridge_type = bridge_type + self.res_num = res_num + self.res_partner = res_partner + Structure.res = res_num + self.i = indices[0] + self.j = indices[1] + +class Helix(Structure): + + def __init__(self, residues, res_num, helix_type): + self.residues = residues + self.res_num = res_num + Structure.res = res_num + self.helix_type = helix_type + + +class Residue: + def __init__(self, atoms_list, indice): + self.atoms = {} + for atom in atoms_list: + self.atoms[atom.atom_name] = atom + self.resid = atom.res_seq_nb + self.res_name = atom.res_name + self.res_letter = self.get_amino_letter(atom.res_name) + self.chain_id = atom.chain_id + self.insertion_code = atom.insertion_code + self.indice = indice + + def get_amino_letter(self, res_name): + code3 = ['ALA', 'ARG', 'ASN', 'ASP', 'CYS', 'GLU', 'GLN', 'GLY', + 'HIS', 'ILE', 'LEU', 'LYS', 'MET', 'PHE', 'PRO', 'SER', + 'THR', 'TRP', 'TYR', 'VAL'] + + + code1 = ['A','R','N','D','C','E','Q','G','H','I','L','K','M','F','P', + 'S','T','W','Y','V'] + + return code1[code3.index(res_name)] + + def h_bond(self, res2): + if("H" not in res2.atoms.keys()): + return(False) + # dimensionnal factor, in kcal/mole + f = 332 + # partial charges + q1 = 0.42 + q2 = 0.20 + # distance between O-N atoms, in angströms + r_ON = self.atoms["O"].dist_atoms(res2.atoms["N"]) + # distance between C-H atoms, in angströms + r_CH = self.atoms["C"].dist_atoms(res2.atoms["H"]) + # distance between O-H atoms, in angströms + r_OH = self.atoms["O"].dist_atoms(res2.atoms["H"]) + # distance between C-N atoms, in angströms + r_CN = self.atoms["C"].dist_atoms(res2.atoms["N"]) + # electrostatic interaction energy, in kcal/mole + E = q1*q2*(1/r_ON + 1/r_CH - 1/r_OH - 1/r_CN)*f + return(E) + + def get_turns(self, residues): + """ + Get all the turns from a specific residue. + """ + turns = {} + i = residues.index(self) + k = 0 + for j in range(3,6): + if(i+j=2 and i 70): + return [angle, 'S'] + return [angle, ''] + else: + return [360.0, ''] + + def get_bridges(self, residues): + bridges = {} + bridge = {} + strongest_bridge = {} + i = residues.index(self) + if(i >= 1 and i < len(residues)-4): + E_min = 0 + for j in range(i+2,len(residues)-1): + # select triplet with the minimal energy + + if(residues[i-1].h_bond(residues[j])<-0.5 + and residues[j].h_bond(residues[i+1])<-0.5): + bridge = {'res1':residues[i-1].h_bond(residues[j]), + 'res2':residues[j].h_bond(residues[i+1]), + 'ipos':residues[i].resid, + 'jpos':residues[j].resid, + 'btype':"para"} + + if(residues[j-1].h_bond(residues[i])<-0.5 + and residues[i].h_bond(residues[j+1])<-0.5): + bridge = {'res1':residues[j-1].h_bond(residues[i]), + 'res2':residues[i].h_bond(residues[j+1]), + 'ipos':residues[i].resid, + 'jpos':residues[j].resid, + 'btype':"para"} + + if(residues[i].h_bond(residues[j])<-0.5 + and residues[j].h_bond(residues[i])<-0.5): + bridge = {'res1':residues[i].h_bond(residues[j]), + 'res2':residues[j].h_bond(residues[i]), + 'ipos':residues[i].resid, + 'jpos':residues[j].resid, + 'btype':"anti"} + + if(residues[i-1].h_bond(residues[j+1])<-0.5 + and residues[j-1].h_bond(residues[i+1])<-0.5): + bridge = {'res1':residues[i-1].h_bond(residues[j+1]), + 'res2':residues[j-1].h_bond(residues[i+1]), + 'ipos':residues[i].resid, + 'jpos':residues[j].resid, + 'btype':"anti"} + + if(bridge): + if(bridge['res1']+bridge['res2']0): + return(bridges[coord_bridge[0]]) + else: + return(False) + + def get_helix(self, residues): + """ + Return if there is an helix at a given residue, + as well as its type. + """ + i = residues.index(self) + # if there are no turns or it is the first residue, skip + if i == 0: + return False + + if(self.get_turns(residues) and residues[i-1].get_turns(residues)): + #print(self.get_turns(residues).turn_type,"- HELIX at", residues[i].indice) + return(self.get_turns(residues).turn_type, residues[i].indice) + return(False) + + def get_ladder(self, residues): + #ladders = {} + i = residues.index(self) + if i != 0: + if self.get_bridges(residues): + if (residues[i-1].get_bridges(residues)): + local_bridge = self.get_bridges(residues) + consec_bridge = residues[i-1].get_bridges(residues) + if local_bridge.bridge_type == consec_bridge.bridge_type: + #print("ladder", consec_bridge.res_num, local_bridge.res_num) + ladder = {'start':consec_bridge.res_num, + 'end':local_bridge.res_num, + 'bridges':[consec_bridge, local_bridge]} + return ladder + return False + + def get_tco(self, residues): + i = residues.index(self) + if(i!=0): + res2 = residues[i-1] + CO_res1 = vector_from_pos(self.atoms["C"].coords, + self.atoms["O"].coords) + CO_res2 = vector_from_pos(res2.atoms["C"].coords, + res2.atoms["O"].coords) + angle = vector_angles(CO_res1, CO_res2) + else: + angle = math.pi/2 + return(math.cos(angle)) + + def get_chirality(self, residues): + i = residues.index(self) + if (i >=1 and i < len(residues)-2): + chirality = {} + angle = calc_dihedral(residues[i-1].atoms["CA"].coords, + residues[i].atoms["CA"].coords, + residues[i+1].atoms["CA"].coords, + residues[i+2].atoms["CA"].coords) + + if(angle>0 and angle<=180): + sign="+" + + if(angle<=0 and angle > -180): + sign="-" + + else: + angle = 360.0 + sign = '' + + return [angle, sign] + + def get_phi_psi(self, residues): + i = residues.index(self) + if(i==0): + phi = 360.0 + else: + phi = calc_dihedral(residues[i-1].atoms["C"].coords, + residues[i].atoms["N"].coords, + residues[i].atoms["CA"].coords, + residues[i].atoms["C"].coords) + if(i==len(residues)-1): + psi = 360.0 + else: + psi = calc_dihedral(residues[i].atoms["N"].coords, + residues[i].atoms["CA"].coords, + residues[i].atoms["C"].coords, + residues[i+1].atoms["N"].coords) + + return((phi, psi)) +def get_bridges(residues): + bridges = {} + bridge = {} + strongest_bridge = {} + for i in range(1,len(residues)-4): + E_min = 0 + for j in range(i+2,len(residues)-1): + # select triplet with the minimal energy + + if(residues[i-1].h_bond(residues[j])<-0.5 + and residues[j].h_bond(residues[i+1])<-0.5): + bridge = {'res1':residues[i-1].h_bond(residues[j]), + 'res2':residues[j].h_bond(residues[i+1]), + 'ipos':residues[i].resid, + 'jpos':residues[j].resid, + 'i':residues[i].indice, + 'j':residues[j].indice, + 'btype':"para"} + + if(residues[j-1].h_bond(residues[i])<-0.5 + and residues[i].h_bond(residues[j+1])<-0.5): + bridge = {'res1':residues[j-1].h_bond(residues[i]), + 'res2':residues[i].h_bond(residues[j+1]), + 'ipos':residues[i].resid, + 'jpos':residues[j].resid, + 'i':residues[i].indice, + 'j':residues[j].indice, + 'btype':"para"} + + if(residues[i].h_bond(residues[j])<-0.5 + and residues[j].h_bond(residues[i])<-0.5): + bridge = {'res1':residues[i].h_bond(residues[j]), + 'res2':residues[j].h_bond(residues[i]), + 'ipos':residues[i].resid, + 'jpos':residues[j].resid, + 'i':residues[i].indice, + 'j':residues[j].indice, + 'btype':"anti"} + + if(residues[i-1].h_bond(residues[j+1])<-0.5 + and residues[j-1].h_bond(residues[i+1])<-0.5): + bridge = {'res1':residues[i-1].h_bond(residues[j+1]), + 'res2':residues[j-1].h_bond(residues[i+1]), + 'ipos':residues[i].resid, + 'jpos':residues[j].resid, + 'i':residues[i].indice, + 'j':residues[j].indice, + 'btype':"anti"} + + if(bridge): + if(bridge['res1']+bridge['res2']0): + return(bridges) + else: + return(False) + +def get_ladders(bridges, residues): + ladders = {} + i = 1 + while i < len(residues): + k = 1 + if i in bridges.keys(): + temp_bridges = [bridges[i]] + while ((i+k in bridges.keys()) and + (bridges[i].bridge_type == bridges[i+k].bridge_type)): + temp_bridges.append(bridges[i+k]) + k+=1 + if k>1: + #print("ladder", bridges[i].res_num, bridges[i+k-1].res_num) + ladders[i] = {'start':bridges[i].res_num, + 'end':bridges[i+k-1].res_num, + 'bridges':temp_bridges, + 'i':i, + 'j':i+k-1} + i+=k-1 + else: + i+=1 + return ladders + +def connected_ladders(ladd_1, ladd_2): + links = [] + for bridge in ladd_1['bridges']: + if bridge.res_partner in res_list(ladd_2): + return ladd_2 + + return False + +def connected_ladders2(ladd_1, ladd_2): + links = [] + for bridge in ladd_1['bridges']: + if bridge.res_partner in res_list(ladd_2): + return([ladd_1['i'], ladd_1['j'], bridge.i, bridge.j, + ladd_2['i'], ladd_2['j']]) + #return ladd_2 + + return False + + +def get_sheets(ladders): + """ + Bridges between ladders. + Check if 1 bridge between one ladder and one or more other ladders. + Iterate over all residues of one ladder and check if bridge with other residues + of the other ladders. + """ + ladds = [ elem for elem in ladders.values() ] + sheets = {} + + corresp = {} + for ladd1 in ladds: + for ladd2 in ladds: + if connected_ladders(ladd1, ladd2)!=False: + corresp_list = [ elem for elem in corresp.keys() ] + if ladd1['i'] not in corresp_list and ladd2['i'] not in corresp_list: + ind = len(sheets.keys()) + sheets[ind] = [] + sheets[ind].append(ladd1) + sheets[ind].append(ladd2) + corresp[ladd1['i']] = ind + corresp[ladd2['i']] = ind + elif ladd2 not in corresp_list and ladd1 in corresp_list: + sheets[corresp[ladd1['i']]].append(ladd2) + corresp[ladd2['i']] = corresp[ladd1['i']] + elif ladd1 not in corresp_list and ladd2 in corresp_list: + sheets[corresp[ladd2['i']]].append(ladd1) + corresp[ladd1['i']] = corresp[ladd2['i']] + + return sheets + +def get_sheets2(ladders): + """ + Bridges between ladders. + Check if 1 bridge between one ladder and one or more other ladders. + Iterate over all residues of one ladder and check if bridge with other residues + of the other ladders. + """ + sheets = {} + for ladder in ladders: + for ladd2 in ladders: + if connected_ladders(ladders[ladder], ladders[ladd2]): + bridge_i = connected_ladders2(ladders[ladder], ladders[ladd2])[2] + bridge_j = connected_ladders2(ladders[ladder], ladders[ladd2])[3] + print("ladder",ladders[ladder]['i'], ladders[ladder]['j'],"bridge",bridge_i, bridge_j, + "ladder 2",ladders[ladd2]['i'], ladders[ladd2]['j']) + + +def res_list(ladder): + # TODO : method in ladder class + l=[] + for i in range(ladder['i'], ladder['j']): + l.append(i) + return(l) + +def build_turns_patterns(residues): + turns_3 = {} + turns_4 = {} + turns_5 = {} + for i,res in enumerate(residues): + turn = residues[i].get_turns(residues) + if(turn): + for k in range(turn.turn_type): + if turn.turn_type == 3: + turns_3[i+1+k] = turn.turn_type + if turn.turn_type == 4: + turns_4[i+1+k] = turn.turn_type + if turn.turn_type == 5: + turns_5[i+1+k] = turn.turn_type + return[turns_3, turns_4, turns_5] + +def build_helix_patterns(residues): + helix_3 = {} + helix_4 = {} + helix_5 = {} + for i,res in enumerate(residues): + helix = residues[i].get_helix(residues) + if(helix): + helix_type = residues[i].get_helix(residues)[0] + helix_pos = residues[i].get_helix(residues)[1] + #print("TYPE", helix_type) + for k in range(helix_type): + if helix_type == 3: + helix_3[i+1+k] = "G" + if helix_type == 4: + helix_4[i+1+k] = "H" + if helix_type == 5: + helix_5[i+1+k] = "I" + #print(helix_3) + return[helix_3, helix_4, helix_5] + +def print_helix_pattern(residues, res, helix): + i = residues.index(res)+1 + if i in helix.keys(): + return (helix[i]) + else: + return(' ') + +def print_turn_pattern(residues, res, turns): + i = residues.index(res)+1 + if i in turns.keys() and not i-1 in turns.keys(): + return(">") + elif i in turns.keys() and i-1 in turns.keys(): + return(turns[i]) + elif i not in turns.keys() and i-1 in turns.keys(): + return("<") + else: + return(' ') + +def raw_ladders_print(ladders, ): + for ladd1 in ladders: + ladd_1 = ladders[ladd1] + for ladd2 in ladders: + ladd_2 = ladders[ladd2] + for bridge in ladd_1['bridges']: + if bridge.j in res_list(ladd_2): + print("ladder", ladd_1['i'],"-",ladd_1['j'], "|", bridge.i, "...", bridge.j, "| ladder", ladd_2['i'], ladd_2['j']) + diff --git a/src/atom.py b/src/chem.py~ similarity index 95% rename from src/atom.py rename to src/chem.py~ index c3db1b2..5b9c05a 100644 --- a/src/atom.py +++ b/src/chem.py~ @@ -1,5 +1,7 @@ import math -from structure import * +#from structure import * +from geom import * + class Atom: def dist_atoms(self, atom2): @@ -21,6 +23,35 @@ class Atom: self.coord_z = coordinates[2] self.coords = coordinates +class Structure: + def __init__(self, res): + self.res = res + +class Turn(Structure): + + def __init__(self, turn_type, res): + self.turn_type = turn_type + Structure.res = res + +class Bridge(Structure): + + def __init__(self, bridge_type, res_num, res_partner, indices): + self.bridge_type = bridge_type + self.res_num = res_num + self.res_partner = res_partner + Structure.res = res_num + self.i = indices[0] + self.j = indices[1] + +class Helix(Structure): + + def __init__(self, residues, res_num, helix_type): + self.residues = residues + self.res_num = res_num + Structure.res = res_num + self.helix_type = helix_type + + class Residue: def __init__(self, atoms_list, indice): self.atoms = {} diff --git a/src/dssp.py b/src/dssp.py index 3f1c9a2..84dbb3c 100644 --- a/src/dssp.py +++ b/src/dssp.py @@ -1,104 +1,86 @@ import pdb -import atom import sys from structure import * -from atom import * +from chem import * +import argparse -def print_dssp(): +from datetime import date - empt_line=" " * 145 +def print_dssp(pdb_file, verbose=False, raw_ladd=False): + # gets export date + export_date = date.today().strftime("%d/%m/%Y") + + ################### OUTPUT #################### + # print DSSP-style formatted header of PDB file + print("################### DSSP OUTPUT ####################\n####", + "Exported on",export_date, + "\n####################################################") + for elem in pdb_file.get_header() : + print(pdb_file.get_header()[elem], end="") + # Get preliminary data for print loop + residues = pdb_file.residues + bridges = get_bridges(residues) + ladders = get_ladders(bridges, residues) + sheets = get_sheets(ladders) + helix = build_helix_patterns(residues) + # iterating over residues + # print header of table + print(' # RESIDUE AA STRUCTURE BP1 TCO KAPPA ALPHA ' \ + 'PHI PSI X-CA Y-CA Z-CA') + for i,res in enumerate(residues): + #res.get_turns(residues, turns) + kappa = res.get_bends(residues)[0] + bend_symbol = res.get_bends(residues)[1] + t_co = res.get_tco(residues) + alpha = res.get_chirality(residues)[0] + chirality = res.get_chirality(residues)[1] + phi = res.get_phi_psi(residues)[0] + psi = res.get_phi_psi(residues)[1] + x_ca = res.atoms["CA"].coord_x + y_ca = res.atoms["CA"].coord_y + z_ca = res.atoms["CA"].coord_z + turns = build_turns_patterns(residues) + turn_3 = print_turn_pattern(residues, res, turns[0]) + turn_4 = print_turn_pattern(residues, res, turns[1]) + turn_5 = print_turn_pattern(residues, res, turns[2]) + helix_3 = print_helix_pattern(residues, res, helix[0]) + helix_4 = print_helix_pattern(residues, res, helix[1]) + helix_5 = print_helix_pattern(residues, res, helix[2]) + if i in bridges.keys(): + bp1 = bridges[i].j + else: + bp1 = 0 + + print('{0: >5}{1: >5}{2: >2}{3: >2}{4: >2}{5: >2}{6: >2}{7: >2}'\ + '{8: >2}{9: >2}{10: >2}{11: >1}{12: >5}{13: >7}{14: >7}' \ + '{15: >7}{16: >7}{17: >7}{18: >7}{19: >7}'\ + '{20: >7}'.format(i+1, res.resid, + res.chain_id, + res.res_letter, + helix_3, helix_4, + helix_5, turn_3, + turn_4, turn_5, bend_symbol, + chirality, bp1,round(t_co, 3), + round(kappa,1), round(alpha, 1), + round(phi, 1), round(psi,1), + round(x_ca,1), round(y_ca,1), + round(z_ca,1))) + if(raw_ladd): + raw_ladders_print(ladders) + +parser = argparse.ArgumentParser() if(len(sys.argv)<2): print("Not enough arguments! Run with --help to learn more about proper" "call structure and parameters.") else: + parser.add_argument("file", type=str, help="input file in PDB format") + parser.add_argument("-v","--verbose", help="add ladders verbose output", + action="store_true") + args = parser.parse_args() pdb_file = pdb.PDBFile(sys.argv[1]) - - #turns = get_turns(pdb_file.residues) - #get_helix(pdb_file.residues, turns) - #get_bends(pdb_file.residues) - # get_chirality(pdb_file.residues) - - - # print("NBRIDGES",len(bridges)) - - # bridges = get_bridges(pdb_file.residues) - - # get_bonds(pdb_file.residues) - - - # e_min = 0 - # for res in residues: - # ene = residues[1].h_bond(res) - # if ene <= e_min: - # e_min = ene - # best_res = res - #print(residues[1].resid, best_res.resid,e_min) - - #get_phi_psi(residues) - #print(residues[2].atoms, residues[0].resid, residues[1].resid) - #print(get_TCO(residues[2],residues[3])) - - -# for i,ladder in enumerate(ladders.values()): -# print(chr(65+i)) - -# for ind, sheet in sheets.items(): -# print(ind, sheet) -# for ladder in sheet: -# print(ladder['start'], ladder['end']) - - -################### OUTPUT #################### -# print DSSP-style formatted header of PDB file -for elem in pdb_file.get_header() : - print(pdb_file.get_header()[elem], end="") - -# Get preliminary data for print loop -residues = pdb_file.residues -bridges = get_bridges(residues) -ladders = get_ladders(bridges, residues) -sheets = get_sheets(ladders) -helix = build_helix_patterns(residues) -# iterating over residues - - -for i,res in enumerate(residues): - #res.get_turns(residues, turns) - kappa = res.get_bends(residues)[0] - bend_symbol = res.get_bends(residues)[1] - t_co = res.get_tco(residues) - alpha = res.get_chirality(residues)[0] - chirality = res.get_chirality(residues)[1] - phi = res.get_phi_psi(residues)[0] - psi = res.get_phi_psi(residues)[1] - x_ca = res.atoms["CA"].coord_x - y_ca = res.atoms["CA"].coord_y - z_ca = res.atoms["CA"].coord_z - turns = build_turns_patterns(residues) - turn_3 = print_turn_pattern(residues, res, turns[0]) - turn_4 = print_turn_pattern(residues, res, turns[1]) - turn_5 = print_turn_pattern(residues, res, turns[2]) - helix_3 = print_helix_pattern(residues, res, helix[0]) - helix_4 = print_helix_pattern(residues, res, helix[1]) - helix_5 = print_helix_pattern(residues, res, helix[2]) - if i in bridges.keys(): - bp1 = bridges[i].j + if(args.verbose): + print_dssp(pdb_file, raw_ladd=True) else: - bp1 = 0 - print(i+1, res.resid, res.chain_id, res.res_letter, helix_3, helix_4, helix_5, turn_3, turn_4, turn_5, bend_symbol, chirality, bp1,round(t_co, 3), round(kappa,1), - round(alpha, 1), round(phi, 1), round(psi,1), round(x_ca,1), round(y_ca,1), round(z_ca,1)) - - -for ladd1 in ladders: - ladd_1 = ladders[ladd1] - for ladd2 in ladders: - ladd_2 = ladders[ladd2] - for bridge in ladd_1['bridges']: - if bridge.j in res_list(ladd_2): - print("ladder", ladd_1['i'],"-",ladd_1['j'], "|", bridge.i, "...", bridge.j, "| ladder", ladd_2['i'], ladd_2['j']) - -# for ind, sheet in sheets.items(): -# print(ind, sheet) -# for ladder in sheet: -# print(ladder['start'], ladder['end']) + print_dssp(pdb_file) + diff --git a/src/geom.py~ b/src/geom.py~ deleted file mode 100644 index 4e96251..0000000 --- a/src/geom.py~ +++ /dev/null @@ -1,66 +0,0 @@ - -def vector_from_pos(a, b): - xAB = b[0]-a[0] - yAB = b[1]-a[1] - zAB = b[2]-a[2] - coord_AB = [xAB,yAB,zAB] - return coord_AB - -def vector_norm(v): - norm = math.sqrt(v[0]**2 + v[1]**2 + v[2]**2) - return norm - -def dot_product(v1,v2): - dot_product = v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2] - return(dot_product) - -def position_vector(c): - vector = vector_from_pos([0,0,0],c) - return vector - -def vectors_substr(v1, v2): - return ([v1[0]-v2[0], - v1[1]-v2[1], - v1[2]-v2[2]]) - -def vector_angles(v1,v2): - dot_prod = dot_product(v1,v2) - norm_v1 = vector_norm(v1) - norm_v2 = vector_norm(v2) - term = dot_prod/(abs(norm_v1)*abs(norm_v2)) - rad_angle = math.acos(term) - return rad_angle - -def calc_dihedral(u1, u2, u3, u4): - """ Calculate dihedral angle method. From bioPython.PDB - (adapted to np.array) - Calculate the dihedral angle between 4 vectors - representing 4 connected points. The angle is in - [-pi, pi]. - - Adapted function of dihedral_angle_from_coordinates.py - by Eric Alcaide. - Source : https://gist.github.com/EricAlcaide - URL : https://gist.github.com/EricAlcaide/30d6bfdb8358d3a57d010c9a501fda56 - """ - #convert coords to numpy arrays - u1 = np.array(u1) - u2 = np.array(u2) - u3 = np.array(u3) - u4 = np.array(u4) - - a1 = u2 - u1 - a2 = u3 - u2 - a3 = u4 - u3 - - v1 = np.cross(a1, a2) - v1 = v1 / (v1 * v1).sum(-1)**0.5 - v2 = np.cross(a2, a3) - v2 = v2 / (v2 * v2).sum(-1)**0.5 - porm = np.sign((v1 * a3).sum(-1)) - rad = np.arccos((v1*v2).sum(-1) / ((v1**2).sum(-1) * (v2**2).sum(-1))**0.5) - if not porm == 0: - rad = rad * porm - - deg_angle = rad*(180/math.pi) - return deg_angle diff --git a/src/pdb.py b/src/pdb.py index cb1b216..ce616aa 100644 --- a/src/pdb.py +++ b/src/pdb.py @@ -1,5 +1,5 @@ # custom imports -from atom import * +from chem import * # import libs import sys import pymol @@ -67,9 +67,7 @@ class PDBFile: # hydrogens should represent in average 50% of total atoms... # We use 30% threshold... if(count_h/len(temp_atoms)<0.30): - #if(output_pdb==None): - print("Need to add hydrogens ! If you want the modified PDB file, " - "please use the -o output.pdb argument") + #Need to add hydrogens self.add_hydrogens(filename) def check_hydrogens(self, atoms):