Super

Some code examples.

A protein backbone superimposer.

The following will superimpose the proteins contained in the pdb files passed as arguments to the program.

For a Quick&Dirty(tm), throw-away program, most errors can just be ignored, so we have left only the error-check that ensures the superposition itself was successful. It is possible that the best rotation matrix that superimposes the structures turns out to be a reflexion. This is the case when the structures are enantiomers of each other. In this case the function Super will return and error and no structure.

If the function Super is called with one or two empty atom lists it will attempt to use the whole corresponding molecules and return an error if the atom numbers of both structures don't match.

Note that the program can be made to handle non-protein molecules just by changing the criteria to include indexes in superlist, and perhaps using XyzRead/XyzWrite instead of their PDB counterparts. The Super function is, thus, general and flexible. Note that it doesnt perform any protein sequence alignment, and will fail if the two proteins have a different number of residues. The addition of a specific protein superimposing function that does sequence alignment is expected to happen at some point.

Note that this program also uses the scu (silly chemistry utillities) library, that can be installed with: go get github.com/rmera/scu

    package main
    
    import (
    	"fmt"
        "github.com/rmera/scu"
    	"github.com/rmera/gochem"
    	"os"
    	"strings"
    )
    
//This will Superimpose the backbone of two proteins.
func main() {
	backbone := []string{"C", "CA", "O", "N"} //The PDB name of the atoms in the backbone.
	mol1, _ := chem.PdbRead(os.Args[1], true) //true means that we try to read the symbol from the PDB file.
	mol2, _ := chem.PdbRead(os.Args[2], true)
	mols := []*chem.Molecule{mol1, mol2}
	superlist := make([][]int, 2, 2)
	//We collect the atoms that are part of the backbone.
	for molnumber, mol := range mols {
		for atomindex, atom := range mol.Atoms {
			if scu.IsIn(atom.Name, backbone) {
				superlist[molnumber] = append(superlist[molnumber], atomindex)
			}
		}
	}
	mol1.Coords[0], err = chem.Super(mol1.Coords[0], mol2.Coords[0], superlist[0], superlist[1])
	if err != nil {
		panic(err.Error())
	}
	newname := strings.Replace(os.Args[1], ".pdb", "_super.pdb", -1)
	chem.PdbWrite(mol1, newname)
}

A more proper version of the previous program, with error checking, is slightly more complex, but still very readable and concise:

    package main
    
    import (
    	"fmt"
        "github.com/rmera/scu"
    	"github.com/rmera/gochem"
    	"os"
    	"strings"
    )
    
//This will Superimpose the backbone of two proteins.
func main() {
	backbone := []string{"C", "CA", "O", "N"} //The PDB name of the atoms in the backbone.
	if len(os.Args) < 3 || !strings.Contains(os.Args[1], ".pdb") || !strings.Contains(os.Args[2], ".pdb") {
		fmt.Printf("Usage: %s  mol1.pdb mol2.pdb\nNo caps in the pdb extension allowed.\n", os.Args[0])
		os.Exit(1)
	}
	mol1, err := chem.PdbRead(os.Args[1], true) //true means that we try to read the symbol from the PDB file.
	mol2, err2 := chem.PdbRead(os.Args[2], true)
	if err != nil || err2 != nil {
		panic("Unable to open input files!")
	}
	mols := []*chem.Molecule{mol1, mol2}
	superlist := make([][]int, 2, 2)
	//We collect the atoms that are part of the backbone.
	for molnumber, mol := range mols {
		for atomindex, atom := range mol.Atoms {
			if scu.IsIn(atom.Name, backbone) {
				superlist[molnumber] = append(superlist[molnumber], atomindex)
			}
		}
	}
	mol1.Coords[0], err = chem.Super(mol1.Coords[0], mol2.Coords[0], superlist[0], superlist[1])
	if err != nil {
		panic(err.Error())
	}
	newname := strings.Replace(os.Args[1], ".pdb", "_super.pdb", -1)
	chem.PdbWrite(mol1, newname)
}