Python Molecule.filter方法代码示例

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import filter [as 别名]
    def test_disulfideWithInsertion(self):
        from htmd.molecule.molecule import Molecule
        from htmd.builder.solvate import solvate
        from htmd.home import home
        from htmd.util import tempname, assertSameAsReferenceDir
        import os
        import numpy as np

        # Use pre-prepared files so we can tell whether the error is in prepare or in build
        # Inputs are reference outputs of proteinprepare.
        preparedInputDir = home(dataDir='test-proteinprepare')

        pdb = '3PTB'

        print('Building {}'.format(pdb))
        inFile = os.path.join(preparedInputDir, pdb, "{}-prepared.pdb".format(pdb))
        mol = Molecule(inFile)
        mol.filter('protein')  # Fix for bad proteinPrepare hydrogen placing

        np.random.seed(1)  # Needed for ions
        smol = solvate(mol)
        topos = ['top/top_all36_prot.rtf', 'top/top_water_ions.rtf']
        params = ['par/par_all36_prot_mod.prm', 'par/par_water_ions.prm']

        smol.insertion[smol.resid == 42] = 'A'  # Adding an insertion to test that disulfide bonds with insertions work
        tmpdir = tempname()
        _ = build(smol, topo=topos, param=params, outdir=tmpdir)
        compareDir = home(dataDir=os.path.join('test-charmm-build', '3PTB_insertion'))
        assertSameAsReferenceDir(compareDir, tmpdir)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import filter [as 别名]
    def test_build(self):
        from htmd.molecule.molecule import Molecule
        from htmd.builder.solvate import solvate
        from htmd.home import home
        from htmd.util import tempname, assertSameAsReferenceDir
        import os
        import numpy as np

        # Use pre-prepared files so we can tell whether the error is in prepare or in build
        # Inputs are reference outputs of proteinprepare.
        preparedInputDir = home(dataDir='test-proteinprepare')

        pdbids = ['3PTB', '1A25', '1GZM', '1U5U']
        for pdb in pdbids:
            with self.subTest(pdb=pdb):
                print('Building {}'.format(pdb))
                inFile = os.path.join(preparedInputDir, pdb, "{}-prepared.pdb".format(pdb))
                mol = Molecule(inFile)
                mol.filter('protein')  # Fix for bad proteinPrepare hydrogen placing

                np.random.seed(1)  # Needed for ions
                smol = solvate(mol)
                topos = ['top/top_all36_prot.rtf', 'top/top_water_ions.rtf']
                params = ['par/par_all36_prot_mod.prm', 'par/par_water_ions.prm']
                tmpdir = tempname()
                _ = build(smol, topo=topos, param=params, outdir=tmpdir)

                compareDir = home(dataDir=os.path.join('test-charmm-build', pdb))
                assertSameAsReferenceDir(compareDir, tmpdir)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import filter [as 别名]
