contactmaps 3.14 compatible; black/isort

Author: mimakaev

examples/customIntegrators/corr_noise.py

@@ -22,6 +22,7 @@
 >>> python corr_noise.py [gpuid]
 
 """
+
 import os
 import sys
 import time

examples/example/example.py

@@ -5,7 +5,6 @@
 In this simulation, a simple polymer chain of 10,000 monomers is simulated.
 """
 
-
 import os
 import sys
 

examples/loopExtrusion/directionalModel/flagshipNormLifetime.py

@@ -3,26 +3,25 @@
 
 matplotlib.use("Agg")
 
+import ctypes
+import multiprocessing as mp
 import os
-import numpy as np
+
 import matplotlib.pyplot as plt
+import numpy as np
 import pandas as pd
-import ctypes
-from mirnylib.plotting import nicePlot
-import multiprocessing as mp
 import pyximport
-
+from mirnylib.plotting import nicePlot
 from openmmlib import polymerutils
-
 from openmmlib.polymerutils import scanBlocks
 
 pyximport.install()
+import sys
+from contextlib import closing
+
 from mirnylib.h5dict import h5dict
+from mirnylib.numutils import coarsegrain, completeIC, zoomArray
 from mirnylib.systemutils import fmap, setExceptionHook
-from mirnylib.numutils import coarsegrain, completeIC
-from mirnylib.numutils import zoomArray
-from contextlib import closing
-import sys
 from smcTranslocator import smcTranslocatorDirectional
 
 filename = "/net/levsha/share/nezar/ctcf_sites/GM12878.ctcf_narrowPeak.loj.encodeMotif.rad21.txt"
@@ -70,9 +69,7 @@ def averageContacts(contactFunction, inValues, N, **kwargs):
 
     filenameChunks = [filenames[i::nproc] for i in range(nproc)]
 
-    with closing(
-        mp.Pool(processes=nproc, initializer=init, initargs=sharedARrays_)
-    ) as p:
+    with closing(mp.Pool(processes=nproc, initializer=init, initargs=sharedARrays_)) as p:
         p.map(worker, filenameChunks)
 
 
@@ -201,9 +198,7 @@ def doSim(i):
     hicdata = np.clip(hicdata, 0, np.percentile(hicdata, 99.99))
     hicdata /= np.mean(np.sum(hicdata, axis=1))
 
-    fmap(
-        doSim, range(30), n=20
-    )  # number of threads to use.  On a 20-core machine I use 20.
+    fmap(doSim, range(30), n=20)  # number of threads to use.  On a 20-core machine I use 20.
 
     arr = coarsegrain(arr, 20)
     arr = np.clip(arr, 0, np.percentile(arr, 99.9))
@@ -249,9 +244,10 @@ def calculateAverageLoop():
 
 
 def doPolymerSimulation(steps, dens, stiff, folder):
+    import time
+
     from openmmlib.openmmlib import Simulation
     from openmmlib.polymerutils import grow_rw
-    import time
 
     SMCTran = initModel()
 
@@ -322,7 +318,5 @@ def doPolymerSimulation(steps, dens, stiff, folder):
     steps=5000,
     stiff=2,
     dens=0.2,
-    folder="flagshipLifetime{1}Mu{2}_try={0}".format(
-        sys.argv[1], sys.argv[2], sys.argv[3]
-    ),
+    folder="flagshipLifetime{1}Mu{2}_try={0}".format(sys.argv[1], sys.argv[2], sys.argv[3]),
 )

examples/loopExtrusion/directionalModel/showContactmap.py

@@ -1,23 +1,21 @@
-import matplotlib
-
-# matplotlib.use("Agg")
-
-from mirnylib.plotting import nicePlot
 import os
 import pickle
-from openmmlib import contactmaps
-from mirnylib.numutils import zoomArray
-from openmmlib import polymerutils
 
+import matplotlib
 import matplotlib.pyplot as plt
 import numpy as np
-from mirnylib.h5dict import h5dict
+import pandas as pd
 from mirnylib.genome import Genome
-from mirnylib.numutils import completeIC, coarsegrain
+from mirnylib.h5dict import h5dict
+from mirnylib.numutils import coarsegrain, completeIC, zoomArray
+from mirnylib.plotting import nicePlot
 from mirnylib.systemutils import setExceptionHook
