Merge pull request #76 from robomics/feature/stripepy-plot

Initial implementation of stripepy plot

Author: rea1991

.github/workflows/build-dockerfile.yml

@@ -8,6 +8,7 @@ on:
   push:
     branches: [main]
     paths:
+      - ".github/workflows/cache-test-datasets.yml"
       - ".github/workflows/build-dockerfile.yml"
       - "src/**"
       - "test/**"
@@ -22,6 +23,7 @@ on:
 
   pull_request:
     paths:
+      - ".github/workflows/cache-test-datasets.yml"
       - ".github/workflows/build-dockerfile.yml"
       - "src/**"
       - "test/**"
@@ -41,17 +43,14 @@ defaults:
   run:
     shell: bash
 
-env:
-  TEST_DATASET_ID: 4DNFI9GMP2J8
-  TEST_DATASET_URL: "https://zenodo.org/records/14283922/files/4DNFI9GMP2J8.stripepy.mcool?download=1"
-  TEST_DATASET_MD5: "a17d08460c03cf6c926e2ca5743e4888"
-  TEST_RESULT_FILE_V1_ID: results_4DNFI9GMP2J8_v1.hdf5
-  TEST_RESULT_FILE_V1_URL: "https://zenodo.org/records/14283922/files/results_4DNFI9GMP2J8_v1.hdf5?download=1"
-  TEST_RESULT_FILE_V1_MD5: "632b2a7a6e5c1a24dc3635710ed68a80"
-
 jobs:
+  cache-test-datasets:
+    name: Cache test dataset
+    uses: paulsengroup/StripePy/.github/workflows/cache-test-datasets.yml@main
+
   build-dockerfile:
     name: Build Dockerfile
+    needs: [cache-test-datasets]
     runs-on: ubuntu-latest
     permissions:
       contents: "read"
@@ -75,32 +74,13 @@ jobs:
       - name: Install build requirements
         run: pip install -r requirements.txt
 
-      - name: Cache test dataset
-        id: cache-dset
+      - name: Restore test dataset
         uses: actions/cache/restore@v4
         with:
-          key: test-dataset
-          path: ${{ env.TEST_DATASET_ID }}.mcool
-
-      - name: Download test dataset
-        if: steps.cache-dset.outputs.cache-hit != 'true'
-        run: |
-          mkdir -p test/data/
-          curl -L '${{ env.TEST_DATASET_URL }}' -o 'test/data/${{ env.TEST_DATASET_ID }}.mcool'
-          curl -L '${{ env.TEST_RESULT_FILE_V1_URL }}' -o 'test/data/${{ env.TEST_RESULT_FILE_V1_ID }}'
-
-          echo "${{ env.TEST_DATASET_MD5 }}  test/data/${{ env.TEST_DATASET_ID }}.mcool" > checksum.md5
-          echo "${{ env.TEST_RESULT_FILE_V1_MD5 }}  test/data/${{ env.TEST_RESULT_FILE_V1_ID }}" >> checksum.md5
-          md5sum -c checksum.md5
-
-      - name: Save test dataset
-        if: steps.cache-dset.outputs.cache-hit != 'true'
-        uses: actions/cache/save@v4
-        with:
-          key: test-dataset
-          path: |
-            test/data/${{env.TEST_DATASET_ID }}.mcool
-            test/data/results*.hdf5
+          key: ${{ needs.cache-test-datasets.outputs.cache-key }}
+          path: test/data/
+          fail-on-cache-miss: true
+          enableCrossOsArchive: true
 
       - name: Generate build args
         id: build-args

.github/workflows/ci.yml

@@ -8,6 +8,7 @@ on:
   push:
     branches: [main]
     paths:
+      - ".github/workflows/cache-test-datasets.yml"
       - ".github/workflows/ci.yml"
       - "src/**"
       - "test/**"
@@ -16,6 +17,7 @@ on:
 
   pull_request:
     paths:
+      - ".github/workflows/cache-test-datasets.yml"
       - ".github/workflows/ci.yml"
       - "src/**"
       - "test/**"