def _loadMolecules(lipids, files):
    from htmd.rotationmatrix import rotationMatrix
    # Create Molecules
    for l in lipids:
        randidx = np.random.randint(len(files[l.resname]))
        mol = Molecule(files[l.resname][randidx])
        mol.filter('not water', _logger=False)
        if l.xyz[2] < 0:
            mol.rotateBy(rotationMatrix([1, 0, 0], np.deg2rad(180)))  # Rotate the lower leaflet lipids upside down
        l.mol = mol
        l.rot = np.random.random() * 360 - 180  # Random starting rotation

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import filter [as 别名]
def opm(pdb):
    """Download a molecule from the OPM.

    Removes DUM atoms.

    Parameters
    ----------
    pdb: str
        The 4-letter PDB code

    Returns
    -------
    mol: Molecule
        The oriented molecule

    thickness: float
        The bilayer thickness (both layers)

    Examples
    --------
    >>> mol, thickness = opm("1z98")
    >>> mol.numAtoms
    7902
    >>> thickness
    28.2

    """

    from htmd.molecule.support import string_to_tempfile
    from htmd.molecule.molecule import Molecule
    # http://opm.phar.umich.edu/pdb/1z98.pdb
    r = requests.get("http://opm.phar.umich.edu/pdb/{:s}.pdb".format(pdb.lower()))

    if r.status_code != 200:
        raise NameError('PDB code not found in the OPM database')

    tempfile = string_to_tempfile(r.content.decode('ascii'), "pdb")
    mol = Molecule(tempfile)
    mol.filter("not resname DUM")

    # Assuming the half-thickness is the last word in the first line
    # REMARK      1/2 of bilayer thickness:   14.1
    f = open(tempfile)
    h = f.readline()
    f.close()
    os.unlink(tempfile)

    hs = h.split()
    thickness = 2.0 * float(hs[-1])

    return mol, thickness

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import filter [as 别名]
    def testWithoutProteinPrepare(self):
        from htmd.builder.solvate import solvate
        # Test without proteinPrepare
        pdbids = ['3PTB']
        # pdbids = ['3PTB', '1A25', '1GZM', '1U5U']
        for pid in pdbids:
            np.random.seed(1)
            mol = Molecule(pid)
            mol.filter('protein')
            smol = solvate(mol)
            ffs = defaultFf()
            tmpdir = os.path.join(self.testDir, 'withoutProtPrep', pid)
            _ = build(smol, ff=ffs, outdir=tmpdir)

            refdir = home(dataDir=join('test-amber-build', 'nopp', pid))
            TestAmberBuild._compareResultFolders(refdir, tmpdir, pid)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import filter [as 别名]
    def _analyse(self, mol, pdb, rtf, prm, traj, ftraj):
        t = Molecule(pdb)
        t.read(traj)
        t.filter('not water')
        t.write(ftraj)
        m = FFMolecule(filename=mol, rtf=rtf, prm=prm)
        m.read(ftraj)
        torsions = m.getRotatableDihedrals()
        # For each torsion
        for i in range(len(torsions)):
            # Create title
            title = '{}-{}-{}-{}'.format(m.name[torsions[i][0]], m.name[torsions[i][1]], m.name[torsions[i][2]],
                                         m.name[torsions[i][3]])
            # Measure
            (r, theta) = self._measure_torsion(torsions[i], m.coords)

            self._plot_scatter(r, theta, title)
            self._plot_hist(theta, title)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import filter [as 别名]
    def testWithProteinPrepare(self):
        from htmd.builder.preparation import proteinPrepare
        from htmd.builder.solvate import solvate
        # Test with proteinPrepare
        pdbids = ['3PTB']
        # pdbids = ['3PTB', '1A25', '1GZM']  # '1U5U' out because it has AR0 (no parameters)
        for pid in pdbids:
            np.random.seed(1)
            mol = Molecule(pid)
            mol.filter('protein')
            mol = proteinPrepare(mol)
            mol.filter('protein')  # Fix for bad proteinPrepare hydrogen placing
            smol = solvate(mol)
            ffs = defaultFf()
            tmpdir = os.path.join(self.testDir, 'withProtPrep', pid)
            _ = build(smol, ff=ffs, outdir=tmpdir)

            refdir = home(dataDir=join('test-amber-build', 'pp', pid))
            TestAmberBuild._compareResultFolders(refdir, tmpdir, pid)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import filter [as 别名]
    def test_customDisulfideBonds(self):
        from htmd.builder.solvate import solvate
        # Test without proteinPrepare
        pdbids = ['1GZM', ]
        for pid in pdbids:
            np.random.seed(1)
            mol = Molecule(pid)
            mol.filter('protein')
            smol = solvate(mol)
            ffs = defaultFf()
            disu = [['segid 1 and resid 110', 'segid 1 and resid 187'],
                    ['segid 0 and resid 110', 'segid 0 and resid 187']]
            tmpdir = os.path.join(self.testDir, 'withoutProtPrep', pid)
            _ = build(smol, ff=ffs, outdir=tmpdir, disulfide=disu)

            refdir = home(dataDir=join('test-amber-build', 'nopp', pid))
            TestAmberBuild._compareResultFolders(refdir, tmpdir, pid)

            np.random.seed(1)
            tmpdir = os.path.join(self.testDir, 'withoutProtPrep', pid)
            _ = build(smol, ff=ffs, outdir=tmpdir)

            refdir = home(dataDir=join('test-amber-build', 'nopp', pid))
            TestAmberBuild._compareResultFolders(refdir, tmpdir, pid)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import filter [as 别名]
    def test_customDisulfideBonds(self):
        from htmd.molecule.molecule import Molecule
        from htmd.builder.solvate import solvate
        from htmd.home import home
        from htmd.util import tempname, assertSameAsReferenceDir
        import os
        import numpy as np

        # Use pre-prepared files so we can tell whether the error is in prepare or in build
        # Inputs are reference outputs of proteinprepare.
        preparedInputDir = home(dataDir='test-proteinprepare')

        pdb = '1GZM'
        inFile = os.path.join(preparedInputDir, pdb, "{}-prepared.pdb".format(pdb))
        mol = Molecule(inFile)
        mol.filter('protein')  # Fix for bad proteinPrepare hydrogen placing

        np.random.seed(1)  # Needed for ions
        smol = solvate(mol)
        topos = ['top/top_all36_prot.rtf', 'top/top_water_ions.rtf']
        params = ['par/par_all36_prot_mod.prm', 'par/par_water_ions.prm']
        disu = [['segid 1 and resid 110', 'segid 1 and resid 187'], ['segid 0 and resid 110', 'segid 0 and resid 187']]
        tmpdir = tempname()
        _ = build(smol, topo=topos, param=params, outdir=tmpdir, disulfide=disu)

        compareDir = home(dataDir=os.path.join('test-charmm-build', pdb))
        assertSameAsReferenceDir(compareDir, tmpdir)

        # TODO: Remove this after deprecation
        from htmd.builder.builder import DisulfideBridge
        np.random.seed(1)  # Needed for ions
        disu = [DisulfideBridge('1', 110, '1', 187), DisulfideBridge('0', 110, '0', 187)]
        tmpdir = tempname()
        _ = build(smol, topo=topos, param=params, outdir=tmpdir, disulfide=disu)
        compareDir = home(dataDir=os.path.join('test-charmm-build', pdb))
        assertSameAsReferenceDir(compareDir, tmpdir)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import filter [as 别名]