+from openmmlib import contactmaps, polymerutils
 from openmmlib.contactmapManager import averageContacts
-import pandas as pd
-from mirnylib.numutils import coarsegrain
+
+# matplotlib.use("Agg")
+
+
 
 setExceptionHook()
 
@@ -33,8 +31,8 @@ def __init__(self, i, forw, rev, blocks, steps):
         import pyximport
 
         pyximport.install()
-        from smcTranslocatorDirectional import smcTranslocator
         import numpy as np
+        from smcTranslocatorDirectional import smcTranslocator
 
         N = len(forw)
         birthArray = np.zeros(N, dtype=np.double) + 0.1
@@ -365,16 +363,8 @@ def showCmapNew():
                 if end > 75000:
                     continue
                 plt.xlim([st, end])
-                plt.savefig(
-                    os.path.join(
-                        "heatmaps", "{0}_st={1}_end={2}_r=2.png".format(fname, st, end)
-                    )
-                )
-                plt.savefig(
-                    os.path.join(
-                        "heatmaps", "{0}_st={1}_end={2}_r=2.pdf".format(fname, st, end)
-                    )
-                )
+                plt.savefig(os.path.join("heatmaps", "{0}_st={1}_end={2}_r=2.png".format(fname, st, end)))
+                plt.savefig(os.path.join("heatmaps", "{0}_st={1}_end={2}_r=2.pdf".format(fname, st, end)))
         plt.clf()
 
     plt.show()

polychrom/contactmaps.py

@@ -281,6 +281,10 @@ def worker(x):
                 return
 
 
+def identity(x):
+    return x
+
+
 def averageContacts(contactIterator, inValues, N, **kwargs):
     """
     A main workhorse for averaging contacts on multiple cores into one shared contact
@@ -352,7 +356,7 @@ def next(self):  # actual realization of the self.next method
     contactBlock = kwargs.get("contactBlock", 5000000)
     classInitArgs = kwargs.get("classInitArgs", [])
     classInitKwargs = kwargs.get("classInitKwargs", {})
-    contactProcessing = kwargs.get("contactProcessing", lambda x: x)
+    contactProcessing = kwargs.get("contactProcessing", identity)
     finalSize = N * (N + 1) // 2
     boundaries = np.linspace(0, finalSize, bucketNum + 1, dtype=int)
     chunks = zip(boundaries[:-1], boundaries[1:])

polychrom/contrib/integrators.py

@@ -44,6 +44,7 @@
 integrators use 0.1 or below).
 
 """
+
 import numpy as np
 import openmm as mm
 from openmmtools import utils

polychrom/hdf5_format.py

@@ -91,6 +91,7 @@
 `/path/to/the/trajectory/blocks_x-y.h5::42` automatically
 
 """
+
 import glob
 import os
 import warnings

polychrom/legacy/forces.py

@@ -1,9 +1,10 @@
+import os
+import pickle
+
 import numpy as np
 import simtk.openmm as openmm
-import simtk.unit.nanometer as nm
 import simtk.unit as units
-import os
-import pickle
+import simtk.unit.nanometer as nm
 
 """
 This is a collection of old forces that are likely no longer used
@@ -186,7 +187,6 @@ def lennard_jones_force(
     sigmaRep=None,
     sigmaAttr=None,
 ):
