diff --git a/nireports/interfaces/__init__.py b/nireports/interfaces/__init__.py
index e190ef32..5391f075 100644
--- a/nireports/interfaces/__init__.py
+++ b/nireports/interfaces/__init__.py
@@ -24,11 +24,12 @@
 
 from nireports.interfaces.fmri import FMRISummary
 from nireports.interfaces.mosaic import PlotContours, PlotMosaic, PlotSpikes
-from nireports.interfaces.nuisance import CompCorVariancePlot, ConfoundsCorrelationPlot
+from nireports.interfaces.nuisance import CompCorVariancePlot, ConfoundsCorrelationPlot, RaincloudPlot
 
 __all__ = (
     "CompCorVariancePlot",
     "ConfoundsCorrelationPlot",
+    "RaincloudPlot",
     "FMRISummary",
     "PlotContours",
     "PlotMosaic",
diff --git a/nireports/interfaces/nuisance.py b/nireports/interfaces/nuisance.py
index 734767f4..d88d5858 100644
--- a/nireports/interfaces/nuisance.py
+++ b/nireports/interfaces/nuisance.py
@@ -30,9 +30,13 @@
     isdefined,
     traits,
 )
+<<<<<<< HEAD
 from nipype.utils.filemanip import fname_presuffix
 
 from nireports.reportlets.nuisance import confounds_correlation_plot
+=======
+from nireports.reportlets.nuisance import confounds_correlation_plot, plot_raincloud
+>>>>>>> 80753e5 (ENH: Add raincloud plot capabilities)
 from nireports.reportlets.xca import compcor_variance_plot
 
 
@@ -142,3 +146,70 @@ def _run_interface(self, runtime):
             ignore_initial_volumes=self.inputs.ignore_initial_volumes,
         )
         return runtime
