Словарь ключ из файла pdb

Question

Я пытаюсь просмотреть файл.pdb, вычислить расстояние между альфа-атомами углерода из разных остатков в цепочках A и B белкового комплекса, а затем сохранить расстояние в словаре вместе с идентификатором цепи и номером остатка.

Например, если первый альфа-углерод ("CA") найден в остатке 100 в цепи A, а тот, к которому он привязан, находится в остатке 123 в цепи B I, то мой словарь должен выглядеть примерно так: d={(A, 100):[B, 123, distance_between_atoms]}

from Bio.PDB.PDBParser import PDBParser
parser=PDBParser()
struct = parser.get_structure("1trk", "1trk.pdb")

def getAlphaCarbons(chain):
    vec = []
    for residue in chain:
        for atom in residue:
            if atom.get_name() == "CA":
                vec = vec + [atom.get_vector()]
    return vec

def dist(a,b):
    return (a-b).norm()


chainA = struct[0]['A']
chainB = struct[0]['B']

vecA = getAlphaCarbons(chainA)
vecB = getAlphaCarbons(chainB)

t={}
model=struct[0]

for model in struct:
    for chain in model:
        for residue in chain:
            for a in vecA:
               for b in vecB:
                if dist(a,b)<=8:
                    t={(chain,residue):[(a, b, dist(a, b))]}

     break
print t  

Он запускает программу целую вечность, и мне пришлось прервать запуск (я сделал где-то бесконечный цикл??)

Я также пытался сделать это:

t = {i:[((a, b, dist(a,b)) for a in vecA) for b in vecB if dist(a, b) <= 8] for i in chainA}
print t

Но это печатает информацию об остатках в следующем формате:

<Residue PHE het=  resseq=591 icode= >: []    

Он не печатает ничего, что связано с расстоянием.

Большое спасибо, надеюсь, все понятно.

Answer 1

Настоятельно рекомендую использовать библиотеки C при расчете расстояний. Я использую mdtraj для такого рода вещей, и он работает намного быстрее, чем все циклы for в BioPython.

Чтобы получить все пары альфа-углеродов:

import mdtraj as md
def get_CA_pairs(self,pdbfile):
  traj = md.load_pdb(pdbfile)
  topology = traj.topology 
  CA_index = ([atom.index for atom in topology.atoms if (atom.name == 'CA')])
  pairs=list(itertools.combinations(CA_index,2))
return pairs

Тогда для быстрого расчета расстояний:

  def get_distances(self,pdbfile,pairs):
  #returns list of resid1, resid2,distances between CA-CA
   traj = md.load_pdb(pdbfile)
   pairs=self.get_CA_pairs(pdbfile)
   dist=md.compute_distances(traj,pairs)
#make dictionary you desire. 
   dict=dict(zip(CA, pairs))
  return dict

Это включает в себя все альфа-углерода. В mdtraj также должен быть идентификатор цепочки, чтобы выбрать CA из каждой цепочки.