Signal plot highlightning and bug corrections.

- Fixed issue where the option -p would not plot the KK layout in layout.py.
- get_kk_plot() has now an argument to draw a circle of radious neighbourhood.
signal_plot() is now working. Argument -m is a flag to select highlightning
of the signal plots. If True, only the neighbouring bins from the point of
interest will be highlighted (independently of the quadrant and significance
of those bins, but only if the point of interest is significant). If False,
all significant regions that correspond to the quadrant selected with -q will
be highlighted (default: False).
- Fixed MetalociObject.kk_coords() that was saving moran_index instead of
bin_index. That was requiring some stupid workaround in the plot code,
but not anymore.
- Fixed signal not parsing correctly in load_region_signal() when provided
a signal string as an argument instead of a file with signal names.
- Changed file format for bed output in figure.py.

metaloci/VERSION

@@ -1 +1 @@
-0.2.5
+0.2.8

metaloci/plot/plot.py

@@ -14,13 +14,12 @@
 from scipy.ndimage import rotate
 from scipy.stats import linregress, zscore
 from shapely.geometry import Point
-from shapely.geometry.multipolygon import MultiPolygon
 
 from metaloci import mlo
 from metaloci.misc import misc
 
 
-def get_kk_plot(mlobject: mlo.MetalociObject, restraints: bool = True):
+def get_kk_plot(mlobject: mlo.MetalociObject, restraints: bool = True, neighbourhood: float = None):
  """
  Generate Kamada-Kawai plot from pre-calculated restraints.
 
@@ -36,7 +35,6 @@ def get_kk_plot(mlobject: mlo.MetalociObject, restraints: bool = True):
  matplotlib.pyplot.figure.Figure
  Kamada-Kawai layout plot object.
  """
-
  PLOTSIZE = 10
  POINTSIZE = PLOTSIZE * 5
 
@@ -48,7 +46,7 @@ def get_kk_plot(mlobject: mlo.MetalociObject, restraints: bool = True):
  options = {"node_size": 50, "edge_color": "black", "linewidths": 0.1, "width": 0.05}
 
  if restraints == True:
