User:Remig/plico/remapNT

From Jmol
< User:Remig‎ | plico
Revision as of 17:22, 12 April 2016 by Remig (talk | contribs)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to navigation Jump to search

RemapNT allows you to change the chain ID, atom numbers and/or residue numbers of a polynucleotide chain by mouse actions. It also calculates group values [nucleotide names (DU, A, etc.)]. Finally it prints the resultant 1 char string to the console.

When you click on a polynucleotide chain, it gives the current chain ID, residue, residue number and atom number of the most 5'ward atom in that chain. You may then edit any value (except residue). Remap then remaps the entire chain based on those values by conventional increments and identifies each nucleotide residue.

Note that it will also remove all waters, ligands, hydrogens and %B alternates when any chain is updated.

RemapNT is a member of the Plico suite of protein folding tools described here. It may be installed and accessed as a macro with the file:

#   remapNT - Jmol script by Ron Mignery
#   v1.7  beta     4/12/2016 -require latest common includes
#
#   Calculate or change polynucleotide chain, atom number, residue names
#    and/or residue numbers and print the resultant 1 char string
#
gBusy = false
kRemapNT=2

function find_5_prime(pIdx) {
    while (pIdx >= 0) {

        # Find C3'
        var c3pIdx = -1
        select {atomIndex=pIdx}
        var cSet = {selected}
        for (var i = 1; i <= cSet.size; i++) {
            var ocSet = connected(cSet[i])
            for (var j = 1; j <= ocSet.size; j++) {
                var occSet = connected(ocSet[j]) and not {atomIndex=pIdx}
                for (var k = 1; k <= occSet.size; k++) {
                    if (connected(occSet[k]).size > 2) {
                        c3pIdx = 0 + occSet[k].atomIndex
                    }
                }
            }
        }

        if (c3pIdx < 0) {
            return pIdx
        }
        # Find C4'
        var endIdx = -1
        cSet = connected({atomIndex=c3pIdx}) and not {oxygen}
        for (var i = 1; i <= cSet.size; i++) {
            var ocSet = connected(cSet[i]) and {oxygen}
            for (var j = 1; j <= ocSet.size; j++) {
                if (connected(ocSet[j]).size > 1) {
                    if (endIdx = -1) {
                        endIdx = 0 + cSet[i].atomIndex
                    }
                    else {
                        var Ox = connected(ocSet[j]) and not {atomindex=@{cset[i].atomIndex}}
                        if (connected(Ox) > 1) {
                            endIdx = 0 + cSet[i].atomIndex
                        }
                    
                    }
                }
            }
        }

        # Find C5'
        var c5idx = -1
        cSet = (connected({atomIndex=endIdx}) and not {atomIndex=c3pIdx}
            and not {oxygen})
        if (cSet) {
            c5idx = 0 + cSet[1].atomIndex
        }

        # Find O5'
        var o5idx = -1
        cSet = connected({atomIndex=c5idx}) and {oxygen}
        if (cSet) {
            o5idx = 0 + cSet[1].atomIndex
        }

        if (o5idx < 0) {
            return c5idx
        }

        # Find P
        pIdx = -1
        cSet = connected({atomIndex=o5idx}) and {phosphorus}
        if (cSet) {
            Pidx =  0 + cSet[1].atomIndex
        }

        if (Pidx < 0) {
            return o5idx
        }
    }

    return -1
}

function find_p_idx(idx) {
    select {atomIndex=idx}
    var cSet = {selected}
    while (cSet) {
        for (var i = 1; i <= cSet.size; i++) {
            if (cSet[i].element == "P") {
                return cSet[i].atomIndex
            }
        }
        cSet = connected({selected}) and not {selected} and not {hydrogen}
        select {selected} or cset
    }
    return -1
}

# cSet, s, and newGreek are arrays and thus passed by reference
function ring_common(cSet, nIdx, s, newGreek, nextGreek) {
    if (cSet.size > 2) {
        print format("Unrecognized structure with set %s", cSet)
    }
    var oldGreek = 0 + newGreek[1]
    newGreek[1] = nextGreek
    for (var i = 1; i <= cSet.size; i++) {
        var ccSet = connected(cSet[i])
        if (ccSet.size == 1) {
            if (cSet[i].element == ccSet[1].element) {
                s[2] = i
                s[1] = ((i > 1) ? 1 : 2)
                newGreek[i] = 1 + nextGreek
                newGreek[s[1]] = nextGreek
                return cSet[s[1]].atomIndex
            }
            else {
                s[1] = i
                s[2] = ((i > 1) ? 1 : 2)
                newGreek[i] = oldGreek
                newGreek[s[2]] = nextGreek
                return cSet[s[2]].atomIndex
            }
        }
    }
    return cSet[1].atomIndex
}