@@ -62,17 +64,18 @@ jobs:
 
             return { include: includes }
 
+  cache-test-datasets:
+    name: Cache test datasets
+    uses: paulsengroup/StripePy/.github/workflows/cache-test-datasets.yml@main
+
   ci:
     name: CI
-    needs: [matrix-factory]
+    needs: [cache-test-datasets, matrix-factory]
     strategy:
       fail-fast: false
       matrix: ${{ fromJson(needs.matrix-factory.outputs.matrix) }}
     runs-on: ${{ matrix.os }}
 
-    env:
-      TEST_DATASET_ID: 4DNFI9GMP2J8
-
     steps:
       - uses: actions/checkout@v4
         with:
@@ -93,35 +96,13 @@ jobs:
           stripepy --help
           stripepy --version
 
-      - name: Cache test dataset
-        id: cache-dset
+      - name: Restore test dataset
         uses: actions/cache/restore@v4
         with:
-          key: test-dataset
-          path: |
-            test/data/${{env.TEST_DATASET_ID }}.mcool
-            test/data/results*.hdf5
-
-      - name: Download test dataset
-        if: steps.cache-dset.outputs.cache-hit != 'true'
-        run: |
-          mkdir test/data/
-          stripepy download \
-            --name ${{ env.TEST_DATASET_ID }} \
-            --output test/data/${{ env.TEST_DATASET_ID }}.mcool
-
-          stripepy download \
-            --name __results_v1 \
-            --output test/data/results_4DNFI9GMP2J8_v1.hdf5
-
-      - name: Save test dataset
-        if: steps.cache-dset.outputs.cache-hit != 'true'
-        uses: actions/cache/save@v4
-        with:
-          key: test-dataset
-          path: |
-            test/data/${{env.TEST_DATASET_ID }}.mcool
-            test/data/results*.hdf5
+          key: ${{ needs.cache-test-datasets.outputs.cache-key }}
+          path: test/data/
+          fail-on-cache-miss: true
+          enableCrossOsArchive: true
 
       - name: Run unit tests
         run: |

.gitignore

@@ -168,6 +168,6 @@ cython_debug/
 # Custom entries
 .DS_Store
 output/
-*.log
+test/data
 *.log*
 *.bak*

README.md

@@ -70,6 +70,7 @@ StripePy is organized into a few subcommands:
 
 - `stripepy call`: run the stripe detection algorithm and store the identified stripes in a `.hdf5` file.
 - `stripepy view`: take the `result.hdf5` file generated by `stripepy call` and extract stripes in BEDPE format.
+- `stripepy plot`: generate various kinds of plots to inspect the stripes identified by `stripepy call`.
 - `stripepy download`: download a minified sample dataset suitable to quickly test StripePy.
 
 ### Walkthrough
@@ -1108,6 +1109,44 @@ Running the last cell will display a HiGlass window embedded in the Jupyter note
 
 ![HiGlass window](docs/assets/4DNFI9GMP2J8_chr2_156mbp_higlass_view.png)
 