+
+
+class _RaincloudPlotInputSpec(BaseInterfaceInputSpec):
+    data_file = File(
+        exists=True, mandatory=True, desc="File containing the data"
+    )
+    out_file = traits.Either(
+        None, File, value=None, usedefault=True, desc="Path to save plot"
+    )
+    group_label = traits.Str(
+        "group_label",
+        mandatory=True,
+        desc="Group name of interest",
+    )
+    feature = traits.Str(
+        "feature",
+        mandatory=True,
+        desc="Feature of interest",
+    )
+    palette = traits.Str(
+        "Set2",
+        usedefault=True,
+        desc="Color palette name",
+    )
+    orient = traits.Str(
+        "v",
+        usedefault=True,
+        desc="Orientation",
+    )
+    density = traits.Bool(
+        True,
+        usedefault=True,
+        desc="``True`` to plot the density",
+    )
+
+
+class _RaincloudPlotOutputSpec(TraitedSpec):
+    out_file = File(exists=True, desc="Path to saved plot")
+
+
+class RaincloudPlot(SimpleInterface):
+    """Plot a raincloud of values."""
+
+    input_spec = _RaincloudPlotInputSpec
+    output_spec = _RaincloudPlotOutputSpec
+
+    def _run_interface(self, runtime):
+        if self.inputs.out_file is None:
+            self._results["out_file"] = fname_presuffix(
+                self.inputs.data_file,
+                suffix="_raincloud.svg",
+                use_ext=False,
+                newpath=runtime.cwd,
+            )
+        else:
+            self._results["out_file"] = self.inputs.out_file
+        plot_raincloud(
+            data_file=self.inputs.data_file,
+            group_label=self.inputs.group_label,
+            feature=self.inputs.feature,
+            palette=self.inputs.palette,
+            orient=self.inputs.orient,
+            density=self.inputs.density,
+            output_file=self._results["out_file"],
+            **kwargs,
+        )
+        return runtime
diff --git a/nireports/reportlets/modality/dwi.py b/nireports/reportlets/modality/dwi.py
index 117b16f5..b1e9216b 100644
--- a/nireports/reportlets/modality/dwi.py
+++ b/nireports/reportlets/modality/dwi.py
@@ -405,3 +405,13 @@ def plot_gradients(
         plt.suptitle(title)
 
     return ax
+
+
+def plot_tissue_values(pvms, fa):
+
+    # Plot FA values in voxels for GM, WM, CSF as boxplots
+    # ToDo
+    # Define colors for CST, GM, WM or read them from somehwere
+
+    df = pd.DataFrame(pvms, columns=["group, score"])
+    return plot_raincloud(df)
diff --git a/nireports/reportlets/nuisance.py b/nireports/reportlets/nuisance.py
index 444104ec..9e186852 100644
--- a/nireports/reportlets/nuisance.py
+++ b/nireports/reportlets/nuisance.py
@@ -30,6 +30,7 @@
 import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
+import pandas as pd
 import seaborn as sns
 from matplotlib.backends.backend_pdf import FigureCanvasPdf as FigureCanvas
 from matplotlib.colorbar import ColorbarBase
@@ -626,7 +627,6 @@ def confoundplot(
     cutoff=None,
     ylims=None,
 ):
-    import seaborn as sns
 
     # Define TR and number of frames
     notr = False
@@ -856,8 +856,6 @@ def confounds_correlation_plot(
     output_file: :obj:`str`
         The file where the figure is saved.
     """
-    import pandas as pd
-    import seaborn as sns
 
     confounds_data = pd.read_table(confounds_file)
 
@@ -932,3 +930,158 @@ def confounds_correlation_plot(
         figure = None
         return output_file
     return [ax0, ax1], gs
+
+
+def plot_raincloud(
+    data_file,
+    group_name,
+    feature,
+    palette="Set2",
+    orient="v",
+    density=True,
+    figure=None,
+    output_file=None,
+    **kwargs,
+):
+    """
+    Generate a raincloud plot with the data points corresponding to the
+    ``feature`` value contained in the data file.
+
+    Parameters
+    ----------
+    data_file : :obj:`str`
+        File containing the data of interest.
+    figure : :obj:`matplotlib.pyplot.figure` or None
+        Existing figure on which to plot.
+    group_name : :obj:`str`
+        The group name of interest to be plot.
+    feature : :obj:`str`
+        The feature of interest to be plot.
+    palette : :obj:`str`, optional
+        Color palette name provided to :func:`sns.stripplot`.
+    orient : :obj:`str`, optional
+        Plot orientation (``v`` or ``h``).
+    density : :obj:`bool`, optional
+        ``True`` to plot the density of the data points.
+    output_file : :obj:`str` or :obj:`None`
+        Path where the output figure should be saved. If this is not defined,
+        then the plotting axes will be returned instead of the saved figure
+        path.
+    kwargs : :obj:`dict`
+        Extra args given to :func:`sns.violinplot`, :func:`sns.stripplot` and
+        :func:`sns.boxplot`.
+
+    Returns
+    -------
+    axes and gridspec
+        Plotting axes and gridspec. Returned only if ``output_file`` is ``None``.
+    output_file : :obj:`str`
+        The file where the figure is saved.
+    """
+
+    # ToDo
+    # Think how data will come or why not accept the df directly ?
+    # df = pd.read_table(data_file)
+    df = pd.read_csv(data_file, sep=r"[\t\s]+", engine="python")
+
+    if figure is None:
+        plt.figure(figsize=(7, 5))
+
+    gs = GridSpec(1, 1)
+    ax = plt.subplot(gs[0, 0])
+
+    # ToDo
+    # Make all the below go to the kwargs. Note that they take different kwargs
+    # so we will need 3 dictionaries: one fro violinplot, one for stripplot, and
+    # the other for the boxplot, and one for the general style ?
+    sns.set(style="whitegrid", font_scale=2)
+
+    x = feature
+    y = group_name
+    # Swap x/y if the requested orientation is vertical
+    if orient == "v":
+        x = group_name
+        y = feature
+
+    # Plot the violin plots
+    if density:
+        ax = sns.violinplot(
+            x=x,
+            y=y,
+            data=df,
+            hue=group_name,
+            dodge=False,
+            palette=palette,
+            scale="width",
+            inner=None,
+            orient=orient,
+        )
+
+        # Cut half of the violins
+        for violin in ax.collections:
+            bbox = violin.get_paths()[0].get_extents()
+            x0, y0, width, height = bbox.bounds
+            width_denom = 2
+            height_denom = 1
+            if orient == "h":
+                width_denom = 1
+                height_denom = 2
+            violin.set_clip_path(
+                plt.Rectangle(
+                    (x0, y0),
+                    width/width_denom,
+                    height/height_denom,
+                    transform=ax.transData)
+            )
+
+    # Add boxplots
+    width = 0.15
+    sns.boxplot(
+        x=x,
+        y=y,
+        data=df,
+        color="black",
+        width=width,
+        zorder=10,
+        showcaps=True,
+        boxprops={"facecolor": "none", "zorder": 10},
+        showfliers=True,
+        whiskerprops={"linewidth": 2, "zorder": 10},
+        saturation=1,
+        orient=orient,
+        ax=ax,
+    )
+
+    old_len_collections = len(ax.collections)
+
+    # Plot the data points as dots
+    sns.stripplot(
+        x=x,
+        y=y,
+        hue=group_name,
+        data=df,
+        palette=palette,
+        edgecolor="white",
+        size=3,
+        jitter=0.1,
+        zorder=0,
+        orient=orient,
+        ax=ax,
+    )
+
+    # Offset the dots that would be otherwise shadowed by the violins
+    if density:
+        offset = np.array([width, 0])
+        if orient == "h":
+            offset = np.array([0, width])
+        for dots in ax.collections[old_len_collections:]:
+            dots.set_offsets(dots.get_offsets() + offset)
+
+    if output_file is not None:
+        figure = plt.gcf()
+        figure.savefig(output_file, bbox_inches="tight")
+        plt.close(figure)
+        figure = None
+        return output_file
+
+    return ax, gs
diff --git a/nireports/tests/test_dwi.py b/nireports/tests/test_dwi.py
index 6e43fea4..7d189c8f 100644
--- a/nireports/tests/test_dwi.py
+++ b/nireports/tests/test_dwi.py
@@ -27,7 +27,7 @@
 import pytest
 from matplotlib import pyplot as plt
 
-from nireports.reportlets.modality.dwi import plot_dwi, plot_gradients
+from nireports.reportlets.modality.dwi import plot_dwi, plot_gradients, plot_tissue_values
 
 
 def test_plot_dwi(tmp_path, testdata_path, outdir):
@@ -69,3 +69,27 @@ def test_plot_gradients(tmp_path, testdata_path, dwi_btable, outdir):
 
     if outdir is not None:
         plt.savefig(outdir / f"{dwi_btable}.svg", bbox_inches="tight")
+
+
+def test_plot_tissue_values():
+
+    # Plot FA values in voxels for GM, WM, CSF as raincloud
+    # df = pd.DataFrame(pvms, columns=["group, score"])
+    data_file = ""
+    group_name = "tissue"
+    feature = "FA"
+    figure = None
+    output_file = None
+    kwargs = {}
+
+    plot_tissue_values(
+        data_file,
+        group_name,
+        feature,
+        figure=figure,
+        output_file=output_file,
+        **kwargs,
+    )
+
+    print("bat")
+>>>>>>> 80753e5 (ENH: Add raincloud plot capabilities)
diff --git a/nireports/tests/test_interfaces.py b/nireports/tests/test_interfaces.py
index d361dd6b..0f2e3eae 100644
--- a/nireports/tests/test_interfaces.py
+++ b/nireports/tests/test_interfaces.py
@@ -27,7 +27,7 @@
 
 import pytest
 
-from nireports.interfaces.nuisance import CompCorVariancePlot, ConfoundsCorrelationPlot
+from nireports.interfaces.nuisance import CompCorVariancePlot, ConfoundsCorrelationPlot, RaincloudPlot
 
 
 def _smoke_test_report(report_interface, artifact_name):
@@ -56,3 +56,24 @@ def test_ConfoundsCorrelationPlot(datadir, ignore_initial_volumes):
         ignore_initial_volumes=ignore_initial_volumes,
     )
     _smoke_test_report(cc_rpt, f"confounds_correlation_{ignore_initial_volumes}.svg")
+
+
+# @pytest.mark.parametrize("orient", ["h", "v"])
+# @pytest.mark.parametrize("density", (True, False))
+def test_RaincloudPlot(datadir,):
+    """Raincloud plot report test"""
+    data_file = os.path.join(datadir, "raincloud_test.tsv")
+    group_label =
+    feature = 
+    palette =
+    orient =
+    density =
+    rc_rpt = RaincloudPlot(
+        data_file=data_file,
+        group_label=group_label,
+        feature=feature,
+        palette=palette,
+        orient=orient,
+        density=density,
+    )
+    _smoke_test_report(rc_rpt, f"raincloud_{orient}_{density}.svg")
diff --git a/nireports/tests/test_reportlets.py b/nireports/tests/test_reportlets.py
index 671334e0..e720c791 100644
--- a/nireports/tests/test_reportlets.py
+++ b/nireports/tests/test_reportlets.py
@@ -35,7 +35,7 @@
 
 from nireports.reportlets.modality.func import fMRIPlot
 from nireports.reportlets.mosaic import plot_mosaic
-from nireports.reportlets.nuisance import plot_carpet
+from nireports.reportlets.nuisance import plot_carpet, plot_raincloud
 from nireports.reportlets.surface import cifti_surfaces_plot
 from nireports.reportlets.xca import compcor_variance_plot, plot_melodic_components
 from nireports.tools.timeseries import cifti_timeseries as _cifti_timeseries
@@ -368,3 +368,56 @@ def test_mriqc_plot_mosaic_2(tmp_path, testdata_path, outdir):
         maxrows=12,
         annotate=True,
     )
+
+
+# @pytest.mark.parametrize("orient", ["h", "v"])
+# @pytest.mark.parametrize("density", (True, False))
+def test_plot_raincloud():
+
+    values_label = "value"
+    group_label = "group"
+    data_file = "data.tsv"
+
+    def _generate_random_data():
+
+        rng = np.random.default_rng(1234)
+
+        # Create some random data in the [0.3, 1.0] interval
+        n_samples = 250
+        min_val = 0.3
+        max_val = 1.0
+        range_size = (max_val - min_val)
+        values_half1 = rng.random(n_samples) * range_size + min_val
+        # Create some random data in the [0.0, 0.6] interval
+        min_val = 0.0
+        max_val = 0.6
+        range_size = (max_val - min_val)
+        values_half2 = rng.random(n_samples) * range_size + min_val
+
+        values = np.hstack((values_half1, values_half2))
+
+        df = pd.DataFrame(values, columns=[values_label])
+
+        # Categorize values according to the values into group1 and group2
+        df[group_label] = np.where(df[values_label].between(min_val, max_val), "group1", "group2")
+
+        df.to_csv(data_file, sep="\t")
+
+    _generate_random_data()
+
+    palette = "Set2"
+    orient = "v"
+    density = True
+    kwargs = {}
+    output_file = "raincloud.png"
+
+    plot_raincloud(
+        data_file,
+        group_label,
+        values_label,
+        palette=palette,
+        orient=orient,
+        density=density,
+        output_file=output_file,
+        **kwargs,
+    )