-
     """
     Adds a lennard-jones force, that allows for mutual attraction.
     This is the slowest force out of all repulsive.
@@ -372,7 +372,6 @@ def exclude_sphere(sim_object, r=5, position=(0, 0, 0)):
 
 
 def attraction_to_the_core(sim_object, k, r0, coreParticles=[]):
-
     """Attracts a subset of particles to the core,
     repells the rest from the core"""
 

polychrom/param_units.py

@@ -1,4 +1,4 @@
-""" 
+"""
 Simulation parameters and the Rouse model
 -----------------------------------------
 
@@ -9,19 +9,19 @@
 to fit experimental data on chromatin dynamics, it provides some intuition for how
 to guess polychrom parameters from first principles.
 
-For example, the mass of a particle is automatically set to 100 amu and 
+For example, the mass of a particle is automatically set to 100 amu and
 the collision rate is set depending on the integrator. For Brownian integrators,
 setting the collision rate to 2.0 works well, whereas for Langevin integrators,
 it is best to set the collision rate to 0.001-0.01. For loop extrusion simulations,
 it is typically set to 0.1. These values have been determined empirically by group
 members. For Brownian dynamics, for example, using a smaller collision rate leads
-to integration failures. By default the temperature is set to 300K. 
+to integration failures. By default the temperature is set to 300K.
 
 In any case, the monomer diffusion coefficient naturally follows as
 :math:`D=k_B T / m \zeta`, where :math:`\zeta` is the collision rate. So once the
 mass and collision rate are set, there is no control over the choice of D unless
 you use a custom Brownian integrator (see integrators module) that directly takes
-in D as a parameter. However, even then, the units of D are in terms of 
+in D as a parameter. However, even then, the units of D are in terms of
 :math:`k_B T / m \zeta`.
 
 The other arbitrary parameters are the monomer radius, set to sim.conlen = 1 nm
@@ -31,26 +31,26 @@
 :math:`k = 2k_B T / x^2`. For a Rouse chain, :math:`k = 3 k_B T / y^2`, where
 :math:`y` is the standard deviation of the bond extension. The Rouse model of DNA
 is secretely a wormlike chain at shorter length scales; thus, :math:`y` should
-be set to the end-to-end distance of the underlying WLC, which is defined as 
+be set to the end-to-end distance of the underlying WLC, which is defined as
 :math:`y = \sqrt{L_0 b}`, where :math:`L_0` is the length of DNA per monomer and
 :math:`b` is the Kuhn length. These relations imply that the bondWiggleDistance
-should be set to :math:`x = \sqrt{2L_0 b / 3}`. 
+should be set to :math:`x = \sqrt{2L_0 b / 3}`.
 
 Note that all length scales in
 polychrom are in terms of sim.conlen = 1 nm. So :math:`L_0` and :math:`b` should
 also be in nanometers. It is not obvious how to convert from nanometers of
 chromatin to basepairs. One way of doing it is using the formula
 n_basepairs = (n_nm / 0.34 nm/bp) * (1 + 146/<L>), where <L> is the average linker
-length in the cell. For example, for human T cells, <L> = 50 bp; so b=40nm of 
+length in the cell. For example, for human T cells, <L> = 50 bp; so b=40nm of
 cumulative linker length would translate to 469 bp of chromatin, where we account
 for the buried DNA in nucleosomes. We can also invert this relation to get
 n_nm = n_basepairs/(1 + 146/<L>) * 0.34 nm/bp. So if we would like each monomer
 to represent 2 kilobases, this would translate to 174 nm of cumulative linker length.
 Using these values as an example, the resulting bondWiggleDistance would be
 sqrt(2Lb/3) = 68 nm. We then divide by the size of a monomer to understand what the
 bondWiggleDistance would be in terms of sim.conlen. So if we posit that the diameter
-of a bead is equal to 34 nm, the bondWiggleDistance would be close to 2.0. This is 
-a much larger number than 0.1, which is what is used as default in polychrom! 
+of a bead is equal to 34 nm, the bondWiggleDistance would be close to 2.0. This is
+a much larger number than 0.1, which is what is used as default in polychrom!
 
 Generally, the more coarse-grained the simulation is, the more flexible the springs
 should be and the larger the bondWiggleDistance should be, since there is more DNA
@@ -61,15 +61,14 @@
 
 For a self avoiding polymer there are 2 main dimensionless numbers to keep in mind.
 One is Ddt/b^2 and the other is a/b, where a is the radius of a monomer and b is the
-Kuhn length. Recall that D and dt (timestep) are set arbitrarily based on 
-computational convenience, and a is set to 1 nm by default. 
+Kuhn length. Recall that D and dt (timestep) are set arbitrarily based on
+computational convenience, and a is set to 1 nm by default.
 Thus, even if length scales and time scales can be
 rescaled at the end of the simulation to match experimental data, these ratios should
 be decided on beforehand and preserved. Most people set the rest length of the spring
-to be the diameter of the monomer = 1 nm. 
+to be the diameter of the monomer = 1 nm.
 """
 
-
 import numpy as np
 from simtk import unit
 

polychrom/polymer_analyses.py

@@ -37,8 +37,8 @@
 
 import numpy as np
 import pandas as pd
-from scipy.spatial import cKDTree
 from scipy.ndimage import gaussian_filter1d
+from scipy.spatial import cKDTree
 
 try:
     from . import _polymer_math