diff --git a/app/Main.hs b/app/Main.hs
--- a/app/Main.hs
+++ b/app/Main.hs
@@ -78,6 +78,8 @@
                                        <?> "([1] | INT) Number of eigenvectors to use while clustering with kmeans. Takes from the first eigenvector. Recommended to start at 2 and work up from there if needed."
                        , minModularity  :: Maybe Double
                                        <?> "([0] | DOUBLE) Minimum modularity to be over to continue recursion."
+                       , numRuns        :: Maybe Int
+                                       <?> "([Nothing] | INT) Number of runs for permutation test at each split. If present, returns p-value for each split in the tree."
                        , eigenGroup     :: Maybe String
                                        <?> "([SignGroup] | KMeansGroup) Whether to group the eigenvector using the sign or kmeans while clustering. While the default is sign, kmeans may be more accurate (but starting points are arbitrary)."
                        , separateComponents :: Bool
@@ -117,6 +119,7 @@
         minSize'            = fmap MinSize . unHelpful . minSize $ opts
         numEigen'           = fmap NumEigen . unHelpful . numEigen $ opts
         minModularity'      = fmap Q . unHelpful . minModularity $ opts
+        numRuns'            = unHelpful . numRuns $ opts
         eigenGroup'         =
           maybe SignGroup (readOrErr "Cannot read --eigen-group")
             . unHelpful
@@ -143,40 +146,40 @@
             Dense -> do
                 (items, mat) <- readDenseAdjMatrix decodeOpt stdin
 