+## Generating plots
+
+StripePy comes with a `plot` subcommand that can be used to generate various kinds of plots.
+
+`stripepy plot` supports the following subcommands:
+
+- `contact-map` (`cm`): plot stripes and other features over the Hi-C matrix
+- `pseudodistribution` (`pd`): plot the pseudo-distribution over the given region of interest
+- `stripe-hist` (`hist`): generate and plot the histograms showing the distribution of the stripe heights and widths
+
+`stripepy cm` takes as input a Hi-C matrix in `.cool`, `.mcool`, or `.hic` format, and optionally the `.hdf5` file generated by `stripepy call` (this parameter is mandatory when highlighting stripes or stripe seeds).
+
+`stripepy pd` and `stripepy hist` do not require the Hi-C matrix file, and require the `.hdf5` file generated by `stripepy call` instead.
+
+All three subcommands support specifying a region of interest through the `--region` option.
+When the commands are run without specifying the region of interest, `stripepy cm` and `stripepy pd` will generate plots for a random 2.5 Mbp region, while `stripepy hist` will generate histograms using data from the entire genome.
+
+Example usage:
+
+```bash
+# Plot the pseudo-distribution over a region of interest
+stripepy plot pd results.hdf5 /tmp/pseudodistribution.png --region chr2:120100000-122100000
+
+# Plot the histograms using genome-wide data
+stripepy plot hist results.hdf5 /tmp/stripe_hist_gw.png
+
+# Plot the Hi-C matrix
+stripepy plot cm 4DNFI9GMP2J8.mcool 10000 /tmp/matrix.png
+
+# Plot the Hi-C matrix higlighting the stripe seeds
+stripepy plot cm 4DNFI9GMP2J8.mcool 10000 /tmp/matrix_with_seeds.png --stripepy-hdf5 results.hdf5 --highlight-seeds
+
+# Plot the Hi-C matrix higlighting the architectural stripes
+stripepy plot cm 4DNFI9GMP2J8.mcool 10000 /tmp/matrix_with_stripes.png --stripepy-hdf5 results.hdf5 --highlight-stripes
+```
+
+Some example plots generated with `stripepy plot` can be found in file `stripepy-plot-test-images.tar.xz` from [doi.org/10.5281/zenodo.14283921](https://doi.org/10.5281/zenodo.14283921)
+
 ## Getting help
 
 For any issues regarding StripePy installation, walkthrough, and output interpretation please open a [discussion](https://github.com/paulsengroup/StripePy/discussions) on GitHub.

src/stripepy/IO.py

@@ -695,43 +695,6 @@ def format_ticks(ax, x=True, y=True, rotate=True):
         ax.tick_params(axis="x", rotation=45)
 
 
-def HiC(I, RoI, plot_in_bp=False, output_folder=None, file_name=None, title=None, compactify=False):
-    """
-    :param I:                   Hi-C matrix to be plotted as image and saved
-    :param  RoI:                refers to the Region of Interest [RoI[0], RoI[1]]x[RoI[2], RoI[3]]
-                                (e.g., in genomic coordinates)
-    :param plot_in_bp:          if True, labels are set in base pairs w.r.t. the genomic interval in RoI;
-                                if False, labels are set in pixel coordinates
-    :param output_folder:       path to folder where to save the image
-    :param file_name:           name of the file to be created
-    :param title:               title to give to the image
-    :param compactify:          if False, it adds axes ticks, color bars
-    :return:                    -
-    """
-    if output_folder is None or file_name is None:
-        return
-
-    fig, ax = plt.subplots(1, 1)
-    img = ax.matshow(I, vmax=np.amax(I), extent=(RoI[0], RoI[1], RoI[3], RoI[2]), cmap=fruit_punch)
-    if plot_in_bp:
-        format_ticks(ax)
-
-    if compactify is True:
-        plt.axis("off")
-    else:
-        # plt.colorbar(img)
-        if title is not None:
-            fig.suptitle(title)
-
-    fig.set_dpi(256)
-    plt.axis("scaled")
-    fig.tight_layout()
-
-    plt.savefig(f"{output_folder}/{file_name}", bbox_inches="tight")
-
-    plt.close()
-
-
 def pseudodistrib(
     pseudo_distrib,
     IoI,
@@ -743,6 +706,7 @@ def pseudodistrib(
     title=None,
     display=False,
 ):
+    # TODO remove
     """
     :param pseudo_distrib:          1D ndarray representing a scalar function sampled over a uniform mesh
     :param  IoI:                    refers to the Interval of Interest [IoI[0], IoI[1]] (e.g., in genomic coordinates)
@@ -808,6 +772,7 @@ def pseudodistrib(
 def pseudodistrib_and_HIoIs(
     pseudo_distrib, IoIs, resolution, colors=None, output_folder=None, file_name=None, title=None, display=False
 ):
+    # TODO remove
     """
     :param pseudo_distrib:          1D ndarray representing a scalar function sampled over a uniform mesh
     :param  IoIs:                   list of lists, where the innermost lists are pairs of coordinates; the first pair
@@ -824,7 +789,6 @@ def pseudodistrib_and_HIoIs(
     :param display:                 if False, it does not display the plot
     :return:                        -
     """
-
     fig, ax = plt.subplots(1, 1)
     for IoI, color in zip(IoIs, colors):
         ax.plot(
@@ -834,21 +798,17 @@ def pseudodistrib_and_HIoIs(
             linewidth=0.5,
             linestyle="solid",
         )
-
     ax.xaxis.set_major_formatter(EngFormatter("b"))
     if no_frills_in_images is False:
         if title is not None:
             fig.suptitle(title)
-
     ax.set_xlabel("genomic coordinates (bp)")
     ax.set_ylabel("pseudo-distribution")
     fig.tight_layout()
     ax.grid(True)
     # plt.axis('scaled')
-
     if output_folder is not None and file_name is not None:
         plt.savefig(f"{output_folder}/{file_name}", bbox_inches="tight")
-
     if display:
         plt.show()
     plt.close()
@@ -865,7 +825,7 @@ def HiC_and_sites(
     file_name=None,
     title=None,
     display=False,
-):
+):  # TODO remove
     """
     :param I:                  Hi-C matrix to be plotted as image and saved
     :param sites:              list of locations of interest
@@ -935,6 +895,7 @@ def HiC_and_HIoIs(
     title=None,
     display=False,
 ):
+    # TODO: remove
     """
     :param I:                   Hi-C matrix to be plotted as image and saved
     :param HIoIs:               list of lists, where the innermost lists are pairs of elements
@@ -1036,6 +997,7 @@ def plot_stripes(
     title=None,
     display=False,
 ):
+    # TODO remove
     """
     :param I:                  Hi-C matrix to be plotted as image and saved
     :param LT_HIoIs:           Horizontal Intervals of Interest (lower-triangular part)
@@ -1280,22 +1242,3 @@ def plot_stripes_and_peaks(
     if display:
         plt.show()
     plt.close()
-
-
-def save_candidates_bedpe(HIoIs, VIoIs, resolution, chr, output_folder, file_name):
-    """
-    :param HIoIs:              Horizontal Intervals of Interest
-    :param VIoIs:              Vertical Intervals of Interest
-    :param resolution:         resolution
-    :param chr:                chromosome
-    :param output_folder:      path to folder where to save the image
-    :param file_name:          name of the file to be created
-    :return:                   -
-    """
-
-    with open(f"{output_folder}/{file_name}", "w") as f:
-        for HIoI, VIoI in zip(HIoIs, VIoIs):
-            f.write(
-                f"{chr}\t{resolution * HIoI[0]}\t{resolution * HIoI[1]}\t"
-                f"{chr}\t{resolution * VIoI[0]}\t{resolution * VIoI[1]}\n"
-            )

src/stripepy/cli/call.py

@@ -34,14 +34,14 @@ def run(
     start_global_time = time.time()
 
     # Data loading:
-    f, chr_starts, chr_ends, bp_lengths = others.cmap_loading(configs_input["contact-map"], configs_input["resolution"])
+    f, chr_starts, chr_ends, bp_lengths = others.cmap_loading(configs_input["contact_map"], configs_input["resolution"])
 
     # Remove existing folders:
     # configs_output["output_folder"] = (
-    #     f"{configs_output['output_folder']}/{configs_input['contact-map'].stem}/{configs_input['resolution']}"
+    #     f"{configs_output['output_folder']}/{configs_input['contact_map'].stem}/{configs_input['resolution']}"
     # )
     configs_output["output_folder"] = (
-        configs_output["output_folder"] / configs_input["contact-map"].stem / str(configs_input["resolution"])
+        configs_output["output_folder"] / configs_input["contact_map"].stem / str(configs_input["resolution"])
     )
     IO.remove_and_create_folder(configs_output["output_folder"], configs_output["force"])