def ConvertSequence2StructuralAlignment( src1, src2, source=None, format="plain", check_residues = 1):

    ca1 = GetPdbCoordinates( src1, renumber = 1)

    if len(ca1) == 0:
        raise "no coordinates found for %s" % src1

    ca2 = GetPdbCoordinates( src2, renumber = 1 )

    if len(ca2) == 0:
        raise "no coordinates found for %s" % src2

    if string.lower(format) not in ("plain",):
        raise "unknown alignment format %s" % format

    if source:
        lines = open(source, "r").readlines()
    else:
        lines = sys.stdin.readlines()

    ## replace gap characters
    lines = map(lambda x: re.sub( "\s", "", string.replace(x, ".", "-")), lines)
    if not lines:
        raise ValueError, "alignment is empty"

    lali = len(lines[0])

    current1 = 0
    current2 = 0

    index1 = 0
    index2 = 0

    output = []

    alignment = []

    for x in range(0, lali):

        res1 = lines[0][x]
        res2 = lines[1][x]

        if res1 != "-": current1+=1
        if res2 != "-": current2+=1

        try:
            while (ca1[index1][0] < current1): index1 += 1
            while (ca2[index2][0] < current2): index2 += 1                    
        except IndexError:
            break

        if res1 == "-" or res2 == "-":
            continue

        (i1, aa1, x1, y1, z1) = ca1[index1]
        (i2, aa2, x2, y2, z2) = ca2[index2]        

        if check_residues:
            if aa1 != res1:
                sys.stderr.write("# mismatch in 1:%s at residue alignment %i(%s) -> structure %i(%s)\n" %\
                                 (source, current1, res1, index1, aa1))
            if aa2 != res2:
                sys.stderr.write("# mismatch in 2:%s at residue %i(%s) -> %i(%s)\n" %\
                                 (source, current2, res2, index2, aa2))

        alignment.append( (x1, y1, z1, x2, y2, z2, 1) )

    return alignment
calculate a structural alignment from two pdb files.
dataset/ETHPy150Open CGATOxford/cgat/legacy/PdbTools.py/ConvertSequence2StructuralAlignment