-                let cluster items = clusteringTreeToGenericClusteringTree
+                let cluster items = fmap clusteringTreeToGenericClusteringTree
                                   . HD.hierarchicalSpectralClusterAdj
                                       eigenGroup'
                                       (fmap unNumEigen numEigen')
                                       (fmap unMinSize minSize')
                                       minModularity'
+                                      numRuns'
                                       items
 
-                return $
-                    if separateComponents'
-                        then Multiple
-                           . fmap (uncurry cluster)
-                           . getComponentMatsItems items
-                           $ mat
-                        else Single $ cluster items mat
+                if separateComponents'
+                    then fmap Multiple
+                       . mapM (uncurry cluster)
+                       . getComponentMatsItems items
+                       $ mat
+                    else Single <$> cluster items mat
             Sparse -> do
                 -- (items, mat) <- readEigenSparseAdjMatrix decodeOpt stdin
                 (items, mat) <- readSparseAdjMatrix decodeOpt stdin
 
-                let cluster items = clusteringTreeToGenericClusteringTree
+                let cluster items = fmap clusteringTreeToGenericClusteringTree
                                   . HS.hierarchicalSpectralClusterAdj
                                       eigenGroup'
                                       (fmap unNumEigen numEigen')
                                       (fmap unMinSize minSize')
                                       minModularity'
+                                      numRuns'
                                       items
 
-                return $
-                    if separateComponents'
-                        then Multiple
-                           . fmap (uncurry cluster)
-                           . getComponentMatsItems items
-                           $ mat
-                        else Single $ cluster items mat
+                if separateComponents'
+                    then fmap Multiple
+                       . mapM (uncurry cluster)
+                       . getComponentMatsItems items
+                       $ mat
+                    else Single <$> cluster items mat
 
             (Premade file) ->
                 fmap (either error Single . A.eitherDecode) . B.readFile $ file
diff --git a/hierarchical-spectral-clustering.cabal b/hierarchical-spectral-clustering.cabal
--- a/hierarchical-spectral-clustering.cabal
+++ b/hierarchical-spectral-clustering.cabal
@@ -1,6 +1,6 @@
 cabal-version: >=1.10
 name: hierarchical-spectral-clustering
-version: 0.4.1.3
+version: 0.5.0.0
 license: GPL-3
 license-file: LICENSE
 copyright: 2019 Gregory W. Schwartz
@@ -43,6 +43,7 @@
         managed >=1.0.6,
         modularity >=0.2.1.1,
         mtl >=2.2.2,
+        mwc-random >=0.13.6.0,
         safe >=0.3.17,
         sparse-linear-algebra >=0.3.1,
         spectral-clustering >=0.3.1.2,
diff --git a/src/Math/Clustering/Hierarchical/Spectral/Dense.hs b/src/Math/Clustering/Hierarchical/Spectral/Dense.hs
--- a/src/Math/Clustering/Hierarchical/Spectral/Dense.hs
+++ b/src/Math/Clustering/Hierarchical/Spectral/Dense.hs
@@ -53,37 +53,48 @@
                             -> Maybe NumEigen
                             -> Maybe Int
                             -> Maybe Q
+                            -> Maybe Runs
                             -> Items a
                             -> Either FeatureMatrix B
-                            -> ClusteringTree a
-hierarchicalSpectralCluster eigenGroup normFlag numEigenMay minSizeMay minModMay initItems initMat =
+                            -> IO (ClusteringTree a)
+hierarchicalSpectralCluster eigenGroup normFlag numEigenMay minSizeMay minModMay runsMay initItems initMat =
     go initItems initB
   where
     initB = either (getB normFlag) id $ initMat
     minMod   = fromMaybe (Q 0) minModMay
     minSize  = fromMaybe 1 minSizeMay
     numEigen = fromMaybe 1 numEigenMay
-    go :: Items a -> B -> ClusteringTree a
+    go :: Items a -> B -> IO (ClusteringTree a)
     go !items !b =
         if (H.rows $ unB b) > 1
             && hasMultipleClusters clusters
             && ngMod > minMod
             && H.rows (unB left) >= minSize
             && H.rows (unB right) >= minSize
-            then
-                Node { rootLabel = vertex
-                     , subForest = [ go (subsetVector items leftIdxs) left
-                                   , go (subsetVector items rightIdxs) right
-                                   ]
-                     }
-
-            else
-                Node {rootLabel = vertex, subForest = []}
+            then do
+                goLeft <- go (subsetVector items leftIdxs) left
+                goRight <- go (subsetVector items rightIdxs) right
+                vertex <- getVertex
+                return $ Node { rootLabel = vertex
+                              , subForest = [goLeft, goRight]
+                              }
+            else do
+                vertex <- getVertex
+                return $ Node {rootLabel = vertex, subForest = []}
       where
-        vertex      = ClusteringVertex
-                        { _clusteringItems = items
-                        , _ngMod = ngMod
-                        }
+        getVertex = do
+          pVal <- mapM (\ x -> permutationTest
+                                (\l m -> unQ $ getBModularity l m)
+                                x
+                                clusters
+                                b
+                       )
+                  runsMay
+          return $ ClusteringVertex
+                 { _clusteringItems = items
+                 , _ngMod = ngMod
+                 , _pVal = pVal
+                 }
         clusters :: H.Vector Double
         clusters = spectralClustering eigenGroup b
         spectralClustering :: EigenGroup -> B -> H.Vector Double
@@ -104,43 +115,56 @@
                             -> Maybe NumEigen
                             -> Maybe Int
                             -> Maybe Q
+                            -> Maybe Runs
                             -> Items a
                             -> AdjacencyMatrix
-                            -> ClusteringTree a
-hierarchicalSpectralClusterAdj !eigenGroup !numEigenMay !minSizeMay !minModMay !items !adjMat =
+                            -> IO (ClusteringTree a)
+hierarchicalSpectralClusterAdj !eigenGroup !numEigenMay !minSizeMay !minModMay !runsMay !items !adjMat =
     if H.rows adjMat > 1
         && hasMultipleClusters clusters
         && ngMod > minMod
         && H.rows left >= minSize
         && H.rows right >= minSize
-        then
-            Node { rootLabel = vertex
-                 , subForest =
-                    [ hierarchicalSpectralClusterAdj
+        then do
+            goLeft <- hierarchicalSpectralClusterAdj
                         eigenGroup
                         numEigenMay
                         minSizeMay
                         minModMay
+                        runsMay
                         (subsetVector items leftIdxs)
                         left
-                    , hierarchicalSpectralClusterAdj
+            goRight <- hierarchicalSpectralClusterAdj
                         eigenGroup
                         numEigenMay
                         minSizeMay
                         minModMay
+                        runsMay
                         (subsetVector items rightIdxs)
                         right
-                    ]
-                 }
-        else
-            Node {rootLabel = vertex, subForest = []}
+            vertex <- getVertex
+            return $ Node { rootLabel = vertex , subForest = [goLeft, goRight] }
+        else do
+            vertex <- getVertex
+            return $ Node {rootLabel = vertex, subForest = []}
   where
     minMod      = fromMaybe (Q 0) minModMay
     minSize     = fromMaybe 1 minSizeMay
     numEigen    = fromMaybe 1 numEigenMay
-    vertex      = ClusteringVertex { _clusteringItems = items
-                                   , _ngMod = ngMod
-                                   }
+    getVertex = do
+      pVal <- mapM
+                (\ x -> permutationTest
+                        (\l m -> unQ $ getModularity l m)
+                        x
+                        clusters
+                        adjMat
+                )
+                runsMay
+      return $ ClusteringVertex
+             { _clusteringItems = items
+             , _ngMod = ngMod
+             , _pVal = pVal
+             }
     clusters    = spectralClustering eigenGroup adjMat
     spectralClustering :: EigenGroup -> AdjacencyMatrix -> H.Vector Double
     spectralClustering SignGroup   = spectralClusterNorm
diff --git a/src/Math/Clustering/Hierarchical/Spectral/Sparse.hs b/src/Math/Clustering/Hierarchical/Spectral/Sparse.hs
--- a/src/Math/Clustering/Hierarchical/Spectral/Sparse.hs
+++ b/src/Math/Clustering/Hierarchical/Spectral/Sparse.hs
@@ -29,6 +29,7 @@
 import qualified Data.Sparse.Common as S
 import qualified Data.Vector as V
 import qualified Data.Vector.Storable as VS
+import qualified Data.Vector.Unboxed as U
 import qualified Numeric.LinearAlgebra.Sparse as S
 
 -- Local
@@ -55,37 +56,48 @@
                             -> Maybe NumEigen
                             -> Maybe Int
                             -> Maybe Q
+                            -> Maybe Runs
                             -> Items a
                             -> Either FeatureMatrix B
-                            -> ClusteringTree a
-hierarchicalSpectralCluster eigenGroup normFlag numEigenMay minSizeMay minModMay initItems initMat =
+                            -> IO (ClusteringTree a)
+hierarchicalSpectralCluster eigenGroup normFlag numEigenMay minSizeMay minModMay runsMay initItems initMat =
     go initItems initB
   where
     initB = either (getB normFlag) id $ initMat
     minMod   = fromMaybe (Q 0) minModMay
     minSize  = fromMaybe 1 minSizeMay
     numEigen = fromMaybe 1 numEigenMay
-    go :: Items a -> B -> ClusteringTree a
+    go :: Items a -> B -> IO (ClusteringTree a)
     go !items !b =
         if (S.nrows $ unB b) > 1
             && hasMultipleClusters clusters
             && ngMod > minMod
             && S.nrows (unB left) >= minSize
             && S.nrows (unB right) >= minSize
-            then
-                Node { rootLabel = vertex
-                     , subForest = [ go (subsetVector items leftIdxs) left
-                                   , go (subsetVector items rightIdxs) right
-                                   ]
-                     }
-
-            else
-                Node {rootLabel = vertex, subForest = []}
+            then do
+                goLeft <- go (subsetVector items leftIdxs) left
+                goRight <- go (subsetVector items rightIdxs) right
+                vertex <- getVertex
+                return $ Node { rootLabel = vertex
+                              , subForest = [goLeft, goRight]
+                              }
+            else do
+                vertex <- getVertex
+                return $ Node {rootLabel = vertex, subForest = []}
       where
-        vertex      = ClusteringVertex
-                        { _clusteringItems = items
-                        , _ngMod = ngMod
-                        }
+        getVertex = do
+          pVal <- mapM (\ x -> permutationTestSparse
+                                (\l m -> unQ $ getBModularity l m)
+                                x
+                                clusters
+                                b
+                       )
+                  runsMay
+          return $ ClusteringVertex
+                 { _clusteringItems = items
+                 , _ngMod = ngMod
+                 , _pVal = pVal
+                 }
         clusters :: S.SpVector Double
         clusters = spectralClustering eigenGroup b
         spectralClustering :: EigenGroup -> B -> S.SpVector Double
@@ -118,36 +130,47 @@
                                -> Maybe NumEigen
                                -> Maybe Int
                                -> Maybe Q
+                               -> Maybe Runs
                                -> Items a
                                -> AdjacencyMatrix
-                               -> ClusteringTree a
-hierarchicalSpectralClusterAdj eigenGroup numEigenMay minSizeMay minModMay initItems initMat =
+                               -> IO (ClusteringTree a)
+hierarchicalSpectralClusterAdj eigenGroup numEigenMay minSizeMay minModMay runsMay initItems initMat =
     go initItems initMat
   where
     minMod      = fromMaybe (Q 0) minModMay
     minSize     = fromMaybe 1 minSizeMay
     numEigen    = fromMaybe 1 numEigenMay
-    go :: Items a -> AdjacencyMatrix -> ClusteringTree a
+    go :: Items a -> AdjacencyMatrix -> IO (ClusteringTree a)
     go !items !mat =
         if S.nrows mat > 1
             && hasMultipleClusters clusters
             && ngMod > minMod
             && S.nrows left >= minSize
             && S.nrows right >= minSize
-            then
-                Node { rootLabel = vertex
-                     , subForest = [ go (subsetVector items leftIdxs) left
-                                   , go (subsetVector items rightIdxs) right
-                                   ]
-                     }
-
-            else
-                Node {rootLabel = vertex, subForest = []}
+            then do
+                goLeft <- go (subsetVector items leftIdxs) left
+                goRight <- go (subsetVector items rightIdxs) right
+                vertex <- getVertex
+                return $ Node { rootLabel = vertex
+                              , subForest = [ goLeft, goRight]
+                              }
+            else do
+                vertex <- getVertex
+                return $ Node {rootLabel = vertex, subForest = []}
       where
-        vertex      = ClusteringVertex
-                        { _clusteringItems = items
-                        , _ngMod = ngMod
-                        }
+        getVertex = do
+          pVal <- mapM (\ x -> permutationTestSparse
+                                (\l m -> unQ $ getModularity l m)
+                                x
+                                clusters
+                                mat
+                       )
+                  runsMay
+          return $ ClusteringVertex
+                 { _clusteringItems = items
+                 , _ngMod = ngMod
+                 , _pVal = pVal
+                 }
         clusters :: S.SpVector Double
         clusters = spectralClustering eigenGroup mat
         spectralClustering :: EigenGroup -> AdjacencyMatrix -> S.SpVector Double
diff --git a/src/Math/Clustering/Hierarchical/Spectral/Test.hs b/src/Math/Clustering/Hierarchical/Spectral/Test.hs
--- a/src/Math/Clustering/Hierarchical/Spectral/Test.hs
+++ b/src/Math/Clustering/Hierarchical/Spectral/Test.hs
@@ -86,6 +86,7 @@
                     Nothing
                     Nothing
                     Nothing
+                    Nothing
                     exampleItems
                     (Left $ exampleMatrix 3 exampleData)
 
@@ -95,6 +96,7 @@
                     (Just 2)
                     Nothing
                     Nothing
+                    Nothing
                     exampleItems
                     (Left $ exampleMatrix 3 exampleData)
 
@@ -103,6 +105,7 @@
                         Nothing
                         Nothing
                         Nothing
+                        Nothing
                         exampleItems
                         adjacencyExample
 
@@ -111,6 +114,7 @@
                         (Just 2)
                         Nothing
                         Nothing
+                        Nothing
                         exampleItems
                         adjacencyExample
 
@@ -136,6 +140,7 @@
                           Nothing
                           Nothing
                           Nothing
+                          Nothing
                           exampleItems
                           denseAdjacencyExample
 
@@ -144,6 +149,7 @@
                           (Just 2)
                           Nothing
                           Nothing
+                          Nothing
                           exampleItems
                           denseAdjacencyExample
 
@@ -153,6 +159,7 @@
                         Nothing
                         Nothing
                         Nothing
+                        Nothing
                         exampleItems
                         (Left denseFeatureExample)
 
@@ -160,6 +167,7 @@
                         KMeansGroup
                         True
                         (Just 2)
+                        Nothing
                         Nothing
                         Nothing
                         exampleItems
diff --git a/src/Math/Clustering/Hierarchical/Spectral/Types.hs b/src/Math/Clustering/Hierarchical/Spectral/Types.hs
--- a/src/Math/Clustering/Hierarchical/Spectral/Types.hs
+++ b/src/Math/Clustering/Hierarchical/Spectral/Types.hs
@@ -22,6 +22,7 @@
     , getClusterItemsGenericTree
     , Q (..)
     , NumEigen (..)
+    , Runs (..)
     ) where
 
 -- Remote
@@ -42,17 +43,20 @@
 type ClusteringTree a = Tree (ClusteringVertex a)
 type GenericClusteringTree a = Tree (GenericClusteringVertex a)
 type NumEigen = Int
+type Runs = Int
 
 data EigenGroup = SignGroup | KMeansGroup deriving (Read, Show, Generic)
 
 data ClusteringVertex a = ClusteringVertex
     { _clusteringItems  :: !(Items a)
     , _ngMod            :: !Q
+    , _pVal           :: !(Maybe Double)
     } deriving (Eq, Ord, Read, Show, Generic)
 
 data GenericClusteringVertex a = GenericClusteringVertex
     { _item :: !(Maybe (Items a))
     , _distance :: !(Maybe Double)
+    , _pValue :: !(Maybe Double)
     } deriving (Eq, Ord, Read, Show, Generic)
 
 -- | Convert a ClusteringTree to a Dendrogram. Modularity is the distance.
@@ -95,6 +99,7 @@
       Node { rootLabel = ( GenericClusteringVertex
                              { _item = Just $ _clusteringItems n
                              , _distance = Nothing
+                             , _pValue = Nothing
                              }
                          )
            , subForest = []
@@ -103,6 +108,7 @@
       Node { rootLabel = ( GenericClusteringVertex
                             { _item = Nothing
                             , _distance = Just . unQ . _ngMod $ n
+                            , _pValue = _pVal n
                             }
                          )
            , subForest = fmap go xs
@@ -115,7 +121,7 @@
 -- | Gather clusters (leaves) from the tree.
 getClusterItemsTree :: ClusteringTree a -> [Items a]
 getClusterItemsTree = fmap _clusteringItems . leaves
-                    
+
 -- | Gather clusters (leaves) from the generic tree.
 getClusterItemsGenericTree :: GenericClusteringTree a -> [Items a]
 getClusterItemsGenericTree = catMaybes . fmap _item . leaves
diff --git a/src/Math/Clustering/Hierarchical/Spectral/Utility.hs b/src/Math/Clustering/Hierarchical/Spectral/Utility.hs
--- a/src/Math/Clustering/Hierarchical/Spectral/Utility.hs
+++ b/src/Math/Clustering/Hierarchical/Spectral/Utility.hs
@@ -8,12 +8,20 @@
 
 module Math.Clustering.Hierarchical.Spectral.Utility
   ( subsetVector
+  , permutationTest
+  , permutationTestSparse
   ) where
 
 -- Remote
+import Data.List (genericLength)
 import Data.Maybe (fromMaybe)
 import qualified Data.Foldable as F
+import qualified Data.Sparse.Common as S
 import qualified Data.Vector as V
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Generic as G
+import qualified System.Random.MWC as R
+import qualified System.Random.MWC.Distributions as R
 
 -- Local
 
@@ -26,3 +34,40 @@
                       $ xs V.!? i
                       ) : acc
                     ) []
+
+
+-- | Permutation test returning p-value.
+permutationTest :: (G.Vector v a)
+                => (v a -> b -> Double)
+                -> Int
+                -> v a
+                -> b
+                -> IO Double
+permutationTest f runs labs samples = do
+  let obs = f labs samples
+      run r = do
+        g <- R.initialize . V.singleton $ fromIntegral r
+        ls <- R.uniformShuffle labs g
+        return $ f ls samples
+
+  n <- genericLength . filter (>= obs) <$> mapM run [1..runs]
+
+  return $ n / fromIntegral runs
+
+-- | Permutation test returning p-value from sparse vector.
+permutationTestSparse :: (S.SpVector Double -> b -> Double)
+                -> Int
+                -> S.SpVector Double
+                -> b
+                -> IO Double
+permutationTestSparse f runs labsInit samples = do
+  let obs = f labsInit samples
+      labs = U.fromList $ S.toDenseListSV labsInit
+      run r = do
+        g <- R.initialize . V.singleton $ fromIntegral r
+        ls <- S.vr . U.toList <$> R.uniformShuffle labs g
+        return $ f ls samples
+
+  n <- genericLength . filter (>= obs) <$> mapM run [1..runs]
+
+  return $ n / fromIntegral runs