Typically one performs some kind of clustering.

I personally often use hierarchical clustering based on the Z-transformed log2-normalized counts of all differentially-expressed genes. For this you first transform the normalized counts (e.g. from DESeq2 or edgeR) from your DEGs to log2-space.

log2.counts <- log2(norm.counts+1)

Then you transform to Z-score:

Z.counts <- t(scale(t(log2.counts)))

Eventually you cluster, e.g. with ward.D2 method:

hclust.counts <- hclust(d = dist(Z.counts), method = "ward.D2")

Plot a heatmap using the ComplexHeatmap package from Biconductor. There are plenty of way to make this heatmap nicer and more custom, please check the ComplexHeatmap documentation. This heatmap will give you a first glance over the patterns in your data. Check the row_order function to extract the genes that belong to each cluster, as explained in the package documentation.

library(ComplexHeatmap)
htmp <- Heatmap(matrix        = Z.counts, 
                cluster_rows  = hclust.counts)

Below the relevant code with some dummy data assuming y was the log2-normalized and Z-transformed counts, so that you get an idea how this works. The heatmap here is of course not meaningful because the count data are just random numbers.

y <- matrix( rpois(1000, lambda=5), nrow=200 )
h <- hclust(dist(y))
htmp <- Heatmap(matrix = y, 
                cluster_rows = h) 
draw(htmp)

With actual data, here from one of my projects, this could look like this:

The degPatterns() function of the DEGreport package sounds like it could fit the bill. See the vignette here.