-
+ 
  nx.draw(
  mlobject.kk_graph,
  mlobject.kk_nodes,
@@ -69,8 +67,15 @@ def get_kk_plot(mlobject: mlo.MetalociObject, restraints: bool = True):
  g.annotate(f" {mlobject.chrom}:{mlobject.start}", (xs[0], ys[0]), size=8)
  g.annotate(f" {mlobject.chrom}:{mlobject.end}", (xs[len(xs) - 1], ys[len(ys) - 1]), size=8)
 
+ poi_x, poi_y = xs[mlobject.poi], ys[mlobject.poi]
+
+ if neighbourhood:
+
+ circle = plt.Circle((poi_x, poi_y), neighbourhood, color="red", fill=False, lw=1, zorder=3)
+ plt.gca().add_patch(circle)
+
  sns.scatterplot(
- x=[xs[mlobject.poi]], y=[ys[mlobject.poi]], s=POINTSIZE * 1.5, ec="lime", fc="none", zorder=3
+ x=[poi_x], y=[poi_y], s=POINTSIZE * 1.5, ec="lime", fc="none", zorder=4
  )
 
  return kk_plt
@@ -207,17 +212,14 @@ def get_gaudi_signal_plot(mlobject: mlo.MetalociObject, lmi_geometry: pd.DataFra
  matplotlib figure containing the Gaudi signal plot.
  """
 
- poi = lmi_geometry.loc[lmi_geometry["moran_index"] == mlobject.poi, "bin_index"].iloc[0]
-
  cmap = "PuOr_r"
  min_value = lmi_geometry.signal.min()
  max_value = lmi_geometry.signal.max()
 
  gsp, ax = plt.subplots(figsize=(12, 10), subplot_kw={"aspect": "equal"})
  lmi_geometry.plot(column="signal", cmap=cmap, linewidth=2, edgecolor="white", ax=ax)
 
- sns.scatterplot(x=[lmi_geometry.X[poi]], y=[lmi_geometry.Y[poi]], s=50, ec="none", fc="lime")
-
+ sns.scatterplot(x=[lmi_geometry.X[mlobject.poi]], y=[lmi_geometry.Y[mlobject.poi]], s=50, ec="none", fc="lime")
  sm = plt.cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=min_value, vmax=max_value))
  sm.set_array([1, 2, 3, 4])
 
@@ -262,7 +264,6 @@ def get_gaudi_type_plot(
  gtp : matplotlib.pyplot.figure.Figure
  matplotlib figure containing the Gaudi type plot.
  """
- poi = lmi_geometry.loc[lmi_geometry["moran_index"] == mlobject.poi, "bin_index"].iloc[0]
 
  legend_elements = [
  Line2D([0], [0], marker="o", color="w", markerfacecolor=colors[1], label="HH", markersize=20),
@@ -272,15 +273,15 @@ def get_gaudi_type_plot(
  ]
 
  cmap = LinearSegmentedColormap.from_list("Custom cmap", [colors[nu] for nu in colors], len(colors))
- alpha = [1.0 if quad <= signipval else 0.3 for quad in lmi_geometry.LMI_pvalue]
+ alpha = [1.0 if pval <= signipval else 0.3 for pval in lmi_geometry.LMI_pvalue]
 
  gtp, ax = plt.subplots(figsize=(12, 10), subplot_kw={"aspect": "equal"})
  lmi_geometry.plot(column="moran_quadrant", cmap=cmap, alpha=alpha, linewidth=2, edgecolor="white", ax=ax)
  plt.axis("off")
 
  sns.scatterplot(
- x=[lmi_geometry.X[poi]],
- y=[lmi_geometry.Y[poi]],
+ x=[lmi_geometry.X[mlobject.poi]],
+ y=[lmi_geometry.Y[mlobject.poi]],
  s=50,
  ec="none",
  fc="lime",
@@ -368,13 +369,102 @@ def get_lmi_scatterplot(
 
  return scatter, r_value_scat, p_value_scat
 
+# def signal_bed(
+# mlobject: mlo.MetalociObject,
+# lmi_geometry: pd.DataFrame,
+# neighbourhood: float,
+# quadrants: list = [1, 3],
+# signipval: float = 0.05,
+# ):
+# """
+# signal_bed _summary_
+
+# Parameters
+# ----------
+# mlobject : mlo.MetalociObject
+# _description_
+# lmi_geometry : pd.DataFrame
+# _description_
+# neighbourhood : float
+# _description_
+# quadrants : list, optional
+# _description_, by default [1, 3]
+# signipval : float, optional
+# _description_, by default 0.05
+
+# Returns
+# -------
+# _type_
+# _description_
+# """
+
+# poi_distance = mlobject.kk_distances[mlobject.poi]
+
+# # Select the polygons that are in the quadrant of interest and are significative.
+# ml_indexes = lmi_geometry[
+# (lmi_geometry.moran_quadrant.isin(quadrants)) & (lmi_geometry.LMI_pvalue <= signipval)
+# ].bin_index.values.tolist()
+
+# # Make a big polygon from the small poligons that are significant.
+# metalocis = lmi_geometry[lmi_geometry.bin_index.isin(ml_indexes)].unary_union
+
+# if metalocis and metalocis.geom_type == "Polygon":
+
+# metalocis = MultiPolygon([metalocis]) # Need a multipolygon in order for the code to work.
+
+# poi_point = Point(
+# (lmi_geometry[lmi_geometry.moran_index == mlobject.poi].X, lmi_geometry[lmi_geometry.moran_index == mlobject.poi].Y)
+# )
+
+# metalocis_bed = []
+
+# bed_data = defaultdict(list)
+
+# try:
+
+# for metaloci in metalocis.geoms:
+
+# metaloci = gpd.GeoSeries(metaloci)
+
+# if poi_point.within(metaloci[0]):
+
+# for _, row_ml in lmi_geometry.iterrows():
+
+# adjacent_point = Point((row_ml.X, row_ml.Y))
+
+# if adjacent_point.within(metaloci[0]):
+
+# metalocis_bed.append(row_ml.bin_index)
+
+# # Add close particles
+# metalocis_bed.sort()
+# close_bins = [i for i, distance in enumerate(poi_distance) if distance <= neighbourhood / 2]
+# metalocis_bed = np.sort(list(set(close_bins + metalocis_bed)))
+
+# for point in metalocis_bed:
+
+# bed_data["chr"].append(lmi_geometry.bin_chr[point])
+# bed_data["start"].append(lmi_geometry.bin_start[point])
+# bed_data["end"].append(lmi_geometry.bin_end[point])
+# bed_data["bin"].append(point)
+
+# except:
+
+# pass
+
+# bed_data = pd.DataFrame(bed_data)
+
+# return bed_data, metalocis_bed
+
+
 
 def signal_bed(
  mlobject: mlo.MetalociObject,
  lmi_geometry: pd.DataFrame,
  neighbourhood: float,
  quadrants: list = [1, 3],
  signipval: float = 0.05,
+ metaloci_only: bool = False,
 ):
  """
  signal_bed _summary_
@@ -391,72 +481,66 @@ def signal_bed(
  _description_, by default [1, 3]
  signipval : float, optional
  _description_, by default 0.05
+ metaloci_only : bool, optional
+ _description_, by default False
 
  Returns
  -------
  _type_
  _description_
  """
 
- poi_distance = mlobject.kk_distances[mlobject.poi]
+ metalocis_bed = []
 
- # Select the polygons that are in the quadrant of interest and are significative.
- ml_indexes = lmi_geometry[
- (lmi_geometry.moran_quadrant.isin(quadrants)) & (lmi_geometry.LMI_pvalue <= signipval)
- ].bin_index.values.tolist()
+ bed_data = defaultdict(list)
 
- # Make a big polygon from the small poligons that are significant.
- metalocis = lmi_geometry[lmi_geometry.bin_index.isin(ml_indexes)].unary_union
+ if metaloci_only == False:
 
- if metalocis and metalocis.geom_type == "Polygon":
+  for _, row_ml in lmi_geometry.iterrows():
 
- metalocis = MultiPolygon([metalocis]) # Need a multipolygon in order for the code to work.
+  if row_ml.LMI_pvalue <= signipval and row_ml.moran_quadrant in quadrants:
 
- poi_point = Point(
- (lmi_geometry[lmi_geometry.bin_index == mlobject.poi].X, lmi_geometry[lmi_geometry.bin_index == mlobject.poi].Y)
- )
+ metalocis_bed.append(row_ml.bin_index)
 
- metalocis_bed = []
+ else:
 
- bed_data = defaultdict(list)
+  poi_row = lmi_geometry[lmi_geometry.bin_index == mlobject.poi].squeeze()
 
- try:
+  # print(lmi_geometry[lmi_geometry.bin_index == mlobject.poi])
 
- for metaloci in metalocis.geoms:
+ if poi_row.LMI_pvalue <= signipval and poi_row.moran_quadrant in quadrants:
 
- metaloci = gpd.GeoSeries(metaloci)
+ poi_point = Point(
+ (lmi_geometry[lmi_geometry.bin_index == mlobject.poi].X, lmi_geometry[lmi_geometry.bin_index == mlobject.poi].Y)
+ )
 
- if poi_point.within(metaloci[0]):
+ poi_point = Point(lmi_geometry.loc[lmi_geometry["bin_index"] == mlobject.poi, ["X", "Y"]].iloc[0])
 
-  for _, row_ml in lmi_geometry.iterrows():
+ for _, row_ml in lmi_geometry.iterrows():
 
-  adjacent_point = Point((row_ml.X, row_ml.Y))
+ adjacent_point = Point((row_ml.X, row_ml.Y))
 
-  if adjacent_point.within(metaloci[0]):
+ if adjacent_point.distance(poi_point) <= neighbourhood:
 
-  metalocis_bed.append(row_ml.bin_index)
+ metalocis_bed.append(row_ml.bin_index)
 
- # Add close particles
- metalocis_bed.sort()
- close_bins = [i for i, distance in enumerate(poi_distance) if distance <= neighbourhood / 2]
- metalocis_bed = np.sort(list(set(close_bins + metalocis_bed)))
+ metalocis_bed.sort()
 
-  for point in metalocis_bed:
+ # print(metalocis_bed)
 
- bed_data["chr"].append(lmi_geometry.bin_chr[point])
- bed_data["start"].append(lmi_geometry.bin_start[point])
- bed_data["end"].append(lmi_geometry.bin_end[point])
- bed_data["bin"].append(point)
+ for point in metalocis_bed:
 
- except:
+ bed_data["chr"].append(lmi_geometry.bin_chr[point])
+ bed_data["start"].append(lmi_geometry.bin_start[point])
+ bed_data["end"].append(lmi_geometry.bin_end[point])
+ bed_data["bin"].append(point)
+ bed_data["quadrant"].append(lmi_geometry.moran_quadrant[point])
 
- pass
 
  bed_data = pd.DataFrame(bed_data)
 
  return bed_data, metalocis_bed
 
-
 def signal_plot(mlobject: mlo.MetalociObject, lmi_geometry: pd.DataFrame, metalocis, bins_sig, coords_sig):
 
  sig_plt = plt.figure(figsize=(10, 1.5))

metaloci/spatial_stats/lmi.py

@@ -18,7 +18,6 @@
 from metaloci import mlo
 from metaloci.misc import misc
 
-
 warnings.filterwarnings("ignore", category=ShapelyDeprecationWarning)
 
 
@@ -85,8 +84,8 @@ def construct_voronoi(mlobject: mlo.MetalociObject, buffer: float):
 
  df_geometry["bin_index"].append(x)
  df_geometry["moran_index"].append(i)
- df_geometry["X"].append(mlobject.kk_coords[i][0])
- df_geometry["Y"].append(mlobject.kk_coords[i][1])
+ df_geometry["X"].append(mlobject.kk_coords[x][0])
+ df_geometry["Y"].append(mlobject.kk_coords[x][1])
  df_geometry["geometry"].append(geometry_data.loc[i, "geometry"])
 
  df_geometry = pd.DataFrame(df_geometry)
@@ -197,9 +196,16 @@ def load_region_signals(mlobject: mlo.MetalociObject, signal_data: dict, signal_
  """
 
  # Read signal file. Will only process the signals present in this list.
- with open(signal_file) as signals_handler:
 
- signal_types = [line.rstrip() for line in signals_handler]
+ if os.path.isfile(signal_file):
+
+ with open(signal_file) as signals_handler:
+
+ signal_types = [line.rstrip() for line in signals_handler]
+
+ else:
+
+ signal_types = [signal_file]
 
  # region_signal = signal_data[mlobject.chrom][
  # (signal_data[mlobject.chrom]["start"] >= int(mlobject.start / mlobject.resolution) * mlobject.resolution)