# Bound to ALT-LEFT-CLICK by plico_remap_nt
function remap_nt_cargo_mb() {
    var idx =_atomPicked
    if ({atomIndex=idx}.element == "H") {
        idx = connected({atomIndex=idx})[1].atomIndex
    }
    delete {hydrogen}
    delete {hoh}
    delete %B
    delete ligands
    connect

    remap_nt( idx, false, 0)
    
    set echo TOP LEFT
    echo @gEcho
    background ECHO yellow
    refresh
    print_1c_chain( newChain)
}
        
function remap_nt(idx, auto, base) {

    # If P can be found
    var pIdx = find_p_idx(idx)
    var isValid = false
    var newResno = 1
    var newChain = "A"
    var newAtomno = 1
    var t5idx = -1
    if (pIdx >= 0) {

        t5idx = find_5_prime(pIdx)
        if (t5idx >= 0) {
            var f = {atomIndex=t5idx}.file
            var m = {atomIndex=t5idx}.model
            newResno = {atomIndex=t5idx}.resno
            newChain = {atomIndex=t5idx}.chain
            newAtomno = {atomIndex=t5idx}.atomno
            select {thisModel}
            color {selected} @gScheme
            select {(chain=newChain) and thisModel}
            color {selected} @gAltScheme
            refresh

            if (auto) {
                newResno = base
                isValid = true
                newChain = gChain1
            }
            else {
            
                # Prompt for new designators
                var p = prompt(("Enter 5\'-terminal atom designator as\n"
                    + "   <resno>:<chain>#<atomno>"),
                    format("%d:%s#%d", newResno, newChain, newAtomno))%0
                # If valid
                if (p != "null") {
                    var iColon = p.find(":")
                    if (iColon > 0) {
                        var iHash = p.find("#")
                        if (iHash > 0) {
                            newResno = 0 + (p[1][iColon-1])
                            newChain = p[iColon+1][iHash-1]
                            newAtomno = 0 + (p[iHash+1][0])
                            if ((newResno > 0)
                                and (newChain.size == 1)
                                and (newAtomno > 0)) {
                                    isValid = true
                            }
                        }
                    }
                    if (not isValid) {
                        prompt ("Entry not valid!")
                    }
                }
            }
        }
    }

    if (isValid) {
        background ECHO pink
        refresh

        # Build inline pdb file
        var ls = "data \"append remapNT\"\n"    # global PDB atom record
        var rs = ""

        select {atomIndex=t5idx}
        var cSet = {selected}
        var nextAtomName = {atomIndex=t5idx}.element
        var newGroup = "UNK"
        var newGreek = array("", "", "", "")
        var nIdx = t5idx
        var c1pIdx = -1
        var o3pIdx = -1
        var stopIdx = -1
        var endIdx = -1
        var isRNA = false
        var first = true
        var psu = false
        while (cSet) {
            var s = array(1, 2, 3, 4)
            var iKeep = -1
            var iDrop = -1
            switch( nextAtomName) {
            case "O" :
                newGreek[1] = (first ? "5\'" : "P3")
                nextAtomName = (first ? "C5\'" : "P")
                nIdx = cSet[1].atomIndex
                break
            case "P" :
                newGreek[1] = ""
                nextAtomName = "OP"
                nIdx = cSet[1].atomIndex
                break
            case "OP" :
                var oc5set = ({})
                for (var i = 1; i <= cSet.size; i++) {
                    newGreek[i] = ""
                    if (connected(cSet[i]).size > 1) {
                        s[cSet.size] = i
                        newGreek[cSet.size] = "5\'"
                        oc5set = connected(cSet[i]) and {carbon}
                        nIdx = cSet[i].atomIndex
                    }
                }
                var isP1 = false
                if (cSet.size > 3) {
                    newGreek[1] = "P1"
                    newGreek[2] = "P2"
                    newGreek[3] = "P3"
                    s[1] = 1
                    s[2] = 2
                    s[3] = 3
                }                
                else {
                    for (var i = 1; i <= cSet.size; i++) {
                        if (i != s[cSet.size]) {
                            if ((isP1 == false)
                                and (abs(angle(cSet[i], {atomIndex=pIdx},
                                cSet[s[cSet.size]], oc5set[1])) < 90.0)) {
                                s[1] = i
                                newGreek[1] = "P1"
                                isP1 = true
                            }
                            else {
                                s[2] = i
                                newGreek[2] = "P2"
                            }
                        }
                    }
                }
                #nIdx = pIdx
                nextAtomName = "C5\'"
                break
            case "C5\'" :
                nIdx = cSet[1].atomIndex
                newGreek[1] = "5\'"
                nextAtomName = "C4\'"
                break
            case "C4\'" :
                nIdx = cSet[1].atomIndex
                newGreek[1] = "4\'"
                nextAtomName = "C3\'"
                break
            case "C3\'" :
                for (var i = 1; i <= cSet.size; i++) {
                    if (cSet[i].element == "O") {
                        s[1] = i
                        newGreek[1] = "4\'"
                        cSet[i].selected = false
                        stopIdx = cSet[i].atomIndex
                    }
                    else {
                        s[2] = i
                        newGreek[2] = "3\'"
                        nIdx = cSet[i].atomIndex

                    }
                }
                nextAtomName = "O3\'"
                break
            case "O3\'" :
                for (var i = 1; i <= cSet.size; i++) {
                    if (cSet[i].element == "O") {
                        s[1] = i
                        newGreek[1] = "3\'"
                        o3pIdx = cSet[i].atomIndex
                    }
                    else {
                        s[2] = i
                        newGreek[2] = "2\'"
                        nIdx = cSet[i].atomIndex
                    }
                }
                nextAtomName = "C2\'"
                break
            case "C2\'" :
                pIdx = -1
                for (var i = 1; i <= cSet.size; i++) {
                    if (cSet[i].element == "P") {
                        pIdx = cSet[i].atomIndex
                        cSet = cSet and not cSet[i]
                        break
                    }
                }
                for (var i = 1; i <= cSet.size; i++) {
                    if (cSet[i].element == "O") {
                        s[1] = i
                        newGreek[1] = "2\'"
                        isRNA = true
                    }
                    else {
                        s[2] = i
                        newGreek[2] = "1\'"
                        c1pIdx = cSet[i].atomIndex
                        nIdx = cSet[i].atomIndex
                    }
                }
                nextAtomName = "C1\'"
                break
            case "C1\'" :
                for (var i = 1; i <= cSet.size; i++) {
                    if (cSet[i].element == "N") {
                        iKeep = i
                        nIdx = cSet[i].atomIndex
                    }
                    else if ((cSet[i].element == "C") and
                        ((connected(cSet[i]) and {oxygen}).size == 0)) { #PSU
                        psu = true
                        iKeep = i
                        nIdx = cSet[i].atomIndex
                    }
                    else {
                        cSet[i].selected = false
                    }
                }
                cSet = cSet[iKeep]
                var ccSet = connected(cSet[1]) and not {atomIndex=c1pIdx}
                newGreek[1] = "9"
                nextAtomName = "N9u"
                newGroup = "PU"
                for (var j = 1; j <= ccSet.size; j++) {
                    if ((connected(ccSet[j]) and {oxygen}) > 0) {
                        newGreek[1] = "1"
                        nextAtomName = "N1y"
                        newGroup = "PY"
                    }
                }
                break
            case "N1y" :
                for (var i = 1; i <= cSet.size; i++) {
                    if (connected(cSet[i]) > 2) {
                        iKeep = i
                        nIdx = cSet[i].atomIndex
                    }
                    else {
                        stopIdx = cSet[i].atomIndex
                    }
                }
                cSet = cSet[iKeep]
                newGreek[1] = "2"
                nextAtomName = "C2"
                break
            case "N9u" :
                # Find N-C-N-C-N
                for (var i = 1; i <= cSet.size; i++) {
                    var n1atom = (connected(cSet[i]) and {nitrogen}
                        and not {atomIndex=nIdx})
                    var c2set = connected(n1atom) and {carbon} and not cSet[i]
                    for (var j = 1; j <= c2set.size; j++) {
                        if ((connected(c2set[j]) and {nitrogen}) > 1) {
                            iDrop = i
                        }
                    }
                }
                stopIdx = cSet[iDrop].atomIndex
                cSet = cSet and not cSet[iDrop]
                nIdx = cSet[1].atomIndex
                newGreek[1] = "8"
                nextAtomName = "C8"
                break
            case "C8" :
                nIdx = ring_common( cSet, nIdx, s, newGreek, "7")
                nextAtomName = "N7"
                break
            case "N7" :
                nIdx = ring_common( cSet, nIdx, s, newGreek, "5")
                nextAtomName = "C5"
                break
            case "C5" :
                if (isRNA and (newGroup == "DU ")) {
                    var c5set = {atomIndex=nIdx} or connected({atomIndex=nIdx})
                    if (angle(c5set[1], c5set[2], c5set[3]) < 114.0) {
                        newGroup = "D  "
                    }
                }
                nIdx = ring_common( cSet, nIdx, s, newGreek, "6")
                if ((newGroup == "DU ") and (cSet.size > 1)) {
                    newGroup = "DT "
                }
                nextAtomName = "C6"
                break
            case "C6" :
                if (newGroup == "PU") {
                    nIdx = ring_common( cSet, nIdx, s, newGreek, "1")
                    newGroup = ((cSet[1].element == "O") ? "DG " : "DA ")
                    nextAtomName = "N1"
                }
                else {
                    if (psu) {
                        psu = false
                        newGroup = "DU "
                    }
                    cSet = ({})
                }
                break
            case "N1" :
                if (connected({atomIndex=nIdx}).size > 2) { # YG
                    newGroup = "X  "
                }
                nIdx = ring_common( cSet, nIdx, s, newGreek, "2")
                nextAtomName = "C2"
                break
            case "C2" :
                nIdx = ring_common( cSet, nIdx, s, newGreek, "3")
                nextAtomName = "N3"
                stopIdx = -1
                break
            case "N3" :
                nIdx = ring_common( cSet, nIdx, s, newGreek, "4")
                nextAtomName = "C4"
                break
            case "C4" :
                if (newGroup != "PY") {
                    cSet = ({})
                }
                else {
                    nIdx = ring_common( cSet, nIdx, s, newGreek, "5")
                    newGroup = ((cSet[1].element == "N") ? "DC " : "DU ")
                    nextAtomName = "C5"
                }
                break
            }
            first = false
            for (var i = 1; i <= cSet.size; i++) {
                rs += format("ATOM  %5d  %-4sUNK ", newAtomNo,
                    (cSet[i].element + newGreek[s[i]]))
                rs += format("%s%4d    %8.3f", newChain, newResno, cSet[i].x)
                rs += format("%8.3f%8.3f\n", cSet[i].y, cSet[i].z)
                newAtomno++
            }

            cSet = (connected(cSet and not {atomIndex=stopIdx})
                and not cSet and not {atomIndex=stopIdx} and not {atomIndex=endIdx})
            endIdx = nIdx

            if (cSet.size == 0) {
                if (isRNA) {
                    newGroup = (newGroup.replace("DA ","A  ").replace("DG ","G  ")
                        .replace("DC ","C  ").replace("DT ","T  ").replace("DU ","U  "))
                }
                ls += rs.replace("UNK", newGroup)
                rs = ""
                newResno++

                if (pIdx >= 0) {
                    cSet = {atomIndex=pIdx}
                    nextAtomName = "P"
                    newGroup = "UNK"
                    newGreek[1] = ""
                    c1pIdx = -1
                    stopIdx = o3pIdx
                    endIdx = -1
                    isRNA = false
                }
                else {
                    break
                }
            }
        } # endwhile

        # Replace chain with new chain
        cset = {atomIndex=idx}
        select cSet
        while (cSet) {
            cSet = connected({selected}) and not {selected}
            select {selected} or cSet
        }
        delete {selected}

        ls += "end \"append remapNT\""
        gAppendNew = appendNew
        set appendNew false
        script inline @{ls}
        set appendNew gAppendNew        
    }
    else {
        color {selected} @gScheme
    }
    print_1c_chain( newChain)
}