def opm(pdb, keep=False, keepaltloc='A'):
    """Download a molecule from the OPM.

    Removes DUM atoms.

    Parameters
    ----------
    pdb: str
        The 4-letter PDB code
    keep: bool
        If False, removes the DUM atoms. If True, it keeps them.

    Returns
    -------
    mol: Molecule
        The oriented molecule

    thickness: float or None
        The bilayer thickness (both layers)

    Examples
    --------
    >>> mol, thickness = opm("1z98")
    >>> mol.numAtoms
    7902
    >>> thickness
    28.2
    >>> _, thickness = opm('4u15')
    >>> thickness is None
    True

    """

    from htmd.molecule.support import string_to_tempfile
    from htmd.molecule.molecule import Molecule
    # http://opm.phar.umich.edu/pdb/1z98.pdb
    r = requests.get("http://opm.phar.umich.edu/pdb/{:s}.pdb".format(pdb.lower()))

    if r.status_code != 200:
        raise NameError('PDB code not found in the OPM database')

    tempfile = string_to_tempfile(r.content.decode('ascii'), "pdb")
    mol = Molecule(tempfile, keepaltloc=keepaltloc)
    if not keep:
        mol.filter("not resname DUM")

    # Assuming the half-thickness is the last word in the matched line
    # REMARK      1/2 of bilayer thickness:   14.1
    tmp = open(tempfile)
    pattern = re.compile('^REMARK.+thickness')

    match = None
    for line in tmp:
        if re.match(pattern, line):
            match = line
    tmp.close()
    os.unlink(tempfile)

    if not match:
        thickness = None
    else:
        split_line = match.split()
        thickness = 2.0 * float(split_line[-1])

    return mol, thickness

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import filter [as 别名]
    import os

    # No arguments - quick travis test
    if len(sys.argv) == 1:
        tryp_op, prepData = proteinPrepare(Molecule("3PTB"), returnDetails=True)
        d = prepData.data
        assert d.protonation[d.resid == 40].iloc[0] == 'HIE'
        assert d.protonation[d.resid == 57].iloc[0] == 'HIP'
        assert d.protonation[d.resid == 91].iloc[0] == 'HID'

    # Long test
    elif sys.argv[1] == "long-test":
        pdbids = ['3PTB', '1A25', '1GZM', '1U5U']
        for pdb in pdbids:
            mol = Molecule(pdb)
            mol.filter("protein")
            mol_op, prepData = proteinPrepare(mol, returnDetails=True)
            mol_op.write("./{}-prepared.pdb".format(pdb))
            prepData.data.to_csv("./{}-prepared.csv".format(pdb), float_format="%.2f")

            compareDir = htmd.home.home(dataDir=os.path.join('test-proteinprepare', pdb))
            htmd.util.assertSameAsReferenceDir(compareDir)
        import doctest
        doctest.testmod()

    # Stand-alone executable: prepare the PDB given as argument
    else:
        mol = Molecule(sys.argv[1])
        mol.filter("protein")
        mol_op, prepData = proteinPrepare(mol, returnDetails=True)
        mol_op.write("./prepared.pdb")

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import filter [as 别名]
    def cart2pol(x, y):
        # Cartesian to polar coordinates. Rho is the rotation angle
        rho = np.sqrt(x ** 2 + y ** 2)
        phi = np.arctan2(y, x)
        return rho, phi

    angle, _ = cart2pol(xa, ya)
    return angle + np.radians(45)


