On May 21, 2014, at 3:54 AM, Michał Kadlof <ramidas@gmail.com> wrote:
Hello,
I would like to ask if there is any quick way to make aligment (with secondary structure score included) beetwen two homologous protein structures, and save them in Stockholm format (information about missing residues is needed). I have houndreds of such pairs, and it would be very helpfull to do it by script for headless Chimera.
Hi Michał, Well, depending on your definition of "quick", the answer is a qualified yes. The first thing would be to look through this primer for looping through data files and doing things with them in Chimera: http://www.cgl.ucsf.edu/chimera/docs/ProgrammersGuide/basicPrimer.html The second thing to know is that you can superimpose two structures (with secondary structures considered) with the "matchmaker" (or "mm") command. Type "help matchmaker" in Chimera's command line to see the help page with the command's options. The tricky part is running Match->Align to generate an alignment based on the superposition without using the graphical interface, and then saving the results in Stockholm format. You will need to use the Python interface to Match->Align, basically: from StructSeqAlign import makeAlignment seqs = makeAlignment(chains, makeMAV=False) class fakeMAV: fileAttrs = {} from MultAlignViewer.formatters import saveStockholm f = open("output.sto", "w") saveStockholm.save(f, fakeMAV(), seqs, {}) f.close() Note that makeAlignment() has additional optional arguments for changing the alignment parameters from their defaults. If you look at the code I think they're reasonably self-explanatory (e.g. circular=True/False). So the main question is what is the 'chains' argument to the makeAlignment() call? Well, it's a list of the chains you want in the alignment. Let's say you want to match chain A from model 0 to chain B of model 1. You'd need code like this (which assumes you have just those two models open): from chimera import openModels, Molecule models = openModels.list(modelTypes=[Molecule]) if models[0].id == 0: m0, m1 = models else: m1, m0 = models chains = [m0.sequence('A'), m1.sequence('B')] --Eric Eric Pettersen UCSF Computer Graphics Lab http://www.cgl.ucsf.edu