function print_1c_chain(iChain) {
    var resmin = get_resno_min(iChain)
    var resmax = get_resno_max(iChain)
    var rchar = (({(resno=resmin) and (chain=iChain) and thisModel}.group[0].size > 1) ? "" : "R")
    var lcAtoms = (within(3.1, false, {(resno=resmin) and (chain=iChain)
        and thisModel and base})
        and not {(resno=resmin) and (chain=iChain) and thisModel})
    var chain2 = ""
    var schar = "S"
    if (lcAtoms) {
        chain2 = lcAtoms[1].chain
        if (((rchar == "R") and (lcAtoms[1].group.size > 1))
            or ((rchar == "") and (lcAtoms[1].group.size == 1))) {
            schar = "M"
        }
        else {
            schar = "S"
            chain2 = ""
        }
    }
    var ls = format("%s%s:%s", iChain, chain2, format("%s%s", rchar, schar))
    for (var i = get_resno_min(iChain); i <= get_resno_max(iChain); i++) {
        ls += ({(resno=i) and (chain=iChain) and thisModel}.group[0])[0]
    }
    print ls
}

# Top level of Remap
function plico_remap_nt() {

    # Load common functions if not already
    if (kCommon < 7) {
        script $SCRIPT_PATH$plicoCommon.spt
        if (kCommon < 7) {
            prompt ("A newer version of plicoCommon.SPT is required")
            quit
        }
    }

    gPlico = "RemapNT"
    plico_prelim(false, true)
    
    gEcho = "_____REMAP NT_____|ALT-CLICK=select NT chain|SHIFT-DOUBLE-CLICK=exit"
    echo @gEcho

    set picking ON
    bind "ALT-LEFT-CLICK" "_pickAtom";
    bind "ALT-LEFT-CLICK" "+:remap_nt_cargo_mb";
    bind "SHIFT-DOUBLE" "plico_exit(true)";
    bind "LEFT-CLICK" "+:plico_menu_toggle";
}

# End of REMAPNT.SPT

Contributors

Remig