if __name__ == "__main__":
    from htmd import *
    from os import path

    p = Molecule(path.join(home(), 'data', 'building-protein-membrane', '4dkl.pdb'))
    p.filter('(chain B and protein) or water')
    p = autoSegment(p, 'protein', 'P')
    m = Molecule(path.join(home(), 'data', 'building-protein-membrane', 'membrane.pdb'))
    a = embed(p, m)
    print(np.unique(m.get('segid')))

    mol = Molecule(path.join(home(), 'data', 'building-protein-membrane', '1ITG_clean.pdb'))
    ref = Molecule(path.join(home(), 'data', 'building-protein-membrane', '1ITG.pdb'))
    mol = autoSegment(mol, sel='protein')
    assert np.all(mol.segid == ref.segid)

    mol = Molecule(path.join(home(), 'data', 'building-protein-membrane', '3PTB_clean.pdb'))
    ref = Molecule(path.join(home(), 'data', 'building-protein-membrane', '3PTB.pdb'))
    mol = autoSegment(mol, sel='protein')
    assert np.all(mol.segid == ref.segid)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import filter [as 别名]
        sschar[i][sschar[i] == 'G'] = 'H'
        sschar[i][sschar[i] == 'I'] = 'H'
        sschar[i][sschar[i] == 'B'] = 'E'
        sschar[i][sschar[i] == 'b'] = 'E'
        sschar[i][sschar[i] == 'S'] = 'C'
        sschar[i][sschar[i] == 'T'] = 'C'
        ssnum[i, np.where(sschar[i] == 'C')[0]] = 0
        ssnum[i, np.where(sschar[i] == 'E')[0]] = 1
        ssnum[i, np.where(sschar[i] == 'H')[0]] = 2
    return ssnum

if __name__ == "__main__":
    from htmd.molecule.molecule import Molecule
    import numpy as np
    mol = Molecule('2HBB')  # NTL9
    mol.filter('protein')
    metr = MetricSecondaryStructure()
    data = metr.project(mol)
    ref = np.array([[0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0,
                     2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 2, 2, 2, 2,
                     2, 2, 2, 2, 2, 2, 0]], dtype=np.int32)
    assert np.array_equal(data, ref), 'MetricSecondaryStructure assertion failed'
    metr = MetricSecondaryStructure('resid 5 to 49')
    data = metr.project(mol)
    ref = np.array([[0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2,
                    2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 0]], dtype=np.int32)
    assert np.array_equal(data, ref), 'MetricSecondaryStructure assertion failed'

    mol = Molecule('3PTB')  # Contains insertion which used to be a problem
    mol.filter('protein')
    metr = MetricSecondaryStructure()

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import filter [as 别名]
                files.append(os.path.basename(f))

        match, mismatch, error = filecmp.cmpfiles(tmpdir, compare, files, shallow=False)
        if len(mismatch) != 0 or len(error) != 0 or len(match) != len(files):
            raise RuntimeError(
                'Different results produced by amber.build for test {} between {} and {} in files {}.'.format(pid, compare, tmpdir, mismatch))

        for f in deletefiles:
            os.remove(f)

    # Test with proteinPrepare
    pdbids = ['3PTB', '1A25', '1GZM']  # '1U5U' out because it has AR0 (no parameters)
    for pid in pdbids:
        np.random.seed(1)
        mol = Molecule(pid)
        mol.filter('protein')
        mol = proteinPrepare(mol)
        mol.filter('protein')  # Fix for bad proteinPrepare hydrogen placing
        smol = solvate(mol)
        ffs = ['leaprc.lipid14', 'leaprc.ff14SB', 'leaprc.gaff']
        tmpdir = tempname()
        bmol = build(smol, ff=ffs, outdir=tmpdir)

        refdir = home(dataDir=os.path.join('test-amber-build', pid))
        _compareResultFolders(refdir, tmpdir, pid)
        shutil.rmtree(tmpdir)

    # Test without proteinPrepare
    pdbids = ['3PTB', '1A25', '1GZM', '1U5U']
    for pid in pdbids:
        np.random.seed(1)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import filter [as 别名]
for ext in _MDTRAJ_SAVERS:
    if ext not in _WRITERS:
        _WRITERS[ext] = MDTRAJwrite


if __name__ == '__main__':
    from htmd.home import home
    from htmd.molecule.molecule import Molecule, mol_equal
    from htmd.util import tempname
    import numpy as np
    import os
    testfolder = home(dataDir='metricdistance')
    mol = Molecule(os.path.join(testfolder, 'filtered.pdb'))
    mol.coords = np.tile(mol.coords, (1, 1, 2))
    mol.filter('protein and resid 1 to 20')
    mol.boxangles = np.ones((3, 2), dtype=np.float32) * 90
    mol.box = np.ones((3, 2), dtype=np.float32) * 15
    mol.step = np.arange(2)
    mol.time = np.arange(2) * 1E5

    for ext in _WRITERS:
        tmp = tempname(suffix='.'+ext)
        mol.write(tmp)
        print('Can write {} files'.format(ext))


    # from difflib import Differ
    # d = Differ()
    #
    # with open(tmp, 'rb') as f1, open(a, 'rb') as f2: