Commit ba104221 authored by Florian RICHOUX's avatar Florian RICHOUX
Browse files

Merge branch 'release/tcbb'

parents 8f13a246 4121879c
This project aims at discovering unmet descriptions of Protein-Protein Interactions using a large-scale approach.
\ No newline at end of file
# What is DeepPPI?
DeepPPI is a deep learning project to predict protein-protein interactions based on protein sequence only.
It is declined into PyTorch models (deprecated) and Keras models.
# How to use our Keras models?
In the 'keras' folder, run the main.py with Python. For instance, to print helps, type
`python main.py -h`
To train a model, you must give both a training set and a validation set. For instance
`python main.py -train ../data/small_1166_train.txt -val ../data/small_1166_val.txt`
In our IEEE/ACM Transactions on Computational Biology and Bioinformatics / arXiv paper, our datasets can be found in:
..* text files in `data/mirror` for our regular dataset
..* text files in `data/mirror/double` for our strict dataset
The two models from these papers can be found at `keras/models/fc2_20_2dense.py` and `keras/models/lstm32_3conv3_2dense_shared.py`. Trained weights of these models are in `keras/weights/`.
To load weights and test a model on a test set, the command is the following:
`python main.py -load <weights_file> -test <test_file>`
\ No newline at end of file
'''
(c) 2010-2011 Thomas Holder, MPI for Developmental Biology
Python parser for AAindex: Amino Acid Index Database
http://www.genome.jp/aaindex/
PyMOL commands:
aaindex2b
pmf
'''
from __future__ import print_function
from __future__ import absolute_import
import sys
import os
_aaindex = dict()
_pymol_auto_arg_update = lambda: None
def search(pattern, searchtitle=True, casesensitive=False):
'''
Search for pattern in description and title (optional) of all records and
return matched records as list. By default search case insensitive.
'''
whatcase = lambda i: i
if not casesensitive:
pattern = pattern.lower()
whatcase = lambda i: i.lower()
matches = []
for record in _aaindex.values():
if pattern in whatcase(record.desc) or searchtitle and pattern in whatcase(record.title):
matches.append(record)
return matches
def grep(pattern):
'''
Search for pattern in title and description of all records (case
insensitive) and print results on standard output.
'''
for record in search(pattern):
print(record)
class Record:
'''
Amino acid index (AAindex) Record
'''
aakeys = 'ARNDCQEGHILKMFPSTWYV'
def __init__(self):
self.key = None
self.desc = ''
self.ref = ''
self.authors = ''
self.title = ''
self.journal = ''
self.correlated = dict()
self.index = dict()
self.comment = ''
def extend(self, row):
i = len(self.index)
for x in row:
self.index[self.aakeys[i]] = x
i += 1
def get(self, aai, aaj=None, d=None):
assert aaj is None
return self.index.get(aai, d)
def __getitem__(self, aai):
return self.get(aai)
def median(self):
x = sorted(filter(None, self.index.values()))
half = len(x) // 2
if len(x) % 2 == 1:
return x[half]
return (x[half - 1] + x[half]) / 2.0
def __str__(self):
desc = self.desc.replace('\n', ' ').strip()
return '%s(%s: %s)' % (self.__class__.__name__, self.key, desc)
class MatrixRecord(Record):
'''
Matrix record for mutation matrices or pair-wise contact potentials
'''
def __init__(self):
Record.__init__(self)
self.index = []
self.rows = dict()
self.cols = dict()
def extend(self, row):
self.index.append(row)
def _get(self, aai, aaj):
i = self.rows[aai]
j = self.cols[aaj]
return self.index[i][j]
def get(self, aai, aaj, d=None):
try:
return self._get(aai, aaj)
except:
pass
try:
return self._get(aaj, aai)
except:
return d
def __getitem__(self, aaij):
return self.get(aaij[0], aaij[1])
def median(self):
x = []
for y in self.index:
x.extend(filter(None, y))
x.sort()
if len(x) % 2 == 1:
return x[len(x) // 2]
return sum(x[len(x) // 2 - 1:len(x) // 2 + 1]) / 2.0
def get(key):
'''
Get record for key
'''
if len(_aaindex) == 0:
init()
return _aaindex[key]
def _float_or_None(x):
if x == 'NA' or x == '-':
return None
return float(x)
def init(path=None, index='13'):
'''
Read in the aaindex files. You need to run this (once) before you can
access any records. If the files are not within the current directory,
you need to specify the correct directory path. By default all three
aaindex files are read in.
'''
index = str(index)
if path is None:
for path in [os.path.split(__file__)[0], '.', cmd.get('fetch_path')]:
if os.path.exists(os.path.join(path, 'aaindex' + index[0])):
break
print('path =', path, file=sys.stderr)
if '1' in index:
_parse(path + '/aaindex1', Record)
if '2' in index:
_parse(path + '/aaindex2', MatrixRecord)
if '3' in index:
_parse(path + '/aaindex3', MatrixRecord)
_pymol_auto_arg_update()
def init_from_file(filename, type=Record):
_parse(filename, type)
def _parse(filename, rec, quiet=True):
'''
Parse aaindex input file. `rec` must be `Record` for aaindex1 and
`MarixRecord` for aaindex2 and aaindex3.
'''
if not os.path.exists(filename):
if sys.version_info[0] < 3:
from urllib import urlretrieve
else:
from urllib.request import urlretrieve
url = 'ftp://ftp.genome.jp/pub/db/community/aaindex/' + os.path.split(filename)[1]
print('Downloading "%s"' % (url))
filename = urlretrieve(url, filename)[0]
print('Saved to "%s"' % (filename))
f = open(filename)
current = rec()
lastkey = None
for line in f:
key = line[0:2]
if key[0] == ' ':
key = lastkey
if key == '//':
_aaindex[current.key] = current
current = rec()
elif key == 'H ':
current.key = line[2:].strip()
elif key == 'R ':
current.ref += line[2:]
elif key == 'D ':
current.desc += line[2:]
elif key == 'A ':
current.authors += line[2:]
elif key == 'T ':
current.title += line[2:]
elif key == 'J ':
current.journal += line[2:]
elif key == '* ':
current.comment += line[2:]
elif key == 'C ':
a = line[2:].split()
for i in range(0, len(a), 2):
current.correlated[a[i]] = float(a[i + 1])
elif key == 'I ':
a = line[1:].split()
if a[0] != 'A/L':
current.extend(map(_float_or_None, a))
elif list(Record.aakeys) != [i[0] for i in a] + [i[-1] for i in a]:
print('Warning: wrong amino acid sequence for', current.key)
else:
try:
assert list(Record.aakeys[:10]) == [i[0] for i in a]
assert list(Record.aakeys[10:]) == [i[2] for i in a]
except:
print('Warning: wrong amino acid sequence for', current.key)
elif key == 'M ':
a = line[2:].split()
if a[0] == 'rows':
if a[4] == 'rows':
a.pop(4)
assert a[3] == 'cols' and len(a) == 6
i = 0
for aa in a[2]:
current.rows[aa] = i
i += 1
i = 0
for aa in a[5]:
current.cols[aa] = i
i += 1
else:
current.extend(map(_float_or_None, a))
elif not quiet:
print('Warning: line starts with "%s"' % (key))
lastkey = key
########## PYMOL ###########
# from Bio.SCOP.Raf import to_one_letter_code
# See also http://www.pymolwiki.org/index.php/Aa_codes
to_one_letter_code = {'PAQ': 'Y', 'AGM': 'R', 'ILE': 'I', 'PR3': 'C',
'GLN': 'Q', 'DVA': 'V', 'CCS': 'C', 'ACL': 'R', 'GLX': 'Z', 'GLY': 'G',
'GLZ': 'G', 'DTH': 'T', 'OAS': 'S', 'C6C': 'C', 'NEM': 'H', 'DLY': 'K',
'MIS': 'S', 'SMC': 'C', 'GLU': 'E', 'NEP': 'H', 'BCS': 'C', 'ASQ': 'D',
'ASP': 'D', 'SCY': 'C', 'SER': 'S', 'LYS': 'K', 'SAC': 'S', 'PRO': 'P',
'ASX': 'B', 'DGN': 'Q', 'DGL': 'E', 'MHS': 'H', 'ASB': 'D', 'ASA': 'D',
'NLE': 'L', 'DCY': 'C', 'ASK': 'D', 'GGL': 'E', 'STY': 'Y', 'SEL': 'S',
'CGU': 'E', 'ASN': 'N', 'ASL': 'D', 'LTR': 'W', 'DAR': 'R', 'VAL': 'V',
'CHG': 'A', 'TPO': 'T', 'CLE': 'L', 'GMA': 'E', 'HAC': 'A', 'AYA': 'A',
'THR': 'T', 'TIH': 'A', 'SVA': 'S', 'MVA': 'V', 'SAR': 'G', 'LYZ': 'K',
'BNN': 'A', '5HP': 'E', 'IIL': 'I', 'SHR': 'K', 'HAR': 'R', 'FME': 'M',
'PYX': 'C', 'ALO': 'T', 'PHI': 'F', 'ALM': 'A', 'PHL': 'F', 'MEN': 'N',
'TPQ': 'A', 'GSC': 'G', 'PHE': 'F', 'ALA': 'A', 'MAA': 'A', 'MET': 'M',
'UNK': 'X', 'LEU': 'L', 'ALY': 'K', 'SET': 'S', 'GL3': 'G', 'TRG': 'K',
'CXM': 'M', 'TYR': 'Y', 'SCS': 'C', 'DIL': 'I', 'TYQ': 'Y', '3AH': 'H',
'DPR': 'P', 'PRR': 'A', 'CME': 'C', 'IYR': 'Y', 'CY1': 'C', 'TYY': 'Y',
'HYP': 'P', 'DTY': 'Y', '2AS': 'D', 'DTR': 'W', 'FLA': 'A', 'DPN': 'F',
'DIV': 'V', 'PCA': 'E', 'MSE': 'M', 'MSA': 'G', 'AIB': 'A', 'CYS': 'C',
'NLP': 'L', 'CYQ': 'C', 'HIS': 'H', 'DLE': 'L', 'CEA': 'C', 'DAL': 'A',
'LLP': 'K', 'DAH': 'F', 'HMR': 'R', 'TRO': 'W', 'HIC': 'H', 'CYG': 'C',
'BMT': 'T', 'DAS': 'D', 'TYB': 'Y', 'BUC': 'C', 'PEC': 'C', 'BUG': 'L',
'CYM': 'C', 'NLN': 'L', 'CY3': 'C', 'HIP': 'H', 'CSO': 'C', 'TPL': 'W',
'LYM': 'K', 'DHI': 'H', 'MLE': 'L', 'CSD': 'A', 'HPQ': 'F', 'MPQ': 'G',
'LLY': 'K', 'DHA': 'A', 'DSN': 'S', 'SOC': 'C', 'CSX': 'C', 'OMT': 'M',
'DSP': 'D', 'PTR': 'Y', 'TRP': 'W', 'CSW': 'C', 'EFC': 'C', 'CSP': 'C',
'CSS': 'C', 'SCH': 'C', 'OCS': 'C', 'NMC': 'G', 'SEP': 'S', 'BHD': 'D',
'KCX': 'K', 'SHC': 'C', 'C5C': 'C', 'HTR': 'W', 'ARG': 'R', 'TYS': 'Y',
'ARM': 'R', 'DNP': 'A'}
def aaindex2b(key='KYTJ820101', selection='(all)', quiet=0, var='b'):
'''
DESCRIPTION
"aaindex" looks up the Amino Acid Index from
http://www.genome.jp/aaindex/
for the given key and assignes b-factors to the given selection. Unknown
residues get the average index value assigned.
USAGE
aaindex2b [key [, selection]]
ARGUMENTS
key = string: Key of AAindex entry
selection = string: atoms to assign b-factors {default: (all)}
EXAMPLE
# Hydropathy index by Kyte-Doolittle
aaindex2b KYTJ820101
spectrumany b, white yellow forest
show surface
'''
from pymol import cmd, stored
entry = get(key)
median = entry.median()
if int(quiet) != 0:
print(entry.desc.strip())
def lookup(resn):
one_letter = to_one_letter_code.get(resn, 'X')
value = entry.get(one_letter)
if value is None:
return median
return value
stored.aaindex = lookup
cmd.alter(selection, var + '=stored.aaindex(resn)')
def pmf(key, cutoff=7.0, selection1='(name CB)', selection2='', state=1, quiet=1):
'''
DESCRIPTION
Potential of Mean Force
ARGUMENTS
key = string: aaindex key
cutoff = float: distance cutoff {default: 7.0}
cutoff = (float, float): distance shell
selection1 = string: atom selection {default: (name CB)}
selection2 = string: atom selection {default: selection1}
NOTES
Does also support a list of keys and a list of cutoffs to deal with
multiple distance shells.
EXAMPLES
# distance dependent c-beta contact potentials
pmf SIMK990101, 5, /2x19//A//CB
pmf SIMK990102, [5, 7.5], /2x19//A//CB
pmf [SIMK990101, SIMK990102, SIMK990103], [0, 5, 7.5, 10], /2x19//A//CB
# interface potential
sidechaincenters 2x19_scc, 2x19
pmf KESO980102, 7.0, /2x19_scc//A, /2x19_scc//B
distance /2x19_scc//A, /2x19_scc//B, cutoff=7.0
'''
from pymol import cmd, stored
from chempy import cpv
if cmd.is_string(key):
if key.lstrip().startswith('['):
key = cmd.safe_alpha_list_eval(key)
else:
key = [key]
if cmd.is_string(cutoff):
cutoff = eval(cutoff)
if not cmd.is_sequence(cutoff):
cutoff = [cutoff]
if len(cutoff) == len(key):
cutoff = [0.0] + list(cutoff)
if len(cutoff) != len(key) + 1:
print('Error: Number of keys and number of cutoffs inconsistent')
return
state = int(state)
quiet = int(quiet)
if len(selection2) == 0:
selection2 = selection1
if not quiet and len(key) > 1:
print('Distance shells:')
for i in range(len(key)):
print('%s %.1f-%.1f' % (key[i], cutoff[i], cutoff[i + 1]))
idmap = dict()
cmd.iterate_state(state, '(%s) or (%s)' % (selection1, selection2),
'idmap[model,index] = [(resn,name),(x,y,z)]', space={'idmap': idmap})
twoN = cmd.count_atoms(selection1) + cmd.count_atoms(selection2)
pairs = cmd.find_pairs(selection1, selection2, cutoff=max(cutoff),
state1=state, state2=state)
if len(pairs) == 0:
print('Empty pair list')
return 0.0
matrix = list(map(get, key))
for i in matrix:
assert isinstance(i, MatrixRecord)
i_list = list(range(len(key)))
u_sum = 0
count = 0
for id1, id2 in pairs:
a1 = idmap[id1]
a2 = idmap[id2]
r = cpv.distance(a1[1], a2[1])
for i in i_list:
if cutoff[i] <= r and r < cutoff[i + 1]:
try:
aa1 = to_one_letter_code[a1[0][0]]
aa2 = to_one_letter_code[a2[0][0]]
u_sum += matrix[i].get(aa1, aa2)
count += 1
except:
print('Failed for', a1[0], a2[0])
value = float(u_sum) / twoN
if not quiet:
print('PMF: %.4f (%d contacts, %d residues)' % (value, count, twoN))
return value
try:
from pymol import cmd
cmd.extend('aaindex2b', aaindex2b)
cmd.extend('pmf', pmf)
def pymol_auto_arg_update():
aaindexkey_sc = cmd.Shortcut(_aaindex.keys())
cmd.auto_arg[0].update({
'aaindex2b': [aaindexkey_sc, 'aaindexkey', ', '],
'pmf': [aaindexkey_sc, 'aaindexkey', ', '],
})
cmd.auto_arg[1].update({
'aaindex2b': [cmd.selection_sc, 'selection', ''],
})
cmd.auto_arg[2].update({
'pmf': [cmd.selection_sc, 'selection', ''],
})
cmd.auto_arg[3].update({
'pmf': [cmd.selection_sc, 'selection', ''],
})
_pymol_auto_arg_update = pymol_auto_arg_update
except:
pass
# vi: ts=4:sw=4:smarttab:expandtab
This diff is collapsed.
This source diff could not be displayed because it is too large. You can view the blob instead.
This diff is collapsed.
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import os, sys
#sys.path.append('/comptes/teletchea-s/venv_specifiques/PPIpredict/lib/python3.4/site-packages/')
#sys.path.append('/home/flo/share/Developpement/PPIpredict/lib/python3.4/')
# Utilisation de Bio.PDB
# http://biopython.org/wiki/The_Biopython_Structural_Bioinformatics_FAQ
from Bio.PDB import *
### Specify path for aaindex init othwerwise python cmd is called...
import aaindex as aaindex
aaindex.init(os.getcwd(),index='13')
# TODO FR : normaliser les données pour les coordonnées 3D
# TODO FR : piste 1 -> normaliser en fonction de la taille de la maille cristallo : difficile car on doit tenir compte du nmbre de structures dans la maille
# TODO FR : piste 2 -> normalisation magique avec numpy -> http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.norm.html
# TODO ST : retourner un vecteur par acide aminé (entier de 0 à 19)
# TODO ST : retourner un vecteur par type d'atome (0 à 4 (ou 5 ?))
# TODO ST : vecteur interaction séparé du vecteur data
# TODO ST : meilleur méthode de séparation des acides aminés à l'interface
# voir l'option NeighborSearch ou http://www.pymolwiki.org/index.php/InterfaceResidues
# TODO ST : mieux séparer les chaines pour être sûr d'avoir deux protéines différentes, en interaction 1:1
one_letter ={'VAL':'V', 'ILE':'I', 'LEU':'L', 'GLU':'E', 'GLN':'Q',
'ASP':'D', 'ASN':'N', 'HIS':'H', 'TRP':'W', 'PHE':'F', 'TYR':'Y',
'ARG':'R', 'LYS':'K', 'SER':'S', 'THR':'T', 'MET':'M', 'ALA':'A',
'GLY':'G', 'PRO':'P', 'CYS':'C'}
def get_residues_in_contact(model, aa):
ca1=aa['CA']
distance=11
for c in model:
if c != aa.get_parent():
for residue in c:
ca2=residue['CA']
ca2name=residue.get_resname()
distance = ca2-ca1
if distance < 10:
# print ("La distance entre {} et {} est : {}".format(c, ca2name, distance))
return ca2name, 1
else:
return ca2name, 0
def get_amino_acid_parameters(AA):
aa=aaindex.to_one_letter_code[AA]
# TODO : sélection arbitraire des sources
# bien entendu, à raffiner ...
# http://www.genome.jp/aaindex/AAindex/list_of_indices
# http://www.genome.jp/aaindex/AAindex/Appendix
parameters={}
parameters['JOND750101']='Hydrophobicity (Jones, 1975)'
parameters['KRIW790103']='Side chain volume (Krigbaum-Komoriya, 1979)'
parameters['CHOP780201']='Normalized frequency of alpha-helix (Chou-Fasman, 1978b)'
parameters['CHOP780202']='Normalized frequency of beta-sheet (Chou-Fasman, 1978b)'
parameters['DAYM780101']='Amino acid composition (Dayhoff et al., 1978b)'
parameters['CHAM820101']='Polarizability parameter (Charton-Charton, 1982)'
v=[]
for p in parameters.keys():
param=aaindex.get(p)
try:
v.append(param.get(aa))
except Exception as e:
v.append(0.0)
# print ("Le vecteur trouvé est {}".format(l))
return v
def output_old_version(structure):
vect_coords=[]
for model in structure:
for chain in model:
for residue in chain:
r = residue.get_resname()
interaction=0
s=''
# if r != 'HOH' and r != ' MG':
if r in one_letter.keys():
try:
s,interaction=get_residues_in_contact(model, residue)
# Phase 2 : uniquement AA1 AA2 contact ou pas
# vect_tout=[r,s,interaction]
# output_vecteurs(vect_tout)
except Exception as e:
# print ("L'erreur est %s" % e)
pass
vect_aa=get_amino_acid_parameters(r)
for atom in residue:
print r , atom.get_name(), atom.get_coord()
coords=atom.get_coord()
# Ne garder que le type d'atome ou son nom ?
# print ("%d %s %s %5.2f %5.2f %5.2f" % (interaction, r, atom.get_name()[0], coords[0],coords[1], coords[2]))
# Phase 3
# vect_coords=[interaction, r, atom.get_name()[0], coords[0],coords[1], coords[2]]
# vect_tout=vect_coords+vect_aa
# output_vecteurs(vect_tout)
def output_vecteurs(mon_vecteur):
print (mon_vecteur)
def main():
# p = MMCIFParser()
# structure = p.get_structure('demo', '1a2y.cif')
p=PDBParser()
structure = p.get_structure('demo', '1wq1.pdb')
# Structures à considérer pour le jeu de données initial
# - 1tv6
# - 1e27
# - 1pys
# - 1tnu (mais A,B)
# - 2cc1 (mais A,B)
# - 1m2t